Problem 396: Weak Goodstein Sequence

Project Euler Problem 396 Solution

EulerSolve provides an optimized solution for Project Euler Problem 396, Weak Goodstein Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For the weak Goodstein sequence starting from $n$, write the current value in base $b$, keep the same digit string, reinterpret it in base $b+1$, subtract $1$, and then increase the base. If $g_0(n)=n$, then $G(n)$ is the least step count for which the iterated process reaches $0$. The program computes $$\sum_{n=1}^{15} G(n) \pmod{10^9}.$$ Mathematical Approach Step 1: One Weak Goodstein Step If the current base is $b$ and $$x=\sum_{i=0}^{k} d_i b^i,\qquad 0\le d_i \lt b,$$ then one weak Goodstein step replaces $x$ by $$T_b(x)=\sum_{i=0}^{k} d_i (b+1)^i - 1.$$ So the sequence is $$g_{t+1}(n)=T_{t+2}(g_t(n)),\qquad g_0(n)=n,$$ because the starting base is $2$, then $3$, then $4$, and so on. Step 2: Exact Small Values Used as Anchors The implementations compute $G(n)$ exactly for $0 \le n \lt 8$ by direct simulation. The checkpoint values are $$G(0)=0,\quad G(1)=1,\quad G(2)=3,\quad G(3)=5,\quad G(4)=21,\quad G(5)=61,\quad G(6)=381,\quad G(7)=2045.$$ In particular, the code verifies $$\sum_{n=1}^{7} G(n)=2517.$$ These exact values are the only non-modular inputs needed for the larger cases $8 \le n \le 15$....

Detailed mathematical approach

Problem Summary

For the weak Goodstein sequence starting from $n$, write the current value in base $b$, keep the same digit string, reinterpret it in base $b+1$, subtract $1$, and then increase the base. If $g_0(n)=n$, then $G(n)$ is the least step count for which the iterated process reaches $0$. The program computes

$$\sum_{n=1}^{15} G(n) \pmod{10^9}.$$

Mathematical Approach

Step 1: One Weak Goodstein Step

If the current base is $b$ and

$$x=\sum_{i=0}^{k} d_i b^i,\qquad 0\le d_i \lt b,$$

then one weak Goodstein step replaces $x$ by

$$T_b(x)=\sum_{i=0}^{k} d_i (b+1)^i - 1.$$

So the sequence is

$$g_{t+1}(n)=T_{t+2}(g_t(n)),\qquad g_0(n)=n,$$

because the starting base is $2$, then $3$, then $4$, and so on.

Step 2: Exact Small Values Used as Anchors

The implementations compute $G(n)$ exactly for $0 \le n \lt 8$ by direct simulation. The checkpoint values are

$$G(0)=0,\quad G(1)=1,\quad G(2)=3,\quad G(3)=5,\quad G(4)=21,\quad G(5)=61,\quad G(6)=381,\quad G(7)=2045.$$

In particular, the code verifies

$$\sum_{n=1}^{7} G(n)=2517.$$

These exact values are the only non-modular inputs needed for the larger cases $8 \le n \le 15$.

Step 3: Reducing $n$ to a Pure Cube State

For $8 \le n \le 15$, write

$$n=8+r=2^3+r,\qquad 0\le r\lt 8.$$

Let

$$t=G(r).$$

The lower three binary positions evolve exactly like the weak Goodstein sequence for $r$, while the leading binary digit at position $3$ stays above that lower part. After those $t$ steps, the tail has vanished, the current base is $t+2$, and the state is the pure cube

$$g_t(n)=(t+2)^3.$$

So the whole problem reduces to counting how many steps remain from a number of the form $b^3$ when the current base is $b$, with

$$b=t+2.$$

Step 4: The Square-Block Lemma

Now consider a state whose current base is $b$ and whose value is $a b^2$, so its base-$b$ digits are $a,0,0$. The local digit dynamics produce a clean block length:

$$Q(b)=(b+1)\left(2^{b+1}-1\right).$$

A direct digit chase shows that one such square block lowers the leading coefficient by $1$ and replaces the base by

$$b'=(b+1)2^{b+1}-1.$$

In other words, from the state $a b^2$ in base $b$, after exactly $Q(b)$ further steps the process reaches

