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