Problem 267: Billionaire

View on Project Euler

Project Euler Problem 267 Solution

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

Problem Summary A gambler starts with capital \(1\). In each of \(n\) fair tosses, the gambler bets the same fixed fraction \(f\) of the current capital. A win triples the stake portion, so total capital is multiplied by \((1+2f)\); a loss destroys the stake portion, so total capital is multiplied by \((1-f)\). For the original Project Euler problem, \(n=1000\) and the target is \(T=10^9\). The C++ solution keeps both quantities as parameters, so the same derivation works for any \(n\ge 1\) and \(T>1\). Mathematical Approach 1. Wealth After Exactly \(k\) Wins If exactly \(k\) of the \(n\) tosses are wins, then the order of wins and losses does not matter, because multiplication is commutative. The final wealth is therefore $$W_k(f)=(1+2f)^k(1-f)^{n-k}.$$ So for each possible win count \(k\), the problem is a one-variable optimization over the admissible range \(0\le f<1\). 2. Why the Code Optimizes the Logarithm Directly maximizing \(W_k(f)\) is equivalent to maximizing its logarithm, because \(\ln(x)\) is strictly increasing on \(x>0\). The code therefore works with $$g_k(f)=\ln W_k(f)=k\ln(1+2f)+(n-k)\ln(1-f).$$ This removes huge intermediate numbers and turns the product into a sum. It also makes the derivative easy to analyze. 3....

Detailed mathematical approach

Problem Summary

We start with £1 and toss a fair coin \(1000\) times. Before every toss we bet the same fixed proportion \(f\) of our current capital, where \(0\le f\le1\). Heads returns the bet together with twice the bet as a win; tails loses the bet. In the statement's example \(f=1/4\): we bet £0.25, heads brings us to £1.5, we then bet £0.375, and tails on the second toss leaves £1.125. We choose \(f\) so that the probability of having at least £1,000,000,000 after the \(1000\) tosses is as large as possible, and must give this largest probability, computed exactly and rounded to \(12\) digits after the decimal point. [1]

The whole solution in one sentence: the final capital depends only on how many tosses show heads, and more heads never hurt, so every \(f\) succeeds exactly when the number of heads reaches a threshold; the best \(f\) is one with the smallest threshold, which a short calculation with derivatives finds, and the answer is the probability that \(1000\) fair tosses show at least that many heads.

Mathematical Approach

Throughout, \(n\) is the number of tosses and \(T\gt1\) the target in pounds, so \(n=1000\) and \(T=10^9\) in the problem, while smaller values serve as examples that can be checked by hand.

1. Only the number of heads matters

If the capital before a toss is \(C\), the bet is \(fC\). Heads adds \(2fC\), so the capital becomes \((1+2f)\,C\); tails takes away \(fC\), so it becomes \((1-f)\,C\). Every toss multiplies the capital by one of two fixed factors, and since the order of the factors in a product does not matter, \(k\) heads and \(n-k\) tails end with

$$W_k(f)=(1+2f)^k\,(1-f)^{n-k}$$

pounds, whatever the order of the tosses. For \(f=1/4\) the factors are \(3/2\) and \(3/4\): heads and then tails gives \(1\to3/2\to9/8\), the £1.125 of the statement, and tails and then heads gives \(1\to3/4\to9/8\) as well. The example already shows the key feature of the game: one head and one tail together leave us richer than before. In general they multiply the capital by \((1+2f)(1-f)=1+f-2f^2\), which exceeds \(1\) for \(0\lt f\lt1/2\), so the target may be reachable with fewer heads than tails.

2. More heads never hurt

For \(0\lt f\lt1\), turning one tail into a head multiplies the final capital by

$$\frac{W_{k+1}(f)}{W_k(f)}=\frac{1+2f}{1-f}\gt1,$$

so \(W_k(f)\) grows strictly with \(k\); for \(f=0\) it is constant, and for \(f=1\) it is \(0\) until \(k=n\). Hence, for every fixed \(f\), the numbers of heads that reach the target form an unbroken block \(k(f),k(f)+1,\dots,n\), where the threshold \(k(f)\) is the smallest \(k\) with \(W_k(f)\ge T\); if not even \(n\) heads suffice, we put \(k(f)=n+1\).

Let \(K\) be the number of heads. Each of the \(2^n\) sequences of heads and tails has probability \(2^{-n}\), and the number of sequences with exactly \(k\) heads is the binomial coefficient \(\binom nk\), the number of ways to choose the \(k\) tosses that show heads; for example \(\binom42=6\). So \(K\) has the binomial distribution [3]

