# Pragmas.
use strict;

# Modules.
use bignum;
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;
        }
}


print gamma_lanczos(25.78), "\n";
system("pause");