Problem 121: Disc Game Prize Fund

View on Project Euler

Project Euler Problem 121 Solution

EulerSolve provides an optimized solution for Project Euler Problem 121, Disc Game Prize Fund, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The game starts with one blue disc and one red disc in the bag. On turn \(t\), the player draws one disc uniformly at random, notes its color, returns it, and then an extra red disc is added. Therefore the bag contains exactly one blue disc and \(t\) red discs on turn \(t\), so $$P(\text{blue on turn }t)=\frac{1}{t+1},\qquad P(\text{red on turn }t)=\frac{t}{t+1}.$$ There are \(15\) turns in total, and the player wins exactly when blue appears more often than red. Since \(15\) is odd, that means at least \(8\) blue draws. If the total prize fund is \(M\) pounds, including the original stake, then the banker does not lose money on average precisely when \(M \cdot P(\text{win}) \le 1\). The task is to compute the largest integer \(M\) satisfying that inequality. Mathematical Approach The crucial quantity is not the order of the draws, but how many blue draws have occurred after each turn. That turns the problem into a small distribution update over the possible blue counts. The right state variable Let \(B_t\) be the number of blue draws after the first \(t\) turns, and let $$P_t(b)=P(B_t=b).$$ Only the values \(0 \le b \le t\) are possible....

Detailed mathematical approach

Problem Summary

A bag starts with one red disc and one blue disc. On each turn the player takes a disc at random and its colour is noted; the disc is returned and an extra red disc is added. The player pays £1 and wins if more blue than red discs have been taken at the end. A winner receives the prize fund, a whole number of pounds that includes the £1 stake. For four turns the winning probability is \(11/120\), and £10 is the largest prize fund for which the banker does not expect a loss. We must find this maximum prize fund for a game of fifteen turns. [1]

The whole solution in one sentence: before turn \(t\) the bag holds one blue and \(t\) red discs, so the turns are independent with blue probability \(1/(t+1)\); a table that grows by one turn at a time counts the winning games exactly, and the prize fund is the whole part of one divided by the winning probability.

Mathematical Approach

Throughout, \(n\) is the number of turns and \(B_n\) the number of blue discs taken in them.

1. The bag before each turn

Every turn returns the drawn disc and adds a red one, whatever colour was drawn. So before turn \(t\) the bag holds one blue disc and \(t\) red discs, and the blue disc is taken with probability

$$p_t=\frac{1}{t+1},\qquad 1-p_t=\frac{t}{t+1}.$$

The bag never depends on the earlier draws, so the turns are independent, and the probability of a colour sequence is the product of the probabilities of its turns: for three turns, blue, blue, red has probability \(\frac12\cdot\frac13\cdot\frac34=\frac18\). A count of successes in independent trials with different probabilities, such as \(B_n\), has a Poisson binomial distribution. [2]

2. Who wins

With \(b\) blue discs there are \(n-b\) red ones, and the player wins when \(b\gt n-b\), that is

$$b\ge\lfloor n/2\rfloor+1,$$

at least \(3\) blue discs in four turns and \(8\) in fifteen; a tie, possible when \(n\) is even, does not win. By linearity of expectation [5] a game yields on average \(p_1+p_2+\dots+p_n\) blue discs, \(77/60\approx1.28\) for four turns and only about \(2.38\) for fifteen, so a win is rare.

3. Counting equally likely outcomes

To avoid fractions we tell the discs apart. On turn \(t\) each of the \(t+1\) discs is equally likely, so a game is a list of \(n\) taken discs, and all \(2\cdot3\cdots(n+1)=(n+1)!\) lists are equally likely. A colour sequence with red discs on the turns of a set \(R\) comes from \(\prod_{t\in R}t\) lists, because on a red turn \(t\) any of the \(t\) red discs may be taken and on a blue turn there is one choice. Hence the winning probability is

$$P_n=\frac{W_n}{(n+1)!},$$

where \(W_n\) is the number of winning lists. Four turns by hand: the player needs \(3\) or \(4\) blue discs. Four blue discs come from one list, and a single red disc on turn \(t\) from \(t\) lists. So \(W_4=1+(1+2+3+4)=11\) of the \(2\cdot3\cdot4\cdot5=120\) lists win, the \(11/120\) of the statement. Likewise three turns need \(2\) blue discs, and \(W_3=1+(1+2+3)=7\) of the \(24\) lists win.

4. A table that grows by one turn

Let \(A_t(b)\) be the number of lists of the first \(t\) turns with exactly \(b\) blue discs. The disc of turn \(t\) is red in \(t\) ways, and then the earlier turns had \(b\) blue discs, or it is the blue disc, and then they had \(b-1\). Therefore

$$A_t(b)=t\,A_{t-1}(b)+A_{t-1}(b-1),$$

