(* Copyright (c) 2008 Gary A. Huber, Howard Hughes Medical Institute 
   See the file COPYRIGHT for copying permission
*)


(* Approximate the definite integral of f from a to b by Romberg's method.
    eps is the desired accuracy.
*)

let isum f nb ne = 
  let sum = ref 0.0 in
    for i = nb to ne do
      sum := !sum +. (f i)
    done;
    !sum
;;

let rec pow x n = 
  if n = 0 then 
    1
  else (
    if (n mod 2) = 0 then
      let xnh = pow x (n/2) in
	xnh*xnh
    else
      x * (pow x (n-1))
  )
;;

let integral ~f ~low:a ~high:b ~accuracy:eps = 

  let rec loop n rnm = 
    let h = (b -. a)/.(float (pow 2 n)) in
    let rn = Array.init (n+1) (fun i -> nan) in 
      rn.(0) <- 0.5*.rnm.(0) +. h*.(isum (fun k -> f (a +. (float (2*k-1))*.h)) 1 (pow 2 (n-1)));
      for m = 1 to n do
	rn.(m) <- rn.(m-1) +. (rn.(m-1) -. rnm.(m-1))/.( (float ((pow 4 m) - 1)));
      done;

      if (abs_float ((rnm.(n-1) -. rn.(n))/.rn.(n))) < eps then
	rn.(n)
      else
	loop (n+1) rn
  in

  let r0 = [| 0.5 *. (b -. a) *. (f(a) +. f(b))|]  in 
    loop 1 r0
;;


(**********************************************************************)
(* Test code *)
(*
Printf.printf "%f\n" 
  (integral (fun x -> let x2 = x*.x in x2*.x2*.x2) 0.0 1.0 1.0e-8);;
*)

(*
let pi = 3.1415926;;

Printf.printf "%f\n"  (integral (fun t -> 2.0/.(sqrt pi)*.(exp (-.t*.t))) 0.0  1.0 1.0e-8);;
*)