Problem 968: 5D Summation

Project Euler Problem 968 Solution

EulerSolve provides an optimized solution for Project Euler Problem 968, 5D Summation, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each $n=0,1,\dots,99$, the problem builds a 10-component bound vector from the recurrence $$a_0=1,\qquad a_1=7,\qquad a_m \equiv 7a_{m-1}+a_{m-2}^2 \pmod{10^9+7}.$$ The block $(a_{10n},a_{10n+1},\dots,a_{10n+9})$ is attached to the 10 edges of the complete graph $K_5$ in the order $$ (0,1),(0,2),(0,3),(0,4),(1,2),(1,3),(1,4),(2,3),(2,4),(3,4). $$ For one such block $U=(U_e)_{e\in E(K_5)}$, we must evaluate $$P(U)=\sum_{\substack{x_0,\dots,x_4\ge 0\\ x_u+x_v\le U_{uv}\ \forall\,0\le u<v\le 4}} 2^{x_0}3^{x_1}5^{x_2}7^{x_3}11^{x_4}\pmod{10^9+7},$$ and the final answer is $$\sum_{n=0}^{99} P(a_{10n},a_{10n+1},\dots,a_{10n+9}) \pmod{10^9+7},$$ with the ten numbers interpreted in the fixed edge order above. The direct five-dimensional summation is only realistic for tiny bounds; the real inputs are large, so the solution rewrites the inequalities bit by bit and performs a digit DP on carry states. Mathematical Approach The key observation is that every constraint is of the form $x_u+x_v\le U_{uv}$. In binary, such an inequality can be checked locally: at bit $k$, only the two chosen bits, the $k$-th bit of the bound, and one small carry value matter. Once that is done simultaneously on all 10 edges, the original 5D sum becomes a finite-state weighted DP....

Detailed mathematical approach

Problem Summary

For each $n=0,1,\dots,99$, the problem builds a 10-component bound vector from the recurrence

$$a_0=1,\qquad a_1=7,\qquad a_m \equiv 7a_{m-1}+a_{m-2}^2 \pmod{10^9+7}.$$

The block $(a_{10n},a_{10n+1},\dots,a_{10n+9})$ is attached to the 10 edges of the complete graph $K_5$ in the order

$$ (0,1),(0,2),(0,3),(0,4),(1,2),(1,3),(1,4),(2,3),(2,4),(3,4). $$

For one such block $U=(U_e)_{e\in E(K_5)}$, we must evaluate

$$P(U)=\sum_{\substack{x_0,\dots,x_4\ge 0\\ x_u+x_v\le U_{uv}\ \forall\,0\le u<v\le 4}} 2^{x_0}3^{x_1}5^{x_2}7^{x_3}11^{x_4}\pmod{10^9+7},$$

and the final answer is

$$\sum_{n=0}^{99} P(a_{10n},a_{10n+1},\dots,a_{10n+9}) \pmod{10^9+7},$$

with the ten numbers interpreted in the fixed edge order above. The direct five-dimensional summation is only realistic for tiny bounds; the real inputs are large, so the solution rewrites the inequalities bit by bit and performs a digit DP on carry states.

Mathematical Approach

The key observation is that every constraint is of the form $x_u+x_v\le U_{uv}$. In binary, such an inequality can be checked locally: at bit $k$, only the two chosen bits, the $k$-th bit of the bound, and one small carry value matter. Once that is done simultaneously on all 10 edges, the original 5D sum becomes a finite-state weighted DP.

The weighted sum on the 10 edges of $K_5$

The five variables are the vertices of $K_5$, and every edge $\{u,v\}$ contributes one upper bound $U_{uv}$. The term attached to a feasible point $(x_0,\dots,x_4)$ is multiplicative:

$$2^{x_0}3^{x_1}5^{x_2}7^{x_3}11^{x_4}=\prod_{i=0}^{4} p_i^{x_i},\qquad (p_0,p_1,p_2,p_3,p_4)=(2,3,5,7,11).$$

So the job is not just to count feasible 5-tuples, but to sum a geometric weight over all of them.

The 100 independent instances generated by the recurrence

The outer problem does not use one bound vector but 100 of them. The recurrence produces 1000 values, and each consecutive block of 10 values becomes one complete set of pairwise bounds. Because each block is independent of the others once the sequence has been generated, the total answer is simply the sum of 100 values of $P(U)$.

