GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
gamln.f
Go to the documentation of this file.
1  FUNCTION gamln(Z,IERR)
2 C***BEGIN PROLOGUE GAMLN
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 GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
12 C Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES
13 C GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION
14 C G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN. THE FUNCTION WAS MADE AS
15 C PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE
16 C 10 DIGITS IN A WORD, RLN=AMAX1(-ALOG10(R1MACH(4)),0.5E-18)
17 C LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.
18 C
19 C SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100
20 C VALUES IS USED FOR SPEED OF EXECUTION.
21 C
22 C DESCRIPTION OF ARGUMENTS
23 C
24 C INPUT
25 C Z - REAL ARGUMENT, Z.GT.0.0E0
26 C
27 C OUTPUT
28 C GAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z
29 C IERR - ERROR FLAG
30 C IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
31 C IERR=1, Z.LE.0.0E0, NO COMPUTATION
32 C
33 C***REFERENCES COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
34 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
35 C***ROUTINES CALLED I1MACH,R1MACH
36 C***END PROLOGUE GAMLN
37 C
38  INTEGER i, i1m, k, mz, nz, ierr, i1mach
39  REAL cf, con, fln, fz, gln, rln, s, tlg, trm, tst, t1, wdtol, z,
40  * zdmy, zinc, zm, zmin, zp, zsq
41  REAL r1mach
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.00000000000000000e+00, 0.00000000000000000e+00,
49  5 6.93147180559945309e-01, 1.79175946922805500e+00,
50  6 3.17805383034794562e+00, 4.78749174278204599e+00,
51  7 6.57925121201010100e+00, 8.52516136106541430e+00,
52  8 1.06046029027452502e+01, 1.28018274800814696e+01,
53  9 1.51044125730755153e+01, 1.75023078458738858e+01,
54  a 1.99872144956618861e+01, 2.25521638531234229e+01,
55  b 2.51912211827386815e+01, 2.78992713838408916e+01,
56  c 3.06718601060806728e+01, 3.35050734501368889e+01,
57  d 3.63954452080330536e+01, 3.93398841871994940e+01,
58  e 4.23356164607534850e+01, 4.53801388984769080e+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.84711813518352239e+01, 5.16066755677643736e+01,
64  5 5.47847293981123192e+01, 5.80036052229805199e+01,
65  6 6.12617017610020020e+01, 6.45575386270063311e+01,
66  7 6.78897431371815350e+01, 7.12570389671680090e+01,
67  8 7.46582363488301644e+01, 7.80922235533153106e+01,
68  9 8.15579594561150372e+01, 8.50544670175815174e+01,
69  a 8.85808275421976788e+01, 9.21361756036870925e+01,
70  b 9.57196945421432025e+01, 9.93306124547874269e+01,
71  c 1.02968198614513813e+02, 1.06631760260643459e+02,
72  d 1.10320639714757395e+02, 1.14034211781461703e+02,
73  e 1.17771881399745072e+02, 1.21533081515438634e+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.25317271149356895e+02, 1.29123933639127215e+02,
79  5 1.32952575035616310e+02, 1.36802722637326368e+02,
80  6 1.40673923648234259e+02, 1.44565743946344886e+02,
81  7 1.48477766951773032e+02, 1.52409592584497358e+02,
82  8 1.56360836303078785e+02, 1.60331128216630907e+02,
83  9 1.64320112263195181e+02, 1.68327445448427652e+02,
84  a 1.72352797139162802e+02, 1.76395848406997352e+02,
85  b 1.80456291417543771e+02, 1.84533828861449491e+02,
86  c 1.88628173423671591e+02, 1.92739047287844902e+02,
87  d 1.96866181672889994e+02, 2.01009316399281527e+02,
88  e 2.05168199482641199e+02, 2.09342586752536836e+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.13532241494563261e+02, 2.17736934113954227e+02,
94  5 2.21956441819130334e+02, 2.26190548323727593e+02,
95  6 2.30439043565776952e+02, 2.34701723442818268e+02,
96  7 2.38978389561834323e+02, 2.43268849002982714e+02,
97  8 2.47572914096186884e+02, 2.51890402209723194e+02,
98  9 2.56221135550009525e+02, 2.60564940971863209e+02,
99  a 2.64921649798552801e+02, 2.69291097651019823e+02,
100  b 2.73673124285693704e+02, 2.78067573440366143e+02,
101  c 2.82474292687630396e+02, 2.86893133295426994e+02,
102  d 2.91323950094270308e+02, 2.95766601350760624e+02,
103  e 3.00220948647014132e+02, 3.04686856765668715e+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.09164193580146922e+02, 3.13652829949879062e+02,
107  3 3.18152639620209327e+02, 3.22663499126726177e+02,
108  4 3.27185287703775217e+02, 3.31717887196928473e+02,
109  5 3.36261181979198477e+02, 3.40815058870799018e+02,
110  6 3.45379407062266854e+02, 3.49954118040770237e+02,
111  7 3.54539085519440809e+02, 3.59134205369575399e+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.33333333333333333e-02, -2.77777777777777778e-03,
117  4 7.93650793650793651e-04, -5.95238095238095238e-04,
118  5 8.41750841750841751e-04, -1.91752691752691753e-03,
119  6 6.41025641025641026e-03, -2.95506535947712418e-02,
120  7 1.79644372368830573e-01, -1.39243221690590112e+00,
121  8 1.34028640441683920e+01, -1.56848284626002017e+02,
122  9 2.19310333333333333e+03, -3.61087712537249894e+04,
123  a 6.91472268851313067e+05, -1.52382215394074162e+07,
124  b 3.82900751391414141e+08, -1.08822660357843911e+10,
125  c 3.47320283765002252e+11, -1.23696021422692745e+13,
126  d 4.88788064793079335e+14, -2.13203339609193739e+16/
127 C
128 C LN(2*PI)
129  DATA con / 1.83787706640934548e+00/
130 C
131 C***FIRST EXECUTABLE STATEMENT GAMLN
132  ierr=0
133  IF (z.LE.0.0e0) GO TO 70
134  IF (z.GT.101.0e0) GO TO 10
135  nz = int(z)
136  fz = z - float(nz)
137  IF (fz.GT.0.0e0) GO TO 10
138  IF (nz.GT.100) GO TO 10
139  gamln = gln(nz)
140  RETURN
141  10 CONTINUE
142  wdtol = r1mach(4)
143  wdtol = amax1(wdtol,0.5e-18)
144  i1m = i1mach(11)
145  rln = r1mach(5)*float(i1m)
146  fln = amin1(rln,20.0e0)
147  fln = amax1(fln,3.0e0)
148  fln = fln - 3.0e0
149  zm = 1.8000e0 + 0.3875e0*fln
150  mz = int(zm) + 1
151  zmin = float(mz)
152  zdmy = z
153  zinc = 0.0e0
154  IF (z.GE.zmin) GO TO 20
155  zinc = zmin - float(nz)
156  zdmy = z + zinc
157  20 CONTINUE
158  zp = 1.0e0/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 (abs(trm).LT.tst) GO TO 40
168  s = s + trm
169  30 CONTINUE
170  40 CONTINUE
171  IF (zinc.NE.0.0e0) GO TO 50
172  tlg = alog(z)
173  gamln = z*(tlg-1.0e0) + 0.5e0*(con-tlg) + s
174  RETURN
175  50 CONTINUE
176  zp = 1.0e0
177  nz = int(zinc)
178  DO 60 i=1,nz
179  zp = zp*(z+float(i-1))
180  60 CONTINUE
181  tlg = alog(zdmy)
182  gamln = zdmy*(tlg-1.0e0) - alog(zp) + 0.5e0*(con-tlg) + s
183  RETURN
184 C
185 C
186  70 CONTINUE
187  ierr=1
188  RETURN
189  END
gamln
function gamln(Z, IERR)
Definition: gamln.f:2
i1mach
integer function i1mach(i)
Definition: i1mach.f:23
abs
static T abs(T x)
Definition: pr-output.cc:1678
r1mach
real function r1mach(i)
Definition: r1mach.f:23