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^{12}= 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?

Congratulations! The first 100 Project Euler posts were amazing, and I’m looking forward to the next 100.

Now go celebrate! You deserve it!

Thank you. I am amazed that I have kept on doing it for so long. However, I am still enjoying it not least due to the interesting discussions that keeps going on here on the site.

Yes, congratulations. Very good job with those problems :D. Good luck on another 100.

Thanks. For some reason I doubt that the next 100 will be solved and posted linearly, but lets see. I have already learned a million tricks and mathematical techniques from the first 100 problems, so who know what I might utilize to solve the next 100.

Line 2 contains an error, as n should = 21.

In addition, the output posted shows the result using the erroneous value for n, but produces correct results once the value is fixed.

Thanks for noticing. Kristian has now fixed it.

Which I also tried to write, but apparently that comment was eaten by the comment trolls. Anyway Bjarki summed it up just fine. And thanks for noticing.

Hi,

Google brought me here after i failed to solve this problem. I’m not great at mathematics, but created a loop in R and got another number: 707106781187

I validate the outcome using:

calculateprobability = function (blueballs, totalballs) (

(blueballs/totalballs)*((blueballs-1)/(totalballs-1))

)

This gives me:

> calculateprobability (707106781187, 10^12)

[1] 0.5

and

> calculateprobability (756872327473, 10^12)

[1] 0.5728557201

Isn’t that strange? What am I doing wrong?

Sorry for the long response time. Without knowing the inner workings of R, that sounds to me like a problem with floating point operations.

That is brilliant.. thank u so much i am not good in mathematics.. you have helped me a lot when i was in a trouble with maths..

You are welcome, I am just happy that people learn from the posts.

Now i am late withnthe response :-), thanks for answer. Indeed looks like floating point problem. But this problem is also a bit to har for me. I will start with the less complex prblems.

Cheers and happy new year,

Pieter

I came up with the answer:

707106783028.0*707106783027.0/1000000002604.0/1000000002603.0

which seems to be right. But gets the big cross on Euler. I started at 10**12 and increased one at a time and solved a quadratic equation for the number of blue discs. Do you think this is a valid answer?

Seems like a rounding error to me. Have you tried checking it with integer arithmetic, such as biginteger and checking that the numerator is half the size of the denominator.

You’re right. Thanks.

Hi

I wonder why I don’t obtain the correct value, I’m probably missing something.

Assume we let where is your calculated value. If I use the calculated value in the equation below (with help of wolframalpha) I get

.

According to wolframalpha the correct value is

which is approx

Thanks for your help and once again, great blog!

/Robin

P.S (please edit the latex code in the post if I messed it up, and remove this row – thanks!)

Hi,

First off all I wanna thank you for your great work and explanation of all Project euler problems.

Second, for this problem 100, I dont want to discuss the method but I want an explanation why my answer is wrong :

I solved this problem by brute force and i obtained :

– 707106788769 blue discs

– 292893221954 red discs

– total = 1000000010723

the total number of discs >= 10^12 and my answer is much smaller than yours !! Do I missed something in this problem ???

Thank you

Are you sure that you use integer values, because the most likely error is that you are using floating point and then get a rounding error

Actually I wrote a function that return the total number of discs (y) from the desired blue discs (x) with these relations :

2x(x-1) = y(y-1)

2x^2 – 2x = y^2 – y

2x^2 – 2x = (y – 1/2)^2 -1/4

sqrt(2x^2 – 2x + 1/4) = y – 1/2

sqrt(2x^2 – 2x + 1/4) + 1/2 = y

So the first value of X that gives an integer Y > 10^12 will be the right answer. This is my python code :

def equation (x):

return (math.sqrt(2*x**2 – 2*x + 1/4) + 1/2)

x = 7.071067 * 10**11

result = 0

while result == 0 :

y = equation(x)

if y >= 10**12 and int(y)==y :

print (x)

break

i=i+1

Even if the math.sqrt() function returns a float value, I am checking this case later (“if y >= 10**12 and int(y)==y”).

I took a very different approach. After playing around for a bit I started looking for a pattern. Take the sequence of solutions n[i]. It turns out the sequence n[i+1]/n[i] is a monotonically increasing converging sequence. The ratio is about 5.82. Thus having found a value of n the next solution will be for a value n[i+1] just over 5.82*n[i].

Implementation was to take n[0] = 21 and start the search at n = 119 looking for 2(b-1)b == (n-1)n where b = int(2^0.5 * n) + 1. This recognizes the to get (b/n)((b-1)/(n-1)) to be 1/2 one of the ratios has to be just above 0.5^0.5 and one just below that number.

Having found a case where the ration is 1/2 the next search starts at n[i]*n[i]/[n[i-1].

The algorithm ended up looking at about 170 different values of n altogether before arriving at the solution.