LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
Code Releases Activity

Optimized gamma_q

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-01 23:00:22 +01:00
parent 91f1dacaf4
commit 9365e69488
1 changed files with 7 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Utils/Util.ml
View File

@ -15,7 +15,9 @@ let factmax = 150
sha, Tokyo, 1991) p.227 [in Japanese] *) sha, Tokyo, 1991) p.227 [in Japanese] *)
(* Incomplete gamma function (* 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 = let rec p_gamma a x loggamma_a =
if x >= 1. +. a then 1. -. q_gamma a x loggamma_a if x >= 1. +. a then 1. -. q_gamma a x loggamma_a
else if x = 0. then 0. else if x = 0. then 0.
@ -29,8 +31,7 @@ let rec p_gamma a x loggamma_a =
in in
let r0 = exp (a *. log x -. x -. loggamma_a) /. a in let r0 = exp (a *. log x -. x -. loggamma_a) /. a in
pg_loop min_float r0 r0 1. 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 = and q_gamma a x loggamma_a =
if x < 1. +. a then 1. -. p_gamma a x loggamma_a if x < 1. +. a then 1. -. p_gamma a x loggamma_a
else else
@ -38,10 +39,11 @@ let rec p_gamma a x loggamma_a =
if k > 1000. then failwith "q_gamma could not converge." if k > 1000. then failwith "q_gamma could not converge."
else if prev = res then res else if prev = res then res
else else
let k_inv = 1. /. k in
let la, lb = 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 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 let prev, res = res, res +. w /. (la *. lb) in
qg_loop prev res la lb w (k +. 1.) qg_loop prev res la lb w (k +. 1.)
in in