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