Problem 809: Rational Recurrence Relation
View on Project EulerProject Euler Problem 809 Solution
EulerSolve provides an optimized solution for Project Euler Problem 809, Rational Recurrence Relation, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The implementations ultimately evaluate the stabilized residue of the tower \(2, 2^2, 2^{2^2}, \dots\) modulo \(10^{15}\), and then report that residue shifted by \(3\): $$T(10^{15}) - 3 \pmod{10^{15}}.$$ Here \(T(m)\) denotes the eventual residue of a sufficiently tall tower of 2s modulo \(m\). The main difficulty is that powers of 2 behave very differently on the odd part of the modulus and on its power-of-two part, so the computation separates those two pieces and handles them in different ways. Mathematical Approach Define the finite tower by \(a_1 = 2\) and \(a_{n+1} = 2^{a_n}\). For each modulus \(m\), the residues \(a_n \bmod m\) stabilize for large \(n\); call the stable residue \(T(m)\). The implementations compute \(T(m)\) by descending through the odd part of the modulus via Euler's totient function. Step 1: Turn the tower into a modular fixed point Once the residues stabilize, one more exponentiation produces the same class again, so the limit satisfies $$T(m) \equiv 2^{T(m)} \pmod{m}.$$ This identity is not used as a brute-force equation solver. Its real purpose is to show that a very tall tower can be reconstructed from the behavior of its exponent modulo a smaller modulus....
Detailed mathematical approach
Problem Summary
The implementations ultimately evaluate the stabilized residue of the tower \(2, 2^2, 2^{2^2}, \dots\) modulo \(10^{15}\), and then report that residue shifted by \(3\):
$$T(10^{15}) - 3 \pmod{10^{15}}.$$
Here \(T(m)\) denotes the eventual residue of a sufficiently tall tower of 2s modulo \(m\). The main difficulty is that powers of 2 behave very differently on the odd part of the modulus and on its power-of-two part, so the computation separates those two pieces and handles them in different ways.
Mathematical Approach
Define the finite tower by \(a_1 = 2\) and \(a_{n+1} = 2^{a_n}\). For each modulus \(m\), the residues \(a_n \bmod m\) stabilize for large \(n\); call the stable residue \(T(m)\). The implementations compute \(T(m)\) by descending through the odd part of the modulus via Euler's totient function.
Step 1: Turn the tower into a modular fixed point
Once the residues stabilize, one more exponentiation produces the same class again, so the limit satisfies
$$T(m) \equiv 2^{T(m)} \pmod{m}.$$
This identity is not used as a brute-force equation solver. Its real purpose is to show that a very tall tower can be reconstructed from the behavior of its exponent modulo a smaller modulus.
Step 2: Split the modulus into a power of 2 and an odd part
Write
$$m = 2^s q, \qquad q \text{ odd}.$$
For a sufficiently tall tower, the outermost value is \(2^E\) with an exponent \(E\) much larger than \(s\), so the residue is automatically divisible by \(2^s\). Therefore
$$T(m) \equiv 0 \pmod{2^s}.$$
That means the residue can be written as
$$T(m) \equiv 2^s u \pmod{m}$$
for some class \(u \pmod{q}\). The problem is now reduced to finding this odd-part factor \(u\).
Step 3: Reduce the exponent modulo \(\varphi(q)\)
Because \(q\) is odd, \(\gcd(2,q)=1\), so Euler's theorem applies:
$$2^{k+\varphi(q)} \equiv 2^k \pmod{q}.$$
Thus the exponent only matters modulo \(\varphi(q)\). But that exponent is again a very tall tower of 2s, so its stable residue modulo \(\varphi(q)\) is exactly \(T(\varphi(q))\). Hence, modulo \(q\),
$$T(m) \equiv 2^{T(\varphi(q))} \pmod{q}.$$
Substituting \(T(m) \equiv 2^s u\) gives
$$2^s u \equiv 2^{T(\varphi(q))} \pmod{q}.$$
Step 4: Isolate the odd part and allow negative exponents
Since 2 is invertible modulo the odd number \(q\), we may divide by \(2^s\):
$$u \equiv 2^{T(\varphi(q)) - s} \pmod{q}.$$
If \(T(\varphi(q)) \ge s\), this is an ordinary modular power. If \(T(\varphi(q)) \lt s\), the exponent becomes negative, so we interpret it as
$$2^{-r} \equiv (2^{-1})^r \pmod{q}.$$
This is why the implementations explicitly compute the modular inverse of 2 on the odd modulus whenever the corrected exponent is negative.
Step 5: Reassemble the residue modulo \(m\)
Once \(u\) is known, the stabilized tower value is simply
$$T(m) \equiv 2^s u \pmod{m}.$$
For pure powers of two we have \(q=1\), so the answer is \(T(2^s)=0\). Otherwise the recursion continues with the strictly smaller modulus \(\varphi(q)\), which quickly reaches the base case
$$T(1)=0.$$
Worked Example: \(m = 40\)
Take \(m = 40 = 2^3 \cdot 5\). Then \(s=3\), \(q=5\), and
$$\varphi(q) = \varphi(5) = 4.$$
Because \(4\) is a power of two, a tall enough tower is divisible by \(4\), so
$$T(4)=0.$$
Therefore
$$u \equiv 2^{T(4)-3} = 2^{-3} \pmod{5}.$$
The inverse of \(2\) modulo \(5\) is \(3\), so
$$u \equiv 3^3 \equiv 27 \equiv 2 \pmod{5}.$$
Hence
$$T(40) \equiv 2^3 \cdot 2 = 16 \pmod{40}.$$
A direct check confirms the fixed-point relation:
$$2^{16} = 65536 \equiv 16 \pmod{40}.$$
How the Code Works
The C++, Python, and Java implementations all follow the same recursion. They memoize the stabilized residue for each modulus appearing on the totient chain, starting from the base value \(T(1)=0\). For the current modulus they remove all factors of 2, compute the odd part \(q\), factor \(q\) to obtain \(\varphi(q)\), recurse to find \(T(\varphi(q))\), and then correct the exponent by subtracting the number of removed factors of 2.
The odd-part contribution is evaluated as a modular power of 2. When the corrected exponent is negative, the implementation first finds the modular inverse of 2 modulo \(q\) with the extended Euclidean algorithm and then raises that inverse to the needed positive power. Finally it restores the factor \(2^s\), reduces modulo the original modulus, and after the top-level call at \(10^{15}\) it subtracts \(3\) modulo \(10^{15}\).
Complexity Analysis
Let \(m_0 = m\), and at level \(i\) write \(m_i = 2^{s_i} q_i\) with \(q_i\) odd. The recursion then moves to \(m_{i+1} = \varphi(q_i)\) until it reaches 1. The depth is the length of this chain. At each level, the direct implementation computes \(\varphi(q_i)\) by trial division up to \(\sqrt{q_i}\), and performs one modular exponentiation costing \(O(\log q_i)\) modular multiplications. Therefore the running time is
$$O\left(\sum_i \sqrt{q_i} + \sum_i \log q_i\right),$$
while the memo table stores one residue per chain value, so the memory use is \(O(k)\), where \(k\) is the chain length. For the actual target modulus the chain shrinks very quickly, so the practical runtime is small.
Footnotes and References
- Problem page: https://projecteuler.net/problem=809
- Euler's totient function: Wikipedia - Euler's totient function
- Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse
- Tetration and infinite power towers: Wikipedia - Tetration
- Chinese remainder theorem: Wikipedia - Chinese remainder theorem
Problem 809 source code
C++
#include <cstdint>
#include <iostream>
#include <unordered_map>
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
static u64 pow_mod(u64 a, u64 e, u64 mod) {
if (mod == 1) return 0;
u64 r = 1 % mod;
a %= mod;
while (e > 0) {
if (e & 1ULL) r = static_cast<u64>((static_cast<u128>(r) * a) % mod);
a = static_cast<u64>((static_cast<u128>(a) * a) % mod);
e >>= 1ULL;
}
return r;
}
static u64 inv_mod(u64 a, u64 mod) {
i64 t = 0, nt = 1;
i64 r = static_cast<i64>(mod), nr = static_cast<i64>(a % mod);
while (nr != 0) {
const i64 q = r / nr;
const i64 tt = t - q * nt;
t = nt;
nt = tt;
const i64 rr = r - q * nr;
r = nr;
nr = rr;
}
if (t < 0) t += static_cast<i64>(mod);
return static_cast<u64>(t);
}
static u64 pow_mod_signed_2(i64 exp, u64 mod) {
if (mod == 1) return 0;
if (exp >= 0) return pow_mod(2, static_cast<u64>(exp), mod);
const u64 inv2 = inv_mod(2, mod);
return pow_mod(inv2, static_cast<u64>(-exp), mod);
}
static u64 phi(u64 n) {
if (n == 0) return 0;
u64 result = n;
if ((n & 1ULL) == 0) {
result -= result / 2;
while ((n & 1ULL) == 0) n >>= 1ULL;
}
for (u64 p = 3; p * p <= n; p += 2) {
if (n % p != 0) continue;
result -= result / p;
while (n % p == 0) n /= p;
}
if (n > 1) result -= result / n;
return result;
}
static u64 tower2_fixpoint_mod(u64 mod, std::unordered_map<u64, u64>& memo) {
auto it = memo.find(mod);
if (it != memo.end()) return it->second;
u64 odd = mod;
i64 twos = 0;
while ((odd & 1ULL) == 0) {
odd >>= 1ULL;
++twos;
}
const u64 t = phi(odd);
const i64 w = (t == 1) ? 0 : static_cast<i64>(tower2_fixpoint_mod(t, memo)) - twos;
const u64 odd_part = pow_mod_signed_2(w, odd);
const u64 ans = (odd_part << twos) % mod;
memo.emplace(mod, ans);
return ans;
}
int main() {
const u64 MOD = 1'000'000'000'000'000ULL;
std::unordered_map<u64, u64> memo;
memo.reserve(64);
memo.emplace(1, 0);
const u64 top = tower2_fixpoint_mod(MOD, memo);
std::cout << (top + MOD - 3) % MOD << '\n';
return 0;
}
Python
import sys
sys.setrecursionlimit(2000)
def inv_mod(a, mod):
t, nt = 0, 1
r, nr = mod, a % mod
while nr != 0:
q = r // nr
t, nt = nt, t - q * nt
r, nr = nr, r - q * nr
if t < 0:
t += mod
return t
def pow_mod_signed_2(exp, mod):
if mod == 1:
return 0
if exp >= 0:
return pow(2, exp, mod)
inv2 = inv_mod(2, mod)
return pow(inv2, -exp, mod)
def phi(n):
if n == 0:
return 0
result = n
if n % 2 == 0:
result -= result // 2
while n % 2 == 0:
n //= 2
p = 3
while p * p <= n:
if n % p == 0:
result -= result // p
while n % p == 0:
n //= p
p += 2
if n > 1:
result -= result // n
return result
memo = {1: 0}
def tower2_fixpoint_mod(mod):
if mod in memo:
return memo[mod]
odd = mod
twos = 0
while odd % 2 == 0:
odd //= 2
twos += 1
t = phi(odd)
w = 0 if t == 1 else tower2_fixpoint_mod(t) - twos
odd_part = pow_mod_signed_2(w, odd)
ans = (odd_part << twos) % mod
memo[mod] = ans
return ans
def solve():
MOD = 1000000000000000
top = tower2_fixpoint_mod(MOD)
return str((top + MOD - 3) % MOD)
if __name__ == "__main__":
print(solve())
Java
import java.util.HashMap;
public class Euler809 {
static long powMod(long a, long e, long mod) {
if (mod == 1)
return 0;
long r = 1 % mod;
a %= mod;
while (e > 0) {
if ((e & 1L) == 1L) {
r = java.math.BigInteger.valueOf(r).multiply(java.math.BigInteger.valueOf(a))
.mod(java.math.BigInteger.valueOf(mod)).longValue();
}
a = java.math.BigInteger.valueOf(a).multiply(java.math.BigInteger.valueOf(a))
.mod(java.math.BigInteger.valueOf(mod)).longValue();
e >>= 1L;
}
return r;
}
static long invMod(long a, long mod) {
long t = 0, nt = 1;
long r = mod, nr = a % mod;
while (nr != 0) {
long q = r / nr;
long tt = t - q * nt;
t = nt;
nt = tt;
long rr = r - q * nr;
r = nr;
nr = rr;
}
if (t < 0)
t += mod;
return t;
}
static long powModSigned2(long exp, long mod) {
if (mod == 1)
return 0;
if (exp >= 0)
return powMod(2, exp, mod);
long inv2 = invMod(2, mod);
return powMod(inv2, -exp, mod);
}
static long phi(long n) {
if (n == 0)
return 0;
long result = n;
if ((n & 1L) == 0L) {
result -= result / 2;
while ((n & 1L) == 0L) {
n >>= 1L;
}
}
for (long p = 3; p * p <= n; p += 2) {
if (n % p != 0)
continue;
result -= result / p;
while (n % p == 0) {
n /= p;
}
}
if (n > 1) {
result -= result / n;
}
return result;
}
static HashMap<Long, Long> memo = new HashMap<>();
static long tower2FixpointMod(long mod) {
Long cached = memo.get(mod);
if (cached != null)
return cached;
long odd = mod;
long twos = 0;
while ((odd & 1L) == 0L) {
odd >>= 1L;
twos++;
}
long t = phi(odd);
long w = (t == 1) ? 0 : tower2FixpointMod(t) - twos;
long oddPart = powModSigned2(w, odd);
long ans = (oddPart << twos) % mod;
memo.put(mod, ans);
return ans;
}
public static String solve() {
memo.put(1L, 0L);
long MOD = 1000000000000000L;
long top = tower2FixpointMod(MOD);
return Long.toString((top + MOD - 3) % MOD);
}
public static void main(String[] args) {
System.out.println(solve());
}
}