This is an important structural simplification: all difficult mathematics is concentrated in evaluating one generic quantity $P(U)$. After that, the final Euler sum is just an outer aggregation over 100 instances.

An edge-wise carry invariant

Fix one edge $e=\{u,v\}$. Write

$$x_i=\sum_{k\ge 0} b_{i,k}2^k,\qquad b_{i,k}\in\{0,1\}.$$

After bits $0,1,\dots,k-1$ have been decided, define the remaining upper parts

$$y_i^{(k)}=\left\lfloor \frac{x_i}{2^k}\right\rfloor.$$

For the edge $e$, keep a nonnegative integer $c_{e,k}$ such that the unfinished higher bits must satisfy

$$y_u^{(k)}+y_v^{(k)}+c_{e,k}\le \left\lfloor \frac{U_e}{2^k}\right\rfloor.$$

This carry is the amount of excess forced upward by the lower bits already chosen. At the start, no lower bits have been processed, so $c_{e,0}=0$ for every edge.

The local transition rule

Let $u_{e,k}\in\{0,1\}$ be the $k$-th binary digit of $U_e$. Since

$$y_u^{(k)}=b_{u,k}+2y_u^{(k+1)},\qquad y_v^{(k)}=b_{v,k}+2y_v^{(k+1)},$$

substituting into the invariant gives

$$b_{u,k}+b_{v,k}+c_{e,k}+2\bigl(y_u^{(k+1)}+y_v^{(k+1)}\bigr)\le u_{e,k}+2\left\lfloor \frac{U_e}{2^{k+1}}\right\rfloor.$$

The smallest carry that makes the same statement true one level higher is therefore

$$c_{e,k+1}=\max\left(0,\left\lceil \frac{b_{u,k}+b_{v,k}+c_{e,k}-u_{e,k}}{2}\right\rceil\right).$$

This is the core recurrence of the whole method. It says that an edge transition depends only on the two chosen vertex bits, the incoming carry on that edge, and the current bit of the bound.

Why only three carry values are needed

For one edge, the pair bit sum $b_{u,k}+b_{v,k}$ is in $\{0,1,2\}$, the bound bit $u_{e,k}$ is in $\{0,1\}$, and if the incoming carry is already in $\{0,1,2\}$, then

$$b_{u,k}+b_{v,k}+c_{e,k}-u_{e,k}\in\{-1,0,1,2,3,4\}.$$

Applying the formula above shows that the next carry also lies in $\{0,1,2\}$. So each edge has exactly three possible carry states, and the full DP state space is

$$\{0,1,2\}^{10},\qquad |\{0,1,2\}^{10}|=3^{10}=59049.$$

One 5-bit choice updates all 10 edges at once

At bit $k$, choose the 5 binary digits $(b_{0,k},\dots,b_{4,k})$. It is convenient to regard them as a 5-bit mask $m\in\{0,1\}^5$. That mask immediately determines the pair sums $b_{u,k}+b_{v,k}$ on all 10 edges, so it also determines the next carry vector once the current state and the bound bits are known.

If the carry vector at level $k$ is $c^{(k)}\in\{0,1,2\}^{10}$, then for a fixed bound vector $U$ the transition map $T_k$ is defined edgewise by

$$\bigl(T_k(c^{(k)},m)\bigr)_e=\max\left(0,\left\lceil \frac{s_e(m)+c^{(k)}_e-u_{e,k}}{2}\right\rceil\right),$$

where $s_e(m)\in\{0,1,2\}$ is the pair bit sum on edge $e$ induced by the mask.

The weight factor splits by bit

The exponential weight also factorizes over the binary digits. Since $x_i=\sum_k b_{i,k}2^k$,

$$\prod_{i=0}^{4} p_i^{x_i}=\prod_{k\ge 0}\prod_{i=0}^{4} p_i^{\,b_{i,k}2^k}.$$

So one mask $m$ chosen at bit $k$ contributes the multiplicative factor

$$W_k(m)=\prod_{i=0}^{4}\left(p_i^{2^k}\right)^{b_{i,k}}\pmod{10^9+7}.$$

This is why the implementations precompute the 32 mask weights for each bit position: once those numbers are available, each DP transition only performs a table lookup and one modular multiplication.

Worked example: one edge transition

