How to compute a^^b mod m?

Question

I have to compute efficiently a^^b mod m for large values of a,b,m<2^32
where ^^ is the tetration operator: 2^^4=2^(2^(2^2))

m is not a prime number and not a power of ten.

Can you help?

Answer 1

To be clear, a^^b is not the same thing as a^b, it is the exponential tower a^(a^(a^...^a)) where there are b copies of a, also known as tetration. Let T(a,b) = a^^b so T(a,1) = a and T(a,b) = a^T(a,b-1).

To compute T(a,b) mod m = a^T(a,b-1) mod m, we want to compute a power of a mod m with an extremely large exponent. What you can use is that modular exponentiation is preperiodic with preperiod length at most the greatest power of a prime in the prime factorization of m, which is at most log_2 m, and the period length divides phi(m), where phi(m) is Euler's totient function. In fact, the period length divides Carmichael's function of m, lambda(m). So,

a^k mod m = a^(k+phi(m)) mod m as long as k>log_2 m.

Be careful that a is not necessarily relatively prime to m (or later, to phi(m), phi(phi(m)), etc.). If it were, you could say that a^k mod m = a^(k mod phi(m)) mod m. However, this is not always true when a and m are not relatively prime. For example, phi(100) = 40, and 2^1 mod 100 = 2, but 2^41 mod 100 = 52. You can reduce large exponents to congruent numbers mod phi(m) that are at least log_2 m, so you can say that 2^10001 mod 100 = 2^41 mod 100 but you can't reduce that to 2^1 mod 100. You could define a mod m [minimum x] or use min + mod(a-min,m) as long as a>min.

If T(a,b-1) > [log_2 m], then

a^T(a,b-1) mod m = a^(T(a,b-1) mod phi(m) [minimum [log_2 m]])

otherwise just calculate a^T(a,b-1) mod m.

Recursively calculate this. You can replace phi(m) with lambda(m).

It doesn't take very long to compute the prime factorization of a number under 2^32 since you can determine the prime factors in at most 2^16 = 65,536 trial divisions. Number-theoretic function like phi and lambda are easily expressed in terms of the prime factorization.

At each step, you will need to be able to calculate modular powers with small exponents.

You end up calculating powers mod phi(m), then powers mod phi(phi(m)), then powers mod phi(phi(phi(m))), etc. It doesn't take that many iterations before the iterated phi function is 1, which means you reduce everything to 0, and you no longer get any change by increasing the height of the tower.

Here is an example, of a type that is included in high school math competitions where the competitors are supposed to rediscover this and execute it by hand. What are the last two digits of 14^^2016?

14^^2016 mod 100 
= 14^T(14,2015) mod 100
= 14^(T(14,2015) mod lambda(100) [minimum 6]) mod 100
= 14^(T(14,2015 mod 20 [minimum 6]) mod 100

T(14,2015) mod 20 
= 14^T(14,2014) mod 20
= 14^(T(14,2014) mod 4 [minimum 4]) mod 20

T(14,2014) mod 4
= 14^T(14,2013) mod 4
= 14^(T(14,2013 mod 2 [minimum 2]) mod 4

T(14,2013) mod 2
= 14^T(14,2012) mod 2
= 14^(T(14,2012 mod 1 [minimum 1]) mod 2
= 14^(1) mod 2
= 14 mod 2
= 0

T(14,2014) mod 4 
= 14^(0 mod 2 [minimum 2]) mod 4
= 14^2 mod 4
= 0

T(14,2015) mod 20
= 14^(0 mod 4 [minimum 4]) mod 20 
= 14^4 mod 20
= 16

T(14,2016) mod 100
= 14^(16 mod 20 [minimum 6]) mod 100
= 14^16 mod 100
= 36

So, 14^14^14^...^14 ends in the digits ...36.

Answer 2

I struggled a bit in understanding the accepted answer, & the other stack-overflow question linked by Kaveh helped.

The full solution uses a variety of ideas:

1. Exponentiation by Squaring.
2. Phi function computation from prime factors.
3. The fact that exponentiation is (pre-)periodic with the (loose) guarantee that values greater than phi(m) will be a part of the cycle.

Adding a full simplified working rust implementation for anyone looking for more clarity:

fn phi(mut x: u64) -> u64 {
    let mut c = 2;
    let mut s = x;
    while x > 1 {
        if x % c == 0 {
            s -= s / c;
        }
        while x % c == 0 {
            x /= c;
        }
        c += 1;
    }
    s
}

fn mod_exp(b: u64, p: u64, m: u64) -> u64 {
    if p == 0 {
        1
    } else if p & 1 == 1 {
        (b * mod_exp(b, p - 1, m)) % m
    } else {
        mod_exp((b * b) % m, p >> 1, m)
    }
}

fn mod_tetr(b: u64, p: u64, m: u64) -> u64 {
    if b % m == 0 {
        0
    } else if p == 1 {
        b
    } else {
        mod_exp(b, phi(m) + mod_tetr(b, p - 1, phi(m)), m)
    }
}