Round primes

and gammas

Weekly challenge 167 — 30 May 2022

Week 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

The content of this website which has been created by

Peter Campbell Smith is hereby placed in the public domain

Peter Campbell Smith is hereby placed in the public domain