Problem 286: Scoring Probabilities

Project Euler Problem 286 Solution

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

Problem Summary There are $N=50$ shots. On shot $x$, the hit probability depends on an unknown parameter $q$: $$p_x(q)=1-\frac{x}{q},\qquad x=1,2,\dots,50.$$ We must find the value of $q$ for which the probability of getting exactly $20$ hits is $$0.02.$$ The final numeric Project Euler answer is intentionally omitted here. Mathematical Approach 1) Basic constraints on $q$. To keep every hit probability between $0$ and $1$, we need $$q>50.$$ The expected number of hits is $$\mu(q)=\sum_{x=1}^{50}\left(1-\frac{x}{q}\right)=50-\frac{1275}{q}.$$ So even when $q\downarrow 50$, the mean is already above $24.5$. Since the target score is only $20$, we are searching for a left-tail probability. 2) This is a Poisson-binomial distribution. The 50 shots are independent, but the success probabilities are different. Therefore the total number of hits is not binomial; it is a Poisson-binomial random variable. A useful generating-function view is $$G(z)=\prod_{x=1}^{50}\bigl((1-p_x)+p_x z\bigr).$$ The coefficient of $z^s$ in $G(z)$ is exactly the probability of finishing with $s$ hits. 3) Exact DP recurrence. Let $P_x(s)$ be the probability of having exactly $s$ hits after the first $x$ shots. For shot $x$, either we miss and stay at $s$, or we hit and come from $s-1$....

Detailed mathematical approach

Problem Summary

In every practice run Barbara takes one shot from each of the distances $x=1,2,\dots,50$. A shot from distance $x$ scores one point with probability $1-x/q$, where $q\gt50$ is a real constant, so each of these probabilities lies strictly between $0$ and $1$. The shots are independent of each other, which the statement takes for granted. Her records say that she scores a total of exactly $20$ points with probability precisely $2\%$. We must find $q$ and give it rounded to $10$ decimal places. ^[1]

The whole solution in one sentence: for any trial value of $q$, adding the shots one at a time gives the exact chance of scoring $20$ points; this chance falls steadily as $q$ grows, so halving an interval that contains the right $q$ again and again pins $q$ down far beyond ten decimal places.

Mathematical Approach

Throughout, $p_x=1-x/q$ is the chance of scoring from distance $x$, $S$ is the number of points in one run, and $F(q)=\Pr[S=20]$. We must solve $F(q)=0.02$ with $q\gt50$.

1. Three shots by hand

Keep only the distances $1$, $2$, $3$ and take $q=5$. The three shots score with probabilities $4/5$, $3/5$, $2/5$ and miss with probabilities $1/5$, $2/5$, $3/5$. By independence, a pattern of hits and misses has the product of its factors as its probability, and a score collects all patterns with that many hits:

$0$ points, only miss–miss–miss: $\frac15\cdot\frac25\cdot\frac35=\frac{6}{125}$;
$1$ point, hit–miss–miss, miss–hit–miss or miss–miss–hit: $\frac{24+9+4}{125}=\frac{37}{125}$;
$2$ points, with the miss at the third, second or first shot: $\frac{36+16+6}{125}=\frac{58}{125}$;
$3$ points, only hit–hit–hit: $\frac45\cdot\frac35\cdot\frac25=\frac{24}{125}$.

The four probabilities add up to $1$, as they must, and the most likely score is $2$. For the real problem such a list is hopeless: $20$ points out of $50$ shots can happen in $\binom{50}{20}=47\,129\,212\,243\,960$ patterns.

2. The distribution of the score

The number of successes in independent yes-or-no trials with possibly different success probabilities has the Poisson binomial distribution. ^[2] Its probability generating function is

$$G(z)=\prod_{x=1}^{50}\bigl(1-p_x+p_x z\bigr),$$

and the coefficient of $z^k$ in $G(z)$ is $\Pr[S=k]$. Indeed, multiplying out the product means choosing from every factor either the miss term $1-p_x$ or the hit term $p_x z$; each choice is one pattern as in Section 1, its coefficient is the probability of that pattern, and the power of $z$ counts its hits. For the three shots:

