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