|
| 1 | +(* |
| 2 | + * Solution to Project Euler problem 100 |
| 3 | + * Copyright (c) Project Nayuki. All rights reserved. |
| 4 | + * |
| 5 | + * https://www.nayuki.io/page/project-euler-solutions |
| 6 | + * https://github.com/nayuki/Project-Euler-solutions |
| 7 | + *) |
| 8 | + |
| 9 | + |
| 10 | +(* |
| 11 | + * Suppose the box has b blue discs and r red discs. |
| 12 | + * The probability of taking 2 blue discs is [b / (b + r)] * [(b - 1) / (b + r - 1)], |
| 13 | + * which we want to be equal to 1/2. Rearrange the equation: |
| 14 | + * [b(b - 1)] / [(b + r)(b + r - 1)] = 1 / 2. |
| 15 | + * 2b(b - 1) = (b + r)(b + r - 1). |
| 16 | + * 2b^2 - 2b = b^2 + br - b + br + r^2 - r. |
| 17 | + * b^2 - b = r^2 + 2br - r. |
| 18 | + * b^2 - (2r + 1)b + (r - r^2) = 0. |
| 19 | + * Apply the quadratic equation to solve for b: |
| 20 | + * b = [(2r + 1) +/- sqrt((2r + 1)^2 - 4(r - r^2))] / 2 |
| 21 | + * = r + [1 +/- sqrt(8r^2 + 1)]/2 |
| 22 | + * = r + [sqrt(8r^2 + 1) + 1]/2. (Discard the minus solution because it would make b < r) |
| 23 | + * |
| 24 | + * For b to be an integer, we need sqrt(8r^2 + 1) to be odd, and also 8r^2 + 1 be a perfect square. |
| 25 | + * Assume 8y^2 + 1 = x^2 for some integer x > 0. |
| 26 | + * We can see this is in fact a Pell's equation: x^2 - 8y^2 = 1. |
| 27 | + * |
| 28 | + * Suppose we have the solution (x0, y0) such that x0 > 0 and x0 is as small as possible. |
| 29 | + * This is called the fundamental solution, and all other solutions be derived from it (proven elsewhere). |
| 30 | + * Suppose (x0, y0) and (x1, y1) are solutions. Then we have: |
| 31 | + * x0^2 - 8*y0^2 = 1. |
| 32 | + * (x0 + y0*sqrt(8))(x0 - y0*sqrt(8)) = 1. |
| 33 | + * (x1 + y1*sqrt(8))(x1 - y1*sqrt(8)) = 1. (Similarly) |
| 34 | + * Multiply them together: |
| 35 | + * [(x0 + y0*sqrt(8))(x0 - y0*sqrt(8))][(x1 + y1*sqrt(8))(x1 - y1*sqrt(8))] = 1 * 1. |
| 36 | + * [(x0 + y0*sqrt(8))(x1 + y1*sqrt(8))][(x0 - y0*sqrt(8))(x1 - y1*sqrt(8))] = 1. |
| 37 | + * [x0*x1 + x0*y1*sqrt(8) + x1*y0*sqrt(8) + 8y0*y1][x0*x1 - x0*y1*sqrt(8) - x1*y0*sqrt(8) + 8y0*y1] = 1. |
| 38 | + * [(x0*x1 + 8y0*y1) + (x0*y1 + x1*y0)*sqrt(8)][(x0*x1 + 8y0*y1) - (x0*y1 + x1*y0)*sqrt(8)] = 1. |
| 39 | + * (x0*x1 + 8y0*y1)^2 - 8*(x0*y1 + x1*y0)^2 = 1. |
| 40 | + * Therefore (x0*x1 + 8y0*y1, x0*y1 + x1*y0) is also a solution. |
| 41 | + * By inspection, the fundamental solution is (3, 1). |
| 42 | + *) |
| 43 | + |
| 44 | +(* Fundamental solution *) |
| 45 | +x0 = 3; |
| 46 | +y0 = 1; |
| 47 | + |
| 48 | +(* Current solution *) |
| 49 | +x = x0; |
| 50 | +y = y0; (* An alias for the number of red discs *) |
| 51 | + |
| 52 | +While[True, |
| 53 | + (* Check if this solution is acceptable *) |
| 54 | + sqrt = Sqrt[y^2 * 8 + 1]; |
| 55 | + blue = (sqrt + 1) / 2 + y; |
| 56 | + If[OddQ[sqrt] && blue + y > 10^12, |
| 57 | + ans = blue; |
| 58 | + Break[]]; |
| 59 | + nextx = x * x0 + y * y0 * 8; |
| 60 | + nexty = x * y0 + y * x0; |
| 61 | + x = nextx; |
| 62 | + y = nexty;] |
| 63 | +ans |
0 commit comments