Problem 414: Kaprekar Constant

Project Euler Problem 414 Solution

EulerSolve provides an optimized solution for Project Euler Problem 414, Kaprekar Constant, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each base $B=6k+3$ with $2 \le k \le 300$, we consider every integer $i$ with $0 \lt i \lt B^5$, write it as a 5-digit base-$B$ string with leading zeros allowed, sort its digits in decreasing and increasing order, subtract, and repeat. For these bases, the 5-digit routine has a unique non-zero Kaprekar constant $C_B$. Define $s_B(i)=0$ when $i=C_B$ or when the 5 digits of $i$ are all equal, and otherwise let $s_B(i)$ be the number of iterations needed to reach $C_B$. Then $$S(B)=\sum_{0 \lt i \lt B^5} s_B(i),$$ and the goal is the last $18$ digits of $$\sum_{k=2}^{300} S(6k+3).$$ Mathematical Approach Step 1: Compress Every 5-Digit State to Two Gaps Let the sorted digits of the current 5-digit string be $$x_0 \le x_1 \le x_2 \le x_3 \le x_4.$$ The descending and ascending numbers are $$N_{\downarrow}=x_4 B^4 + x_3 B^3 + x_2 B^2 + x_1 B + x_0,$$ $$N_{\uparrow}=x_0 B^4 + x_1 B^3 + x_2 B^2 + x_3 B + x_4.$$ The Kaprekar subtraction is therefore $$N_{\downarrow}-N_{\uparrow}=(x_4-x_0)(B^4-1)+(x_3-x_1)(B^3-B).$$ So the next value depends only on the two outer gaps $$u=x_4-x_0,\qquad v=x_3-x_1,$$ with $0 \le v \le u \le B-1$. Define $$R_B(u,v)=u(B^4-1)+vB(B^2-1).$$ To advance one step, the implementation writes $R_B(u,v)$ in base $B$, sorts its five digits, and recomputes the new pair $(u',v')$....

Detailed mathematical approach

Problem Summary

Write a number with a fixed number of digits (padding with leading zeros), sort the digits once downwards and once upwards, and subtract the smaller arrangement from the larger. In base ten with four digits this Kaprekar step drives almost everything to $6174$: from $0837$ the statement finds $8730-0378=8352$, then $8532-2358=6174$, and $6174$ reproduces itself since $7641-1467=6174$. Four equal digits give $0$ instead. Repeating the step until one of the two happens is the Kaprekar routine. ^[1] ^[4]

This problem asks the same with five digits in a base $b$. A Kaprekar constant need not exist then: the routine must end in a fixed point or a repeating cycle, ^[2] and both a cycle and several different constants are possible. For five digits and $b=6t+3\ne9$ one does exist, and the statement names $(10,4,14,9,5)_{15}$ and $(14,6,20,13,7)_{21}$. Write $C_b$ for it. The function $s_b(i)$ is $0$ when $i=C_b$ and when the five base-$b$ digits of $i$ are equal, and otherwise is the number of iterations needed to reach $C_b$. Then

$$S(b)=\sum_{0\lt i\lt b^5}s_b(i),$$

with $S(15)=5274369$ and $S(111)=400668930299$ given. We want the last $18$ digits of $\sum_{k=2}^{300}S(6k+3)$, a sum over the $299$ bases $15,21,\dots,1803$. ^[1]

The whole solution in one sentence: one Kaprekar step looks only at two differences between sorted digits, so each base collapses from $b^5$ starting values to about $b^2/2$ states, and $S(b)$ becomes a weighted sum of distances in that small graph.

Mathematical Approach

1. One step, done by hand

Fix $b=15$ and start at $i=3874$, whose five base-$15$ digits are $(0,1,2,3,4)$. Sorted downwards they read $(4,3,2,1,0)_{15}=213090$; sorted upwards they read $3874$ again. The step therefore sends $3874$ to $213090-3874$, which is $209216=(4,1,14,12,11)_{15}$. Nine further steps reach $523040=(10,4,14,9,5)_{15}$, the constant $C_{15}$, so $s_{15}(3874)=10$. Base $15$ already has $759374$ starting values to simulate one by one, and base $1803$ has about $1.9\cdot10^{16}$. We must see what a step really depends on.

2. A step sees only two gaps

Let the sorted digits of the current value be $d_0\le d_1\le d_2\le d_3\le d_4$. The two arrangements are

$$\begin{aligned}A&=d_4b^4+d_3b^3+d_2b^2+d_1b+d_0,\\ B&=d_0b^4+d_1b^3+d_2b^2+d_3b+d_4,\end{aligned}$$

