Micropuzzles is a charming book from the early 1980s by a chap called J.J. Clessa. It is a sort of forerunner of Project Euler, containing little puzzles aimed at computer hobbyists who have acquired those newfangled micros and are looking for something challenging to do with their machines.
The first main puzzle was quite interesting.
In coming up with a solution, I needed to generate Pythogorean triples, which are a series of 3 numbers, x,y and z such that x^2 + y^2 == z^2. They are Pythagorean in the sense that their relationship is the same as the sides of a right-angle triangle.
The puzzle goes like this
I’d like you to find the smallest-area Pythogorean triangle whose perimeter is a perfect square and whose area is a perfect cube.
Supporting functions
My implementation language was Haskell. Because yucky floating point numbers were involved I first needed an isInt function that would tell me if numbers were whole. My first attempt went like this.
isInt x = x == (fromIntegral $ round x)
That proved the source of some annoying glitches (thanks IEEE-754!) when I was working out if things were perfect cubes. That is because, to take 8 as an example
Prelude> 8 ** 3
512.0
Prelude> 512.0 ** (1/3)
7.999999999999999
So I refactored. Maybe it should be called isDamnNearInt
now.
isInt x = x - (fromIntegral $ ceiling x) < 0.000001
With a working way to detect if something was an integer, I made myself isSquare and isCube functions.
isCube x = isInt $ x ** (1/3)
isSquare = isInt . sqrt
Now I could take a shot at implementing a solution. As it turns out there are shitloads of ways to generate Pythagorean triples.
Take one: the Euclidian formula
I had gone down the route of using Euclid’s m/n formula — take any m and n with m>n and your x,y and z are 2mn, m^2 - n^2 and m^2 + n^2. Here that is in Haskell.
pTriplets i = [(x,y,z) | m <- [2..i]
, n <- [1..(m-1)]
, let x = m^2 - n^2
, let y = 2 * m * n
, let z = m^2 + n^2
]
But there’s a problem. This way of doing it doesn’t generate all possible Pythogorean triples.
It generates all primitive Pythagorean triples — like 3,4,5 — and many of the multiples — 6,8,10 is 3,4,5 with each number doubled, so it is called a multiple. But Euclid’s formula does not generate all of the multiples. And, importantly, it doesn’t generate the one we are after!
The brute-force approach
Back to the drawing board. I tried something much more obvious (at least to me). Just iterate through all x and y up to n and take the square root of the sum of their squares as z. Note the calls to isInt, isCube and isSquare. This method actually finds the solution.
pTriplets' n = [(x,y,z) | x <- [1..n]
, y <- [x..n]
, let z = sqrt (x^2 + y^2)
, isInt z
, isCube ((x*y)/2)
, isSquare (x+y+z)
]
And this remarkably straightforward implementation works as expected if you give it an n of 200 or so (the exact number would give the puzzle away).
Dickson’s method
What happened, though, was that I started reading about the other ways of generating Pythagorean triples. Many of them are similar to the Euclidian formula in that they generate only the primitive Pythagorean triples. One of the more straightforward methods that does generate them all is Dicksons’ method, which was first described by Leonard Eugene Dickson in 1920.
There is a nice proof that it generates all the triples, due to Josef Rukavicka in 2013.
Interestingly, Dickson’s method gave a significant performance boost, sometimes up to doubling the speed at which I could get the answer.
pTriplets'' n = [(x,y,z) | s <- [1..n]
, t <- [s..n]
, let r = sqrt (2 * s * t)
, let x = r + s
, let y = r + t
, let z = r + s + t
, isInt r
, isCube ((x*y)/2)
, isSquare (x+y+z)
]
Not too shabby.
Main> pTriplets' 192
[(ANSWER REDACTED)]
(0.12 secs, 68105184 bytes)
Main> pTriplets'' 192
[(ANSWER REDACTED)]
(0.05 secs, 47123304 bytes)
What did I learn?
Two things, I think.
- Floating point numbers will always catch you out even when you are looking out for them to catch you out.
- There’s more than one way to generate Pythagorean triples and some are quicker than others.