$$\begin{aligned}&\Bigl(\tfrac15+\tfrac45z\Bigr)\Bigl(\tfrac25+\tfrac35z\Bigr)\Bigl(\tfrac35+\tfrac25z\Bigr)\\&\qquad=\tfrac{1}{125}\bigl(6+37z+58z^2+24z^3\bigr).\end{aligned}$$

The mean of $S$ is the sum of the $p_x$. ^[2] Since $1+2+\dots+50=1275$,

$$\mathbb E[S]=\sum_{x=1}^{50}\Bigl(1-\frac{x}{q}\Bigr)=50-\frac{1275}{q},$$

which exceeds $24.5$ for every $q\gt50$: a score of $20$ is always below average.

3. Adding the shots one at a time

Let $P_x(s)$ be the probability of exactly $s$ points after the first $x$ shots, starting from $P_0(0)=1$ and $P_0(s)=0$ for $s\ne0$. Exactly $s$ points after $x$ shots arise in two ways that exclude each other: $s$ points after $x-1$ shots and a miss at shot $x$, or $s-1$ points and a hit. Shot $x$ is independent of the earlier shots, so the law of total probability gives

$$\begin{aligned}P_x(s)&=(1-p_x)\,P_{x-1}(s)\\&\quad+p_x\,P_{x-1}(s-1),\end{aligned}$$

where $P_{x-1}(-1)=0$; in the terms of Section 2, this multiplies the polynomial of the first $x-1$ shots by the factor of shot $x$. For the three shots the rows are $\tfrac15,\tfrac45$, then $\tfrac{2}{25},\tfrac{11}{25},\tfrac{12}{25}$, and then the four values of Section 1. For instance, $\tfrac{11}{25}=\tfrac25\cdot\tfrac45+\tfrac35\cdot\tfrac15$: a point at the first shot followed by a miss, or a miss followed by a point. In the same way $\tfrac{58}{125}=\tfrac35\cdot\tfrac{12}{25}$${}+\tfrac25\cdot\tfrac{11}{25}$. Shot $x$ produces $x+1$ new values, so the $50$ shots need $1275$ updates with two products each and $50$ single products for $s=0$, instead of $2^{50}$ patterns. This is the direct convolution algorithm for the Poisson binomial distribution. ^[2]

4. The chance as a function of q

Exactly $20$ points means exactly $30$ misses. Adding the probabilities over all sets $A$ of $30$ missed distances gives

$$F(q)=\sum_{|A|=30}\ \prod_{x\in A}\frac{x}{q}\prod_{x\notin A}\Bigl(1-\frac{x}{q}\Bigr).$$

Every term is $q^{-30}$ times a polynomial in $1/q$ of degree $20$, so $F$ is a polynomial in $1/q$ of degree at most $50$; in particular it is continuous and differentiable for $q\gt0$. For the three shots, one point means two misses, and the same kind of sum gives $\Pr[S=1]=11/q^2-18/q^3$, which is $55/125-18/125=37/125$ at $q=5$. The recurrence of Section 3, run in $80$-digit arithmetic, gives these values:

$F(q)$ tends to $0.0412276961\ldots$ as $q$ decreases to $50$, where the last shot could no longer score;
$F(52)\approx0.0240157724$ and $F(53)\approx0.0180966946$;
$F(100)\approx2.5\cdot10^{-8}$, and $F(q)\to0$ as $q\to\infty$, because every $p_x$ tends to $1$ and the score to $50$.

5. The chance falls as q grows

Claim. $F$ is strictly decreasing for $q\gt50$. Together with Section 4 this shows that $F(q)=0.02$ has exactly one solution with $q\gt50$, and that it lies between $52$ and $53$.

One shot at a time. Fix a distance $x$ and let $S_x$ be the score of the other $49$ shots. Conditioning on shot $x$ as in Section 3,

$$\begin{aligned}F(q)&=p_x\Pr[S_x=19]\\&\quad+(1-p_x)\Pr[S_x=20].\end{aligned}$$

