1 /*
2 * Implementation of the true gamma function (as opposed to lgamma)
3 * for 128-bit long double.
4 *
5 * Copyright (c) 2006-2024, Arm Limited.
6 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
7 */
8
9 /*
10 * This module implements the float128 gamma function under the name
11 * tgamma128. It's expected to be suitable for integration into system
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
14 * handling and optimize the initial process of extracting the
15 * exponent, which is done here by simple and portable (but
16 * potentially slower) methods.
17 */
18
19 #include <float.h>
20 #include <math.h>
21 #include <stdbool.h>
22 #include <stddef.h>
23
24 /* Only binary128 format is supported. */
25 #if LDBL_MANT_DIG == 113
26
27 #include "tgamma128.h"
28
29 #define lenof(x) (sizeof(x)/sizeof(*(x)))
30
31 /*
32 * Helper routine to evaluate a polynomial via Horner's rule
33 */
poly(const long double * coeffs,size_t n,long double x)34 static long double poly(const long double *coeffs, size_t n, long double x)
35 {
36 long double result = coeffs[--n];
37
38 while (n > 0)
39 result = (result * x) + coeffs[--n];
40
41 return result;
42 }
43
44 /*
45 * Compute sin(pi*x) / pi, for use in the reflection formula that
46 * relates gamma(-x) and gamma(x).
47 */
sin_pi_x_over_pi(long double x)48 static long double sin_pi_x_over_pi(long double x)
49 {
50 int quo;
51 long double fracpart = remquol(x, 0.5L, &quo);
52
53 long double sign = 1.0L;
54 if (quo & 2)
55 sign = -sign;
56 quo &= 1;
57
58 if (quo == 0 && fabsl(fracpart) < 0x1.p-58L) {
59 /* For numbers this size, sin(pi*x) is so close to pi*x that
60 * sin(pi*x)/pi is indistinguishable from x in float128 */
61 return sign * fracpart;
62 }
63
64 if (quo == 0) {
65 return sign * sinl(pi*fracpart) / pi;
66 } else {
67 return sign * cosl(pi*fracpart) / pi;
68 }
69 }
70
71 /* Return tgamma(x) on the assumption that x >= 8. */
tgamma_large(long double x,bool negative,long double negadjust)72 static long double tgamma_large(long double x,
73 bool negative, long double negadjust)
74 {
75 /*
76 * In this range we compute gamma(x) as x^(x-1/2) * e^-x * K,
77 * where K is a correction factor computed as a polynomial in 1/x.
78 *
79 * (Vaguely inspired by the form of the Lanczos approximation, but
80 * I tried the Lanczos approximation itself and it suffers badly
81 * from big cancellation leading to loss of significance.)
82 */
83 long double t = 1/x;
84 long double p = poly(coeffs_large, lenof(coeffs_large), t);
85
86 /*
87 * To avoid overflow in cases where x^(x-0.5) does overflow
88 * but gamma(x) does not, we split x^(x-0.5) in half and
89 * multiply back up _after_ multiplying the shrinking factor
90 * of exp(-(x-0.5)).
91 *
92 * Note that computing x-0.5 and (x-0.5)/2 is exact for the
93 * relevant range of x, so the only sources of error are pow
94 * and exp themselves, plus the multiplications.
95 */
96 long double powhalf = powl(x, (x-0.5L)/2.0L);
97 long double expret = expl(-(x-0.5L));
98
99 if (!negative) {
100 return (expret * powhalf) * powhalf * p;
101 } else {
102 /*
103 * Apply the reflection formula as commented below, but
104 * carefully: negadjust has magnitude less than 1, so it can
105 * turn a case where gamma(+x) would overflow into a case
106 * where gamma(-x) doesn't underflow. Not only that, but the
107 * FP format has greater range in the tiny domain due to
108 * denormals. For both reasons, it's not good enough to
109 * compute the positive result and then adjust it.
110 */
111 long double ret = 1 / ((expret * powhalf) * (x * negadjust) * p);
112 return ret / powhalf;
113 }
114 }
115
116 /* Return tgamma(x) on the assumption that 0 <= x < 1/32. */
tgamma_tiny(long double x,bool negative,long double negadjust)117 static long double tgamma_tiny(long double x,
118 bool negative, long double negadjust)
119 {
120 /*
121 * For x near zero, we use a polynomial approximation to
122 * g = 1/(x*gamma(x)), and then return 1/(g*x).
123 */
124 long double g = poly(coeffs_tiny, lenof(coeffs_tiny), x);
125 if (!negative)
126 return 1.0L / (g*x);
127 else
128 return g / negadjust;
129 }
130
131 /* Return tgamma(x) on the assumption that 0 <= x < 2^-113. */
tgamma_ultratiny(long double x,bool negative,long double negadjust)132 static long double tgamma_ultratiny(long double x, bool negative,
133 long double negadjust)
134 {
135 /* On this interval, gamma can't even be distinguished from 1/x,
136 * so we skip the polynomial evaluation in tgamma_tiny, partly to
137 * save time and partly to avoid the tiny intermediate values
138 * setting the underflow exception flag. */
139 if (!negative)
140 return 1.0L / x;
141 else
142 return 1.0L / negadjust;
143 }
144
145 /* Return tgamma(x) on the assumption that 1 <= x <= 2. */
tgamma_central(long double x)146 static long double tgamma_central(long double x)
147 {
148 /*
149 * In this central interval, our strategy is to finding the
150 * difference between x and the point where gamma has a minimum,
151 * and approximate based on that.
152 */
153
154 /* The difference between the input x and the minimum x. The first
155 * subtraction is expected to be exact, since x and min_hi have
156 * the same exponent (unless x=2, in which case it will still be
157 * exact). */
158 long double t = (x - min_x_hi) - min_x_lo;
159
160 /*
161 * Now use two different polynomials for the intervals [1,m] and
162 * [m,2].
163 */
164 long double p;
165 if (t < 0)
166 p = poly(coeffs_central_neg, lenof(coeffs_central_neg), -t);
167 else
168 p = poly(coeffs_central_pos, lenof(coeffs_central_pos), t);
169
170 return (min_y_lo + p * (t*t)) + min_y_hi;
171 }
172
tgamma128(long double x)173 long double tgamma128(long double x)
174 {
175 /*
176 * Start by extracting the number's sign and exponent, and ruling
177 * out cases of non-normalized numbers.
178 *
179 * For an implementation integrated into a system libm, it would
180 * almost certainly be quicker to do this by direct bitwise access
181 * to the input float128 value, using whatever is the local idiom
182 * for knowing its endianness.
183 *
184 * Integration into a system libc may also need to worry about
185 * setting errno, if that's the locally preferred way to report
186 * math.h errors.
187 */
188 int sign = signbit(x);
189 int exponent;
190 switch (fpclassify(x)) {
191 case FP_NAN:
192 return x+x; /* propagate QNaN, make SNaN throw an exception */
193 case FP_ZERO:
194 return 1/x; /* divide by zero on purpose to indicate a pole */
195 case FP_INFINITE:
196 if (sign) {
197 return x-x; /* gamma(-inf) has indeterminate sign, so provoke an
198 * IEEE invalid operation exception to indicate that */
199 }
200 return x; /* but gamma(+inf) is just +inf with no error */
201 case FP_SUBNORMAL:
202 exponent = -16384;
203 break;
204 default:
205 frexpl(x, &exponent);
206 exponent--;
207 break;
208 }
209
210 bool negative = false;
211 long double negadjust = 0.0L;
212
213 if (sign) {
214 /*
215 * Euler's reflection formula is
216 *
217 * gamma(1-x) gamma(x) = pi/sin(pi*x)
218 *
219 * pi
220 * => gamma(x) = --------------------
221 * gamma(1-x) sin(pi*x)
222 *
223 * But computing 1-x is going to lose a lot of accuracy when x
224 * is very small, so instead we transform using the recurrence
225 * gamma(t+1)=t gamma(t). Setting t=-x, this gives us
226 * gamma(1-x) = -x gamma(-x), so we now have
227 *
228 * pi
229 * gamma(x) = ----------------------
230 * -x gamma(-x) sin(pi*x)
231 *
232 * which relates gamma(x) to gamma(-x), which is much nicer,
233 * since x can be turned into -x without rounding.
234 */
235 negadjust = sin_pi_x_over_pi(x);
236 negative = true;
237 x = -x;
238
239 /*
240 * Now the ultimate answer we want is
241 *
242 * 1 / (gamma(x) * x * negadjust)
243 *
244 * where x is the positive value we've just turned it into.
245 *
246 * For some of the cases below, we'll compute gamma(x)
247 * normally and then compute this adjusted value afterwards.
248 * But for others, we can implement the reciprocal operation
249 * in this formula by _avoiding_ an inversion that the
250 * sub-case was going to do anyway.
251 */
252
253 if (negadjust == 0) {
254 /*
255 * Special case for negative integers. Applying the
256 * reflection formula would cause division by zero, but
257 * standards would prefer we treat this error case as an
258 * invalid operation and return NaN instead. (Possibly
259 * because otherwise you'd have to decide which sign of
260 * infinity to return, and unlike the x=0 case, there's no
261 * sign of zero available to disambiguate.)
262 */
263 return negadjust / negadjust;
264 }
265 }
266
267 /*
268 * Split the positive domain into various cases. For cases where
269 * we do the negative-number adjustment the usual way, we'll leave
270 * the answer in 'g' and drop out of the if statement.
271 */
272 long double g;
273
274 if (exponent >= 11) {
275 /*
276 * gamma of any positive value this large overflows, and gamma
277 * of any negative value underflows.
278 */
279 if (!negative) {
280 long double huge = 0x1p+12288L;
281 return huge * huge; /* provoke an overflow */
282 } else {
283 long double tiny = 0x1p-12288L;
284 return tiny * tiny * negadjust; /* underflow, of the right sign */
285 }
286 } else if (exponent >= 3) {
287 /* Negative-number adjustment happens inside here */
288 return tgamma_large(x, negative, negadjust);
289 } else if (exponent < -113) {
290 /* Negative-number adjustment happens inside here */
291 return tgamma_ultratiny(x, negative, negadjust);
292 } else if (exponent < -5) {
293 /* Negative-number adjustment happens inside here */
294 return tgamma_tiny(x, negative, negadjust);
295 } else if (exponent == 0) {
296 g = tgamma_central(x);
297 } else if (exponent < 0) {
298 /*
299 * For x in [1/32,1) we range-reduce upwards to the interval
300 * [1,2), using the inverse of the normal recurrence formula:
301 * gamma(x) = gamma(x+1)/x.
302 */
303 g = tgamma_central(1+x) / x;
304 } else {
305 /*
306 * For x in [2,8) we range-reduce downwards to the interval
307 * [1,2) by repeated application of the recurrence formula.
308 *
309 * Actually multiplying (x-1) by (x-2) by (x-3) and so on
310 * would introduce multiple ULPs of rounding error. We can get
311 * better accuracy by writing x = (k+1/2) + t, where k is an
312 * integer and |t|<1/2, and expanding out the obvious factor
313 * (x-1)(x-2)...(x-k+1) as a polynomial in t.
314 */
315 long double mult;
316 int i = x;
317 if (i == 2) { /* x in [2,3) */
318 mult = (x-1);
319 } else {
320 long double t = x - (i + 0.5L);
321 switch (i) {
322 /* E.g. for x=3.5+t, we want
323 * (x-1)(x-2) = (2.5+t)(1.5+t) = 3.75 + 4t + t^2 */
324 case 3:
325 mult = 3.75L+t*(4.0L+t);
326 break;
327 case 4:
328 mult = 13.125L+t*(17.75L+t*(7.5L+t));
329 break;
330 case 5:
331 mult = 59.0625L+t*(93.0L+t*(51.50L+t*(12.0L+t)));
332 break;
333 case 6:
334 mult = 324.84375L+t*(570.5625L+t*(376.250L+t*(
335 117.5L+t*(17.5L+t))));
336 break;
337 case 7:
338 mult = 2111.484375L+t*(4033.5L+t*(3016.1875L+t*(
339 1140.0L+t*(231.25L+t*(24.0L+t)))));
340 break;
341 }
342 }
343
344 g = tgamma_central(x - (i-1)) * mult;
345 }
346
347 if (!negative) {
348 /* Positive domain: return g unmodified */
349 return g;
350 } else {
351 /* Negative domain: apply the reflection formula as commented above */
352 return 1.0L / (g * x * negadjust);
353 }
354 }
355
356 #endif
357