.TITLE GAMMA ; .MACRO .GAMMA P,Y ;THE GAMMA FUNCTION IS DESCRIBED IN H. WAYLAND, DIFFERENTIAL EQUATIONS ;APPLIED TO SCIENCE AND ENGINEERING, VAN NOSTRAND, NEW YORK, 1957, P. 164. ;THE GAMMA FUNCTION IS A GENERALIZATION OF THE FACTORIAL FUNCTION. IT IS ;DEFINED FOR X > 0: ; GAMMA(P)=INTEGRAL FROM 0 TO INFINITY OF EXP(-V)* V**(P-1) DV ;THM 1: GAMMA(P)=(P-1)*GAMMA(P-1) ;THM 2: GAMMA(P)=(P-1)! (FACTORIAL), IF P IS POSITIVE INTEGER ;THM 3: GAMMA(0.5) = SQRT(PI) ;GAMMA FUNCTION IS TABULATED IN CD HODGMAN, RC WEAST & SM SELBY, HANDBOOK OF ;CHEMISTRY & PHYSICS, CHEMICAL RUBBER PUBL, CLEVELAND, 1958, P 278 ;THE COMPUTER ALGORITHM WAS TAKEN FROM C. HASTINGS, APPROXIMATIONS FOR ;DIGITAL COMPUTERS, PRINCETON UNIV. PRESS, PRINCETON, NJ, 1955, P. 155. ;GAMMA(1+X) ~= A[0] + A1[1]*X + A[2]*X**2 + A[3]*X**3 + A[4]*X**4 + A[5]*X**5 ;WHERE THE As ARE GIVEN AT THE END OF THE PROGRAM. OVER THE RANGE 0 <= X <= 1, ;THIS POWER SERIES GIVES GAMMA ACCURATE TO WITHIN +-0.00005. .GLOBL .GAMMA,GAMMA,.CALL,.POLY .MCALL .TLQ,.MOVF,.GAMMA GAMMA: .GAMMA @2(R5),-(SP) ;F=GAMMA(X,Y) F:=GAMMA(X,Y); MOV (SP)+,R0 MOV (SP)+,R1 RETURN A0=R0 A1=R1 Y: .BLKW 2 P: .BLKW 2 .GAMMA: STF A0,-(SP) STF A1,-(SP) LDF P,A0 CMPF #1E-30,A0 CFCC BLT S1 .TLQ <.GAMMA NOT DEFINED FOR NONPOSITIVE NUMBER> JMP .CALL S1: CMPF #33,A0 CFCC BGE S2 .TLQ <.GAMMA NUMBER GT 33> JMP .CALL S2: CMPF #1,A0 CFCC BLT S3 ;CALCULATE FOR P < 1 ADDF #1,A0 ;FROM THM 1: GAMMA(P)=GAMMA(P+1)/P CALL QGAMMA SUBF #1,A0 DIVF A0,A1 BR RET S3: CALL QGAMMA RET: STF A1,Y LDF (SP)+,A1 LDF (SP)+,A0 RETURN ;REENTRANT SUBROUTINE. A0 IS P. GAMMA(P) IS RETURNED IN A1. QGAMMA: CMPF #2,A0 ;POWER SERIES ONLY VALID FOR 1<=P<2 CFCC BGT SERIES SUBF #1,A0 CALL GAMMA MULF A0,A1 ;GAMMA(P)=(P-1)* GAMMA(P-1) ADDF #1,A0 RETURN SERIES: SUBF #1,A0 ;NOW A0=X=P-1 STF A0,.POLY-4 MOV #5,.POLY-6 MOV #A,.POLY-10 CALL .POLY LDF .POLY-14,A1 ADDF #1,A0 RETURN A: .FLT2 1, -0.5748646, 0.9512363, -0.6998588, 0.4245549, -0.1010678 .END