(** Constants *)
let pi = acos (-1.)
let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)
let factmax = 150


(** Generalized Boys function. Uses GSL's incomplete Gamma function.
    maxm : Maximum total angular momentum
*)
let boys_function ~maxm t =
  begin
    if t = 0. then
      Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1))
    else
      let incomplete_gamma ~alpha x = 
        Gsl.Sf.gamma alpha *. (  Gsl.Sf.gamma_inc_P alpha x )
      in
      let t_inv = 
        1. /. t
      in
      let factor =
        Array.make (maxm+1) (0.5, sqrt t_inv);
      in
      for i=1 to maxm
      do
        let (dm, f) = factor.(i-1) in
        factor.(i) <- (dm +. 1., f *. t_inv);
      done;
      Array.map (fun (dm, f) -> (incomplete_gamma dm t ) *. 0.5 *. f) factor 
  end



let fact_memo =
  let rec aux accu_l accu = function
  | 0 ->  aux [1.] 1. 1
  | i when (i = factmax) ->
          let x = (float_of_int factmax) *. accu in
          List.rev (x::accu_l)
  | i ->  let x = (float_of_int i) *. accu in
          aux (x::accu_l) x (i+1)
  in
  aux [] 0. 0
  |> Array.of_list



(** Factorial function. 
    @raise Invalid_argument for negative or arguments >100. *)
let fact = function
  | i when (i < 0) ->
    raise (Invalid_argument "Argument of factorial should be non-negative")
  | i when (i > 150) ->
    raise (Invalid_argument "Result of factorial is infinite")
  | i -> fact_memo.(i)



(** Integer powers of floats *)
let rec pow a = function
  | 0 -> 1.
  | 1 -> a
  | 2 -> a *. a
  | 3 -> a *. a *. a
  | -1 -> 1. /. a
  | n when (n<0) -> pow (1./.a) (-n)
  | n ->
    let b = pow a (n / 2) in
    b *. b *. (if n mod 2 = 0 then 1. else a)
;;