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 a positive integer $n$ the $n$th weak Goodstein sequence begins with $g_1=n$, and for $k\gt 1$ the term $g_k$ is built from $g_{k-1}$ in three moves: write $g_{k-1}$ in base $k$, read that same digit string as a base $k+1$ number, and subtract $1$. The sequence stops as soon as a term is $0$. For $n=6$ it runs $\{6,11,17,25,\dots\}$: $6=110_2$ and $110_3=12$, so $g_2=11$, and Section 1 continues the walk. Every such sequence can be shown to terminate. Writing $G(n)$ for the number of nonzero elements of the $n$th sequence, the statement supplies $G(2)=3$, $G(4)=21$, $G(6)=381$ and the total $2517$ for $\sum G(n)$ over $1\le n\lt 8$. We must find the last nine digits of $\sum G(n)$ over $1\le n\lt 16$. ^[1]

The whole solution in one sentence: every step raises the base by exactly one, so $G(n)$ is just the final base minus $2$, and following how the base grows while the four binary digits of $n$ are emptied one position at a time turns $G(n)$ into the map $x\mapsto x\,2^{x}$ iterated a few thousand times, which we evaluate modulo $10^9$ by carrying every exponent modulo a short chain of shrinking moduli.

Mathematical Approach

1. One step is a subtraction in the new base

The values soon grow far too large to store, but the digit strings stay short, never more than four positions in this problem, even though the base and the individual digits become enormous. Suppose the base-$k$ digits of $g_{k-1}$ are $c_m,\dots,c_1,c_0$, so that

$$g_{k-1}=\sum_{i=0}^{m}c_i\,k^{i},\qquad 0\le c_i\le k-1.$$

Every $c_i$ is at most $k-1\lt k+1$, so the same string is also legal in base $k+1$. Reading it there and subtracting $1$ therefore amounts to schoolbook subtraction of $1$ from that string in the new base: find the lowest $j$ with $c_j\ne 0$, replace $c_j$ by $c_j-1$, and fill the digits below it, all $0$, with the largest base-$k+1$ digit, namely $k$.

So we may forget the values and keep only a state: a digit string with its base. One step raises the base by one and subtracts $1$ from the string in the raised base, the borrow filling with the new base minus one, not the old one — the easiest mistake to make here. For $n=6$ the states are $(1,1,0)$ at base $2$, then $(1,0,2)$ at base $3$, then $(1,0,1)$ at base $4$, then $(1,0,0)$ at base $5$, which read as $6,\,102_3=11,\,101_4=17,\,100_5=25$, exactly the sequence in the statement.

2. Why every sequence stops

Two observations settle termination. First, the number of positions never grows: if the top position is $m$ and the base is $b$, every digit is at most $b-1$, so

$$\begin{aligned}&\sum_{i=0}^{m}(b-1)(b+1)^{i}\\&\qquad=\frac{b-1}{b}\big((b+1)^{m+1}-1\big)\end{aligned}$$

is smaller than $(b+1)^{m+1}$, and subtracting $1$ only lowers it, so the new value still fits in positions $0$ to $m$.

Second, read the string as the vector $(c_m,\dots,c_0)$, top position first. The step leaves every digit above $j$ untouched and strictly lowers $c_j$, so the vector strictly decreases lexicographically. Lexicographic order on $\mathbb{N}^{m+1}$ admits no infinite strictly decreasing chain, by induction on the number of positions: one coordinate cannot decrease forever, and between two drops of the top coordinate the remaining $m$ positions supply only finitely many steps. So the string reaches all zeros after finitely many steps, which is termination.

The classical phrasing attaches the ordinal $\omega^{m}c_m+\dots+\omega c_1+c_0$, with the coefficients on the right because ordinal multiplication is not commutative; comparing these Cantor normal forms is the lexicographic comparison above, and they stay below $\omega^{\omega}$. ^[3] The genuine Goodstein sequences, which rewrite the exponents too, need ordinals up to $\varepsilon_0$ and are famously unprovable in Peano arithmetic. ^[2]

