Detailed mathematical approach

Problem Summary

The local solutions model the pie-eating process through an expectation function \(E(x)\). Here \(x \ge 1\) is the current scaled state, one step is spent immediately, and then the state is multiplied by a random factor \(W \in [0,1]\) whose density is \(2(1-w)\). The process stops as soon as the state falls below 1, so the boundary condition is \(E(u)=0\) for \(u \lt 1\).

Therefore the programming task is not to simulate random trials, but to derive and evaluate the exact formula for \(E(x)\), especially at the default input \(x=40\).

All three implementations use the same closed form:

$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$

Mathematical Approach

Step 1: Write the Renewal Equation

Condition on the first random multiplier \(W\). We always spend one step now, and after that the expected remaining work is \(E(xW)\). Averaging over the density \(2(1-w)\) gives the renewal equation used in the C++ comments:

$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$

Because \(E(xw)=0\) whenever \(xw \lt 1\), only the interval \(w \ge 1/x\) contributes when \(x \ge 1\). So the same equation can be rewritten as

$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$

Step 2: Remove the Inner Scaling

Set \(t=xw\). Then \(dw=dt/x\) and \(w=t/x\), which converts the integral into a Volterra-type expression with fixed lower limit:

$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$

This form is much easier to differentiate because the unknown function now appears as \(E(t)\) instead of \(E(xw)\).

Step 3: Differentiate to Obtain an ODE

Define

$I(x)=\int_1^x (x-t)\,E(t)\,dt.$

The integral equation becomes

$x^2\bigl(E(x)-1\bigr)=2I(x).$

Differentiate once:

$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$

Differentiate a second time:

$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$

The \(E(x)\) terms cancel, leaving the differential equation

$x^2E''(x)+4xE'(x)=2.$

Step 4: Solve the Differential Equation

Let \(Y(x)=E'(x)\). Then

$x^2Y'(x)+4xY(x)=2.$

After division by \(x^2\), this is a first-order linear equation:

$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$

The integrating factor is \(x^4\), so

$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$

Integrating gives

$x^4Y(x)=\frac{2}{3}x^3+C_1,$

hence

$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$

Integrating once more,

$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$

The logarithm is the dominant term, while the homogeneous part contributes the \(x^{-3}\) correction seen in the code.

Step 5: Fix the Constants from the Boundary Data

At \(x=1\), the integral part vanishes, so

$E(1)=1.$

Now evaluate the once-differentiated identity at \(x=1\):

$2\bigl(E(1)-1\bigr)+E'(1)=0.$

Since \(E(1)=1\), this yields

$E'(1)=0.$

Substitute into the general solution:

$0=E'(1)=\frac{2}{3}+C_1,$

so \(C_1=-\frac{2}{3}\). Then

$1=E(1)=0+\frac{2}{9}+C_2,$

which gives \(C_2=\frac{7}{9}\). Therefore

$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$

Step 6: Checkpoints and Final Evaluation

The C++ program verifies the formula at three checkpoints, and they match exactly:

$E(1)=1.0000000000,$

$E(2)=1.2676536759,$

$E(7.5)=2.1215732071.$

For the default Project Euler input, the same formula gives

$E(40)=3.2370342194.$

As \(x\) grows, the correction term decays rapidly, so the expectation behaves like

$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$

How the Code Works

The C++ solution stores the mathematics in expected_steps(x) and evaluates the closed form with long double. It also supports an optional --x=... argument and a --skip-checkpoints switch, while the helper nearly_equal validates the known values \(E(1)\), \(E(2)\), and \(E(7.5)\).

The Python and Java versions intentionally strip the implementation down to the essential formula. Each computes the same expression and formats the value at \(x=40\) with 10 digits after the decimal point. None of the three languages performs simulation or numerical integration at runtime.

Complexity Analysis

Once the closed form is known, evaluation needs one logarithm and a constant number of arithmetic operations. Therefore the running time is \(O(1)\) and the memory usage is \(O(1)\). This is far better than Monte Carlo sampling or iterative quadrature, both of which would need many repeated evaluations to approach the same precision.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=394
2. Renewal theory overview: Wikipedia - Renewal theory
3. Integral equation background: Wikipedia - Integral equation
4. Volterra integral equations: Wikipedia - Volterra integral equation
5. First-order linear differential equations: Wikipedia - Linear differential equation