1

I tried to implement an algorithm to calculate power towers modulo m. Below the procedure tower should calculate 2^3^...^14^15 (mod m) and tower2 should calculate 15^14^...^3^2 (mod m). But for m = 163 , tower2 produces a wrong answer. I found out that a immediate result is 0 and the procedure does not get this. Can anyone fix the error ?

The procedure powmod is implemented and works perfectly :

powmod(basis,exponent,modul)={if(exponent==0,hilf=1);if(exponent>0,bin=binary(exponent);hilf=basis;hilf=hilf-truncate(hilf/modul)*modul;for(stelle=2,length(bin),hilf=hilf^2;if(bin[stelle]==1,hilf=hilf*basis);hilf=hilf-truncate(hilf/modul)*modul));hilf}

? tower

%19 = (p,q,r)->if(q==0,hilf=1);if(q==1,hilf=p);if(q==2,hilf=powmod(p,p,r));if(q>
2,x=[];for(j=1,q,x=concat(x,r);r=eulerphi(r));hilf=14^15;forstep(j=13,2,-1,r=x[j
-1];if(r>=2,hilf=powmod(j,hilf,r);w=factorint(r);w=component(w,2);while(hilf<vec
max(w),hilf=hilf+r))));component(Mod(hilf,r),2)

? tower2

%20 = (p,q,r)->if(q==0,hilf=1);if(q==1,hilf=p);if(q==2,hilf=powmod(p,p,r));if(q>
2,x=[];for(j=1,q,x=concat(x,r);r=eulerphi(r));hilf=3^2;forstep(j=13,2,-1,r=x[j-1
];if(r>=2,hilf=powmod(17-j,hilf,r);w=factorint(r);w=component(w,2);while(hilf<ve
cmax(w),hilf=hilf+r))));component(Mod(hilf,r),2)
?
max taldykin
  • 12,459
  • 5
  • 45
  • 64
Peter
  • 214
  • 1
  • 10
  • 1
    is it possible to make the code a little bit more readable? – max taldykin Apr 23 '14 at 06:15
  • possible duplicate of [Is there an algorithm known for power towers modulo a number managing all cases?](http://stackoverflow.com/questions/21305782/is-there-an-algorithm-known-for-power-towers-modulo-a-number-managing-all-cases) – max taldykin Apr 23 '14 at 06:18

1 Answers1

0

The reason your code doesn't work is that you (recursively) compute x^n (mod r) as x^(n mod phi(r)) and this isn't true unless gcd(x,r) = 1.

Also, you don't need powmod since Mod(basis,modul)^expo is built-in. Here's a general possibility :

\\ x[1]^(x[2]^( ...^ x[#x])) mod m, assuming x[i] > 1 for all i
tower(x, m) =
{ my(f = factor(m), P = f[,1], E = f[,2]);
  chinese(vector(#P, i, towerp(x, P[i], E[i])));
}

towerp(x, p, e) =
{ my(q = p^e, i, t, v);

  if (#x == 0, return (Mod(1, q)));
  if (#x == 1, return (Mod(x[1], q)));
  if (v = valuation(x[1], p),
    t = x[#x]; i = #x;
    while (i > 1,
      if (t >= e, return (Mod(0, q)));
      t = x[i]^t; i--);
    if (t * v >= e, return (Mod(0, q)));
    return (Mod(x[1], q)^t);
  );
  Mod(x[1], q)^lift(tower(x[2..#x], (p-1)*p^e));
}

? tower([2..15], 163)
%1 = Mod(162, 163)
? tower(Vecrev([2..15]), 163)
%2 = Mod(16, 163)
K.B.
  • 861
  • 5
  • 14