let find_pi tol =
  let pi = ref 0.0 in
  let rec find_pi a b t p =
    if abs_float (a -. b) < tol
    then !pi
    else (
      print_endline "iteration";
      let a_n1 = (a +. b) /. 2.0 in
      let b_n1 = sqrt (a *. b) in
      let t_n1 = t -. (p *. ((a -. a_n1) ** 2.0)) in
      let p_n1 = 2.0 *. p in
      pi := ((a_n1 +. b_n1) ** 2.0) /. (4.0 *. t_n1);
      find_pi a_n1 b_n1 t_n1 p_n1
    ) in
  find_pi 1.0 (1.0 /. (sqrt 2.0)) 0.25 1.0

let () =
  print_endline (string_of_float (find_pi (epsilon_float)))