and subtracting removes the middle digit completely:

$$A-B=(d_4-d_0)(b^4-1)+(d_3-d_1)(b^3-b).$$

Only two numbers survive. Put $p=d_4-d_0$ and $q=d_3-d_1$; then $0\le q\le p\le b-1$. Define

$$X(p,q)=p(b^4-1)+q(b^3-b).$$

In Section 1, $p=4-0=4$ and $q=3-1=2$, and indeed $X(4,2)=4\cdot 50624+2\cdot 3360$, which is $209216$.

Call $(p,q)$ a state. There are only $b(b+1)/2$ states instead of $b^5$ starting values, and values sharing a state share their successor: in base $15$, $34=(0,0,0,2,4)_{15}$, $62=(0,0,0,4,2)_{15}$ and $259=(0,0,1,2,4)_{15}$ all have state $(4,2)$ and all step to $209216$.

A by-product: $b^4-1$ and $b^3-b$ are multiples of $b-1$, so every value a step produces is divisible by $b-1$ ^[3] — in base $15$, by $14$, and $523040=14\cdot 37360$.

3. The digits of the result, in closed form

Writing $X(p,q)$ in base $b$ needs no sorting at all. Expand

$$X(p,q)=pb^4+qb^3+0\cdot b^2-qb-p$$

and clear the negative coefficients from the units place upwards, assuming $1\le q\le p\le b-1$. The units coefficient $-p$ becomes the digit $b-p$ and borrows $1$; the $-q-1$ in the $b$ place becomes $b-1-q$ and borrows $1$; the $-1$ in the $b^2$ place becomes $b-1$ and borrows $1$; the $b^3$ coefficient becomes $q-1$ and the borrowing stops. Hence

$$X(p,q)=(p,\,q-1,\,b-1,\,b-1-q,\,b-p)_b.$$

For $q=0$ and $p\ge1$ the same computation gives $(p-1,b-1,b-1,b-1,b-p)_b$. Only the state $(0,0)$ is left out, and there $X=0$ has five zero digits. The example confirms the formula: $p=4$, $q=2$, $b=15$ predict $(4,1,14,12,11)_{15}$, exactly the digits found by hand.

Two conclusions follow. For $p\ge1$ every digit listed is at most $b-1$ and one equals $b-1$, so the largest digit of $X(p,q)$ is $b-1$; and the successor state, written $T(p,q)$, comes from sorting five explicitly known quantities, so one transition is constant work.

4. Where the constant comes from

A value is left alone by the routine exactly when it equals $X(p,q)$ for its own state and $T(p,q)=(p,q)$. Different states give different values, since $(b-1)(b^3-b)\lt b^4-1$ means the $q$ part cannot undo a change in $p$; a nonzero fixed value and a nonzero fixed state are thus the same thing, and the case analysis below finds all of them.

Let $1\le q\le p\le b-1$ and let $E$ be the digit multiset $\{p,q-1,b-1,b-1-q,b-p\}$. Its maximum is $b-1$, so the new first gap is $b-1-\min E$ and $T(p,q)=(p,q)$ forces $\min E=b-1-p$. The minimum is not $p$, since $q-1\lt p$, nor $b-p$, since $b-p\gt b-1-p$. Two cases remain.

$b-1-q=b-1-p$, that is $q=p$. Then $\min E=b-1-p\le p-1$ gives $b\le 2p$, and sorting $E$ makes the new second gap $1$ if $b=2p$ and $2p-b$ if $b\lt 2p$; equating either to $q=p$ forces $b=2$ or $p=b$, so nothing survives for $b\ge3$.
$q-1=b-1-p$, that is $p+q=b$. Then $b-1-q=p-1$ and $b-p=q$, so $E=\{q-1,q,p-1,p,b-1\}$. If $p\ge q+2$ this list is already sorted, the new second gap is $p-q$, and $T(p,q)=(p,q)$ means $p-q=q$; together with $p+q=b$ this gives $q=b/3$ and $p=2b/3$. In the leftover cases $p=q$ and $p=q+1$, sorting $E$ makes the new second gap $1$, so a fixed state needs $q=1$; that gives $b=2$ and $b=3$, and base $3$ has $p=2$, $q=1=b/3$.

The case $q=0$ is quicker: $\{p-1,b-p,b-1,b-1,b-1\}$ has $b-1$ as its second largest element, so a vanishing second gap needs a fourth copy of $b-1$, forcing $p=1$; but then the new first gap is $b-1\ne1$ for $b\ge3$.