$$\Pr[K=k]=\binom{n}{k}2^{-n},\qquad k=0,1,\dots,n.$$

The probability of becoming a billionaire with the fraction \(f\) is therefore \(\Pr[K\ge k(f)]\).

3. The best fraction has the smallest threshold

Every \(\Pr[K=k]\) is positive, so \(\Pr[K\ge k]\) becomes strictly smaller each time \(k\) grows by one. Maximising the probability over \(f\) is therefore the same as minimising the threshold \(k(f)\). For every \(k\) let

$$M_k=\max_{0\le f\le1}W_k(f)$$

be the largest capital that \(k\) heads can produce, and let \(k_{\min}\) be the smallest \(k\) with \(M_k\ge T\). By Section 2, \(W_{k+1}(f)\ge W_k(f)\) for every \(f\), so \(M_{k+1}\ge M_k\): once a count of heads reaches the target, every larger count does too. Then two facts hold:

  • every \(f\) has \(k(f)\ge k_{\min}\), because \(W_k(f)\ge T\) for some \(k\lt k_{\min}\) would give \(M_k\ge W_k(f)\ge T\); so no \(f\) succeeds with a probability above \(\Pr[K\ge k_{\min}]\);
  • a fraction \(f\) at which \(W_{k_{\min}}\) takes its maximum \(M_{k_{\min}}\ge T\) has \(k(f)=k_{\min}\), so by Section 2 it succeeds with probability exactly \(\Pr[K\ge k_{\min}]\).

So the requested chance is the binomial tail

$$P=\Pr[K\ge k_{\min}]=2^{-n}\sum_{k=k_{\min}}^{n}\binom{n}{k}.$$

The optimal \(f\) itself is never needed: every \(f\) with \(W_{k_{\min}}(f)\ge T\) is optimal, and all of them give the same \(P\). They form an interval, because \(W_{k_{\min}}\) rises up to its maximiser and falls after it (Section 4); it has positive length unless \(M_{k_{\min}}=T\) exactly, and then the maximiser of \(W_{k_{\min}}\) is the only optimal fraction. If no \(k\) reaches the target, which happens only for \(T\gt3^n\), the answer is \(0\).

4. The best fraction for a fixed number of heads

To find \(M_k\) we maximise the logarithm of \(W_k(f)\), which has the same maximiser because \(\ln\) is increasing:

$$g_k(f)=k\ln(1+2f)+(n-k)\ln(1-f).$$

For \(k\lt n\) and \(0\le f\lt1\) its derivative

$$g_k'(f)=\frac{2k}{1+2f}-\frac{n-k}{1-f}$$

