Problem 436: Unfair Wager

View on Project Euler

Project Euler Problem 436 Solution

EulerSolve provides an optimized solution for Project Euler Problem 436, Unfair Wager, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Louise and Julie draw independent \(U(0,1)\) numbers and keep a common running total. Louise stops when the total first exceeds \(1\) and records her final draw \(x\). Julie then continues from that same total, stops when the total first exceeds \(2\), and records her final draw \(y\). The goal is to compute $$\Pr(y \gt x).$$ The implementations evaluate this probability exactly and obtain $$\boxed{\Pr(y \gt x)=\frac{1+14e-5e^2}{4}\approx 0.5276662759.}$$ So Julie is favored, but only slightly. Mathematical Approach The solution has three clean pieces: first derive the distribution of the final draw for crossing a threshold in \([0,1]\), then derive the joint density of Louise's last draw and overshoot above \(1\), and finally integrate the conditional probability that Julie's last draw is larger. Step 1: Distribution of the Last Draw for a Residual Threshold For \(0 \le t \le 1\) and \(0 \le z \le 1\), define $$G(t,z)=\Pr(L_t \le z).$$ Condition on the first new draw \(u \sim U(0,1)\). If \(u \gt t\), then the threshold is crossed immediately and the last draw is exactly \(u\). This contributes $$\Pr(t \lt u \le z)=\max(0,z-t).$$ If \(u \le t\), then after subtracting \(u\) from the remaining distance, the problem restarts with threshold \(t-u\)....

Detailed mathematical approach

Problem Summary

Louise and Julie draw independent \(U(0,1)\) numbers and keep a common running total. Louise stops when the total first exceeds \(1\) and records her final draw \(x\). Julie then continues from that same total, stops when the total first exceeds \(2\), and records her final draw \(y\). The goal is to compute

$$\Pr(y \gt x).$$

The implementations evaluate this probability exactly and obtain

$$\boxed{\Pr(y \gt x)=\frac{1+14e-5e^2}{4}\approx 0.5276662759.}$$

So Julie is favored, but only slightly.

Mathematical Approach

The solution has three clean pieces: first derive the distribution of the final draw for crossing a threshold in \([0,1]\), then derive the joint density of Louise's last draw and overshoot above \(1\), and finally integrate the conditional probability that Julie's last draw is larger.

Step 1: Distribution of the Last Draw for a Residual Threshold

For \(0 \le t \le 1\) and \(0 \le z \le 1\), define

$$G(t,z)=\Pr(L_t \le z).$$

Condition on the first new draw \(u \sim U(0,1)\).

If \(u \gt t\), then the threshold is crossed immediately and the last draw is exactly \(u\). This contributes

$$\Pr(t \lt u \le z)=\max(0,z-t).$$

If \(u \le t\), then after subtracting \(u\) from the remaining distance, the problem restarts with threshold \(t-u\). Averaging over \(u\) gives the renewal equation

$$G(t,z)=\int_0^t G(t-u,z)\,du+\max(0,z-t).$$

With the substitution \(v=t-u\), this becomes

$$G(t,z)=\int_0^t G(v,z)\,dv+\max(0,z-t).$$

Differentiate with respect to \(t\).

For \(0 \le t \lt z\), the boundary term is \(z-t\), so

$$\frac{\partial G}{\partial t}=G-1,\qquad G(0,z)=z.$$

The solution is

$$G(t,z)=1-(1-z)e^t \qquad (0 \le t \lt z \le 1).$$

For \(0 \le z \le t \le 1\), the boundary term is \(0\), so

$$\frac{\partial G}{\partial t}=G.$$

Using continuity at \(t=z\), we obtain

$$G(t,z)=e^t\bigl(e^{-z}-1+z\bigr)\qquad (0 \le z \le t \le 1).$$

Thus the last-draw CDF is

$$G(t,z)= \begin{cases} 1-(1-z)e^t, & 0 \le t \lt z \le 1,\\[4pt] e^t\bigl(e^{-z}-1+z\bigr), & 0 \le z \le t \le 1. \end{cases}$$

It also satisfies the expected boundary values \(G(0,z)=z\), \(G(t,0)=0\), and \(G(t,1)=1\).

