Problem 470: Super Ramvok

Project Euler Problem 470 Solution

EulerSolve provides an optimized solution for Project Euler Problem 470, Super Ramvok, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each dimension $d$ and roll cost $c$, the problem first asks for the best expected net profit of the ordinary Ramvok game on a fixed non-empty subset $A\subseteq\{1,\dots,d\}$. Those profits are then averaged over all subsets of the same size, and the averaged values are used in a second process whose state depends only on the current subset size. The resulting quantity is denoted $S(d,c)$, and the final target is $$F(n)=\sum_{d=4}^{n}\sum_{c=0}^{n} S(d,c).$$ The implementation verifies two checkpoints along the way: $$R(\{1,2,3,4\},0.2)=2.65,\qquad S(6,1)=208.3.$$ Mathematical Approach Fix $d\ge 4$, a cost $c$, and a non-empty subset $$A=\{v_1,\dots,v_k\}\subseteq\{1,\dots,d\},\qquad v_1<\dots<v_k.$$ The exact method implemented by the programs has three layers: optimize the ordinary game on one subset, average over all subsets with the same cardinality, then solve a one-dimensional Bellman system for the super-game. Step 1: Finite-Horizon Value on a Fixed Subset Let $h_t$ be the best expected gross reward when the available values are exactly the elements of $A$ and at most $t$ further rolls are allowed. The recurrence used by the implementation is $$h_0=0,\qquad h_t=\frac{1}{k}\sum_{i=1}^{k}\max(v_i,h_{t-1}).$$ This recurrence is an optimal-stopping equation. On the first roll we see one of the $v_i$, each with probability $1/k$....

Detailed mathematical approach

Problem Summary

For each dimension $d$ and roll cost $c$, the problem first asks for the best expected net profit of the ordinary Ramvok game on a fixed non-empty subset $A\subseteq\{1,\dots,d\}$. Those profits are then averaged over all subsets of the same size, and the averaged values are used in a second process whose state depends only on the current subset size. The resulting quantity is denoted $S(d,c)$, and the final target is

$$F(n)=\sum_{d=4}^{n}\sum_{c=0}^{n} S(d,c).$$

The implementation verifies two checkpoints along the way:

$$R(\{1,2,3,4\},0.2)=2.65,\qquad S(6,1)=208.3.$$

Mathematical Approach

Fix $d\ge 4$, a cost $c$, and a non-empty subset

$$A=\{v_1,\dots,v_k\}\subseteq\{1,\dots,d\},\qquad v_1<\dots<v_k.$$

The exact method implemented by the programs has three layers: optimize the ordinary game on one subset, average over all subsets with the same cardinality, then solve a one-dimensional Bellman system for the super-game.

Step 1: Finite-Horizon Value on a Fixed Subset

Let $h_t$ be the best expected gross reward when the available values are exactly the elements of $A$ and at most $t$ further rolls are allowed. The recurrence used by the implementation is

$$h_0=0,\qquad h_t=\frac{1}{k}\sum_{i=1}^{k}\max(v_i,h_{t-1}).$$

This recurrence is an optimal-stopping equation. On the first roll we see one of the $v_i$, each with probability $1/k$. After seeing $v_i$, we either stop immediately and keep $v_i$, or continue and replace it by the continuation value $h_{t-1}$. Therefore the contribution of outcome $v_i$ is $\max(v_i,h_{t-1})$, and averaging over all outcomes gives the displayed formula.

Step 2: Convert Gross Reward into Net Profit

If each roll costs $c$, then using exactly $t$ allowed rolls gives expected net profit

$$h_t-c t.$$

Hence the subset profit is

$$R(A,c)=\max_{t\ge 0}\bigl(h_t-c t\bigr).$$

The code exploits two useful bounds. If $c=0$, the value is simply the largest available reward, so $R(A,0)=\max(A)$. If $c\ge 1$ and $A\subseteq\{1,\dots,d\}$, then $h_t\le \max(A)\le d$, so for $t>d$ we have

$$h_t-c t\le d-t<0.$$

