Euler Project problems can be addictive. This one gave me fits.

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)×(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^12 = 1,000,000,000,000 discs in total, determine the number of blue discs that the box would contain.

I like to use Mathematica to work on these problems, but sometimes “a sufficiently advanced technology is indistinguishable from magic”. Here’s the magical solution:

FindInstance[2 x (x - 1) ==  y(y - 1) && x > 0 && y > 10^12, {x, y}, Integers]

This returned the solution {{x -> 756872327473, y -> 1070379110497}} but I have no idea how it was found.

So here’s a different approach. Consider the the two sequences {2 x (x – 1) } and { x (x – 1) }. We can just hunt for common elements using brute force, as follows:

foo @ n_ := n (n - 1)
foo2 @ n_ := 2 n (n - 1)

hunt[steps_] := Module[{j, k, x, y},
    j = 1; k = 1;
    x = foo@j; y = foo2@ j;
    step = 1;
    While[step < steps,
      If[x < y, x  = foo@j ; j += 1];
      If[x > y, y = foo2@ k; k += 1];
      If[x == y, Print @{x, j - 1, k - 1}; x = foo@j ; j += 1];
      step += 1;
      ]
    ]

Running hunt@100 produces

{0, 0, 0}
{0, 1, 0}
{12, 4, 3}
{420, 21, 15}
{14280, 120, 85}
{485112, 697, 493}

You can read off the solutions to the examples given in the problem statement. The trouble is that the brute force approach does not even begin to scale to exploring out around 10^12. So we turn to the On-Line Encyclopedia of Integer Sequences, also known as Sloane’s. Entering the initial terms 3, 15, 85, 493 immediately gets a hit.

Sloane’s provides the formula a(n)=6*a(n-1)-a(n-2)-2, with a(0)=1, a(1)=3. I don’t know where this formula came from, but I think the second approach was more enlightening than the initial magic approach.

A final line of approach. The initial setup is (x/y) (x-1)/(y-1) = 1/2. A bit of rearrangement gives 2 = y/x (y-1)/(x-1). In particular, y/x and (y-1)/(x-1) are rational approximations of the square root of two. I understand that continued fractions can be used to attack this kind of problem, and plan to look into that as time allows.