Problem 970: Kangaroo Hopping over Sixes

Project Euler Problem 970 Solution

EulerSolve provides an optimized solution for Project Euler Problem 970, Kangaroo Hopping over Sixes, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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|.$$

$$-\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

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

Problem 970 source code

C++

#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using Dec = boost::multiprecision::cpp_dec_float_100;

std::string extract_non_six_digits_from_fraction(Dec fractional_value, int wanted) {
    std::string out;
    out.reserve(static_cast<std::size_t>(wanted));

    Dec f = fractional_value;
    if (f < 0) {
        f = f - boost::multiprecision::floor(f);
    }

    for (int step = 0; static_cast<int>(out.size()) < wanted && step < 100000; ++step) {
        f *= 10;
        int digit = static_cast<int>(boost::multiprecision::floor(f + Dec("1e-70")));
        if (digit < 0) {
            digit = 0;
        }
        if (digit > 9) {
            digit = 9;
        }
        f -= Dec(digit);
        if (digit != 6) {
            out.push_back(static_cast<char>('0' + digit));
        }
    }

    return out;
}

std::string non_six_digits_large_n(std::int64_t n, int wanted) {
    // Dominant nonzero pole of 1 / (s - 1 + e^{-s}), i.e. 1 + W_1(-1/e).
    const Dec alpha("-2.0888430156130438559570867167749475005456937410367296732391125442446071101031945");
    const Dec beta("7.4614892856542545569061166121864153345090949932022092409344113914118766543223747");
    const Dec den = alpha * alpha + beta * beta;

    const Dec two_pi = 2 * boost::math::constants::pi<Dec>();
    const Dec raw = beta * Dec(n);
    const Dec theta = raw - boost::multiprecision::floor(raw / two_pi) * two_pi;
    const Dec c = boost::multiprecision::cos(theta);
    const Dec s = boost::multiprecision::sin(theta);

    const Dec envelope = 2 * (c * alpha + s * beta) / den;
    const Dec sign = (envelope >= 0 ? Dec(1) : Dec(-1));

    const Dec log10_eps =
        alpha * Dec(n) / boost::multiprecision::log(Dec(10)) +
        boost::multiprecision::log10(boost::multiprecision::abs(envelope));

    const std::int64_t shift = static_cast<std::int64_t>(
        boost::multiprecision::floor(-log10_eps).convert_to<long long>());
    const Dec frac = -log10_eps - Dec(shift);
    const Dec t = sign * boost::multiprecision::pow(Dec(10), -frac);

    const Dec base = Dec(2) / 3 + t;
    const int q = static_cast<int>(boost::multiprecision::floor(base));
    Dec r = base - q;
    if (r < 0) {
        r += 1;
    }

    // Last few digits before the shifted position: ...666666 + q.
    std::vector<int> tail(40, 6);
    int carry = q;
    for (int i = static_cast<int>(tail.size()) - 1; i >= 0 && carry != 0; --i) {
        const int value = tail[static_cast<std::size_t>(i)] + carry;
        if (value >= 0) {
            tail[static_cast<std::size_t>(i)] = value % 10;
            carry = value / 10;
        } else {
            const int borrow = (-value + 9) / 10;
            tail[static_cast<std::size_t>(i)] = value + 10 * borrow;
            carry = -borrow;
        }
    }

    std::string out;
    out.reserve(static_cast<std::size_t>(wanted));

    for (int d : tail) {
        if (d != 6) {
            out.push_back(static_cast<char>('0' + d));
            if (static_cast<int>(out.size()) == wanted) {
                return out;
            }
        }
    }

    std::string rest = extract_non_six_digits_from_fraction(r, wanted - static_cast<int>(out.size()));
    out += rest;
    return out;
}

