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
12  * maths libraries under the standard name tgammal, if long double is
13  * 128-bit. Such a library will probably want to check the error
34 static long double poly(const long double *coeffs, size_t n, long double x)  in poly()
36     long double result = coeffs[--n];  in poly()
39         result = (result * x) + coeffs[--n];  in poly()
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()
51     long double fracpart = remquol(x, 0.5L, &quo);  in sin_pi_x_over_pi()
53     long double sign = 1.0L;  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()
72 static long double tgamma_large(long double x,  in tgamma_large()
73                                 bool negative, long double negadjust)  in tgamma_large()
76      * In this range we compute gamma(x) as x^(x-1/2) * e^-x * K,  in tgamma_large()
83     long double t = 1/x;  in tgamma_large()
84     long double p = poly(coeffs_large, lenof(coeffs_large), t);  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()
106          * where gamma(-x) doesn't underflow. Not only that, but the  in tgamma_large()
107          * FP format has greater range in the tiny domain due to  in tgamma_large()
111         long double ret = 1 / ((expret * powhalf) * (x * negadjust) * p);  in tgamma_large()
117 static long double tgamma_tiny(long double x,  in tgamma_tiny()
118                                bool negative, long double negadjust)  in tgamma_tiny()
124     long double g = poly(coeffs_tiny, lenof(coeffs_tiny), 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()
133                                     long double negadjust)  in tgamma_ultratiny()
146 static long double tgamma_central(long double x)  in tgamma_central()
158     long double t = (x - min_x_hi) - min_x_lo;  in tgamma_central()
164     long double p;  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()
177      * out cases of non-normalized numbers.  in tgamma128()
197             return x-x; /* gamma(-inf) has indeterminate sign, so provoke an  in tgamma128()
202         exponent = -16384;  in tgamma128()
206         exponent--;  in tgamma128()
211     long double negadjust = 0.0L;  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()
237         x = -x;  in tgamma128()
250          * sub-case was going to do anyway.  in tgamma128()
269      * we do the negative-number adjustment the usual way, we'll leave  in tgamma128()
272     long double g;  in tgamma128()
280             long double huge = 0x1p+12288L;  in tgamma128()
283             long double tiny = 0x1p-12288L;  in tgamma128()
287         /* Negative-number adjustment happens inside here */  in tgamma128()
289     } else if (exponent < -113) {  in tgamma128()
290         /* Negative-number adjustment happens inside here */  in tgamma128()
292     } else if (exponent < -5) {  in tgamma128()
293         /* Negative-number adjustment happens inside here */  in tgamma128()
299          * For x in [1/32,1) we range-reduce upwards to the interval  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()
313          * (x-1)(x-2)...(x-k+1) as a polynomial in t.  in tgamma128()
315         long double mult;  in tgamma128()
318             mult = (x-1);  in tgamma128()
320             long double t = x - (i + 0.5L);  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()