So for $b\ge3$ a nonzero fixed point exists if and only if $3\mid b$, and it is then unique. With $m=b/3$,

$$C_b=X(2m,m)=m(b^2-1)(2b^2+b+2),$$

and by Section 3 its digits are $(2m,m-1,3m-1,2m-1,m)_b$: for $b=15$ that is $(10,4,14,9,5)_{15}=523040$ and for $b=21$ it is $(14,6,20,13,7)_{21}=2787400$, the values printed in the statement, while $b=3$ gives $184$ and $b=1803$ gives $12705969088986184$.

5. Which bases work, and why $9$ is left out

A nonzero fixed point is not enough; the statement guarantees that every starting value whose five digits are not all equal reaches it. Section 4 already discards every base not divisible by $3$, and a second argument discards every even base $b=2h$, $h\ge2$. Apply Section 3 to the state $(h,h-1)$: its digits are $(h,h-2,2h-1,h,h)_b$, with gaps $(h+1,0)$. Apply it again to $(h+1,0)$: the digits are $(h,2h-1,2h-1,2h-1,h-1)_b$, with gaps $(h,h-1)$. Those two states form a two-cycle, so the routine never settles; in base $12$ they are $(6,5)$ and $(7,0)$, and $132990$ and $145145$ alternate for ever.

Only bases $b\equiv3\pmod 6$ are left, and base $9$ still fails: it has the fixed point $41520$ but also a cycle of length five, which $41$ of its $45$ states run into. ^[3]

$$\begin{aligned}34960\to40080\to55360\\ \to49520\to42240\to34960.\end{aligned}$$

That is why the statement excludes $b=9$. What remains matches the statement: for $2\le b\le 60$, the bases in which every value whose digits are not all equal reaches one nonzero constant are precisely $3,15,21,27,33,39,45,51$ and $57$. For the $299$ bases the problem uses, all $163\,760\,805$ states were checked: every state with $p\ge1$ reaches $(2m,m)$, no cycle longer than one state occurs, and the only self-loops are $(0,0)$ and $(2m,m)$.

6. How many numbers share a state

Let $M(p,q)$ be the number of five-digit base-$b$ strings of state $(p,q)$. Write the sorted digits as $a+f_j$ with $f_0=0$, so $a$ is the smallest digit; then $f_4=p$, the offsets rise weakly, and $f_3-f_1=q$. The smallest digit is free as long as $a+p\le b-1$, leaving $b-p$ choices and splitting off a factor:

$$M(p,q)=(b-p)\,W(p,q),$$

where $W(p,q)$ counts offset patterns with their arrangements. A multiset of five digits with multiplicities $k_1,\dots,k_r$ has $5!/(k_1!\cdots k_r!)$ arrangements, ^[6] so a pattern contributes $120$, $60$, $30$, $20$, $10$, $5$, or $1$ when all five digits agree. Three disjoint cases exhaust the offsets.

$q=0$ forces $f_1=f_2=f_3=c$. The choices $c=0$ and $c=p$ give a fourfold digit, $5$ arrangements each; the $p-1$ inner values give a threefold digit, $20$ each. So $W(p,0)=10+20(p-1)$ for $p\ge1$, and $W(0,0)=1$.
$q=p\ge1$ forces $f_1=0$ and $f_3=p$. The choices $f_2=0$ and $f_2=p$ give $10$ each and the $p-1$ inner ones $30$, so $W(p,p)=20+30(p-1)$.
For $0\lt q\lt p$, either $f_1=0$, or $f_3=p$, or neither — and the first two cannot hold together, as that would mean $q=p$. If $f_1=0$ then $f_3=q$ and $f_2$ runs over $0,\dots,q$, contributing $20+30+60(q-1)$; the reflection $f\mapsto p-f$ carries this family onto the case $f_3=p$, which contributes the same. If neither holds, $f_1$ runs over $1,\dots,p-1-q$ and $f_2$ over $f_1,\dots,f_1+q$, contributing $60+60+120(q-1)$ for each $f_1$.

Adding the three parts of the last case,

$$\begin{aligned}W(p,q)&=100+120(q-1)\\&\quad+120q(p-1-q).\end{aligned}$$

No special case for $p=q+1$ is needed: the factor $p-1-q$ is then $0$. As a check, $W(4,2)=100+120+240=460$ and $M(4,2)=11\cdot 460=5060$ in base $15$, so $5060$ values take exactly the $10$ steps found in Section 1.