The distribution of $S_x$ does not involve $p_x$, so if only $p_x$ changes, $F$ changes at the rate $\Pr[S_x=19]-\Pr[S_x=20]$. As $q$ grows, $p_x=1-x/q$ grows at the rate $x/q^2$, and the chain rule gives

$$\begin{aligned}F'(q)&=\sum_{x=1}^{50}\frac{x}{q^2}\,d_x,\\ d_x&=\Pr[S_x=19]-\Pr[S_x=20].\end{aligned}$$

It is therefore enough to show that every $d_x$ is negative: for the other $49$ shots, $20$ points must be more likely than $19$.

The shape of a Poisson binomial distribution. Let $T$ count the successes of $n$ independent trials with success probabilities $r_1,\dots,r_n$ strictly between $0$ and $1$, and let $a_i=r_i/(1-r_i)$ be their odds. Taking the factor $\prod_i(1-r_i)$ out of the probability of every pattern gives $\Pr[T=k]=c\,e_k$, where $c=\prod_i(1-r_i)$ and $e_k$, the elementary symmetric polynomial, is the sum of all products of $k$ different odds, with $e_0=1$. For the three shots the odds are $4$, $3/2$, $2/3$ and $c=6/125$; then $e_1=37/6$, $e_2=29/3$ and $e_3=4$ give the values $37/125$, $58/125$ and $24/125$ of Section 1. Newton's inequalities for nonnegative numbers state ^[3]

$$\Bigl(\frac{e_k}{\binom nk}\Bigr)^2\ge\frac{e_{k-1}}{\binom n{k-1}}\cdot\frac{e_{k+1}}{\binom n{k+1}},$$

and since every $e_j$ is positive and $\binom nk^2\gt\binom n{k-1}\binom n{k+1}$, it follows for $1\le k\le n-1$ that

$$\Pr[T=k]^2\gt\Pr[T=k-1]\,\Pr[T=k+1].$$

So the ratio $\Pr[T=k+1]/\Pr[T=k]$ strictly decreases as $k$ grows: the probabilities rise strictly while this ratio exceeds $1$ and fall strictly once it is below $1$. In particular they rise strictly up to the smallest most likely value $m$: if the ratio were at most $1$ at some $k\lt m$, every later ratio would be below $1$, and $\Pr[T=m]$ could not exceed $\Pr[T=k]$. Where $m$ lies is settled by Darroch's rule: the most likely value of $T$ is $\lfloor\mu\rfloor$ or $\lceil\mu\rceil$, where $\mu$ is the mean of $T$. ^[4] For the three shots, $\mu=9/5$ and the most likely score is $2=\lceil\mu\rceil$.

Applying this to $S_x$. The mean of $S_x$ is $49-(1275-x)/q$, which for $q\gt50$ is larger than $49-1274/50=23.52$. So the most likely value of $S_x$ is at least $23$, and $\Pr[S_x=19]\lt\Pr[S_x=20]$ for every $x$. Every $d_x$ is negative, so $F'(q)\lt0$, which proves the claim.

6. Bisection

Because $F$ is continuous and strictly decreasing, the bisection method finds the root reliably. ^[5] It keeps an interval $[L,H]$ with $F(L)\gt0.02\ge F(H)$, so that the root lies inside, and tests the midpoint $M=(L+H)/2$: if $F(M)\gt0.02$, the root lies to the right of $M$ and $M$ becomes the new $L$; otherwise $M$ becomes the new $H$. Each step halves the interval at the cost of one evaluation of $F$. Starting from $[50,100]$, the first steps are:

$M=75$, $62.5$, $56.25$, $53.125$: each time $F(M)\lt0.02$, so $H$ moves down to $53.125$;
$M=51.5625$: now $F(M)\gt0.02$, so $L=51.5625$.

After $k$ steps the interval has width $50/2^k$, and $43$ steps already bring it below $10^{-11}$. Newton's method can converge faster close to the root, but it needs the derivative and a good start.

7. Is double precision accurate enough?

Write $u=2^{-53}\approx1.1\cdot10^{-16}$ for the largest relative error of one rounding in double precision. ^[6] We bound the error of the computed $F$ for every $q$ between $52$ and $53$, and then its effect on the root.