Step 2: Joint Density of Louise's Last Draw and Overshoot

Let \(s\) be the running total immediately before Louise's last draw. Then \(0 \lt s \lt 1\). If this state is reached after exactly \(n\) previous draws, the density of the sum of \(n\) independent \(U(0,1)\) variables on the interval \([0,1]\) is the first branch of the Irwin-Hall distribution:

$$f_n(s)=\frac{s^{n-1}}{(n-1)!}\qquad (0 \le s \le 1,\ n \ge 1).$$

Summing over all possible numbers of previous draws gives the density of the pre-crossing state:

$$r(s)=\sum_{n=1}^{\infty}\frac{s^{n-1}}{(n-1)!}=e^s \qquad (0 \le s \le 1).$$

Now let \(x\) be Louise's last draw and let \(w\) be the overshoot above \(1\). Then

$$s+x=1+w,\qquad s=1+w-x.$$

Because \(0 \le s \le 1\) and \(0 \le x \le 1\), the admissible region is

$$0 \le w \le x \le 1.$$

Given \(s\), the last draw is uniform on the interval \((1-s,1)\), so its density is simply \(1\) there. Therefore the joint density of \((x,w)\) is

$$f(x,w)=e^{1+w-x}\qquad (0 \le w \le x \le 1).$$

A quick normalization check confirms that this is correct:

$$\int_0^1\int_w^1 e^{1+w-x}\,dx\,dw=1.$$

Step 3: Condition Julie's Phase on Louise's Terminal State

After Louise stops, the common running total is \(1+w\). Julie therefore needs only an additional amount

$$t=1-w$$

to exceed \(2\). Conditioned on \((x,w)\), Julie's phase is a fresh threshold-crossing problem with threshold \(1-w\). Hence

$$\Pr(y \gt x \mid x,w)=1-G(1-w,x).$$

Integrating over Louise's terminal state gives

$$\Pr(y \gt x)=\int_0^1\int_w^1 e^{1+w-x}\bigl(1-G(1-w,x)\bigr)\,dx\,dw.$$

Step 4: Split the Domain Along \(x+w=1\)

The formula for \(G\) depends on whether \(1-w\) is smaller or larger than \(x\).

If \(x+w \gt 1\), then \(1-w \lt x\), so the first branch applies:

$$1-G(1-w,x)=(1-x)e^{1-w}.$$

Multiplying by the joint density gives

$$e^{1+w-x}\bigl(1-G(1-w,x)\bigr)=(1-x)e^{2-x}.$$

If \(x+w \le 1\), then \(1-w \ge x\), so the second branch applies:

$$1-G(1-w,x)=1-e^{1-w}\bigl(e^{-x}-1+x\bigr).$$

After multiplication by \(e^{1+w-x}\), the integrand becomes

$$e^{1+w-x}-e^{2-2x}+(1-x)e^{2-x}.$$

Therefore we split the integral along the line \(x=1-w\).

Step 5: Evaluate the Integral Exactly

For \(0 \le w \le \tfrac12\), the line \(x=1-w\) lies inside the interval \([w,1]\). For \(w \ge \tfrac12\), the whole interval satisfies \(x+w \gt 1\). Thus

$$\begin{aligned} \Pr(y \gt x)=&\int_0^{1/2}\left[\int_w^{1-w}\left(e^{1+w-x}-e^{2-2x}+(1-x)e^{2-x}\right)\,dx\right]dw\\ &+\int_0^{1/2}\left[\int_{1-w}^{1}(1-x)e^{2-x}\,dx\right]dw\\ &+\int_{1/2}^{1}\left[\int_w^{1}(1-x)e^{2-x}\,dx\right]dw. \end{aligned}$$

The inner integrals simplify to

$$I_1(w)=2e-we^{2-w}-\frac{e^{2w}}{2}-\frac{e^{2-2w}}{2}\qquad \left(0 \le w \le \frac12\right),$$

$$I_2(w)=e-we^{2-w}\qquad \left(\frac12 \le w \le 1\right).$$

So we only need two one-dimensional integrals:

$$A=\int_0^{1/2} I_1(w)\,dw=\frac14+e-\frac{5e^2}{4}+\frac{3e^{3/2}}{2},$$

