diff --git a/Makefile b/Makefile
index 92b2ee1..6bc9375 100644
--- a/Makefile
+++ b/Makefile
@@ -6,7 +6,7 @@ PKGS=
 OCAMLCFLAGS="-g -warn-error A"
 OCAMLOPTFLAGS="opt -O3 -nodynlink -remove-unused-arguments -rounds 16 -inline 100 -inline-max-unroll 100"
 OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -Is $(INCLUDE_DIRS) -ocamlopt $(OCAMLOPTFLAGS)
-MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) 
+MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) Utils/math_functions.c
 MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) 
 MLFILES= $(wildcard */*.ml)  $(wildcard *.ml) 
 MLIFILES=$(wildcard */*.mli) $(wildcard *.mli) 
diff --git a/Utils/Util.ml b/Utils/Util.ml
index e81b8d0..018ae15 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -1,5 +1,9 @@
-include Constants
+(** Functions from libm *)
 external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed]  [@@noalloc]
+external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed]  [@@noalloc]
+external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed]  [@@noalloc]
+
+include Constants
 
 let factmax = 150
 
@@ -13,9 +17,6 @@ let rec boys_function ~maxm t =
     begin
       if t = 0. then [| 1. |] else
       let sq_t = sqrt t in
-      (*
-      [| (sq_pi_over_two /. sq_t) *.  Gsl.Sf.erf sq_t |]
-      *)
       [| (sq_pi_over_two /. sq_t) *.  erf_float sq_t |]
     end
   | _ ->
@@ -24,7 +25,7 @@ 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 = 
-          Gsl.Sf.gamma alpha *. (  Gsl.Sf.gamma_inc_P alpha x )
+          (gamma_float alpha) *. (  Gsl.Sf.gamma_inc_P alpha x )
         in
         let t_inv = 
           1. /. t
diff --git a/Utils/math_functions.c b/Utils/math_functions.c
index 4109aab..1ee4939 100644
--- a/Utils/math_functions.c
+++ b/Utils/math_functions.c
@@ -13,3 +13,25 @@ CAMLprim double erf_float(double x)
   return erf(x);
 }
 
+
+CAMLprim value erfc_float_bytecode(value x)
+{
+  return copy_double(erfc(Double_val(x)));
+}  
+
+CAMLprim double erfc_float(double x)
+{
+  return erfc(x);
+}
+
+
+CAMLprim value gamma_float_bytecode(value x)
+{
+  return copy_double(tgamma(Double_val(x)));
+}  
+
+CAMLprim double gamma_float(double x)
+{
+  return tgamma(x);
+}
+