include Constants 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) ;; (** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *) let chop f g = if (abs_float f) < cutoff then 0. else f *. (g ())