To be honest, problem 121 of Project Euler scared me a whole lot, since my probability theory is not very strong. However, I do believe that I came up with quite a nice solution after all. The problem reads

A bag contains one red disc and one blue disc. In a game of chance a player takes a disc at random and its colour is noted. After each turn the disc is returned to the bag, an extra red disc is added, and another disc is taken at random.

The player pays £1 to play and wins if they have taken more blue discs than red discs at the end of the game.

If the game is played for four turns, the probability of a player winning is exactly 11/120, and so the maximum prize fund the banker should allocate for winning in this game would be £10 before they would expect to incur a loss. Note that any payout will be a whole number of pounds and also includes the original £1 paid to play the game, so in the example given the player actually wins £9.

Find the maximum prize fund that should be allocated to a single game in which fifteen turns are played.

At first I tried to solve it by making Monte Carlo simulations, but I kept getting different results. Which shows me that I would need to go very far up in the number of simulations. And that in turn took a long time. So I had to find another approach.

Since the events of drawing a disc are independent, I tried to draw it as scenario trees, which seems to be a valid approach. So this is the solution I will present to you.

## Solving for the small example

Let us start out by solving for the case of four turns. First of all there is a total of 5! = 120 outcomes. Since for one drawing there are 2 = 2! outcomes, for two drawings there are 2*3 = 6 = 3! outcomes and so on. So for n drawings there are (n+1)! outcomes.

If we draw once where we have 1 red and 1 blue disc, we get a scenario tree looking like

Which means in 1 of the 2! cases there is drawn 1 blue, and in the other case there is drawn 1 red. Pretty simple so far. If we add another drawing the tree becomes

Starting from the left, we see that there if we have already drawn a blue, there is one scenario more where we will draw another blue giving us a total of 1 scenario where we will draw 2 blues. Then it becomes interesting. We can get to the center of this line from two places. If we have drawn a blue the first time we can draw a red this time. Since there are 2 red discs there are 2 scenarios out of three where we will draw a red. This gives a total of 2 scenarios that will end us in the center if we draw a blue in the first drawing. However, we can also draw a red disc in first drawing. And from there we have 1 scenario where we draw a blue. And the last I will let you figure out. This means for a 2 drawings we have 1/6 chance of winning.

Let us add another layer

The leftmost is trivial to explain, so let us jump to the case where we end up with 2 reds and 1 blue. In second drawing there are 3 scenarios in which we end up with 1 blue and 1 red disc. From there we have 3 scenarios which will give us yet another red. That gives us 3*3=9 scenarios that will land us in having 1 blue and 2 red. Add to that another 2 scenarios from the fact that we can have 2 reds and then draw a blue.

Soooo now we only miss the last layer

My guess is that you can follow the logic now. So how do we see if we have won? The leftmost is the number of scenarios where we have drawn only blues. The second is the scenarios where we have 3 blue and 1 red. The center one we have 2 blue and 2 red. So it is only the 2 leftmost cases where we win. Which gives us a total of 11 scenarios in which we win.

Based on the scenario tree we can make an implementation which will calculate the number of scenarios. I have made the c# implementation as this

int limit = 15; long[] outcomes = new long[limit+1]; outcomes[limit] = 1; outcomes[limit-1] = 1; for (int i = 2; i <= limit; i++) { for(int j = 0; j < outcomes.Length-1; j++){ outcomes[j] = outcomes[j + 1]; } outcomes[limit] = 0; for (int j = outcomes.Length-1; j > 0; j--) { outcomes[j] += outcomes[j - 1] * i; } } long positive = 0; for (int i = 0; i < limit / 2+1; i++) { positive += outcomes[i]; }

All we need to do now is to calculate the total number of scenarios which is 16! in our case. We can make this calculation in 0.0021ms on my computer. So even though there might be prettier implementations out there, I don’t think you can find one which is much faster.

The result is

There are 9219406943 positive outcomes out of 20922789888000 This gives a prize allocation of 2269

## Wrapping up

The outcomes of the above scenarios are the same as the coefficients for the polynomial (x+1)(x+2)(x+3)(x+4). I am not sure what you can use the information for, but who knows.

You can find the source code here.

These kind of problems quite often come up in programming contests (for an example, see problem D: Darts from SWERC 2009). It wasn’t until quite recently I started feeling a little comfortable solving them. The process usually consists of formulating the problem recursively and then computing the formula using dynamic programming.

