%PROCESS MARGINS(1,100); (NOUNDERFLOW): GAMMA: PROC(Z, MG, DG) OPTIONS (REORDER); /*----------------------------------------------------------------- Procedure GAMMA computes Gamma function ?(Z) for complex argument Z. Result ?(Z) outputs in "normalized" form (see below) without reduction mantissa to zero order. For values Z such that ABS(Z)12 (**) Stirlirg asymptotic formula is used ([1], 6.1.37, [2], table 2). For transfer from left halfplane to right one done by reflection formula ([1], 6.1.17). Values Z not satisfying to condition ((*) and (**) are tansferred to one of theese areas by addition or substraction appropiate natural number, and final result are calculated by recurrence formula [1],6.1.16. So, computation of ?(Z) for -25 <= RE(Z) <= 25, -1 <= IM(Z) <= 1 reduces to power series, and of others to Sirling formula. NOTE! In returned mantissa MG L=lg(DG) digits are lost For ABS(L) < 1 error no more few unit of 15th digit. Used procedures: TEXP. ENTRY declaration: DCL GAMMAP ENTRY(CPLX(18), CPLX(18), FLOAT(18)); Parameters: Z CPLX(18) - given complex argument; MG CPLX(18) - computed mantissa of Gamma function; DG FLOAT(18) - computed order (degree of 10) of Gamma function Result: ?(Z)=MG*(10**DG). References: 1.Handbook of Mathematical Functions. Eds M. Abramovitz & I.A. Stegun, NBS, 1964 2.J.W.WRENCH,JR. CONCERNING TWO SERIES FOR THE GAMMA FUNCTION, MATH. COMP. 1968,V 22,N 103, P.617-626. Version: 09.01.90 M.G. Belkina, M.Kh. Zimnov ---------------------------------------------------------------------*/ DCL Z CPLX float(18), /* Given argument */ MG CPLX float(18), /* Mantissa of Gamma function */ DG FLOAT(18), /* Order of Gamma function */ (K, L, N, T) INIT(0) BIN FIXED(31), (AZ, H, P, RZ, IZ) INIT(0.0q0) FLOAT(18), (S INIT(2.30258509299404568q0), /* S=LN(10) */ PI INIT(3.14159265358979324q0), LZ INIT(9.18938533204672742q-1) /*LZ=LN(SQRT(2*PI))*/) FLOAT(18), (EC, SS, Z1, ZS) INIT(0.0q0) CPLX float(18), (ABS, FLOOR, IMAG, LOG, REAL, SIGN, TRUNC) BUILTIN, SYSPRINT FILE PRINT, TEXP ENTRY(CPLX(18),FLOAT(18),CPLX(18)); /* ON ERROR BEGIN; PUT DATA; STOP; END;*/ AZ=ABS(Z); RZ=REAL(Z); IZ=IMAG(Z); T=SIGN(RZ); IF IMAG(Z)=0.0q0 & RZ<0.0q0 & RZ=TRUNC(RZ) | AZ=0.0q0 | AZ>1.0q10 THEN DO; IF IMAG(Z)=0.0q0 & RZ<0.0q0 & RZ=TRUNC(RZ) | AZ=0.0q0 THEN PUT SKIP EDIT('GAMMA: AT Z=', Z, 'GAMMA FUNCTION HAS SINGLE ', 'POLE WITH RESIDUAL EQUAL TO ', '((-1)**Z)/(-FACTORIAL(-Z)).') (A, C(e(24,15,16)), SKIP, A, A, A); ELSE PUT SKIP EDIT('GAMMA: FOR GIVEN ARGUMENT Z=', Z, ' PROCEDURE DOES NOT WORK.') (A, C(e(24,15,16)), A); PUT SKIP EDIT('GAMMA: VALUES MG=1+1I, DG=1E75 ARE RETURNED.')(A); MG=1.0q0+1.0q0I; DG=1.0q75; GO TO FIN; END; ZS=Z; /* Computation by power series for 1/?(Z) */ IF ABS(RZ)<=2.5q1 & ABS(IZ)<=1.0q0 THEN DO; N=FLOOR(ABS(RZ)); ZS=ZS-N*T; IF ZS=0.0q0 THEN DO; ZS=ZS+1.0q0; N=N-1; END; MG=ZS*(1.0q0 +ZS*( 5.772156649015329q-1 + ZS*(-6.558780715202538q-1 +ZS*(-4.200263503409523q-2 + ZS*( 1.665386113822915q-1 +ZS*(-4.219773455554434q-2 + ZS*(-9.621971527876974q-3 +ZS*( 7.218943246663099q-3 + ZS*(-1.165167591859065q-3 +ZS*(-2.152416741149510q-4 + ZS*( 1.28050282388116q-4 +ZS*(-2.0134854780788q-5 + ZS*(-1.250493482143q-6 +ZS*( 1.133027231982q-6 + ZS*(-2.05633841698q-7 +ZS*( 6.11609510448q-9 + ZS*( 5.0020076445q-9 +ZS*(-1.181274579q-9 + ZS*( 1.043426712q-10 +ZS*( 7.78226344q-12 + ZS*(-3.6968056326q-12 +ZS*( 5.1003703q-13 + ZS*(-2.058326q-14 +ZS*(-5.34812q-15 + ZS*( 1.22678q-15 +ZS*(-1.18126q-16 + ZS*( 1.187q-18 +ZS*( 1.412q-18 + ZS*(-2.299q-19 +ZS*( 1.71q-20)))))))))))))))))))))))))))))); IF ABS(MG)>1.0q-70 THEN DG=0.0q0; ELSE DO; MG=MG*1.0q10; DG=-10.0q0; END; MG=1.0q0/MG; DO K=0 TO N-1; IF RZ>=0.0q0 THEN MG=MG*(ZS+K); ELSE MG=MG/(ZS-K-1.0q0); END /*K*/; GO TO FIN; END; /* End of computation by power series for 1/?(Z) */ ELSE DO; /* Computation bty Stirling formula */ IF RZ<0.0q0 THEN DO; L=-1; /* transition to right halfplane */ ZS=1.0q0-Z; END; IF AZ<1.2q1 THEN DO; N=12-FLOOR(AZ); /* Argument reduction to form convenient for*/ ZS=ZS+N; /* Stirling formula using */ END; Z1=ZS; DO WHILE(ABS(Z1)>1.0q1); /* Z scaling */ Z1=Z1*1.0q-1; K=K+1; END; MG=LZ+(ZS-0.5q0)*(LOG(Z1)+K*S)-ZS; CALL TEXP(MG,DG,MG); Z1=1.0q0/ZS; /* Stirling formula */ MG=MG*( 1.000000000000000q00+ Z1*( 8.333333333333333q-2+ Z1*( 3.472222222222222q-3+ Z1*(-2.681327160493827q-3+ Z1*(-2.294720936213992q-4+ Z1*( 7.840392217200666q-4+ Z1*( 6.972813758365858q-5+ Z1*(-5.921664373536939q-4+ Z1*(-5.171790908260592q-5+ Z1*( 8.394987206720873q-4+ Z1*( 7.204895416020011q-5+ Z1*(-1.914438498565478q-3+ Z1*(-1.625162627839158q-4+ Z1*( 6.403362833808070q-3+ Z1*( 5.401647678926045q-4+ Z1*(-2.952788094569912q-2+ Z1*(-2.481743600264998q-3+ Z1*( 1.795401170612349q-1+ Z1* 1.505611304002642q-2)))))))))))))))))); DO K=1 TO N; MG=MG/(ZS-K); /* Recurrence formula */ END; IF L=-1 THEN DO; /* Reflection formula */ CALL TEXP(PI*Z*1.0q0I,P,EC); H=10.0q0**(-2.0q0*ABS(P)); IF P > 0.0q0 THEN SS=EC-H/EC; ELSE SS=EC*H-1.0q0/EC; MG=2.0q0I*PI/(SS*MG); DG=-ABS(P)-DG; END; END; /* End of computation bty Stirling formula */ FIN: END GAMMA; /* 1st release Jan. 26, 1987 */ TEXP: PROC(W,PW,EW) OPTIONS (REORDER); /* EXP(W)=EW*(10**PW) computation with extended accyracy */ DCL(W,EW) CPLX float(18), PW FLOAT(18), (S INIT(2.302585092994046E0), S1 INIT(2.302585E0), S2 INIT(9.29940456E-8), S3 INIT(8.40179915E-17), /* S=S1+S2+S3=LN(10).*/ PI2 INIT(6.283185307179586E0), PI21 INIT(6.283185000000000E0), PI22 INIT(3.071795860000000E-7), PI23 INIT(4.769252870000000E-16), P1) FLOAT(16); DCL (EXP,FLOOR,REAL,IMAG) BUILTIN; PW=FLOOR(REAL(W)/S); P1=FLOOR(IMAG(W)/PI2); EW=EXP(W-(S1*PW)-(S2*PW)-(S3*PW)- (PI21*P1*1E0I)-(PI22*P1*1E0I)-(PI23*P1*1E0I)); DO WHILE(ABS(EW)>=10E0); EW=EW*1E-1; PW=PW+1E0; END; DO WHILE(ABS(EW)<1E0 & ABS(EW)^=0E0); EW=EW*1E1; PW=PW-1E0; END; END TEXP;