The weights. The only subtractions are $1-x/q$ and $1$ minus the computed hit probability. Since $x/q\lt1$, the computed hit and miss probabilities are off by less than $2u$ and $3u$. By Section 5, changing the two weights of shot $x$ changes $F$ by these changes times probabilities of the other $49$ shots, which are at most $1$, so the $50$ shots contribute less than $250u$.
The updates. Each update adds two nonnegative products, so no cancellation occurs, and a table entry collects at most two roundings per shot. After $50$ shots its relative error is at most about $100u$, which is less than $3u$ for a probability below $0.03$.
The root. So the computed $F$ is off by less than $3\cdot10^{-14}$. Since $|F'(q)|\ge0.005$ between $52$ and $53$, there $F(M)$ differs from $0.02$ by at least $0.005$ times the distance from $M$ to the root, and outside this interval by more than $0.0019$. A comparison $F(M)\gt0.02$ can therefore come out wrong only for an $M$ within $3\cdot10^{-14}/0.005=6\cdot10^{-12}$ of the root, and the final endpoints, neighbouring doubles $2^{-47}\approx7.1\cdot10^{-15}$ apart, lie within about $6\cdot10^{-12}$ of it.

The exact root lies more than $10^{-11}$ away from the nearest point at which the rounded tenth decimal would change, so the printed digits are correct. The actual errors are far smaller than these bounds: at the root the computed $F$ is off by about $2.4\cdot10^{-18}$, and the result of the bisection differs from the root computed with $80$ digits by $3.5\cdot10^{-15}$.

How the Code Works

The C++ program ^[7] keeps the parameters of the problem and a flag for its built-in checks in Options. parse_arguments changes them with --shots=, --target-score=, --target-probability= and --skip-checkpoints. It requires at least one shot, a target score from $0$ to the number of shots and a target probability strictly between $0$ and $1$; otherwise, or on an unknown argument, the program ends with exit code 1.

exact_score_probability implements Section 3 with a single table dp of type std::array<long double, 64>. It sets dp[0] to $1$ and, for every shot $x$, computes hit as $1-x/q$ and miss as $1$ minus hit. The inner loop runs $s=x,x-1,\dots,1$ and replaces dp[s] by dp[s]*miss + dp[s-1]*hit; afterwards dp[0] is multiplied by miss. The downward order lets one table replace the rows of Section 3: when dp[s] is overwritten, dp[s-1] still holds its value for $x-1$ shots. All scores up to $x$ are kept, although only dp[target_score] is returned.

solve implements Section 6. Starting from low at $50$ and high at $100$, it doubles high while the probability there is still above the target, which never happens for the default problem, then halves the interval $220$ times and returns the midpoint. run_checkpoints checks $F(52)\gt0.02\gt F(53)$, solves the default problem and checks that the probability at the solution differs from $0.02$ by at most $10^{-13}$; if a check fails, the program ends with exit code 2. Finally main solves the problem with the chosen options, converts the result to double and prints it with $10$ decimals.

These checks do not cover everything the method needs, and Section 5 is proved for the default parameters only. The table has $64$ entries, but any number of shots is accepted: --shots=64 accesses entries past its end, as AddressSanitizer reports, and the digit loop of parse_int_after_prefix can overflow an int. The lower end $50$ of the search is fixed, so --shots=60 returns a $q$ below $60$, where some $1-x/q$ are negative, and nothing checks that a solution exists: --target-probability=0.9 prints 50.0000000000.

The type long double is the 80-bit extended format with GCC on x86 processors, but the same as double with Microsoft Visual C++ and on 64-bit ARM under Windows and macOS. ^[8] In the optimised Clang and GCC builds on 64-bit ARM used here, a multiplication and the following addition are also fused into one rounding, which changes the last bits of some probabilities. Both effects only shrink the bounds of Section 7, so the printed result is the same. Python and Java have the same probability function (in Java exactScoreProbability) and bisection in double precision, without checks; Python's solve takes the parameters as optional arguments, and Java formats with Locale.US to print a decimal point in every locale.

Complexity and Verification

