Problem 423: Consecutive Die Throws

Project Euler Problem 423 Solution

EulerSolve provides an optimized solution for Project Euler Problem 423, Consecutive Die Throws, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each $n$, the counting formula used by the implementations is $$C(n)=6\sum_{k=0}^{\pi(n)} \binom{n-1}{k}5^{\,n-1-k},$$ where $\pi(n)$ is the number of primes not exceeding $n$. The final target is the cumulative sum $$S(L)=\sum_{n=1}^{L} C(n)\pmod{10^9+7},\qquad L=50{,}000{,}000.$$ A direct evaluation of every truncated binomial sum would be far too slow, so the solution derives a constant-time transition from $n$ to $n+1$. Mathematical Approach Step 1: Isolate the Summand Define $$B_n(k)=\binom{n-1}{k}5^{\,n-1-k},\qquad t_n=\pi(n).$$ Then the inner truncated sum is $$P_n=\sum_{k=0}^{t_n} B_n(k),$$ so the quantity contributed at level $n$ is simply $$C(n)=6P_n.$$ The whole problem is therefore reduced to updating $P_n$ efficiently as $n$ increases. Step 2: Pascal's Identity Gives the Row Transition Using $$\binom{n}{k}=\binom{n-1}{k}+\binom{n-1}{k-1},$$ we obtain $$B_{n+1}(k)=\binom{n}{k}5^{\,n-k}=5B_n(k)+B_n(k-1).$$ This recurrence is the key structural fact. It says that every term in the next row depends only on two neighboring terms from the current row, so there is no need to rebuild the full truncated sum from scratch. Step 3: Only the Boundary Matters The cutoff $t_n=\pi(n)$ changes by at most one when $n$ increases....

Detailed mathematical approach

Problem Summary

For each $n$, the counting formula used by the implementations is

$$C(n)=6\sum_{k=0}^{\pi(n)} \binom{n-1}{k}5^{\,n-1-k},$$

where $\pi(n)$ is the number of primes not exceeding $n$. The final target is the cumulative sum

$$S(L)=\sum_{n=1}^{L} C(n)\pmod{10^9+7},\qquad L=50{,}000{,}000.$$

A direct evaluation of every truncated binomial sum would be far too slow, so the solution derives a constant-time transition from $n$ to $n+1$.

Mathematical Approach

Step 1: Isolate the Summand

Define

$$B_n(k)=\binom{n-1}{k}5^{\,n-1-k},\qquad t_n=\pi(n).$$

Then the inner truncated sum is

$$P_n=\sum_{k=0}^{t_n} B_n(k),$$

so the quantity contributed at level $n$ is simply

$$C(n)=6P_n.$$

The whole problem is therefore reduced to updating $P_n$ efficiently as $n$ increases.

Step 2: Pascal's Identity Gives the Row Transition

Using

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

we obtain

$$B_{n+1}(k)=\binom{n}{k}5^{\,n-k}=5B_n(k)+B_n(k-1).$$

This recurrence is the key structural fact. It says that every term in the next row depends only on two neighboring terms from the current row, so there is no need to rebuild the full truncated sum from scratch.

Step 3: Only the Boundary Matters

The cutoff $t_n=\pi(n)$ changes by at most one when $n$ increases. Therefore it is enough to keep:

$$P_n=\sum_{k=0}^{t_n} B_n(k),\qquad U_n=B_n(t_n),\qquad V_n=B_n(t_n+1).$$

These are, respectively, the truncated sum, the last included term, and the first excluded term. The initial state is immediate from $n=1$:

$$t_1=\pi(1)=0,\qquad P_1=1,\qquad U_1=1,\qquad V_1=0.$$

Step 4: Transition When $n+1$ Is Composite

If $n+1$ is composite, then $t_{n+1}=t_n=t$. Summing the row transition up to $k=t$ gives

$$\begin{aligned} P_{n+1} &=\sum_{k=0}^{t}\left(5B_n(k)+B_n(k-1)\right) \\ &=5\sum_{k=0}^{t}B_n(k)+\sum_{k=0}^{t-1}B_n(k) \\ &=6P_n-U_n. \end{aligned}$$

The new boundary terms are

$$U_{n+1}=5U_n+B_n(t-1),\qquad V_{n+1}=5V_n+U_n.$$

