/********************************************************************************
*  This PL/1 subroutine calculates the Gamma function using a direct numerical 
*   integration (adaptive stepwise 16-point Gauss quadrature) of Euler's integral.
*
*   Γ(z)=0∫∞dt tz-1.e-t
*
*
*  Input parameters are:
*    KV: desired relative accuracy
*    XD: integration stepwidth (recommended =5)
*    Z: argument of gamma function (complex)
*
*  Output parameter is:
*    G: value of Gamma function (complex)
*
*   Note that this program uses extended accuracy because of the possibility
*   that the integrand oscillates rapidly and the integral becomes very small
*   as a result. Normal (double or single) accuracy may not produce correct
*   results in these cases.
*   Note also that this may not be the computationally most efficient method 
*   to calculate the Gamma function. Other methods may be much faster, but 
*   are not as flexible with regard to the desired accuracy for instance.
*
*   This program has worked well and correct for me in the past (using IBM PL/1 
*   compilers), but I can't guarantee this for all compilers, especially if the 
*   code is translated into a different language.
* 
*   If there are any problems, please let me know and I mention the issue here.
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, April 2007
*****************************************************************************/

 GAMMA: PROCEDURE (KV,XD,Z,G) ; 
 DEFAULT RANGE(A:Z) COMPLEX FLOAT VALUE (DEC FLOAT(33)) ;
 DCL QG(16) INIT
 (9.894009349916499325961541735E-01, 2.71524594117540948517805725E-02,
 9.445750230732325760779884155E-01, 6.22535239386478928628438370E-02,
 8.656312023878317438804678977E-01, 9.51585116824927848099251076E-02,
 7.554044083550030338951011948E-01, 1.246289712555338720524762822E-01,
 6.178762444026437484466717640E-01, 1.495959888165767320815017305E-01,
 4.580167776572273863424194430E-01, 1.691565193950025381893120790E-01,
 2.816035507792589132304605015E-01, 1.826034150449235888667636680E-01,
 9.50125098376374401853193354E-02, 1.894506104550684962853967232E-01)
 REAL DEC FLOAT(33) ;
 DCL I REAL DEC FIXED(5,0) ;
 DCL (Z,G) COMPLEX DEC FLOAT(33) ;
 DCL (CONVGC,KV,RZ, JJ,XL,XU,XD,EIZ,PI2 INIT(3.141592653589793E00),IZ,
     QYMR,QYMI,GR,GI,CONVGCI,CONVGCR,QYR,QYI,DI,XL0) REAL DEC FLOAT(33) ;

 ON UNDERFLOW ; 
 IZ=ABS(IMAG(Z)) ; RZ=REAL(Z) ; 
 IF IZ>0 THEN EIZ=EXP(PI2/IZ); ELSE EIZ=1E10/KV ;
 QG=QG/2 ;

 XL0=MAX(RZ-1,1) ;  XU=XL0 ;
 XD=XD/(FLOOR(RZ/32)+1) ;

 G=0  ; 
 CONVGC=1 ; JJ=1 ;DI=1; 

 INT:
 DO WHILE(CONVGC>=KV) ;    /* Integration loop */
 IF DI=1 THEN DO ;   /* Integration into positive direction from real maximum  */
 XL=XU ;
 XU=XL+MIN(XL*EIZ-XL,XD) ;
 END ;
 ELSE DO ;   /* Integration into negative direction from real maximum  */
 XU=XL ;     
 XL=MAX(XU/EIZ,XU-XD) ; 
 END ;
 CALL DQG16 ;
 G=G+QY ;
 GR=REAL(G) ;  QYR=REAL(QY) ;
 GI=IMAG(G) ;  QYI=IMAG(QY) ;
 QYMR=ABS(GR)/JJ ;
 QYMI=ABS(GI)/JJ ;
 IF QYMR>0 THEN CONVGCR=ABS(QYR/QYMR);
 IF QYMI>0 THEN CONVGCI=ABS(QYI/QYMI);
 CONVGC=MAX(CONVGCR,CONVGCI) ;
 JJ=JJ+1 ;
 END ;

 IF DI=1 THEN DO ;   
 CONVGC=1 ; DI=-1 ; XL=XL0 ;
 GO TO INT ;   
 END ;


 /* Internal Subroutines */

 FCT:PROCEDURE(X) ;
 DCL X REAL DEC FLOAT(33) ;
 RETURN(X**(Z-1)*EXP(-X))   ;
 END FCT ;

 DQG16: PROCEDURE ;
 A=0.5*(XU+XL) ;
 B=XU-XL ;
 LY=0 ; 
 DO I=1 TO 15 BY 2 ;
 C=QG(I)*B ;
 LY=LY+QG(I+1)*(FCT(A+C)+FCT(A-C)) ;
 END ;
 QY=LY*B ;
 END DQG16 ;

 ENDE: END GAMMA ;