is strictly decreasing, since the first fraction never grows and the second one grows with \(f\); in other words \(g_k''\lt0\), so \(g_k\) is strictly concave. [4] Setting \(g_k'(f)=0\) gives \(2k(1-f)=(n-k)(1+2f)\), that is

$$f_k^{*}=\frac{3k-n}{2n}.$$

Left of \(f_k^{*}\) the derivative is positive and right of it negative, so \(g_k\) rises up to \(f_k^{*}\) and falls after it. Three cases cover every \(k\):

  • \(3k\le n\): then \(f_k^{*}\le0\), \(g_k\) falls on the whole interval, and the best choice is \(f=0\) with \(M_k=1\lt T\). So few heads can never succeed, whatever the bet.
  • \(n/3\lt k\lt n\): then \(0\lt f_k^{*}\lt1\), where the second inequality is equivalent to \(k\lt n\), and the maximum is at \(f_k^{*}\).
  • \(k=n\): then \(W_n(f)=(1+2f)^n\) grows with \(f\), so \(M_n=3^n\), reached at \(f=1\).

In the middle case \(1+2f_k^{*}=3k/n\) and \(1-f_k^{*}=3(n-k)/(2n)\), which gives the closed form

$$M_k=\Bigl(\frac{3k}{n}\Bigr)^{k}\Bigl(\frac{3(n-k)}{2n}\Bigr)^{n-k},$$

and with the convention \(0^0=1\) it also gives \(M_n=3^n\). For a whole-number target \(T\), as in the problem, the test \(M_k\ge T\) can even be made exactly, by comparing the integers \((3k)^k(3(n-k))^{n-k}\) and \(T\,n^k(2n)^{n-k}\); Section 9 uses this.

5. The connection with the Kelly criterion

Write \(p=k/n\) for the proportion of heads. Then \(g_k(f)/n\) equals \(p\ln(1+2f)+(1-p)\ln(1-f)\), the expected logarithm of the factor by which one round multiplies the capital, for a bet that gains twice the stake with probability \(p\) and loses the stake otherwise. For \(p\gt1/3\) the fraction that maximises it,

$$f_k^{*}=p-\frac{1-p}{2},$$

is the Kelly bet for that wager; [2] for \(p\le1/3\) this formula is not positive and the best choice is not to bet, as in Section 4. A gambler who knew in advance that exactly \(k\) tosses will show heads would thus bet the Kelly fraction for \(p=k/n\). For the real coin, \(p=1/2\), Kelly's fraction is \(1/4\), the very fraction of the statement's example. It maximises the long-run growth rate, not the chance of reaching one fixed target; Section 7 shows in an example that a smaller bet can do better.

6. Worked example: four tosses and a target of 4

Take \(n=4\) and \(T=4\), one of the checks built into the C++ program. The counts \(k=0\) and \(k=1\) satisfy \(3k\le4\), so they never succeed. For the others:

  • \(k=2\): \(f_2^{*}=1/4\) and \(M_2=(3/2)^2(3/4)^2=81/64\), about \(1.27\), less than \(4\);
  • \(k=3\): \(f_3^{*}=5/8\) and \(M_3=(9/4)^3\cdot3/8=2187/512\), about \(4.27\), at least \(4\);
  • \(k=4\): \(M_4=3^4=81\).

So \(k_{\min}=3\), and the answer for this small game is \(\Pr[K\ge3]=(4+1)/16=5/16\). The optimal fractions are those with \(W_3(f)=(1+2f)^3(1-f)\ge4\); numerically they run from \(1/2\) to about \(0.7338\). The left end is exact: with \(f=1/2\), three heads and one tail give \(2^3\cdot\tfrac12=4\), exactly the target. Such ties are why the statement insists on exact computation, and "at least" counts this one as a success. Every other \(f\) needs four heads or fails altogether, so it succeeds with probability \(1/16\) or \(0\).

The other built-in check, \(n=2\) and \(T=2\), is even shorter: \(M_1=(3/2)(3/4)=9/8\lt2\), so only two heads can succeed, and the answer is \(1/4\).

7. A larger example: one hundred tosses and a target of 10

For \(n=100\) and \(T=10\) the closed form gives \(M_{43}\approx7.54\) and \(M_{44}\approx11.61\), so \(k_{\min}=44\): as Section 1 predicted, the target can be reached with fewer heads than tails. The best probability is \(\Pr[K\ge44]\approx0.903326\), and the optimal fractions form the interval from about \(0.1195\) to \(0.2009\) around \(f_{44}^{*}=0.16\). The probability as a function of \(f\) is a staircase; a few of its steps:

  • \(f=0.05\): threshold \(51\), probability about \(0.4602\);
  • \(f=0.10\), and also Kelly's \(f=0.25\): threshold \(45\), probability about \(0.8644\);
  • \(f=0.15\) and \(f=0.20\): threshold \(44\), probability about \(0.9033\);
  • \(f=0.30\): threshold \(46\), probability about \(0.8159\);
  • \(f=0.40\): threshold \(49\), probability about \(0.6178\).

A small bet makes the capital grow too slowly, a large bet loses too much on every tail, and the best bets lie in between. With exactly \(50\) heads, Kelly's \(f=1/4\) ends with \((9/8)^{50}\approx361\) pounds, more than any other fraction, yet it needs one head more than the optimal fractions to reach \(10\).

8. Adding up the binomial tail

The probabilities \(P_k=\binom nk2^{-n}\) are best computed one from the other. Since \(\binom{n}{k+1}=\binom nk\cdot\frac{n-k}{k+1}\), [3]

$$P_0=2^{-n},\qquad P_{k+1}=P_k\cdot\frac{n-k}{k+1}.$$

For \(n=4\) this gives \(1/16\), \(4/16\), \(6/16\), \(4/16\), \(1/16\), and the tail from \(k=3\) is again \(5/16\). The recurrence avoids huge intermediate numbers: \(1000!\) has \(2568\) digits, far beyond the largest double, about \(1.8\cdot10^{308}\), whereas every \(P_k\) for \(n=1000\) lies between \(P_0=2^{-1000}\approx9.33\cdot10^{-302}\) and \(P_{500}\approx0.0252\). The start value is close to the bottom of the double range: the smallest positive normal number is \(2^{-1022}\approx2.2\cdot10^{-308}\), the smallest subnormal one \(2^{-1074}\). [5] For \(n\ge1075\) the number \(2^{-n}\) is at most half of \(2^{-1074}\) and rounds to zero, and so does every \(P_k\) computed from it.

9. Exact decisions and rounding

The programs use floating-point numbers, so two steps need care. The first is the test \(M_k\ge T\), which the programs make by comparing \(g_k(f_k^{*})\) with \(\ln T\approx20.723\). For \(n=1000\) and \(T=10^9\), the difference \(g_k(f_k^{*})-\ln T\) is about \(+0.311\) at \(k=k_{\min}\) and about \(-0.106\) at \(k=k_{\min}-1\), while the double-precision values of \(g_k(f_k^{*})\) at these two \(k\) are off by less than \(2\cdot10^{-14}\). So the scan cannot pick a wrong threshold here, as the exact integer test of Section 4 confirms. Without such a margin the logarithms can fail: for \(n=8\) and \(T=M_3=553584375/536870912\), which three heads reach exactly at \(f=1/16\), the true threshold is \(3\) and the answer \(219/256\), but in double precision the rounded \(g_3(f_3^{*})\) falls just below the rounded \(\ln T\), so the scan reports the threshold \(4\) and \(163/256\). The second step is the rounding to \(12\) digits: the tail computed in double precision by left-to-right addition differs from the exact fraction by about \(4\cdot10^{-16}\), and the exact value lies more than \(2\cdot10^{-13}\) away from the nearest point where the twelfth digit would change.

How the Code Works

The C++ program [7] keeps its parameters in Options: \(1000\) tosses, the target \(10^9\) and a flag for the built-in checks. parse_arguments changes them with --tosses=, --target= and --skip-checkpoints, and rejects fewer than one toss or a target that is not above \(1\); but the digit loop of its helper parse_int_after_prefix has no overflow check, so an oversized toss count is misread (with GCC and Clang, 4294967298 becomes 2). minimum_wins_needed carries out Sections 3 and 4: it scans \(k=0,1,\dots,n\), takes \(n\ln3\) for \(k=n\), and otherwise computes \(f_k^{*}\) and takes \(g_k(f_k^{*})\) when \(0\lt f_k^{*}\lt1\) and \(0\) when \(f_k^{*}\le0\); it returns the first \(k\) whose value reaches \(\ln T\), or \(n+1\) if none does. binomial_tail_probability implements Section 8: it returns \(1\) for a threshold of at most \(0\) and \(0\) for one above \(n\); otherwise it fills a vector with \(P_0,\dots,P_n\), starting from std::ldexp(1.0L, -n), and adds \(P_{k_{\min}},\dots,P_n\) in increasing order of \(k\). solve_probability chains the two functions. Before the main computation, run_checkpoints checks the two examples of Section 6, \(1/4\) for \(n=T=2\) and \(5/16\) for \(n=T=4\), with a tolerance of \(10^{-15}\) and stops the program with exit code 2 if one fails. The result is converted to double and printed with \(12\) decimals.

The C++ program computes in long double. With GCC on x86 processors this is an 80-bit format with a far wider range, but with Microsoft Visual C++, and on 64-bit ARM under Windows and macOS, it is the same as double. [6] On such systems the underflow of Section 8 applies: a run with --tosses=1100, for example, prints 0.000000000000. The default \(n=1000\) is not affected.

The Python and Java versions have the same three functions (in Java binomialTailProbability and so on), compute in double and have \(n=1000\) and \(T=10^9\) built in, without options or checks. Python starts with math.ldexp(1.0, -n) and adds the tail with sum over a slice; Java starts with Math.pow(2.0, -n) and formats with Locale.US, so that the decimal separator is a point in every locale.

Complexity and Verification

The scan evaluates at most \(n+1\) values with two logarithms each, and the tail needs \(n\) multiplications and divisions and at most \(n+1\) additions, so the time is \(O(n)\). The vector of probabilities takes \(O(n)\) memory; keeping only the current \(P_k\) would need \(O(1)\). For \(n=1000\) each program finishes within a tenth of a second on a desktop computer, most of it start-up time.

The C++ program checks itself on the two examples of Section 6. For this article the following checks were also made; they are not part of the published solutions.

  • A brute force with exact fractions follows the capital toss by toss through all \(2^n\) sequences, for \(401\) equally spaced fractions \(f\) and the fractions \(f_k^{*}\). For every \(n\le10\) and every target \(T\in\{2,3,5,10,30,100,1000\}\), \(70\) cases in all, it finds the best probability given by Section 3.
  • The floating-point scan of minimum_wins_needed finds the threshold of the exact integer test of Section 4 for every \(n\le1000\) with \(T=10^3\) and with \(T=10^9\), and for every \(n\le300\) with \(T\in\{2,3,10,100,10^6\}\), \(3500\) pairs in all.
  • For \(n=1000\) and \(T=10^9\), the C++ program (built with GCC and with Clang) and the Python and Java programs print the same twelve digits: the exact fraction, correctly rounded.
  • With --tosses and --target the C++ program reproduces the examples of Sections 6 and 7, for instance 0.903326047752 for \(n=100\) and \(T=10\).

Footnotes and References

The references below supply the background facts used above. The reduction to a threshold number of heads, the closed form of \(M_k\) and the checks of exactness are derived explicitly in this article.

  1. Project Euler 267 — Billionaire. The official statement: starting with £1, a fixed proportion f of the capital is bet on each of 1000 fair coin tosses; heads returns double the bet, tails loses it. With f chosen to maximise the chance of ending with at least £1,000,000,000, the task asks for that chance, rounded to 12 digits after the decimal point.
  2. Wikipedia — Kelly criterion. For a bet that gains the fraction g of the stake with probability p and loses the fraction l with probability q = 1 − p, maximising the expected logarithm of wealth, p·log(1 + fg) + q·log(1 − fl), gives the Kelly bet f = p/l − q/g; with g = 2 and l = 1 this is the fraction p − q/2 of Section 5.
  3. Wikipedia — Binomial distribution. The number of successes in n independent trials with success probability p equals k with probability C(n, k) p^k (1 − p)^(n − k); consecutive probabilities have the ratio (n − k)p / ((k + 1)(1 − p)), which for p = 1/2 gives the recurrence of Section 8.
  4. Wikipedia — Concave function. A twice-differentiable function with negative second derivative is strictly concave, and any local maximum of a concave function is also a global maximum.
  5. Wikipedia — Double-precision floating-point format. The format's largest normal number is about 1.8 × 10^308, its smallest positive normal number 2^−1022 ≈ 2.2 × 10^−308 and its smallest positive subnormal number 2^−1074 ≈ 4.9 × 10^−324.
  6. Wikipedia — Long double. With the GNU C Compiler, long double is the 80-bit extended precision format on x86 processors; Microsoft Visual C++ makes it a synonym for double, and so does 64-bit ARM (AArch64) under Windows and macOS.
  7. C++ — Euler267.cpp. The linked immutable C++ revision contains minimum_wins_needed, binomial_tail_probability and the checks of run_checkpoints. It is the implementation record for this article.

Problem 267 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

struct Options {
    int tosses = 1000;
    long double target = 1.0e9L;
    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;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_ld_after_prefix(const std::string& arg, const std::string& prefix, long double& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    char* end_ptr = nullptr;
    value = std::strtold(tail.c_str(), &end_ptr);
    return end_ptr != nullptr && *end_ptr == '\0';
}

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, "--tosses=", options.tosses) ||
            parse_ld_after_prefix(arg, "--target=", options.target)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.tosses >= 1 && options.target > 1.0L;
}

