Task — Gamma function

Implement the subroutine gamma() using the Lanczos approximation method.

Useful discussion by Ryan Thompson here.

Examples

Example 1:
gamma(3) = 1.99
gamma(5) = 24
gamma(7) = 719.99

Analysis

It's just about 50 years since I completed my PhD and I don't recall using a gamma function since then. I do, though, remember that for an integer j, gamma(j) is (j-1)! But that doesn't help much.

Fortunately the kind people at Wikipedia have provided some Python code to implement the Lanczos approximation, and I was able to translate that easily into Perl, handling the complex numbers with Math::Complex. I'd never used Math::Complex before, and I must say it's very easy to use - compared, for example, to BigInt.

Now you could complain that I haven't really answered the question because I didn't compute the parameters, but at least I've produced working code, and in the real world, that's what counts.

Try it 

Script

#!/usr/bin/perl

# Peter Campbell Smith - 2022-05-30
# PWC 167 task 2

use v5.26;
use strict;
use warnings;
use utf8;
use Math::Complex;

my (@p, $EPSILON, $j);

@p = (676.5203681218851,
    -1259.1392167224028,
      771.32342877765313,
     -176.61502916214059,
       12.507343278686905,
       -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7);

$EPSILON = 1e-7;

# some test cases
for $j (3, 5, 7, 1.9, 2.1, 50, 0.6, 0.5, 0.4) {
    say qq[gamma($j) = ] . gamma($j);
}

sub drop_imag {
    
    # removes the imaginary part of arg if < 10**-7
    my $z = shift;
    if (abs(Im($z)) <= $EPSILON) {
        $z = Re($z);
    }
    return $z;
}

sub gamma {
    
    # implements Lanczos approximation method
    my ($i, $t, $x, $y, $z);
    $z = shift;
    
    # reflects value if arg < 0.5
    if (Re($z) < 0.5) {
        $y = pi / (sin(pi * $z) * gamma(1 - $z));
        
    # the basic algorithm
    } else {
        $z -= 1;
        $x = 0.99999999999980993;
        for $i (0 .. scalar(@p) - 1) {
            $x += $p[$i] / ($z + $i + 1);
        }
        $t = $z + scalar(@p) - 0.5;
        $y = sqrt(2 * pi) * $t ** ($z + 0.5) * exp(-$t) * $x;
    }
    return drop_imag($y);
}

Output

gamma(3) = 2
gamma(5) = 24
gamma(7) = 720.000000000002
gamma(1.9) = 0.961765831907388
gamma(2.1) = 1.04648584685356
gamma(50) = 6.08281864034254e+62
gamma(0.6) = 1.48919224881282
gamma(0.5) = 1.77245385090552
gamma(0.4) = 2.21815954375769