diff --git a/Utils/Util.ml b/Utils/Util.ml
index 3a7e37d..38c4335 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -15,7 +15,9 @@ let factmax = 150
   sha, Tokyo, 1991) p.227 [in Japanese] *)
 
 (* Incomplete gamma function
-    1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  *)
+  p:  1 / Gamma(a) * Int_0^x   exp(-t) t^(a-1) dt  
+  q:  1 / Gamma(a) * Int_x^inf 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.
@@ -29,8 +31,7 @@ let rec p_gamma a x loggamma_a =
     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
@@ -38,10 +39,11 @@ let rec p_gamma a x loggamma_a =
         if k > 1000. then failwith "q_gamma could not converge."
         else if prev = res then res
         else
+          let k_inv = 1. /. k in
           let la, lb =
-            lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) /. k
+            lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) *. k_inv
           in
-          let w = w *. (k -. 1. -. a) /. k in
+          let w = w *. (k -. 1. -. a) *. k_inv in
           let prev, res = res, res +. w /. (la *. lb) in
           qg_loop prev res la lb w (k +. 1.)
       in