3. The step count is the final base minus two

Each step raises the base by exactly $1$, and $g_1$ is read at base $2$. If the string first becomes all zeros at base $\beta(n)$, exactly $\beta(n)-2$ steps were taken, and the nonzero terms are the ones before that zero, so

$$G(n)=\beta(n)-2.$$

So we never count steps, only track the base. The brute-force run gives the final bases $3,5,7,23,63,383,2047$ for $n=1,\dots,7$, hence $G(n)=1,3,5,21,61,381,2045$, which match the statement's three checkpoints and sum to $2517$. Both sequences are catalogued. ^[4] ^[5]

Every formula below is cleaner if we carry the shifted base $C=B+1$ rather than $B$: base $B$ has the $B$ digits $0$ to $B-1$, and the shift turns the map $B\mapsto 2B+1$ of Section 5 into $C\mapsto 2C$.

4. Emptying the units digit

Start from a state with digit $c_0$ at position $0$ and base $b$. While $c_0\gt 0$ the subtraction touches nothing else, lowering $c_0$ by one and raising the base, so after $c_0$ steps position $0$ holds $0$ and the base is $b+c_0$. Since $n\lt 16$ its binary digits $c_3,c_2,c_1,c_0$ are $0$ or $1$, and the starting base is $2$, so this stage leaves

$$C_1=3+c_0.$$

5. Position one: the shifted base doubles

Now let the state be $(a,0)$ at base $B$, meaning digit $a\ge 1$ at position $1$ and $0$ at position $0$. One step must borrow: the state becomes $(a-1,B)$ at base $B+1$, because the new base is $B+1$ and the fill digit is $B$. Section 4 then empties position $0$ in $B$ further steps, ending at base $2B+1$ with the state $(a-1,0)$. In terms of $C=B+1$ one such round is simply

$$C\longmapsto 2C,$$

and it lowers the digit at position $1$ by one. After $a$ rounds position $1$ is empty and $C$ has been multiplied by $2^{a}$, at a cost of $(B+1)(2^{a}-1)$ steps, the geometric sum $C+2C+\dots+2^{a-1}C$. For our starting states this gives

$$C_2=2^{c_1}C_1.$$

Two examples to check by hand: $n=2$ has $c_0=0$, $c_1=1$, so $C_2=6$ and $G(2)=3$; $n=3$ has $c_0=c_1=1$, so $C_2=8$ and $G(3)=5$. Both agree with the simulation.

6. Position two: the map $x\mapsto x\,2^{x}$

Next let the state be $(a,0,0)$ at base $B$, with $a\ge 1$ at position $2$. The borrow step turns it into $(a-1,B,B)$ at base $B+1$. Emptying position $0$, whose digit is $B$, then raises the base to $2B+1$ by Section 4, so the shifted base has grown from $C=B+1$ to $2C$. Position $1$ still holds $B=C-1$, so Section 5 multiplies that $2C$ by a further $2^{C-1}$. The net effect of one position-two round is

$$\begin{aligned}C&\longmapsto 2^{C-1}\cdot 2C\\&=C\,2^{C}=u(C),\\ u(x)&=x\,2^{x},\end{aligned}$$

and it lowers the digit at position $2$ by one. Writing $u^{(r)}$ for $u$ applied $r$ times, emptying position $2$ leaves

$$C_3=u^{(c_2)}(C_2).$$

One round costs $u(C)-1-B=C(2^{C}-1)$ steps. For $n=4$ we get $C_2=3$ and $C_3=u(3)=24$, hence $G(4)=21$; for $n=6$, $C_2=6$ and $C_3=u(6)=384$, hence $G(6)=381$ — the statement's own checkpoints. A digit-by-digit simulation confirms the round at full size too: from $n=8$ the base reaches $23$ after $21$ steps with the state $(2,0,0)$, and the next position-two round takes exactly $24\cdot(2^{24}-1)=402653160$ further steps, ending at base $402653183=u(24)-1$.

