(* Project Euler problem 75 *)
let mlimit = (int_of_float(750000. ** 0.5));;
let limit = 1500000;;
 
let rec gcd a b =
  if a = b then a
  else
    if a < b then gcd a (b - a) else gcd (a - b) b;;
 
let perim m n =  (* For use with m and n coprimes with m > n  *)
  let a = m * m - n * n in
  let b = 2 * m * n in
  let c = m * m + n * n in
  a + b + c ;;
 
let is_co_prime m n = gcd m n = 1;;
 
let answer =
  let rec aux m acc =
    if m = mlimit then (List.sort compare acc)
    else
    let rec aux2 n accu=
      if n >= m then aux (m + 1) accu
      else
        if ((n + m) mod 2 = 1) && is_co_prime m n then
        let rec aux3 p accum =
          if p > limit then aux2 (n + 1) accum
          else
            let p2 = 2 * p in
            aux3 p2 (p::accum)
         in aux3 (perim m n) accu
         else aux2 (n + 1) accu
    in aux2 1 acc
  in aux 2 [];;
 
let prob75 =
  let rec aux lst acc =
    match lst with
    | [] -> acc
    | hd::tl when (List.exists (fun x -> x = hd) tl) -> aux (List.filter (fun x -> x <> hd) tl) acc
    | hd::tl -> aux tl (acc + 1)
  in aux answer 0;;
 
print_int prob75; print_newline();;

KCogUHJvamVjdCBFdWxlciBwcm9ibGVtIDc1ICopCmxldCBtbGltaXQgPSAoaW50X29mX2Zsb2F0KDc1MDAwMC4gKiogMC41KSk7OwpsZXQgbGltaXQgPSAxNTAwMDAwOzsKCmxldCByZWMgZ2NkIGEgYiA9CiAgaWYgYSA9IGIgdGhlbiBhCiAgZWxzZQogICAgaWYgYSA8IGIgdGhlbiBnY2QgYSAoYiAtIGEpIGVsc2UgZ2NkIChhIC0gYikgYjs7CgpsZXQgcGVyaW0gbSBuID0gICgqIEZvciB1c2Ugd2l0aCBtIGFuZCBuIGNvcHJpbWVzIHdpdGggbSA+IG4gICopCiAgbGV0IGEgPSBtICogbSAtIG4gKiBuIGluCiAgbGV0IGIgPSAyICogbSAqIG4gaW4KICBsZXQgYyA9IG0gKiBtICsgbiAqIG4gaW4KICBhICsgYiArIGMgOzsKCmxldCBpc19jb19wcmltZSBtIG4gPSBnY2QgbSBuID0gMTs7CgpsZXQgYW5zd2VyID0KICBsZXQgcmVjIGF1eCBtIGFjYyA9CiAgICBpZiBtID0gbWxpbWl0IHRoZW4gKExpc3Quc29ydCBjb21wYXJlIGFjYykKICAgIGVsc2UKICAgIGxldCByZWMgYXV4MiBuIGFjY3U9CiAgICAgIGlmIG4gPj0gbSB0aGVuIGF1eCAobSArIDEpIGFjY3UKICAgICAgZWxzZQogICAgICAgIGlmICgobiArIG0pIG1vZCAyID0gMSkgJiYgaXNfY29fcHJpbWUgbSBuIHRoZW4KICAgICAgICBsZXQgcmVjIGF1eDMgcCBhY2N1bSA9CiAgICAgICAgICBpZiBwID4gbGltaXQgdGhlbiBhdXgyIChuICsgMSkgYWNjdW0KICAgICAgICAgIGVsc2UKICAgICAgICAgICAgbGV0IHAyID0gMiAqIHAgaW4KICAgICAgICAgICAgYXV4MyBwMiAocDo6YWNjdW0pCiAgICAgICAgIGluIGF1eDMgKHBlcmltIG0gbikgYWNjdQogICAgICAgICBlbHNlIGF1eDIgKG4gKyAxKSBhY2N1CiAgICBpbiBhdXgyIDEgYWNjCiAgaW4gYXV4IDIgW107OwoKbGV0IHByb2I3NSA9CiAgbGV0IHJlYyBhdXggbHN0IGFjYyA9CiAgICBtYXRjaCBsc3Qgd2l0aAogICAgfCBbXSAtPiBhY2MKICAgIHwgaGQ6OnRsIHdoZW4gKExpc3QuZXhpc3RzIChmdW4geCAtPiB4ID0gaGQpIHRsKSAtPiBhdXggKExpc3QuZmlsdGVyIChmdW4geCAtPiB4IDw+IGhkKSB0bCkgYWNjCiAgICB8IGhkOjp0bCAtPiBhdXggdGwgKGFjYyArIDEpCiAgaW4gYXV4IGFuc3dlciAwOzsKICAKcHJpbnRfaW50IHByb2I3NTsgcHJpbnRfbmV3bGluZSgpOzs=