LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2025-01-03 01:55:40 +01:00
Code Releases Activity

Fixed overflow in Boys

Browse Source
This commit is contained in:
Anthony Scemama 2020-03-27 17:58:50 +01:00
parent 4269327046
commit e7aa4c6ac7
2 changed files with 11 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Notebooks/SpVec.ipynb
View File

@ -1746,7 +1746,7 @@
"metadata": { "metadata": {
"celltoolbar": "Raw Cell Format", "celltoolbar": "Raw Cell Format",
"kernelspec": { "kernelspec": {
"display_name": "OCaml 4.07.1", "display_name": "OCaml default",
"language": "OCaml", "language": "OCaml",
"name": "ocaml-jupyter" "name": "ocaml-jupyter"
}, },

20
Utils/Util.ml
View File

@ -209,29 +209,29 @@ let boys_function ~maxm t =
let result = let result =
Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1))
in in
(* let power_t_inv = (maxm+maxm+1) in
assert (abs_float t > 1.e-10); try
*)
if t <> 0. then
begin
let fmax = let fmax =
let t_inv = sqrt (1. /. t) in let t_inv = sqrt (1. /. t) in
let n = float_of_int maxm in let n = float_of_int maxm in
let dm = 0.5 +. n in let dm = 0.5 +. n in
let f = (pow t_inv (maxm+maxm+1) ) in let f = (pow t_inv power_t_inv ) in
match classify_float f with match classify_float f with
| FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f | FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f
| FP_zero | FP_zero
| FP_subnormal -> 0. | FP_subnormal -> 0.
| _ -> invalid_arg "zero_m overflow" | _ ->
Printf.eprintf "t=%e\n" t;
Printf.eprintf "maxm=%d\n" maxm;
raise Exit
in in
let emt = exp (-. t) in let emt = exp (-. t) in
result.(maxm) <- fmax; result.(maxm) <- fmax;
for n=maxm-1 downto 0 do for n=maxm-1 downto 0 do
result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n) result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n)
done done;
end; result
result with Exit -> result
end end