One evaluation of $F$ for $n$ shots makes $n(n+1)/2$ two-term updates, $1275$ for $n=50$; stopping each inner loop at the target score would cut them to $810$. A bisection with $I$ steps and no doubling makes $I+1$ evaluations, so the time is $O(In^2)$; the method needs $O(n)$ memory, and the programs use a fixed table of $64$ numbers. With its checks, the C++ program evaluates $F$ $445$ times, $567\,375$ updates in all, and the Python and Java programs $221$ times; on the ARM computer used here each finishes within a few hundredths of a second. In these runs only the first $52$ of the $220$ halvings change the interval; in an emulation of the 80-bit format, the first $64$.

Besides the self-checks of run_checkpoints, the following checks were made for this article; they are not part of the published solutions.

All numbers of the three-shot example were recomputed with exact fractions.
For $716$ values of $q$ ($50.000001$, every $0.05$ up to $69.95$, every $3$ from $70$ to $1000$, and $10^4,\dots,10^8$), the formula for $F'(q)$ of Section 5, with $40$ digits, agrees with a central difference and is negative; every $S_x$ has its most likely value at $24$ or higher, and the ratio $\Pr[S_x=20]/\Pr[S_x=19]$ is never below $1.61$.
On a grid of step $0.01$ from $52$ to $53$, $|F'(q)|$ falls steadily from $0.0067$ to $0.0052$, which gives the bound of Section 7. An $80$-digit bisection gives the root that Section 7 compares with the double-precision result of the Python program; the same bisection with a $64$-bit significand, emulated with the mpmath library, differs from it by $3.6\cdot10^{-18}$.
The C++ program, built with Clang and with GCC on a 64-bit ARM computer, and the Python and Java programs print the same ten decimals.

Footnotes and References

The references below supply the background facts used above. The polynomial form of $F$, the proof that $F$ decreases and the accuracy estimates are derived explicitly in this article.

Project Euler 286 — Scoring Probabilities. The official statement: from distance x Barbara scores with probability 1 − x/q for a real constant q greater than 50; shooting from x = 1, 2, …, 50 she scores exactly 20 points with probability precisely 2%, and the task asks for q rounded to 10 decimal places.
Wikipedia — Poisson binomial distribution. The distribution of a sum of independent Bernoulli trials that need not be identically distributed: its probability generating function is the product of the factors 1 − pⱼ + pⱼz, its mean is the sum of the pⱼ, and the direct convolution algorithm, which adds one trial at a time, is numerically stable.
Wikipedia — Newton's inequalities. For non-negative real numbers a₁, …, aₙ with elementary symmetric polynomials eₖ and elementary symmetric means Sₖ = eₖ / C(n, k), the inequality Sₖ₋₁ · Sₖ₊₁ ≤ Sₖ² holds.
Alexander Gnedin — Cross Modality of the Extended Binomial Sums (arXiv:2408.06477). Section 5.1 recalls Darroch's rule (J. N. Darroch, Ann. Math. Statist. 35 (1964), 1317–1321): the mode of a sum of independent Bernoulli variables is the mean when the mean is an integer, and otherwise the integer just below or just above the mean, or both.
Wikipedia — Bisection method. If a continuous function has opposite signs at the ends of an interval, the intermediate value theorem gives a root inside; each step keeps the half in which the sign changes, halving the width.
Wikipedia — Double-precision floating-point format. A double has a 53-bit significand; in the range from 2ⁿ to 2ⁿ⁺¹ neighbouring doubles are 2ⁿ⁻⁵² apart, and the maximum relative error of rounding to the nearest double is 2⁻⁵³.
C++ — Euler286.cpp. The linked immutable C++ revision contains exact_score_probability, solve and the checks of run_checkpoints. It is the implementation record for this article.
Wikipedia — Long double. With GCC, long double is the 80-bit extended precision format on x86 processors; it is a synonym for double with Microsoft Visual C++ and on 64-bit ARM under Windows and macOS.

Problem 286 source code

C++

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

