removed GSL dependency

Utils/Util.ml

@@ -8,6 +8,48 @@ include Constants
 let factmax = 150
 
 
+
+
+(*reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
+  (New Algorithm handbook in C language) (Gijyutsu hyouron
+  sha, Tokyo, 1991) p.227 [in Japanese] *)
+
+(* Incomplete gamma function
+    1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  *)
+let rec p_gamma a x loggamma_a =
+  if x >= 1. +. a then 1. -. q_gamma a x loggamma_a
+  else if x = 0. then 0.
+  else
+    let rec pg_loop prev res term k =
+      if k > 1000. then failwith "p_gamma could not converge."
+      else if prev = res then res
+      else
+	let term = term *. x /. (a +. k) in
+	pg_loop res (res +. term) term (k +. 1.)
+    in
+    let r0 =  exp (a *. log x -. x -. loggamma_a) /. a in
+    pg_loop min_float r0 r0 1.
+  (* Incomplete gamma function
+     1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  *)
+  and q_gamma a x loggamma_a =
+    if x < 1. +. a then 1. -. p_gamma a x loggamma_a
+    else
+      let rec qg_loop prev res la lb w k =
+        if k > 1000. then failwith "q_gamma could not converge."
+        else if prev = res then res
+        else
+          let la, lb =
+            lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) /. k
+          in
+          let w = w *. (k -. 1. -. a) /. k in
+          let prev, res = res, res +. w /. (la *. lb) in
+          qg_loop prev res la lb w (k +. 1.)
+      in
+      let w = exp (a *. log x -. x -. loggamma_a) in
+      let lb = (1. +. x -. a) in
+      qg_loop min_float (w /. lb) 1. lb w 2.0
+
+
 (** Generalized Boys function. Uses GSL's incomplete Gamma function.
     maxm : Maximum total angular momentum
 *)
@@ -25,7 +67,10 @@ let rec boys_function ~maxm t =
         Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1))
       else
         let incomplete_gamma ~alpha x = 
-          (gamma_float alpha) *. (  Gsl.Sf.gamma_inc_P alpha x )
+(*          (gamma_float alpha) *. (  Gsl.Sf.gamma_inc_P alpha x )
+*)
+            let gf = gamma_float alpha in
+            gf *. p_gamma alpha x (log gf)
         in
         let t_inv = 
           1. /. t