# Micropuzzles and generating Pythagorean triples

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.

1. Floating point numbers will always catch you out even when you are looking out for them to catch you out.
2. There’s more than one way to generate Pythagorean triples and some are quicker than others.

• Be respectful. You may want to read the comment guidelines before posting.
• You can use Markdown syntax to format your comments. You can only use level 5 and 6 headings.
• You can add class="your language" to code blocks to help highlight.js highlight them correctly.

Privacy note: This form will forward your IP address, user agent and referrer to the Akismet, StopForumSpam and Botscout spam filtering services. I don’t log these details. Those services will. I do log everything you type into the form. Full privacy statement.

• #870

I learned a lot following along with your discussion and code. I know it's been a few years since you wrote this blog entry, but wanted to add some comments to correctly guide other readers.

Your isInt function should be changed, one suggestion: isInt x = abs (x - (fromIntegral \$ round x)) < 0.000001

Ceiling will always round UP, so for example ceiling 1.01 gives 2. The abs function ensures that things work whether the "nearly" integer value is slightly above or below its neighboring integer.

Regarding your Euclidean approach: it is incorrect that the method only produces primitives. For example m=3, n=1 produces (6,8,10). However, it is true that the triplets produced this way is a strict subset of all the triplets, and that the solution to the problem is NOT contained in this subset.

Final note: Floats can be avoided in finding the triplets. At a cost of speed loss, the algorithm can search to find hypotenuse lengths satisfying x^2+y^2==z^2, rather than calculating z. Of course that doesn't help in locating those meeting the perimeter and area restriction.