1 /*-
2  * SPDX-License-Identifier: BSD-3-Clause
15  * 3. Neither the name of the University nor the names of its contributors
72 static const volatile double tiny = 1e-300;
77  * equal-ripples:
79  * log(G(x)) ~= (x-0.5)*(log(x)-1) + 0.5(log(2*pi)-1) + 1/x*P(1/(x*x))
81  * Keep extra precision in multiplying (x-.5)(log(x)-1), to avoid
82  * premature round-off.
84  * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
88  * The following is a decomposition of 0.5 * (log(2*pi) - 1) into the
93 ln2pi_hiu = LD80C(0xd680000000000000,  -2,  4.18945312500000000000e-01L),
94 ln2pi_lou = LD80C(0xe379b414b596d687, -18, -6.77929532725821967032e-06L);
99     Pa0u = LD80C(0xaaaaaaaaaaaaaaaa,  -4,  8.33333333333333333288e-02L),
100     Pa1u = LD80C(0xb60b60b60b5fcd59,  -9, -2.77777777777776516326e-03L),
101     Pa2u = LD80C(0xd00d00cffbb47014, -11,  7.93650793635429639018e-04L),
102     Pa3u = LD80C(0x9c09c07c0805343e, -11, -5.95238087960599252215e-04L),
103     Pa4u = LD80C(0xdca8d31f8e6e5e8f, -11,  8.41749082509607342883e-04L),
104     Pa5u = LD80C(0xfb4d4289632f1638, -10, -1.91728055205541624556e-03L),
105     Pa6u = LD80C(0xd15a4ba04078d3f8,  -8,  6.38893788027752396194e-03L),
106     Pa7u = LD80C(0xe877283110bcad95,  -6, -2.83771309846297590312e-02L),
107     Pa8u = LD80C(0x8da97eed13717af8,  -3,  1.38341887683837576925e-01L),
108     Pa9u = LD80C(0xf093b1c1584e30ce,  -2, -4.69876818515470146031e-01L);
133 	u.a -= 1;  in large_gam()
135 	/* Split (x - 0.5) in high and low parts. */  in large_gam()
136 	x -= 0.5L;  in large_gam()
138 	xlo = x - xhi;  in large_gam()
140 	/* Compute  t = (x-.5)*(log(x)-1) in extra precision. */  in large_gam()
149 	u.b = thi - u.a;  in large_gam()
156  * [1.066.., 2.066..] accurate to 4.25e-19.
161     a0_hiu = LD80C(0xe2b6e4153a57746c,  -1, 8.85603194410888700265e-01L),
162     a0_lou = LD80C(0x851566d40f32c76d, -66, 1.40907742727049706207e-20L);
167 P0u = LD80C(0xdb629fb9bbdc1c1d,    -2,  4.28486815855585429733e-01L),
168 P1u = LD80C(0xe6f4f9f5641aa6be,    -3,  2.25543885805587730552e-01L),
169 P2u = LD80C(0xead1bd99fdaf7cc1,    -6,  2.86644652514293482381e-02L),  variable
170 P3u = LD80C(0x9ccc8b25838ab1e0,    -8,  4.78512567772456362048e-03L),
171 P4u = LD80C(0x8f0c4383ef9ce72a,    -9,  2.18273781132301146458e-03L),
172 P5u = LD80C(0xe732ab2c0a2778da,   -13,  2.20487522485636008928e-04L),
173 P6u = LD80C(0xce70b27ca822b297,   -16,  2.46095923774929264284e-05L),
174 P7u = LD80C(0xa309e2e16fb63663,   -19,  2.42946473022376182921e-06L),
175 P8u = LD80C(0xaf9c110efb2c633d,   -23,  1.63549217667765869987e-07L),
176 Q1u = LD80C(0xd4d7422719f48f15,    -1,  8.31409582658993993626e-01L),
177 Q2u = LD80C(0xe13138ea404f1268,    -5, -5.49785826915643198508e-02L),
178 Q3u = LD80C(0xd1c6cc91989352c0,    -4, -1.02429960435139887683e-01L),
179 Q4u = LD80C(0xa7e9435a84445579,    -7,  1.02484853505908820524e-02L),
180 Q5u = LD80C(0x83c7c34db89b7bda,    -8,  4.02161632832052872697e-03L),
181 Q6u = LD80C(0xbed06bf6e1c14e5b,   -11, -7.27898206351223022157e-04L),
182 Q7u = LD80C(0xef05bf841d4504c0,   -18,  7.12342421869453515194e-06L),
183 Q8u = LD80C(0xf348d08a1ff53cb1,   -19,  3.62522053809474067060e-06L);
186 #define	P2	(P2u.e)
216 	tlo = (z - thi) + c;  in ratfun_gam()
223 	tlo += (q - thi);  in ratfun_gam()
227 	r.b = p - r.a;  in ratfun_gam()
232 	r.b = ((a0_hi - r.a) + thi) + tlo;  in ratfun_gam()
246   xm1u = LD80C(0xec5b0c6ad7c7edc3, -2, 4.61632144968362341254e-01L);
250     left = -0.3955078125;	/* left boundary for rat. approx */
258 	y = x - 1;  in small_gam()
261 		yy = ratfun_gam(y - x0, 0);  in small_gam()
266 	yy.a = r.a - 1;  in small_gam()
267 	y = y - 1 ;  in small_gam()
268 	r.b = yy.b = y - yy.a;  in small_gam()
271 	for (ym1 = y - 1; ym1 > left + x0; y = ym1--, yy.a--) {  in small_gam()
275 		r.b += (t - r.a);  in small_gam()
279 	yy = ratfun_gam(y - x0, 0);  in small_gam()
295 		d = (t + x) * (x - t);  in smaller_gam()
298 		xlo = x - xhi;  in smaller_gam()
301 		t = 1 - x0;  in smaller_gam()
303 		d = 1 - x0;  in smaller_gam()
304 		d -= t;  in smaller_gam()
309 		xlo = x - xhi;  in smaller_gam()
310 		t = x - x0;  in smaller_gam()
311 		d = - x0 - t;  in smaller_gam()
317 	r.a -= d * xhi;  in smaller_gam()
318 	r.a -= d * xlo;  in smaller_gam()
342 		return ((x - x) / zero);  in neg_gam()
344 	z = y - x;  in neg_gam()
346 		z = 1 - z;  in neg_gam()
350 		sgn = -1;  in neg_gam()
355 		z = cospil(0.5 - z);  in neg_gam()
357 	/* Special case: G(1-x) = Inf; G(x) may be nonzero. */  in neg_gam()
358 	if (x < -1753) {  in neg_gam()
360 		if (x < -1760)  in neg_gam()
364 		return (sgn < 0 ? -y : y);  in neg_gam()
368 	y = 1 - x;  in neg_gam()
369 	if (1 - y == x)  in neg_gam()
371 	else		/* 1-x is inexact */  in neg_gam()
372 		y = - x * tgammal(-x);  in neg_gam()
374 	if (sgn < 0) y = -y;  in neg_gam()
378  * xmax comes from lgamma(xmax) - emax * log(2) = 0.
387 static const double iota = 0x1p-116;
409 	if (x > -iota) {  in tgammal()
411 			u.a = 1 - tiny;	/* raise inexact */  in tgammal()
416 		RETURNI(x - x);		/* x is NaN or -Inf */  in tgammal()