/* Implementation of gamma function according to ISO C.
- Copyright (C) 1997-2016 Free Software Foundation, Inc.
+ Copyright (C) 1997-2019 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
- <http://www.gnu.org/licenses/>. */
+ <https://www.gnu.org/licenses/>. */
#include <math.h>
+#include <math-narrow-eval.h>
#include <math_private.h>
+#include <fenv_private.h>
+#include <math-underflow.h>
#include <float.h>
/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
{
/* Adjust into the range for using exp (lgamma). */
*exp2_adj = 0;
- double n = __ceil (x - 1.5);
+ double n = ceil (x - 1.5);
double x_adj = x - n;
double eps;
double prod = __gamma_product (x_adj, 0, n, &eps);
{
/* Adjust into the range for applying Stirling's
approximation. */
- double n = __ceil (12.0 - x);
+ double n = ceil (12.0 - x);
x_adj = math_narrow_eval (x + n);
x_eps = (x - (x_adj - n));
prod = __gamma_product (x_adj - n, x_eps, n, &eps);
starting by computing pow (X_ADJ, X_ADJ) with a power of 2
factored out. */
double exp_adj = -eps;
- double x_adj_int = __round (x_adj);
+ double x_adj_int = round (x_adj);
double x_adj_frac = x_adj - x_adj_int;
int x_adj_log2;
double x_adj_mant = __frexp (x_adj, &x_adj_log2);
double ret = (__ieee754_pow (x_adj_mant, x_adj)
* __ieee754_exp2 (x_adj_log2 * x_adj_frac)
* __ieee754_exp (-x_adj)
- * __ieee754_sqrt (2 * M_PI / x_adj)
+ * sqrt (2 * M_PI / x_adj)
/ prod);
exp_adj += x_eps * __ieee754_log (x_adj);
double bsum = gamma_coeff[NCOEFF - 1];
__ieee754_gamma_r (double x, int *signgamp)
{
int32_t hx;
- u_int32_t lx;
+ uint32_t lx;
double ret;
EXTRACT_WORDS (hx, lx, x);
return 1.0 / x;
}
if (__builtin_expect (hx < 0, 0)
- && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
+ && (uint32_t) hx < 0xfff00000 && rint (x) == x)
{
/* Return value for integer x < 0 is NaN with invalid exception. */
*signgamp = 0;
}
else
{
- double tx = __trunc (x);
- *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1;
+ double tx = trunc (x);
+ *signgamp = (tx == 2.0 * trunc (tx / 2.0)) ? -1 : 1;
if (x <= -184.0)
/* Underflow. */
ret = DBL_MIN * DBL_MIN;
{
if (*signgamp < 0)
{
- ret = math_narrow_eval (-__copysign (DBL_MAX, ret) * DBL_MAX);
+ ret = math_narrow_eval (-copysign (DBL_MAX, ret) * DBL_MAX);
ret = -ret;
}
else
- ret = math_narrow_eval (__copysign (DBL_MAX, ret) * DBL_MAX);
+ ret = math_narrow_eval (copysign (DBL_MAX, ret) * DBL_MAX);
return ret;
}
else if (ret == 0)
{
if (*signgamp < 0)
{
- ret = math_narrow_eval (-__copysign (DBL_MIN, ret) * DBL_MIN);
+ ret = math_narrow_eval (-copysign (DBL_MIN, ret) * DBL_MIN);
ret = -ret;
}
else
- ret = math_narrow_eval (__copysign (DBL_MIN, ret) * DBL_MIN);
+ ret = math_narrow_eval (copysign (DBL_MIN, ret) * DBL_MIN);
return ret;
}
else