GNU Octave  4.2.1
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
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
e
i e
Definition: data.cc:2724
d
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
Definition: lo-slatec-proto.h:109
a
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
to
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
abs
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
b
b
Definition: cellfun.cc:398
gamln
function gamln(Z, IERR)
Definition: gamln.f:1
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8