GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dgamln.f
Go to the documentation of this file.
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
double precision function d1mach(i)
Definition: d1mach.f:23
double precision function dgamln(Z, IERR)
Definition: dgamln.f:2
integer function i1mach(i)
Definition: i1mach.f:23