long double binomial_tail_probability(const int n, const int k_min) {
    if (k_min <= 0) {
        return 1.0L;
    }
    if (k_min > n) {
        return 0.0L;
    }

    std::vector<long double> prob(static_cast<std::size_t>(n + 1), 0.0L);
    prob[0] = std::ldexp(1.0L, -n);  // 2^-n
    for (int k = 0; k < n; ++k) {
        prob[static_cast<std::size_t>(k + 1)] =
            prob[static_cast<std::size_t>(k)] * static_cast<long double>(n - k) /
            static_cast<long double>(k + 1);
    }

    long double sum = 0.0L;
    for (int k = k_min; k <= n; ++k) {
        sum += prob[static_cast<std::size_t>(k)];
    }
    return sum;
}

int minimum_wins_needed(const int n, const long double target) {
    const long double ln_target = std::log(target);
    for (int k = 0; k <= n; ++k) {
        long double best_gain = -1.0e300L;
        if (k == n) {
            best_gain = static_cast<long double>(n) * std::log(3.0L);
        } else {
            const long double f_star =
                (3.0L * static_cast<long double>(k) - static_cast<long double>(n)) /
                (2.0L * static_cast<long double>(n));
            if (f_star > 0.0L && f_star < 1.0L) {
                best_gain =
                    static_cast<long double>(k) * std::log(1.0L + 2.0L * f_star) +
                    static_cast<long double>(n - k) * std::log(1.0L - f_star);
            } else if (f_star <= 0.0L) {
                best_gain = 0.0L;
            }
        }
        if (best_gain >= ln_target) {
            return k;
        }
    }
    return n + 1;
}