Suppose that on one particular edge the incoming carry is $c_{e,k}=1$, the bound bit is $u_{e,k}=0$, and the current mask sets $b_{u,k}=1$ and $b_{v,k}=0$. Then

$$c_{e,k+1}=\max\left(0,\left\lceil \frac{1+0+1-0}{2}\right\rceil\right)=1.$$

So the excess created by the lower bits is still large enough that one unit must be carried to the next binary level. If instead both endpoint bits were 1, then

$$c_{e,k+1}=\max\left(0,\left\lceil \frac{1+1+1-0}{2}\right\rceil\right)=2,$$

which shows why carry value 2 really can occur and must be represented in the state space.

As a concrete full-instance checkpoint, when all ten bounds are equal to 2, the digit DP gives

$$P(2,2,\dots,2)=7120,$$

exactly matching brute force. The method is not approximate; it is an exact reformulation of the original summation.

The DP recurrence and the terminal condition

Let $D_k(c)$ be the total weight of all assignments to the first $k$ bits whose carry vector after those bits is $c$. Then

$$D_{k+1}(c')=\sum_{\substack{c\in\{0,1,2\}^{10}\\ m\in\{0,1\}^5\\ T_k(c,m)=c'}} D_k(c)\,W_k(m)\pmod{10^9+7},$$

with initial condition $D_0(\mathbf 0)=1$ and $D_0(c)=0$ for $c\ne \mathbf 0$.

The implementations process 31 bit positions. That is enough because every bound produced by the recurrence lies below $10^9+7<2^{30}$, and one extra zero bit safely flushes any residual carry. A 5-tuple $(x_0,\dots,x_4)$ is feasible exactly when all carries are zero after the final round, so

$$P(U)=D_{31}(\mathbf 0).$$

How the Code Works

Shared preprocessing

The C++, Python, and Java implementations all rely on the same mathematical ingredients. They precompute the 1000-term recurrence sequence, decode every one of the $3^{10}$ carry states into its 10 ternary digits, build the pair-bit table for all 32 masks on the 10 edges, and precompute the 32 weight factors for each of the 31 bit positions.

The Python implementation is a thin wrapper around the compiled solver, so the heavy computation is the same carry-based digit DP used by the compiled versions. The Java implementation reproduces the same recurrence and transition logic directly.

Evaluating one bound vector