starting from \(A_0(0)=1\), the one empty list, with entries of row \(t-1\) outside \(0\le b\le t-1\) counted as \(0\). For example, \(A_2(1)=2\cdot A_1(1)+A_1(0)=3\): one blue disc in two turns is blue and then one of the two red discs, or red and then blue. Filling a table from smaller results, each computed once, is called dynamic programming. [6] The first rows, with \(b=0,1,\dots,t\) from left to right, are

$$\begin{aligned} t=0&:\ 1\\ t=1&:\ 1,\ 1\\ t=2&:\ 2,\ 3,\ 1\\ t=3&:\ 6,\ 11,\ 6,\ 1\\ t=4&:\ 24,\ 50,\ 35,\ 10,\ 1 \end{aligned}$$

Row \(t\) sums to \((t+1)!\), since every list has some number of blue discs, and the winning entries give \(W_3=6+1=7\) and \(W_4=10+1=11\), as above. In general

$$W_n=\sum_{b=\lfloor n/2\rfloor+1}^{n}A_n(b).$$

Dividing by \((t+1)!\) gives the same rule for the probabilities \(P(B_t=b)=A_t(b)/(t+1)!\), the form used by the Python program:

$$\begin{aligned} P(B_t=b)&=(1-p_t)\,P(B_{t-1}=b)\\ &\quad+p_t\,P(B_{t-1}=b-1). \end{aligned}$$

5. A product and the Stirling numbers

For the polynomial \(F_t(x)=\sum_b A_t(b)\,x^b\) the recurrence reads \(F_t(x)=(t+x)\,F_{t-1}(x)\): the term \(t\) is a red disc and keeps the power of \(x\), the term \(x\) is the blue disc and raises it. Hence

$$F_n(x)=(1+x)(2+x)\cdots(n+x),$$

for instance

$$\begin{aligned} &(1+x)(2+x)(3+x)(4+x)\\ &\quad=24+50x+35x^2+10x^3+x^4. \end{aligned}$$

Since \((n+1)!=2\cdot3\cdots(n+1)\), dividing \(F_n(x)\) by \((n+1)!\) divides each factor \(t+x\) by \(t+1\) and turns it into \((1-p_t)+p_t\,x\), the probability generating function of one turn. The product, whose coefficient of \(x^b\) is \(P(B_n=b)\), is the generating function of \(B_n\), since the generating function of a sum of independent variables is the product of theirs. [3]

A permutation of \(m\) objects splits them into cycles, found by following an object to its image, and so on, until it returns; swapping \(1\) and \(3\) and fixing \(2\) gives the two cycles \((1\,3)(2)\). The unsigned Stirling numbers of the first kind \(c(m,k)\) count the permutations of \(m\) objects with \(k\) cycles and are the coefficients of \(x^k\) in \(x(x+1)\cdots(x+m-1)\), [4] so \(A_n(b)=c(n+1,b+1)\); for \(m=5\) and \(k=1,\dots,5\) they are \(24, 50, 35, 10, 1\). Their recurrence \(c(t+1,k)=t\,c(t,k)+c(t,k-1)\) matches the game: a new object either forms a cycle of its own, in one way like the blue disc, or follows one of the \(t\) earlier objects in its cycle, in \(t\) ways like the red discs. So fifteen turns are won exactly as often as a random permutation of \(16\) objects has at least \(9\) cycles.

6. From the probability to the prize fund

The banker receives £1 and pays the prize fund \(M\), stake included, with probability \(P_n\), so the expected gain per game is \(1-M\,P_n\). [5] There is no expected loss exactly when \(M\le1/P_n\), and as \(M\) is a whole number,

$$M_n=\left\lfloor\frac{(n+1)!}{W_n}\right\rfloor.$$

For four turns \(120/11\approx10.91\) gives \(M_4=10\): the expected payout is then \(10\cdot\frac{11}{120}=\frac{110}{120}\) pounds, less than the stake, while \(M=11\) would give \(\frac{121}{120}\). The boundary case \(M\,P_n=1\) is not a loss, and it occurs: \(P_1=1/2\) and \(P_2=1/6\) give \(M_1=2\) and \(M_2=6\). Exact arithmetic keeps such a quotient from being rounded to just below a whole number.

7. Small values

For \(n=1,\dots,8\) the winning probabilities are \(\frac12\), \(\frac16\), \(\frac{7}{24}\), \(\frac{11}{120}\), \(\frac{101}{720}\), \(\frac{197}{5040}\), \(\frac{2311}{40320}\) and \(\frac{5119}{362880}\), and the prize funds are \(2\), \(6\), \(3\), \(10\), \(7\), \(25\), \(17\) and \(70\). They zigzag: \(2k+1\) turns need \(k+1\) blue discs, as many as \(2k\) turns, but give one more draw to get them, while \(2k+2\) turns need one blue disc more.

