71 lines
2.6 KiB
Perl
71 lines
2.6 KiB
Perl
use strict;
|
|
use warnings;
|
|
use constant pi => 4*atan2(1, 1);
|
|
use constant e => exp(1);
|
|
|
|
# Normally would be: use Math::MPFR
|
|
# but this will use it if it's installed and ignore otherwise
|
|
my $have_MPFR = eval { require Math::MPFR; Math::MPFR->import(); 1; };
|
|
|
|
sub Gamma {
|
|
my $z = shift;
|
|
my $method = shift // 'lanczos';
|
|
if ($method eq 'lanczos') {
|
|
use constant g => 9;
|
|
$z < .5 ? pi / sin(pi * $z) / Gamma(1 - $z, $method) :
|
|
sqrt(2* pi) *
|
|
($z + g - .5)**($z - .5) *
|
|
exp(-($z + g - .5)) *
|
|
do {
|
|
my @coeff = qw{
|
|
1.000000000000000174663
|
|
5716.400188274341379136
|
|
-14815.30426768413909044
|
|
14291.49277657478554025
|
|
-6348.160217641458813289
|
|
1301.608286058321874105
|
|
-108.1767053514369634679
|
|
2.605696505611755827729
|
|
-0.7423452510201416151527e-2
|
|
0.5384136432509564062961e-7
|
|
-0.4023533141268236372067e-8
|
|
};
|
|
my ($sum, $i) = (shift(@coeff), 0);
|
|
$sum += $_ / ($z + $i++) for @coeff;
|
|
$sum;
|
|
}
|
|
} elsif ($method eq 'taylor') {
|
|
$z < .5 ? Gamma($z+1, $method)/$z :
|
|
$z > 1.5 ? ($z-1)*Gamma($z-1, $method) :
|
|
do {
|
|
my $s = 0; ($s *= $z-1) += $_ for qw{
|
|
0.00000000000000000002 -0.00000000000000000023 0.00000000000000000141
|
|
0.00000000000000000119 -0.00000000000000011813 0.00000000000000122678
|
|
-0.00000000000000534812 -0.00000000000002058326 0.00000000000051003703
|
|
-0.00000000000369680562 0.00000000000778226344 0.00000000010434267117
|
|
-0.00000000118127457049 0.00000000500200764447 0.00000000611609510448
|
|
-0.00000020563384169776 0.00000113302723198170 -0.00000125049348214267
|
|
-0.00002013485478078824 0.00012805028238811619 -0.00021524167411495097
|
|
-0.00116516759185906511 0.00721894324666309954 -0.00962197152787697356
|
|
-0.04219773455554433675 0.16653861138229148950 -0.04200263503409523553
|
|
-0.65587807152025388108 0.57721566490153286061 1.00000000000000000000
|
|
}; 1/$s;
|
|
}
|
|
} elsif ($method eq 'stirling') {
|
|
no warnings qw(recursion);
|
|
$z < 100 ? Gamma($z + 1, $method)/$z :
|
|
sqrt(2*pi*$z)*($z/e + 1/(12*e*$z))**$z / $z;
|
|
} elsif ($method eq 'MPFR') {
|
|
my $result = Math::MPFR->new();
|
|
Math::MPFR::Rmpfr_gamma($result, Math::MPFR->new($z), 0);
|
|
$result;
|
|
} else { die "unknown method '$method'" }
|
|
}
|
|
|
|
for my $method (qw(MPFR lanczos taylor stirling)) {
|
|
next if $method eq 'MPFR' && !$have_MPFR;
|
|
printf "%10s: ", $method;
|
|
print join(' ', map { sprintf "%.12f", Gamma($_/3, $method) } 1 .. 10);
|
|
print "\n";
|
|
}
|