For a single instance $U$, the implementation first extracts the 31 binary columns of the 10 bounds. Then it runs the DP with two rolling layers: the current layer stores the values $D_k(c)$, and the next layer accumulates $D_{k+1}(c')$. The process starts from the all-zero carry vector and tries all 32 masks at each step.

A practical optimization is to track only active states, meaning states whose current DP value is nonzero. In the worst case the state space has 59049 elements, but in many layers only a fraction of them are reachable, so this sparse traversal saves time without changing the recurrence.

Validation and outer aggregation

The reference implementation checks the DP against brute force on small cases. In particular, it verifies the exact values $P(2,\dots,2)=7120$ and $P(1,2,\dots,10)=799809376$, and it also tests random small bound vectors. Those checks confirm that the carry recurrence reproduces the original five-dimensional sum exactly.

After that, the 100 problem instances are independent: for each $n$, take the block $(a_{10n},\dots,a_{10n+9})$, evaluate $P(U)$, and add it to the total modulo $10^9+7$. The compiled implementations parallelize this outer loop because each block can be solved independently.

Complexity Analysis

For one bound vector, the worst-case DP cost is

$$O\!\left(31\cdot 3^{10}\cdot 2^5\cdot 10\right),$$

because there are 31 bit levels, at most $3^{10}$ carry states, 32 masks per state, and 10 edges to update in one transition. The actual runtime is usually lower because the active-state optimization avoids scanning states with zero weight.

The memory usage is $O(3^{10})$ for the two rolling DP layers, plus the precomputed state-decoding, mask, and weight tables. The outer summation over 100 instances multiplies the total amount of work by 100, but those 100 calls are independent and therefore parallel-friendly.

A naive brute-force approach would require iterating over five nested loops up to bounds near $10^9$, which is completely infeasible. The digit-DP succeeds because it replaces huge numeric ranges by a tiny per-edge carry alphabet $\{0,1,2\}$.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=968
Complete graph: Wikipedia - Complete graph
Dynamic programming: Wikipedia - Dynamic programming
Binary number system: Wikipedia - Binary number
Adder and carry propagation: Wikipedia - Adder
Modular exponentiation: Wikipedia - Modular exponentiation

Problem 968 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <random>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u32 MOD = 1'000'000'007U;
constexpr int VARS = 5;
constexpr int EDGES = 10;
constexpr int BITS = 31;
constexpr std::array<int, EDGES> POW3 = {
    1, 3, 9, 27, 81, 243, 729, 2187, 6561, 19683,
};
constexpr int STATE_COUNT = 59049;  // 3^10
constexpr std::array<std::array<int, 2>, EDGES> EDGE_ENDPOINTS = {{
    {0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 2},
    {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4},
}};
constexpr std::array<u32, VARS> BASES = {2U, 3U, 5U, 7U, 11U};

std::array<std::array<std::uint8_t, EDGES>, STATE_COUNT> CARRY_DIGITS{};
std::array<std::array<std::uint8_t, EDGES>, 32> PAIR_BITS{};
std::array<std::array<u32, 32>, BITS> BIT_FACTORS{};

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

void build_tables() {
    for (int s = 0; s < STATE_COUNT; ++s) {
        int t = s;
        for (int e = 0; e < EDGES; ++e) {
            CARRY_DIGITS[static_cast<std::size_t>(s)][static_cast<std::size_t>(e)] =
                static_cast<std::uint8_t>(t % 3);
            t /= 3;
        }
    }

    for (int mask = 0; mask < 32; ++mask) {
        for (int e = 0; e < EDGES; ++e) {
            const int u = EDGE_ENDPOINTS[static_cast<std::size_t>(e)][0];
            const int v = EDGE_ENDPOINTS[static_cast<std::size_t>(e)][1];
            const int bit_u = (mask >> u) & 1;
            const int bit_v = (mask >> v) & 1;
            PAIR_BITS[static_cast<std::size_t>(mask)][static_cast<std::size_t>(e)] =
                static_cast<std::uint8_t>(bit_u + bit_v);
        }
    }

    std::array<std::array<u32, BITS>, VARS> powers{};
    for (int i = 0; i < VARS; ++i) {
        powers[static_cast<std::size_t>(i)][0] = BASES[static_cast<std::size_t>(i)] % MOD;
        for (int k = 1; k < BITS; ++k) {
            const u64 v = powers[static_cast<std::size_t>(i)][static_cast<std::size_t>(k - 1)];
            powers[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)] =
                static_cast<u32>((v * v) % MOD);
        }
    }

    for (int k = 0; k < BITS; ++k) {
        for (int mask = 0; mask < 32; ++mask) {
            u64 value = 1U;
            for (int i = 0; i < VARS; ++i) {
                if ((mask >> i) & 1) {
                    value = (value * powers[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)]) % MOD;
                }
            }
            BIT_FACTORS[static_cast<std::size_t>(k)][static_cast<std::size_t>(mask)] =
                static_cast<u32>(value);
        }
    }
}

u32 compute_p(const std::array<u32, EDGES>& bounds) {
    std::array<std::array<std::uint8_t, EDGES>, BITS> bound_bits{};
    for (int k = 0; k < BITS; ++k) {
        for (int e = 0; e < EDGES; ++e) {
            bound_bits[static_cast<std::size_t>(k)][static_cast<std::size_t>(e)] =
                static_cast<std::uint8_t>((bounds[static_cast<std::size_t>(e)] >> k) & 1U);
        }
    }

    std::vector<u32> cur(static_cast<std::size_t>(STATE_COUNT), 0U);
    std::vector<u32> nxt(static_cast<std::size_t>(STATE_COUNT), 0U);
    std::vector<int> active;
    std::vector<int> next_active;
    active.reserve(4096);
    next_active.reserve(4096);
    active.push_back(0);
    cur[0] = 1U;

    for (int k = 0; k < BITS; ++k) {
        next_active.clear();
        const auto& bits = bound_bits[static_cast<std::size_t>(k)];

        for (int state : active) {
            const u32 state_value = cur[static_cast<std::size_t>(state)];
            if (state_value == 0U) {
                continue;
            }
            cur[static_cast<std::size_t>(state)] = 0U;

            const auto& carry = CARRY_DIGITS[static_cast<std::size_t>(state)];

            for (int mask = 0; mask < 32; ++mask) {
                const auto& pair = PAIR_BITS[static_cast<std::size_t>(mask)];
                int next_state = 0;
                for (int e = 0; e < EDGES; ++e) {
                    const int t = static_cast<int>(pair[static_cast<std::size_t>(e)]) +
                                  static_cast<int>(carry[static_cast<std::size_t>(e)]) -
                                  static_cast<int>(bits[static_cast<std::size_t>(e)]);
                    const int next_carry = (t <= 0) ? 0 : ((t + 1) >> 1);
                    next_state += next_carry * POW3[static_cast<std::size_t>(e)];
                }

                const u64 add = (static_cast<u64>(state_value) *
                                 BIT_FACTORS[static_cast<std::size_t>(k)][static_cast<std::size_t>(mask)]) %
                                MOD;
                u32& dst = nxt[static_cast<std::size_t>(next_state)];
                if (dst == 0U) {
                    next_active.push_back(next_state);
                }
                const u64 sum = static_cast<u64>(dst) + add;
                dst = static_cast<u32>(sum >= MOD ? sum - MOD : sum);
            }
        }

        active.swap(next_active);
        cur.swap(nxt);
    }

    return cur[0];
}

