Problem 765: Trillionaire

Problem 765 source code

C++

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

namespace {

long double optimal_success_probability(int rounds, long double p_win, long double target) {
    if (rounds < 0) {
        throw std::invalid_argument("rounds must be non-negative");
    }
    if (!(p_win > 0.0L && p_win < 1.0L)) {
        throw std::invalid_argument("p_win must be in (0, 1)");
    }
    if (target <= 1.0L) {
        return 1.0L;
    }

    const long double budget_q = 1.0L / target;

    // q[w] = Q(W = w), where Q is the fair-coin measure (0.5 / 0.5).
    std::vector<long double> q(rounds + 1, 0.0L);
    q[0] = std::ldexp(1.0L, -rounds);
    for (int w = 0; w < rounds; ++w) {
        q[w + 1] = q[w] * static_cast<long double>(rounds - w) / static_cast<long double>(w + 1);
    }

    // p[w] = P(W = w), where P uses the true win probability p_win.
    // Using p[w] = q[w] * (2(1-p))^rounds * (p/(1-p))^w keeps numbers stable.
    std::vector<long double> p(rounds + 1, 0.0L);
    const long double base = std::powl(2.0L * (1.0L - p_win), rounds);
    const long double ratio = p_win / (1.0L - p_win);
    long double ratio_pow = 1.0L;
    for (int w = 0; w <= rounds; ++w) {
        p[w] = q[w] * base * ratio_pow;
        ratio_pow *= ratio;
    }

    // Neyman-Pearson on leaf likelihood ratios: include all largest w, then possibly
    // a fractional part of one boundary layer to use Q-budget exactly.
    long double used_q = 0.0L;
    long double success_p = 0.0L;
    int boundary_w = -1;

    for (int w = rounds; w >= 0; --w) {
        const long double candidate_q = used_q + q[w];
        if (candidate_q <= budget_q + 1e-30L) {
            used_q = candidate_q;
            success_p += p[w];
        } else {
            boundary_w = w;
            break;
        }
    }

    if (boundary_w == -1) {
        return 1.0L;
    }

    long double frac = (budget_q - used_q) / q[boundary_w];
    frac = std::clamp(frac, 0.0L, 1.0L);
    success_p += frac * p[boundary_w];

    return success_p;
}

void validate() {
    const long double tol = 1e-15L;

    // 1 round, target 2: all-in gives success exactly with one win.
    {
        long double got = optimal_success_probability(1, 0.6L, 2.0L);
        long double want = 0.6L;
        if (std::fabsl(got - want) > tol) {
            throw std::runtime_error("Validation failed: N=1, T=2");
        }
    }

    // 2 rounds, target 4: need two net doublings, probability p^2.
    {
        long double got = optimal_success_probability(2, 0.6L, 4.0L);
        long double want = 0.36L;
        if (std::fabsl(got - want) > tol) {
            throw std::runtime_error("Validation failed: N=2, T=4");
        }
    }

    // If target <= starting capital, success is guaranteed.
    {
        long double got = optimal_success_probability(1000, 0.6L, 1.0L);
        long double want = 1.0L;
        if (std::fabsl(got - want) > tol) {
            throw std::runtime_error("Validation failed: T<=1 guarantee");
        }
    }
}

}  // namespace

