Branch data Line data Source code
1 : : /*-
2 : : * Copyright (c) 1992, 1993
3 : : * The Regents of the University of California. All rights reserved.
4 : : *
5 : : * Redistribution and use in source and binary forms, with or without
6 : : * modification, are permitted provided that the following conditions
7 : : * are met:
8 : : * 1. Redistributions of source code must retain the above copyright
9 : : * notice, this list of conditions and the following disclaimer.
10 : : * 2. Redistributions in binary form must reproduce the above copyright
11 : : * notice, this list of conditions and the following disclaimer in the
12 : : * documentation and/or other materials provided with the distribution.
13 : : * 3. Neither the name of the University nor the names of its contributors
14 : : * may be used to endorse or promote products derived from this software
15 : : * without specific prior written permission.
16 : : *
17 : : * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
18 : : * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 : : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 : : * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
21 : : * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 : : * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 : : * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 : : * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 : : * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 : : * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 : : * SUCH DAMAGE.
28 : : */
29 : :
30 : : /* @(#)gamma.c 8.1 (Berkeley) 6/4/93 */
31 : : #include "cdefs-compat.h"
32 : : //__FBSDID("$FreeBSD: src/lib/msun/bsdsrc/b_tgamma.c,v 1.10 2008/02/22 02:26:51 das Exp $");
33 : :
34 : : /*
35 : : * This code by P. McIlroy, Oct 1992;
36 : : *
37 : : * The financial support of UUNET Communications Services is greatfully
38 : : * acknowledged.
39 : : */
40 : :
41 : : #include <float.h>
42 : : #include <openlibm_math.h>
43 : :
44 : : #include "mathimpl.h"
45 : :
46 : : /* METHOD:
47 : : * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
48 : : * At negative integers, return NaN and raise invalid.
49 : : *
50 : : * x < 6.5:
51 : : * Use argument reduction G(x+1) = xG(x) to reach the
52 : : * range [1.066124,2.066124]. Use a rational
53 : : * approximation centered at the minimum (x0+1) to
54 : : * ensure monotonicity.
55 : : *
56 : : * x >= 6.5: Use the asymptotic approximation (Stirling's formula)
57 : : * adjusted for equal-ripples:
58 : : *
59 : : * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x))
60 : : *
61 : : * Keep extra precision in multiplying (x-.5)(log(x)-1), to
62 : : * avoid premature round-off.
63 : : *
64 : : * Special values:
65 : : * -Inf: return NaN and raise invalid;
66 : : * negative integer: return NaN and raise invalid;
67 : : * other x ~< 177.79: return +-0 and raise underflow;
68 : : * +-0: return +-Inf and raise divide-by-zero;
69 : : * finite x ~> 171.63: return +Inf and raise overflow;
70 : : * +Inf: return +Inf;
71 : : * NaN: return NaN.
72 : : *
73 : : * Accuracy: tgamma(x) is accurate to within
74 : : * x > 0: error provably < 0.9ulp.
75 : : * Maximum observed in 1,000,000 trials was .87ulp.
76 : : * x < 0:
77 : : * Maximum observed error < 4ulp in 1,000,000 trials.
78 : : */
79 : :
80 : : static double neg_gam(double);
81 : : static double small_gam(double);
82 : : static double smaller_gam(double);
83 : : static struct Double large_gam(double);
84 : : static struct Double ratfun_gam(double, double);
85 : :
86 : : /*
87 : : * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval
88 : : * [1.066.., 2.066..] accurate to 4.25e-19.
89 : : */
90 : : #define LEFT -.3955078125 /* left boundary for rat. approx */
91 : : #define x0 .461632144968362356785 /* xmin - 1 */
92 : :
93 : : #define a0_hi 0.88560319441088874992
94 : : #define a0_lo -.00000000000000004996427036469019695
95 : : #define P0 6.21389571821820863029017800727e-01
96 : : #define P1 2.65757198651533466104979197553e-01
97 : : #define P2 5.53859446429917461063308081748e-03
98 : : #define P3 1.38456698304096573887145282811e-03
99 : : #define P4 2.40659950032711365819348969808e-03
100 : : #define Q0 1.45019531250000000000000000000e+00
101 : : #define Q1 1.06258521948016171343454061571e+00
102 : : #define Q2 -2.07474561943859936441469926649e-01
103 : : #define Q3 -1.46734131782005422506287573015e-01
104 : : #define Q4 3.07878176156175520361557573779e-02
105 : : #define Q5 5.12449347980666221336054633184e-03
106 : : #define Q6 -1.76012741431666995019222898833e-03
107 : : #define Q7 9.35021023573788935372153030556e-05
108 : : #define Q8 6.13275507472443958924745652239e-06
109 : : /*
110 : : * Constants for large x approximation (x in [6, Inf])
111 : : * (Accurate to 2.8*10^-19 absolute)
112 : : */
113 : : #define lns2pi_hi 0.418945312500000
114 : : #define lns2pi_lo -.000006779295327258219670263595
115 : : #define Pa0 8.33333333333333148296162562474e-02
116 : : #define Pa1 -2.77777777774548123579378966497e-03
117 : : #define Pa2 7.93650778754435631476282786423e-04
118 : : #define Pa3 -5.95235082566672847950717262222e-04
119 : : #define Pa4 8.41428560346653702135821806252e-04
120 : : #define Pa5 -1.89773526463879200348872089421e-03
121 : : #define Pa6 5.69394463439411649408050664078e-03
122 : : #define Pa7 -1.44705562421428915453880392761e-02
123 : :
124 : : static const double zero = 0., one = 1.0, tiny = 1e-300;
125 : :
126 : : OLM_DLLEXPORT double
127 : 28 : tgamma(x)
128 : : double x;
129 : : {
130 : : struct Double u;
131 : :
132 [ + + ]: 28 : if (isgreaterequal(x, 6)) {
133 [ + - ]: 2 : if(x > 171.63)
134 : 2 : return (x / zero);
135 : 0 : u = large_gam(x);
136 : 0 : return(__exp__D(u.a, u.b));
137 [ + + ]: 26 : } else if (isgreaterequal(x, 1.0 + LEFT + x0))
138 : 6 : return (small_gam(x));
139 [ + + ]: 20 : else if (isgreater(x, 1.e-17))
140 : 8 : return (smaller_gam(x));
141 [ + + ]: 12 : else if (isgreater(x, -1.e-17)) {
142 [ - + ]: 4 : if (x != 0.0)
143 : 0 : u.a = one - tiny; /* raise inexact */
144 : 4 : return (one/x);
145 [ + + ]: 8 : } else if (!isfinite(x))
146 : 4 : return (x - x); /* x is NaN or -Inf */
147 : : else
148 : 4 : return (neg_gam(x));
149 : : }
150 : : /*
151 : : * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
152 : : */
153 : : static struct Double
154 : 0 : large_gam(x)
155 : : double x;
156 : : {
157 : : double z, p;
158 : : struct Double t, u, v;
159 : :
160 : 0 : z = one/(x*x);
161 : 0 : p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7))))));
162 : 0 : p = p/x;
163 : :
164 : 0 : u = __log__D(x);
165 : 0 : u.a -= one;
166 : 0 : v.a = (x -= .5);
167 : 0 : TRUNC(v.a);
168 : 0 : v.b = x - v.a;
169 : 0 : t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */
170 : 0 : t.b = v.b*u.a + x*u.b;
171 : : /* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */
172 : 0 : t.b += lns2pi_lo; t.b += p;
173 : 0 : u.a = lns2pi_hi + t.b; u.a += t.a;
174 : 0 : u.b = t.a - u.a;
175 : 0 : u.b += lns2pi_hi; u.b += t.b;
176 : 0 : return (u);
177 : : }
178 : : /*
179 : : * Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.)
180 : : * It also has correct monotonicity.
181 : : */
182 : : static double
183 : 6 : small_gam(x)
184 : : double x;
185 : : {
186 : : double y, ym1, t;
187 : : struct Double yy, r;
188 : 6 : y = x - one;
189 : 6 : ym1 = y - one;
190 [ + + ]: 6 : if (y <= 1.0 + (LEFT + x0)) {
191 : 4 : yy = ratfun_gam(y - x0, 0);
192 : 4 : return (yy.a + yy.b);
193 : : }
194 : 2 : r.a = y;
195 : 2 : TRUNC(r.a);
196 : 2 : yy.a = r.a - one;
197 : 2 : y = ym1;
198 : 2 : yy.b = r.b = y - yy.a;
199 : : /* Argument reduction: G(x+1) = x*G(x) */
200 [ + + ]: 4 : for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
201 : 2 : t = r.a*yy.a;
202 : 2 : r.b = r.a*yy.b + y*r.b;
203 : 2 : r.a = t;
204 : 2 : TRUNC(r.a);
205 : 2 : r.b += (t - r.a);
206 : : }
207 : : /* Return r*tgamma(y). */
208 : 2 : yy = ratfun_gam(y - x0, 0);
209 : 2 : y = r.b*(yy.a + yy.b) + r.a*yy.b;
210 : 2 : y += yy.a*r.a;
211 : 2 : return (y);
212 : : }
213 : : /*
214 : : * Good on (0, 1+x0+LEFT]. Accurate to 1ulp.
215 : : */
216 : : static double
217 : 8 : smaller_gam(x)
218 : : double x;
219 : : {
220 : : double t, d;
221 : : struct Double r, xx;
222 [ - + ]: 8 : if (x < x0 + LEFT) {
223 : 0 : t = x, TRUNC(t);
224 : 0 : d = (t+x)*(x-t);
225 : 0 : t *= t;
226 : 0 : xx.a = (t + x), TRUNC(xx.a);
227 : 0 : xx.b = x - xx.a; xx.b += t; xx.b += d;
228 : 0 : t = (one-x0); t += x;
229 : 0 : d = (one-x0); d -= t; d += x;
230 : 0 : x = xx.a + xx.b;
231 : : } else {
232 : 8 : xx.a = x, TRUNC(xx.a);
233 : 8 : xx.b = x - xx.a;
234 : 8 : t = x - x0;
235 : 8 : d = (-x0 -t); d += x;
236 : : }
237 : 8 : r = ratfun_gam(t, d);
238 : 8 : d = r.a/x, TRUNC(d);
239 : 8 : r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b;
240 : 8 : return (d + r.a/x);
241 : : }
242 : : /*
243 : : * returns (z+c)^2 * P(z)/Q(z) + a0
244 : : */
245 : : static struct Double
246 : 14 : ratfun_gam(z, c)
247 : : double z, c;
248 : : {
249 : : double p, q;
250 : : struct Double r, t;
251 : :
252 : 14 : q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8)))))));
253 : 14 : p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4)));
254 : :
255 : : /* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */
256 : 14 : p = p/q;
257 : 14 : t.a = z, TRUNC(t.a); /* t ~= z + c */
258 : 14 : t.b = (z - t.a) + c;
259 : 14 : t.b *= (t.a + z);
260 : 14 : q = (t.a *= t.a); /* t = (z+c)^2 */
261 : 14 : TRUNC(t.a);
262 : 14 : t.b += (q - t.a);
263 : 14 : r.a = p, TRUNC(r.a); /* r = P/Q */
264 : 14 : r.b = p - r.a;
265 : 14 : t.b = t.b*p + t.a*r.b + a0_lo;
266 : 14 : t.a *= r.a; /* t = (z+c)^2*(P/Q) */
267 : 14 : r.a = t.a + a0_hi, TRUNC(r.a);
268 : 14 : r.b = ((a0_hi-r.a) + t.a) + t.b;
269 : 14 : return (r); /* r = a0 + t */
270 : : }
271 : :
272 : : static double
273 : 4 : neg_gam(x)
274 : : double x;
275 : : {
276 : 4 : int sgn = 1;
277 : : struct Double lg, lsine;
278 : : double y, z;
279 : :
280 : 4 : y = ceil(x);
281 [ + + ]: 4 : if (y == x) /* Negative integer. */
282 : 2 : return ((x - x) / zero);
283 : 2 : z = y - x;
284 [ - + ]: 2 : if (z > 0.5)
285 : 0 : z = one - z;
286 : 2 : y = 0.5 * y;
287 [ + - ]: 2 : if (y == ceil(y))
288 : 2 : sgn = -1;
289 [ - + ]: 2 : if (z < .25)
290 : 0 : z = sin(M_PI*z);
291 : : else
292 : 2 : z = cos(M_PI*(0.5-z));
293 : : /* Special case: G(1-x) = Inf; G(x) may be nonzero. */
294 [ - + ]: 2 : if (x < -170) {
295 [ # # ]: 0 : if (x < -190)
296 : 0 : return ((double)sgn*tiny*tiny);
297 : 0 : y = one - x; /* exact: 128 < |x| < 255 */
298 : 0 : lg = large_gam(y);
299 : 0 : lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */
300 : 0 : lg.a -= lsine.a; /* exact (opposite signs) */
301 : 0 : lg.b -= lsine.b;
302 : 0 : y = -(lg.a + lg.b);
303 : 0 : z = (y + lg.a) + lg.b;
304 : 0 : y = __exp__D(y, z);
305 [ # # ]: 0 : if (sgn < 0) y = -y;
306 : 0 : return (y);
307 : : }
308 : 2 : y = one-x;
309 [ + - ]: 2 : if (one-y == x)
310 : 2 : y = tgamma(y);
311 : : else /* 1-x is inexact */
312 : 0 : y = -x*tgamma(-x);
313 [ + - ]: 2 : if (sgn < 0) y = -y;
314 : 2 : return (M_PI / (y*z));
315 : : }
316 : :
317 : : #if (LDBL_MANT_DIG == 53)
318 : : openlibm_weak_reference(tgamma, tgammal);
319 : : #endif
|