7. Distance to the constant

Every state has exactly one successor, so the states form a functional graph, in which each vertex has out-degree one. ^[5] For $p\ge1$, let $D(p,q)$ be the number of further steps that lead from $X(p,q)$ to $C_b$:

$$D(p,q)=\begin{cases}0,&T(p,q)=(p,q),\\ 1+D\bigl(T(p,q)\bigr),&T(p,q)\ne(p,q).\end{cases}$$

Section 5 guarantees that this recursion terminates, because every path from a state with $p\ge1$ ends at the self-loop $(2m,m)$. The state $(0,0)$ loops onto itself too, so the same recursion returns $0$ for it; Section 8 never needs that value. Keeping each computed value in a table — memoisation ^[8] — costs one transition per state in total.

8. Assembling $S(b)$

Each $i$ with $0\lt i\lt b^5$ falls into one of three groups. The $b-1$ values with five equal digits have state $(0,0)$ and $s_b(i)=0$ by definition; $C_b$ gives $0$ as well; every other $i$ spends one step reaching $X(p,q)$ and then $D(p,q)$ more, so $s_b(i)=1+D(p,q)$. Counting those first steps gives $b^5-1-(b-1)-1=b^5-b-1$, and grouping the rest by state gives

$$S(b)=b^5-b-1+\sum_{p=1}^{b-1}\sum_{q=0}^{p}M(p,q)\,D(p,q).$$

The constant takes care of itself: its state $(2m,m)$ has depth $0$, so the double sum gives it nothing, and it was already removed from the first-step count. In base $15$ that state holds $M(10,5)=14900$ values, and the $14899$ other than $C_{15}$ itself are exactly those with $s_{15}(i)=1$.

9. Base $3$ from start to finish

The smallest admissible base is $b=3$, where the whole machine fits on a few lines. There are six states; for each, the list gives $X(p,q)$ with its digits, the successor $T(p,q)$, the depth $D$ and the multiplicity $M$.

$(0,0)$: $X=0=(0,0,0,0,0)_3$, $T=(0,0)$, $D=0$, $M=3$.
$(1,0)$: $X=80=(0,2,2,2,2)_3$, $T=(2,0)$, $D=3$, $M=20$.
$(1,1)$: $X=104=(1,0,2,1,2)_3$, $T=(2,1)$, $D=1$, $M=40$.
$(2,0)$: $X=160=(1,2,2,2,1)_3$, $T=(1,1)$, $D=2$, $M=30$.
$(2,1)$: $X=184=(2,0,2,1,1)_3$, $T=(2,1)$, $D=0$, $M=100$.
$(2,2)$: $X=208=(2,1,2,0,1)_3$, $T=(2,1)$, $D=1$, $M=50$.

The multiplicities $3,20,40,30,100,50$ add up to $243=3^5$, as they must. The only self-loop with $p\ge1$ is $(2,1)$, so $C_3=184$, matching the closed form of Section 4 and the published table of five-digit Kaprekar cycles. ^[3] The weighted depths give

$$20\cdot3+40\cdot1+30\cdot2+100\cdot0+50\cdot1=210,$$

and $3^5-3-1=239$, so $S(3)=239+210=449$. Simulating all $242$ starting values one by one also gives $449$: $3$ have $s_3=0$, and $99$, $90$, $30$, $20$ have $s_3=1,2,3,4$ — the states $(2,1)$ without the constant itself, then $(1,1)$ and $(2,2)$, then $(2,0)$, then $(1,0)$. The longest run starts at $i=1$, of state $(1,0)$: it goes $1\to80\to160\to104\to184$, four steps, as $1+D(1,0)=4$ predicts.

10. The outer sum and the last $18$ digits

The $299$ values $S(6k+3)$ do not interact, so the outer sum is a plain addition — but a large one. Already $S(1803)=6334961735971658529$ has $19$ digits, and $S(15)=5274369$, $S(21)=34289199$, $S(27)=132630799$ show how fast they grow. Since only the last $18$ digits are wanted, everything may be reduced modulo $10^{18}$ as it is accumulated: congruences are compatible with addition and multiplication, ^[7] and Section 8 uses nothing else.

How the Code Works

The C++ program ^[10] follows the ten sections directly. parse_arguments reads the optional flags --k-from=, --k-to= and --skip-checkpoints, accepting them only when the lower end of the range is at least $2$ and the upper end is not below it. That bound keeps base $9$ out of reach: $k=1$ would give $b=9$, whose cycle would make the recursion of Section 7 run for ever.