Since $t=0$ yields profit $0$, no horizon larger than $d$ can be optimal for the integer costs used in the main sum. That is why the exact search stops at $t=d$.

Step 3: Average over Subset Size

For fixed $d$ and $c$, define the average subset profit at level $k$ by

$$a_k(c)=\frac{1}{\binom{d}{k}}\sum_{\substack{A\subseteq\{1,\dots,d\}\\ |A|=k}} R(A,c),\qquad 1\le k\le d.$$

This is the key compression step. The super-game does not need the identity of the subset once we average over all $k$-subsets; it only needs the current size $k$. All detailed subset structure is absorbed into the scalar quantity $a_k(c)$.

Step 4: Bellman Equation for the Super-Game

Let $x_k$ denote the expected total super-value when the current subset has size $k$. The linear system in the implementations shows that the next transition depends only on whether the next uniformly chosen label is already present or absent. From a state of size $k<d$:

$$\Pr(k\to k-1)=\frac{k}{d},\qquad \Pr(k\to k+1)=1-\frac{k}{d}=\frac{d-k}{d}.$$

Conditioning on that next step gives

$$x_k=a_k(c)+\frac{k}{d}x_{k-1}+\left(1-\frac{k}{d}\right)x_{k+1},\qquad 1\le k<d.$$

The state $k=0$ is an implicit zero boundary, so no separate unknown is needed there. At the top state $k=d$, the next move must go to $d-1$, hence

$$x_d=a_d(c)+x_{d-1}.$$

Rearranging these equations produces a tridiagonal linear system, and the desired super-value for the pair $(d,c)$ is

$$S(d,c)=x_d.$$

Step 5: Final Aggregation

Once $S(d,c)$ is known for every $d$ and every integer cost $0\le c\le n$, the problem asks for their double sum:

$$F(n)=\sum_{d=4}^{n}\sum_{c=0}^{n} S(d,c).$$

This is exactly the quantity accumulated by the implementations before the final rounding step.

Worked Example: $A=\{1,2,3,4\}$ and $c=0.2$

Here $k=4$. The recurrence gives

$$h_1=\frac{1+2+3+4}{4}=2.5,$$

$$h_2=\frac{\max(1,2.5)+\max(2,2.5)+\max(3,2.5)+\max(4,2.5)}{4}=3,$$

$$h_3=\frac{3+3+3+4}{4}=3.25,$$

$$h_4=\frac{3.25+3.25+3.25+4}{4}=3.4375.$$

The corresponding profits are

$$h_1-0.2=2.3,\qquad h_2-0.4=2.6,\qquad h_3-0.6=2.65,\qquad h_4-0.8=2.6375.$$

So the optimum is attained at $t=3$, proving

$$R(\{1,2,3,4\},0.2)=2.65.$$

The second checkpoint, $S(6,1)=208.3$, verifies that the averaging stage and the tridiagonal solve are also assembled correctly.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical pipeline. For each $d$, they enumerate every non-empty subset of $\{1,\dots,d\}$ with a bitmask. The chosen bits already determine the subset in sorted order, so no expensive reordering is needed.

For each subset, the implementation computes the best net profit for every integer cost $c=0,1,\dots,n$. The $c=0$ case is handled immediately by taking the largest element of the subset. For $c\ge 1$, the implementation advances the recurrence for $h_t$ only up to $t=d$, updating the best value $h_t-c t$ for all costs in parallel.

These profits are accumulated in buckets indexed by subset size $k$. After all subsets have been processed, each bucket sum is divided by the number of subsets of that size, producing the averages $a_k(c)$.

Then, for each fixed cost $c$, the implementation builds the tridiagonal system corresponding to the Bellman equations above and solves it by forward elimination followed by back substitution. The last component gives $S(d,c)$. Finally it sums all $S(d,c)$ over $d=4,\dots,n$ and $c=0,\dots,n$, and rounds the total to the nearest integer.

Complexity Analysis