u32 brute_p(const std::array<u32, EDGES>& bounds) {
    const int ub_a = static_cast<int>(std::min({bounds[0], bounds[1], bounds[2], bounds[3]}));
    const int ub_b = static_cast<int>(std::min({bounds[0], bounds[4], bounds[5], bounds[6]}));
    const int ub_c = static_cast<int>(std::min({bounds[1], bounds[4], bounds[7], bounds[8]}));
    const int ub_d = static_cast<int>(std::min({bounds[2], bounds[5], bounds[7], bounds[9]}));
    const int ub_e = static_cast<int>(std::min({bounds[3], bounds[6], bounds[8], bounds[9]}));

    std::vector<u32> p2(static_cast<std::size_t>(ub_a + 1), 1U);
    std::vector<u32> p3(static_cast<std::size_t>(ub_b + 1), 1U);
    std::vector<u32> p5(static_cast<std::size_t>(ub_c + 1), 1U);
    std::vector<u32> p7(static_cast<std::size_t>(ub_d + 1), 1U);
    std::vector<u32> p11(static_cast<std::size_t>(ub_e + 1), 1U);
    for (int i = 1; i <= ub_a; ++i) {
        p2[static_cast<std::size_t>(i)] = static_cast<u32>((2ULL * p2[static_cast<std::size_t>(i - 1)]) % MOD);
    }
    for (int i = 1; i <= ub_b; ++i) {
        p3[static_cast<std::size_t>(i)] = static_cast<u32>((3ULL * p3[static_cast<std::size_t>(i - 1)]) % MOD);
    }
    for (int i = 1; i <= ub_c; ++i) {
        p5[static_cast<std::size_t>(i)] = static_cast<u32>((5ULL * p5[static_cast<std::size_t>(i - 1)]) % MOD);
    }
    for (int i = 1; i <= ub_d; ++i) {
        p7[static_cast<std::size_t>(i)] = static_cast<u32>((7ULL * p7[static_cast<std::size_t>(i - 1)]) % MOD);
    }
    for (int i = 1; i <= ub_e; ++i) {
        p11[static_cast<std::size_t>(i)] = static_cast<u32>((11ULL * p11[static_cast<std::size_t>(i - 1)]) % MOD);
    }

    u32 total = 0U;
    for (int a = 0; a <= ub_a; ++a) {
        for (int b = 0; b <= ub_b; ++b) {
            if (static_cast<u32>(a + b) > bounds[0]) {
                continue;
            }
            for (int c = 0; c <= ub_c; ++c) {
                if (static_cast<u32>(a + c) > bounds[1] || static_cast<u32>(b + c) > bounds[4]) {
                    continue;
                }
                for (int d = 0; d <= ub_d; ++d) {
                    if (static_cast<u32>(a + d) > bounds[2] || static_cast<u32>(b + d) > bounds[5] ||
                        static_cast<u32>(c + d) > bounds[7]) {
                        continue;
                    }
                    for (int e = 0; e <= ub_e; ++e) {
                        if (static_cast<u32>(a + e) > bounds[3] || static_cast<u32>(b + e) > bounds[6] ||
                            static_cast<u32>(c + e) > bounds[8] || static_cast<u32>(d + e) > bounds[9]) {
                            continue;
                        }
                        u64 term = p2[static_cast<std::size_t>(a)];
                        term = (term * p3[static_cast<std::size_t>(b)]) % MOD;
                        term = (term * p5[static_cast<std::size_t>(c)]) % MOD;
                        term = (term * p7[static_cast<std::size_t>(d)]) % MOD;
                        term = (term * p11[static_cast<std::size_t>(e)]) % MOD;
                        total = static_cast<u32>(total + term);
                        if (total >= MOD) {
                            total -= MOD;
                        }
                    }
                }
            }
        }
    }

    return total;
}