How the Code Works

The C++ program [7] starts in main. parse_arguments accepts --turns=N, by default \(15\), digits only and at least \(1\), and --skip-checkpoints; unless the latter is given, run_checkpoints requires prize_fund(4) to be \(10\), and the program exits with code \(2\) otherwise. The result of prize_fund is printed as a long long, since the standard streams cannot print 128-bit integers.

prize_fund keeps the table of Section 4 in a std::vector of signed 128-bit integers, i128, an alias of __int128_t, indexed by the number of blue discs. Starting from dp[0] = 1, each turn \(t\) fills a fresh vector next, whose entry blue collects dp[blue] * t, a red disc, when blue is at most \(t-1\), and dp[blue - 1], the blue disc, when blue is positive. Then denominator becomes \(2\cdot3\cdots(n+1)\), the entries from turns / 2 + 1 to turns add up to win_numerator \(=W_n\), and denominator / win_numerator is \(M_n\), because integer division discards the fractional part. [8] For fifteen turns no number exceeds \(16!\approx2.1\cdot10^{13}\), a 45-bit number. No upper limit for \(N\) is enforced, however: from \(N=33\) on, \((N+1)!\ge34!\gt2^{127}\) overflows the signed 128-bit type, which is undefined behaviour, so no output is guaranteed; the digit parser does not guard its int against overflow either.

The Python file calls solve(), which starts from dp[0] = Fraction(1) and stores the probabilities themselves as exact Fraction values in lowest terms. [8] Every nonzero state b sends dp[b] * p_red to new_dp[b] and dp[b] * p_blue to new_dp[b + 1], where p_red is \(t/(t+1)\) and p_blue is \(1/(t+1)\), and the program prints int(1 / win_prob), the whole part of \(1/P_n\). Java's main pushes whole BigInteger numbers forward in the same way from dp[0] = BigInteger.ONE, forms denom and winNum like the C++ program and prints denom.divide(winNum), which truncates like Java's integer division. Both fix the number of turns at \(15\).

Complexity and Verification

Before turn \(t\) only the \(t\) states \(b\le t-1\) are reachable, each with two contributions, so \(n\) turns make \(n(n+1)\) contributions, \(240\) for fifteen turns: \(O(n^2)\) arithmetic operations in all, with two rows of \(n+1\) numbers stored. All numbers, including the numerators and denominators of the fractions, are at most \((n+1)!\) and have \(O(n\log n)\) bits. Listing all \(2^n\) colour sequences would also be quick for \(n=15\), but it costs about \(n\,2^n\) steps.

The C++ program checks prize_fund(4) \(=10\) unless --skip-checkpoints is given; the Python and Java programs have no checks. Separately from the published code, we verified that

  • the three programs print the same value for fifteen turns, equal to \(M_{15}\) from the brute-force \(W_{15}\) of the third item;
  • enumerating all lists of taken discs for \(n\le8\), up to \(9!=362880\) lists, gives the table's \(W_n\), including \(7\) of \(24\) and \(11\) of \(120\);
  • enumerating all \(2^n\) colour sequences with their weights \(\prod_{t\in R}t\) gives the table's \(W_n\) for every \(n\le20\);
  • counting the permutations of up to \(8\) objects by their cycles reproduces the rows \(A_n\) for \(n\le7\);
  • the C++ option --turns=N is exact for every \(N\le32\); for \(N\) from \(33\) to \(40\) our build printed wrong values, as expected after the overflow.

Footnotes and References

The references below supply the background facts used above. The counting of lists, the recurrence for \(A_t(b)\), the prize-fund formula and the small values are derived explicitly in this article.

  1. Project Euler 121 — Disc Game Prize Fund. The official statement: one red and one blue disc at the start, the drawn disc returned and a red disc added after each turn, a £1 stake and a win with more blue than red discs; for four turns it gives 11/120 and £10, stake included, and it asks for the maximum prize fund for fifteen turns.
  2. Wikipedia — Poisson binomial distribution. The distribution of a sum of independent Bernoulli trials that need not be identically distributed. Its generating function is the product of the factors 1 − pⱼ + pⱼz, and its direct convolution algorithm adds one trial at a time, like the table of Section 4.
  3. Wikipedia — Probability-generating function. For a random variable X with values 0, 1, 2, … it is G(z) = E[zˣ]; the generating function of a sum of independent random variables is the product of theirs.
  4. Wikipedia — Stirling numbers of the first kind. The unsigned numbers count permutations by cycles, are the coefficients of the rising factorial x(x + 1)⋯(x + n − 1) and satisfy [n + 1, k] = n·[n, k] + [n, k − 1], where the new object forms a cycle of its own or is inserted in one of n ways; for n = 5 and k = 1, …, 5 they are 24, 50, 35, 10, 1.
  5. Wikipedia — Expected value. With finitely many outcomes the expectation is the sum of the outcomes weighted by their probabilities, and it is linear: E[X + Y] = E[X] + E[Y].
  6. Wikipedia — Dynamic programming. A problem is broken into simpler sub-problems that are each solved once and stored; the bottom-up form solves the smaller ones first, as the rows with fewer turns are filled first here.
  7. C++ — Euler121.cpp. The linked immutable revision contains parse_arguments, prize_fund and run_checkpoints. It is the implementation record for this article.
  8. Exact arithmetic in the three languages: by the C++ working draft, [expr.mul], / on integers discards any fractional part of the quotient, and by [expr.pre] an arithmetic result outside the range of its type is undefined behaviour; Python's fractions module keeps every Fraction in lowest terms; Java's BigInteger holds immutable arbitrary-precision integers whose operations mimic Java's integer operators.

