/src/r-source/src/nmath/gamma_cody.c
Line | Count | Source |
1 | | /* From http://www.netlib.org/specfun/gamma Fortran translated by f2c,... |
2 | | * ------------------------------##### Martin Maechler, ETH Zurich |
3 | | * |
4 | | *=========== was part of ribesl (Bessel I(.)) |
5 | | *=========== ~~~~~~ |
6 | | */ |
7 | | |
8 | | // used in bessel_i.c and bessel_j.c, hidden if possible. |
9 | | |
10 | | #include "nmath.h" |
11 | | |
12 | | attribute_hidden double Rf_gamma_cody(double x) |
13 | 0 | { |
14 | | /* ---------------------------------------------------------------------- |
15 | | |
16 | | This routine calculates the GAMMA function for a float argument X. |
17 | | Computation is based on an algorithm outlined in reference [1]. |
18 | | The program uses rational functions that approximate the GAMMA |
19 | | function to at least 20 significant decimal digits. Coefficients |
20 | | for the approximation over the interval (1,2) are unpublished. |
21 | | Those for the approximation for X >= 12 are from reference [2]. |
22 | | The accuracy achieved depends on the arithmetic system, the |
23 | | compiler, the intrinsic functions, and proper selection of the |
24 | | machine-dependent constants. |
25 | | |
26 | | ******************************************************************* |
27 | | |
28 | | Error returns |
29 | | |
30 | | The program returns the value XINF for singularities or |
31 | | when overflow would occur. The computation is believed |
32 | | to be free of underflow and overflow. |
33 | | |
34 | | Intrinsic functions required are: |
35 | | |
36 | | INT, DBLE, EXP, LOG, REAL, SIN |
37 | | |
38 | | |
39 | | References: |
40 | | [1] "An Overview of Software Development for Special Functions", |
41 | | W. J. Cody, Lecture Notes in Mathematics, 506, |
42 | | Numerical Analysis Dundee, 1975, G. A. Watson (ed.), |
43 | | Springer Verlag, Berlin, 1976. |
44 | | |
45 | | [2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968. |
46 | | |
47 | | Latest modification: October 12, 1989 |
48 | | |
49 | | Authors: W. J. Cody and L. Stoltz |
50 | | Applied Mathematics Division |
51 | | Argonne National Laboratory |
52 | | Argonne, IL 60439 |
53 | | ----------------------------------------------------------------------*/ |
54 | | |
55 | | /* ---------------------------------------------------------------------- |
56 | | Mathematical constants |
57 | | ----------------------------------------------------------------------*/ |
58 | 0 | const static double sqrtpi = .9189385332046727417803297; /* == ??? */ |
59 | | |
60 | | /* ******************************************************************* |
61 | | |
62 | | Explanation of machine-dependent constants |
63 | | |
64 | | beta - radix for the floating-point representation |
65 | | maxexp - the smallest positive power of beta that overflows |
66 | | XBIG - the largest argument for which GAMMA(X) is representable |
67 | | in the machine, i.e., the solution to the equation |
68 | | GAMMA(XBIG) = beta**maxexp |
69 | | XINF - the largest machine representable floating-point number; |
70 | | approximately beta**maxexp |
71 | | EPS - the smallest positive floating-point number such that 1.0+EPS > 1.0 |
72 | | XMININ - the smallest positive floating-point number such that |
73 | | 1/XMININ is machine representable |
74 | | |
75 | | Approximate values for some important machines are: |
76 | | |
77 | | beta maxexp XBIG |
78 | | |
79 | | CRAY-1 (S.P.) 2 8191 966.961 |
80 | | Cyber 180/855 |
81 | | under NOS (S.P.) 2 1070 177.803 |
82 | | IEEE (IBM/XT, |
83 | | SUN, etc.) (S.P.) 2 128 35.040 |
84 | | IEEE (IBM/XT, |
85 | | SUN, etc.) (D.P.) 2 1024 171.624 |
86 | | IBM 3033 (D.P.) 16 63 57.574 |
87 | | VAX D-Format (D.P.) 2 127 34.844 |
88 | | VAX G-Format (D.P.) 2 1023 171.489 |
89 | | |
90 | | XINF EPS XMININ |
91 | | |
92 | | CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 |
93 | | Cyber 180/855 |
94 | | under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 |
95 | | IEEE (IBM/XT, |
96 | | SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 |
97 | | IEEE (IBM/XT, |
98 | | SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 |
99 | | IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 |
100 | | VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 |
101 | | VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 |
102 | | |
103 | | ******************************************************************* |
104 | | |
105 | | ---------------------------------------------------------------------- |
106 | | Machine dependent parameters |
107 | | ---------------------------------------------------------------------- |
108 | | */ |
109 | | |
110 | |
|
111 | 0 | const static double xbig = 171.624; |
112 | | /* ML_POSINF == const double xinf = 1.79e308;*/ |
113 | | /* DBL_EPSILON = const double eps = 2.22e-16;*/ |
114 | | /* DBL_MIN == const double xminin = 2.23e-308;*/ |
115 | | |
116 | | /*---------------------------------------------------------------------- |
117 | | Numerator and denominator coefficients for rational minimax |
118 | | approximation over (1,2). |
119 | | ----------------------------------------------------------------------*/ |
120 | 0 | const static double p[8] = { |
121 | 0 | -1.71618513886549492533811, |
122 | 0 | 24.7656508055759199108314,-379.804256470945635097577, |
123 | 0 | 629.331155312818442661052,866.966202790413211295064, |
124 | 0 | -31451.2729688483675254357,-36144.4134186911729807069, |
125 | 0 | 66456.1438202405440627855 }; |
126 | 0 | const static double q[8] = { |
127 | 0 | -30.8402300119738975254353, |
128 | 0 | 315.350626979604161529144,-1015.15636749021914166146, |
129 | 0 | -3107.77167157231109440444,22538.1184209801510330112, |
130 | 0 | 4755.84627752788110767815,-134659.959864969306392456, |
131 | 0 | -115132.259675553483497211 }; |
132 | | /*---------------------------------------------------------------------- |
133 | | Coefficients for minimax approximation over (12, INF). |
134 | | ----------------------------------------------------------------------*/ |
135 | 0 | const static double c[7] = { |
136 | 0 | -.001910444077728,8.4171387781295e-4, |
137 | 0 | -5.952379913043012e-4,7.93650793500350248e-4, |
138 | 0 | -.002777777777777681622553,.08333333333333333331554247, |
139 | 0 | .0057083835261 }; |
140 | | |
141 | | /* Local variables */ |
142 | 0 | int i, n; |
143 | 0 | int parity;/*logical*/ |
144 | 0 | double fact, xden, xnum, y, z, yi, res, sum, ysq; |
145 | |
|
146 | 0 | parity = (0); |
147 | 0 | fact = 1.; |
148 | 0 | n = 0; |
149 | 0 | y = x; |
150 | 0 | if (y <= 0.) { |
151 | | /* ------------------------------------------------------------- |
152 | | Argument is negative |
153 | | ------------------------------------------------------------- */ |
154 | 0 | y = -x; |
155 | 0 | yi = trunc(y); |
156 | 0 | res = y - yi; |
157 | 0 | if (res != 0.) { |
158 | 0 | if (yi != trunc(yi * .5) * 2.) |
159 | 0 | parity = (1); |
160 | 0 | fact = -M_PI / sinpi(res); |
161 | 0 | y += 1.; |
162 | 0 | } else { |
163 | 0 | return(ML_POSINF); |
164 | 0 | } |
165 | 0 | } |
166 | | /* ----------------------------------------------------------------- |
167 | | Argument is positive |
168 | | -----------------------------------------------------------------*/ |
169 | 0 | if (y < DBL_EPSILON) { |
170 | | /* -------------------------------------------------------------- |
171 | | Argument < EPS |
172 | | -------------------------------------------------------------- */ |
173 | 0 | if (y >= DBL_MIN) { |
174 | 0 | res = 1. / y; |
175 | 0 | } else { |
176 | 0 | return(ML_POSINF); |
177 | 0 | } |
178 | 0 | } else if (y < 12.) { |
179 | 0 | yi = y; |
180 | 0 | if (y < 1.) { |
181 | | /* --------------------------------------------------------- |
182 | | EPS < argument < 1 |
183 | | --------------------------------------------------------- */ |
184 | 0 | z = y; |
185 | 0 | y += 1.; |
186 | 0 | } else { |
187 | | /* ----------------------------------------------------------- |
188 | | 1 <= argument < 12, reduce argument if necessary |
189 | | ----------------------------------------------------------- */ |
190 | 0 | n = (int) y - 1; |
191 | 0 | y -= (double) n; |
192 | 0 | z = y - 1.; |
193 | 0 | } |
194 | | /* --------------------------------------------------------- |
195 | | Evaluate approximation for 1. < argument < 2. |
196 | | ---------------------------------------------------------*/ |
197 | 0 | xnum = 0.; |
198 | 0 | xden = 1.; |
199 | 0 | for (i = 0; i < 8; ++i) { |
200 | 0 | xnum = (xnum + p[i]) * z; |
201 | 0 | xden = xden * z + q[i]; |
202 | 0 | } |
203 | 0 | res = xnum / xden + 1.; |
204 | 0 | if (yi < y) { |
205 | | /* -------------------------------------------------------- |
206 | | Adjust result for case 0. < argument < 1. |
207 | | -------------------------------------------------------- */ |
208 | 0 | res /= yi; |
209 | 0 | } else if (yi > y) { |
210 | | /* ---------------------------------------------------------- |
211 | | Adjust result for case 2. < argument < 12. |
212 | | ---------------------------------------------------------- */ |
213 | 0 | for (i = 0; i < n; ++i) { |
214 | 0 | res *= y; |
215 | 0 | y += 1.; |
216 | 0 | } |
217 | 0 | } |
218 | 0 | } else { |
219 | | /* ------------------------------------------------------------- |
220 | | Evaluate for argument >= 12., |
221 | | ------------------------------------------------------------- */ |
222 | 0 | if (y <= xbig) { |
223 | 0 | ysq = y * y; |
224 | 0 | sum = c[6]; |
225 | 0 | for (i = 0; i < 6; ++i) { |
226 | 0 | sum = sum / ysq + c[i]; |
227 | 0 | } |
228 | 0 | sum = sum / y - y + sqrtpi; |
229 | 0 | sum += (y - .5) * log(y); |
230 | 0 | res = exp(sum); |
231 | 0 | } else { |
232 | 0 | return(ML_POSINF); |
233 | 0 | } |
234 | 0 | } |
235 | | /* ---------------------------------------------------------------------- |
236 | | Final adjustments and return |
237 | | ----------------------------------------------------------------------*/ |
238 | 0 | if (parity) |
239 | 0 | res = -res; |
240 | 0 | if (fact != 1.) |
241 | 0 | res = fact / res; |
242 | 0 | return res; |
243 | 0 | } |
244 | | |