long double solve_probability(const int tosses, const long double target) {
    const int k_min = minimum_wins_needed(tosses, target);
    return binomial_tail_probability(tosses, k_min);
}

bool run_checkpoints() {
    // For two tosses and target 2.0, requiring exactly two wins gives probability 1/4.
    {
        const long double p = solve_probability(2, 2.0L);
        if (std::fabsl(p - 0.25L) > 1e-15L) {
            std::cerr << "Checkpoint failed for tosses=2 target=2.0" << '\n';
            return false;
        }
    }
    // For four tosses and target 4.0, at least three wins are required at optimum.
    {
        const long double p = solve_probability(4, 4.0L);
        const long double expected = 5.0L / 16.0L;
        if (std::fabsl(p - expected) > 1e-15L) {
            std::cerr << "Checkpoint failed for tosses=4 target=4.0" << '\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;
    }
    const long double ans = solve_probability(options.tosses, options.target);
    std::cout << std::fixed << std::setprecision(12) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

import math

def binomial_tail_probability(n, k_min):
    if k_min <= 0:
        return 1.0
    if k_min > n:
        return 0.0
        
    prob = [0.0] * (n + 1)
    # math.ldexp is available in math package but pow(2.0, -n) is also fine
    prob[0] = math.ldexp(1.0, -n)
    
    for k in range(n):
        prob[k + 1] = prob[k] * (n - k) / (k + 1)
        
    return sum(prob[k_min : n + 1])

def minimum_wins_needed(n, target):
    ln_target = math.log(target)
    for k in range(n + 1):
        best_gain = -1.0e300
        if k == n:
            best_gain = n * math.log(3.0)
        else:
            f_star = (3.0 * k - n) / (2.0 * n)
            if 0.0 < f_star < 1.0:
                best_gain = k * math.log(1.0 + 2.0 * f_star) + (n - k) * math.log(1.0 - f_star)
            elif f_star <= 0.0:
                best_gain = 0.0
                
        if best_gain >= ln_target:
            return k
            
    return n + 1

def solve_probability(tosses, target):
    k_min = minimum_wins_needed(tosses, target)
    return binomial_tail_probability(tosses, k_min)

def solve():
    ans = solve_probability(1000, 1.0e9)
    return f"{ans:.12f}"

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

Java

import java.util.Locale;

public class Euler267 {
    static double binomialTailProbability(int n, int kMin) {
        if (kMin <= 0)
            return 1.0;
        if (kMin > n)
            return 0.0;

        double[] prob = new double[n + 1];
        prob[0] = Math.pow(2.0, -n);

        for (int k = 0; k < n; k++) {
            prob[k + 1] = prob[k] * (n - k) / (k + 1);
        }

        double sum = 0.0;
        for (int k = kMin; k <= n; k++) {
            sum += prob[k];
        }
        return sum;
    }

    static int minimumWinsNeeded(int n, double target) {
        double lnTarget = Math.log(target);
        for (int k = 0; k <= n; k++) {
            double bestGain = -1.0e300;
            if (k == n) {
                bestGain = n * Math.log(3.0);
            } else {
                double fStar = (3.0 * k - n) / (2.0 * n);
                if (fStar > 0.0 && fStar < 1.0) {
                    bestGain = k * Math.log(1.0 + 2.0 * fStar) + (n - k) * Math.log(1.0 - fStar);
                } else if (fStar <= 0.0) {
                    bestGain = 0.0;
                }
            }
            if (bestGain >= lnTarget) {
                return k;
            }
        }
        return n + 1;
    }

    static double solveProbability(int tosses, double target) {
        int kMin = minimumWinsNeeded(tosses, target);
        return binomialTailProbability(tosses, kMin);
    }

    public static String solve() {
        double ans = solveProbability(1000, 1.0e9);
        return String.format(Locale.US, "%.12f", ans);
    }

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