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
dgamma.f
Go to the documentation of this file.
1 *DECK DGAMMA
2  DOUBLE PRECISION FUNCTION dgamma (X)
3 C***BEGIN PROLOGUE DGAMMA
4 C***PURPOSE Compute the complete Gamma function.
5 C***LIBRARY SLATEC (FNLIB)
6 C***CATEGORY C7A
7 C***TYPE DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
8 C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
9 C***AUTHOR Fullerton, W., (LANL)
10 C***DESCRIPTION
11 C
12 C DGAMMA(X) calculates the double precision complete Gamma function
13 C for double precision argument X.
14 C
15 C Series for GAM on the interval 0. to 1.00000E+00
16 C with weighted error 5.79E-32
17 C log weighted error 31.24
18 C significant figures required 30.00
19 C decimal places required 32.05
20 C
21 C***REFERENCES (NONE)
22 C***ROUTINES CALLED D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
23 C***REVISION HISTORY (YYMMDD)
24 C 770601 DATE WRITTEN
25 C 890531 Changed all specific intrinsics to generic. (WRB)
26 C 890911 Removed unnecessary intrinsics. (WRB)
27 C 890911 REVISION DATE from Version 3.2
28 C 891214 Prologue converted to Version 4.0 format. (BAB)
29 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
30 C 920618 Removed space from variable name. (RWC, WRB)
31 C***END PROLOGUE DGAMMA
32  DOUBLE PRECISION x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax,
33  1 xmin, y, d9lgmc, dcsevl, d1mach
34  LOGICAL first
35 C
36  SAVE gamcs, pi, sq2pil, ngam, xmin, xmax, dxrel, first
37  DATA gamcs( 1) / +.8571195590 9893314219 2006239994 2 d-2 /
38  DATA gamcs( 2) / +.4415381324 8410067571 9131577165 2 d-2 /
39  DATA gamcs( 3) / +.5685043681 5993633786 3266458878 9 d-1 /
40  DATA gamcs( 4) / -.4219835396 4185605010 1250018662 4 d-2 /
41  DATA gamcs( 5) / +.1326808181 2124602205 8400679635 2 d-2 /
42  DATA gamcs( 6) / -.1893024529 7988804325 2394702388 6 d-3 /
43  DATA gamcs( 7) / +.3606925327 4412452565 7808221722 5 d-4 /
44  DATA gamcs( 8) / -.6056761904 4608642184 8554829036 5 d-5 /
45  DATA gamcs( 9) / +.1055829546 3022833447 3182350909 3 d-5 /
46  DATA gamcs( 10) / -.1811967365 5423840482 9185589116 6 d-6 /
47  DATA gamcs( 11) / +.3117724964 7153222777 9025459316 9 d-7 /
48  DATA gamcs( 12) / -.5354219639 0196871408 7408102434 7 d-8 /
49  DATA gamcs( 13) / +.9193275519 8595889468 8778682594 0 d-9 /
50  DATA gamcs( 14) / -.1577941280 2883397617 6742327395 3 d-9 /
51  DATA gamcs( 15) / +.2707980622 9349545432 6654043308 9 d-10 /
52  DATA gamcs( 16) / -.4646818653 8257301440 8166105893 3 d-11 /
53  DATA gamcs( 17) / +.7973350192 0074196564 6076717535 9 d-12 /
54  DATA gamcs( 18) / -.1368078209 8309160257 9949917230 9 d-12 /
55  DATA gamcs( 19) / +.2347319486 5638006572 3347177168 8 d-13 /
56  DATA gamcs( 20) / -.4027432614 9490669327 6657053469 9 d-14 /
57  DATA gamcs( 21) / +.6910051747 3721009121 3833697525 7 d-15 /
58  DATA gamcs( 22) / -.1185584500 2219929070 5238712619 2 d-15 /
59  DATA gamcs( 23) / +.2034148542 4963739552 0102605193 2 d-16 /
60  DATA gamcs( 24) / -.3490054341 7174058492 7401294910 8 d-17 /
61  DATA gamcs( 25) / +.5987993856 4853055671 3505106602 6 d-18 /
62  DATA gamcs( 26) / -.1027378057 8722280744 9006977843 1 d-18 /
63  DATA gamcs( 27) / +.1762702816 0605298249 4275966074 8 d-19 /
64  DATA gamcs( 28) / -.3024320653 7353062609 5877211204 2 d-20 /
65  DATA gamcs( 29) / +.5188914660 2183978397 1783355050 6 d-21 /
66  DATA gamcs( 30) / -.8902770842 4565766924 4925160106 6 d-22 /
67  DATA gamcs( 31) / +.1527474068 4933426022 7459689130 6 d-22 /
68  DATA gamcs( 32) / -.2620731256 1873629002 5732833279 9 d-23 /
69  DATA gamcs( 33) / +.4496464047 8305386703 3104657066 6 d-24 /
70  DATA gamcs( 34) / -.7714712731 3368779117 0390152533 3 d-25 /
71  DATA gamcs( 35) / +.1323635453 1260440364 8657271466 6 d-25 /
72  DATA gamcs( 36) / -.2270999412 9429288167 0231381333 3 d-26 /
73  DATA gamcs( 37) / +.3896418998 0039914493 2081663999 9 d-27 /
74  DATA gamcs( 38) / -.6685198115 1259533277 9212799999 9 d-28 /
75  DATA gamcs( 39) / +.1146998663 1400243843 4761386666 6 d-28 /
76  DATA gamcs( 40) / -.1967938586 3451346772 9510399999 9 d-29 /
77  DATA gamcs( 41) / +.3376448816 5853380903 3489066666 6 d-30 /
78  DATA gamcs( 42) / -.5793070335 7821357846 2549333333 3 d-31 /
79  DATA pi / 3.1415926535 8979323846 2643383279 50 d0 /
80  DATA sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
81  DATA first /.true./
82 C***FIRST EXECUTABLE STATEMENT DGAMMA
83  IF (first) THEN
84  ngam = initds(gamcs, 42, 0.1*REAL(D1MACH(3)) )
85 C
86  CALL dgamlm(xmin, xmax)
87  dxrel = sqrt(d1mach(4))
88  ENDIF
89  first = .false.
90 C
91  y = abs(x)
92  IF (y.GT.10.d0) go to 50
93 C
94 C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND. REDUCE INTERVAL AND FIND
95 C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
96 C
97  n = x
98  IF (x.LT.0.d0) n = n - 1
99  y = x - n
100  n = n - 1
101  dgamma = 0.9375d0 + dcsevl(2.d0*y-1.d0, gamcs, ngam)
102  IF (n.EQ.0) RETURN
103 C
104  IF (n.GT.0) go to 30
105 C
106 C COMPUTE GAMMA(X) FOR X .LT. 1.0
107 C
108  n = -n
109  IF (x .EQ. 0.d0) CALL xermsg('SLATEC', 'DGAMMA', 'X IS 0', 4, 2)
110  IF (x .LT. 0.0 .AND. x+n-2 .EQ. 0.d0) CALL xermsg('SLATEC',
111  + 'DGAMMA', 'X IS A NEGATIVE INTEGER', 4, 2)
112  IF (x .LT. (-0.5d0) .AND. abs((x-aint(x-0.5d0))/x) .LT. dxrel)
113  + CALL xermsg('SLATEC', 'DGAMMA',
114  + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',
115  + 1, 1)
116 C
117  DO 20 i=1,n
118  dgamma = dgamma/(x+i-1 )
119  20 CONTINUE
120  RETURN
121 C
122 C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
123 C
124  30 DO 40 i=1,n
125  dgamma = (y+i) * dgamma
126  40 CONTINUE
127  RETURN
128 C
129 C GAMMA(X) FOR ABS(X) .GT. 10.0. RECALL Y = ABS(X).
130 C
131  50 IF (x .GT. xmax) CALL xermsg('SLATEC', 'DGAMMA',
132  + 'X SO BIG GAMMA OVERFLOWS', 3, 2)
133 C
134  dgamma = 0.d0
135  IF (x .LT. xmin) CALL xermsg('SLATEC', 'DGAMMA',
136  + 'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
137  IF (x.LT.xmin) RETURN
138 C
139  dgamma = exp((y-0.5d0)*log(y) - y + sq2pil + d9lgmc(y) )
140  IF (x.GT.0.d0) RETURN
141 C
142  IF (abs((x-aint(x-0.5d0))/x) .LT. dxrel) CALL xermsg('SLATEC',
143  + 'DGAMMA',
144  + 'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
145 C
146  sinpiy = sin(pi*y)
147  IF (sinpiy .EQ. 0.d0) CALL xermsg('SLATEC', 'DGAMMA',
148  + 'X IS A NEGATIVE INTEGER', 4, 2)
149 C
150  dgamma = -pi/(y*sinpiy*dgamma)
151 C
152  RETURN
153  END
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2