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
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:1
i e
Definition: data.cc:2724
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
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
double precision function dgamln(Z, IERR)
Definition: dgamln.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
b
Definition: cellfun.cc:398