Generating All Pythagorean Triples

As per ProjectEuler 75, we would like to generate all integer triples (a, b, c) such that \(a^2 + b^2 = c^2\) with \(a + b + c \leq N\)

These triples are called Pythagorean triples. There is a wide variety of algorithms to generate these as discussed on Wikipedia. We use Euclid’s formula to generate primitive triples, i.e. triples such that the greatest common divisor of a, b, and, c is 1. Primitive triples are uniquely generated by iterating on two integers m and n, such that \(0 < n < m\), \(gcd(m, n) = 1\), and \(m - n\) is odd. Then, \(a = 2mn\), \(b = m^2 - n^2\), and \(c = m^2 + n^2\) define a Pythagorean triple.

Note that c is always the largest element of the triple, but both a and b can be the smallest element depending on the relative values of m and n. This generator can be implemented with the following simple ocaml code:

let rec gcd x y = if y = 0 then x else gcd y (x mod y)
let fold_prim init f max_s =
  let rec loop acc m n =
    let m2 = m * m in
    let mn = m * n in
    if m <= n || max_s < 2 * (m2 + mn) then
      let m = m + 1 in
      if max_s < 2 * (m2 + m) then acc
      else loop acc m (1 + m land 1)
    else if gcd m n &lt;&gt; 1 then loop acc m (n+2)
    else
      let n2 = n*n in
      let acc = f acc (m2 - n2) (2 * mn) (m2 + n2) in
      loop acc m (n+2)
  in
  loop init 2 1

The fold_prim function takes as input an initial element init, a function f that will be applied to all primitive Pythagorean triple (a, b, c) such that \(a + b + c \leq max_s\). It returns:

f (... f (f init a1 b1 c1) a2 b2 c2 ... ) aN bN cN

And so the type of fold_int is:

fold_int: `init:`a -> f:(`a -> int -> int -> int -> `a) -> int -> `a

Generating primitive triples in a unique way, i.e. each triple only corresponds to one pair (m, n) and appears once, gives us a way to generate all triples in a unique way too: we only have to iterate over some positive integer k and return (k.a, k.b, k.c) for every positive k and every primitive triple (a, b, c).

We could easily apply this to count the number of Pythagorean triples, including non-primitive ones, below some limit. For example using the code below returns the number of triples below \(10^8\). It runs in ~2s on my modestly powered ChromeBook.

let () =
  let s = 100_000_000 in
  let nb =
    fold_prim 0 (fun acc a b c -> acc + max_n / (a + b + c)) s
  in
  Format.printf "%d\n" nb

Note that in the implementation of fold_prim, we generate all the necessary pairs (m, n) and filter the non coprime ones using Euclid’s algorithm for GCD. This is quite inefficient and could be improved using a Stern-Brocot tree for example.