Problem 444: The Roundtable Lottery

View on Project Euler

Project Euler Problem 444 Solution

EulerSolve provides an optimized solution for Project Euler Problem 444, The Roundtable Lottery, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The round-table lottery is reduced, in the implementations, to an expectation kernel and its repeated cumulative sums. The basic quantity is $$E(n)=H_n=\sum_{j=1}^{n}\frac{1}{j},$$ where \(H_n\) is the \(n\)-th harmonic number. From that kernel we define $$S_0(N)=E(N),\qquad S_r(N)=\sum_{m=1}^{N} S_{r-1}(m)\quad(r\ge 1).$$ The required output is \(S_{20}(10^{14})\) in scientific notation with 10 significant digits. The implementations also use the checkpoints \(E(111)=5.2912\) and \(S_3(100)=5.983679014\times 10^5\). Mathematical Approach Step 1: The expectation kernel is harmonic After the probabilistic part of the problem is simplified, the only primitive quantity that remains is the harmonic number. That is why the first validation step is simply the direct evaluation of $$E(n)=H_n.$$ So the task is no longer a simulation over many random eliminations. It becomes a summation problem built on reciprocal terms. Step 2: Repeated summation gives a weighted reciprocal sum Unrolling the recurrence for \(S_r(N)\) shows how often each reciprocal term \(1/j\) appears after \(k\) cumulative sums....

Detailed mathematical approach

Problem Summary

The round-table lottery is reduced, in the implementations, to an expectation kernel and its repeated cumulative sums. The basic quantity is

$$E(n)=H_n=\sum_{j=1}^{n}\frac{1}{j},$$

where \(H_n\) is the \(n\)-th harmonic number. From that kernel we define

$$S_0(N)=E(N),\qquad S_r(N)=\sum_{m=1}^{N} S_{r-1}(m)\quad(r\ge 1).$$

The required output is \(S_{20}(10^{14})\) in scientific notation with 10 significant digits. The implementations also use the checkpoints \(E(111)=5.2912\) and \(S_3(100)=5.983679014\times 10^5\).

Mathematical Approach

Step 1: The expectation kernel is harmonic

After the probabilistic part of the problem is simplified, the only primitive quantity that remains is the harmonic number. That is why the first validation step is simply the direct evaluation of

$$E(n)=H_n.$$

So the task is no longer a simulation over many random eliminations. It becomes a summation problem built on reciprocal terms.

Step 2: Repeated summation gives a weighted reciprocal sum

Unrolling the recurrence for \(S_r(N)\) shows how often each reciprocal term \(1/j\) appears after \(k\) cumulative sums. For fixed \(j\), it contributes once for every weakly increasing chain

$$j\le n_1\le n_2\le \cdots \le n_k\le N.$$

The number of such chains is the stars-and-bars coefficient

$$\binom{N-j+k}{k}.$$

Therefore the \(k\)-fold cumulative sum can be written as one explicit sum:

$$S_k(N)=\sum_{j=1}^{N}\frac{1}{j}\binom{N-j+k}{k}.$$

This identity already explains why the final algorithm never needs to enumerate anything near \(10^{14}\).

Step 3: Closed form for the cumulative transform

Introduce

$$T_k(N)=\binom{N+k}{k}\bigl(H_{N+k}-H_k\bigr).$$

For \(k=0\), this gives \(T_0(N)=H_N=S_0(N)\). Now compare first differences. Using

$$\binom{N+k}{k}=\binom{N+k-1}{k}+\binom{N+k-1}{k-1},$$

$$H_{N+k}=H_{N+k-1}+\frac{1}{N+k},$$

and

$$\frac{1}{N+k}\binom{N+k}{k}=\frac{1}{k}\binom{N+k-1}{k-1},$$

we obtain

$$T_k(N)-T_k(N-1)=\binom{N+k-1}{k-1}\bigl(H_{N+k-1}-H_{k-1}\bigr)=T_{k-1}(N).$$

The repeated sums satisfy the same recurrence:

$$S_k(N)-S_k(N-1)=S_{k-1}(N),\qquad S_k(0)=0.$$

Since \(S_0(N)=T_0(N)\), induction on \(k\) gives the exact identity

$$\boxed{S_k(N)=\binom{N+k}{k}\bigl(H_{N+k}-H_k\bigr).}$$

Step 4: Numerical checks and final target

The closed form reproduces the validation values used by the implementations:

$$E(111)=H_{111}\approx 5.2912,$$

$$S_3(100)=\binom{103}{3}\bigl(H_{103}-H_3\bigr)\approx 5.983679014\times 10^5.$$

For the required parameters, the same formula yields

$$S_{20}(10^{14})\approx 1.200856722\times 10^{263}.$$

Step 5: Harmonic numbers for huge arguments

Only the argument \(N+k\) is enormous. Small harmonic numbers such as \(H_k\) are summed directly, while the large argument is evaluated with the Euler-Maclaurin expansion

