Project Euler #27

#load "nums.cma"

open Big_int;;

type triple = { a: big_int; 
                b: big_int; 
                len: big_int };;

let is_prime n =
  if lt_big_int n zero_big_int then
    false
  else
    let lim = sqrt_big_int n in
    let rec check_mod i l =
      if gt_big_int i  l then
        true
      else
        if eq_big_int (mod_big_int n i) zero_big_int then
          false
        else
          true && (check_mod (succ_big_int i) l)
    in
    check_mod (big_int_of_int 2) lim;;
      
let primes n =
  let rec inner_primes i =
    if gt_big_int i n then
      []
    else
      if is_prime i then
        [i] @ inner_primes (succ_big_int i)
      else
        inner_primes (succ_big_int i)
  in
  inner_primes (big_int_of_int 2);;

let make_form a b = 
  fun n -> 
    add_big_int 
      (mult_big_int n n) 
      (add_big_int (mult_big_int a n) b);;

let sequence_length f =
  let rec advance_sequence n =
    if not (is_prime (f n)) then
      zero_big_int
    else
      add_big_int 
        unit_big_int (advance_sequence (succ_big_int n))
  in
  advance_sequence (big_int_of_int 0);;

let triple_max a b =
  if gt_big_int a.len b.len then
    a
  else
    b;;

let find_longest_sequence i1 i2 n =
  let rec make_lengths_list i =
    if gt_big_int i i2 then
      []
    else
      [{ a=i; b=n; len=sequence_length (make_form i n); }] @ 
      (make_lengths_list (succ_big_int i))
  in
  List.fold_left (triple_max) 
    { a=zero_big_int; b=zero_big_int; len=zero_big_int; } 
    (make_lengths_list i1);;

let bs = primes (big_int_of_int 1000);;
let a1 = big_int_of_string "-999";;
let a2 = big_int_of_string "999";;

let max_ab = List.fold_left (triple_max) 
    { a=zero_big_int; b=zero_big_int; len=zero_big_int; } 
    (List.map (find_longest_sequence a1 a2) bs);;

Printf.printf "%s\n" 
    (string_of_big_int (mult_big_int max_ab.a max_ab.b));;