1*DECK DGAMMA
2      DOUBLE PRECISION FUNCTION DGAMMA (X)
3C***BEGIN PROLOGUE  DGAMMA
4C***PURPOSE  Compute the complete Gamma function.
5C***LIBRARY   SLATEC (FNLIB)
6C***CATEGORY  C7A
7C***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
8C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
9C***AUTHOR  Fullerton, W., (LANL)
10C***DESCRIPTION
11C
12C DGAMMA(X) calculates the double precision complete Gamma function
13C for double precision argument X.
14C
15C Series for GAM        on the interval  0.          to  1.00000E+00
16C                                        with weighted error   5.79E-32
17C                                         log weighted error  31.24
18C                               significant figures required  30.00
19C                                    decimal places required  32.05
20C
21C***REFERENCES  (NONE)
22C***ROUTINES CALLED  D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
23C***REVISION HISTORY  (YYMMDD)
24C   770601  DATE WRITTEN
25C   890531  Changed all specific intrinsics to generic.  (WRB)
26C   890911  Removed unnecessary intrinsics.  (WRB)
27C   890911  REVISION DATE from Version 3.2
28C   891214  Prologue converted to Version 4.0 format.  (BAB)
29C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
30C   920618  Removed space from variable name.  (RWC, WRB)
31C***END PROLOGUE  DGAMMA
32      DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX,
33     1  XMIN, Y, D9LGMC, DCSEVL, D1MACH
34      LOGICAL FIRST
35C
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./
82C***FIRST EXECUTABLE STATEMENT  DGAMMA
83      IF (FIRST) THEN
84         NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACH(3)) )
85C
86         CALL DGAMLM (XMIN, XMAX)
87         DXREL = SQRT(D1MACH(4))
88      ENDIF
89      FIRST = .FALSE.
90C
91      Y = ABS(X)
92      IF (Y.GT.10.D0) GO TO 50
93C
94C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND
95C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
96C
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
103C
104      IF (N.GT.0) GO TO 30
105C
106C COMPUTE GAMMA(X) FOR X .LT. 1.0
107C
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)
116C
117      DO 20 I=1,N
118        DGAMMA = DGAMMA/(X+I-1 )
119 20   CONTINUE
120      RETURN
121C
122C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
123C
124 30   DO 40 I=1,N
125        DGAMMA = (Y+I) * DGAMMA
126 40   CONTINUE
127      RETURN
128C
129C GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
130C
131 50   IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DGAMMA',
132     +   'X SO BIG GAMMA OVERFLOWS', 3, 2)
133C
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
138C
139      DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y) )
140      IF (X.GT.0.D0) RETURN
141C
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)
145C
146      SINPIY = SIN (PI*Y)
147      IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA',
148     +   'X IS A NEGATIVE INTEGER', 4, 2)
149C
150      DGAMMA = -PI/(Y*SINPIY*DGAMMA)
151C
152      RETURN
153      END
154