10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-07-04 18:35:50 +02:00

gamma function from libm

This commit is contained in:
Anthony Scemama 2018-02-01 22:39:23 +01:00
parent daf96d567b
commit d93729deff
3 changed files with 29 additions and 6 deletions

View File

@ -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)

View File

@ -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

View File

@ -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);
}