30 lines
937 B
Text
30 lines
937 B
Text
INT g = 7;
|
|
[]REAL p = []REAL(0.99999999999980993, 676.5203681218851, -1259.1392167224028,
|
|
771.32342877765313, -176.61502916214059, 12.507343278686905,
|
|
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)[@0];
|
|
|
|
PROC complex gamma = (COMPL in z)COMPL: (
|
|
# Reflection formula #
|
|
COMPL z := in z;
|
|
IF re OF z < 0.5 THEN
|
|
pi / (complex sin(pi*z)*complex gamma(1-z))
|
|
ELSE
|
|
z -:= 1;
|
|
COMPL x := p[0];
|
|
FOR i TO g+1 DO x +:= p[i]/(z+i) OD;
|
|
COMPL t := z + g + 0.5;
|
|
complex sqrt(2*pi) * t**(z+0.5) * complex exp(-t) * x
|
|
FI
|
|
);
|
|
|
|
OP ** = (COMPL z, p)COMPL: ( z=0|0|complex exp(complex ln(z)*p) );
|
|
PROC factorial = (COMPL n)COMPL: complex gamma(n+1);
|
|
|
|
FORMAT compl fmt = $g(-16, 8)"⊥"g(-10, 8)$;
|
|
|
|
test:(
|
|
printf(($q"factorial(-0.5)**2="f(compl fmt)l$, factorial(-0.5)**2));
|
|
FOR i TO 9 DO
|
|
printf(($q"factorial("d")="f(compl fmt)l$, i, factorial(i)))
|
|
OD
|
|
)
|