12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061 |
- #!/usr/bin/perl
- # Algorithm from Wikipedia:
- # https://en.wikipedia.org/wiki/Lanczos_approximation#Simple_implementation
- use 5.020;
- use strict;
- use warnings;
- use lib qw(../lib);
- use Math::GComplex qw(:overload real imag);
- use experimental qw(signatures lexical_subs);
- my $pi = log(-1) * -i;
- sub gamma($z) {
- my $epsilon = 0.0000001;
- my sub withinepsilon($x) {
- abs($x - abs($x)) <= $epsilon;
- }
- state $p = [
- 676.5203681218851, -1259.1392167224028,
- 771.32342877765313, -176.61502916214059,
- 12.507343278686905, -0.13857109526572012,
- 9.9843695780195716e-6, 1.5056327351493116e-7,
- ];
- my $result;
- if (real($z) < 0.5) {
- $result = ($pi / (sin($pi * $z) * gamma(1 - $z)));
- }
- else {
- $z -= 1;
- my $x = 0.99999999999980993;
- while (my ($i, $pval) = each @$p) {
- $x += $pval / ($z + $i + 1);
- }
- my $t = ($z + @$p - 0.5);
- $result = (sqrt($pi * 2) * $t**($z + 0.5) * exp(-$t) * $x);
- }
- withinepsilon(imag($result)) ? real($result) : $result;
- }
- foreach my $i (0.5, 4, 5, 6, 30, 40, 50) {
- printf("gamma(%3s) =~ %s\n", real($i), gamma($i));
- }
- __END__
- gamma(0.5) =~ 1.77245385090552
- gamma( 4) =~ 6
- gamma( 5) =~ 24
- gamma( 6) =~ 120
- gamma( 30) =~ 8.84176199373977e+30
- gamma( 40) =~ 2.0397882081197e+46
- gamma( 50) =~ 6.08281864034252e+62
|