GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
xmax
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1224
log
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
d1mach
double precision function d1mach(i)
Definition: d1mach.f:20
abs
static T abs(T x)
Definition: pr-output.cc:1696
xmin
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1233