Simpson's Rule in OCaml

(*
   An implementation of Simpson's Rule in OCaml
   Cf.
http://www.ccs.neu.edu/home/dorai/t-y-scheme/t-y-scheme-Z-H-22.html#node_sec_C.1
*)

let pi = 4.0 *. atan 1.0;;

let even n = (n mod 2) = 0;;
let odd n = not (even n);;

let rec integrate_simpson f a b n =
  if odd n then integrate_simpson f a b (n + 1)
  else
    let h = ((b-.a)/.(float_of_int n)) in
    let sum_every_other_ordinate_starting_from x0 num_ordinates =
      let rec loop x i r =
	if i >= num_ordinates then r else
	loop (x+.h*.2.) (i+1) (r+.(f x))
      in
      loop x0 0 0.
    in
    (((f a +. f b)
	+.(2.*.sum_every_other_ordinate_starting_from (a+.h*.2.) ((n/2)-1))
	+.(4.*.sum_every_other_ordinate_starting_from (a+.h)      (n/2)))
       *.(1./.3.)
       *.h);;

let phi x =
  (1. /. (sqrt (2. *. pi)))
    *.  exp (-. 0.5 *. x *. x);;

integrate_simpson phi 0. 1. 10;; (* 0.3413 *)
integrate_simpson phi 0. 2. 10;; (* 0.4772 *)
integrate_simpson phi 0. 3. 10;; (* 0.4987 *)