Problem 394: Eating Pie

Project Euler Problem 394 Solution

EulerSolve provides an optimized solution for Project Euler Problem 394, Eating Pie, with C++, Python, Java, and a step-by-step mathematical explanation.

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

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

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

Problem 394 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>

namespace {

struct Options {
    long double x = 40.0L;
    bool run_checkpoints = true;
};

bool parse_ld_after_prefix(const std::string& arg, const std::string& prefix, long double& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stold(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_ld_after_prefix(arg, "--x=", options.x)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.x >= 1.0L;
}

long double expected_steps(const long double x) {
    // Closed form from the renewal equation:
    // E(x) = 1 + 2 * integral_0^1 (1-w) E(xw) dw (with E(u)=0 for u<1),
    // yielding E(x) = 7/9 + (2/3) ln x + 2/(9 x^3).
    return 7.0L / 9.0L + (2.0L / 3.0L) * std::log(x) + 2.0L / (9.0L * x * x * x);
}

bool nearly_equal(const long double a, const long double b, const long double rel_tol = 1e-12L) {
    const long double scale = std::max(std::fabsl(a), std::fabsl(b));
    if (scale == 0.0L) {
        return true;
    }
    return std::fabsl(a - b) <= rel_tol * scale;
}

bool run_checkpoints() {
    if (!nearly_equal(expected_steps(1.0L), 1.0L, 1e-14L)) {
        std::cerr << "Checkpoint failed: E(1)\n";
        return false;
    }

    if (!nearly_equal(expected_steps(2.0L), 1.2676536759L, 1e-10L)) {
        std::cerr << "Checkpoint failed: E(2)\n";
        return false;
    }

    if (!nearly_equal(expected_steps(7.5L), 2.1215732071L, 1e-10L)) {
        std::cerr << "Checkpoint failed: E(7.5)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(expected_steps(options.x)) << '\n';
    return 0;
}

Python

import math

def expected_steps(x):
    return 7.0 / 9.0 + (2.0 / 3.0) * math.log(x) + 2.0 / (9.0 * x * x * x)

def solve():
    return "{:.10f}".format(expected_steps(40.0))

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

Java

public class Euler394 {
    private static double expectedSteps(double x) {
        return 7.0 / 9.0 + (2.0 / 3.0) * Math.log(x) + 2.0 / (9.0 * x * x * x);
    }

    public static String solve() {
        return String.format("%.10f", expectedSteps(40.0)).replace(',', '.');
    }

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