Evaluate modular power towers

Question

Given two numbers n and m, evaluate the infinite power tower:

n^(n+1)^(n+2)^(n+3)^(n+4)^... mod m

Keep in mind that ^ is right-associative. So 2^3^4 = 2^(3^4). Now how can you possibly assign a value to an infinite sequence of right-associative operators?

Define f(n,m,i) as the power tower containing the first i terms of the infinite power tower. Then there is some constant C such that for every i > C, f(n,m,i) = f(n,m,C). So you could say the infinite power tower converges on a certain value. We're interested in that value.

---

Your program must be able to compute n = 2017, m = 10^10 in under 10 seconds on a reasonable modern PC. That is, you should implement an actual algorithm, no bruteforcing.

You may assume that n < 2³⁰ and m < 2⁵⁰ for the numerical limits in your programming language, but your algorithm must theoretically work for any size n, m. However your program must be correct for inputs within these size limits, intermediate value overflows are not excused if the inputs are within these limits.

Examples:

2, 10^15
566088170340352

4, 3^20
4

32, 524287
16

Answer 1

Pyth, 23 bytes

M&tG.^HsgBu-G/GH{PGGhHG

Defines a function g, taking m and n in that order.

Try it online

How it works

M&tG.^HsgBu-G/GH{PGGhHG
M                         def g(G, H):
 &tG                        0 if G == 1, else …
                 PG         prime factors of G
                {           deduplicate that
          u-G/GH   G        reduce that on lambda G,H:G-G/H, starting at G
                              (this gives the Euler totient φ(G))
        gB          hH      bifurcate: two-element list [that, g(that, H + 1)]
       s                    sum
    .^H               G     H^that mod G

Python 2, 109 76 bytes

import sympy
def g(n,m):j=sympy.totient(m);return m-1and pow(n,j+g(n+1,j),m)

Try it online!

Why it works

We use the following generalization of Euler’s theorem.

Lemma. n^2φ(m) ≡ n^φ(m) (mod m) for all n (whether or not n is coprime to m).

Proof: For all prime powers p^k dividing m,

• If p divides n, then because φ(m) ≥ φ(p^k) = p^{k − 1}(p − 1) ≥ 2^{k − 1} ≥ k, we have n^2φ(m) ≡ 0 ≡ n^φ(m) (mod p^k).
• Otherwise, because φ(p^k) divides φ(m), Euler’s theorem gives n^2φ(m) ≡ 1 ≡ n^φ(m) (mod p^k).

Therefore, n^2φ(m) ≡ n^φ(m) (mod m).

Corollary. If k ≥ φ(m), then n^k ≡ n^{φ(m) + (k mod φ(m))} (mod m).

Proof: If k ≥ 2φ(m), the lemma gives n^k = n^2φ(m)n^{k - 2φ(m)} ≡ n^φ(m)n^{k - 2φ(m)} = n^{k - φ(m)} (mod m) and we repeat until the exponent is less than 2φ(m).

Answer 2

Haskell, 156 bytes

(?) takes two Integers and returns an Integer, use as (10^10)?2017 (reversed order compared to OP.)

1?n=0
m?n=n&m$m#2+m#2?(n+1)
1#_=1
n#p|m<-until((<2).gcd p)(`div`p)n=m#(p+1)*1`max`div(n*p-n)(p*m)
(_&_)0=1
(x&m)y|(a,b)<-divMod y 2=mod(x^b*(mod(x*x)m&m)a)m

Try it online! (I put the cases to test in the header this time, since they use exponentiation notation.)

Curiously the slowest test case isn't the one with a speed limit (that's almost instant), but the 524287 ? 32 one, because 524287 is a much larger prime than appears in the factors of the other test cases.

How it works

• (x&m)y is x^y `mod` m, or power mod, using exponentiation by squaring.
• n#p is the Euler totient function of n, assuming n has no smaller prime factors than p.
  • m is n with all p factors divided out.
  • If there are k such factors, the totient should itself get a corresponding factor (p-1)*p^(k-1), which is calculated as div(n*p-n)(p*m).
  • 1`max`... handles the case where n wasn't actually divisible by p, which makes the other argument of max equal to 0.
• The main function m?n uses that when y is large enough, n^y `mod` m is the same as n^(t+(y`mod`t)) `mod` m, when t is the totient of m. (The t+ is needed for those prime factors n and m have in common, which all get maximized.)
• The algorithm stops because iterated totient functions eventually hit 1.

Answer 3

Mathematica, 55 bytes

n_~f~1=0;n_~f~m_:=PowerMod[n,(t=EulerPhi@m)+f[n+1,t],m]

Examples:

In[1]:= n_~f~1=0;n_~f~m_:=PowerMod[n,(t=EulerPhi@m)+f[n+1,t],m]

In[2]:= f[2, 10^15]

Out[2]= 566088170340352

In[3]:= f[4, 3^20]

Out[3]= 4

In[4]:= f[32, 524287]

Out[4]= 16

In[5]:= f[2017, 10^10]

Out[5]= 7395978241

Answer 4

Pari/GP, 59 bytes

f(n,m)=if(m==1,0,lift(Mod(n,m)^((t=eulerphi(m))+f(n+1,t))))

Try it online!