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

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. The 20-State Automaton

To decide whether another digit may be appended, only two pieces of information matter: the last digit $d\in\{0,\dots,9\}$ and whether the current run length is $1$ or $2$. This gives $20$ states $s=(d,r)$.

From $(d,1)$, every next digit is legal. If we append $d$, the next state is $(d,2)$; if we append $x\ne d$, the next state is $(x,1)$. From $(d,2)$, the digit $d$ is forbidden, while each $x\ne d$ leads to $(x,1)$. This is exactly what build_automaton constructs in all three language versions.

3. Tail Moments and the Linear System

Fix a state $s$, and let $W_s$ be the set of all finite valid tails that may follow $s$, including the empty tail $\epsilon$. For $m\ge 0$, define

$$H_s^{(m)}=\sum_{w\in W_s}\frac{\operatorname{val}(w)^m}{10^{(m+1)|w|}},$$

with the convention that the empty tail contributes $1$ when $m=0$ and $0$ otherwise. If a non-empty tail is written as $w=xu$, where the first appended digit is $x$ and the remaining tail $u$ starts from the next state $t$, then

$$\operatorname{val}(xu)=x\cdot 10^{|u|}+\operatorname{val}(u).$$

Applying the binomial theorem gives the recurrence

$$H_s^{(m)}=\delta_{m,0}+\frac{1}{10^{m+1}}\sum_{(x,t)\in T(s)}\sum_{j=0}^{m}\binom{m}{j}x^{m-j}H_t^{(j)},$$

where $T(s)$ is the set of labeled outgoing transitions from state $s$. The term $j=m$ involves the unknown moments of the same order, so we move it to the left:

$$H_s^{(m)}-\frac{1}{10^{m+1}}\sum_{(x,t)\in T(s)}H_t^{(m)} =\delta_{m,0}+\frac{1}{10^{m+1}}\sum_{(x,t)\in T(s)}\sum_{j=0}^{m-1}\binom{m}{j}x^{m-j}H_t^{(j)}.$$

For each fixed $m$, this is a $20\times20$ linear system. Because the right-hand side uses only moments of order $\lt m$, the code solves the systems in increasing order of $m$, precomputing binomial coefficients and digit powers and then using Gaussian elimination with partial pivoting.

4. From Moments to the Harmonic Sum

Choose a valid decimal prefix $p$ and let $s$ be the state determined by its last digit and current run length. Every admissible integer that begins with $p$ can be written as

$$p\cdot 10^{|w|}+\operatorname{val}(w),\qquad w\in W_s.$$

Its entire contribution is therefore

$$R(p,s)=\sum_{w\in W_s}\frac{1}{p\cdot 10^{|w|}+\operatorname{val}(w)}.$$

Since $0\le \operatorname{val}(w) \lt 10^{|w|}$, we have $\operatorname{val}(w)/(p10^{|w|}) \lt 1/p$. For the actual computation the code uses prefixes with at least two or three digits, so this ratio is comfortably below $1$. Hence

$$\frac{1}{p\cdot 10^{|w|}+\operatorname{val}(w)} =\sum_{m\ge 0}(-1)^m\frac{\operatorname{val}(w)^m}{p^{m+1}10^{(m+1)|w|}}.$$

Summing over all valid tails gives the key identity

$$R(p,s)=\sum_{m\ge 0}(-1)^m\frac{H_s^{(m)}}{p^{m+1}}.$$

This is the formula used by estimate_series_value. Small admissible integers are summed directly, and the rest are grouped by prefix; each prefix contributes an alternating series whose coefficients are the precomputed moments.

5. Code-Level Checkpoints

The C++ implementation verifies two facts before printing the final value. First, among $1\le n\le 1200$, exactly $20$ numbers contain three equal consecutive digits. Second, splitting the same infinite sum with prefix_digits = 2 and with prefix_digits = 3 gives answers that differ by less than $10^{-15}$. These checks confirm both the automaton and the prefix-tail decomposition.

How the Code Works

has_three_equal_consecutive tests the forbidden pattern directly. build_automaton creates the 20 labeled states. solve_moments computes $H_s^{(m)}$ for $0\le m\le 16$ using high-precision arithmetic: cpp_dec_float_100 in C++, Decimal in Python, and BigDecimal in Java. Finally, estimate_series_value adds exact reciprocals below $10^{D-1}$, scans the admissible $D$-digit prefixes, determines their terminal state, and accumulates the truncated alternating expansion for $R(p,s)$. The three solutions differ only in syntax, not in mathematics.

Complexity Analysis

Let $M$ be the largest moment index, with $M=16$ in the repository. For each $m$, one dense $20\times20$ linear system is solved, so the dominant algebraic work is $O(M\cdot 20^3)$. The prefix phase scans $9\cdot 10^{D-1}$ candidate $D$-digit prefixes, which is tiny for $D=3$. Memory usage is $O(M\cdot 20)$ for the table of moments, plus $O(20^2)$ temporary storage for elimination.

Footnotes and References

Problem page: https://projecteuler.net/problem=368
Kempner series: https://en.wikipedia.org/wiki/Kempner_series
Finite-state machine: https://en.wikipedia.org/wiki/Finite-state_machine
Binomial theorem: https://en.wikipedia.org/wiki/Binomial_theorem
Gaussian elimination: https://en.wikipedia.org/wiki/Gaussian_elimination

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