bool run_checkpoints() {
    const Dec e = boost::multiprecision::exp(Dec(1));
    const Dec h2 = e * e - e;
    const Dec h3 = e * e * e - 2 * e * e + e / 2;

    const std::string d2 =
        extract_non_six_digits_from_fraction(h2 - boost::multiprecision::floor(h2), 8);
    const std::string d3 =
        extract_non_six_digits_from_fraction(h3 - boost::multiprecision::floor(h3), 8);

    if (d2 != "70774270") {
        std::cerr << "Checkpoint failed for H(2): got " << d2 << '\n';
        return false;
    }
    if (d3 != "55395558") {
        std::cerr << "Checkpoint failed for H(3): got " << d3 << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    std::cout << non_six_digits_large_n(1'000'000, 8) << '\n';
    return 0;
}

Python

import decimal
decimal.getcontext().prec = 120

def solve():
    D = decimal.Decimal
    pi = D('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807')
    ln10 = D('2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834')

    alpha = D('-2.0888430156130438559570867167749475005456937410367296732391125442446071101031945')
    beta = D('7.4614892856542545569061166121864153345090949932022092409344113914118766543223747')
    den = alpha*alpha + beta*beta
    n = 1000000; wanted = 8

    two_pi = 2*pi; raw = beta*D(n)
    theta = raw - (raw/two_pi).to_integral_value(rounding=decimal.ROUND_FLOOR)*two_pi
    # cos/sin via Taylor series
    def cos_d(x):
        s = D(1); t = D(1)
        for i in range(1, 80):
            t *= -x*x/D(2*i*(2*i-1)); s += t
        return s
    def sin_d(x):
        s = x; t = x
        for i in range(1, 80):
            t *= -x*x/D((2*i)*(2*i+1)); s += t
        return s
    c = cos_d(theta); s = sin_d(theta)
    envelope = 2*(c*alpha + s*beta)/den
    sign = D(1) if envelope >= 0 else D(-1)

    log10_eps = alpha*D(n)/ln10
    def log10_d(x):
        return x.ln()/ln10
    log10_eps += log10_d(abs(envelope))

    shift = int((-log10_eps).to_integral_value(rounding=decimal.ROUND_FLOOR))
    frac = -log10_eps - D(shift)
    # 10^(-frac)
    def pow10_d(x):
        return (x*ln10).exp()
    t = sign * pow10_d(-frac)

    base = D(2)/3 + t
    q = int(base.to_integral_value(rounding=decimal.ROUND_FLOOR))
    r = base - q
    if r < 0: r += 1

    tail = [6]*40
    carry = q
    for i in range(len(tail)-1, -1, -1):
        if carry == 0: break
        value = tail[i] + carry
        if value >= 0:
            tail[i] = value % 10; carry = value // 10
        else:
            borrow = (-value+9)//10
            tail[i] = value + 10*borrow; carry = -borrow

    out = ''
    for d in tail:
        if d != 6:
            out += str(d)
            if len(out) == wanted: return out

    # Extract from fractional part
    f = r
    while len(out) < wanted:
        f *= 10
        digit = int(f.to_integral_value(rounding=decimal.ROUND_FLOOR))
        digit = max(0, min(9, digit))
        f -= D(digit)
        if digit != 6: out += str(digit)

    return out[:wanted]

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

Java

import java.math.*;

public class Euler970 {
    public static String solve() {
        MathContext MC = new MathContext(120);
        BigDecimal PI = new BigDecimal(
                "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807",
                MC);
        BigDecimal LN10 = new BigDecimal(
                "2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834",
                MC);
        BigDecimal ALPHA = new BigDecimal(
                "-2.0888430156130438559570867167749475005456937410367296732391125442446071101031945", MC);
        BigDecimal BETA = new BigDecimal(
                "7.4614892856542545569061166121864153345090949932022092409344113914118766543223747", MC);
        BigDecimal den = ALPHA.multiply(ALPHA, MC).add(BETA.multiply(BETA, MC), MC);
        int n = 1000000;
        int wanted = 8;
        BigDecimal TWO_PI = PI.multiply(new BigDecimal(2), MC);
        BigDecimal raw = BETA.multiply(new BigDecimal(n), MC);
        BigDecimal theta = raw.subtract(raw.divideToIntegralValue(TWO_PI).multiply(TWO_PI, MC), MC);

        BigDecimal c = cosD(theta, MC), s = sinD(theta, MC);
        BigDecimal envelope = new BigDecimal(2).multiply(c.multiply(ALPHA, MC).add(s.multiply(BETA, MC), MC), MC)
                .divide(den, MC);
        BigDecimal sign = envelope.signum() >= 0 ? BigDecimal.ONE : BigDecimal.ONE.negate();

        BigDecimal log10eps = ALPHA.multiply(new BigDecimal(n), MC).divide(LN10, MC);
        log10eps = log10eps
                .add(envelope.abs().multiply(LN10, MC).stripTrailingZeros().equals(BigDecimal.ZERO) ? BigDecimal.ZERO
                        : logBD(envelope.abs(), MC).divide(LN10, MC), MC);

        int shift = log10eps.negate().setScale(0, RoundingMode.FLOOR).intValue();
        BigDecimal frac = log10eps.negate().subtract(new BigDecimal(shift));
        BigDecimal t = sign.multiply(expBD(frac.negate().multiply(LN10, MC), MC), MC); // 10^(-frac)

        BigDecimal base = new BigDecimal(2).divide(new BigDecimal(3), MC).add(t, MC);
        int q = base.setScale(0, RoundingMode.FLOOR).intValue();
        BigDecimal r = base.subtract(new BigDecimal(q));
        if (r.signum() < 0)
            r = r.add(BigDecimal.ONE);

        int[] tail = new int[40];
        java.util.Arrays.fill(tail, 6);
        int carry = q;
        for (int i = tail.length - 1; i >= 0 && carry != 0; i--) {
            int value = tail[i] + carry;
            if (value >= 0) {
                tail[i] = value % 10;
                carry = value / 10;
            } else {
                int borrow = (-value + 9) / 10;
                tail[i] = value + 10 * borrow;
                carry = -borrow;
            }
        }

        StringBuilder out = new StringBuilder();
        for (int d : tail) {
            if (d != 6) {
                out.append(d);
                if (out.length() == wanted)
                    return out.toString();
            }
        }
        BigDecimal f = r;
        while (out.length() < wanted) {
            f = f.multiply(BigDecimal.TEN);
            int digit = f.setScale(0, RoundingMode.FLOOR).intValue();
            digit = Math.max(0, Math.min(9, digit));
            f = f.subtract(new BigDecimal(digit));
            if (digit != 6)
                out.append(digit);
        }
        return out.substring(0, wanted);
    }

    static BigDecimal cosD(BigDecimal x, MathContext mc) {
        BigDecimal s = BigDecimal.ONE, t = BigDecimal.ONE;
        for (int i = 1; i < 80; i++) {
            t = t.negate().multiply(x, mc).multiply(x, mc).divide(new BigDecimal(2 * i * (2 * i - 1)), mc);
            s = s.add(t, mc);
        }
        return s;
    }

    static BigDecimal sinD(BigDecimal x, MathContext mc) {
        BigDecimal s = x, t = x;
        for (int i = 1; i < 80; i++) {
            t = t.negate().multiply(x, mc).multiply(x, mc).divide(new BigDecimal(2 * i * (2 * i + 1)), mc);
            s = s.add(t, mc);
        }
        return s;
    }

    static BigDecimal logBD(BigDecimal x, MathContext mc) {
        BigDecimal result = new BigDecimal(Math.log(x.doubleValue()), mc);
        for (int i = 0; i < 5; i++) {
            BigDecimal ex = expBD(result, mc);
            BigDecimal diff = x.subtract(ex, mc).divide(ex, mc);
            result = result.add(diff, mc);
        }
        return result;
    }

    static BigDecimal expBD(BigDecimal x, MathContext mc) {
        BigDecimal sum = BigDecimal.ONE, term = BigDecimal.ONE;
        for (int i = 1; i <= 100; i++) {
            term = term.multiply(x, mc).divide(new BigDecimal(i), mc);
            sum = sum.add(term, mc);
            if (term.abs().compareTo(new BigDecimal("1e-110")) < 0)
                break;
        }
        return sum;
    }

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