March 2014 update

This commit is contained in:
Ingy döt Net 2014-04-02 16:56:35 +00:00
parent 09687c4926
commit a25938f123
1846 changed files with 21876 additions and 5203 deletions

View file

@ -1,17 +1,22 @@
import std.stdio, std.mathspecial;
import std.stdio, std.math, std.mathspecial;
real myGamma(in real x) pure nothrow {
real taylorGamma(in real x) pure nothrow {
static immutable real[30] table = [
1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108,
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
-0.00021524167411495097, 0.00012805028238811619, -0.00002013485478078824,
-0.00000125049348214267, 0.00000113302723198170, -0.00000020563384169776,
0.00000000611609510448, 0.00000000500200764447, -0.00000000118127457049,
0.00000000010434267117, 0.00000000000778226344, -0.00000000000369680562,
0.00000000000051003703, -0.00000000000002058326, -0.00000000000000534812,
0.00000000000000122678, -0.00000000000000011813, 0.00000000000000000119,
0.00000000000000000141, -0.00000000000000000023, 0.00000000000000000002];
0x1p+0, 0x1.2788cfc6fb618f4cp-1,
-0x1.4fcf4026afa2dcecp-1, -0x1.5815e8fa27047c8cp-5,
0x1.5512320b43fbe5dep-3, -0x1.59af103c340927bep-5,
-0x1.3b4af28483e214e4p-7, 0x1.d919c527f60b195ap-8,
-0x1.317112ce3a2a7bd2p-10, -0x1.c364fe6f1563ce9cp-13,
0x1.0c8a78cd9f9d1a78p-13, -0x1.51ce8af47eabdfdcp-16,
-0x1.4fad41fc34fbb2p-20, 0x1.302509dbc0de2c82p-20,
-0x1.b9986666c225d1d4p-23, 0x1.a44b7ba22d628acap-28,
0x1.57bc3fc384333fb2p-28, -0x1.44b4cedca388f7c6p-30,
0x1.cae7675c18606c6p-34, 0x1.11d065bfaf06745ap-37,
-0x1.0423bac8ca3faaa4p-38, 0x1.1f20151323cd0392p-41,
-0x1.72cb88ea5ae6e778p-46, -0x1.815f72a05f16f348p-48,
0x1.6198491a83bccbep-50, -0x1.10613dde57a88bd6p-53,
0x1.5e3fee81de0e9c84p-60, 0x1.a0dc770fb8a499b6p-60,
-0x1.0f635344a29e9f8ep-62, 0x1.43d79a4b90ce8044p-66];
immutable real y = x - 1.0L;
real sm = table[$ - 1];
@ -20,9 +25,38 @@ real myGamma(in real x) pure nothrow {
return 1.0L / sm;
}
real lanczosGamma(real z) pure nothrow {
// Coefficients used by the GNU Scientific Library.
// http://en.wikipedia.org/wiki/Lanczos_approximation
enum g = 7;
static immutable real[9] table =
[ 0.99999_99999_99809_93,
676.52036_81218_851,
-1259.13921_67224_028,
771.32342_87776_5313,
-176.61502_91621_4059,
12.50734_32786_86905,
-0.13857_10952_65720_12,
9.98436_95780_19571_6e-6,
1.50563_27351_49311_6e-7];
// Reflection formula.
if (z < 0.5L) {
return PI / (sin(PI * z) * lanczosGamma(1 - z));
} else {
z -= 1;
real x = table[0];
foreach (immutable i; 1 .. g + 2)
x += table[i] / (z + i);
immutable real t = z + g + 0.5L;
return sqrt(2 * PI) * t ^^ (z + 0.5L) * exp(-t) * x;
}
}
void main() {
foreach (immutable i; 1 .. 11) {
immutable real x = i / 3.0L;
writefln("%f: %20.19e %20.19e", x, myGamma(x), gamma(x));
writefln("%f: %20.19e %20.19e %20.19e", x,
x.taylorGamma, x.lanczosGamma, x.gamma);
}
}

View file

@ -0,0 +1,15 @@
8 constant (gamma-shift)
: (mortici) ( f1 -- f2)
-1 s>f f+ 1 s>f
fover 271828183e-8 f* 12 s>f f* f/
fover 271828183e-8 f/ f+
fover f** fswap
628318530e-8 f* fsqrt f* \ 2*pi
;
: gamma ( f1 -- f2)
fdup f0< >r fdup f0= r> or abort" Gamma less or equal to zero"
fdup (gamma-shift) s>f f+ (mortici) fswap
1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
;

View file

@ -0,0 +1,44 @@
implement Lanczos7;
include "sys.m"; sys: Sys;
include "draw.m";
include "math.m"; math: Math;
lgamma, exp, pow, sqrt: import math;
Lanczos7: module {
init: fn(nil: ref Draw->Context, nil: list of string);
};
init(nil: ref Draw->Context, nil: list of string)
{
sys = load Sys Sys->PATH;
math = load Math Math->PATH;
# We ignore some floating point exceptions:
math->FPcontrol(0, Math->OVFL|Math->UNFL);
ns : list of real = -0.5 :: 0.1 :: 0.5 :: 1.0 :: 1.5 :: 2.0 :: 3.0 :: 10.0 :: 140.0 :: 170.0 :: nil;
sys->print("%5s %24s %24s\n", "x", "math->lgamma", "lanczos7");
while(ns != nil) {
x := hd ns;
ns = tl ns;
# math->lgamma returns a tuple.
(i, r) := lgamma(x);
g := real i * exp(r);
sys->print("%5.1f %24.16g %24.16g\n", x, g, lanczos7(x));
}
}
lanczos7(z: real): real
{
t := z + 6.5;
x := 0.99999999999980993 +
676.5203681218851/z -
1259.1392167224028/(z+1.0) +
771.32342877765313/(z+2.0) -
176.61502916214059/(z+3.0) +
12.507343278686905/(z+4.0) -
0.13857109526572012/(z+5.0) +
9.9843695780195716e-6/(z+6.0) +
1.5056327351493116e-7/(z+7.0);
return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;
}

View file

@ -1,18 +1,51 @@
sub Γ($x) {
1 / <
0.00000_00000_00000_00002 -0.00000_00000_00000_00023 0.00000_00000_00000_00141
0.00000_00000_00000_00119 -0.00000_00000_00000_11813 0.00000_00000_00001_22678
-0.00000_00000_00005_34812 -0.00000_00000_00020_58326 0.00000_00000_00510_03703
-0.00000_00000_03696_80562 0.00000_00000_07782_26344 0.00000_00001_04342_67117
-0.00000_00011_81274_57049 0.00000_00050_02007_64447 0.00000_00061_16095_10448
-0.00000_02056_33841_69776 0.00000_11330_27231_98170 -0.00000_12504_93482_14267
-0.00002_01348_54780_78824 0.00012_80502_82388_11619 -0.00021_52416_74114_95097
-0.00116_51675_91859_06511 0.00721_89432_46663_09954 -0.00962_19715_27876_97356
-0.04219_77345_55544_33675 0.16653_86113_82291_48950 -0.04200_26350_34095_23553
-0.65587_80715_20253_88108 0.57721_56649_01532_86061 1.00000_00000_00000_00000
>.reduce: * *($x-1) + *
sub Γ(\z, :$method = 'lanczos') {
given $method {
when 'lanczos' {
constant g = 9;
z < .5 ?? pi/ sin(pi * z) / Γ(1 - z, :$method) !!
sqrt(2*pi) *
(z + g - 1/2)**(z - 1/2) *
exp(-(z + g - 1/2)) *
[+] <
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
> Z* 1, map 1/(z + *), 0..*
}
when 'stirling' {
z < 10 ?? Γ(z + 1, :$method) / z !!
sqrt(2*pi*z)*(z/e + 1/(12*z*e))**z / z
}
when 'taylor' {
z > 3/2 ?? (z - 1)*Γ(z - 1, :$method) !!
z < 1/2 ?? Γ(z + 1, :$method)/z !!
1 / <
0.00000_00000_00000_00002 -0.00000_00000_00000_00023 0.00000_00000_00000_00141
0.00000_00000_00000_00119 -0.00000_00000_00000_11813 0.00000_00000_00001_22678
-0.00000_00000_00005_34812 -0.00000_00000_00020_58326 0.00000_00000_00510_03703
-0.00000_00000_03696_80562 0.00000_00000_07782_26344 0.00000_00001_04342_67117
-0.00000_00011_81274_57049 0.00000_00050_02007_64447 0.00000_00061_16095_10448
-0.00000_02056_33841_69776 0.00000_11330_27231_98170 -0.00000_12504_93482_14267
-0.00002_01348_54780_78824 0.00012_80502_82388_11619 -0.00021_52416_74114_95097
-0.00116_51675_91859_06511 0.00721_89432_46663_09954 -0.00962_19715_27876_97356
-0.04219_77345_55544_33675 0.16653_86113_82291_48950 -0.04200_26350_34095_23553
-0.65587_80715_20253_88108 0.57721_56649_01532_86061 1.00000_00000_00000_00000
>.reduce: * *(z-1) + *
}
default { die "unknow method '$method'" }
}
}
 
for 1/3, 2/3 ... 10/3 {
printf "Γ(%.5f)  %.14f\n", $_, Γ($_);
for <lanczos taylor stirling> -> $method {
# a bug in fmt forces us to round numbers before displaying them
say "{$method.fmt("%10s")}: ",
map { Γ($_, :$method).round(1e-12).fmt("%.12f") }, (1/3, 2/3 ... 10/3);
}

View file

@ -1,44 +1,61 @@
my @coeff = reverse <DATA>;
use strict;
use warnings;
use constant pi => 4*atan2(1, 1);
use constant e => exp(1);
sub Gamma {
my $x = shift;
my $s;
$s = $s*($x-1) + $_ for @coeff;
return 1 / $s;
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;
} else { die "unknown method '$method'" }
}
for (1 .. 10) {
$_ /= 3;
printf "Gamma(%.5f) ≈ %.14e\n", $_, Gamma($_);
for my $method (qw(lanczos taylor stirling)) {
printf "%10s: ", $method;
print join(' ', map { sprintf "%.12f", Gamma($_/3, $method) } 1 .. 10);
print "\n";
}
__DATA__
1.00000000000000000000
0.57721566490153286061
-0.65587807152025388108
-0.04200263503409523553
0.16653861138229148950
-0.04219773455554433675
-0.00962197152787697356
0.00721894324666309954
-0.00116516759185906511
-0.00021524167411495097
0.00012805028238811619
-0.00002013485478078824
-0.00000125049348214267
0.00000113302723198170
-0.00000020563384169776
0.00000000611609510448
0.00000000500200764447
-0.00000000118127457049
0.00000000010434267117
0.00000000000778226344
-0.00000000000369680562
0.00000000000051003703
-0.00000000000002058326
-0.00000000000000534812
0.00000000000000122678
-0.00000000000000011813
0.00000000000000000119
0.00000000000000000141
-0.00000000000000000023
0.00000000000000000002

View file

@ -1,70 +1,68 @@
/*REXX pgm calculates GAMMA using Taylor series coefficients, ≈80 digits*/
/*The GAMMA function symbol is the Greek cap letter: Γ */
numeric digits 82 /*able to use extended precision.*/
parse arg y . /*allow specification of Γ arg.*/
/*either show a range or a */
do j=word(y 1,1) to word(y 9,1) /* single gamma value(s). */
say 'gamma('j") =" gamma(j) /*compute gamma of J and display.*/
end /*j*/
exit /*stick a fork in it, we're done.*/
/*The GAMMA function symbol is the Greek capital letter: Γ */
numeric digits 84 /*able to use extended precision.*/
parse arg y z . /*allow specification of Γ arg.*/
/*either show a range or a ··· */
do j=word(y 1,1) to word(z y 9,1) /* ··· single gamma value(s). */
say 'gamma('j") =" gamma(j) /*compute gamma of J and display.*/
end /*j*/
exit /*stick a fork in it, we're done.*/
/*───────────────────────────────────GAMMA subroutine=──────────────────*/
gamma: procedure; parse arg x; xm=x-1; sum=0
#.1= 1 /*#.2 is the Euler-Mascheroni constant.*/
/*coefficients thanks to: Arne Fransén and Staffan Wrigge.*/
#.2= 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467
#.3=-0.65587807152025388107701951514539048127976638047858434729236244568387083835372210
#.4=-0.04200263503409523552900393487542981871139450040110609352206581297618009687597599
#.5= 0.16653861138229148950170079510210523571778150224717434057046890317899386605647425
#.6=-0.04219773455554433674820830128918739130165268418982248637691887327545901118558900
#.7=-0.00962197152787697356211492167234819897536294225211300210513886262731167351446074
#.8= 0.00721894324666309954239501034044657270990480088023831800109478117362259497415854
#.9=-0.00116516759185906511211397108401838866680933379538405744340750527562002584816653
#.10=-0.00021524167411495097281572996305364780647824192337833875035026748908563946371678
#.11= 0.00012805028238811618615319862632816432339489209969367721490054583804120355204347
#.12=-0.00002013485478078823865568939142102181838229483329797911526116267090822918618897
#.13=-0.00000125049348214267065734535947383309224232265562115395981534992315749121245561
#.14= 0.00000113302723198169588237412962033074494332400483862107565429550539546040842730
#.15=-0.00000020563384169776071034501541300205728365125790262933794534683172533245680371
#.16= 0.00000000611609510448141581786249868285534286727586571971232086732402927723507435
#.17= 0.00000000500200764446922293005566504805999130304461274249448171895337887737472132
#.18=-0.00000000118127457048702014458812656543650557773875950493258759096189263169643391
#.19= 0.00000000010434267116911005104915403323122501914007098231258121210871073927347588
#.20= 0.00000000000778226343990507125404993731136077722606808618139293881943550732692987
#.21=-0.00000000000369680561864220570818781587808576623657096345136099513648454655443000
#.22= 0.00000000000051003702874544759790154813228632318027268860697076321173501048565735
#.23=-0.00000000000002058326053566506783222429544855237419746091080810147188058196444349
#.24=-0.00000000000000534812253942301798237001731872793994898971547812068211168095493211
#.25= 0.00000000000000122677862823826079015889384662242242816545575045632136601135999606
#.26=-0.00000000000000011812593016974587695137645868422978312115572918048478798375081233
#.27= 0.00000000000000000118669225475160033257977724292867407108849407966482711074006109
#.28= 0.00000000000000000141238065531803178155580394756670903708635075033452562564122263
#.29=-0.00000000000000000022987456844353702065924785806336992602845059314190367014889830
#.30= 0.00000000000000000001714406321927337433383963370267257066812656062517433174649858
#.31= 0.00000000000000000000013373517304936931148647813951222680228750594717618947898583
#.32=-0.00000000000000000000020542335517666727893250253513557337960820379352387364127301
#.33= 0.00000000000000000000002736030048607999844831509904330982014865311695836363370165
#.34=-0.00000000000000000000000173235644591051663905742845156477979906974910879499841377
#.35=-0.00000000000000000000000002360619024499287287343450735427531007926413552145370486
#.36= 0.00000000000000000000000001864982941717294430718413161878666898945868429073668232
#.37=-0.00000000000000000000000000221809562420719720439971691362686037973177950067567580
#.38= 0.00000000000000000000000000012977819749479936688244144863305941656194998646391332
#.39= 0.00000000000000000000000000000118069747496652840622274541550997151855968463784158
#.40=-0.00000000000000000000000000000112458434927708809029365467426143951211941179558301
#.41= 0.00000000000000000000000000000012770851751408662039902066777511246477487720656005
#.42=-0.00000000000000000000000000000000739145116961514082346128933010855282371056899245
#.43= 0.00000000000000000000000000000000001134750257554215760954165259469306393008612196
#.44= 0.00000000000000000000000000000000004639134641058722029944804907952228463057968680
#.45=-0.00000000000000000000000000000000000534733681843919887507741819670989332090488591
#.46= 0.00000000000000000000000000000000000032079959236133526228612372790827943910901464
#.47=-0.00000000000000000000000000000000000000444582973655075688210159035212464363740144
#.48=-0.00000000000000000000000000000000000000131117451888198871290105849438992219023663
#.49= 0.00000000000000000000000000000000000000016470333525438138868182593279063941453996
#.50=-0.00000000000000000000000000000000000000001056233178503581218600561071538285049997
#.51= 0.00000000000000000000000000000000000000000026784429826430494783549630718908519485
#.52= 0.00000000000000000000000000000000000000000002424715494851782689673032938370921241
/*coefficients thanks to: Arne Fransén and Staffan Wrigge.*/
#.1 = 1 /* [↓] #.2 is the Euler-Mascheroni constant.*/
#.2 = 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467
#.3 = -0.65587807152025388107701951514539048127976638047858434729236244568387083835372210
#.4 = -0.04200263503409523552900393487542981871139450040110609352206581297618009687597599
#.5 = 0.16653861138229148950170079510210523571778150224717434057046890317899386605647425
#.6 = -0.04219773455554433674820830128918739130165268418982248637691887327545901118558900
#.7 = -0.00962197152787697356211492167234819897536294225211300210513886262731167351446074
#.8 = 0.00721894324666309954239501034044657270990480088023831800109478117362259497415854
#.9 = -0.00116516759185906511211397108401838866680933379538405744340750527562002584816653
#.10 = -0.00021524167411495097281572996305364780647824192337833875035026748908563946371678
#.11 = 0.00012805028238811618615319862632816432339489209969367721490054583804120355204347
#.12 = -0.00002013485478078823865568939142102181838229483329797911526116267090822918618897
#.13 = -0.00000125049348214267065734535947383309224232265562115395981534992315749121245561
#.14 = 0.00000113302723198169588237412962033074494332400483862107565429550539546040842730
#.15 = -0.00000020563384169776071034501541300205728365125790262933794534683172533245680371
#.16 = 0.00000000611609510448141581786249868285534286727586571971232086732402927723507435
#.17 = 0.00000000500200764446922293005566504805999130304461274249448171895337887737472132
#.18 = -0.00000000118127457048702014458812656543650557773875950493258759096189263169643391
#.19 = 0.00000000010434267116911005104915403323122501914007098231258121210871073927347588
#.20 = 0.00000000000778226343990507125404993731136077722606808618139293881943550732692987
#.21 = -0.00000000000369680561864220570818781587808576623657096345136099513648454655443000
#.22 = 0.00000000000051003702874544759790154813228632318027268860697076321173501048565735
#.23 = -0.00000000000002058326053566506783222429544855237419746091080810147188058196444349
#.24 = -0.00000000000000534812253942301798237001731872793994898971547812068211168095493211
#.25 = 0.00000000000000122677862823826079015889384662242242816545575045632136601135999606
#.26 = -0.00000000000000011812593016974587695137645868422978312115572918048478798375081233
#.27 = 0.00000000000000000118669225475160033257977724292867407108849407966482711074006109
#.28 = 0.00000000000000000141238065531803178155580394756670903708635075033452562564122263
#.29 = -0.00000000000000000022987456844353702065924785806336992602845059314190367014889830
#.30 = 0.00000000000000000001714406321927337433383963370267257066812656062517433174649858
#.31 = 0.00000000000000000000013373517304936931148647813951222680228750594717618947898583
#.32 = -0.00000000000000000000020542335517666727893250253513557337960820379352387364127301
#.33 = 0.00000000000000000000002736030048607999844831509904330982014865311695836363370165
#.34 = -0.00000000000000000000000173235644591051663905742845156477979906974910879499841377
#.35 = -0.00000000000000000000000002360619024499287287343450735427531007926413552145370486
#.36 = 0.00000000000000000000000001864982941717294430718413161878666898945868429073668232
#.37 = -0.00000000000000000000000000221809562420719720439971691362686037973177950067567580
#.38 = 0.00000000000000000000000000012977819749479936688244144863305941656194998646391332
#.39 = 0.00000000000000000000000000000118069747496652840622274541550997151855968463784158
#.40 = -0.00000000000000000000000000000112458434927708809029365467426143951211941179558301
#.41 = 0.00000000000000000000000000000012770851751408662039902066777511246477487720656005
#.42 = -0.00000000000000000000000000000000739145116961514082346128933010855282371056899245
#.43 = 0.00000000000000000000000000000000001134750257554215760954165259469306393008612196
#.44 = 0.00000000000000000000000000000000004639134641058722029944804907952228463057968680
#.45 = -0.00000000000000000000000000000000000534733681843919887507741819670989332090488591
#.46 = 0.00000000000000000000000000000000000032079959236133526228612372790827943910901464
#.47 = -0.00000000000000000000000000000000000000444582973655075688210159035212464363740144
#.48 = -0.00000000000000000000000000000000000000131117451888198871290105849438992219023663
#.49 = 0.00000000000000000000000000000000000000016470333525438138868182593279063941453996
#.50 = -0.00000000000000000000000000000000000000001056233178503581218600561071538285049997
#.51 = 0.00000000000000000000000000000000000000000026784429826430494783549630718908519485
#.52 = 0.00000000000000000000000000000000000000000002424715494851782689673032938370921241
do j=52 by -1 to 1
sum = sum * xm + #.j
end /*j*/
end /*j*/
return 1/sum