For a fixed $d$, there are $2^d-1$ non-empty subsets. Processing one subset takes $O(d)$ time to decode its elements, $O(dk)$ time to advance the recurrence over $t=1,\dots,d$, and $O(dn)$ time to update all integer costs, where $k\le d$. Thus one subset costs $O(d(d+n))$ in the worst case, and the dominant cost for that dimension is

$$O\bigl(2^d\,d(d+n)\bigr).$$

The tridiagonal solve contributes only $O(nd)$ time for that same $d$, which is much smaller than the subset enumeration. Summing over $d=4,\dots,n$ gives an overall exponential exact algorithm, dominated by the largest dimension; for $d$ on the order of $n$, this is roughly $O(2^n n^2)$. The working memory is $O(d(n+1))$ for the level averages plus $O(d)$ for the linear-system arrays, with additional copies only when parallel workers are used.

Footnotes and References

Problem page: https://projecteuler.net/problem=470
Optimal stopping and Bellman recurrences: Wikipedia — Bellman equation
Tridiagonal linear systems: Wikipedia — Tridiagonal matrix algorithm
Birth-death chains: Wikipedia — Birth-death process

Problem 470 source code

C++

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

namespace {

using u64 = std::uint64_t;

constexpr unsigned kDefaultN = 20U;
constexpr unsigned kMaxPracticalN = 24U;

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

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

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

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--n=", parsed_unsigned)) {
            options.n = parsed_unsigned;
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

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

