Неправильный результат в PARI-реализации

Question

Я попытался реализовать алгоритм для расчета мощности башни по модулю m. Ниже процедуры башня должна вычислить 2^3^...^14^15 (мод м), а башня2 должна вычислить 15^14^...^3^2 (мод м). Но для m = 163 tower2 дает неправильный ответ. Я обнаружил, что немедленный результат равен 0, и процедура не получает это. Кто-нибудь может исправить ошибку?

Процедура powmod реализована и отлично работает:

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)
?

Answer 1

Причина, по которой ваш код не работает, заключается в том, что вы (рекурсивно) вычисляете x ^ n (mod r) как x ^ (nmod phi (r)), и это неверно, если gcd(x,r) = 1.

Кроме того, вам не нужен powmod, так как Mod(basis,modul)^expoвстроен. Вот общая возможность:

      \\ 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)