bool run_checkpoints() {
    {
        std::array<u32, EDGES> bounds{};
        bounds.fill(2U);
        if (compute_p(bounds) != 7120U) {
            std::cerr << "Checkpoint failed: P(2,...,2) != 7120" << '\n';
            return false;
        }
    }
    {
        const std::array<u32, EDGES> bounds = {1U, 2U, 3U, 4U, 5U, 6U, 7U, 8U, 9U, 10U};
        if (compute_p(bounds) != 799809376U) {
            std::cerr << "Checkpoint failed: P(1..10) mismatch" << '\n';
            return false;
        }
    }
    {
        std::mt19937 rng(968U);
        std::uniform_int_distribution<int> dist(0, 4);
        for (int it = 0; it < 20; ++it) {
            std::array<u32, EDGES> bounds{};
            for (int e = 0; e < EDGES; ++e) {
                bounds[static_cast<std::size_t>(e)] = static_cast<u32>(dist(rng));
            }
            const u32 fast = compute_p(bounds);
            const u32 brute = brute_p(bounds);
            if (fast != brute) {
                std::cerr << "Checkpoint failed on random small case" << '\n';
                return false;
            }
        }
    }
    return true;
}

u32 solve() {
    std::vector<u32> a(1000U, 0U);
    a[0] = 1U;
    a[1] = 7U;
    for (int i = 2; i < 1000; ++i) {
        const u64 term = (7ULL * a[static_cast<std::size_t>(i - 1)] +
                          static_cast<u64>(a[static_cast<std::size_t>(i - 2)]) *
                              a[static_cast<std::size_t>(i - 2)]) %
                         MOD;
        a[static_cast<std::size_t>(i)] = static_cast<u32>(term);
    }

    std::vector<u32> q(100U, 0U);
    std::atomic<int> next_idx{0};

    unsigned workers = std::thread::hardware_concurrency();
    if (workers == 0U) {
        workers = 4U;
    }
    workers = std::min(workers, 100U);

    std::vector<std::thread> pool;
    pool.reserve(workers);
    for (unsigned t = 0; t < workers; ++t) {
        pool.emplace_back([&]() {
            while (true) {
                const int n = next_idx.fetch_add(1, std::memory_order_relaxed);
                if (n >= 100) {
                    break;
                }
                std::array<u32, EDGES> bounds{};
                for (int e = 0; e < EDGES; ++e) {
                    bounds[static_cast<std::size_t>(e)] =
                        a[static_cast<std::size_t>(10 * n + e)];
                }
                q[static_cast<std::size_t>(n)] = compute_p(bounds);
            }
        });
    }
    for (auto& th : pool) {
        th.join();
    }

    u32 total = 0U;
    for (u32 value : q) {
        total = static_cast<u32>(total + value);
        if (total >= MOD) {
            total -= MOD;
        }
    }
    return total;
}

}  // namespace

