GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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