diff --git a/Basis/ERI.ml b/Basis/ERI.ml index b7324a2..7e13d3f 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -19,11 +19,19 @@ type t = (float, float64_elt, fortran_layout) Bigarray.Genarray.t let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = let exp_pq = 1. /. expo_pq_inv in let t = norm_pq_sq *. exp_pq in - boys_function ~maxm t - |> Array.mapi (fun m fm -> - two_over_sq_pi *. (if m mod 2 = 0 then fm else -.fm) *. - (pow exp_pq m) *. (sqrt exp_pq) - ) + let f = two_over_sq_pi *. (sqrt exp_pq) in + let result = boys_function ~maxm t in + let rec aux accu k = function + | 0 -> result.(k) <- result.(k) *. accu + | l -> + begin + result.(k) <- result.(k) *. accu; + let new_accu = -. accu *. exp_pq in + aux new_accu (k+1) (l-1) + end + in + aux f 0 maxm; + result (** Compute all the integrals of a contracted class *)