我想在 PARI/GP 中实现计算
a_1 ^ a_2 ^ ... ^ a_n (mod m)
它管理所有情况,尤其是在 phi 链中出现高权力的情况。
有谁知道这样的实现?
这是一种使用中国余数来确保模数是主要幂的可能性。这简化了在 gcd(x,m) 不为 1 的痛苦情况下 x^n mod m 的计算。代码假设 a_i > 1;大多数代码检查 p^a_1^a_2^...^a_n 对于素数 p 是否为 0 mod (p^e),同时避免溢出。
\\ 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[^1], (p-1)*p^e));
}
例如
? 5^(4^(3^2)) % 163 \\ direct computation, wouldn't scale
%1 = 158
? tower([5,4,3,2], 163)
%2 = Mod(158, 163)