Problem 368: A Kempner-like Series

Project Euler Problem 368 Solution

EulerSolve provides an optimized solution for Project Euler Problem 368, A Kempner-like Series, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $\mathcal A$ be the set of positive integers whose decimal expansion never contains a block $ddd$ of three equal consecutive digits. The required value is the convergent Kempner-like series $$S=\sum_{n\in\mathcal A}\frac{1}{n}.$$ The local C++, Python, and Java programs all implement the same method and print $S \approx 253.6135092068$. The real task is therefore not brute force, but how to sum the infinite tail of admissible numbers in a controlled way. Mathematical Approach 1. Why the Series Converges The forbidden pattern can be counted with the same state machine used by the code. Let $a_k$ be the number of length-$k$ digit strings with no three equal consecutive digits, allowing leading zeroes for counting purposes. Split them into strings ending with run length $1$ and run length $2$, say $u_k$ and $v_k$. Appending a different digit gives $u_{k+1}=9(u_k+v_k)$, while appending the same digit is possible only from run length $1$, so $v_{k+1}=u_k$. Hence $$a_{k+1}=9a_k+9a_{k-1}.$$ The dominant root of $x^2-9x-9=0$ is $$\rho=\frac{9+\sqrt{117}}{2}\lt 10.$$ Therefore the number of admissible $k$-digit denominators is $O(\rho^k)$, so the total contribution of all $k$-digit terms is $O((\rho/10)^k)$. Since $\rho/10\lt 1$, the series converges exponentially fast by digit length. 2....

Detailed mathematical approach

Problem Summary

The harmonic series $1+\frac12+\frac13+\cdots$ diverges. Dropping every term whose denominator contains the digit $9$ leaves the Kempner series, which converges to about $22.9206766193$. Here $1/n$ is dropped instead whenever the decimal expansion of $n$ contains three or more equal consecutive digits. Writing $\mathcal A$ for the surviving positive integers, we must evaluate

$$S=\sum_{n\in\mathcal A}\frac1n$$

and give the result rounded to $10$ digits after the decimal point. The statement supplies one checkpoint: of the first $1200$ terms exactly $20$ are omitted — the reciprocals of $111,222,\dots,999$, then $1000$, $1110$ and $1111$ through $1119$. ^[1]

The whole solution in one sentence: a twenty-state machine remembers only the last digit and whether it has just been repeated, the surviving numbers are grouped by their leading three digits, and the infinite tail behind each group is replaced by a short alternating series whose coefficients come from seventeen linear systems of twenty unknowns.

Mathematical Approach

Two obstacles stand in the way: the series converges far too slowly for direct summation, and $\mathcal A$ is defined by digits rather than by arithmetic. The second cures the first, because a digit condition can be tracked by a finite machine, and a finite machine turns the infinite tail into finite linear algebra.

1. What the rule forbids

A run of four or more equal digits contains a block of three equal digits, so forbidding "three or more equal consecutive digits" is the same as forbidding any block $ddd$; a direct test confirms the two agree for every $n\lt300\,000$. A member of $\mathcal A$ therefore never carries a run longer than two, so reading from the left we need never remember a run length above $2$. Only runs matter: $1212$ and $9119$ are in $\mathcal A$, $1222$ and $45551$ are not.

2. How many numbers survive, and why the sum converges

Let $A_d$ count the $d$-digit members of $\mathcal A$ (no leading zero), split by the run that ends them: $U_d$ ending in a run of one digit, $V_d$ in a run of two, so $A_d=U_d+V_d$. Appending a digit different from the last is always allowed, in $9$ ways, and starts a new run; appending a copy is allowed only from a run of length one. Hence $U_{d+1}=9A_d$ and $V_{d+1}=U_d$, and since $U_d=9A_{d-1}$ for $d\ge2$,

$$A_{d+1}=9A_d+9A_{d-1}\qquad(d\ge2),$$

with $A_1=9$ and $A_2=90$, giving $891$, $8829$, $87480$, $866781$; the $891$ three-digit members are the $900$ three-digit numbers minus $111,\dots,999$. The sequence is catalogued with this interpretation and recurrence. ^[3]

The characteristic equation $x^2=9x+9$ has dominant root

$$\rho=\frac{9+\sqrt{117}}2=9.9083269132\ldots,$$

so $A_d$ grows like $\rho^{\,d}$, and the ratios $A_{d+1}/A_d$ settle on $9.908327$. Every $d$-digit number is at least $10^{d-1}$, so the contribution $c_d$ of the $d$-digit members obeys

$$c_d\le\frac{A_d}{10^{\,d-1}}=O\!\left(\left(\frac{\rho}{10}\right)^{\!d}\right).$$

Because $\rho\lt10$ the bounds form a convergent geometric series, so $S$ is finite. ^[6] This is Kempner's argument of 1914, which grouped the no-nine numbers by digit length with ratio $9/10$. ^[2] Here the ratio is $\rho/10=0.9908326913\ldots$, and that one number is the whole difficulty.

3. Why summation alone cannot work

These $c_d$, and the partial sums over all members below $10^{d}$, are rounded to nine decimals:

$d=1$: $c_1=2.828968254$, partial sum $2.828968254$;
$d=2$: $c_2=2.348409264$, partial sum $5.177377518$;
$d=3$: $c_3=2.281607142$, partial sum $7.458984660$;
$d=4$: $c_4=2.257036275$, partial sum $9.716020935$;
$d=5,6$: $c_5=2.235924310$, $c_6=2.215390644$;
$d=7$: $c_7=2.195077387$, partial sum $16.362413276$.