$$\begin{aligned} H_n = {}& \log n + \gamma + \frac{1}{2n} - \frac{1}{12n^2} + \frac{1}{120n^4} - \frac{1}{252n^6} + \frac{1}{240n^8} \\ &- \frac{5}{660n^{10}} + \frac{691}{32760n^{12}} - \frac{1}{12n^{14}} + \frac{3617}{8160n^{16}} \\ &- \frac{43867}{14364n^{18}} + \frac{174611}{6600n^{20}} + O(n^{-22}). \end{aligned}$$

At \(n=10^{14}+20\), the neglected tail is far below the requested 10-digit output accuracy.

How the Code Works

The C++, Python, and Java implementations all use the same plan. They compute the binomial factor through the short product

$$\binom{N+k}{k}=\prod_{i=1}^{k}\frac{N+i}{i},$$

which avoids factorial overflow and costs only \(O(k)\) operations because \(k=20\) is fixed. They then evaluate \(H_{N+k}\) with high-precision arithmetic, subtract the directly computed \(H_k\), multiply the two factors, and format the answer in normalized scientific notation.

Complexity Analysis

For the target instance, the binomial product takes \(O(k)\) arithmetic steps, the exact evaluation of \(H_k\) also takes \(O(k)\), and the large harmonic number uses a constant number of asymptotic correction terms. Since \(k=20\) is fixed, the practical cost is \(O(1)\) time and \(O(1)\) memory. The decisive improvement is that the method never iterates up to \(N=10^{14}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=444
  2. Harmonic number: Wikipedia - Harmonic number
  3. Euler-Maclaurin formula: Wikipedia - Euler-Maclaurin formula
  4. Hockey-stick identity: Wikipedia - Hockey-stick identity
  5. Graham, Knuth, Patashnik, Concrete Mathematics, chapter on sums and binomial identities.

Problem 444 source code

C++

#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>

namespace {

using u64 = std::uint64_t;
using Big = boost::multiprecision::cpp_dec_float_100;

struct Options {
    u64 n = 100'000'000'000'000ULL;
    u64 k = 20ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(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_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        if (parse_u64_after_prefix(arg, "--k=", options.k)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1ULL;
}

Big harmonic_small(const u64 n) {
    Big h = 0;
    for (u64 i = 1; i <= n; ++i) {
        h += Big(1) / Big(i);
    }
    return h;
}

Big harmonic_large(const u64 n) {
    if (n <= 1'000'000ULL) {
        return harmonic_small(n);
    }

    const Big x = Big(n);
    const Big inv = Big(1) / x;
    const Big inv2 = inv * inv;
    const Big inv4 = inv2 * inv2;
    const Big inv6 = inv4 * inv2;
    const Big inv8 = inv4 * inv4;
    const Big inv10 = inv8 * inv2;
    const Big inv12 = inv10 * inv2;
    const Big inv14 = inv12 * inv2;
    const Big inv16 = inv14 * inv2;
    const Big inv18 = inv16 * inv2;
    const Big inv20 = inv18 * inv2;

    return log(x) + boost::math::constants::euler<Big>() + inv / 2 - inv2 / 12 + inv4 / 120 -
           inv6 / 252 + inv8 / 240 - Big(5) * inv10 / 660 + Big(691) * inv12 / 32760 -
           inv14 / 12 + Big(3617) * inv16 / 8160 - Big(43867) * inv18 / 14364 +
           Big(174611) * inv20 / 6600;
}

Big binomial_n_plus_k_choose_k(const u64 n, const u64 k) {
    Big c = 1;
    for (u64 i = 1; i <= k; ++i) {
        c *= Big(n + i);
        c /= Big(i);
    }
    return c;
}

Big compute_s_k(const u64 n, const u64 k) {
    const Big comb = binomial_n_plus_k_choose_k(n, k);
    const Big h_big = harmonic_large(n + k);
    const Big h_k = harmonic_small(k);
    return comb * (h_big - h_k);
}

std::string normalize_scientific(const std::string& s) {
    const std::size_t pos = s.find('e');
    if (pos == std::string::npos) {
        return s;
    }

    const std::string mantissa = s.substr(0, pos);
    std::string exponent = s.substr(pos + 1);

    bool negative = false;
    if (!exponent.empty() && (exponent[0] == '+' || exponent[0] == '-')) {
        negative = (exponent[0] == '-');
        exponent.erase(exponent.begin());
    }

    while (exponent.size() > 1U && exponent[0] == '0') {
        exponent.erase(exponent.begin());
    }

    return mantissa + (negative ? "e-" : "e") + exponent;
}

std::string format_scientific_10sig(const Big& value) {
    std::ostringstream oss;
    oss << std::scientific << std::setprecision(9) << value;
    return normalize_scientific(oss.str());
}

bool run_checkpoints() {
    std::ostringstream e111;
    e111 << std::fixed << std::setprecision(4) << harmonic_small(111ULL);
    if (e111.str() != "5.2912") {
        std::cerr << "Checkpoint failed: E(111)\n";
        return false;
    }

    const std::string s3 = format_scientific_10sig(compute_s_k(100ULL, 3ULL));
    if (s3 != "5.983679014e5") {
        std::cerr << "Checkpoint failed: S3(100)\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 << format_scientific_10sig(compute_s_k(options.n, options.k)) << '\n';
    return 0;
}

Python

from decimal import Decimal, getcontext
import math

def solve():
    getcontext().prec = 120
    n = 100000000000000
    k = 20

    euler_gamma = Decimal('0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495')

    def harmonic_large(n_):
        if n_ <= 1000000:
            h = Decimal(0)
            for i in range(1, n_ + 1): h += Decimal(1) / Decimal(i)
            return h
        x = Decimal(n_)
        inv = Decimal(1) / x
        inv2 = inv * inv
        inv4 = inv2 * inv2
        inv6 = inv4 * inv2
        inv8 = inv4 * inv4
        inv10 = inv8 * inv2
        inv12 = inv10 * inv2
        inv14 = inv12 * inv2
        inv16 = inv14 * inv2
        inv18 = inv16 * inv2
        inv20 = inv18 * inv2
        return (x.ln() + euler_gamma + inv/2 - inv2/12 + inv4/120 - inv6/252
                + inv8/240 - 5*inv10/660 + 691*inv12/32760 - inv14/12
                + 3617*inv16/8160 - 43867*inv18/14364 + 174611*inv20/6600)

    def harmonic_small(n_):
        h = Decimal(0)
        for i in range(1, n_ + 1): h += Decimal(1) / Decimal(i)
        return h

    comb = Decimal(1)
    for i in range(1, k + 1):
        comb *= Decimal(n + i)
        comb /= Decimal(i)
    result = comb * (harmonic_large(n + k) - harmonic_small(k))

    # Format as scientific with 10 sig digits
    s = f'{result:.9e}'
    # Normalize exponent
    parts = s.split('e')
    mantissa = parts[0]
    exp_str = parts[1]
    exp_val = int(exp_str)
    return f'{mantissa}e{exp_val}'

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

Java

import java.math.*;

public class Euler444 {
    public static String solve() {
        MathContext mc = new MathContext(120);
        long n = 100000000000000L;
        int k = 20;
        BigDecimal EG = new BigDecimal(
                "0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495",
                mc);
        BigDecimal x = new BigDecimal(n + k, mc);
        BigDecimal inv = BigDecimal.ONE.divide(x, mc);
        BigDecimal inv2 = inv.multiply(inv, mc);
        BigDecimal inv4 = inv2.multiply(inv2, mc), inv6 = inv4.multiply(inv2, mc), inv8 = inv4.multiply(inv4, mc);
        BigDecimal inv10 = inv8.multiply(inv2, mc), inv12 = inv10.multiply(inv2, mc), inv14 = inv12.multiply(inv2, mc);
        BigDecimal inv16 = inv14.multiply(inv2, mc), inv18 = inv16.multiply(inv2, mc), inv20 = inv18.multiply(inv2, mc);
        BigDecimal lnx = ln(x, mc);
        BigDecimal H = lnx.add(EG, mc).add(inv.divide(new BigDecimal(2), mc), mc)
                .subtract(inv2.divide(new BigDecimal(12), mc), mc).add(inv4.divide(new BigDecimal(120), mc), mc)
                .subtract(inv6.divide(new BigDecimal(252), mc), mc).add(inv8.divide(new BigDecimal(240), mc), mc)
                .subtract(new BigDecimal(5).multiply(inv10, mc).divide(new BigDecimal(660), mc), mc)
                .add(new BigDecimal(691).multiply(inv12, mc).divide(new BigDecimal(32760), mc), mc)
                .subtract(inv14.divide(new BigDecimal(12), mc), mc)
                .add(new BigDecimal(3617).multiply(inv16, mc).divide(new BigDecimal(8160), mc), mc)
                .subtract(new BigDecimal(43867).multiply(inv18, mc).divide(new BigDecimal(14364), mc), mc)
                .add(new BigDecimal(174611).multiply(inv20, mc).divide(new BigDecimal(6600), mc), mc);
        BigDecimal Hk = BigDecimal.ZERO;
        for (int i = 1; i <= k; i++)
            Hk = Hk.add(BigDecimal.ONE.divide(new BigDecimal(i), mc), mc);
        BigDecimal comb = BigDecimal.ONE;
        for (int i = 1; i <= k; i++) {
            comb = comb.multiply(new BigDecimal(n + i), mc);
            comb = comb.divide(new BigDecimal(i), mc);
        }
        BigDecimal result = comb.multiply(H.subtract(Hk, mc), mc);
        String s = String.format("%.9e", result);
        String[] parts = s.split("e");
        int exp = Integer.parseInt(parts[1].startsWith("+") ? parts[1].substring(1) : parts[1]);
        return parts[0] + "e" + exp;
    }

    static BigDecimal ln(BigDecimal x, MathContext mc) {
        // Use log10 conversion: ln(x) = log10(x) / log10(e)
        double approx = Math.log(x.doubleValue());
        BigDecimal result = new BigDecimal(approx, mc);
        // Newton refinement: ln(x) via exp
        for (int i = 0; i < 5; i++) {
            BigDecimal ex = exp(result, mc);
            BigDecimal diff = x.subtract(ex, mc).divide(ex, mc);
            result = result.add(diff, mc);
        }
        return result;
    }

    static BigDecimal exp(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());
    }
}