    if (options.n < 4U) {
        std::cerr << "--n must be >= 4.\n";
        return false;
    }
    if (options.n > kMaxPracticalN) {
        std::cerr << "--n too large for this exact subset-enumeration implementation (max "
                  << kMaxPracticalN << ").\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const u64 workload_units) {
    constexpr u64 kMinUnitsForParallel = 200'000ULL;
    if (!allow_multithreading || workload_units < kMinUnitsForParallel) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

long double ramvok_profit_real_cost(const int* values, const int k, const long double cost) {
    if (k <= 0) {
        return 0.0L;
    }

    const int max_value = values[k - 1];
    if (cost <= 0.0L) {
        return static_cast<long double>(max_value);
    }

    long double h_prev = 0.0L;
    long double best = 0.0L;

    for (int t = 1; t <= 10'000; ++t) {
        long double sum = 0.0L;
        for (int i = 0; i < k; ++i) {
            const long double x = static_cast<long double>(values[i]);
            sum += (x > h_prev) ? x : h_prev;
        }

        const long double h_cur = sum / static_cast<long double>(k);
        const long double candidate = h_cur - cost * static_cast<long double>(t);
        if (candidate > best) {
            best = candidate;
        }

        // Upper bound on any future profit at horizon >= t.
        const long double future_upper = static_cast<long double>(max_value) - cost * static_cast<long double>(t);
        if (future_upper <= best + 1e-16L) {
            break;
        }

        h_prev = h_cur;
    }

    return best;
}

void subset_best_profits_integer_costs(const std::uint32_t mask,
                                       const unsigned d,
                                       const unsigned max_cost,
                                       std::vector<long double>& output) {
    int values[24];
    int k = 0;

    for (unsigned bit = 0U; bit < d; ++bit) {
        if ((mask >> bit) & 1U) {
            values[k++] = static_cast<int>(bit + 1U);
        }
    }

    output.assign(max_cost + 1U, 0.0L);
    if (k == 0) {
        return;
    }

    const int max_value = values[k - 1];
    output[0] = static_cast<long double>(max_value);

    // For integer c>=1 and max roll <= d, t>d can never beat t=0 because
    // h_t <= max_value <= d, so h_t - c*t <= d - t < 0.
    const unsigned t_limit = d;

    long double h_prev = 0.0L;
    for (unsigned t = 1U; t <= t_limit; ++t) {
        long double sum = 0.0L;
        for (int i = 0; i < k; ++i) {
            const long double x = static_cast<long double>(values[i]);
            sum += (x > h_prev) ? x : h_prev;
        }

        const long double h_cur = sum / static_cast<long double>(k);
        for (unsigned c = 1U; c <= max_cost; ++c) {
            const long double candidate =
                h_cur - static_cast<long double>(c) * static_cast<long double>(t);
            if (candidate > output[c]) {
                output[c] = candidate;
            }
        }
        h_prev = h_cur;
    }
}

std::vector<long double> compute_average_reward_by_level(const unsigned d,
                                                         const unsigned max_cost,
                                                         const bool allow_multithreading,
                                                         const unsigned requested_threads) {
    const std::uint32_t total_masks = (1U << d);
    const u64 workload = static_cast<u64>(total_masks) - 1ULL;
    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload);

    const std::size_t stride = static_cast<std::size_t>(max_cost + 1U);
    const std::size_t table_size = static_cast<std::size_t>(d + 1U) * stride;

    std::vector<long double> global_sum(table_size, 0.0L);
    std::vector<u64> global_count(static_cast<std::size_t>(d + 1U), 0ULL);

    if (threads == 1U) {
        std::vector<long double> profits;
        for (std::uint32_t mask = 1U; mask < total_masks; ++mask) {
            const unsigned k = static_cast<unsigned>(__builtin_popcount(mask));
            subset_best_profits_integer_costs(mask, d, max_cost, profits);
            ++global_count[k];
            const std::size_t base = static_cast<std::size_t>(k) * stride;
            for (unsigned c = 0U; c <= max_cost; ++c) {
                global_sum[base + c] += profits[c];
            }
        }
    } else {
        std::vector<std::thread> workers;
        workers.reserve(threads);

        std::vector<std::vector<long double>> partial_sum(threads, std::vector<long double>(table_size, 0.0L));
        std::vector<std::vector<u64>> partial_count(threads, std::vector<u64>(d + 1U, 0ULL));

        for (unsigned t = 0U; t < threads; ++t) {
            const std::uint32_t begin =
                1U + static_cast<std::uint32_t>((static_cast<u64>(total_masks - 1U) * t) / threads);
            const std::uint32_t end =
                1U + static_cast<std::uint32_t>((static_cast<u64>(total_masks - 1U) * (t + 1U)) / threads);

            workers.emplace_back([&, t, begin, end]() {
                std::vector<long double> profits;
                for (std::uint32_t mask = begin; mask < end; ++mask) {
                    const unsigned k = static_cast<unsigned>(__builtin_popcount(mask));
                    subset_best_profits_integer_costs(mask, d, max_cost, profits);
                    ++partial_count[t][k];
                    const std::size_t base = static_cast<std::size_t>(k) * stride;
                    for (unsigned c = 0U; c <= max_cost; ++c) {
                        partial_sum[t][base + c] += profits[c];
                    }
                }
            });
        }

        for (auto& worker : workers) {
            worker.join();
        }

        for (unsigned t = 0U; t < threads; ++t) {
            for (unsigned k = 1U; k <= d; ++k) {
                global_count[k] += partial_count[t][k];
            }
            for (std::size_t i = 0; i < table_size; ++i) {
                global_sum[i] += partial_sum[t][i];
            }
        }
    }

    std::vector<long double> average(table_size, 0.0L);
    for (unsigned k = 1U; k <= d; ++k) {
        const u64 cnt = global_count[k];
        const long double inv_cnt = 1.0L / static_cast<long double>(cnt);
        const std::size_t base = static_cast<std::size_t>(k) * stride;
        for (unsigned c = 0U; c <= max_cost; ++c) {
            average[base + c] = global_sum[base + c] * inv_cnt;
        }
    }

    return average;
}

std::vector<long double> solve_super_by_cost(const unsigned d,
                                             const unsigned max_cost,
                                             const std::vector<long double>& average_reward_by_level) {
    const std::size_t stride = static_cast<std::size_t>(max_cost + 1U);
    std::vector<long double> s_by_cost(max_cost + 1U, 0.0L);

    std::vector<long double> lower(d, 0.0L);
    std::vector<long double> diag(d, 0.0L);
    std::vector<long double> upper(d, 0.0L);
    std::vector<long double> rhs(d, 0.0L);

    for (unsigned c = 0U; c <= max_cost; ++c) {
        for (unsigned k = 1U; k <= d; ++k) {
            const unsigned idx = k - 1U;
            const long double a_k =
                average_reward_by_level[static_cast<std::size_t>(k) * stride + c];

            if (k == d) {
                lower[idx] = -1.0L;
                diag[idx] = 1.0L;
                upper[idx] = 0.0L;
                rhs[idx] = a_k;
                continue;
            }

            const long double p = static_cast<long double>(k) / static_cast<long double>(d);
            const long double q = 1.0L - p;

            lower[idx] = (k == 1U) ? 0.0L : -p;
            diag[idx] = 1.0L;
            upper[idx] = -q;
            rhs[idx] = a_k;
        }

        for (unsigned i = 1U; i < d; ++i) {
            const long double factor = lower[i] / diag[i - 1U];
            diag[i] -= factor * upper[i - 1U];
            rhs[i] -= factor * rhs[i - 1U];
        }

        std::vector<long double> x(d, 0.0L);
        x[d - 1U] = rhs[d - 1U] / diag[d - 1U];
        for (int i = static_cast<int>(d) - 2; i >= 0; --i) {
            x[static_cast<std::size_t>(i)] =
                (rhs[static_cast<std::size_t>(i)] -
                 upper[static_cast<std::size_t>(i)] * x[static_cast<std::size_t>(i + 1)]) /
                diag[static_cast<std::size_t>(i)];
        }

        s_by_cost[c] = x[d - 1U];
    }

    return s_by_cost;
}

bool nearly_equal(const long double a, const long double b, const long double eps) {
    return std::fabsl(a - b) <= eps;
}

bool run_checkpoints() {
    {
        const int values[] = {1, 2, 3, 4};
        const long double r = ramvok_profit_real_cost(values, 4, 0.2L);
        if (!nearly_equal(r, 2.65L, 1e-12L)) {
            std::cerr << "Checkpoint failed: R(4,0.2) expected 2.65, got "
                      << std::setprecision(18) << static_cast<double>(r) << '\n';
            return false;
        }
    }

    {
        constexpr unsigned d = 6U;
        constexpr unsigned max_cost = 1U;
        const std::vector<long double> avg =
            compute_average_reward_by_level(d, max_cost, false, 1U);
        const std::vector<long double> s = solve_super_by_cost(d, max_cost, avg);

        if (!nearly_equal(s[1], 208.3L, 1e-10L)) {
            std::cerr << "Checkpoint failed: S(6,1) expected 208.3, got "
                      << std::setprecision(18) << static_cast<double>(s[1]) << '\n';
            return false;
        }
    }

    return true;
}

long double compute_f(const unsigned n,
                      const bool allow_multithreading,
                      const unsigned requested_threads) {
    long double total = 0.0L;

    for (unsigned d = 4U; d <= n; ++d) {
        const std::vector<long double> avg =
            compute_average_reward_by_level(d, n, allow_multithreading, requested_threads);
        const std::vector<long double> s = solve_super_by_cost(d, n, avg);
        for (unsigned c = 0U; c <= n; ++c) {
            total += s[c];
        }
    }

    return total;
}

}  // namespace

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

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

    const long double f_n = compute_f(options.n, options.allow_multithreading, options.requested_threads);
    const long long rounded = std::llround(f_n);

    std::cout << "F(" << options.n << ") = " << std::setprecision(18) << static_cast<double>(f_n)
              << '\n';
    std::cout << "Rounded: " << rounded << '\n';
    std::cout << "Answer: " << rounded << '\n';

    return 0;
}

