64 lines
1.6 KiB
Plaintext
64 lines
1.6 KiB
Plaintext
let e = exp 1.
|
|
let pi = 4. *. atan 1.
|
|
let sqrttwopi = sqrt (2. *. pi)
|
|
|
|
module Lanczos = struct
|
|
(* Lanczos method *)
|
|
(* Coefficients used by the GNU Scientific Library *)
|
|
let g = 7.
|
|
let c = [|0.99999999999980993; 676.5203681218851; -1259.1392167224028;
|
|
771.32342877765313; -176.61502916214059; 12.507343278686905;
|
|
-0.13857109526572012; 9.9843695780195716e-6; 1.5056327351493116e-7|]
|
|
|
|
let rec ag z d =
|
|
if d = 0 then c.(0) +. ag z 1
|
|
else if d < 8 then c.(d) /. (z +. float d) +. ag z (succ d)
|
|
else c.(d) /. (z +. float d)
|
|
|
|
let gamma z =
|
|
let z = z -. 1. in
|
|
let p = z +. g +. 0.5 in
|
|
sqrttwopi *. p ** (z +. 0.5) *. exp (-. p) *. ag z 0
|
|
end
|
|
|
|
module Stirling = struct
|
|
(* Stirling method *)
|
|
let gamma z =
|
|
sqrttwopi /. sqrt z *. (z /. e) ** z
|
|
|
|
end
|
|
|
|
module Stirling2 = struct
|
|
(* Extended Stirling method seen in Abramowitz and Stegun *)
|
|
let d = [|1./.12.; 1./.288.; -139./.51840.; -571./.2488320.|]
|
|
|
|
let rec corr z x n =
|
|
if n < Array.length d - 1 then d.(n) /. x +. corr z (x *. z) (succ n)
|
|
else d.(n) /. x
|
|
|
|
let gamma z = Stirling.gamma z *. (1. +. corr z z 0)
|
|
end
|
|
|
|
let mirror gma z =
|
|
if z > 0.5 then gma z
|
|
else pi /. sin (pi *. z) /. gma (1. -. z)
|
|
|
|
let _ =
|
|
Printf.printf "z\t\tLanczos\t\tStirling\tStirling2\n";
|
|
for i = 1 to 20 do
|
|
let z = float i /. 10. in
|
|
Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
|
|
z
|
|
(mirror Lanczos.gamma z)
|
|
(mirror Stirling.gamma z)
|
|
(mirror Stirling2.gamma z)
|
|
done;
|
|
for i = 1 to 7 do
|
|
let z = 10. *. float i in
|
|
Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
|
|
z
|
|
(Lanczos.gamma z)
|
|
(Stirling.gamma z)
|
|
(Stirling2.gamma z)
|
|
done
|