The observed ratios $c_{d+1}/c_d$ are $0.830129$, $0.971554$, $0.989231$, $0.990646$, $0.990816$, $0.990831$, visibly approaching $q=\rho/10$. Section 2 proves convergence, but an upper bound alone does not prove this ratio tends to $q$; the figures only suggest it. Extrapolating from $c_7$, the tail beyond length $d$ is roughly $c_7q^{\,d-6}/(1-q)$, which drops below $5\cdot10^{-11}$ — the accuracy ten decimals demand — only near $d=3177$. That estimate is no certified cutoff, but it shows no refinement of direct summation closes the gap: the tail must be summed in closed form, with the rigorous bound of Section 10.

4. A twenty-state machine

Read the digits from left to right. To decide whether the next digit is legal, two facts suffice: the last digit $d$, and whether the run ending the string has length $1$ or $2$. These $20$ states $s=(d,r)$ form a deterministic finite automaton, stored as a partial one — an illegal digit simply has no transition, so no rejecting state is kept. ^[4] The transitions are

from $(d,1)$: the digit $d$ leads to $(d,2)$, and each $x\ne d$ leads to $(x,1)$ — ten moves;
from $(d,2)$: the digit $d$ is rejected, and each $x\ne d$ leads to $(x,1)$ — nine moves.

The crucial consequence: for an admissible prefix the state depends only on its last two digits, never on its length or leading digits. Those earlier digits are not ignored; they decide admissibility, and $123$ and $11123$ end alike but only the first survives. All $20$ states occur: of the $891$ admissible three-digit prefixes, $81$ end in two equal digits and land in a run-$2$ state, the other $810$ in a run-$1$ state.

5. Grouping the series by prefix

Fix a prefix length $D$. Every member of $\mathcal A$ with at least $D$ digits has a unique leading $D$-digit block $p$, itself admissible, whose remaining digits form a string $w$ of length $j\ge0$ spelled out by the automaton from $s(p)$. With $v(w)$ the value of $w$ read as a $j$-digit block, leading zeros allowed, the number is $p\cdot10^{\,j}+v(w)$. Let $W_s$ be the admissible continuations from $s$, the empty one included, and put

$$R(p)=\sum_{w\in W_{s(p)}}\frac{1}{p\cdot10^{|w|}+v(w)}.$$

Then the series splits exactly:

$$S=\sum_{\substack{n\in\mathcal A\\ n\lt10^{D-1}}}\frac1n+\sum_{p}R(p),$$

the second sum running over the admissible $D$-digit blocks. The first part is a short finite sum; all the work is in $R(p)$. Two prefixes in the same state have identical continuation sets $W_s$, so the infinite part is done once per state — $20$ times, not once per prefix — though their sums $R(p)$ still differ because $p$ does.

6. Expanding one reciprocal

Since $0\le v(w)\lt10^{\,j}$ we have $v(w)/(p\cdot10^{\,j})\lt1/p\le1$, so the geometric expansion of $1/(1+x)$ applies: ^[6]

$$\frac{1}{p\cdot10^{j}+v}=\sum_{m\ge0}\frac{(-1)^m\,v^{m}}{p^{\,m+1}\,10^{\,j(m+1)}}.$$

The factor $p^{-(m+1)}$ is the only place the prefix appears. Exchanging the two sums needs absolute convergence, and bounding each single term will not do, since infinitely many share one bound; sum the geometric series first. For fixed $w$ the absolute values of its $m$-terms total

$$\begin{aligned}&\frac{1}{p\cdot10^{|w|}}\cdot\frac{1}{1-v/(p\cdot10^{|w|})}\\ &\qquad\le\frac{1}{(p-1)10^{|w|}},\end{aligned}$$

and the continuations of length $j$ from a fixed state obey the same two-run recurrence as Section 2 — $1,10,99,981,9720,\dots$ from a run-$1$ state, $1,9,90,891,8829,\dots$ from a run-$2$ state, both $O(\rho^{\,j})$ — so those bounds are summable over $W_s$. The exchange is legitimate and yields the identity the method rests on. Define the moments of the state $s$ by

$$H^{(m)}_s=\sum_{w\in W_s}\frac{v(w)^{m}}{10^{\,|w|(m+1)}},$$

reading $v(w)^0$ as $1$ for every continuation — this matters for zero-valued ones such as $0$ and $00$, not only the empty one, whose contribution is $1$ at order $0$ and $0$ above. The name is apt: read $v(w)/10^{|w|}$ as a decimal fraction of weight $10^{-|w|}$, and $H^{(0)}_s$ sums the weights, $H^{(m)}_s$ the weighted $m$-th powers. Then

$$R(p)=\sum_{m\ge0}(-1)^m\,\frac{H^{(m)}_s}{p^{\,m+1}},\qquad s=s(p).$$

Recasting a Kempner-type series as a geometrically convergent one whose coefficients are moments computed by recurrence is the shape of the modern treatment of these sums. ^[9]

7. A recurrence for the moments