7. Position three: $u$ iterated $C$ times

Finally let the state be $(1,0,0,0)$ at base $B$, the case we need, since $c_3\le 1$. One step gives

$$\begin{aligned}(B+1)^{3}-1&=B(B+1)^{2}\\&\quad+B(B+1)+B,\end{aligned}$$

so the state becomes $(0,B,B,B)$ at base $B+1$: the top digit is gone and the three positions below are full. Emptying them by Sections 4, 5 and 6 in turn sends $C=B+1$ to $2C$, then to $2^{C-1}\cdot 2C=u(C)$, and then, position $2$ also holding $B=C-1$, through $C-1$ further applications of $u$. So one position-three round is

$$C\longmapsto u^{(C)}(C),$$

with $u$ applied $C$ times — the number of iterations is itself the current value. That self-reference is where the explosive growth comes from.

8. The closed form for $G(n)$

Chaining the four stages gives the formula. With $c_3,c_2,c_1,c_0$ the binary digits of $n\lt 16$, set

$$\begin{aligned}C_1&=3+c_0,\\ C_2&=2^{c_1}C_1,\\ C_3&=u^{(c_2)}(C_2),\end{aligned}$$

$$C_4=\begin{cases}u^{(C_3)}(C_3)&c_3=1\\ C_3&c_3=0\end{cases}$$

and $G(n)=C_4-3$, since the final base is $C_4-1$ and Section 3 subtracts $2$. For $1\le n\le 7$ this gives $C_4=4,6,8,24,64,384,2048$, hence $G(n)=1,3,5,21,61,381,2045$, the simulated values.

Extend the notation by $G(0)=0$, a sequence starting at $0$ having no nonzero term. For $8\le n\le 15$ write $t=G(n-8)$. The lower three positions of $n$ carry the binary digits of $n-8$, and the leading digit cannot be touched until they are empty, so stage three ends with

$$\begin{aligned}C_3&=t+3,\\ G(n)&=u^{(t+3)}(t+3)-3.\end{aligned}$$

For $n=8$ this is $u^{(3)}(3)-3$, and $u(3)=24$, $u(24)=24\cdot 2^{24}=3\cdot 2^{27}$, so

$$G(8)=3\cdot 2^{402653211}-3,$$

exactly the value recorded as the next term of the catalogued step-count sequence. ^[4] Modulo $10^9$ it is $722374141$, the checkpoint the C++ program asserts.

9. The telescoping sum the program evaluates

The solutions never iterate $u$ directly; they accumulate

$$\begin{aligned}H(s,m)&=\sum_{i=0}^{m}s_i\big(2^{s_i}-1\big),\\ s_{i+1}&=u(s_i),\end{aligned}$$

with $s_0=s$, and then return $t+H(t+3,t+2)$. This is the same quantity: each summand is a round cost from Section 6, and

$$s_i\big(2^{s_i}-1\big)=s_i2^{s_i}-s_i=s_{i+1}-s_i,$$

so the sum telescopes to $H(s,m)=s_{m+1}-s_0$. With $s_0=t+3$ and $m=t+2$ there are $t+3$ summands, so $s_{m+1}=u^{(t+3)}(t+3)$ and

$$\begin{aligned}&t+H(t+3,t+2)\\&\qquad=u^{(t+3)}(t+3)-3,\end{aligned}$$

which is the closed form of Section 8. The telescoping also explains the mysterious-looking offsets $+3$ and $+2$ in the code: they are $C_3$ and $C_3-1$.

10. How large these numbers are

For $n=15$ we have $t=2045$, so $s_0=2048$ and $u$ is applied $2048$ times. Already $s_1=2048\cdot 2^{2048}$ has $620$ decimal digits, and $s_2=s_1 2^{s_1}$ has more than $10^{619}$ decimal digits. Storing such integers in full is out of the question, so we work modulo $10^9$ from the start. Only $G(1),\dots,G(7)$ are computed exactly, by the simulation of Section 1.