next_state is Section 2 in five lines: it forms x ${}=p(b^4-1)+q(b^2-1)b$, peels off five base-$b$ digits with % and / into a std::array, sorts them and returns $(d_4-d_0,\,d_3-d_1)$. It recomputes the digits rather than using Section 3; that is still constant work and keeps the definition visible. ways is $M(p,q)$ of Section 6: the literals 120/24, 120/6, 120/2/2 and 120/2 are the arrangement counts $5$, $20$, $30$ and $60$, the three branches are the cases $q=0$, $q=p$ and $0\lt q\lt p$, the guard if (p - 2 >= q) adds the family with neither endpoint doubled, and the result is multiplied by b - p.

dfs_depth is the recursion of Section 7 over a b-by-b table memo initialised to $-1$. It takes a reference to the cell it will fill, recurses on next_state unless the state loops onto itself, and writes each depth once; only the triangle $q\le p$ is used, so half the table is idle. S_of_base loops over $p=1,\dots,b-1$ and $q=0,\dots,p$, accumulating ways times dfs_depth modulo $10^{18}$ with the single product formed in __uint128_t, ^[9] then adds total - b - 1, the first-step count of Section 8. Unless --skip-checkpoints is given, main first calls run_checkpoints, which compares S_of_base(15) with $5274369$ and S_of_base(111) with $400668930299$ and returns false on a mismatch, whereupon main exits with status $2$; both targets are below $10^{18}$, so the comparison is exact. solve_sum adds $S(6k+3)$ over the requested $k$, again modulo $10^{18}$.

The Java program repeats the same four per-base routines — nextState, ways, dfsDepth over an int[b][b], sOfBase — in signed long arithmetic, but has no argument parser, fixes the range at $2\le k\le300$ and runs no checkpoints. One inner detail differs too: instead of one multiplication it adds wMod to a running term exactly d times. That is correct but needlessly slow, and it buys no safety — the largest multiplicity anywhere is $M=26049804100$, and even times the largest depth $619$ it stays near $1.6\cdot10^{13}$, far inside the signed 64-bit range. IntStream.rangeClosed supplies the values of $k$ and mapToLong turns the per-base results into a LongStream, whose reduce(0, (a, b) -> (a + b) % MOD) combines them: on residues below $10^{18}$ that operation is associative, has identity $0$ and cannot overflow a long, as the method requires. ^[9]

The Python file is not a third implementation but a bridge. It takes the problem number from its own filename, locates the C++ source in the sibling solutionsCpp directory, compiles it with clang++ or g++ at -O2 into a cached hidden binary (recompiling only when the source is newer), and parses the output with ANSWER_RE and EQUAL_RE, falling back to the last non-empty line — which is what the C++ program prints. Since should_skip_cpp_checkpoints finds --skip-checkpoints in that source, the bridge always passes the flag, so the checkpoints are skipped on this path.

Complexity and Verification

For one base the C++ program clears a table of $b^2$ cells and then visits each state of the triangle once, at a constant cost per visit: five divisions and a sort of five values in next_state, a few multiplications in ways, one table lookup. That is $\Theta(b^2)$ time and $\Theta(b^2)$ memory, and for $b=1803$ the table holds $1803^2=3\,250\,809$ four-byte entries, about $13$ MB. The Python bridge runs this same program. Java pays more, because its repeated addition costs $D(p,q)$ operations at every state; with

$$A_b=\sum_{p=1}^{b-1}\sum_{q=0}^{p}D(p,q)$$

its time is $\Theta(b^2+A_b)$, and $A_b$ dominates: $A_{1803}=635\,141\,195$ against some $1.6\cdot10^{6}$ states. Over the whole range the number of states is

$$\sum_{k=2}^{300}\frac{(6k+3)(6k+4)}{2}=163\,760\,805,$$

of which the loops visit $163\,760\,506$, one zero state per base being skipped. The C++ cost thus grows like the cube of the largest base — against $b^5\approx1.9\cdot10^{16}$ values for that base alone — while Java adds $\sum_bA_b=48\,125\,267\,720$ further additions; on one current desktop machine C++ needed about $3.3$ seconds of processor work for the whole range and Java about $15.6$, which those additions account for. The recursion of Section 7 descends only the uncached part of a path, so at most $620$ calls are ever active; the longest chain for $b=1803$ has $619$ edges, and the order in which S_of_base drives it peaks at $614$ frames.

Besides the two built-in checkpoints, the following were verified for this article by separate computations.

