3  * for 128-bit long double.
5  * Copyright (c) 2006-2024, Arm Limited.
6  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
13  * 128-bit. Such a library will probably want to check the error
29 #define lenof(x) (sizeof(x)/sizeof(*(x)))  argument
34 static long double poly(const long double *coeffs, size_t n, long double x)  in poly()  argument
36     long double result = coeffs[--n];  in poly()
39         result = (result * x) + coeffs[--n];  in poly()
45  * Compute sin(pi*x) / pi, for use in the reflection formula that
46  * relates gamma(-x) and gamma(x).
48 static long double sin_pi_x_over_pi(long double x)  in sin_pi_x_over_pi()  argument
51     long double fracpart = remquol(x, 0.5L, &quo);  in sin_pi_x_over_pi()
55         sign = -sign;  in sin_pi_x_over_pi()
58     if (quo == 0 && fabsl(fracpart) < 0x1.p-58L) {  in sin_pi_x_over_pi()
59         /* For numbers this size, sin(pi*x) is so close to pi*x that  in sin_pi_x_over_pi()
60          * sin(pi*x)/pi is indistinguishable from x in float128 */  in sin_pi_x_over_pi()
71 /* Return tgamma(x) on the assumption that x >= 8. */
72 static long double tgamma_large(long double x,  in tgamma_large()  argument
76      * In this range we compute gamma(x) as x^(x-1/2) * e^-x * K,  in tgamma_large()
77      * where K is a correction factor computed as a polynomial in 1/x.  in tgamma_large()
83     long double t = 1/x;  in tgamma_large()
87      * To avoid overflow in cases where x^(x-0.5) does overflow  in tgamma_large()
88      * but gamma(x) does not, we split x^(x-0.5) in half and  in tgamma_large()
90      * of exp(-(x-0.5)).  in tgamma_large()
92      * Note that computing x-0.5 and (x-0.5)/2 is exact for the  in tgamma_large()
93      * relevant range of x, so the only sources of error are pow  in tgamma_large()
96     long double powhalf = powl(x, (x-0.5L)/2.0L);  in tgamma_large()
97     long double expret = expl(-(x-0.5L));  in tgamma_large()
105          * turn a case where gamma(+x) would overflow into a case  in tgamma_large()
106          * where gamma(-x) doesn't underflow. Not only that, but the  in tgamma_large()
111         long double ret = 1 / ((expret * powhalf) * (x * negadjust) * p);  in tgamma_large()
116 /* Return tgamma(x) on the assumption that 0 <= x < 1/32. */
117 static long double tgamma_tiny(long double x,  in tgamma_tiny()  argument
121      * For x near zero, we use a polynomial approximation to  in tgamma_tiny()
122      * g = 1/(x*gamma(x)), and then return 1/(g*x).  in tgamma_tiny()
124     long double g = poly(coeffs_tiny, lenof(coeffs_tiny), x);  in tgamma_tiny()
126         return 1.0L / (g*x);  in tgamma_tiny()
131 /* Return tgamma(x) on the assumption that 0 <= x < 2^-113. */
132 static long double tgamma_ultratiny(long double x, bool negative,  in tgamma_ultratiny()  argument
135     /* On this interval, gamma can't even be distinguished from 1/x,  in tgamma_ultratiny()
140         return 1.0L / x;  in tgamma_ultratiny()
145 /* Return tgamma(x) on the assumption that 1 <= x <= 2. */
146 static long double tgamma_central(long double x)  in tgamma_central()  argument
150      * difference between x and the point where gamma has a minimum,  in tgamma_central()
154     /* The difference between the input x and the minimum x. The first  in tgamma_central()
155      * subtraction is expected to be exact, since x and min_hi have  in tgamma_central()
156      * the same exponent (unless x=2, in which case it will still be  in tgamma_central()
158     long double t = (x - min_x_hi) - min_x_lo;  in tgamma_central()
166         p = poly(coeffs_central_neg, lenof(coeffs_central_neg), -t);  in tgamma_central()
173 long double tgamma128(long double x)  in tgamma128()  argument
177      * out cases of non-normalized numbers.  in tgamma128()
188     int sign = signbit(x);  in tgamma128()
190     switch (fpclassify(x)) {  in tgamma128()
192         return x+x; /* propagate QNaN, make SNaN throw an exception */  in tgamma128()
194         return 1/x; /* divide by zero on purpose to indicate a pole */  in tgamma128()
197             return x-x; /* gamma(-inf) has indeterminate sign, so provoke an  in tgamma128()
200         return x;     /* but gamma(+inf) is just +inf with no error */  in tgamma128()
202         exponent = -16384;  in tgamma128()
205         frexpl(x, &exponent);  in tgamma128()
206         exponent--;  in tgamma128()
217          *    gamma(1-x) gamma(x) = pi/sin(pi*x)  in tgamma128()
220          * => gamma(x) = --------------------  in tgamma128()
221          *               gamma(1-x) sin(pi*x)  in tgamma128()
223          * But computing 1-x is going to lose a lot of accuracy when x  in tgamma128()
225          * gamma(t+1)=t gamma(t). Setting t=-x, this gives us  in tgamma128()
226          * gamma(1-x) = -x gamma(-x), so we now have  in tgamma128()
229          *    gamma(x) = ----------------------  in tgamma128()
230          *               -x gamma(-x) sin(pi*x)  in tgamma128()
232          * which relates gamma(x) to gamma(-x), which is much nicer,  in tgamma128()
233          * since x can be turned into -x without rounding.  in tgamma128()
235         negadjust = sin_pi_x_over_pi(x);  in tgamma128()
237         x = -x;  in tgamma128()
242          *    1 / (gamma(x) * x * negadjust)  in tgamma128()
244          * where x is the positive value we've just turned it into.  in tgamma128()
246          * For some of the cases below, we'll compute gamma(x)  in tgamma128()
250          * sub-case was going to do anyway.  in tgamma128()
260              * infinity to return, and unlike the x=0 case, there's no  in tgamma128()
269      * we do the negative-number adjustment the usual way, we'll leave  in tgamma128()
283             long double tiny = 0x1p-12288L;  in tgamma128()
287         /* Negative-number adjustment happens inside here */  in tgamma128()
288         return tgamma_large(x, negative, negadjust);  in tgamma128()
289     } else if (exponent < -113) {  in tgamma128()
290         /* Negative-number adjustment happens inside here */  in tgamma128()
291         return tgamma_ultratiny(x, negative, negadjust);  in tgamma128()
292     } else if (exponent < -5) {  in tgamma128()
293         /* Negative-number adjustment happens inside here */  in tgamma128()
294         return tgamma_tiny(x, negative, negadjust);  in tgamma128()
296         g = tgamma_central(x);  in tgamma128()
299          * For x in [1/32,1) we range-reduce upwards to the interval  in tgamma128()
301          * gamma(x) = gamma(x+1)/x.  in tgamma128()
303         g = tgamma_central(1+x) / x;  in tgamma128()
306          * For x in [2,8) we range-reduce downwards to the interval  in tgamma128()
309          * Actually multiplying (x-1) by (x-2) by (x-3) and so on  in tgamma128()
311          * better accuracy by writing x = (k+1/2) + t, where k is an  in tgamma128()
313          * (x-1)(x-2)...(x-k+1) as a polynomial in t.  in tgamma128()
316         int i = x;  in tgamma128()
317         if (i == 2) { /* x in [2,3) */  in tgamma128()
318             mult = (x-1);  in tgamma128()
320             long double t = x - (i + 0.5L);  in tgamma128()
322                 /* E.g. for x=3.5+t, we want  in tgamma128()
323                  * (x-1)(x-2) = (2.5+t)(1.5+t) = 3.75 + 4t + t^2 */  in tgamma128()
346         g = tgamma_central(x - (i-1)) * mult;  in tgamma128()
354         return 1.0L / (g * x * negadjust);  in tgamma128()