A non-empty continuation is a first digit $a$ followed by a continuation $w'$ from the state $t$ it leads to, and $v(aw')=a\cdot10^{|w'|}+v(w')$. Raise this to the power $m$ with the binomial theorem ^[5] and divide by $10^{(|w'|+1)(m+1)}$: in the term with exponent $u$ on $v(w')$ the powers of ten combine into $10^{-|w'|(u+1)}\cdot10^{-(m+1)}$, so summing over $w'$ gives $H^{(u)}_t$ and

$$\begin{aligned}H^{(m)}_s&=\delta_{m,0}\\ &+\frac{1}{10^{m+1}}\sum_{(a,t)}\sum_{u=0}^{m}\binom mu a^{\,m-u}H^{(u)}_t,\end{aligned}$$

where $(a,t)$ runs over the legal moves out of $s$, and $\delta_{m,0}$ is $1$ when $m=0$ and $0$ otherwise. The inner term with $u=m$ is $H^{(m)}_t$ itself, so the moments of one order are coupled: a linear system, not an explicit formula per state. Let $h_m$ be the column of the $H^{(m)}_s$ and $M$ the $20\times20$ matrix whose entry $M_{s,t}$ counts the moves from $s$ to $t$. Moving the $u=m$ terms left,

$$\left(I-\frac{M}{10^{m+1}}\right)h_m=b_m,$$

$$\begin{aligned}(b_m)_s&=\delta_{m,0}+\frac{\Sigma_m(s)}{10^{m+1}},\\ \Sigma_m(s)&=\sum_{(a,t)}\sum_{u\lt m}\binom mu a^{\,m-u}H^{(u)}_t.\end{aligned}$$

The right-hand side uses only moments of order below $m$, so the systems are solved for $m=0,1,2,\dots$ in turn. Invertibility needs care: the row sums of $M$ are $9$ and $10$, and the row-sum bound for non-negative matrices ^[7] then permits a spectral radius of exactly $10$ — not enough at $m=0$. Apply that bound to $M^2$, whose row sums count two-move paths: from a run-$1$ state, nine moves reach run-$1$ states with ten continuations each and one reaches a run-$2$ state with nine, giving $9\cdot10+9=99$; a run-$2$ state gives $9\cdot10=90$. Hence

$$\operatorname{spr}(M)^2=\operatorname{spr}(M^2)\le99\lt100,$$

so $\operatorname{spr}(M)\le\sqrt{99}=9.9499\ldots\lt10$, every eigenvalue of $M/10^{m+1}$ is below $1$ in absolute value for all $m\ge0$, and $I-M/10^{m+1}$ is invertible. The exact value is $\rho=9.9083\ldots$, the root of Section 2, but the crude bound already settles it.

8. The case $m=0$, by hand

For $m=0$ the digit $a$ drops out of the right-hand side, so $H^{(0)}_s$ depends only on the run length; put $P=H^{(0)}_{(d,1)}$ and $Q=H^{(0)}_{(d,2)}$. A run-$1$ state has one move into a run-$2$ state and nine into run-$1$ states; a run-$2$ state has nine into run-$1$ states. Hence

$$P=1+\frac{Q+9P}{10},\qquad Q=1+\frac{9P}{10}.$$

Substituting the second into the first gives $0.01P=1.1$, so $P=110$ and $Q=100$ exactly; the program's numerical solve agrees to its working precision. Note how close this came to failing: the coefficient of $P$ on the right is $0.99$, so the system is nearly singular, and indeed $H^{(0)}_s\approx1/(1-\rho/10)=109.08$. That near-singularity is the slow convergence of Section 3: the solve performs in one step a summation term-by-term addition would never finish. For $m\ge1$ the digit powers $a^{m-u}$ enter, the symmetry between digits breaks, and all $20$ unknowns differ: $H^{(1)}_{(0,1)}=54.939011049$, $H^{(1)}_{(7,1)}=54.169780280$, $H^{(1)}_{(4,2)}=49.999456049$.

9. A worked prefix

Take $p=100$. Its last two digits are equal, so $s(100)=(0,2)$ and the next digit may not be $0$ — which is why $1000$ goes while $1001,\dots,1009$ stay. The alternating series has terms

$m=0$: $100/100=1.0000000000000000$;
$m=1$: $-0.0054395060445000$;
$m=2$: $+0.0000365577267675$;
$m=3$: $-0.0000002741651468$;
$m=4$: $+0.0000000021913944$,

whose sum $0.9945967797085151$ gives $R(100)=0.9945967797$ to ten places; stopping one term earlier would have given $0.9945967775$. The same machinery yields $R(123)=0.8907188488$, $R(177)=0.5634864927$, $R(500)=0.1997827118$, $R(998)=0.1101662001$ and $R(1177)=0.0849280735$. Splitting off the empty continuation also gives a check independent of the derivation, $R(p)=1/p+\sum_x R(10p+x)$ over the admissible children of $p$; the Verification section reports how well it holds.

10. Truncating safely

Write $R_M(p)$ for the expansion stopped at order $M$. In exact arithmetic that truncation is the only approximation left, and it is easy to bound: since $v(w)^m\le10^{|w|m}$ term by term, the moments obey $0\le H^{(m)}_s\le H^{(0)}_s\le110$, so

$$\begin{aligned}|R(p)-R_M(p)|&\le110\sum_{m\gt M}p^{-(m+1)}\\ &=\frac{110}{(p-1)p^{\,M+1}}.\end{aligned}$$

With $D=3$ every prefix is at least $100$, so at $M=16$ this is at most $1.1\cdot10^{-34}$ per prefix and $9.9\cdot10^{-32}$ over all $891$ — ample for ten decimals. Against a much longer expansion the discarded corrections amount to $3.7\cdot10^{-35}$ at $M=16$ and $4.7\cdot10^{-31}$ at $M=14$. These are truncation figures only: they do not bound the decimal rounding or the skipped elimination updates described below. A larger $D$ is unnecessary; a smaller one is worse, since at $D=2$ the ratio $v/(p\cdot10^{j})$ rises towards $1/10$, so the same $M$ buys half as many digits.

How the Code Works

All three share the same four pieces, named here as in C++ and Python; Java uses the same names in camel case (hasThreeEqualConsecutive, buildAutomaton and so on). has_three_equal_consecutive decides membership of $\mathcal A$: C++ peels digits with n % 10 keeping the two previous ones, Python and Java scan the decimal string. build_automaton fills the transition table of Section 4, indexing states by state_index(digit, run_len) $=2d+(r-1)$ and omitting the illegal move out of a run-$2$ state, so each state holds ten or nine pairs (digit, to).

solve_moments solves the systems of Section 7 for $m=0,1,\dots,M$. It precomputes Pascal's triangle in binom and the powers $a^{k}$ in digit_pow, then for each $m$ assembles a as $I-M/10^{m+1}$ beside $b_m$, looping over the stored transitions and over $u\lt m$ with the moments already in h. Elimination with partial pivoting follows, ^[8] clearing the pivot column from every other row, above as well as below: this is Gauss–Jordan elimination with no back-substitution step, so the matrix becomes the identity and the last column is the solution. A guard skips a row update whose multiplier is below $10^{-70}$ in absolute value in C++, or $10^{-80}$ in Python and Java — thresholds far above the rounding of $100$-digit arithmetic, so they discard representable values and add a small approximation of their own.

estimate_series_value performs the split of Section 5 for prefix_digits $=D$: it sums $1/n$ directly, in the chosen decimal arithmetic, over the admissible $n\lt10^{D-1}$, then for each admissible $D$-digit $p$ recovers the state from last and prev — the two-digit dependence of Section 4 — and accumulates the alternating series of Section 6, multiplying inv_power by $1/p$ and flipping sign.

The three differ in arithmetic, digit scanning and built-in checks. ^[10] C++ uses Boost's cpp_dec_float_100, Python sets getcontext().prec = 100 for Decimal, and Java carries a $100$-digit MathContext with HALF_UP rounding on BigDecimal, passing mc to its products and quotients: a BigDecimal division called without one throws ArithmeticException when the exact quotient does not terminate. Only C++ self-checks, in run_checkpoints: it recounts the $20$ omissions below $1200$, then sums the series twice, at $D=2$ and $D=3$ with $M=14$, requiring agreement within $10^{-15}$; the main routine recomputes with $M=16$, and --skip-checkpoints suppresses the checks. Java also routes the rounded BigDecimal through doubleValue() before String.format: the thirteen significant digits survive a double here, but the detour is a needless risk that toPlainString would avoid.

Complexity and Verification

Let $K=20$ be the number of states and $M=16$ the highest moment, so $17$ systems are solved. Assembling the right-hand side at order $m$ costs $O(10Km)$, summing to $O(10KM^2)$, and each Gauss–Jordan solve costs $\Theta(K^3)$; the Python solver, instrumented, performs $4349$ multiply-and-subtract updates plus $230$ divisions per order, so the moment phase stays under $10^5$ high-precision operations. The prefix phase scans $9\cdot10^{D-1}$ blocks and spends $M+1$ terms on each admissible one — $891\cdot17=15147$ at $D=3$. Memory is $O(KM)$ for the moments, $O(K^2)$ for the elimination and $O(M^2)$ for Pascal's triangle: a few thousand numbers. The cost does not grow with how far out in $\mathcal A$ the series reaches, though more accuracy needs a larger $M$ and more precision. Indicative timings: C++ $0.02$ s without its checkpoints and $0.64$ s with them, Python $0.04$ s, Java $0.13$ s.

Beyond those two checks, the following were verified for this article.

Brute force below $10^7$ reproduces the statement's $20$ omissions under $1200$, the recurrence of Section 2, the counts $9,90,891,8829,87480,866781$ and every contribution and partial sum of Section 3.
An independent reimplementation in exact rational arithmetic — ordinary fractions, no rounding anywhere, prefix length $D=4$ and cutoff $M=12$ — reproduces the three programs' printed value to all ten decimals, and the three agree character for character.
The algebraic moments $H^{(0)}_s=110$ and $100$ match the numerical solve to its working precision.
The parent–child relation of Section 9 holds across all $891$ admissible three-digit prefixes, largest residual $4.9\cdot10^{-44}$ at $M=20$.
Changing $D$ moves the total only as far as Section 10 allows: at $M=16$ the splittings $D=3,4,5$ agree to $3.7\cdot10^{-35}$, while $D=2$ differs by $6.8\cdot10^{-18}$.

The strongest check uses the statement itself. Adapt the method to the no-nine Kempner series: one state with self-loops labelled $0,\dots,8$, solver dimension one, and the admissibility test and prefix assignment changed to reject any $9$. The moment recurrence is untouched and gives $H^{(0)}=10$, the expected $1/(1-9/10)$. At $D=3$, $M=28$ and $100$ digits the adaptation returns $22.92067661926415034816$ to twenty decimals, matching the published value ^[2] and the statement's $22.9206766193$ after rounding. Reproducing a known constant to twenty decimals tests the moment machinery far better than any internal check, though only the prefix tests above exercise the twenty-state restriction.

One caveat about the C++ checkpoint. At $M=14$ the two splittings differ by $8.0\cdot10^{-16}$, just inside the $10^{-15}$ tolerance; at $M=13$ the difference is $8.8\cdot10^{-15}$ and the check would fail. A margin of $1.25$ ties tolerance and cutoff together more tightly than they appear; $M=16$ there too would restore a wide margin.

Footnotes and References

The references supply the background facts used above. The twenty-state construction, the prefix decomposition, the moment identity and the error bounds are derived in this article.

Project Euler 368 — A Kempner-like Series. The official statement recalls that the harmonic series diverges and that omitting every denominator containing a 9 leaves a series converging to about 22.9206766193, then asks for the value of the series that omits every denominator with three or more equal consecutive digits, rounded to 10 decimal places. Only 20 of its first 1200 terms are omitted, from 1/111 to 1/1119.
Wikipedia — Kempner series. Gives the no-nine sum as 22.92067661926415034816 to twenty decimals, the value cross-checked above, and Kempner's 1914 proof: group the admissible numbers by digit length, note that a d-digit number is at least 10^(d−1), bound the groups geometrically with ratio 9/10 — the argument adapted in Section 2. Both proof and method, it adds, extend to a whole forbidden string of digits.
OEIS A057092. The sequence 1, 9, 90, 891, 8829, 87480, 866781 satisfies a(n) = 9(a(n−1) + a(n−2)), and a comment records that a(n) counts the n-digit integers with no digit repeated three times in a row, with a(3) = 891. These are the A_d of Section 2.
Wikipedia — Deterministic finite automaton. Defines a DFA by its states, input alphabet, transition function, start state and accepting states, moving deterministically along the matching transition as each symbol is read. Section 4 tracks the digit condition this way, with the partial-transition convention the source describes.
Wikipedia — Binomial theorem. For a non-negative integer n, (x + y)^n is the sum over k of C(n, k) x^(n−k) y^k. Section 7 applies it to the m-th power of a·10^|w'| + v(w').
Wikipedia — Geometric series. A geometric series with ratio r converges absolutely to a/(1 − r) exactly when |r| < 1, with terms alternating in sign for negative r. Section 2 uses the criterion for the bound on c_d, Section 6 the expansion with r = −v/(p·10^j).
Wikipedia — Perron–Frobenius theorem. A non-negative square matrix has an eigenvalue, non-negative and at least as large in absolute value as every other, bounded above by its largest row sum. The row sums 9 and 10 of M give only the non-strict bound 10, so Section 7 applies the theorem to M², with row sums 99 and 90, for the strict bound invertibility needs.
Wikipedia — Gaussian elimination. Describes row reduction, its Gauss–Jordan variant, and partial pivoting, which picks the pivot of largest absolute value for numerical stability. Its count of about 2n³/3 operations covers forward elimination plus back substitution; solve_moments reduces all the way to the identity, so its cost is the Θ(n³) of Gauss–Jordan.
J.-F. Burnol, Moments in the exact summation of the curious series of Kempner type (arXiv:2402.08525). Presents Kempner series as geometrically convergent series whose coefficients involve power sums over the allowed digits and the moments of a measure on the unit interval, obtained by recurrence and framed by a priori bounds. It treats forbidden single digits; Sections 6 and 7 derive a state-dependent version for blocks.
Boost.Multiprecision — cpp_dec_float documents the decimal type behind cpp_dec_float_100, Python — decimal the settable precision of the arithmetic context, and Java — BigDecimal that a division without a MathContext throws when the exact quotient does not terminate.

Problem 368 source code

C++

#include <array>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

using Real = boost::multiprecision::cpp_dec_float_100;

constexpr int kStateCount = 20;  // 10 digits * run length {1,2}

int state_index(const int digit, const int run_len) {
    return digit * 2 + (run_len - 1);
}

bool has_three_equal_consecutive(std::uint64_t n) {
    int prev1 = -1;
    int prev2 = -2;
    while (n > 0) {
        const int d = static_cast<int>(n % 10ULL);
        n /= 10ULL;
        if (d == prev1 && d == prev2) {
            return true;
        }
        prev2 = prev1;
        prev1 = d;
    }
    return false;
}

struct Automaton {
    std::array<std::vector<std::pair<int, int>>, kStateCount> trans;
};

Automaton build_automaton() {
    Automaton a;
    for (int d = 0; d <= 9; ++d) {
        for (int r = 1; r <= 2; ++r) {
            const int s = state_index(d, r);
            a.trans[static_cast<std::size_t>(s)].clear();
            for (int x = 0; x <= 9; ++x) {
                if (x == d) {
                    if (r == 2) {
                        continue;
                    }
                    a.trans[static_cast<std::size_t>(s)].push_back({x, state_index(d, 2)});
                } else {
                    a.trans[static_cast<std::size_t>(s)].push_back({x, state_index(x, 1)});
                }
            }
        }
    }
    return a;
}

std::vector<std::vector<Real>> solve_moments(const Automaton& automaton, const int max_m) {
    std::vector<std::vector<Real>> h(static_cast<std::size_t>(kStateCount),
                                     std::vector<Real>(static_cast<std::size_t>(max_m + 1), Real(0)));

    // Binomial coefficients.
    std::vector<std::vector<std::uint64_t>> binom(static_cast<std::size_t>(max_m + 1),
                                                  std::vector<std::uint64_t>(static_cast<std::size_t>(max_m + 1), 0));
    for (int n = 0; n <= max_m; ++n) {
        binom[static_cast<std::size_t>(n)][0] = 1;
        binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(n)] = 1;
        for (int k = 1; k < n; ++k) {
            binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] =
                binom[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k - 1)] +
                binom[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k)];
        }
    }

    std::array<std::array<Real, 32>, 10> digit_pow{};
    for (int d = 0; d <= 9; ++d) {
        digit_pow[static_cast<std::size_t>(d)][0] = Real(1);
        for (int p = 1; p < 32; ++p) {
            digit_pow[static_cast<std::size_t>(d)][static_cast<std::size_t>(p)] =
                digit_pow[static_cast<std::size_t>(d)][static_cast<std::size_t>(p - 1)] * Real(d);
        }
    }

    for (int m = 0; m <= max_m; ++m) {
        const Real c = pow(Real(10), -(m + 1));

        std::vector<std::vector<Real>> a(static_cast<std::size_t>(kStateCount),
                                         std::vector<Real>(static_cast<std::size_t>(kStateCount + 1), Real(0)));

        for (int s = 0; s < kStateCount; ++s) {
            a[static_cast<std::size_t>(s)][static_cast<std::size_t>(s)] = Real(1);

            for (const auto& [digit, to] : automaton.trans[static_cast<std::size_t>(s)]) {
                (void)digit;
                a[static_cast<std::size_t>(s)][static_cast<std::size_t>(to)] -= c;
            }

            Real rhs = (m == 0) ? Real(1) : Real(0);
            for (const auto& [digit, to] : automaton.trans[static_cast<std::size_t>(s)]) {
                for (int u = 0; u < m; ++u) {
                    rhs += c * Real(binom[static_cast<std::size_t>(m)][static_cast<std::size_t>(u)]) *
                           digit_pow[static_cast<std::size_t>(digit)][static_cast<std::size_t>(m - u)] *
                           h[static_cast<std::size_t>(to)][static_cast<std::size_t>(u)];
                }
            }
            a[static_cast<std::size_t>(s)][static_cast<std::size_t>(kStateCount)] = rhs;
        }

        // Gaussian elimination with partial pivoting.
        for (int col = 0; col < kStateCount; ++col) {
            int pivot = col;
            for (int r = col + 1; r < kStateCount; ++r) {
                if (abs(a[static_cast<std::size_t>(r)][static_cast<std::size_t>(col)]) >
                    abs(a[static_cast<std::size_t>(pivot)][static_cast<std::size_t>(col)])) {
                    pivot = r;
                }
            }
            std::swap(a[static_cast<std::size_t>(col)], a[static_cast<std::size_t>(pivot)]);

            const Real pv = a[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)];
            for (int j = col; j <= kStateCount; ++j) {
                a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)] /= pv;
            }

            for (int r = 0; r < kStateCount; ++r) {
                if (r == col) {
                    continue;
                }
                const Real f = a[static_cast<std::size_t>(r)][static_cast<std::size_t>(col)];
                if (abs(f) < Real("1e-70")) {
                    continue;
                }
                for (int j = col; j <= kStateCount; ++j) {
                    a[static_cast<std::size_t>(r)][static_cast<std::size_t>(j)] -=
                        f * a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)];
                }
            }
        }

        for (int s = 0; s < kStateCount; ++s) {
            h[static_cast<std::size_t>(s)][static_cast<std::size_t>(m)] =
                a[static_cast<std::size_t>(s)][static_cast<std::size_t>(kStateCount)];
        }
    }

    return h;
}