ways agrees state by state with a full enumeration of all $b^5$ digit strings for $b=4,7,9,12$, and the multiplicities sum to $b^5$; the closed form $W(p,q)$ agrees with ways for every state of every base $3\le b\le 120$.
The digit formulas of Section 3 hold for every state with $p\ge1$ of every base $3\le b\le 200$; the closed form for $C_b$ holds for every multiple of $3$ up to $2000$; the two-cycle of Section 5 exists in every even base $2h\le 1000$; and for $3\le b\le 600$ a nonzero fixed state exists exactly for the multiples of $3$, always $(2m,m)$.
The state scan behind Section 5 used explicit cycle detection and an iterative depth pass, independent of the recursion in the published code.
An independent Python implementation of Sections 6 to 8 reproduces $S(15)=5274369$, $S(111)=400668930299$ and $S(3)=449$, and direct simulation agrees: all $759374$ base-$15$ values give $5274369$ with longest run $s_{15}(7)=13$, all $242$ base-$3$ values give $449$. The C++ program run with --k-from=2 --k-to=3 prints $39563568=S(15)+S(21)$.

Footnotes and References

The references supply the background facts used above. The two-gap reduction, the digits of $X(p,q)$, the classification of fixed states, the multiplicities $M(p,q)$ and the formula for $S(b)$ are derived in this article.

Project Euler 414 — Kaprekar Constant. The official statement recalls 6174 with the example 0837, states that for five digits and a base b = 6t+3 other than 9 a Kaprekar constant exists, names it in bases 15 and 21, defines s_b(i) and S(b) with the checkpoints S(15) = 5274369 and S(111) = 400668930299, and asks for the last 18 digits of the sum of S(6k+3) over 2 ≤ k ≤ 300.
Wikipedia — Kaprekar's routine. For a fixed base and digit count, every Kaprekar sequence ends in a fixed point or a repeating cycle, and zero is a trivial fixed point in every base, reached only from strings of equal digits when leading zeros are kept. The only base-ten Kaprekar constants are 495 for three digits and 6174 for four, by Prichett, Ludington and Lapenta (1981).
MathWorld — Kaprekar Routine. The two sorted arrangements share a digit sum, hence are congruent modulo b − 1, so every difference is divisible by b − 1. The page also tabulates the nonzero cycles by base and digit count: the five-digit entry is 184 for base 3 and, for base 9, the fixed point 41520 with the cycle (34960, 40080, 55360, 49520, 42240) — the checks of Sections 5 and 9.
Plus Magazine — Y. Nishiyama, "The mysterious number 6174". Credits D. R. Kaprekar of Devlali with the operation, devised in 1949, and reports that every four-digit base-ten number whose digits are not all equal reaches 6174 in at most seven steps.
Wikipedia — Functional graph. A directed graph in which every vertex has out-degree exactly one, equivalent to a function from a set to itself, with a self-loop wherever it fixes a point — the structure of Section 7.
Wikipedia — Multinomial theorem. The multinomial coefficient n! divided by the product of the factorials of the multiplicities is "the number of distinct ways to permute a multiset of n elements"; Section 6 uses n = 5.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition, subtraction and multiplication, which is what allows every partial sum and product to be reduced modulo 10^18.
Wikipedia — Memoization. Storing the result of a call to a pure function so that a repeated call returns it instead of recomputing, trading memory for speed; the depth table of Section 7 is exactly this.
GCC manual — 128-bit integers documents unsigned __int128 as an extension available where the target has a wide enough integer mode. Java — LongStream requires of reduce(identity, op) an associative accumulator and an identity with accumulator.apply(identity, x) equal to x, which addition modulo 10^18 with identity 0 satisfies here.
C++ — Euler414.cpp. The linked immutable revision contains next_state, ways, dfs_depth, S_of_base and run_checkpoints; it is the implementation record for this article. The extra verifications above are separate programs, not packaged in that file.

Problem 414 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 MOD = 1000000000000000000ULL;  // last 18 digits

struct Options {
    int k_from = 2;
    int k_to = 300;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return 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;
        }
        if (parse_int_after_prefix(arg, "--k-from=", options.k_from) ||
            parse_int_after_prefix(arg, "--k-to=", options.k_to)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.k_from >= 2 && options.k_to >= options.k_from;
}

u64 pow_u64(u64 base, int exp) {
    u64 r = 1;
    for (int i = 0; i < exp; ++i) {
        r *= base;
    }
    return r;
}