Python

import sys
from multiprocessing import Pool, cpu_count
import math

def subset_best_profits_integer_costs(mask, d, max_cost):
    values = []
    for bit in range(d):
        if (mask >> bit) & 1:
            values.append(bit + 1)
            
    k = len(values)
    output = [0.0] * (max_cost + 1)
    if k == 0:
        return output
        
    max_value = values[-1]
    output[0] = float(max_value)
    
    t_limit = d
    h_prev = 0.0
    for t in range(1, t_limit + 1):
        s = 0.0
        for x in values:
            if x > h_prev:
                s += x
            else:
                s += h_prev
                
        h_cur = s / k
        for c in range(1, max_cost + 1):
            candidate = h_cur - c * t
            if candidate > output[c]:
                output[c] = candidate
        h_prev = h_cur
        
    return output

def compute_average(args):
    begin, end, d, max_cost = args
    table_size = (d + 1) * (max_cost + 1)
    partial_sum = [0.0] * table_size
    partial_count = [0] * (d + 1)
    
    for mask in range(begin, end):
        k = mask.bit_count()
        profits = subset_best_profits_integer_costs(mask, d, max_cost)
        partial_count[k] += 1
        base = k * (max_cost + 1)
        for c in range(max_cost + 1):
            partial_sum[base + c] += profits[c]
            
    return (partial_sum, partial_count)