11. Powers of two modulo $10^9$

Everything reduces to $2^{s}\bmod 10^9$ when $s$ is known only through residues. Since $10^9=2^9\cdot 5^9$ is not coprime to $2$, Euler's theorem cannot be applied to the whole modulus: the two prime powers must be handled apart and glued together.

Modulo $2^9=512$: if $s\ge 9$ then $2^{s}\equiv 0$. Only small starting values need care, and there $s$ is known exactly.
Modulo $5^9$: here $2$ is invertible, and its powers $2,4,3,1$ modulo $5$ first return to $1$ at the exponent $4=\varphi(5)$, so $2$ is a primitive root modulo $5$. For an odd prime $p$ such a root stays one modulo every $p^k$ unless its $(p-1)$st power is $1$ modulo $p^2$, and $2^4=16\not\equiv 1\pmod{25}$. ^[6]

The multiplicative order of $2$ is therefore the full group order,

$$\operatorname{ord}_{5^k}(2)=\varphi(5^k)=4\cdot 5^{k-1},$$

confirmed by direct computation for $k\le 9$. So $2^{s}\bmod 5^9$ depends only on $s\bmod 4\cdot 5^{8}=1562500$.

Write $p_2$ and $p_5$ for the remainders of $2^{s}$ modulo $512$ and modulo $5^9$. As $512$ and $5^9$ are coprime, exactly one residue modulo $10^9$ has both: every candidate matching $p_2$ has the shape $p_2+512\,v$, and matching $p_5$ then forces $v\equiv(p_5-p_2)\cdot 512^{-1}\pmod{5^9}$, the inverse coming from the extended Euclidean algorithm. ^[7]

12. The chain of moduli

Getting $2^{s}\bmod 10^9$ now needs $s\bmod 4\cdot 5^{8}$. But $s=s_{i}$ is itself $s_{i-1}2^{s_{i-1}}$, so that residue needs $2^{s_{i-1}}\bmod 4\cdot 5^{8}$, which by the same analysis needs $s_{i-1}\bmod 4\cdot 5^{7}$, and so on. Each level drops the $5$-exponent by one, so the recursion is finite and the solutions carry every $s_i$ modulo the ten moduli

$$10^9,\;4\cdot 5^{8},\;4\cdot 5^{7},\;\dots,\;4\cdot 5,\;4,$$

updating all ten with $s_{i+1}\equiv s_i\cdot 2^{s_i}$ each iteration. It bottoms out at $4$, where $2^{s}\equiv 0$ for every $s\ge 2$ and no exponent is needed. The factor $4=2^2$ is kept in every entry so that each modulus is a multiple of the order $4\cdot 5^{k-1}$ the level above it needs.

How the Code Works

All three programs share one structure. goodstein_small performs the digit simulation of Section 1 literally: it peels the base-$b$ digits off by repeated division, rebuilds the number in base $b+1$ by Horner's rule, subtracts $1$ and raises the base. It runs only for $n\lt 8$, where no term exceeds $1023$ and no rebuilt value $1024$. Python and Java fill the cache and add up the residues in solve, work that in C++ sits in main; the cache is g_small, or gSmall in Java, with entry $0$ set to $0$. The snake_case names below are the C++ and Python ones, which Java spells goodsteinSmall, computeGMod, computeHMod, buildModChain and pow2ModComposite.

compute_g_mod implements Section 8: for $n\lt 8$ it returns the cached value, otherwise it sets t = g_small[n-8], s = t+3, m = t+2 and returns $t+H(s,m)$ modulo $10^9$.

build_mod_chain builds the ten ModEntry records of Section 12, each holding its modulus, the exponents pow2 and pow5 of its $2$- and $5$-parts, the modulus period_mod in which its own exponent must be known, and the inverse inv_mod2_in5 for the Chinese remainder step. The head is $10^9$ with period_mod $4\cdot 5^{8}$; the tail $4$ has period_mod zero.