std::pair<int, int> next_state(int p, int q, int b) {
    const u64 b2 = static_cast<u64>(b) * static_cast<u64>(b);
    const u64 b4 = b2 * b2;
    u64 x = static_cast<u64>(p) * (b4 - 1ULL) +
            static_cast<u64>(q) * (b2 - 1ULL) * static_cast<u64>(b);

    std::array<int, 5> d{};
    for (int i = 0; i < 5; ++i) {
        d[static_cast<std::size_t>(i)] = static_cast<int>(x % static_cast<u64>(b));
        x /= static_cast<u64>(b);
    }
    std::sort(d.begin(), d.end());
    return {d[4] - d[0], d[3] - d[1]};
}

i64 ways(int p, int q, int b) {
    i64 t = 0;
    if (q == 0) {
        if (p == 0) {
            t = 1;                        // aaaaa
        } else {
            t += 120 / 24 * 2;            // aaaae, aeeee
            t += 120 / 6 * (p - 1);       // accce
        }
    } else if (p == q) {
        t += 120 / 2 / 2 * (p - 1);       // aacee
        t += 120 / 6 / 2 * 2;             // aaaee, aaeee
    } else {
        t += 120 / 2 * (q - 1);           // abcee
        t += 120 / 2 / 2;                 // abbee
        t += 120 / 6;                     // abeee
        t *= 2;

        if (p - 2 >= q) {
            t += 120 / 2 * (p - 1 - q) * 2;          // abbde, abdde
            t += 120 * (p - 1 - q) * (q - 1);        // abcde
        }
    }
    return static_cast<i64>(b - p) * t;
}

int dfs_depth(int p, int q, int b, std::vector<std::vector<int>>& memo) {
    int& cell = memo[static_cast<std::size_t>(p)][static_cast<std::size_t>(q)];
    if (cell != -1) {
        return cell;
    }
    const auto nxt = next_state(p, q, b);
    if (nxt.first == p && nxt.second == q) {
        cell = 0;
    } else {
        cell = dfs_depth(nxt.first, nxt.second, b, memo) + 1;
    }
    return cell;
}

u64 S_of_base(int b) {
    const u64 total = pow_u64(static_cast<u64>(b), 5);

    std::vector<std::vector<int>> memo(
        static_cast<std::size_t>(b),
        std::vector<int>(static_cast<std::size_t>(b), -1));

    u64 ans = 0ULL;
    for (int p = 1; p <= b - 1; ++p) {
        for (int q = 0; q <= p; ++q) {
            const i64 w = ways(p, q, b);
            const int d = dfs_depth(p, q, b, memo);
            ans = (ans + static_cast<u64>((u128)(w % static_cast<i64>(MOD)) * d % MOD)) % MOD;
        }
    }

    // +1 first-step contribution for all non-special numbers:
    // total numbers minus {0}, minus all-equal positive numbers (b-1), minus Kaprekar constant.
    ans = (ans + total - static_cast<u64>(b) - 1ULL) % MOD;
    return ans;
}

bool run_checkpoints() {
    if (S_of_base(15) != 5274369ULL) {
        std::cerr << "Checkpoint failed: S(15)\n";
        return false;
    }
    if (S_of_base(111) != 400668930299ULL) {
        std::cerr << "Checkpoint failed: S(111)\n";
        return false;
    }
    return true;
}

u64 solve_sum(int k_from, int k_to) {
    u64 ans = 0ULL;
    for (int k = k_from; k <= k_to; ++k) {
        const int b = 6 * k + 3;
        ans = (ans + S_of_base(b)) % MOD;
    }
    return ans;
}

struct SumWorkerArgs {
    int tid = 0;
    int thread_count = 1;
    int k_from = 2;
    int k_to = 300;
    u64 partial = 0ULL;
};

void* solve_sum_worker(void* raw) {
    auto* args = static_cast<SumWorkerArgs*>(raw);
    u64 local = 0ULL;
    for (int k = args->k_from + args->tid; k <= args->k_to; k += args->thread_count) {
        const int b = 6 * k + 3;
        local = (local + S_of_base(b)) % MOD;
    }
    args->partial = local;
    return nullptr;
}

