132 lines
4.2 KiB
Plaintext
132 lines
4.2 KiB
Plaintext
# Coefficients used by the GNU Scientific Library #
|
|
[]LONG REAL p = ( LONG 0.99999 99999 99809 93,
|
|
LONG 676.52036 81218 851,
|
|
-LONG 1259.13921 67224 028,
|
|
LONG 771.32342 87776 5313,
|
|
-LONG 176.61502 91621 4059,
|
|
LONG 12.50734 32786 86905,
|
|
-LONG 0.13857 10952 65720 12,
|
|
LONG 9.98436 95780 19571 6e-6,
|
|
LONG 1.50563 27351 49311 6e-7);
|
|
|
|
PROC lanczos gamma = (LONG REAL in z)LONG REAL: (
|
|
# Reflection formula #
|
|
LONG REAL z := in z;
|
|
IF z < LONG 0.5 THEN
|
|
long pi / (long sin(long pi*z)*lanczos gamma(1-z))
|
|
ELSE
|
|
z -:= 1;
|
|
LONG REAL x := p[1];
|
|
FOR i TO UPB p - 1 DO x +:= p[i+1]/(z+i) OD;
|
|
LONG REAL t = z + UPB p - LONG 1.5;
|
|
long sqrt(2*long pi) * t**(z+LONG 0.5) * long exp(-t) * x
|
|
FI
|
|
);
|
|
|
|
PROC taylor gamma = (LONG REAL x)LONG REAL:
|
|
BEGIN # good for values between 0 and 1 #
|
|
[]LONG REAL a =
|
|
( LONG 1.00000 00000 00000 00000,
|
|
LONG 0.57721 56649 01532 86061,
|
|
-LONG 0.65587 80715 20253 88108,
|
|
-LONG 0.04200 26350 34095 23553,
|
|
LONG 0.16653 86113 82291 48950,
|
|
-LONG 0.04219 77345 55544 33675,
|
|
-LONG 0.00962 19715 27876 97356,
|
|
LONG 0.00721 89432 46663 09954,
|
|
-LONG 0.00116 51675 91859 06511,
|
|
-LONG 0.00021 52416 74114 95097,
|
|
LONG 0.00012 80502 82388 11619,
|
|
-LONG 0.00002 01348 54780 78824,
|
|
-LONG 0.00000 12504 93482 14267,
|
|
LONG 0.00000 11330 27231 98170,
|
|
-LONG 0.00000 02056 33841 69776,
|
|
LONG 0.00000 00061 16095 10448,
|
|
LONG 0.00000 00050 02007 64447,
|
|
-LONG 0.00000 00011 81274 57049,
|
|
LONG 0.00000 00001 04342 67117,
|
|
LONG 0.00000 00000 07782 26344,
|
|
-LONG 0.00000 00000 03696 80562,
|
|
LONG 0.00000 00000 00510 03703,
|
|
-LONG 0.00000 00000 00020 58326,
|
|
-LONG 0.00000 00000 00005 34812,
|
|
LONG 0.00000 00000 00001 22678,
|
|
-LONG 0.00000 00000 00000 11813,
|
|
LONG 0.00000 00000 00000 00119,
|
|
LONG 0.00000 00000 00000 00141,
|
|
-LONG 0.00000 00000 00000 00023,
|
|
LONG 0.00000 00000 00000 00002
|
|
);
|
|
LONG REAL y = x - 1;
|
|
LONG REAL sum := a [UPB a];
|
|
FOR n FROM UPB a - 1 DOWNTO LWB a DO
|
|
sum := sum * y + a [n]
|
|
OD;
|
|
1/sum
|
|
END # taylor gamma #;
|
|
|
|
LONG REAL long e = long exp(1);
|
|
|
|
PROC sterling gamma = (LONG REAL n)LONG REAL:
|
|
( # improves for values much greater then 1 #
|
|
long sqrt(2*long pi/n)*(n/long e)**n
|
|
);
|
|
|
|
PROC factorial = (LONG INT n)LONG REAL:
|
|
(
|
|
IF n=0 OR n=1 THEN 1
|
|
ELIF n=2 THEN 2
|
|
ELSE n*factorial(n-1) FI
|
|
);
|
|
|
|
REF[]LONG REAL fm := NIL;
|
|
|
|
PROC sponge gamma = (LONG REAL x)LONG REAL:
|
|
(
|
|
INT a = 12; # alter to get required precision #
|
|
REF []LONG REAL fm := NIL;
|
|
LONG REAL res;
|
|
|
|
IF fm :=: REF[]LONG REAL(NIL) THEN
|
|
fm := HEAP[0:a-1]LONG REAL;
|
|
fm[0] := long sqrt(LONG 2*long pi);
|
|
FOR k TO a-1 DO
|
|
fm[k] := (((k-1) MOD 2=0) | 1 | -1) * long exp(a-k) *
|
|
(a-k) **(k-LONG 0.5) / factorial(k-1)
|
|
OD
|
|
FI;
|
|
res := fm[0];
|
|
FOR k TO a-1 DO
|
|
res +:= fm[k] / ( x + k )
|
|
OD;
|
|
res *:= long exp(-(x+a)) * (x+a)**(x + LONG 0.5);
|
|
res/x
|
|
);
|
|
|
|
FORMAT real fmt = $g(-real width, real width - 2)$;
|
|
FORMAT long real fmt16 = $g(-17, 17 - 2)$; # accurate to about 16 decimal places #
|
|
|
|
[]STRING methods = ("Genie", "Lanczos", "Sponge", "Taylor","Stirling");
|
|
|
|
printf(($11xg12xg12xg13xg13xgl$,methods));
|
|
|
|
FORMAT sample fmt = $"gamma("g(-3,1)")="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$;
|
|
FORMAT sqr sample fmt = $"gamma("g(-3,1)")**2="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$;
|
|
FORMAT sample exp fmt = $"gamma("g(-3)")="g(-15,11,0)n(UPB methods-1)(","g(-18,14,0))l$;
|
|
|
|
PROC sample = (LONG REAL x)[]LONG REAL:
|
|
(gamma(SHORTEN x), lanczos gamma(x), sponge gamma(x), taylor gamma(x), sterling gamma(x));
|
|
|
|
FOR i FROM 1 TO 20 DO
|
|
LONG REAL x = i / LONG 10;
|
|
printf((sample fmt, x, sample(x)));
|
|
IF i = 5 THEN # insert special case of a half #
|
|
printf((sqr sample fmt,
|
|
x, gamma(SHORTEN x)**2, lanczos gamma(x)**2, sponge gamma(x)**2,
|
|
taylor gamma(x)**2, sterling gamma(x)**2))
|
|
FI
|
|
OD;
|
|
FOR x FROM 10 BY 10 TO 70 DO
|
|
printf((sample exp fmt, x, sample(x)))
|
|
OD
|