Gamma function with the Lanczos approximation (Perl)

From LiteratePrograms

Other implementations: Perl | Python | Ruby

The Lanczos approximation offers a relatively simple way to calculate the gamma function of a complex argument to any desired precision. The precision is determined by a free constant g and a series whose coefficients depend on g and may be computed in advanced if a fixed precision is desired. The following implementation uses the same coefficients as the GNU Scientific Library's gamma function. The typical relative error is about 10^-15 (nearly full double precision accuracy).

Lanczos' series is only valid for numbers in the right half complex plane, so the reflection formula for the gamma function is employed for negative arguments (and also arguments in (0, 0.5), since this should increase accuracy near the singularity at zero.)

<<gamma_lanczos>>=
# Pragmas.
use strict;
# Modules.
use Math::Complex;
# Global variables.
use vars qw($g $lanczos_coef);
$g = 7;
$lanczos_coef = [
        0.99999999999980993,
        676.5203681218851,
        -1259.1392167224028,
        771.32342877765313,
        -176.61502916214059,
        12.507343278686905,
        -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7
];
sub gamma_lanczos {
        my $z = shift;
        $z = Math::Complex->make($z);
        if (Re($z) < 0.5) {
                return pi / (sin(pi * $z) * gamma(1 - $z));
        } else {
                $z -= 1;
                my $x = $lanczos_coef->[0];
                foreach my $i (1 .. $g + 2) {
                        $x += $lanczos_coef->[$i] / ($z + $i);
                }
                my $t = $z + $g + 0.5;
                return sqrt(2 * pi) * $t ** ($z + 0.5) * exp(-$t) * $x;
        }
}

Testing against a known value (Γ(10) = 9! = 362880):

<<gamma_lanczos_test.pl>>=
#!/usr/bin/env perl
gamma_lanczos
print gamma_lanczos(10), "\n";

returns 362880.000000001

Download code