Real estimate_series_value(const std::vector<std::vector<Real>>& h, const int max_m, const int prefix_digits) {
    const std::uint64_t low = 1;
    std::uint64_t high_small = 1;
    for (int i = 0; i < prefix_digits - 1; ++i) {
        high_small *= 10ULL;
    }
    const std::uint64_t high_prefix = high_small * 10ULL;

    Real total = 0;

    for (std::uint64_t n = low; n < high_small; ++n) {
        if (!has_three_equal_consecutive(n)) {
            total += Real(1) / Real(n);
        }
    }

    for (std::uint64_t p = high_small; p < high_prefix; ++p) {
        if (has_three_equal_consecutive(p)) {
            continue;
        }

        const int last = static_cast<int>(p % 10ULL);
        const int prev = static_cast<int>((p / 10ULL) % 10ULL);
        const int run = (last == prev) ? 2 : 1;
        const int s = state_index(last, run);

        const Real inv_p = Real(1) / Real(p);
        Real inv_power = inv_p;
        Real contrib = 0;
        int sign = 1;
        for (int m = 0; m <= max_m; ++m) {
            const Real term = h[static_cast<std::size_t>(s)][static_cast<std::size_t>(m)] * inv_power;
            contrib += (sign > 0) ? term : -term;
            inv_power *= inv_p;
            sign = -sign;
        }
        total += contrib;
    }

    return total;
}

