How to evalute an exponential tower modulo a prime

Question

I want to find a fast algorithm to evaluate an expression like the following, where P is prime.

A ^ B ^ C ^ D ^ E mod P

Example:

(9 ^ (3 ^ (15 ^ (3 ^ 15)))) mod 65537 = 16134

The problem is the intermediate results can grow much too large to handle.

Answer 1

Basically the problem reduces to computing a^T mod m for given a, m and a term T that is ridiulously huge. However, we are able to evaluate T mod n with a given modulus n much faster than T . So we ask: "Is there an integer n, such that a^(T mod n) mod m = a^T mod m?"

Now if a and m are coprime, we know that n = phi(m) fulfills our condition according to Euler's theorem:

  a^T                               (mod m)
= a^((T mod phi(m)) + k * phi(m))   (mod m)    (for some k)
= a^(T mod phi(m)) * a^(k * phi(m)) (mod m)
= a^(T mod phi(m)) * (a^phi(m))^k   (mod m)
= a^(T mod phi(m)) * 1^k            (mod m)
= a^(T mod phi(m))                  (mod m)

If we can compute phi(m) (which is easy to do for example in O(m^(1/2)) or if we know the prime factorization of m), we have reduced the problem to computing T mod phi(m) and a simple modular exponentiation.

What if a and m are not coprime? The situation is not as pleasant as before, since there might not be a valid n with the property a^T mod m = a^(T mod n) mod m for all T. However, we can show that the sequence a^k mod m for k = 0, 1, 2, ... enters a cycle after some point, that is there exist x and C with x, C < m, such that a^y = a^(y + C) for all y >= x.

Example: For a = 2, m = 12, we get the sequence 2^0, 2^1, ... = 1, 2, 4, 8, 4, 8, ... (mod 12). We can see the cycle with parameters x = 2 and C = 2.

We can find the cycle length via brute-force, by computing the sequence elements a^0, a^1, ... until we find two indices X < Y with a^X = a^Y. Now we set x = X and C = Y - X. This gives us an algorithm with O(m) exponentiations per recursion.

What if we want to do better? Thanks to Jyrki Lahtonen from Math Exchange for providing the essentials for the following algorithm!

Let's evaluate the sequence d_k = gcd(a^k, m) until we find an x with d_x = d_{x+1}. This will take at most log(m) GCD computations, because x is bounded by the highest exponent in the prime factorization of m. Let C = phi(m / d_x). We can now prove that a^{k + C} = a^k for all k >= x, so we have found the cycle parameters in O(m^(1/2)) time.

Let's assume we have found x and C and want to compute a^T mod m now. If T < x, the task is trivial to perform with simple modular exponentiation. Otherwise, we have T >= x and can thus make use of the cycle:

  a^T                                     (mod m)
= a^(x + ((T - x) mod C))                 (mod m)
= a^(x + (-x mod C) + (T mod C) + k*C)    (mod m)     (for some k)
= a^(x + (-x mod C) + k*C) * a^(T mod C)  (mod m)
= a^(x + (-x mod C)) * a^(T mod C)        (mod m)

Again, we have reduced the problem to a subproblem of the same form ("compute T mod C") and two simple modular exponentiations.

Since the modulus is reduced by at least 1 in every iteration, we get a pretty weak bound of O(P^(1/2) * min (P, n)) for the runtime of this algorithm, where n is the height of the stack. In practice we should get a lot better, since the moduli are expected to decrease exponentially. Of course this argument is a bit hand-wavy, maybe some more mathematically-inclined person can improve on it.

There are a few edge cases to consider that actually make your life a bit easier: you can stop immediately if m = 1 (the result is 0 in this case) or if a is a multiple of m (the result is 0 as well in this case).

EDIT: It can be shown that x = C = phi(m) is valid, so as a quick and dirty solution we can use the formula

a^T = a^(phi(m) + T mod phi(m))  (mod m)

for T >= phi(m) or even T >= log_2(m).