def compute_average_reward_by_level(d, max_cost):
    total_masks = 1 << d
    threads = max(1, cpu_count())
    
    tasks = []
    chunk = (total_masks - 1 + threads - 1) // threads
    for t in range(threads):
        begin = 1 + t * chunk
        end = 1 + min(total_masks - 1, (t + 1) * chunk)
        if begin < end:
            tasks.append((begin, end, d, max_cost))
            
    with Pool(len(tasks)) as pool:
        results = pool.map(compute_average, tasks)
        
    table_size = (d + 1) * (max_cost + 1)
    global_sum = [0.0] * table_size
    global_count = [0] * (d + 1)
    
    for p_sum, p_count in results:
        for k in range(1, d + 1):
            global_count[k] += p_count[k]
        for i in range(table_size):
            global_sum[i] += p_sum[i]
            
    average = [0.0] * table_size
    for k in range(1, d + 1):
        cnt = global_count[k]
        if cnt == 0: continue
        inv_cnt = 1.0 / cnt
        base = k * (max_cost + 1)
        for c in range(max_cost + 1):
            average[base + c] = global_sum[base + c] * inv_cnt
            
    return average

def solve_super_by_cost(d, max_cost, average_reward_by_level):
    stride = max_cost + 1
    s_by_cost = [0.0] * (max_cost + 1)
    
    for c in range(max_cost + 1):
        lower = [0.0] * d
        diag = [0.0] * d
        upper = [0.0] * d
        rhs = [0.0] * d
        
        for k in range(1, d + 1):
            idx = k - 1
            a_k = average_reward_by_level[k * stride + c]
            
            if k == d:
                lower[idx] = -1.0
                diag[idx] = 1.0
                upper[idx] = 0.0
                rhs[idx] = a_k
                continue
                
            p = float(k) / float(d)
            q = 1.0 - p
            
            lower[idx] = 0.0 if k == 1 else -p
            diag[idx] = 1.0
            upper[idx] = -q
            rhs[idx] = a_k
            
        for i in range(1, d):
            factor = lower[i] / diag[i - 1]
            diag[i] -= factor * upper[i - 1]
            rhs[i] -= factor * rhs[i - 1]
            
        x = [0.0] * d
        x[d - 1] = rhs[d - 1] / diag[d - 1]
        for i in range(d - 2, -1, -1):
            x[i] = (rhs[i] - upper[i] * x[i + 1]) / diag[i]
            
        s_by_cost[c] = x[d - 1]
        
    return s_by_cost

def compute_f(n):
    total = 0.0
    for d in range(4, n + 1):
        avg = compute_average_reward_by_level(d, n)
        s = solve_super_by_cost(d, n, avg)
        for c in range(n + 1):
            total += s[c]
    return total

def solve():
    n = 20
    ans = compute_f(n)
    return str(round(ans))

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler470 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler470.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler470_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler470 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler470 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler470 C++ bridge produced empty output.");
        }
        return parsed;
    }

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