Detailed mathematical approach

Problem Summary

A kangaroo starts at 0, and each hop has an independent length uniformly distributed on \([0,1]\). Let \(H(x)\) denote the expected number of hops needed to reach or pass position \(x\). The task is to determine the first eight digits after the decimal point of \(H(10^6)\) that are not equal to 6.

A direct numerical attack on the defining recurrence is the wrong scale of computation. The important fact is that for large integers \(n\), the fractional part of \(H(n)\) is extraordinarily close to \(0.6666\ldots\), so the real problem is to identify where that long run of sixes first changes and what the first non-6 digits are.

Mathematical Approach

The solution is built from an exact renewal equation, a Laplace transform with a very structured denominator, and a residue calculation at the dominant complex poles.

Conditioning on the first hop

If \(x \le 0\), no hop is required, so \(H(x)=0\). For \(x>0\), after the first hop of length \(U\sim \mathrm{Unif}[0,1]\), the remaining expected number of hops is \(H(x-U)\). Therefore

$H(x)=1+\int_0^1 H(x-u)\,du,$

with the convention \(H(y)=0\) for \(y\le 0\). This is the exact problem-specific object behind all three implementations.

Delay differential equation and exact small values

Differentiating the integral equation gives

$H'(x)=H(x)-H(x-1)\qquad (x>0).$

On \(0<x<1\), the term \(H(x-1)\) vanishes, so \(H'(x)=H(x)\) and \(H(0^+)=1\). Hence

$H(x)=e^x\qquad (0\le x\le 1).$

On \(1<x<2\), this becomes \(H'(x)=H(x)-e^{x-1}\), so

$H(x)=e^x-(x-1)e^{x-1}\qquad (1\le x\le 2).$

In particular,

$H(2)=e^2-e.$

Repeating the same piecewise argument on the next interval yields

$H(3)=e^3-2e^2+\frac{e}{2}.$

These are not ornamental formulas: they are concrete checkpoints for the digit extractor and they show the characteristic exponential-polynomial structure of the exact solution.

Laplace transform and the main term \(2x+\frac23\)

Let

$F(s)=\int_0^\infty H(x)e^{-sx}\,dx.$

Transforming the renewal equation gives

$F(s)=\frac1s+\frac{1-e^{-s}}sF(s),$

hence

$F(s)=\frac{1}{s-1+e^{-s}}.$

The denominator has a double zero at \(s=0\), because

$s-1+e^{-s}=\frac{s^2}{2}-\frac{s^3}{6}+O(s^4).$

So near the origin,

$F(s)=\frac{2}{s^2}+\frac{2}{3s}+O(1).$

After inverting the Laplace transform, this becomes

$H(x)=2x+\frac23+\text{exponentially small oscillating terms}.$

For integer \(n\), the term \(2n\) is integral, so the fractional part of \(H(n)\) is governed by \(2/3\) plus a tiny correction.

The dominant nonzero poles

All remaining terms come from the nonzero roots of

$s-1+e^{-s}=0.$

Rewriting it as \((s-1)e^s=-1\) gives

$s=1+W_k(-e^{-1}),$

where \(W_k\) is a branch of the Lambert \(W\) function. The nearest nonzero roots form a complex-conjugate pair

$\lambda=\alpha+i\beta,\qquad \bar{\lambda}=\alpha-i\beta,$

with

$\alpha\approx -2.08884301561304,\qquad \beta\approx 7.46148928565425.$

If \(p\neq 0\) is a root, then the residue of \(F(s)\) at \(p\) is \(1/p\), because the derivative of the denominator is \(1-e^{-p}=p\) at a root. Therefore the dominant correction at integer \(n\) is

$\varepsilon_n \approx \frac{e^{\lambda n}}{\lambda}+\frac{e^{\bar{\lambda} n}}{\bar{\lambda}}=\frac{2\big(\alpha\cos(\beta n)+\beta\sin(\beta n)\big)}{\alpha^2+\beta^2}e^{\alpha n}.$

Hence

$H(n)=2n+\frac23+\varepsilon_n+\text{smaller terms},$

and the correction decays exponentially because \(\alpha<0\).

Locating the first changed decimal digit

The implementations never expand the full decimal prefix. Instead they compute

$\log_{10}|\varepsilon_n|=\frac{\alpha n}{\ln 10}+\log_{10}\left|\frac{2\big(\alpha\cos(\beta n)+\beta\sin(\beta n)\big)}{\alpha^2+\beta^2}\right|.$

If

$-\log_{10}|\varepsilon_n|=m+\theta,\qquad m=\lfloor -\log_{10}|\varepsilon_n|\rfloor,\quad 0\le \theta<1,$

then

$\varepsilon_n=\sigma 10^{-m-\theta},\qquad \sigma\in\{-1,1\}.$

So the first deviation from \(0.6666\ldots\) occurs around the \(m\)-th decimal position. The code sets \(t=\sigma 10^{-\theta}\), writes

$\frac23+t=q+r,\qquad q=\left\lfloor \frac23+t\right\rfloor,\quad 0\le r<1,$

uses \(q\) as a carry or borrow into a short artificial tail of sixes, and then reads the next digits from the fractional part \(r\). Any digit equal to 6 is skipped, and the first eight remaining digits are the answer.

How the Code Works

The C++, Python, and Java implementations all evaluate the same asymptotic formula at high precision. They store \(\alpha\) and \(\beta\), reduce \(\beta n\) modulo \(2\pi\), and compute the trigonometric envelope that multiplies \(e^{\alpha n}\).

The important numerical trick is that they work with \(\log_{10}|\varepsilon_n|\) rather than with the full tiny number \(\varepsilon_n\). That immediately reveals how far to the right the first non-6 digit appears, without materializing the enormous block of preceding sixes.

After the location step, the implementation creates a short decimal tail initially filled with sixes, applies the possible carry or borrow coming from \(q\), and then continues with the fractional digits of \(r\). Every time a produced digit is 6 it is ignored; every other digit is appended. One implementation also checks the extractor against the exact closed forms for \(H(2)\) and \(H(3)\).

Complexity Analysis

The problem is solved by a constant-size analytic formula, not by iterating up to \(10^6\) in the original recurrence. High-precision transcendental evaluation dominates the setup, but that cost does not grow combinatorially with \(n\).

If \(D\) denotes the number of requested non-6 digits, the extraction phase is \(O(D)\) time and \(O(D)\) space. Here \(D=8\), so the overall computation is effectively constant time and constant memory for practical purposes.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=970
2. Lambert W function: Wikipedia - Lambert W function
3. Laplace transform: Wikipedia - Laplace transform
4. Renewal theory: Wikipedia - Renewal theory
5. Delay differential equation: Wikipedia - Delay differential equation
6. Continuous uniform distribution: Wikipedia - Continuous uniform distribution