引き続きOCaml

自然数rについて、x*x + y * y = rを満たす自然数(x, y)のリストを返す関数
ただし、x ≧ y ≧ 0

# let squares r = 
    let isqrt n = int_of_float (sqrt (float_of_int n)) in
    let y_max = isqrt (r / 2) in
    let rec get_list y list = 
      if y >= y_max then list
      else 
        let x = isqrt (r - y * y) in
        if x * x + y * y = r
          then (x, y) :: get_list (y+1) list
          else get_list (y+1) list
    in get_list 0 [];;

実行例

# squares 36;;
- : (int * int) list = [(6, 0)]
# squares 1024;;
- : (int * int) list = [(32, 0)]
# squares 48612265;;
- : (int * int) list =
[(6972, 59); (6971, 132); (6952, 531); (6948, 581); (6944, 627); (6917, 876);
 (6899, 1008); (6853, 1284); (6789, 1588); (6772, 1659); (6756, 1723);
 (6651, 2092); (6637, 2136); (6576, 2317); (6556, 2373); (6413, 2736);
 (6404, 2757); (6384, 2803); (6259, 3072); (6213, 3164); (6124, 3333);
 (6048, 3469); (5907, 3704); (5832, 3821); (5691, 4028); (5656, 4077);
 (5613, 4136); (5432, 4371); (5243, 4596); (5179, 4668); (5139, 4712);
 (5008, 4851)]