int main() {
    build_tables();
    if (!run_checkpoints()) {
        return 1;
    }
    std::cout << solve() << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler968 {

    static final int MOD = 1000000007;
    static final int VARS = 5;
    static final int EDGES = 10;
    static final int BITS = 31;
    static final int[] POW3 = { 1, 3, 9, 27, 81, 243, 729, 2187, 6561, 19683 };
    static final int STATE_COUNT = 59049;
    static final int[][] EDGE_ENDPOINTS = {
            { 0, 1 }, { 0, 2 }, { 0, 3 }, { 0, 4 }, { 1, 2 },
            { 1, 3 }, { 1, 4 }, { 2, 3 }, { 2, 4 }, { 3, 4 }
    };
    static final int[] BASES = { 2, 3, 5, 7, 11 };

    static byte[][] CARRY_DIGITS = new byte[STATE_COUNT][EDGES];
    static byte[][] PAIR_BITS = new byte[32][EDGES];
    static int[][] BIT_FACTORS = new int[BITS][32];

    static void buildTables() {
        for (int s = 0; s < STATE_COUNT; ++s) {
            int t = s;
            for (int e = 0; e < EDGES; ++e) {
                CARRY_DIGITS[s][e] = (byte) (t % 3);
                t /= 3;
            }
        }

        for (int mask = 0; mask < 32; ++mask) {
            for (int e = 0; e < EDGES; ++e) {
                int u = EDGE_ENDPOINTS[e][0];
                int v = EDGE_ENDPOINTS[e][1];
                int bitU = (mask >> u) & 1;
                int bitV = (mask >> v) & 1;
                PAIR_BITS[mask][e] = (byte) (bitU + bitV);
            }
        }

        int[][] powers = new int[VARS][BITS];
        for (int i = 0; i < VARS; ++i) {
            powers[i][0] = BASES[i] % MOD;
            for (int k = 1; k < BITS; ++k) {
                long v = powers[i][k - 1];
                powers[i][k] = (int) ((v * v) % MOD);
            }
        }

        for (int k = 0; k < BITS; ++k) {
            for (int mask = 0; mask < 32; ++mask) {
                long value = 1;
                for (int i = 0; i < VARS; ++i) {
                    if (((mask >> i) & 1) == 1) {
                        value = (value * powers[i][k]) % MOD;
                    }
                }
                BIT_FACTORS[k][mask] = (int) value;
            }
        }
    }

    static int computeP(int[] bounds) {
        byte[][] boundBits = new byte[BITS][EDGES];
        for (int k = 0; k < BITS; ++k) {
            for (int e = 0; e < EDGES; ++e) {
                boundBits[k][e] = (byte) ((bounds[e] >> k) & 1);
            }
        }

        int[] cur = new int[STATE_COUNT];
        int[] nxt = new int[STATE_COUNT];
        int[] active = new int[STATE_COUNT];
        int[] nextActive = new int[STATE_COUNT];
        int activeCount = 1;
        int nextActiveCount = 0;

        active[0] = 0;
        cur[0] = 1;

        for (int k = 0; k < BITS; ++k) {
            nextActiveCount = 0;
            byte[] bits = boundBits[k];
            int[] bitFactorsK = BIT_FACTORS[k];

            for (int aidx = 0; aidx < activeCount; ++aidx) {
                int state = active[aidx];
                int stateValue = cur[state];
                if (stateValue == 0)
                    continue;
                cur[state] = 0;

                byte[] carry = CARRY_DIGITS[state];

                for (int mask = 0; mask < 32; ++mask) {
                    byte[] pair = PAIR_BITS[mask];
                    int nextState = 0;
                    for (int e = 0; e < EDGES; ++e) {
                        int t = pair[e] + carry[e] - bits[e];
                        int nextCarry = (t <= 0) ? 0 : ((t + 1) >> 1);
                        if (nextCarry != 0) {
                            nextState += nextCarry * POW3[e];
                        }
                    }

                    long add = ((long) stateValue * bitFactorsK[mask]) % MOD;
                    int dst = nxt[nextState];
                    if (dst == 0) {
                        nextActive[nextActiveCount++] = nextState;
                    }
                    dst += add;
                    if (dst >= MOD)
                        dst -= MOD;
                    nxt[nextState] = dst;
                }
            }

            int[] tempActive = active;
            active = nextActive;
            nextActive = tempActive;
            activeCount = nextActiveCount;

            int[] tempCur = cur;
            cur = nxt;
            nxt = tempCur;
        }

        return cur[0];
    }

    public static String solve() {
        buildTables();

        int[] a = new int[1000];
        a[0] = 1;
        a[1] = 7;
        for (int i = 2; i < 1000; ++i) {
            long term = (7L * a[i - 1] + (long) a[i - 2] * a[i - 2]) % MOD;
            a[i] = (int) term;
        }

        int total = IntStream.range(0, 100).parallel().map(n -> {
            int[] bounds = new int[EDGES];
            for (int e = 0; e < EDGES; ++e) {
                bounds[e] = a[10 * n + e];
            }
            return computeP(bounds);
        }).reduce(0, (acc, val) -> {
            int sum = acc + val;
            return sum >= MOD ? sum - MOD : sum;
        });

        return Integer.toString(total);
    }

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