pow2_mod_composite realises Section 11 for one entry: it reads the residue of $s$ at that entry's period_mod, raises $2$ to it modulo the $5$-part, takes the $2$-part from the exact value of $s$ when that is known and below pow2 and otherwise sets it to $0$, then glues the two by the Chinese remainder formula; a pure power-of-two entry returns at once. Only C++ spells the binary exponentiation out, in mod_pow, ^[8] with egcd and mod_inv for the inverse; Python defers to the built-in pow, Java to BigInteger.modPow and modInverse.

compute_h_mod runs the loop of Section 9. It keeps the vector residues of $s_i$ modulo every chain entry, adds residues[0] * (p_main - 1) to the running sum, then fills a second vector next_residues by multiplying each residue by $2^{s_i}$ in its own modulus, so every exponent read during one refresh belongs to the same $s_i$; the last iteration adds its summand and stops. Alongside it carries exact_s and the flag exact_known, tracking $s_i$ exactly while it stays below about $1000$ — at most two iterations. That matters only for the $2^9$ part: once the flag is cleared, every later $s_i$ provably exceeds $1000$ and hence $9$, so $2^{s_i}\equiv 0\pmod{2^9}$ is correct.

The language differences are in the arithmetic types. A product of two residues below $10^9$ stays below $10^{18}$ and so still fits a signed $64$-bit integer, yet C++ casts modular products to __int128 and Java routes them, with modPow and modInverse, through BigInteger — belt and braces rather than necessity. Python needs neither, its integers being unbounded. Only the C++ file has self-checks: run_checkpoints verifies $G(2)=3$, $G(4)=21$, $G(6)=381$ and the total $2517$, compares compute_h_mod(2,2) against brute_h_mod_small, which forms $H(2,2)$ exactly with Boost cpp_int and reduces only the completed sum, ^[9] and asserts $H(3,2)\equiv 722374141$; parse_arguments accepts --skip-checkpoints to bypass them.

Complexity and Verification

The exact phase costs $\sum_{n=1}^{7}G(n)=2517$ digit steps on numbers of at most three digits, so it is free. The modular phase runs $s_0=G(n-8)+3$ iterations for each $n$ from $8$ to $15$, that is $3+4+6+8+24+64$${}+384+2048=2541$ in total, the longest run being $2048$. Every summand costs one power at the head modulus, and every iteration but the last of a run also refreshes all ten residues, nine of them by a modular exponentiation (the modulus-$4$ entry needs none): counting the calls gives $2533$ refreshes and $25338$ exponentiations. No reduced exponent exceeds $21$ bits, so square-and-multiply spends at most $42$ modular multiplications on each, ^[8] a few hundred thousand in all on numbers below $10^9$. Memory is $O(1)$: ten modulus records and two residue vectors. Built with g++ -O2, the C++ binary finishes in about $0.01$ seconds, and all three programs print the same nine digits. One pitfall survives: the task wants the last nine digits while the programs print a residue modulo $10^9$, so a residue of fewer than nine digits would need its leading zeros restored by hand.

