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:

1 2 |
<em><u>B</u> * <u>B-1</u> = 0.5 T T-1</em> |

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.

Dim A As Double, B As Double

Dim TEMP As Double

Dim Answer As Double, T As Single

T = Timer

A = 1

B = 1

Do Until Answer > Sqr(0.5) * 10 ^ 12

TEMP = A

A = A + 2 * B

B = TEMP + B

Answer = (B + 1) / 2 ‘the number of Blue discs

Loop

Debug.Print Answer; ” Time:”; Timer – T

End Sub

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

Michael – very neat. At first I couldn’t see why precision would be a problem, but a quick trial verified that it was. I then set up a brute force solution using the decimal data type, which gives 28 significant figures, but is way too slow unless you know the answer already.

FWIW, here is the code using decimals:

Dim B As Variant, n As Variant, Prob As Variant, ans(1 To 3, 1 To 1) As Double

Dim b_2 As Variant, n_2 As Variant, Remain As Variant

ans(3, 1) = Timer

n = Min_discs – 1

Do

n = CDec(Int(n + 1))

B = CDec(Int(n * 0.5 ^ 0.5) – 1)

Do

B = B + 1

Prob = ((B / n) * ((B – 1) / (n – 1)))

Loop While Prob – 0.5 .LT. 0

b_2 = CDec(B * (B – 1))

n_2 = CDec(n * (n – 1))

Remain = n_2 – (2 * b_2)

Loop While Abs(Remain) .GT. 0

ans(1, 1) = B

ans(2, 1) = n

ans(3, 1) = Timer – ans(3, 1)

P_100 = ans

End Function

I noticed that the solutions are a factor of 5.8284 (sqrt(2) + 1)^2 apart. So, using the “don’t try this” method of adding one to T or B actually works if you can get a good guess for the start. Values for T are 21, 120, 697, 4060, 23661, …