namespace {

struct Options {
    int shots = 50;
    int target_score = 20;
    long double target_probability = 0.02L;
    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, "--shots=", options.shots) ||
            parse_int_after_prefix(arg, "--target-score=", options.target_score) ||
            parse_ld_after_prefix(arg, "--target-probability=", options.target_probability)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.shots >= 1 && options.target_score >= 0 &&
           options.target_score <= options.shots &&
           options.target_probability > 0.0L && options.target_probability < 1.0L;
}

long double exact_score_probability(const int shots, const int target_score, const long double q) {
    std::array<long double, 64> dp{};
    dp.fill(0.0L);
    dp[0] = 1.0L;
    for (int x = 1; x <= shots; ++x) {
        const long double hit = 1.0L - static_cast<long double>(x) / q;
        const long double miss = 1.0L - hit;
        for (int s = x; s >= 1; --s) {
            dp[static_cast<std::size_t>(s)] =
                dp[static_cast<std::size_t>(s)] * miss +
                dp[static_cast<std::size_t>(s - 1)] * hit;
        }
        dp[0] *= miss;
    }
    return dp[static_cast<std::size_t>(target_score)];
}

long double solve(const Options& options) {
    long double low = 50.0L;
    long double high = 100.0L;
    while (exact_score_probability(options.shots, options.target_score, high) > options.target_probability) {
        high *= 2.0L;
    }

    for (int it = 0; it < 220; ++it) {
        const long double mid = (low + high) / 2.0L;
        if (exact_score_probability(options.shots, options.target_score, mid) > options.target_probability) {
            low = mid;
        } else {
            high = mid;
        }
    }
    return (low + high) / 2.0L;
}

bool run_checkpoints() {
    const long double p52 = exact_score_probability(50, 20, 52.0L);
    const long double p53 = exact_score_probability(50, 20, 53.0L);
    if (!(p52 > 0.02L && p53 < 0.02L)) {
        std::cerr << "Checkpoint failed for bracketing around q=52..53" << '\n';
        return false;
    }
    const Options options;
    const long double q = solve(options);
    const long double prob = exact_score_probability(options.shots, options.target_score, q);
    if (std::fabsl(prob - options.target_probability) > 1e-13L) {
        std::cerr << "Checkpoint failed for solved q probability consistency" << '\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 answer = solve(options);
    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

def exact_score_probability(shots, target_score, q):
    dp = [0.0] * 64
    dp[0] = 1.0
    for x in range(1, shots + 1):
        hit = 1.0 - float(x) / q
        miss = 1.0 - hit
        for s in range(x, 0, -1):
            dp[s] = dp[s] * miss + dp[s - 1] * hit
        dp[0] *= miss
    return dp[target_score]

def solve(shots=50, target_score=20, target_probability=0.02):
    low = 50.0
    high = 100.0
    while exact_score_probability(shots, target_score, high) > target_probability:
        high *= 2.0
        
    for _ in range(220):
        mid = (low + high) / 2.0
        if exact_score_probability(shots, target_score, mid) > target_probability:
            low = mid
        else:
            high = mid
            
    return f"{(low + high) / 2.0:.10f}"

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

Java

import java.util.Locale;

public class Euler286 {
    static double exactScoreProbability(int shots, int targetScore, double q) {
        double[] dp = new double[64];
        dp[0] = 1.0;
        for (int x = 1; x <= shots; ++x) {
            double hit = 1.0 - (double) x / q;
            double miss = 1.0 - hit;
            for (int s = x; s >= 1; --s) {
                dp[s] = dp[s] * miss + dp[s - 1] * hit;
            }
            dp[0] *= miss;
        }
        return dp[targetScore];
    }

    public static String solve() {
        int shots = 50;
        int targetScore = 20;
        double targetProbability = 0.02;

        double low = 50.0;
        double high = 100.0;
        while (exactScoreProbability(shots, targetScore, high) > targetProbability) {
            high *= 2.0;
        }

        for (int it = 0; it < 220; ++it) {
            double mid = (low + high) / 2.0;
            if (exactScoreProbability(shots, targetScore, mid) > targetProbability) {
                low = mid;
            } else {
                high = mid;
            }
        }

        return String.format(Locale.US, "%.10f", (low + high) / 2.0);
    }

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