#!/usr/bin/perl

# 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;
        }
}

for($i=0; $i<=100; $i++)
{
  $z = Math::Complex->make(1.0,$i/10.0);
  print "z: $z\n";
  $g = gamma_lanczos($z);
  print "g: $g\n";
  $lg = log($g);
  print "lg: $lg\n";
  $relg = Re($lg);
  $imlg = Im($lg);
  printf "%.1f %.8e %.8e\n",$i/10,$relg,$imlg;
}