Lines Matching +full:range +full:- +full:double
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)
36 long double result = coeffs[--n];
39 result = (result * x) + coeffs[--n];
46 * relates gamma(-x) and gamma(x).
48 static long double sin_pi_x_over_pi(long double x)
51 long double fracpart = remquol(x, 0.5L, &quo);
53 long double sign = 1.0L;
55 sign = -sign;
58 if (quo == 0 && fabsl(fracpart) < 0x1.p-58L) {
72 static long double tgamma_large(long double x,
73 bool negative, long double negadjust)
76 * In this range we compute gamma(x) as x^(x-1/2) * e^-x * K,
83 long double t = 1/x;
84 long double p = poly(coeffs_large, lenof(coeffs_large), t);
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
90 * of exp(-(x-0.5)).
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
96 long double powhalf = powl(x, (x-0.5L)/2.0L);
97 long double expret = expl(-(x-0.5L));
106 * where gamma(-x) doesn't underflow. Not only that, but the
107 * FP format has greater range in the tiny domain due to
111 long double ret = 1 / ((expret * powhalf) * (x * negadjust) * p);
117 static long double tgamma_tiny(long double x,
118 bool negative, long double negadjust)
124 long double g = poly(coeffs_tiny, lenof(coeffs_tiny), x);
131 /* Return tgamma(x) on the assumption that 0 <= x < 2^-113. */
132 static long double tgamma_ultratiny(long double x, bool negative,
133 long double negadjust)
146 static long double tgamma_central(long double x)
158 long double t = (x - min_x_hi) - min_x_lo;
164 long double p;
166 p = poly(coeffs_central_neg, lenof(coeffs_central_neg), -t);
173 long double tgamma128(long double x)
177 * out cases of non-normalized numbers.
197 return x-x; /* gamma(-inf) has indeterminate sign, so provoke an
202 exponent = -16384;
206 exponent--;
211 long double negadjust = 0.0L;
217 * gamma(1-x) gamma(x) = pi/sin(pi*x)
220 * => gamma(x) = --------------------
221 * gamma(1-x) sin(pi*x)
223 * But computing 1-x is going to lose a lot of accuracy when x
225 * gamma(t+1)=t gamma(t). Setting t=-x, this gives us
226 * gamma(1-x) = -x gamma(-x), so we now have
229 * gamma(x) = ----------------------
230 * -x gamma(-x) sin(pi*x)
232 * which relates gamma(x) to gamma(-x), which is much nicer,
233 * since x can be turned into -x without rounding.
237 x = -x;
250 * sub-case was going to do anyway.
269 * we do the negative-number adjustment the usual way, we'll leave
272 long double g;
280 long double huge = 0x1p+12288L;
283 long double tiny = 0x1p-12288L;
287 /* Negative-number adjustment happens inside here */
289 } else if (exponent < -113) {
290 /* Negative-number adjustment happens inside here */
292 } else if (exponent < -5) {
293 /* Negative-number adjustment happens inside here */
299 * For x in [1/32,1) we range-reduce upwards to the interval
306 * For x in [2,8) we range-reduce downwards to the interval
309 * Actually multiplying (x-1) by (x-2) by (x-3) and so on
313 * (x-1)(x-2)...(x-k+1) as a polynomial in t.
315 long double mult;
318 mult = (x-1);
320 long double t = x - (i + 0.5L);
323 * (x-1)(x-2) = (2.5+t)(1.5+t) = 3.75 + 4t + t^2 */
346 g = tgamma_central(x - (i-1)) * mult;