Round primes

and gammas

Weekly challenge 167 — 30 May 2022

;

Task 2

Implement the subroutine `gamma()`

using the
Lanczos approximation method.

Useful discussion by Ryan Thompson here.

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

Cornelius Lanczos

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.

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

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

Any content of this website which has been created by Peter Campbell Smith is in the public domain