1 /* dgamma.f -- translated by f2c (version 20041007).
2    You must link the resulting object file with libf2c:
3         on Microsoft Windows system, link with libf2c.lib;
4         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5         or, if you install libf2c.a in a standard place, with -lf2c -lm
6         -- in that order, at the end of the command line, as in
7                 cc *.o -lf2c -lm
8         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9 
10                 http://www.netlib.org/f2c/libf2c.zip
11 */
12 
13 /** This routine has been editted to be thread safe **/
14 
15 #define V3P_NETLIB_SRC
16 #include "v3p_netlib.h"
17 
18 /* Table of constant values */
19 
20 static integer c__3 = 3;
21 static integer c__42 = 42;
22 static integer c__4 = 4;
23 static integer c__2 = 2;
24 static integer c__1 = 1;
25 
26 /* DECK DGAMMA */
dgamma_(doublereal * x)27 doublereal dgamma_(doublereal *x)
28 {
29     /* Initialized data */
30 
31     static doublereal gamcs[42] = { .008571195590989331421920062399942,
32             .004415381324841006757191315771652,
33             .05685043681599363378632664588789,
34             -.004219835396418560501012500186624,
35             .001326808181212460220584006796352,
36             -1.893024529798880432523947023886e-4,
37             3.606925327441245256578082217225e-5,
38             -6.056761904460864218485548290365e-6,
39             1.055829546302283344731823509093e-6,
40             -1.811967365542384048291855891166e-7,
41             3.117724964715322277790254593169e-8,
42             -5.354219639019687140874081024347e-9,
43             9.19327551985958894688778682594e-10,
44             -1.577941280288339761767423273953e-10,
45             2.707980622934954543266540433089e-11,
46             -4.646818653825730144081661058933e-12,
47             7.973350192007419656460767175359e-13,
48             -1.368078209830916025799499172309e-13,
49             2.347319486563800657233471771688e-14,
50             -4.027432614949066932766570534699e-15,
51             6.910051747372100912138336975257e-16,
52             -1.185584500221992907052387126192e-16,
53             2.034148542496373955201026051932e-17,
54             -3.490054341717405849274012949108e-18,
55             5.987993856485305567135051066026e-19,
56             -1.027378057872228074490069778431e-19,
57             1.762702816060529824942759660748e-20,
58             -3.024320653735306260958772112042e-21,
59             5.188914660218397839717833550506e-22,
60             -8.902770842456576692449251601066e-23,
61             1.527474068493342602274596891306e-23,
62             -2.620731256187362900257328332799e-24,
63             4.496464047830538670331046570666e-25,
64             -7.714712731336877911703901525333e-26,
65             1.323635453126044036486572714666e-26,
66             -2.270999412942928816702313813333e-27,
67             3.896418998003991449320816639999e-28,
68             -6.685198115125953327792127999999e-29,
69             1.146998663140024384347613866666e-29,
70             -1.967938586345134677295103999999e-30,
71             3.376448816585338090334890666666e-31,
72             -5.793070335782135784625493333333e-32 };
73     static doublereal pi = 3.1415926535897932384626433832795;
74     static doublereal sq2pil = .91893853320467274178032973640562;
75     // static logical first = TRUE_;
76 
77     /* System generated locals */
78     integer i__1;
79     real r__1;
80     doublereal ret_val, d__1, d__2;
81 
82     /* Builtin functions */
83     double sqrt(doublereal), d_int(doublereal *), log(doublereal), exp(
84             doublereal), sin(doublereal);
85 
86     /* Local variables */
87     integer i__, n;
88     doublereal y;
89     /* static */ integer ngam;
90     /* static */ doublereal xmin, xmax, dxrel;
91     extern doublereal d1mach_(integer *), d9lgmc_(doublereal *);
92     extern /* Subroutine */ int dgamlm_(doublereal *, doublereal *);
93     extern doublereal dcsevl_(doublereal *, doublereal *, integer *);
94     extern integer initds_(doublereal *, integer *, real *);
95     extern /* Subroutine */ int xermsg_(const char *, const char *, const char *, integer *,
96             integer *, ftnlen, ftnlen, ftnlen);
97     doublereal sinpiy;
98 
99 /* ***BEGIN PROLOGUE  DGAMMA */
100 /* ***PURPOSE  Compute the complete Gamma function. */
101 /* ***LIBRARY   SLATEC (FNLIB) */
102 /* ***CATEGORY  C7A */
103 /* ***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C) */
104 /* ***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS */
105 /* ***AUTHOR  Fullerton, W., (LANL) */
106 /* ***DESCRIPTION */
107 
108 /* DGAMMA(X) calculates the double precision complete Gamma function */
109 /* for double precision argument X. */
110 
111 /* Series for GAM        on the interval  0.          to  1.00000E+00 */
112 /*                                        with weighted error   5.79E-32 */
113 /*                                         log weighted error  31.24 */
114 /*                               significant figures required  30.00 */
115 /*                                    decimal places required  32.05 */
116 
117 /* ***REFERENCES  (NONE) */
118 /* ***ROUTINES CALLED  D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG */
119 /* ***REVISION HISTORY  (YYMMDD) */
120 /*   770601  DATE WRITTEN */
121 /*   890531  Changed all specific intrinsics to generic.  (WRB) */
122 /*   890911  Removed unnecessary intrinsics.  (WRB) */
123 /*   890911  REVISION DATE from Version 3.2 */
124 /*   891214  Prologue converted to Version 4.0 format.  (BAB) */
125 /*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ) */
126 /*   920618  Removed space from variable name.  (RWC, WRB) */
127 /* ***END PROLOGUE  DGAMMA */
128 
129 /* ***FIRST EXECUTABLE STATEMENT  DGAMMA */
130 
131     // d1mach has been made thread safe, so there is no need for the
132     // statics in determining these values
133 //     if (first) {
134 //      r__1 = (real) d1mach_(&c__3) * .1f;
135 //      ngam = initds_(gamcs, &c__42, &r__1);
136 
137 //      dgamlm_(&xmin, &xmax);
138 //      dxrel = sqrt(d1mach_(&c__4));
139 //     }
140 //     first = FALSE_;
141     r__1 = (real) d1mach_(&c__3) * .1f;
142     ngam = initds_(gamcs, &c__42, &r__1);
143     dgamlm_(&xmin, &xmax);
144     dxrel = sqrt(d1mach_(&c__4));
145 
146     y = abs(*x);
147     if (y > 10.) {
148         goto L50;
149     }
150 
151 /* COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND */
152 /* GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL. */
153 
154     n = (integer) (*x);
155     if (*x < 0.) {
156         --n;
157     }
158     y = *x - n;
159     --n;
160     d__1 = y * 2. - 1.;
161     ret_val = dcsevl_(&d__1, gamcs, &ngam) + .9375;
162     if (n == 0) {
163         return ret_val;
164     }
165 
166     if (n > 0) {
167         goto L30;
168     }
169 
170 /* COMPUTE GAMMA(X) FOR X .LT. 1.0 */
171 
172     n = -n;
173     if (*x == 0.) {
174         xermsg_("SLATEC", "DGAMMA", "X IS 0", &c__4, &c__2, (ftnlen)6, (
175                 ftnlen)6, (ftnlen)6);
176     }
177     if (*x < 0.f && *x + n - 2 == 0.) {
178         xermsg_("SLATEC", "DGAMMA", "X IS A NEGATIVE INTEGER", &c__4, &c__2, (
179                 ftnlen)6, (ftnlen)6, (ftnlen)23);
180     }
181     d__2 = *x - .5;
182     if (*x < -.5 && (d__1 = (*x - d_int(&d__2)) / *x, abs(d__1)) < dxrel) {
183         xermsg_("SLATEC", "DGAMMA", "ANSWER LT HALF PRECISION BECAUSE X TOO "
184                 "NEAR NEGATIVE INTEGER", &c__1, &c__1, (ftnlen)6, (ftnlen)6, (
185                 ftnlen)60);
186     }
187 
188     i__1 = n;
189     for (i__ = 1; i__ <= i__1; ++i__) {
190         ret_val /= *x + i__ - 1;
191 /* L20: */
192     }
193     return ret_val;
194 
195 /* GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0 */
196 
197 L30:
198     i__1 = n;
199     for (i__ = 1; i__ <= i__1; ++i__) {
200         ret_val = (y + i__) * ret_val;
201 /* L40: */
202     }
203     return ret_val;
204 
205 /* GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X). */
206 
207 L50:
208     if (*x > xmax) {
209         xermsg_("SLATEC", "DGAMMA", "X SO BIG GAMMA OVERFLOWS", &c__3, &c__2,
210                 (ftnlen)6, (ftnlen)6, (ftnlen)24);
211     }
212 
213     ret_val = 0.;
214     if (*x < xmin) {
215         xermsg_("SLATEC", "DGAMMA", "X SO SMALL GAMMA UNDERFLOWS", &c__2, &
216                 c__1, (ftnlen)6, (ftnlen)6, (ftnlen)27);
217     }
218     if (*x < xmin) {
219         return ret_val;
220     }
221 
222     ret_val = exp((y - .5) * log(y) - y + sq2pil + d9lgmc_(&y));
223     if (*x > 0.) {
224         return ret_val;
225     }
226 
227     d__2 = *x - .5;
228     if ((d__1 = (*x - d_int(&d__2)) / *x, abs(d__1)) < dxrel) {
229         xermsg_("SLATEC", "DGAMMA", "ANSWER LT HALF PRECISION, X TOO NEAR NE"
230                 "GATIVE INTEGER", &c__1, &c__1, (ftnlen)6, (ftnlen)6, (ftnlen)
231                 53);
232     }
233 
234     sinpiy = sin(pi * y);
235     if (sinpiy == 0.) {
236         xermsg_("SLATEC", "DGAMMA", "X IS A NEGATIVE INTEGER", &c__4, &c__2, (
237                 ftnlen)6, (ftnlen)6, (ftnlen)23);
238     }
239 
240     ret_val = -pi / (y * sinpiy * ret_val);
241 
242     return ret_val;
243 } /* dgamma_ */
244 
245