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