Uniquely generate integer triples (a, b, c) such that 2*a*a + b*b = c*c with a+b+c <= max_s.

pyth2.ml

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 n2 = n * n in
    let mn = m * n in
    if max_s < m2 then acc
    else if max_s < 2*n2 + mn then
      let n = int_of_float (float_of_int m /. 1.415) in
      loop acc (m+1) n
    else
      let acc =
        if gcd m n == 1 then
          if m land 1 = 1 then begin
            if m2 <= 2*n2 && 4*n2 + 2*mn <= max_s then
              f acc (2*mn) (2*n2 -m2) (m2 + 2*n2)
            else acc
          end else
            if m2 <= 2*n2 then
              let m2_2 = m2 / 2 in
              f acc mn (n2 - m2_2) (m2_2 + n2)
            else acc
        else acc
      in
      loop acc m (n + 1)
  in
  loop init 1 1

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