$$B=\int_{1/2}^{1} I_2(w)\,dw=\frac{5e}{2}-\frac{3e^{3/2}}{2}.$$

The \(e^{3/2}\) terms cancel, leaving

$$\Pr(y \gt x)=A+B=\boxed{\frac{1+14e-5e^2}{4}}.$$

Numerically,

$$\Pr(y \gt x)\approx 0.5276662759.$$

How the Code Works

The C++, Python, and Java implementations all use the exact closed form \(\frac{1+14e-5e^2}{4}\) as the final result. The C++ implementation additionally performs independent checks: it verifies the boundary behavior of the last-draw distribution, tests the renewal identity numerically at several sample points, and compares the exact closed form against a separate two-dimensional Simpson integration of the probability integral. The Python and Java implementations simply evaluate the closed form and format the answer to 10 decimal places.

Complexity Analysis

The exact computation is \(O(1)\) time and \(O(1)\) memory because it reduces to a fixed closed-form expression. The numerical verification in the C++ implementation uses Simpson quadrature on a grid of \(N_w\) overshoot slices and \(N_x\) inner points, so that validation path costs \(O(N_wN_x)\) time. Its outer slices are independent, so that validation can be parallelized. Memory usage remains \(O(N_w)\) for the stored inner values.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=436
  2. Renewal theory: Wikipedia — Renewal theory
  3. Irwin-Hall distribution: Wikipedia — Irwin-Hall distribution
  4. Simpson's rule: Wikipedia — Simpson's rule

Problem 436 source code

C++

#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

struct Options {
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    std::uint64_t parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<std::uint64_t>(c - '0');
        if (parsed > static_cast<std::uint64_t>(std::numeric_limits<unsigned>::max())) {
            return false;
        }
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

unsigned pick_thread_count(const unsigned requested) {
    if (requested > 0U) {
        return requested;
    }
    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) {
        hw = 4U;
    }
    return hw;
}

long double cdf_last_draw_crossing_t(const long double t, const long double z) {
    if (z <= 0.0L) {
        return 0.0L;
    }
    if (z >= 1.0L) {
        return 1.0L;
    }
    if (t <= 0.0L) {
        return z;
    }
    if (t < z) {
        return 1.0L - (1.0L - z) * std::exp(t);
    }
    return std::exp(t) * (std::exp(-z) - 1.0L + z);
}

long double exact_probability() {
    // Closed form from integrating the exact joint density of (x, overshoot at 1).
    const long double e = std::exp(1.0L);
    return (1.0L + 14.0L * e - 5.0L * e * e) / 4.0L;
}

long double simpson_integral_cdf_segment(const long double a,
                                         const long double b,
                                         const long double z,
                                         const int n_even) {
    assert(n_even > 0 && (n_even % 2 == 0));
    if (b <= a) {
        return 0.0L;
    }
    const long double h = (b - a) / static_cast<long double>(n_even);
    long double sum = cdf_last_draw_crossing_t(a, z) + cdf_last_draw_crossing_t(b, z);
    for (int i = 1; i < n_even; ++i) {
        const long double v = a + h * static_cast<long double>(i);
        const long double w = (i & 1) ? 4.0L : 2.0L;
        sum += w * cdf_last_draw_crossing_t(v, z);
    }
    return sum * h / 3.0L;
}

long double simpson_integral_cdf(const long double t, const long double z, const int n_even) {
    assert(n_even > 0 && (n_even % 2 == 0));
    if (t <= 0.0L) {
        return 0.0L;
    }

    const long double split = std::min(t, std::max(0.0L, z));
    if (split <= 0.0L || split >= t) {
        return simpson_integral_cdf_segment(0.0L, t, z, n_even);
    }

    return simpson_integral_cdf_segment(0.0L, split, z, n_even) +
           simpson_integral_cdf_segment(split, t, z, n_even);
}

long double probability_integrand(const long double w, const long double x) {
    const long double t = 1.0L - w;
    const long double q_xw = std::exp(1.0L - x + w);  // Joint density of (x, overshoot at 1).
    const long double y_beats_x = 1.0L - cdf_last_draw_crossing_t(t, x);
    return q_xw * y_beats_x;
}