u64 solve_sum_parallel(int k_from, int k_to) {
    const int total_k = k_to - k_from + 1;
    if (total_k <= 1) {
        return solve_sum(k_from, k_to);
    }

    long hw = sysconf(_SC_NPROCESSORS_ONLN);
    if (hw <= 1) {
        return solve_sum(k_from, k_to);
    }

    const int thread_count = std::min<int>(total_k, static_cast<int>(hw));
    std::vector<pthread_t> threads(static_cast<std::size_t>(thread_count));
    std::vector<SumWorkerArgs> args(static_cast<std::size_t>(thread_count));

    for (int t = 0; t < thread_count; ++t) {
        args[static_cast<std::size_t>(t)].tid = t;
        args[static_cast<std::size_t>(t)].thread_count = thread_count;
        args[static_cast<std::size_t>(t)].k_from = k_from;
        args[static_cast<std::size_t>(t)].k_to = k_to;
        args[static_cast<std::size_t>(t)].partial = 0ULL;
    }

    bool create_failed = false;
    int started = 0;
    for (int t = 0; t < thread_count; ++t) {
        if (pthread_create(&threads[static_cast<std::size_t>(t)], nullptr,
                           solve_sum_worker, &args[static_cast<std::size_t>(t)]) != 0) {
            create_failed = true;
            break;
        }
        ++started;
    }

    for (int t = 0; t < started; ++t) {
        pthread_join(threads[static_cast<std::size_t>(t)], nullptr);
    }

    if (create_failed) {
        return solve_sum(k_from, k_to);
    }

    u64 ans = 0ULL;
    for (const auto& a : args) {
        ans = (ans + a.partial) % MOD;
    }
    return ans;
}

}  // 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;
    }

    const u64 answer = solve_sum_parallel(options.k_from, options.k_to);
    std::cout << answer << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.Arrays;
import java.util.stream.IntStream;

public class Euler414 {
    static final long MOD = 1000000000000000000L;

    static long powU64(long base, int exp) {
        long r = 1;
        for (int i = 0; i < exp; ++i) {
            r *= base;
        }
        return r;
    }

    static int[] nextState(int p, int q, int b) {
        long b2 = (long) b * b;
        long b4 = b2 * b2;
        long x = (long) p * (b4 - 1L) + (long) q * (b2 - 1L) * (long) b;

        int[] d = new int[5];
        for (int i = 0; i < 5; ++i) {
            d[i] = (int) (x % b);
            x /= b;
        }
        Arrays.sort(d);
        return new int[] { d[4] - d[0], d[3] - d[1] };
    }

    static long ways(int p, int q, int b) {
        long t = 0;
        if (q == 0) {
            if (p == 0) {
                t = 1;
            } else {
                t += 120 / 24 * 2;
                t += 120 / 6 * (p - 1);
            }
        } else if (p == q) {
            t += 120 / 2 / 2 * (p - 1);
            t += 120 / 6 / 2 * 2;
        } else {
            t += 120 / 2 * (q - 1);
            t += 120 / 2 / 2;
            t += 120 / 6;
            t *= 2;

            if (p - 2 >= q) {
                t += 120 / 2 * (p - 1 - q) * 2;
                t += 120 * (p - 1 - q) * (q - 1);
            }
        }
        return (long) (b - p) * t;
    }

    static int dfsDepth(int p, int q, int b, int[][] memo) {
        if (memo[p][q] != -1) {
            return memo[p][q];
        }
        int[] nxt = nextState(p, q, b);
        if (nxt[0] == p && nxt[1] == q) {
            memo[p][q] = 0;
        } else {
            memo[p][q] = dfsDepth(nxt[0], nxt[1], b, memo) + 1;
        }
        return memo[p][q];
    }

    static long sOfBase(int b) {
        long total = powU64(b, 5);
        int[][] memo = new int[b][b];
        for (int i = 0; i < b; i++) {
            Arrays.fill(memo[i], -1);
        }

        long ans = 0;
        for (int p = 1; p < b; ++p) {
            for (int q = 0; q <= p; ++q) {
                long w = ways(p, q, b);
                int d = dfsDepth(p, q, b, memo);

                long wMod = w % MOD;
                long term = 0;
                for (int i = 0; i < d; i++) {
                    term += wMod;
                    if (term >= MOD)
                        term -= MOD;
                }

                ans += term;
                if (ans >= MOD)
                    ans -= MOD;
            }
        }

        long contribution = (total - b - 1) % MOD;
        if (contribution < 0)
            contribution += MOD;
        ans = (ans + contribution) % MOD;
        return ans;
    }

    static String solve() {
        int kFrom = 2;
        int kTo = 300;

        long ans = IntStream.rangeClosed(kFrom, kTo)
                .parallel()
                .mapToLong(k -> {
                    int b = 6 * k + 3;
                    return sOfBase(b);
                })
                .reduce(0, (a, b) -> (a + b) % MOD);

        return Long.toString(ans);
    }

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