The missing term $B_n(t-1)$ is recovered from a ratio:

$$\frac{B_n(t-1)}{B_n(t)}=\frac{\binom{n-1}{t-1}5^{\,n-t}}{\binom{n-1}{t}5^{\,n-1-t}}=\frac{5t}{n-t},$$

hence

$$B_n(t-1)=U_n\cdot \frac{5t}{n-t}.$$

Step 5: Transition When $n+1$ Is Prime

If $n+1$ is prime, then the cutoff moves to the right: $t_{n+1}=t_n+1=t+1$. The next truncated sum now includes one extra boundary term:

$$P_{n+1}=6P_n+5V_n.$$

The new last included term is

$$U_{n+1}=5V_n+U_n,$$

and the new first excluded term satisfies

$$V_{n+1}=5B_n(t+2)+V_n.$$

Again a ratio gives the unseen term:

$$\frac{B_n(t+2)}{B_n(t+1)}=\frac{\binom{n-1}{t+2}5^{\,n-t-3}}{\binom{n-1}{t+1}5^{\,n-t-2}}=\frac{n-t-2}{5(t+2)},$$

$$B_n(t+2)=V_n\cdot \frac{n-t-2}{5(t+2)}.$$

If $t+2>n-1$, that term is simply $0$, exactly as the implementations handle it.

Step 6: Modular Arithmetic and Final Accumulation

Every step contributes

$$C(n)=6P_n,$$

therefore

$$S(L)\equiv \sum_{n=1}^{L} 6P_n \pmod{10^9+7}.$$

The recurrence contains divisions by $n-t$, $t+2$, and $5$. Since the modulus $10^9+7$ is prime and all these numbers are positive and strictly smaller than the modulus for the target range, each division is replaced by multiplication with a modular inverse.

Sanity Checks

The derived formulas reproduce the exact checkpoint values used by the implementations:

$$C(3)=216,\qquad C(4)=1290,\qquad C(11)=361912500,$$

$$C(24)=4727547363281250000,$$

and for the cumulative sum,

$$S(50)\equiv 832833871 \pmod{10^9+7}.$$

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They first build a primality table up to $L+1$, because only the question “is $n+1$ prime?” decides which recurrence branch to use. They also precompute modular inverses up to $L+2$, which turns the rational boundary ratios into constant-time modular multiplications. A single forward pass then maintains the truncated sum together with the two boundary terms, updates them with the prime or composite formulas above, and adds $6P_n$ to the running answer at each step.

Complexity Analysis

The sieve costs $O(L\log\log L)$ time. Precomputing modular inverses costs $O(L)$ time. The recurrence itself performs only constant-time arithmetic per $n$, so the main loop is $O(L)$. Overall the method is $O(L\log\log L)$ time and $O(L)$ memory, with only $O(1)$ recurrence state beyond the primality and inverse tables.

Footnotes and References

Problem page: https://projecteuler.net/problem=423
Binomial coefficient and Pascal's identity: Wikipedia — Binomial coefficient
Prime-counting function $\pi(n)$: Wikipedia — Prime-counting function
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse

Problem 423 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using u64 = std::uint64_t;
using boost::multiprecision::cpp_int;

constexpr u64 MOD = 1000000007ULL;

struct Options {
    int limit = 50000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& 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::stoi(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_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1ULL;
    u64 cur = base % MOD;
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = static_cast<u64>((__uint128_t)result * cur % MOD);
        }
        cur = static_cast<u64>((__uint128_t)cur * cur % MOD);
        e >>= 1ULL;
    }
    return result;
}

std::vector<bool> sieve_primes(int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n + 1), true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; 1LL * p * p <= n; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int q = p * p; q <= n; q += p) {
            is_prime[static_cast<std::size_t>(q)] = false;
        }
    }
    return is_prime;
}

std::vector<u64> modular_inverses(int n) {
    std::vector<u64> inv(static_cast<std::size_t>(n + 1), 0ULL);
    if (n >= 1) {
        inv[1] = 1ULL;
    }
    for (int i = 2; i <= n; ++i) {
        inv[static_cast<std::size_t>(i)] =
            (MOD - (MOD / static_cast<u64>(i)) * inv[static_cast<std::size_t>(MOD % static_cast<u64>(i))] % MOD) % MOD;
    }
    return inv;
}

