# Evaluate modular power towers

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 < 230 and m < 250 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

• Tip (for contenders): n and m are not guaranteed to be co-prime. – Leaky Nun Jun 3 '17 at 16:28
• 10^10 (and 10^20, and potentially 3^20 for signed integers) is larger than many languages' default integer types. Is it required that this large of an input be supported? – Doorknob Jun 3 '17 at 16:46
• @orlp Does that "yes" include 10^20? Because that doesn't fit in a 64-bit integer, so if you want to require it, I'd suggest pointing it out explicitly, because otherwise you're going to get a lot of invalid answers by people just assuming that 64 bit integers are going to be accurate enough. – Martin Ender Jun 3 '17 at 16:59
• Either way, what is the largest input we need to support? – Martin Ender Jun 3 '17 at 16:59
• @Doorknob I added more lenient limits to the challenge. However your algorithm must theoretically work for any size m, n. – orlp Jun 3 '17 at 17:18

# 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. n2φ(m)nφ(m) (mod m) for all n (whether or not n is coprime to m).

Proof: For all prime powers pk dividing m,

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

Therefore, n2φ(m)nφ(m) (mod m).

Corollary. If k ≥ φ(m), then nknφ(m) + (k mod φ(m)) (mod m).

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

• How does this handle the case where the base and modulo isn't coprime? P.S. sympy has totient function. – orlp Jun 4 '17 at 12:51
• @orlp I’ve added a proof. Not sure how I missed sympy.totient. – Anders Kaseorg Jun 4 '17 at 22:35
• I see now. Nice method! – orlp Jun 5 '17 at 4:20

(?) 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)(divp)n=m#(p+1)*1maxdiv(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).
• 1max... 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+(ymodt)) 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.

# 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


# 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!