$$ (a-1)(b')^2 $$

in base $b'$. For $a=1$ this means termination; for $a=2$ it means one square layer has been removed.

The reason for the closed form is that the lower two digits behave like a carry chain whose countdown lengths are

$$ (b+1),\ 2(b+1),\ 4(b+1),\ \dots,\ 2^b(b+1), $$

and summing this geometric progression gives $Q(b)$.

Step 5: From $b^3$ to the Recurrence Used in Code

Starting from $b^3$ in base $b$, the first square block converts the cube into a square with coefficient $b$: after $Q(b)$ steps the process reaches

$$b\,(b_1)^2,\qquad b_1=(b+1)2^{b+1}-1.$$

Each additional square block removes one more unit from that leading coefficient. It is convenient to write

$$s_0=b+1,\qquad s_{i+1}=s_i 2^{s_i},\qquad b_i=s_i-1.$$

Then

$$Q(b_i)=s_i\left(2^{s_i}-1\right),$$

and the total number of remaining steps from $b^3$ is

$$F_3(b)=\sum_{i=0}^{b} s_i\left(2^{s_i}-1\right).$$

Returning to $n=8+r$, with $t=G(r)$ and $b=t+2$, the program uses

$$H(s_0,m)=\sum_{i=0}^{m} s_i\left(2^{s_i}-1\right),\qquad s_{i+1}=s_i 2^{s_i},$$

where

$$s_0=t+3,\qquad m=t+2.$$

Therefore

$$G(n)\equiv t + H(t+3,t+2)\pmod{10^9},\qquad n\ge 8,$$

which is exactly the formula implemented by compute_g_mod.

Step 6: Why Modular Evaluation Is Necessary

The recurrence $s_{i+1}=s_i 2^{s_i}$ explodes immediately. The solver never forms these integers explicitly after the small exact phase. Instead it computes everything modulo

$$10^9=2^9\cdot 5^9.$$

For a modulus of the form $2^a 5^b$, the $2$-power part of $2^s$ is simple: if $s\ge a$, then

$$2^s\equiv 0 \pmod{2^a}.$$

For the $5$-power part, $2$ is coprime to $5$, so Euler's theorem gives period

$$\varphi(5^b)=4\cdot 5^{b-1}.$$

That is why the code builds the modulus chain

$$10^9,\ 4\cdot 5^8,\ 4\cdot 5^7,\ \dots,\ 20,\ 4,$$

and stores the current $s_i$ modulo every entry in that chain. To recover $2^{s_i}\bmod 2^a5^b$, it computes the residue modulo $5^b$ from $s_i \bmod 4\cdot 5^{b-1}$, computes the residue modulo $2^a$ from the exact small value when needed, and combines the two pieces with the Chinese remainder theorem. This is the job of pow2_mod_composite.

Worked Example: $n=8$

Here $r=0$, so $t=G(0)=0$. Hence

$$b=t+2=2,\qquad s_0=t+3=3,\qquad m=t+2=2.$$

The first block length is

$$s_0(2^{s_0}-1)=3(8-1)=21.$$

After these $21$ steps, the sequence has become $2\cdot 23^2$ in base $23$.

Next,

$$s_1=s_0 2^{s_0}=24,$$

so the second block length is

$$s_1(2^{s_1}-1)=24(2^{24}-1)=402653160.$$

After this block, the state is $402653183^2$ in base $402653183$.

Finally,

$$s_2=s_1 2^{s_1}=402653184,$$

and the last block $s_2(2^{s_2}-1)$ is astronomically large. The modular engine evaluates the total directly and the implementations verify the checkpoint

$$H(3,2)\equiv 722374141 \pmod{10^9}.$$

How the Code Works

All three solution files follow the same structure. goodstein_small simulates the process exactly for $n \lt 8$. compute_g_mod uses those anchor values, applies the reduction $n=8+r$, and then calls compute_h_mod.

compute_h_mod iterates through the block sum $H(s_0,m)$ while keeping only residues for the modulus chain. build_mod_chain precomputes the CRT data, and pow2_mod_composite is the core routine that reconstructs $2^{s_i}$ modulo each composite modulus. The final solve function sums $G(n)$ for $1\le n\le 15$ modulo $10^9$.

Complexity Analysis

The exact simulation phase is tiny because it only evaluates $G(n)$ for $n=1,\dots,7$, and the largest checkpoint is $G(7)=2045$. In the modular phase, each call to compute_h_mod performs $m+1$ updates across a fixed chain of length $10$, so its cost is $O(mL)$ with $L=10$, and its memory use is $O(L)$.

For this problem, $m=t+2$ and $t\in\{0,G(1),\dots,G(7)\}$, so even the largest case has $m=2047$. In practice the whole computation is constant-scale and the modular arithmetic dominates the runtime.

Footnotes and References

Problem page: https://projecteuler.net/problem=396
Goodstein sequences and termination: Wikipedia - Goodstein's theorem
Chinese remainder theorem: Wikipedia - Chinese remainder theorem
Euler's theorem and $\varphi(5^k)=4\cdot 5^{k-1}$: Wikipedia - Euler's theorem

Problem 396 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using BigInt = boost::multiprecision::cpp_int;

constexpr i64 kMod = 1000000000LL;

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

i64 mod_pow(i64 base, i64 exp, i64 mod) {
    i64 result = 1 % mod;
    i64 cur = base % mod;
    i64 e = exp;
    while (e > 0) {
        if (e & 1LL) {
            result = static_cast<i64>((__int128)result * cur % mod);
        }
        cur = static_cast<i64>((__int128)cur * cur % mod);
        e >>= 1LL;
    }
    return result;
}

i64 egcd(i64 a, i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 x1 = 0;
    i64 y1 = 0;
    const i64 g = egcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

i64 mod_inv(i64 a, i64 mod) {
    i64 x = 0;
    i64 y = 0;
    const i64 g = egcd(a, mod, x, y);
    if (g != 1) {
        throw std::runtime_error("Inverse does not exist");
    }
    i64 res = x % mod;
    if (res < 0) {
        res += mod;
    }
    return res;
}

struct ModEntry {
    i64 mod;
    int pow2;
    int pow5;
    i64 period_mod;   // modulus needed for exponent (0 => not needed)
    i64 mod2_part;    // 2^pow2
    i64 mod5_part;    // 5^pow5
    i64 inv_mod2_in5; // inverse of mod2_part modulo mod5_part (if both parts exist)
};

std::vector<ModEntry> build_mod_chain() {
    std::vector<ModEntry> entries;

    auto pow5 = [](int k) -> i64 {
        i64 r = 1;
        for (int i = 0; i < k; ++i) {
            r *= 5;
        }
        return r;
    };

    // Main modulus 10^9 = 2^9 * 5^9, exponent period for odd part: 4*5^8.
    entries.push_back(ModEntry{kMod, 9, 9, 4 * pow5(8), 1LL << 9, pow5(9), 0});

    // Chain: 4*5^8, 4*5^7, ..., 4*5^1, 4.
    for (int e5 = 8; e5 >= 0; --e5) {
        const i64 mod = 4 * pow5(e5);
        const i64 period = (e5 == 0 ? 0 : 4 * pow5(e5 - 1));
        entries.push_back(ModEntry{mod, 2, e5, period, 4, pow5(e5), 0});
    }

    for (ModEntry& entry : entries) {
        if (entry.pow2 > 0 && entry.pow5 > 0) {
            entry.inv_mod2_in5 = mod_inv(entry.mod2_part % entry.mod5_part, entry.mod5_part);
        } else {
            entry.inv_mod2_in5 = 0;
        }
    }

    return entries;
}

i64 pow2_mod_composite(
    const ModEntry& entry,
    const std::vector<i64>& residues,
    const std::vector<ModEntry>& chain,
    const i64 exact_s,
    const bool exact_known) {
    if (entry.mod == 1) {
        return 0;
    }
    if (entry.pow5 == 0) {
        if (entry.pow2 == 0) {
            return 0;
        }
        if (exact_known && exact_s < entry.pow2) {
            return (1LL << exact_s) % entry.mod;
        }
        return 0;
    }

    i64 exp_mod = 0;
    if (entry.period_mod > 0) {
        int idx = -1;
        for (int i = 0; i < static_cast<int>(chain.size()); ++i) {
            if (chain[static_cast<std::size_t>(i)].mod == entry.period_mod) {
                idx = i;
                break;
            }
        }
        if (idx < 0) {
            throw std::runtime_error("Period modulus not in chain");
        }
        exp_mod = residues[static_cast<std::size_t>(idx)];
    }

    const i64 p5 = mod_pow(2, exp_mod, entry.mod5_part);

    if (entry.pow2 == 0) {
        return p5;
    }

    i64 p2 = 0;
    if (exact_known && exact_s < entry.pow2) {
        p2 = (1LL << exact_s) % entry.mod2_part;
    } else {
        p2 = 0;
    }

    // CRT for x ≡ p2 (mod 2^a), x ≡ p5 (mod 5^b).
    const i64 diff = (p5 - p2) % entry.mod5_part;
    const i64 adjusted = (diff + entry.mod5_part) % entry.mod5_part;
    const i64 t = static_cast<i64>((__int128)adjusted * entry.inv_mod2_in5 % entry.mod5_part);
    const i64 x = p2 + static_cast<i64>((__int128)entry.mod2_part * t % entry.mod);
    return x % entry.mod;
}

i64 compute_h_mod(i64 s0, int m) {
    const std::vector<ModEntry> chain = build_mod_chain();
    std::vector<i64> residues(chain.size(), 0);
    for (std::size_t i = 0; i < chain.size(); ++i) {
        residues[i] = s0 % chain[i].mod;
    }

    i64 exact_s = s0;
    bool exact_known = true;
    constexpr i64 kExactLimit = 1000;

    i64 sum_mod = 0;
    for (int i = 0; i <= m; ++i) {
        const i64 p_main = pow2_mod_composite(chain[0], residues, chain, exact_s, exact_known);
        const i64 term = static_cast<i64>((__int128)residues[0] * ((p_main + kMod - 1) % kMod) % kMod);
        sum_mod += term;
        if (sum_mod >= kMod) {
            sum_mod -= kMod;
        }

        if (i == m) {
            break;
        }

        std::vector<i64> next_residues(chain.size(), 0);
        for (std::size_t j = 0; j < chain.size(); ++j) {
            const i64 p = pow2_mod_composite(chain[j], residues, chain, exact_s, exact_known);
            next_residues[j] =
                static_cast<i64>((__int128)residues[j] * p % chain[j].mod);
        }
        residues.swap(next_residues);

        if (exact_known) {
            if (exact_s > 20) {
                exact_known = false;
            } else {
                i64 p2 = 1;
                for (int k = 0; k < exact_s; ++k) {
                    if (p2 > kExactLimit) {
                        break;
                    }
                    p2 <<= 1;
                }
                if (p2 > kExactLimit || exact_s > kExactLimit / std::max<i64>(1, p2)) {
                    exact_known = false;
                } else {
                    exact_s *= p2;
                    if (exact_s > kExactLimit) {
                        exact_known = false;
                    }
                }
            }
        }
    }

    return sum_mod;
}

i64 goodstein_small(const int n) {
    i64 g = n;
    int base = 2;
    i64 count = 0;
    while (g > 0) {
        ++count;

        std::vector<int> digits;
        i64 x = g;
        while (x > 0) {
            digits.push_back(static_cast<int>(x % base));
            x /= base;
        }
        if (digits.empty()) {
            digits.push_back(0);
        }

        i64 next = 0;
        for (int i = static_cast<int>(digits.size()) - 1; i >= 0; --i) {
            next = next * static_cast<i64>(base + 1) + digits[static_cast<std::size_t>(i)];
        }
        g = next - 1;
        ++base;
    }
    return count;
}

i64 compute_g_mod(const int n, const std::array<i64, 8>& g_small) {
    if (n < 8) {
        return g_small[static_cast<std::size_t>(n)] % kMod;
    }
    const int low = n - 8;
    const i64 t = g_small[static_cast<std::size_t>(low)];
    const i64 s = t + 3;
    const int m = static_cast<int>(t + 2);
    const i64 h = compute_h_mod(s, m);
    return (t + h) % kMod;
}

i64 brute_h_mod_small(i64 s0, const int m) {
    BigInt s = s0;
    BigInt total = 0;
    const BigInt MOD = kMod;
    for (int i = 0; i <= m; ++i) {
        BigInt p = BigInt(1) << static_cast<unsigned>(s.convert_to<unsigned>());
        total += s * (p - 1);
        if (i < m) {
            s *= p;
        }
    }
    return static_cast<i64>((total % MOD).convert_to<i64>());
}

bool run_checkpoints() {
    if (goodstein_small(2) != 3) {
        std::cerr << "Checkpoint failed: G(2)\n";
        return false;
    }
    if (goodstein_small(4) != 21) {
        std::cerr << "Checkpoint failed: G(4)\n";
        return false;
    }
    if (goodstein_small(6) != 381) {
        std::cerr << "Checkpoint failed: G(6)\n";
        return false;
    }

    i64 sum_upto_7 = 0;
    for (int n = 1; n < 8; ++n) {
        sum_upto_7 += goodstein_small(n);
    }
    if (sum_upto_7 != 2517) {
        std::cerr << "Checkpoint failed: sum_{1<=n<8} G(n)\n";
        return false;
    }

    // Modular chain check on a still-manageable case with exact brute force.
    const i64 fast_h = compute_h_mod(2, 2);
    const i64 brute_h = brute_h_mod_small(2, 2);
    if (fast_h != brute_h) {
        std::cerr << "Checkpoint failed: H(s,m) modular engine\n";
        return false;
    }

    // Extra nontrivial modular checkpoint: H(3,2) can still be verified exactly.
    if (compute_h_mod(3, 2) != 722374141LL) {
        std::cerr << "Checkpoint failed: H(3,2)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::array<i64, 8> g_small{};
    g_small[0] = 0;
    for (int n = 1; n < 8; ++n) {
        g_small[static_cast<std::size_t>(n)] = goodstein_small(n);
    }

    i64 answer = 0;
    for (int n = 1; n < 16; ++n) {
        answer += compute_g_mod(n, g_small);
        answer %= kMod;
    }

    std::cout << answer << '\n';
    return 0;
}

Python

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def egcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = egcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return g, x, y

def mod_inv(a, mod):
    g, x, y = egcd(a, mod)
    if g != 1:
        raise Exception("Inverse does not exist")
    return x % mod

class ModEntry:
    def __init__(self, mod, pow2, pow5, period_mod, mod2_part, mod5_part):
        self.mod = mod
        self.pow2 = pow2
        self.pow5 = pow5
        self.period_mod = period_mod
        self.mod2_part = mod2_part
        self.mod5_part = mod5_part
        if pow2 > 0 and pow5 > 0:
            self.inv_mod2_in5 = mod_inv(self.mod2_part % self.mod5_part, self.mod5_part)
        else:
            self.inv_mod2_in5 = 0

def build_mod_chain():
    entries = []
    
    entries.append(ModEntry(10**9, 9, 9, 4 * (5**8), 1 << 9, 5**9))
    
    for e5 in range(8, -1, -1):
        mod = 4 * (5**e5)
        period = 0 if e5 == 0 else 4 * (5**(e5 - 1))
        entries.append(ModEntry(mod, 2, e5, period, 4, 5**e5))
        
    return entries

def pow2_mod_composite(entry, residues, chain, exact_s, exact_known):
    if entry.mod == 1:
        return 0
    if entry.pow5 == 0:
        if entry.pow2 == 0:
            return 0
        if exact_known and exact_s < entry.pow2:
            return (1 << exact_s) % entry.mod
        return 0
        
    exp_mod = 0
    if entry.period_mod > 0:
        idx = next(i for i, c in enumerate(chain) if c.mod == entry.period_mod)
        exp_mod = residues[idx]
        
    p5 = mod_pow(2, exp_mod, entry.mod5_part)
    
    if entry.pow2 == 0:
        return p5
        
    p2 = 0
    if exact_known and exact_s < entry.pow2:
        p2 = (1 << exact_s) % entry.mod2_part
        
    diff = (p5 - p2) % entry.mod5_part
    adjusted = (diff + entry.mod5_part) % entry.mod5_part
    t = (adjusted * entry.inv_mod2_in5) % entry.mod5_part
    x = p2 + entry.mod2_part * t
    return x % entry.mod

def compute_h_mod(s0, m):
    kMod = 1000000000
    chain = build_mod_chain()
    residues = [s0 % c.mod for c in chain]
    
    exact_s = s0
    exact_known = True
    kExactLimit = 1000
    
    sum_mod = 0
    for i in range(m + 1):
        p_main = pow2_mod_composite(chain[0], residues, chain, exact_s, exact_known)
        term = (residues[0] * ((p_main + kMod - 1) % kMod)) % kMod
        sum_mod = (sum_mod + term) % kMod
        
        if i == m:
            break
            
        next_residues = [0] * len(chain)
        for j, c in enumerate(chain):
            p = pow2_mod_composite(c, residues, chain, exact_s, exact_known)
            next_residues[j] = (residues[j] * p) % c.mod
            
        residues = next_residues
        
        if exact_known:
            if exact_s > 20:
                exact_known = False
            else:
                p2 = 1 << exact_s
                if p2 > kExactLimit or exact_s > kExactLimit // max(1, p2):
                    exact_known = False
                else:
                    exact_s *= p2
                    if exact_s > kExactLimit:
                        exact_known = False
                        
    return sum_mod

def goodstein_small(n):
    g = n
    base = 2
    count = 0
    while g > 0:
        count += 1
        digits = []
        x = g
        while x > 0:
            digits.append(x % base)
            x //= base
        if not digits:
            digits.append(0)
            
        next_val = 0
        for d in reversed(digits):
            next_val = next_val * (base + 1) + d
            
        g = next_val - 1
        base += 1
        
    return count

def compute_g_mod(n, g_small):
    kMod = 1000000000
    if n < 8:
        return g_small[n] % kMod
    low = n - 8
    t = g_small[low]
    s = t + 3
    m = int(t + 2)
    h = compute_h_mod(s, m)
    return (t + h) % kMod

def solve():
    g_small = [0] * 8
    for n in range(1, 8):
        g_small[n] = goodstein_small(n)
        
    kMod = 1000000000
    ans = 0
    for n in range(1, 16):
        ans = (ans + compute_g_mod(n, g_small)) % kMod
        
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler396 {
    static final long K_MOD = 1000000000L;

    static long modPow(long base, long exp, long mod) {
        return BigInteger.valueOf(base)
                .modPow(BigInteger.valueOf(exp), BigInteger.valueOf(mod))
                .longValue();
    }

    static long modInv(long a, long mod) {
        return BigInteger.valueOf(a).modInverse(BigInteger.valueOf(mod)).longValue();
    }

    static class ModEntry {
        long mod;
        int pow2;
        int pow5;
        long period_mod;
        long mod2_part;
        long mod5_part;
        long inv_mod2_in5;

        ModEntry(long m, int p2, int p5, long pm, long m2p, long m5p) {
            this.mod = m;
            this.pow2 = p2;
            this.pow5 = p5;
            this.period_mod = pm;
            this.mod2_part = m2p;
            this.mod5_part = m5p;

            if (p2 > 0 && p5 > 0) {
                this.inv_mod2_in5 = modInv(m2p % m5p, m5p);
            } else {
                this.inv_mod2_in5 = 0;
            }
        }
    }

    static long pow5(int k) {
        long r = 1;
        for (int i = 0; i < k; i++)
            r *= 5;
        return r;
    }

    static List<ModEntry> buildModChain() {
        List<ModEntry> entries = new ArrayList<>();
        entries.add(new ModEntry(K_MOD, 9, 9, 4 * pow5(8), 1L << 9, pow5(9)));

        for (int e5 = 8; e5 >= 0; --e5) {
            long mod = 4 * pow5(e5);
            long period = (e5 == 0 ? 0 : 4 * pow5(e5 - 1));
            entries.add(new ModEntry(mod, 2, e5, period, 4, pow5(e5)));
        }
        return entries;
    }

    static long pow2ModComposite(ModEntry entry, long[] residues, List<ModEntry> chain, long exact_s,
            boolean exact_known) {
        if (entry.mod == 1)
            return 0;
        if (entry.pow5 == 0) {
            if (entry.pow2 == 0)
                return 0;
            if (exact_known && exact_s < entry.pow2) {
                return (1L << exact_s) % entry.mod;
            }
            return 0;
        }

        long exp_mod = 0;
        if (entry.period_mod > 0) {
            int idx = -1;
            for (int i = 0; i < chain.size(); i++) {
                if (chain.get(i).mod == entry.period_mod) {
                    idx = i;
                    break;
                }
            }
            exp_mod = residues[idx];
        }

        long p5 = modPow(2, exp_mod, entry.mod5_part);

        if (entry.pow2 == 0)
            return p5;

        long p2 = 0;
        if (exact_known && exact_s < entry.pow2) {
            p2 = (1L << exact_s) % entry.mod2_part;
        }

        long diff = (p5 - p2) % entry.mod5_part;
        long adjusted = (diff + entry.mod5_part) % entry.mod5_part;

        long t = BigInteger.valueOf(adjusted)
                .multiply(BigInteger.valueOf(entry.inv_mod2_in5))
                .mod(BigInteger.valueOf(entry.mod5_part))
                .longValue();

        long x = p2 + BigInteger.valueOf(entry.mod2_part)
                .multiply(BigInteger.valueOf(t))
                .mod(BigInteger.valueOf(entry.mod))
                .longValue();

        return x % entry.mod;
    }

    static long computeHMod(long s0, int m) {
        List<ModEntry> chain = buildModChain();
        long[] residues = new long[chain.size()];
        for (int i = 0; i < chain.size(); i++) {
            residues[i] = s0 % chain.get(i).mod;
        }

        long exact_s = s0;
        boolean exact_known = true;
        long kExactLimit = 1000;

        long sum_mod = 0;
        for (int i = 0; i <= m; i++) {
            long p_main = pow2ModComposite(chain.get(0), residues, chain, exact_s, exact_known);
            long term = BigInteger.valueOf(residues[0])
                    .multiply(BigInteger.valueOf((p_main + K_MOD - 1) % K_MOD))
                    .mod(BigInteger.valueOf(K_MOD))
                    .longValue();

            sum_mod = (sum_mod + term) % K_MOD;

            if (i == m)
                break;

            long[] next_residues = new long[chain.size()];
            for (int j = 0; j < chain.size(); j++) {
                long p = pow2ModComposite(chain.get(j), residues, chain, exact_s, exact_known);
                next_residues[j] = BigInteger.valueOf(residues[j])
                        .multiply(BigInteger.valueOf(p))
                        .mod(BigInteger.valueOf(chain.get(j).mod))
                        .longValue();
            }
            residues = next_residues;

            if (exact_known) {
                if (exact_s > 20) {
                    exact_known = false;
                } else {
                    long p2 = 1L << exact_s;
                    if (p2 > kExactLimit || exact_s > kExactLimit / Math.max(1, p2)) {
                        exact_known = false;
                    } else {
                        exact_s *= p2;
                        if (exact_s > kExactLimit) {
                            exact_known = false;
                        }
                    }
                }
            }
        }
        return sum_mod;
    }

    static long goodsteinSmall(int n) {
        long g = n;
        long base = 2;
        long count = 0;

        while (g > 0) {
            count++;
            List<Integer> digits = new ArrayList<>();
            long x = g;
            while (x > 0) {
                digits.add((int) (x % base));
                x /= base;
            }
            if (digits.isEmpty())
                digits.add(0);

            long nextVal = 0;
            for (int i = digits.size() - 1; i >= 0; i--) {
                nextVal = nextVal * (base + 1) + digits.get(i);
            }
            g = nextVal - 1;
            base++;
        }
        return count;
    }

    static long computeGMod(int n, long[] gSmall) {
        if (n < 8)
            return gSmall[n] % K_MOD;
        int low = n - 8;
        long t = gSmall[low];
        long s = t + 3;
        int m = (int) (t + 2);
        long h = computeHMod(s, m);
        return (t + h) % K_MOD;
    }

    static String solve() {
        long[] gSmall = new long[8];
        for (int n = 1; n < 8; n++) {
            gSmall[n] = goodsteinSmall(n);
        }

        long ans = 0;
        for (int n = 1; n < 16; n++) {
            ans = (ans + computeGMod(n, gSmall)) % K_MOD;
        }
        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}