u64 solve_mod(const int limit) {
    const std::vector<bool> is_prime = sieve_primes(limit + 1);
    const std::vector<u64> inv = modular_inverses(limit + 2);
    const u64 inv5 = mod_pow(5ULL, MOD - 2ULL);

    int t = 0;        // t = pi(n)
    u64 p = 1ULL;     // sum_{k=0..t} B_n(k)
    u64 bt = 1ULL;    // B_n(t)
    u64 bt1 = 0ULL;   // B_n(t+1)

    u64 answer = 0ULL;

    for (int n = 1; n <= limit; ++n) {
        const u64 c_n = 6ULL * p % MOD;
        answer += c_n;
        if (answer >= MOD) {
            answer -= MOD;
        }

        if (n == limit) {
            break;
        }

        const bool next_is_prime = is_prime[static_cast<std::size_t>(n + 1)];

        u64 p_next = 0ULL;
        u64 bt_next = 0ULL;
        u64 bt1_next = 0ULL;

        if (!next_is_prime) {
            p_next = (6ULL * p + MOD - bt) % MOD;

            u64 bt_minus_1 = 0ULL;
            if (t > 0) {
                bt_minus_1 = bt * (5ULL * static_cast<u64>(t) % MOD) % MOD;
                bt_minus_1 = bt_minus_1 * inv[static_cast<std::size_t>(n - t)] % MOD;
            }

            bt_next = (5ULL * bt + bt_minus_1) % MOD;
            bt1_next = (5ULL * bt1 + bt) % MOD;
        } else {
            p_next = (6ULL * p + 5ULL * bt1) % MOD;

            u64 bt2 = 0ULL;  // B_n(t+2)
            if (t + 2 <= n - 1) {
                bt2 = bt1 * static_cast<u64>(n - t - 2) % MOD;
                bt2 = bt2 * inv5 % MOD;
                bt2 = bt2 * inv[static_cast<std::size_t>(t + 2)] % MOD;
            }

            bt_next = (5ULL * bt1 + bt) % MOD;
            bt1_next = (5ULL * bt2 + bt1) % MOD;
            ++t;
        }

        p = p_next;
        bt = bt_next;
        bt1 = bt1_next;
    }

    return answer;
}

cpp_int comb_small(int n, int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    const int r = std::min(k, n - k);
    cpp_int ans = 1;
    for (int i = 1; i <= r; ++i) {
        ans *= (n - r + i);
        ans /= i;
    }
    return ans;
}

cpp_int pow_cpp_int(cpp_int base, int exp) {
    cpp_int result = 1;
    int e = exp;
    while (e > 0) {
        if (e & 1) {
            result *= base;
        }
        base *= base;
        e >>= 1;
    }
    return result;
}

cpp_int c_exact(int n) {
    const std::vector<bool> is_prime = sieve_primes(n);
    int pi_n = 0;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            ++pi_n;
        }
    }

    cpp_int sum = 0;
    for (int k = 0; k <= pi_n; ++k) {
        sum += comb_small(n - 1, k) * pow_cpp_int(cpp_int(5), n - 1 - k);
    }
    return sum * 6;
}