For this article the method was re-derived and checked independently. A brute-force simulation of the statement's own rule reproduces $G(2)=3$, $G(4)=21$, $G(6)=381$ and the total $2517$, and on every state those seven sequences visit both claims of Section 2 hold. The round formulas of Sections 4 to 6 were verified exhaustively for all digits and all bases $b\le 8$, and the bottom-three closed form behind Section 7 in all $267$ digit triples with $2\le b\le 8$, each digit between $0$ and $b-1$, whose predicted final base stays below $3\cdot 10^{8}$ — no failure — while the full-scale run quoted in Section 6 confirms the map $u$ at its real size. The telescoping identity of Section 9 was checked exactly for small parameters. Finally, compute_h_mod was matched against an independent exact evaluation of $H(s,m)$ in sixteen cases, which carries $s_i$ as a true integer up to $s_m$ and reduces $s_{m+1}=s_m2^{s_m}$ by one built-in modular exponentiation on that exact exponent; they include $H(2048,1)$ and $H(2045,1)$, whose $s_1$ is the $620$-digit integer the $n=15$ run really meets. It was matched too against a second engine on a different principle: the Carmichael function $\lambda$, with $2^{e}\equiv 2^{(e\bmod\lambda(M))+\lambda(M)}\pmod M$ and no $2/5$ split or Chinese remainder step — a reduction that needs the raised exponent to reach the exponent of $2$ in $M$, which we checked holds at every link of its chain. The two engines agree for every $n$ from $1$ to $15$, and none of this is part of the published solutions.

Footnotes and References

The sources below supply the background facts used above; the digit states, the termination argument, the three round maps, the closed form and the telescoping identity are derived in full here.

Project Euler 396 — Weak Goodstein Sequence. The official statement defines the nth weak Goodstein sequence by g(1) = n and, for k > 1, by writing g(k−1) in base k, reading that string in base k+1 and subtracting 1; it gives 6, 11, 17, 25, … for n = 6, says every such sequence terminates, defines G(n) as the count of nonzero elements, lists G(2) = 3, G(4) = 21, G(6) = 381 and the total 2517 for 1 ≤ n < 8, and asks for the last nine digits of that total for 1 ≤ n < 16.
Wikipedia — Goodstein's theorem. The article treats the hereditary variant, in which the exponents are rewritten in the base as well. It states that every such sequence terminates at 0, proves it by replacing each base by ω and observing that the resulting ordinals strictly decrease, and reports that Kirby and Paris showed the theorem unprovable in Peano arithmetic. Our flat variant needs only the elementary argument of Section 2.
Wikipedia — Ordinal arithmetic. The article gives the Cantor normal form, with its coefficients written to the right of the powers, and the rule for comparing two ordinals in it: compare the exponents and coefficients in turn, and at the first inequality the larger component wins — the lexicographic comparison of digit vectors used in Section 2.
OEIS A266203 — number of steps for the weak Goodstein sequence to reach zero. It lists 0, 1, 3, 5, 21, 61, 381, 2045, which is G(n) for n = 0 to 7 in our notation, and a comment gives the next term as 3·2^402653211 − 3, the value Section 8 derives for G(8).
OEIS A056041 — final base of the weak Goodstein sequence. It lists 2, 3, 5, 7, 23, 63, 383, 2047 for n = 0 to 7, describes the same base-change process, and records the comment a(n) = A266203(n) + 2 — the identity of Section 3, that the step count is the final base minus the starting base 2.
Wikipedia — Primitive root modulo n. The article defines the multiplicative order of a as the lowest power of a congruent to 1, lists 2 among the primitive roots modulo 5 and modulo 25, and states that, for an odd prime p, a primitive root g modulo p remains one modulo every p^k unless g^(p−1) ≡ 1 (mod p²). For g = 2, p = 5 that exception fails, giving the order 4·5^(k−1) of Section 11.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli the article states that exactly one residue below their product has prescribed remainders, and its constructive proof for two moduli builds it from Bézout coefficients found by the extended Euclidean algorithm — the step that glues the two parts of a power of two in Section 11.
Wikipedia — Modular exponentiation. The article describes the right-to-left binary method, or exponentiation by squaring, which squares the base and conditionally multiplies once per bit of the exponent, and gives its running time as O(log exponent) against the O(exponent) of the naive methods — the cost model used above.
Boost.Multiprecision — cpp_int. The page documents cpp_int as a signed integer type of arbitrary precision whose values fall back on dynamic allocation once they outgrow the internal storage, illustrating it by computing 100! exactly. The C++ self-check forms H(2, 2) with it and reduces only the completed sum.

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());
    }
}