22 lines
867 B
Forth
22 lines
867 B
Forth
open System.Numerics
|
|
open System
|
|
|
|
let rec gamma (z: Complex) =
|
|
let mutable z = z
|
|
let lanczosCoefficients = [| 676.520368121885; -1259.1392167224; 771.323428777653; -176.615029162141; 12.5073432786869; -0.13857109526572; 9.98436957801957E-06; 1.50563273514931E-07 |]
|
|
|
|
if z.Real < 0.5 then
|
|
Math.PI / (sin (Math.PI * z) * gamma (1.0 - z))
|
|
else
|
|
let mutable x = Complex.One
|
|
z <- z - 1.0
|
|
|
|
for i = 0 to lanczosCoefficients.Length - 1 do
|
|
x <- x + lanczosCoefficients.[i] / (z + Complex(i, 0) + 1.0)
|
|
|
|
let t = z + float lanczosCoefficients.Length - 0.5
|
|
sqrt (2.0 * Math.PI) * (t ** (z + 0.5)) * exp (-t) * x
|
|
|
|
Seq.iter (fun i -> printfn "Gamma(%f) = %A" i (gamma (Complex(i, 0)))) [ 0 .. 100 ]
|
|
Seq.iter2 (fun i j -> printfn "Gamma(%f + i%f) = %A" i j (gamma (Complex(i, j)))) [ 0 .. 100 ] [ 0 .. 100 ]
|