bool run_checkpoints() {
    int omitted = 0;
    for (int n = 1; n <= 1200; ++n) {
        if (has_three_equal_consecutive(static_cast<std::uint64_t>(n))) {
            ++omitted;
        }
    }
    if (omitted != 20) {
        std::cerr << "Checkpoint failed: omitted terms up to 1200\n";
        return false;
    }

    const Automaton automaton = build_automaton();
    const int max_m = 14;
    const auto h = solve_moments(automaton, max_m);

    const Real s2 = estimate_series_value(h, max_m, 2);
    const Real s3 = estimate_series_value(h, max_m, 3);
    if (abs(s2 - s3) > Real("1e-15")) {
        std::cerr << "Checkpoint failed: prefix split mismatch\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const Automaton automaton = build_automaton();
    const int max_m = 16;
    const auto h = solve_moments(automaton, max_m);
    const Real value = estimate_series_value(h, max_m, 3);

    std::cout << std::fixed << std::setprecision(10) << value << '\n';
    return 0;
}

Python

from decimal import Decimal, getcontext

getcontext().prec = 100

kStateCount = 20

def state_index(digit, run_len):
    return digit * 2 + (run_len - 1)

def has_three_equal_consecutive(n):
    s = str(n)
    for i in range(len(s) - 2):
        if s[i] == s[i+1] == s[i+2]: return True
    return False

def build_automaton():
    trans = [[] for _ in range(kStateCount)]
    for d in range(10):
        for r in (1, 2):
            s = state_index(d, r)
            for x in range(10):
                if x == d:
                    if r == 2: continue
                    trans[s].append((x, state_index(d, 2)))
                else:
                    trans[s].append((x, state_index(x, 1)))
    return trans

def solve_moments(trans, max_m):
    h = [[Decimal(0)] * (max_m + 1) for _ in range(kStateCount)]
    
    binom = [[0] * (max_m + 1) for _ in range(max_m + 1)]
    for n in range(max_m + 1):
        binom[n][0] = 1
        binom[n][n] = 1
        for k in range(1, n):
            binom[n][k] = binom[n-1][k-1] + binom[n-1][k]
            
    digit_pow = [[Decimal(1)] * 32 for _ in range(10)]
    for d in range(10):
        for p in range(1, 32):
            digit_pow[d][p] = digit_pow[d][p - 1] * Decimal(d)
            
    c10 = Decimal(10)
            
    for m in range(max_m + 1):
        c = Decimal(1) / (c10 ** (m + 1))
        
        a = [[Decimal(0)] * (kStateCount + 1) for _ in range(kStateCount)]
        
        for s in range(kStateCount):
            a[s][s] = Decimal(1)
            for digit, to in trans[s]:
                a[s][to] -= c
                
            rhs = Decimal(1) if m == 0 else Decimal(0)
            for digit, to in trans[s]:
                for u in range(m):
                    term = c * Decimal(binom[m][u]) * digit_pow[digit][m - u] * h[to][u]
                    rhs += term
            a[s][kStateCount] = rhs
            
        for col in range(kStateCount):
            pivot = col
            for r in range(col + 1, kStateCount):
                if abs(a[r][col]) > abs(a[pivot][col]):
                    pivot = r
            a[col], a[pivot] = a[pivot], a[col]
            
            pv = a[col][col]
            for j in range(col, kStateCount + 1):
                a[col][j] /= pv
                
            for r in range(kStateCount):
                if r == col: continue
                factor = a[r][col]
                if abs(factor) < Decimal("1e-80"): continue
                for j in range(col, kStateCount + 1):
                    a[r][j] -= factor * a[col][j]
                    
        for s in range(kStateCount):
            h[s][m] = a[s][kStateCount]
            
    return h

def estimate_series_value(h, max_m, prefix_digits):
    low = 1
    high_small = 10 ** (prefix_digits - 1)
    high_prefix = high_small * 10
    
    total = Decimal(0)
    for n in range(low, high_small):
        if not has_three_equal_consecutive(n):
            total += Decimal(1) / Decimal(n)
            
    for p in range(high_small, high_prefix):
        if has_three_equal_consecutive(p): continue
        
        last = p % 10
        prev = (p // 10) % 10
        run = 2 if last == prev else 1
        s = state_index(last, run)
        
        inv_p = Decimal(1) / Decimal(p)
        inv_power = inv_p
        contrib = Decimal(0)
        sign = 1
        for m in range(max_m + 1):
            term = h[s][m] * inv_power
            if sign > 0: contrib += term
            else: contrib -= term
            inv_power *= inv_p
            sign = -sign
            
        total += contrib
        
    return total

def solve():
    trans = build_automaton()
    max_m = 16
    h = solve_moments(trans, max_m)
    val = estimate_series_value(h, max_m, 3)
    return "{:.10f}".format(val)

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

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;

public class Euler368 {
    static final int kStateCount = 20;
    static final MathContext mc = new MathContext(100, RoundingMode.HALF_UP);

    static int stateIndex(int digit, int runLen) {
        return digit * 2 + (runLen - 1);
    }

    static boolean hasThreeEqualConsecutive(long n) {
        String s = Long.toString(n);
        for (int i = 0; i < s.length() - 2; i++) {
            if (s.charAt(i) == s.charAt(i + 1) && s.charAt(i) == s.charAt(i + 2)) {
                return true;
            }
        }
        return false;
    }

    static class Transition {
        int digit;
        int to;

        Transition(int digit, int to) {
            this.digit = digit;
            this.to = to;
        }
    }

    static List<List<Transition>> buildAutomaton() {
        List<List<Transition>> trans = new ArrayList<>();
        for (int i = 0; i < kStateCount; i++)
            trans.add(new ArrayList<>());

        for (int d = 0; d <= 9; d++) {
            for (int r = 1; r <= 2; r++) {
                int s = stateIndex(d, r);
                for (int x = 0; x <= 9; x++) {
                    if (x == d) {
                        if (r == 2)
                            continue;
                        trans.get(s).add(new Transition(x, stateIndex(d, 2)));
                    } else {
                        trans.get(s).add(new Transition(x, stateIndex(x, 1)));
                    }
                }
            }
        }
        return trans;
    }

    static BigDecimal[][] solveMoments(List<List<Transition>> trans, int maxM) {
        BigDecimal[][] h = new BigDecimal[kStateCount][maxM + 1];
        for (int i = 0; i < kStateCount; i++) {
            for (int j = 0; j <= maxM; j++)
                h[i][j] = BigDecimal.ZERO;
        }

        long[][] binom = new long[maxM + 1][maxM + 1];
        for (int n = 0; n <= maxM; n++) {
            binom[n][0] = 1;
            binom[n][n] = 1;
            for (int k = 1; k < n; k++) {
                binom[n][k] = binom[n - 1][k - 1] + binom[n - 1][k];
            }
        }

        BigDecimal[][] digitPow = new BigDecimal[10][32];
        for (int d = 0; d <= 9; d++) {
            digitPow[d][0] = BigDecimal.ONE;
            for (int p = 1; p < 32; p++) {
                digitPow[d][p] = digitPow[d][p - 1].multiply(BigDecimal.valueOf(d), mc);
            }
        }

        BigDecimal c10 = BigDecimal.TEN;

        for (int m = 0; m <= maxM; m++) {
            BigDecimal c = BigDecimal.ONE.divide(c10.pow(m + 1, mc), mc);

            BigDecimal[][] a = new BigDecimal[kStateCount][kStateCount + 1];
            for (int i = 0; i < kStateCount; i++) {
                for (int j = 0; j <= kStateCount; j++)
                    a[i][j] = BigDecimal.ZERO;
            }

            for (int s = 0; s < kStateCount; s++) {
                a[s][s] = BigDecimal.ONE;
                for (Transition t : trans.get(s)) {
                    a[s][t.to] = a[s][t.to].subtract(c, mc);
                }

                BigDecimal rhs = (m == 0) ? BigDecimal.ONE : BigDecimal.ZERO;
                for (Transition t : trans.get(s)) {
                    for (int u = 0; u < m; u++) {
                        BigDecimal term = c.multiply(BigDecimal.valueOf(binom[m][u]), mc)
                                .multiply(digitPow[t.digit][m - u], mc)
                                .multiply(h[t.to][u], mc);
                        rhs = rhs.add(term, mc);
                    }
                }
                a[s][kStateCount] = rhs;
            }

            for (int col = 0; col < kStateCount; col++) {
                int pivot = col;
                for (int r = col + 1; r < kStateCount; r++) {
                    if (a[r][col].abs().compareTo(a[pivot][col].abs()) > 0) {
                        pivot = r;
                    }
                }

                BigDecimal[] temp = a[col];
                a[col] = a[pivot];
                a[pivot] = temp;

                BigDecimal pv = a[col][col];
                for (int j = col; j <= kStateCount; j++) {
                    a[col][j] = a[col][j].divide(pv, mc);
                }

                BigDecimal epsilon = new BigDecimal("1e-80");
                for (int r = 0; r < kStateCount; r++) {
                    if (r == col)
                        continue;
                    BigDecimal factor = a[r][col];
                    if (factor.abs().compareTo(epsilon) < 0)
                        continue;
                    for (int j = col; j <= kStateCount; j++) {
                        a[r][j] = a[r][j].subtract(factor.multiply(a[col][j], mc), mc);
                    }
                }
            }

            for (int s = 0; s < kStateCount; s++) {
                h[s][m] = a[s][kStateCount];
            }
        }
        return h;
    }

    static BigDecimal estimateSeriesValue(BigDecimal[][] h, int maxM, int prefixDigits) {
        long low = 1;
        long highSmall = 1;
        for (int i = 0; i < prefixDigits - 1; i++)
            highSmall *= 10;
        long highPrefix = highSmall * 10;

        BigDecimal total = BigDecimal.ZERO;

        for (long n = low; n < highSmall; n++) {
            if (!hasThreeEqualConsecutive(n)) {
                total = total.add(BigDecimal.ONE.divide(BigDecimal.valueOf(n), mc), mc);
            }
        }

        for (long p = highSmall; p < highPrefix; p++) {
            if (hasThreeEqualConsecutive(p))
                continue;

            int last = (int) (p % 10);
            int prev = (int) ((p / 10) % 10);
            int run = (last == prev) ? 2 : 1;
            int s = stateIndex(last, run);

            BigDecimal invP = BigDecimal.ONE.divide(BigDecimal.valueOf(p), mc);
            BigDecimal invPower = invP;
            BigDecimal contrib = BigDecimal.ZERO;
            int sign = 1;

            for (int m = 0; m <= maxM; m++) {
                BigDecimal term = h[s][m].multiply(invPower, mc);
                if (sign > 0)
                    contrib = contrib.add(term, mc);
                else
                    contrib = contrib.subtract(term, mc);
                invPower = invPower.multiply(invP, mc);
                sign = -sign;
            }
            total = total.add(contrib, mc);
        }
        return total;
    }

    static String solve() {
        List<List<Transition>> trans = buildAutomaton();
        int maxM = 16;
        BigDecimal[][] h = solveMoments(trans, maxM);
        BigDecimal val = estimateSeriesValue(h, maxM, 3);
        return String.format(java.util.Locale.US, "%.10f", val.setScale(10, RoundingMode.HALF_UP).doubleValue());
    }

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