1*DECK GAMLN 2 REAL FUNCTION GAMLN (Z, IERR) 3C***BEGIN PROLOGUE GAMLN 4C***SUBSIDIARY 5C***PURPOSE Compute the logarithm of the Gamma function 6C***LIBRARY SLATEC 7C***CATEGORY C7A 8C***TYPE SINGLE PRECISION (GAMLN-S, DGAMLN-D) 9C***KEYWORDS LOGARITHM OF GAMMA FUNCTION 10C***AUTHOR Amos, D. E., (SNL) 11C***DESCRIPTION 12C 13C GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR 14C Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES 15C GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION 16C G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN. THE FUNCTION WAS MADE AS 17C PORTABLE AS POSSIBLE BY COMPUTING ZMIN FROM THE NUMBER OF BASE 18C 10 DIGITS IN A WORD, RLN=MAX(-ALOG10(R1MACH(4)),0.5E-18) 19C LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY. 20C 21C SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100 22C VALUES IS USED FOR SPEED OF EXECUTION. 23C 24C DESCRIPTION OF ARGUMENTS 25C 26C INPUT 27C Z - REAL ARGUMENT, Z.GT.0.0E0 28C 29C OUTPUT 30C GAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z 31C IERR - ERROR FLAG 32C IERR=0, NORMAL RETURN, COMPUTATION COMPLETED 33C IERR=1, Z.LE.0.0E0, NO COMPUTATION 34C 35C***REFERENCES COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 36C BY D. E. AMOS, SAND83-0083, MAY, 1983. 37C***ROUTINES CALLED I1MACH, R1MACH 38C***REVISION HISTORY (YYMMDD) 39C 830501 DATE WRITTEN 40C 830501 REVISION DATE from Version 3.2 41C 910415 Prologue converted to Version 4.0 format. (BAB) 42C 920128 Category corrected. (WRB) 43C 921215 GAMLN defined for Z negative. (WRB) 44C***END PROLOGUE GAMLN 45C 46 INTEGER I, I1M, K, MZ, NZ, IERR, I1MACH 47 REAL CF, CON, FLN, FZ, GLN, RLN, S, TLG, TRM, TST, T1, WDTOL, Z, 48 * ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ 49 REAL R1MACH 50 DIMENSION CF(22), GLN(100) 51C LNGAMMA(N), N=1,100 52 DATA GLN(1), GLN(2), GLN(3), GLN(4), GLN(5), GLN(6), GLN(7), 53 1 GLN(8), GLN(9), GLN(10), GLN(11), GLN(12), GLN(13), GLN(14), 54 2 GLN(15), GLN(16), GLN(17), GLN(18), GLN(19), GLN(20), 55 3 GLN(21), GLN(22)/ 56 4 0.00000000000000000E+00, 0.00000000000000000E+00, 57 5 6.93147180559945309E-01, 1.79175946922805500E+00, 58 6 3.17805383034794562E+00, 4.78749174278204599E+00, 59 7 6.57925121201010100E+00, 8.52516136106541430E+00, 60 8 1.06046029027452502E+01, 1.28018274800814696E+01, 61 9 1.51044125730755153E+01, 1.75023078458738858E+01, 62 A 1.99872144956618861E+01, 2.25521638531234229E+01, 63 B 2.51912211827386815E+01, 2.78992713838408916E+01, 64 C 3.06718601060806728E+01, 3.35050734501368889E+01, 65 D 3.63954452080330536E+01, 3.93398841871994940E+01, 66 E 4.23356164607534850E+01, 4.53801388984769080E+01/ 67 DATA GLN(23), GLN(24), GLN(25), GLN(26), GLN(27), GLN(28), 68 1 GLN(29), GLN(30), GLN(31), GLN(32), GLN(33), GLN(34), 69 2 GLN(35), GLN(36), GLN(37), GLN(38), GLN(39), GLN(40), 70 3 GLN(41), GLN(42), GLN(43), GLN(44)/ 71 4 4.84711813518352239E+01, 5.16066755677643736E+01, 72 5 5.47847293981123192E+01, 5.80036052229805199E+01, 73 6 6.12617017610020020E+01, 6.45575386270063311E+01, 74 7 6.78897431371815350E+01, 7.12570389671680090E+01, 75 8 7.46582363488301644E+01, 7.80922235533153106E+01, 76 9 8.15579594561150372E+01, 8.50544670175815174E+01, 77 A 8.85808275421976788E+01, 9.21361756036870925E+01, 78 B 9.57196945421432025E+01, 9.93306124547874269E+01, 79 C 1.02968198614513813E+02, 1.06631760260643459E+02, 80 D 1.10320639714757395E+02, 1.14034211781461703E+02, 81 E 1.17771881399745072E+02, 1.21533081515438634E+02/ 82 DATA GLN(45), GLN(46), GLN(47), GLN(48), GLN(49), GLN(50), 83 1 GLN(51), GLN(52), GLN(53), GLN(54), GLN(55), GLN(56), 84 2 GLN(57), GLN(58), GLN(59), GLN(60), GLN(61), GLN(62), 85 3 GLN(63), GLN(64), GLN(65), GLN(66)/ 86 4 1.25317271149356895E+02, 1.29123933639127215E+02, 87 5 1.32952575035616310E+02, 1.36802722637326368E+02, 88 6 1.40673923648234259E+02, 1.44565743946344886E+02, 89 7 1.48477766951773032E+02, 1.52409592584497358E+02, 90 8 1.56360836303078785E+02, 1.60331128216630907E+02, 91 9 1.64320112263195181E+02, 1.68327445448427652E+02, 92 A 1.72352797139162802E+02, 1.76395848406997352E+02, 93 B 1.80456291417543771E+02, 1.84533828861449491E+02, 94 C 1.88628173423671591E+02, 1.92739047287844902E+02, 95 D 1.96866181672889994E+02, 2.01009316399281527E+02, 96 E 2.05168199482641199E+02, 2.09342586752536836E+02/ 97 DATA GLN(67), GLN(68), GLN(69), GLN(70), GLN(71), GLN(72), 98 1 GLN(73), GLN(74), GLN(75), GLN(76), GLN(77), GLN(78), 99 2 GLN(79), GLN(80), GLN(81), GLN(82), GLN(83), GLN(84), 100 3 GLN(85), GLN(86), GLN(87), GLN(88)/ 101 4 2.13532241494563261E+02, 2.17736934113954227E+02, 102 5 2.21956441819130334E+02, 2.26190548323727593E+02, 103 6 2.30439043565776952E+02, 2.34701723442818268E+02, 104 7 2.38978389561834323E+02, 2.43268849002982714E+02, 105 8 2.47572914096186884E+02, 2.51890402209723194E+02, 106 9 2.56221135550009525E+02, 2.60564940971863209E+02, 107 A 2.64921649798552801E+02, 2.69291097651019823E+02, 108 B 2.73673124285693704E+02, 2.78067573440366143E+02, 109 C 2.82474292687630396E+02, 2.86893133295426994E+02, 110 D 2.91323950094270308E+02, 2.95766601350760624E+02, 111 E 3.00220948647014132E+02, 3.04686856765668715E+02/ 112 DATA GLN(89), GLN(90), GLN(91), GLN(92), GLN(93), GLN(94), 113 1 GLN(95), GLN(96), GLN(97), GLN(98), GLN(99), GLN(100)/ 114 2 3.09164193580146922E+02, 3.13652829949879062E+02, 115 3 3.18152639620209327E+02, 3.22663499126726177E+02, 116 4 3.27185287703775217E+02, 3.31717887196928473E+02, 117 5 3.36261181979198477E+02, 3.40815058870799018E+02, 118 6 3.45379407062266854E+02, 3.49954118040770237E+02, 119 7 3.54539085519440809E+02, 3.59134205369575399E+02/ 120C COEFFICIENTS OF ASYMPTOTIC EXPANSION 121 DATA CF(1), CF(2), CF(3), CF(4), CF(5), CF(6), CF(7), CF(8), 122 1 CF(9), CF(10), CF(11), CF(12), CF(13), CF(14), CF(15), 123 2 CF(16), CF(17), CF(18), CF(19), CF(20), CF(21), CF(22)/ 124 3 8.33333333333333333E-02, -2.77777777777777778E-03, 125 4 7.93650793650793651E-04, -5.95238095238095238E-04, 126 5 8.41750841750841751E-04, -1.91752691752691753E-03, 127 6 6.41025641025641026E-03, -2.95506535947712418E-02, 128 7 1.79644372368830573E-01, -1.39243221690590112E+00, 129 8 1.34028640441683920E+01, -1.56848284626002017E+02, 130 9 2.19310333333333333E+03, -3.61087712537249894E+04, 131 A 6.91472268851313067E+05, -1.52382215394074162E+07, 132 B 3.82900751391414141E+08, -1.08822660357843911E+10, 133 C 3.47320283765002252E+11, -1.23696021422692745E+13, 134 D 4.88788064793079335E+14, -2.13203339609193739E+16/ 135C 136C LN(2*PI) 137 DATA CON / 1.83787706640934548E+00/ 138C 139C***FIRST EXECUTABLE STATEMENT GAMLN 140 IERR=0 141 IF (Z.LE.0.0E0) GO TO 70 142 IF (Z.GT.101.0E0) GO TO 10 143 NZ = Z 144 FZ = Z - NZ 145 IF (FZ.GT.0.0E0) GO TO 10 146 IF (NZ.GT.100) GO TO 10 147 GAMLN = GLN(NZ) 148 RETURN 149 10 CONTINUE 150 WDTOL = R1MACH(4) 151 WDTOL = MAX(WDTOL,0.5E-18) 152 I1M = I1MACH(11) 153 RLN = R1MACH(5)*I1M 154 FLN = MIN(RLN,20.0E0) 155 FLN = MAX(FLN,3.0E0) 156 FLN = FLN - 3.0E0 157 ZM = 1.8000E0 + 0.3875E0*FLN 158 MZ = ZM + 1 159 ZMIN = MZ 160 ZDMY = Z 161 ZINC = 0.0E0 162 IF (Z.GE.ZMIN) GO TO 20 163 ZINC = ZMIN - NZ 164 ZDMY = Z + ZINC 165 20 CONTINUE 166 ZP = 1.0E0/ZDMY 167 T1 = CF(1)*ZP 168 S = T1 169 IF (ZP.LT.WDTOL) GO TO 40 170 ZSQ = ZP*ZP 171 TST = T1*WDTOL 172 DO 30 K=2,22 173 ZP = ZP*ZSQ 174 TRM = CF(K)*ZP 175 IF (ABS(TRM).LT.TST) GO TO 40 176 S = S + TRM 177 30 CONTINUE 178 40 CONTINUE 179 IF (ZINC.NE.0.0E0) GO TO 50 180 TLG = ALOG(Z) 181 GAMLN = Z*(TLG-1.0E0) + 0.5E0*(CON-TLG) + S 182 RETURN 183 50 CONTINUE 184 ZP = 1.0E0 185 NZ = ZINC 186 DO 60 I=1,NZ 187 ZP = ZP*(Z+(I-1)) 188 60 CONTINUE 189 TLG = ALOG(ZDMY) 190 GAMLN = ZDMY*(TLG-1.0E0) - ALOG(ZP) + 0.5E0*(CON-TLG) + S 191 RETURN 192C 193C 194 70 CONTINUE 195 GAMLN = R1MACH(2) 196 IERR=1 197 RETURN 198 END 199