Euler Problem 100 asks:

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.

Using R for #Red, B for #Blue and T for total (R+B), the basic equation here is:

Since they can neither equal each other nor both be the square root of 0.5 (or 0.707106781…) we want B/T to be a skosh over the square root of 0.5, and (B-1/T-1) a skosh under. In other words, the answer will be the first integer B greater than Sqr(0.5)*10^12.

If we rearrange the equation, we get

• 2B^2 – 2B = T^2 – T

A single double quadratic equation with two unknowns can be morphed into a Pell equation with either no solutions, or an infinite number of them. Again, Euler wants the first solution with total discs T greater than 10^12. The first hit on a Google search for probability and Pell (here) returns an article exactly on point, with the examples for P=0.5 worked out in the first 11 cases. Dr. Sasser must have written it in a hurry, because it can be confusing. He uses two variables y, and they mean different things. I’m going to try to sort that out. There are also a few typos, but you can figure them out from context. The above equation does morph, as Dr. Sasser shows how, into Pell-form x^2 – Dy^2 = C as

• x^2 – 2y^2 = -1

with y = 2B-1. B, the number of Blue discs, is then (y+1)/2. Once we have the y solutions, we have the B solutions. The solutions for this Pell equation are x + y*Sqr(2) = (1 + sqr(2))^n, for integer n.

Flashing back to algebra, the product of (aw + bz)(cw + dz) is:

• acw^2 + (ad + cb)wz + dbz^2

If c = 1, w = 1, d = 1, z = Sqr(2), and z^2 = 2 corresponding to solution 1*1 + 1*Sqr(2) then rearranging the above gives:

• a + 2b + (a+b)*Sqr(2)

This is of the form x + y*Sqr(2) with x = a + 2b and y = a + b. This gives us a method to take (1 + Sqr(2)) to any integer power n. Recall that we’re looking for the y = a + b that gives B = (y+1)/2 greater than Sqr(0.5)*10^12.

This is the code that does that. It runs in a blink.

The first time I attempted this problem I incremented B or T as necessary to get just above Sqr(0.5)*10^12. Don’t go that way. Floating arithmetic errors give “false positives” for that approach.

There usual angle bracket corrections are in the code.

…mrt