int main() {
    validate();

    const int rounds = 1000;
    const long double p_win = 0.6L;
    const long double target = 1.0e12L;

    long double ans = optimal_success_probability(rounds, p_win, target);

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

Python

import decimal

def optimal_success_probability(rounds, p_win, target):
    decimal.getcontext().prec = 50
    
    if target <= 1:
        return decimal.Decimal(1)
        
    budget_q = decimal.Decimal(1) / decimal.Decimal(str(target))
    
    q = [decimal.Decimal(0)] * (rounds + 1)
    q[0] = decimal.Decimal(1) / (decimal.Decimal(2) ** rounds)
    for w in range(rounds):
        q[w + 1] = q[w] * decimal.Decimal(rounds - w) / decimal.Decimal(w + 1)
        
    p = [decimal.Decimal(0)] * (rounds + 1)
    dp_win = decimal.Decimal(str(p_win))
    one_minus_p = decimal.Decimal(1) - dp_win
    base = (decimal.Decimal(2) * one_minus_p) ** rounds
    ratio = dp_win / one_minus_p
    ratio_pow = decimal.Decimal(1)
    
    for w in range(rounds + 1):
        p[w] = q[w] * base * ratio_pow
        ratio_pow *= ratio
        
    used_q = decimal.Decimal(0)
    success_p = decimal.Decimal(0)
    boundary_w = -1
    eps = decimal.Decimal('1e-30')
    
    for w in range(rounds, -1, -1):
        candidate_q = used_q + q[w]
        if candidate_q <= budget_q + eps:
            used_q = candidate_q
            success_p += p[w]
        else:
            boundary_w = w
            break
            
    if boundary_w == -1:
        return decimal.Decimal(1)
        
    frac = (budget_q - used_q) / q[boundary_w]
    if frac < decimal.Decimal(0): frac = decimal.Decimal(0)
    if frac > decimal.Decimal(1): frac = decimal.Decimal(1)
    
    success_p += frac * p[boundary_w]
    return success_p

def solve():
    ans = optimal_success_probability(1000, 0.6, 1000000000000.0)
    return f"{ans:.10f}"

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

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

public class Euler765 {
    static final MathContext MC = new MathContext(50, RoundingMode.HALF_UP);

    static BigDecimal optimalSuccessProbability(int rounds, double pWin, double target) {
        if (target <= 1.0) {
            return BigDecimal.ONE;
        }

        BigDecimal budgetQ = BigDecimal.ONE.divide(new BigDecimal(Double.toString(target)), MC);

        BigDecimal[] q = new BigDecimal[rounds + 1];
        BigDecimal two = new BigDecimal("2");
        q[0] = BigDecimal.ONE.divide(two.pow(rounds, MC), MC);

        for (int w = 0; w < rounds; ++w) {
            BigDecimal rw = new BigDecimal(rounds - w);
            BigDecimal w1 = new BigDecimal(w + 1);
            q[w + 1] = q[w].multiply(rw, MC).divide(w1, MC);
        }

        BigDecimal pWinDec = new BigDecimal(Double.toString(pWin));
        BigDecimal oneMinusP = BigDecimal.ONE.subtract(pWinDec, MC);
        BigDecimal base = two.multiply(oneMinusP, MC).pow(rounds, MC);
        BigDecimal ratio = pWinDec.divide(oneMinusP, MC);

        BigDecimal[] p = new BigDecimal[rounds + 1];
        BigDecimal ratioPow = BigDecimal.ONE;
        for (int w = 0; w <= rounds; ++w) {
            p[w] = q[w].multiply(base, MC).multiply(ratioPow, MC);
            ratioPow = ratioPow.multiply(ratio, MC);
        }

        BigDecimal usedQ = BigDecimal.ZERO;
        BigDecimal successP = BigDecimal.ZERO;
        int boundaryW = -1;
        BigDecimal eps = new BigDecimal("1e-30");

        for (int w = rounds; w >= 0; --w) {
            BigDecimal candidateQ = usedQ.add(q[w], MC);
            if (candidateQ.compareTo(budgetQ.add(eps, MC)) <= 0) {
                usedQ = candidateQ;
                successP = successP.add(p[w], MC);
            } else {
                boundaryW = w;
                break;
            }
        }

        if (boundaryW == -1) {
            return BigDecimal.ONE;
        }

        BigDecimal frac = budgetQ.subtract(usedQ, MC).divide(q[boundaryW], MC);
        if (frac.compareTo(BigDecimal.ZERO) < 0)
            frac = BigDecimal.ZERO;
        if (frac.compareTo(BigDecimal.ONE) > 0)
            frac = BigDecimal.ONE;

        successP = successP.add(frac.multiply(p[boundaryW], MC), MC);
        return successP;
    }

    public static String solve() {
        BigDecimal ans = optimalSuccessProbability(1000, 0.6, 1000000000000.0);
        return ans.setScale(10, RoundingMode.HALF_UP).toPlainString();
    }

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