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.

I loved the Peter's choice of regex and substr() version. I simply love the story telling style. Keep it up great work.

This review may cover either or both challenges for this week.

I am sorry that the 'Try it' feature is currently working very slowly or not at all owing to some issue with my web hosting provider.

Try it 

Try running the script with any input:

Your $x =

example: 3.14159

⤡

#!/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);
}

18 lines of code

Output from script

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