Project Euler Solution 86: Cuboid Route
Cube
Let us first start with a simpler problem, namely that of a cube. There we have only one side length, so $a = b = c$.
There are six possible ways to draw the shortest path. By symmetry we know that it has to intersect the one edge exactly in the middle. The cube has a lot of symmetries, and so the solutions need to have the same symmetries. Therefore we know that all of them have equal length and also that they need to intersect the edge in the middle. The paths sketched look like this:
I gave each path a name. For instance the one that goes through the $a$-$b$-plane first and then through the $a$-$c$-plane is called $(ba,ac)$. This is the one that we will look into further.
We can write the length of the path as this: $$ l = \sqrt{\left(\frac a2 \right)^2 + a^2} + \sqrt{\left(\frac a2 \right)^2 + a^2} = \sqrt 5 \cdot a \,. $$
Let us prove that this is indeed the shortest possible path. Let us call the intersection point $x$, and we need to have $0 \leq x \leq a$ because this intersection point cannot be outside of the cube. Using this $x$, we have the generalized length of the path as a function of $x$: $$ l(x) = \sqrt{x^2 + a^2} + \sqrt{(a - x)^2 + a^2} \,. $$
We want to find the minimum $x$ which minimizes $l$ and expect it to be $x = a/2$. In order to find the minimum, we will take the derivative of $l$ and find its roots (where it is zero): $$ \frac{\mathrm d l}{\mathrm d x} = \frac{2x}{2 \sqrt{x^2 + a^2}} + \frac{- 2 (a - x)}{2\sqrt{(a - x)^2 + a^2}} = 0 \,. $$
We need to do some algebra next. We move one summand to the other side of the equation. Then we also cancel the factor 2 already. We get this expression: $$ \frac{x}{\sqrt{x^2 + a^2}} = \frac{a - x}{\sqrt{(a - x)^2 + a^2}} \,. $$
Next we take the square of both sides. This brings us to this form: $$ \frac{x^2}{x^2 + a^2} = \frac{(a - x)^2}{(a - x)^2 + a^2} \,. $$
Then we multiply both sides by both denominators and get this equation: $$ x^2 \left( (a-x)^2 + a^2 \right) = (x^2 + a^2) (a-x^2) \,. $$
Expanding both sides completely and cancelling the terms, we end up with just these terms: $$ a^4 - 2 a^3 x = 0 \,. $$
This can be solved for $x$ to yield the solution: $$ x = \frac a2 \,. $$
So this is what I have proposed in the beginning. And now we have seen that is indeed the solution, as well as the way to find $x$ if we hadn't known it from a symmetry argument. These symmetry arguments are very powerful in physics, especially theoretical physics. That's why you will find theoretical physicists (like me) marvel about the beautiful symmetries in nature.
Back to the cuboid
Now that we have derived the shortest path length with the cube, let us go back to the original cuboid.
We will still use the route $(ba, ac)$, even though it is not necessarily the shortest one. But that's not a problem. Once we have the shortest path along that route we can just permute $a$, $b$ and $c$ to find all six routes. We just take the shortest from that and learn which permutation is generally the shortest one.
So the path length along that particular route is expressed like this: $$ l(x) = \sqrt{x^2 + b^2} + \sqrt{(a-x^2) + c^2} \,. $$
Now it is not obvious what this $x$ is supposed to be. We use the same approach and compute the root of the derivative, $$ \frac{\mathrm dl}{\mathrm dx} = 0 \,. $$
After similar manipulations, we end up with this equation: $$ (b-c^2) x^2 - 2ab^2x + a^2 b^2 = 0 \,. $$
We can use the solution formula and end up with the following: $$ x = \frac{ab (b \pm c)}{b^2 - c^2} \,.$$
This can be further simplified with a neat trick. Recall the binomial equations and that $b^2 - c^2$ can be written as $(b+c)(b-c)$. Then we can cancel one of them with the $(b \pm c)$ in the numerator and get $$ x = \frac{ab}{b \mp c} \,. $$
We have the constraint $0 \leq x \leq a$. Rewriting the solution as $$ x = a \frac{b}{b \mp c} $$ we can see that the version with minus will lead to $x > a$ and therefore is ruled out. We can therefore compactly write the solution as $$ x = \frac{ab}{b + c} \,. $$
Putting this back into the path length formula, we get $$ l(x) = \sqrt{\left(\frac{ab}{b+c}\right)^2 + b^2} + \sqrt{\left(a - \frac{ab}{b+c}\right)^2 + c^2} \,. $$
Optimal route
The problem that we have now is that we don't know which of the six routes is the optimal one. We need to try all six permutations.
For this we first implement the path length equation in Python:
def path_length(a: int, b: int, c: int) -> float: return math.sqrt(((a * b) / (b + c)) ** 2 + b**2) + math.sqrt( (a - (a * b) / (b + c)) ** 2 + c**2 )
And then we just compute the path lengths for the six permutations of the side lengths that are given in the problem statement:
{edges: path_length(*edges) for edges in itertools.permutations([3, 5, 6])}
We then get this expression:
{(3, 5, 6): 11.40175425099138, (3, 6, 5): 11.40175425099138, (5, 3, 6): 10.295630140987, (5, 6, 3): 10.295630140987, (6, 3, 5): 10.0, (6, 5, 3): 10.0}
It seems that the $a$ needs to be the longest edge for for this route to be the shortest one. And also it seems that the order of $b$ and $c$ doesn't matter. We can show that indeed the expression is unchanged when $b$ and $c$ are changed. For that one needs to square the expression $l(x)$ and carefully group all terms such that the symmetry is clear.
We can check our finding with more examples by generating random cuboids and checking whether that assumption holds:
for i in range(1000): lengths = [random.random() for j in range(3)] actual = min(path_length(*edges) for edges in itertools.permutations(lengths)) expected = path_length(*sorted(lengths, reverse=True)) assert abs(actual - expected) < 1e-10, (lengths, actual, expected)
And indeed it does. The tolerance in the comparison is needed because the symmetry in $b$ and $c$ doesn't hold numerically.
We therefore have a legitimate way to compute the smallest path length for a given cuboid $(a, b, c)$ when we just permute it such that $a \geq b \geq c$ and then use the formula that we have derived.
Integer path lengths
The problem asks us to find how many cubes with $M \geq a \geq b \geq c$ have a shortest path length which happens to be an integer value. So we need to come up with a way to decide that it is going to be an integer.
Let us look at the formula for the path length again: $$ l = \sqrt{\left(\frac{ab}{b+c}\right)^2 + b^2} + \sqrt{\left(a - \frac{ab}{b+c}\right)^2 + c^2} \,. $$
This is an integer if both radicands are perfect squares. But this is too strict and fails for the example with side lengths (6, 5, 3). We have this: $$ l = \sqrt{\frac{30^2}{8^2} + 5^2} + \sqrt{\left(6-\frac{30^2}{8^2}\right) + 3^2} = \sqrt{\frac{2500}{64}} + \sqrt{\frac{900}{64}} \,. $$
So the radicands are perfect squares, but perfect squares of fractions. We therefore need to expand the expressions into numerator and denominator.
We have to rewrite it like this, as explicit fractions: $$ \sqrt{\frac{a^2 b^2 + b^2 (b+c)^2}{(b+c)^2}} + \sqrt{\frac{\left(a (b+c) - ab\right)^2 + c^2 (b+c)^2}{(b+c)^2}} $$
From this we can derive a better check:
def shortest_path_is_integer(a, b, c) -> bool: numerator_1 = a**2 * b**2 + b**2 * (b + c) ** 2 sqrt_1 = int(math.sqrt(numerator_1)) if sqrt_1**2 != numerator_1: return False numerator_2 = (a * (b + c) - a * b) ** 2 + c**2 * (b + c) ** 2 sqrt_2 = int(math.sqrt(numerator_2)) denominator = b + c return sqrt_2**2 == numerator_2 and (sqrt_1 + sqrt_2) % denominator == 0
And then we can write the solution with that by just enumerating all $a \geq b \geq c$:
def solution() -> int: result = 0 for a in itertools.count(1): print(a, result) for b in range(1, a + 1): for c in range(1, b + 1): if shortest_path_is_integer(a, b, c): result += 1 if result > 1_000_000: return a
This quickly finds the given solution when we cap at 2,000. But with a cap at 2,000,000, it took 07:13 min to solve. There needs to be a more clever way.
After submitting the solution, I was able to look into the forums there. And there are a couple of very good ideas which can be used to make it much simpler.
Simplifying the path length
Our expression for the path length is quite complicated. There is a much easier way to get to it. Let us not take a look at the cuboid in 3D but rather unroll it to a 2D shape. Then we can immediately see that the diagonal has length $$ l = \sqrt{a^2 + (b + c)^2} \,. $$
That is the same as the more complicated expression, one can check that numerically. I tried to check that algebraically but eventually ran out of motivation to find a factorization for my sum of eight terms into a square.
We can do the same algorithm as before, just with a simpler equation. We can also try to improve the test for a square root a bit. Instead of using math.sqrt I use Heron's method with integers to provide the square root of an integer if it was a square.
def integer_root(radicand: int) -> Optional[int]: result = radicand while True: new = (result + radicand // result) // 2 if result == new: if result**2 == radicand: return result else: return None elif new > result: return None else: result = new
The test for the shortest path can then be cast in term of that $b + c$ thing and cached:
@functools.cache def shortest_path_is_integer(a, b_plus_c) -> bool: return integer_root(a**2 + b_plus_c**2) is not None
The main solution function then looks like this:
def solution() -> int: ceiling = 1_000_000 result = 0 for a in itertools.count(1): for b in range(1, a + 1): for c in range(1, b + 1): if shortest_path_is_integer(a, b + c): result += 1 if result > ceiling: return a
We iterate through all the boxes and perform the test. This finishes in 02:55 min, so it is much better than the previous solution.
Just to make sure that this integer root finding is an improvement, I have used our old is_square from Solution 66: Diophantine Equation:
@functools.cache def shortest_path_is_integer(a, b_plus_c) -> bool: return is_square(a**2 + b_plus_c**2)
This takes 02:50 minutes, so there is not that much difference between both implementations.
So this is approach is a bit better than my original one with the complicated path length formula, but it still doesn't fulfil the “one minute rule”. And given the age of the problem, this should work much faster on my modern laptop.
Counting multiplicities
The check whether it is an integer path length only depends on $a$ and $b + c$. There is no real need to check all combinations of $b$ and $c$, we just need to check once for their sum and then multiply it with the number of $(b, c)$ there are such that $M \geq a \geq b \geq c$.
I have just drawn a grid of the numbers $b$ (columns) and $c$ (rows) and filled in their sum. If we have $b \geq c$, I have marked the cell with green. This is for $a = 6$.
We can now see that there is one way to have $b + c = 2$, one way for 3, two ways for 4, and so on. This has a triangular shape, so we can find an easy expression for that.
To make sure that we get this right, let's start with the test. We want to have a multiplicities function that reproduces the above grid.
def test_multiplicity() -> None: expected = {2: 1, 3: 1, 4: 2, 5: 2, 6: 3, 7: 3, 8: 3, 9: 2, 10: 2, 11: 1, 12: 1} for b_plus_c, expected_multiplicity in expected.items(): assert multiplicity(6, b_plus_c) == expected_multiplicity
And then we can implement it:
def multiplicity(a: int, b_plus_c: int) -> int: if b_plus_c <= a + 1: return b_plus_c // 2 else: return (2 * a - b_plus_c + 2) // 2
Using this allows us to write the solution in two loops instead of three loops:
def solution() -> int: ceiling = 1_000_000 result = 0 for a in itertools.count(1): for b_plus_c in range(1, 2 * a + 1): if shortest_path_is_integer(a, b_plus_c): result += multiplicity(a, b_plus_c) if result > ceiling: return a
That finishes in 2512 ms, so it fulfils the “one minute rule” and is fast enough.
Since we only call the function shortest_path_is_integer once with each value of $b + c$, we can take out the cache. Without that cache we can get the solution in 1084 ms, which is even better.
Pythagorean triplets
There was another great hint in the forums: Since we look for integer $l$, we can reframe this into the search for Pythagorean triplets $(a, b+c, l)$.
We have already encountered those in Solution 75: Singular integer right triangles. Using this one can construct the triplets. The problem here is that the used way to generate triplets needs some sort of ceiling for $a$ and $b + c$. We have this $M$, but ideally we would just iterate all $b$ and $c$ before moving to the next $a$, making an explicit ceiling $M$ redundant.
The way that the Pythagorean triplets are generated with the multiplier $k$ means that we have to make $M$ explicit. We therefore have to write a function that computes all the solutions for all cuboids up to $M$. And then we would have to bisect in $M$ to find the one where the number of solutions exceeds 1,000,000.
As I was content with the above solution that runs in around 1 s, I didn't implement this myself.