Project Euler 100: Finding the number of blue discs for which there is 50% chance of taking two blue.

So now we are at Problem 100 at Project Euler, you could call this an anniversary and I would have expected something extra difficult.  It took a while for me to crack the following problem

If a box contains twenty-one coloured discs, composed of fifteen blue discs and six red discs, and two discs were taken at random, it can be seen that the probability of taking two blue discs, P(BB) = (15/21)x(14/20) = 1/2.

The next such arrangement, for which there is exactly 50% chance of taking two blue discs at random, is a box containing eighty-five blue discs and thirty-five red discs.

By finding the first arrangement to contain over 10¹² = 1,000,000,000,000 discs in total, determine the number of blue discs that the box would contain.

But once you realize that a little something you can solve it quite easily.

We want to solve the following equation

where n is the total number of discs and b is the number of blue discs.  The  equation can be rearranged as

Finding an integer solution to this thing? Doesn’t it look like a diophantine equation to you? It has two variables and we are looking for an integer solution. My best guess was to search google for a “quadratic diophantine equation”, and that result turned up Dario Alpern’s Generic two integer variable equation solver.

That site can give you the detailed step to give you a basic solution to the equation, which we already have through the problem statement. As well as an algorithm for all other solutions to the equation. You should really go on to study the solutions for these as I wont go through them. I can’t give you better explanations than he does. The site gives us that

From there on, it is just to iterate solutions until you hit an n which is larger than the minimum value. So the solution to this problem can be implemented in C# as

long b = 15;
long n = 21;
long target = 1000000000000;

while(n < target){
    long btemp = 3 * b + 2 * n - 2;
    long ntemp = 4 * b + 3 * n - 3;

    b = btemp;
    n = ntemp;
}

yep it is that simple…. The code runs in

There are 756872327473 blues
Solution took 0,0006 ms

So a really fast solution for the problem. Alternatively it can be solved in Mathematica with the line of code

Solve[b / (b + r) * (b - 1) / (b + r - 1) == 1 / 2 && b + r > 10^12 && b > 0 && r > 0, {b, r}, Integers] // First

This code runs in 156ms.

Wrapping up

This was the 100th solved problem for me, and I want to make a small fanfare for myself. But I don’t think you want to hear me making that kind of noise.

You can download the source code here.  How did you solve the problem?