Problem 121 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i128 = __int128_t;

struct Options {
    int turns = 15;
    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_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, "--turns=", options.turns)) {
            continue;
        }

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

    return options.turns >= 1;
}

i128 prize_fund(const int turns) {
    std::vector<i128> dp(static_cast<std::size_t>(turns + 1), 0);
    dp[0] = 1;

    for (int t = 1; t <= turns; ++t) {
        std::vector<i128> next(static_cast<std::size_t>(turns + 1), 0);
        for (int blue = 0; blue <= t; ++blue) {
            if (blue <= t - 1) {
                next[static_cast<std::size_t>(blue)] += dp[static_cast<std::size_t>(blue)] * t;
            }
            if (blue > 0) {
                next[static_cast<std::size_t>(blue)] += dp[static_cast<std::size_t>(blue - 1)];
            }
        }
        dp = std::move(next);
    }

    i128 denominator = 1;
    for (int t = 1; t <= turns; ++t) {
        denominator *= (t + 1);
    }

    i128 win_numerator = 0;
    for (int blue = turns / 2 + 1; blue <= turns; ++blue) {
        win_numerator += dp[static_cast<std::size_t>(blue)];
    }

    return denominator / win_numerator;
}

bool run_checkpoints() {
    if (prize_fund(4) != 10) {
        std::cerr << "Checkpoint failed for turns=4" << '\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;
    }

    std::cout << static_cast<long long>(prize_fund(options.turns)) << '\n';
    return 0;
}

Python

# Problem 121: Disc game prize fund
# DP on probability of blue wins over 15 turns.

from fractions import Fraction

def solve():
    turns = 15
    # dp[b] = probability of getting exactly b blue discs
    dp = [Fraction(0)] * (turns + 1)
    dp[0] = Fraction(1)
    for t in range(1, turns + 1):
        new_dp = [Fraction(0)] * (turns + 1)
        p_blue = Fraction(1, t + 1)
        p_red = Fraction(t, t + 1)
        for b in range(t + 1):
            if dp[b] == 0: continue
            new_dp[b] += dp[b] * p_red
            if b + 1 <= turns:
                new_dp[b + 1] += dp[b] * p_blue
        dp = new_dp
    win_prob = sum(dp[b] for b in range(turns // 2 + 1, turns + 1))
    print(int(1 / win_prob))

solve()

Java

import java.math.BigInteger;

public class Euler121 {
    public static void main(String[] args) {
        int turns = 15;
        BigInteger[] dp = new BigInteger[turns + 1];
        for (int i = 0; i <= turns; i++)
            dp[i] = BigInteger.ZERO;
        dp[0] = BigInteger.ONE;
        for (int t = 1; t <= turns; t++) {
            BigInteger[] next = new BigInteger[turns + 1];
            for (int i = 0; i <= turns; i++)
                next[i] = BigInteger.ZERO;
            for (int b = 0; b <= t; b++) {
                if (dp[b].signum() == 0)
                    continue;
                // red: multiply by t (red discs)
                next[b] = next[b].add(dp[b].multiply(BigInteger.valueOf(t)));
                // blue: multiply by 1
                if (b + 1 <= turns)
                    next[b + 1] = next[b + 1].add(dp[b]);
            }
            dp = next;
        }
        BigInteger denom = BigInteger.ONE;
        for (int t = 1; t <= turns; t++)
            denom = denom.multiply(BigInteger.valueOf(t + 1));
        BigInteger winNum = BigInteger.ZERO;
        for (int b = turns / 2 + 1; b <= turns; b++)
            winNum = winNum.add(dp[b]);
        System.out.println(denom.divide(winNum));
    }
}