long double inner_integral_for_w(const long double w, const int n_even_x) {
    assert(n_even_x > 0 && (n_even_x % 2 == 0));
    if (w >= 1.0L) {
        return 0.0L;
    }

    const long double a = w;
    const long double b = 1.0L;
    const long double h = (b - a) / static_cast<long double>(n_even_x);
    long double sum = probability_integrand(w, a) + probability_integrand(w, b);
    for (int i = 1; i < n_even_x; ++i) {
        const long double x = a + h * static_cast<long double>(i);
        const long double wt = (i & 1) ? 4.0L : 2.0L;
        sum += wt * probability_integrand(w, x);
    }
    return sum * h / 3.0L;
}

long double numerical_probability_check(const int n_even_w,
                                        const int n_even_x,
                                        const unsigned thread_count) {
    assert(n_even_w > 0 && (n_even_w % 2 == 0));
    assert(n_even_x > 0 && (n_even_x % 2 == 0));

    std::vector<long double> inner(static_cast<std::size_t>(n_even_w) + 1ULL, 0.0L);
    const long double hw = 1.0L / static_cast<long double>(n_even_w);
    std::atomic<int> next_i(0);

    auto worker = [&]() {
        while (true) {
            const int i = next_i.fetch_add(1, std::memory_order_relaxed);
            if (i > n_even_w) {
                return;
            }
            const long double w = hw * static_cast<long double>(i);
            inner[static_cast<std::size_t>(i)] = inner_integral_for_w(w, n_even_x);
        }
    };

    std::vector<std::thread> pool;
    pool.reserve(thread_count);
    for (unsigned i = 0U; i < thread_count; ++i) {
        pool.emplace_back(worker);
    }
    for (auto& th : pool) {
        th.join();
    }

    long double sum = inner.front() + inner.back();
    for (int i = 1; i < n_even_w; ++i) {
        const long double wt = (i & 1) ? 4.0L : 2.0L;
        sum += wt * inner[static_cast<std::size_t>(i)];
    }
    return sum * hw / 3.0L;
}

void run_checkpoints(const unsigned thread_count) {
    constexpr long double kEps = 1e-12L;

    // Basic boundary checks for the CDF.
    assert(std::fabs(cdf_last_draw_crossing_t(0.0L, 0.37L) - 0.37L) < kEps);
    assert(std::fabs(cdf_last_draw_crossing_t(0.42L, 0.0L) - 0.0L) < kEps);
    assert(std::fabs(cdf_last_draw_crossing_t(0.42L, 1.0L) - 1.0L) < kEps);

    // Renewal equation check: G(t,z) = ∫_0^t G(v,z)dv + max(0, z - t).
    const std::vector<long double> ts = {0.11L, 0.29L, 0.53L, 0.77L, 0.93L};
    const std::vector<long double> zs = {0.07L, 0.19L, 0.41L, 0.68L, 0.95L};
    for (const long double t : ts) {
        for (const long double z : zs) {
            const long double lhs = cdf_last_draw_crossing_t(t, z);
            const long double integral = simpson_integral_cdf(t, z, 2000);
            const long double rhs = integral + ((z > t) ? (z - t) : 0.0L);
            assert(std::fabs(lhs - rhs) < 5e-11L);
        }
    }

    // Independent numerical integration of the full 2D probability.
    const long double p_exact = exact_probability();
    const long double p_numeric = numerical_probability_check(600, 1200, thread_count);
    assert(std::fabs(p_numeric - p_exact) < 2e-9L);
}

}  // namespace

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

    const unsigned thread_count = pick_thread_count(options.requested_threads);
    if (options.run_checkpoints) {
        run_checkpoints(thread_count);
    }

    const long double ans = exact_probability();
    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

import math

def solve():
    e = math.e
    ans = (1.0 + 14.0 * e - 5.0 * e * e) / 4.0
    return f"{ans:.10f}"

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

Java

public class Euler436 {
    public static String solve() {
        double e = Math.E;
        double ans = (1.0 + 14.0 * e - 5.0 * e * e) / 4.0;
        return String.format(java.util.Locale.US, "%.10f", ans);
    }

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