From c666878dcde773ac73a217ac9f919c25a6694a21 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 1 Feb 2018 19:06:29 +0100
Subject: [PATCH] Fast 00

---
 Utils/Constants.ml |  2 ++
 Utils/Util.ml      | 48 +++++++++++++++++++++++++++-------------------
 2 files changed, 30 insertions(+), 20 deletions(-)

diff --git a/Utils/Constants.ml b/Utils/Constants.ml
index 556587a..1f855f7 100644
--- a/Utils/Constants.ml
+++ b/Utils/Constants.ml
@@ -2,6 +2,8 @@ let cutoff = 1.e-15
 
 (** Constants *)
 let pi = acos (-1.)
+let sq_pi = sqrt pi
+let sq_pi_over_two = sq_pi *. 0.5
 let pi_inv = 1. /. pi
 let two_over_sq_pi = 2. /. (sqrt pi)
 
diff --git a/Utils/Util.ml b/Utils/Util.ml
index 4c9e76e..cecca3c 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -7,26 +7,34 @@ let factmax = 150
     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
+  match maxm with
+  | 0 ->
+    begin
+      if t = 0. then [| 1. |] else
+      let sq_t = sqrt t in
+      [| (sq_pi_over_two /. sq_t) *.  Gsl.Sf.erf sq_t |]
+    end
+  | _ ->
+    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