bool run_checkpoints() {
    if (c_exact(3) != 216) {
        std::cerr << "Checkpoint failed: C(3)\n";
        return false;
    }
    if (c_exact(4) != 1290) {
        std::cerr << "Checkpoint failed: C(4)\n";
        return false;
    }
    if (c_exact(11) != cpp_int("361912500")) {
        std::cerr << "Checkpoint failed: C(11)\n";
        return false;
    }
    if (c_exact(24) != cpp_int("4727547363281250000")) {
        std::cerr << "Checkpoint failed: C(24)\n";
        return false;
    }
    if (solve_mod(50) != 832833871ULL) {
        std::cerr << "Checkpoint failed: S(50) mod 1e9+7\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 << solve_mod(options.limit) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000000007
    limit = 50000000

    def mod_pow(base, exp):
        r = 1; base %= MOD
        while exp > 0:
            if exp & 1: r = r * base % MOD
            base = base * base % MOD
            exp >>= 1
        return r

    # Sieve primes
    is_prime = bytearray(b'\x01') * (limit + 2)
    is_prime[0] = is_prime[1] = 0
    for p in range(2, int((limit+1)**0.5) + 1):
        if is_prime[p]:
            for q in range(p*p, limit+2, p): is_prime[q] = 0

    # Modular inverses
    inv = [0] * (limit + 3)
    inv[1] = 1
    for i in range(2, limit + 3):
        inv[i] = (MOD - (MOD // i) * inv[MOD % i] % MOD) % MOD

    inv5 = mod_pow(5, MOD - 2)
    t = 0; p = 1; bt = 1; bt1 = 0
    answer = 0
    for n in range(1, limit + 1):
        c_n = 6 * p % MOD
        answer = (answer + c_n) % MOD
        if n == limit: break
        next_prime = is_prime[n + 1]
        if not next_prime:
            p_next = (6 * p + MOD - bt) % MOD
            bt_m1 = 0
            if t > 0:
                bt_m1 = bt * (5 * t % MOD) % MOD * inv[n - t] % MOD
            bt_next = (5 * bt + bt_m1) % MOD
            bt1_next = (5 * bt1 + bt) % MOD
        else:
            p_next = (6 * p + 5 * bt1) % MOD
            bt2 = 0
            if t + 2 <= n - 1:
                bt2 = bt1 * (n - t - 2) % MOD * inv5 % MOD * inv[t + 2] % MOD
            bt_next = (5 * bt1 + bt) % MOD
            bt1_next = (5 * bt2 + bt1) % MOD
            t += 1
        p, bt, bt1 = p_next, bt_next, bt1_next

    return str(answer)

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

Java

public class Euler423 {
    static final long MOD = 1000000007L;

    static long modPow(long base, long exp) {
        long result = 1;
        long cur = base % MOD;
        long e = exp;
        while (e > 0) {
            if ((e & 1) != 0) {
                result = (result * cur) % MOD;
            }
            cur = (cur * cur) % MOD;
            e >>= 1;
        }
        return result;
    }

    static boolean[] sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        if (n >= 0)
            isPrime[0] = false;
        if (n >= 1)
            isPrime[1] = false;

        for (int p = 2; p * p <= n; ++p) {
            if (!isPrime[p])
                continue;
            for (int q = p * p; q <= n; q += p) {
                isPrime[q] = false;
            }
        }
        return isPrime;
    }

    static int[] modularInverses(int n) {
        int[] inv = new int[n + 1];
        if (n >= 1)
            inv[1] = 1;
        for (int i = 2; i <= n; ++i) {
            inv[i] = (int) ((MOD - (MOD / i) * inv[(int) (MOD % i)] % MOD) % MOD);
        }
        return inv;
    }

    static String solveMod(int limit) {
        boolean[] isPrime = sievePrimes(limit + 1);
        int[] inv = modularInverses(limit + 2);
        long inv5 = modPow(5, MOD - 2);

        int t = 0;
        long p = 1;
        long bt = 1;
        long bt1 = 0;

        long answer = 0;

        for (int n = 1; n <= limit; ++n) {
            long cn = (6 * p) % MOD;
            answer += cn;
            if (answer >= MOD)
                answer -= MOD;

            if (n == limit)
                break;

            boolean nextIsPrime = isPrime[n + 1];

            long pNext, btNext, bt1Next;

            if (!nextIsPrime) {
                pNext = (6 * p + MOD - bt) % MOD;

                long btMinus1 = 0;
                if (t > 0) {
                    btMinus1 = (bt * ((5L * t) % MOD)) % MOD;
                    btMinus1 = (btMinus1 * inv[n - t]) % MOD;
                }

                btNext = (5 * bt + btMinus1) % MOD;
                bt1Next = (5 * bt1 + bt) % MOD;
            } else {
                pNext = (6 * p + 5 * bt1) % MOD;

                long bt2 = 0;
                if (t + 2 <= n - 1) {
                    bt2 = (bt1 * (n - t - 2)) % MOD;
                    bt2 = (bt2 * inv5) % MOD;
                    bt2 = (bt2 * inv[t + 2]) % MOD;
                }

                btNext = (5 * bt1 + bt) % MOD;
                bt1Next = (5 * bt2 + bt1) % MOD;
                t++;
            }

            p = pNext;
            bt = btNext;
            bt1 = bt1Next;
        }

        return Long.toString(answer);
    }

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