After becoming familiar with these problems, this problem became straightforward for me. I created a function [latex]\mathrm{f}(n,i,r,b)[/latex], where [latex]n[/latex] is the total number of rounds, [latex]i[/latex] is the round we are currently on (0 based, so in the first round [latex]i=0[/latex]), [latex]r[/latex] is the number of red discs we’ve drawn, and [latex]b[/latex] is the number of blue discs we’ve drawn. [latex]\mathrm{f}(n,i,r,b)[/latex] then returns the probability of winning the game if we are on round [latex]i[/latex] of [latex]n[/latex] rounds, and have already gotten the specified amount of red and blue discs.

So first, the base case is when [latex]n = i[/latex]. Then we’ve played all the rounds and the game is over. Then we simply check whether [latex]b[/latex] is bigger than [latex]r[/latex]. If [latex]b > r[/latex], we return [latex]1[/latex], since we’ve won the game and therefore the probability of winning is [latex]1[/latex]. If [latex]b \leq r[/latex] then we’ve lost the game and so we return [latex]0[/latex].

When [latex]i < n[/latex] we still have some rounds to play. On the [latex]i[/latex]th round the number of blue discs is of course [latex]1[/latex], the number of red discs is [latex]i + 1[/latex], and the total number of discs is [latex]i + 2[/latex]. Now there are two cases. We can either draw a blue disc with probability [latex]\frac{1}{i+2}[/latex], or draw a red disc with probability [latex]\frac{i+1}{i+2}[/latex]. If we draw a blue disc, then the probability of winning is [latex]\mathrm{f}(n,i+1,r,b+1)[/latex]. If we draw a red disc, then the probability of winning is [latex]\mathrm{f}(n,i+1,r+1,b)[/latex]. Using the multiplication and addition rules, we can conclude that

The probability of winning is then [latex]\mathrm{f}(15,0,0,0)[/latex], but the answer to the problem is [latex]\lfloor{}\frac{1}{\mathrm{f}(15,0,0,0)}\rfloor{}[/latex]. This would then normally be computed using dynamic programming (and there are overlapping subproblems for sure), but Mathematica was able to compute the final value in a couple of milliseconds without any optimization, so I didn’t bother.

I hope I explained the idea clearly, especially because I think there are more Project Euler problems that can be solved with very similar techniques (some of them are also about expected values).

Anyways, I’m glad I finally solved it. I had previously tried solving it, but could never figure it out. Also, now I finally get to mark that MathBlog subscription email about “Project Euler 121” as read, since I’ve kept it unread in my inbox until I could solve it myself (didn’t want to spoil the problem by reading the post here). 🙂

Now it’s just onto the other three problems you’ve solved but I haven’t; 88, 124 and 126.

Nice solution, and yes it is a nice problem to approach in a recursive way. I think our approach is very similar, but the code structure is different. And yours is probably nicer 🙂

Well, at least the code is nice. For the sake of completeness, here it is:

And it runs very fast. My computer can’t even measure the time of it, always returns 0ms. 🙂

But it’s not too surprising. [latex]i, r, b[/latex] are all in the range 0..15 (inclusive), so there are only [latex]16^3 = 4096[/latex] distinct states (and even less than that, since [latex]r+b \leq 15[/latex]. My solution (which does recompute some subproblems) calls the function only 65535 times, so visits each possible state about 16 times.

But that is still much more computation than your solution does.

I just finished a Project Euler problem that can be solved with the exact same technique (along with another technique on top of that). For those interested and maybe want to try it out, it is problem number:

You can prove by induction that the coefficients (a_k) of the polynomial p(x)=(x+1)(x+2)…(x+n) give indeed the number of cases in which we obtain k blue disks. In fact you prove that implicitly.

Here is R code to obtain the probability of winning. I calculate it directly by summing probabilities of all outcomes where eight or more blue discs are drawn.

v <- rep(1,15) / 2:16

k <- lapply(8:15,combn,x=1:15)

sum(

unlist(

lapply(k,function(x) {

apply(x,2,function(y) {

z <- v

z[-y] <- abs(1-z[-y])

prod(z)

})

})

)

)

I solved the problem using straight probability theory, summing up probabilities of choosing each disk. I got the right answer, but it did a lot more operations than your and took about 1.8 ms. Then used your information about the polynomial representation, along with a C++ polynomial class I wrote and found the solution in .04 ms. I like your solution better.