diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-96')
47 files changed, 429 insertions, 1124 deletions
diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c index c33701f11e..c1ffa3e0d9 100644 --- a/sysdeps/ieee754/ldbl-96/e_asinl.c +++ b/sysdeps/ieee754/ldbl-96/e_asinl.c @@ -64,9 +64,12 @@ static const long double one = 1.0L, huge = 1.0e+4932L, - pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L, - pio2_lo = 2.9127320560933561582586004641843300502121E-20L, - pio4_hi = 7.8539816339744830960109718553940894025800E-1L, + pio2_hi = 0x1.921fb54442d1846ap+0L, /* pi/2 rounded to nearest to 64 + bits. */ + pio2_lo = -0x7.6733ae8fe47c65d8p-68L, /* pi/2 - pio2_hi rounded to + nearest to 64 bits. */ + pio4_hi = 0xc.90fdaa22168c235p-4L, /* pi/4 rounded to nearest to 64 + bits. */ /* coefficient for R(x^2) */ diff --git a/sysdeps/ieee754/ldbl-96/e_atan2l.c b/sysdeps/ieee754/ldbl-96/e_atan2l.c deleted file mode 100644 index 209f29bbd2..0000000000 --- a/sysdeps/ieee754/ldbl-96/e_atan2l.c +++ /dev/null @@ -1,125 +0,0 @@ -/* e_atan2l.c -- long double version of e_atan2.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* __ieee754_atan2l(y,x) - * Method : - * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): - * ARG (x+iy) = arctan(y/x) ... if x > 0, - * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * - * Special cases: - * - * ATAN2((anything), NaN ) is NaN; - * ATAN2(NAN , (anything) ) is NaN; - * ATAN2(+-0, +(anything but NaN)) is +-0 ; - * ATAN2(+-0, -(anything but NaN)) is +-pi ; - * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; - * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; - * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; - * ATAN2(+-INF,+INF ) is +-pi/4 ; - * ATAN2(+-INF,-INF ) is +-3pi/4; - * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <math.h> -#include <math_private.h> - -static const long double -tiny = 1.0e-4900L, -zero = 0.0, -pi_o_4 = 7.85398163397448309628202E-01L, /* 0x3FFE, 0xC90FDAA2, 0x2168C235 */ -pi_o_2 = 1.5707963267948966192564E+00L, /* 0x3FFF, 0xC90FDAA2, 0x2168C235 */ -pi = 3.14159265358979323851281E+00L, /* 0x4000, 0xC90FDAA2, 0x2168C235 */ -pi_lo = -5.01655761266833202345176e-20L;/* 0xBFBE, 0xECE675D1, 0xFC8F8CBB */ - -long double -__ieee754_atan2l (long double y, long double x) -{ - long double z; - int32_t k,m,hx,hy,ix,iy; - u_int32_t sx,sy,lx,ly; - - GET_LDOUBLE_WORDS(sx,hx,lx,x); - ix = sx&0x7fff; - lx |= hx & 0x7fffffff; - GET_LDOUBLE_WORDS(sy,hy,ly,y); - iy = sy&0x7fff; - ly |= hy & 0x7fffffff; - if(((2*ix|((lx|-lx)>>31))>0xfffe)|| - ((2*iy|((ly|-ly)>>31))>0xfffe)) /* x or y is NaN */ - return x+y; - if(((sx-0x3fff)|lx)==0) return __atanl(y); /* x=1.0 */ - m = ((sy>>15)&1)|((sx>>14)&2); /* 2*sign(x)+sign(y) */ - - /* when y = 0 */ - if((iy|ly)==0) { - switch(m) { - case 0: - case 1: return y; /* atan(+-0,+anything)=+-0 */ - case 2: return pi+tiny;/* atan(+0,-anything) = pi */ - case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ - } - } - /* when x = 0 */ - if((ix|lx)==0) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny; - - /* when x is INF */ - if(ix==0x7fff) { - if(iy==0x7fff) { - switch(m) { - case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ - case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ - case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ - case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ - } - } else { - switch(m) { - case 0: return zero ; /* atan(+...,+INF) */ - case 1: return -zero ; /* atan(-...,+INF) */ - case 2: return pi+tiny ; /* atan(+...,-INF) */ - case 3: return -pi-tiny ; /* atan(-...,-INF) */ - } - } - } - /* when y is INF */ - if(iy==0x7fff) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny; - - /* compute y/x */ - k = sy-sx; - if(k > 70) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**70 */ - else if(sx>=0x8000&&k<-70) z=0.0; /* |y|/x < -2**70 */ - else z=__atanl(fabsl(y/x)); /* safe to do y/x */ - switch (m) { - case 0: return z ; /* atan(+,+) */ - case 1: { - u_int32_t sz; - GET_LDOUBLE_EXP(sz,z); - SET_LDOUBLE_EXP(z,sz ^ 0x8000); - } - return z ; /* atan(-,+) */ - case 2: return pi-(z-pi_lo);/* atan(+,-) */ - default: /* case 3 */ - return (z-pi_lo)-pi;/* atan(-,-) */ - } -} -strong_alias (__ieee754_atan2l, __atan2l_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c index 0885c8fa8c..477c3a61d3 100644 --- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997,1999,2001,2003,2004,2011 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -19,14 +19,102 @@ #include <math.h> #include <math_private.h> +#include <float.h> +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's + approximation to gamma function. */ + +static const long double gamma_coeff[] = + { + 0x1.5555555555555556p-4L, + -0xb.60b60b60b60b60bp-12L, + 0x3.4034034034034034p-12L, + -0x2.7027027027027028p-12L, + 0x3.72a3c5631fe46aep-12L, + -0x7.daac36664f1f208p-12L, + 0x1.a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708p-8L, + }; + +#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) + +/* Return gamma (X), for positive X less than 1766, in the form R * + 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to + avoid overflow or underflow in intermediate calculations. */ + +static long double +gammal_positive (long double x, int *exp2_adj) +{ + int local_signgam; + if (x < 0.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; + } + else if (x <= 1.5L) + { + *exp2_adj = 0; + return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); + } + else if (x < 7.5L) + { + /* Adjust into the range for using exp (lgamma). */ + *exp2_adj = 0; + long double n = __ceill (x - 1.5L); + long double x_adj = x - n; + long double eps; + long double prod = __gamma_productl (x_adj, 0, n, &eps); + return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) + * prod * (1.0L + eps)); + } + else + { + long double eps = 0; + long double x_eps = 0; + long double x_adj = x; + long double prod = 1; + if (x < 13.0L) + { + /* Adjust into the range for applying Stirling's + approximation. */ + long double n = __ceill (13.0L - x); + x_adj = x + n; + x_eps = (x - (x_adj - n)); + prod = __gamma_productl (x_adj - n, x_eps, n, &eps); + } + /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). + Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, + starting by computing pow (X_ADJ, X_ADJ) with a power of 2 + factored out. */ + long double exp_adj = -eps; + long double x_adj_int = __roundl (x_adj); + long double x_adj_frac = x_adj - x_adj_int; + int x_adj_log2; + long double x_adj_mant = __frexpl (x_adj, &x_adj_log2); + if (x_adj_mant < M_SQRT1_2l) + { + x_adj_log2--; + x_adj_mant *= 2.0L; + } + *exp2_adj = x_adj_log2 * (int) x_adj_int; + long double ret = (__ieee754_powl (x_adj_mant, x_adj) + * __ieee754_exp2l (x_adj_log2 * x_adj_frac) + * __ieee754_expl (-x_adj) + * __ieee754_sqrtl (2 * M_PIl / x_adj) + / prod); + exp_adj += x_eps * __ieee754_logl (x); + long double bsum = gamma_coeff[NCOEFF - 1]; + long double x_adj2 = x_adj * x_adj; + for (size_t i = 1; i <= NCOEFF - 1; i++) + bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; + exp_adj += bsum / x_adj; + return ret + ret * __expm1l (exp_adj); + } +} long double __ieee754_gammal_r (long double x, int *signgamp) { - /* We don't have a real gamma implementation now. We'll use lgamma - and the exp function. But due to the required boundary - conditions we must check some values separately. */ u_int32_t es, hx, lx; GET_LDOUBLE_WORDS (es, hx, lx, x); @@ -43,22 +131,55 @@ __ieee754_gammal_r (long double x, int *signgamp) *signgamp = 0; return x - x; } - if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0) - && ((hx & 0x7fffffff) | lx) != 0) + if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0)) { - /* NaN, return it. */ + /* Positive infinity (return positive infinity) or NaN (return + NaN). */ *signgamp = 0; - return x; + return x + x; } - if (__builtin_expect ((es & 0x8000) != 0, 0) - && x < 0xffffffff && __rintl (x) == x) + if (__builtin_expect ((es & 0x8000) != 0, 0) && __rintl (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; return (x - x) / (x - x); } - /* XXX FIXME. */ - return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); + if (x >= 1756.0L) + { + /* Overflow. */ + *signgamp = 0; + return LDBL_MAX * LDBL_MAX; + } + else if (x > 0.0L) + { + *signgamp = 0; + int exp2_adj; + long double ret = gammal_positive (x, &exp2_adj); + return __scalbnl (ret, exp2_adj); + } + else if (x >= -LDBL_EPSILON / 4.0L) + { + *signgamp = 0; + return 1.0f / x; + } + else + { + long double tx = __truncl (x); + *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1; + if (x <= -1766.0L) + /* Underflow. */ + return LDBL_MIN * LDBL_MIN; + long double frac = tx - x; + if (frac > 0.5L) + frac = 1.0L - frac; + long double sinpix = (frac <= 0.25L + ? __sinl (M_PIl * frac) + : __cosl (M_PIl * (0.5L - frac))); + int exp2_adj; + long double ret = M_PIl / (-x * sinpix + * gammal_positive (-x, &exp2_adj)); + return __scalbnl (ret, -exp2_adj); + } } strong_alias (__ieee754_gammal_r, __gammal_r_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c index 306f92924c..d3152f91e5 100644 --- a/sysdeps/ieee754/ldbl-96/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c @@ -85,10 +85,21 @@ long double __ieee754_hypotl(long double x, long double y) u_int32_t high,low; GET_LDOUBLE_WORDS(exp,high,low,b); if((high|low)==0) return a; - SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0); /* t1=2^16382 */ + SET_LDOUBLE_WORDS(t1, 0x7ffd, 0x80000000, 0); /* t1=2^16382 */ b *= t1; a *= t1; k -= 16382; + GET_LDOUBLE_EXP (ea, a); + GET_LDOUBLE_EXP (eb, b); + if (eb > ea) + { + t1 = a; + a = b; + b = t1; + j = ea; + ea = eb; + eb = j; + } } else { /* scale a and b by 2^9600 */ ea += 0x2580; /* a *= 2^9600 */ eb += 0x2580; /* b *= 2^9600 */ diff --git a/sysdeps/ieee754/ldbl-96/e_ilogbl.c b/sysdeps/ieee754/ldbl-96/e_ilogbl.c deleted file mode 100644 index 0c7d9d5440..0000000000 --- a/sysdeps/ieee754/ldbl-96/e_ilogbl.c +++ /dev/null @@ -1,59 +0,0 @@ -/* s_ilogbl.c -- long double version of s_ilogb.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* ilogbl(long double x) - * return the binary exponent of non-zero x - * ilogbl(0) = FP_ILOGB0 - * ilogbl(NaN) = FP_ILOGBNAN (no signal is raised) - * ilogbl(+-Inf) = INT_MAX (no signal is raised) - */ - -#include <limits.h> -#include <math.h> -#include <math_private.h> - -int __ieee754_ilogbl (long double x) -{ - int32_t es,hx,lx,ix; - - GET_LDOUBLE_EXP(es,x); - es &= 0x7fff; - if(es==0) { - GET_LDOUBLE_WORDS(es,hx,lx,x); - if((hx|lx)==0) - return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */ - else /* subnormal x */ - if(hx==0) { - for (ix = -16415; lx>0; lx<<=1) ix -=1; - } else { - for (ix = -16383; hx>0; hx<<=1) ix -=1; - } - return ix; - } - else if (es<0x7fff) return es-0x3fff; - else if (FP_ILOGBNAN != INT_MAX) - { - GET_LDOUBLE_WORDS(es,hx,lx,x); - if (((hx & 0x7fffffff)|lx) == 0) - /* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */ - return INT_MAX; - } - return FP_ILOGBNAN; -} diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c index 785c0b0676..4c13018aea 100644 --- a/sysdeps/ieee754/ldbl-96/e_j1l.c +++ b/sysdeps/ieee754/ldbl-96/e_j1l.c @@ -203,7 +203,7 @@ __ieee754_y1l (long double x) __sincosl (x, &s, &c); ss = -s - c; cc = s - c; - if (ix < 0x7fe00000) + if (ix < 0x7ffe) { /* make sure x+x not overflow */ z = __cosl (x + x); if ((s * c) > zero) diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c index 58a9107f7d..fa8e27efec 100644 --- a/sysdeps/ieee754/ldbl-96/e_jnl.c +++ b/sysdeps/ieee754/ldbl-96/e_jnl.c @@ -302,7 +302,8 @@ __ieee754_ynl (int n, long double x) if (__builtin_expect ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0), 0)) return x + x; if (__builtin_expect ((ix | i0 | i1) == 0, 0)) - return -HUGE_VALL + x; /* -inf and overflow exception. */ + /* -inf or inf and divide-by-zero exception. */ + return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L; if (__builtin_expect (se & 0x8000, 0)) return zero / (zero * x); sign = 1; diff --git a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c index b72230962d..e18be6ee0c 100644 --- a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c @@ -1,5 +1,5 @@ /* Extended-precision floating point argument reduction. - Copyright (C) 1999-2012 Free Software Foundation, Inc. + Copyright (C) 1999-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision code by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/e_remainderl.c b/sysdeps/ieee754/ldbl-96/e_remainderl.c deleted file mode 100644 index 290e483ae5..0000000000 --- a/sysdeps/ieee754/ldbl-96/e_remainderl.c +++ /dev/null @@ -1,72 +0,0 @@ -/* e_remainderl.c -- long double version of e_remainder.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* __ieee754_remainderl(x,p) - * Return : - * returns x REM p = x - [x/p]*p as if in infinite - * precise arithmetic, where [x/p] is the (infinite bit) - * integer nearest x/p (in half way case choose the even one). - * Method : - * Based on fmod() return x-[x/p]chopped*p exactlp. - */ - -#include <math.h> -#include <math_private.h> - -static const long double zero = 0.0; - - -long double -__ieee754_remainderl(long double x, long double p) -{ - u_int32_t sx,sex,sep,x0,x1,p0,p1; - long double p_half; - - GET_LDOUBLE_WORDS(sex,x0,x1,x); - GET_LDOUBLE_WORDS(sep,p0,p1,p); - sx = sex&0x8000; - sep &= 0x7fff; - sex &= 0x7fff; - - /* purge off exception values */ - if((sep|p0|p1)==0) return (x*p)/(x*p); /* p = 0 */ - if((sex==0x7fff)|| /* x not finite */ - ((sep==0x7fff)&& /* p is NaN */ - ((p0|p1)!=0))) - return (x*p)/(x*p); - - - if (sep<0x7ffe) x = __ieee754_fmodl(x,p+p); /* now x < 2p */ - if (((sex-sep)|(x0-p0)|(x1-p1))==0) return zero*x; - x = fabsl(x); - p = fabsl(p); - if (sep<0x0002) { - if(x+x>p) { - x-=p; - if(x+x>=p) x -= p; - } - } else { - p_half = 0.5*p; - if(x>p_half) { - x-=p; - if(x>=p_half) x -= p; - } - } - GET_LDOUBLE_EXP(sex,x); - SET_LDOUBLE_EXP(x,sex^sx); - return x; -} -strong_alias (__ieee754_remainderl, __remainderl_finite) diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c new file mode 100644 index 0000000000..214745624f --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/gamma_product.c @@ -0,0 +1,46 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013-2014 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + 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/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +double +__gamma_product (double x, double x_eps, int n, double *eps) +{ + long double x_full = (long double) x + (long double) x_eps; + long double ret = x_full; + for (int i = 1; i < n; i++) + ret *= x_full + i; + +#if FLT_EVAL_METHOD != 0 + volatile +#endif + double fret = ret; + *eps = (ret - fret) / fret; + + return fret; +} diff --git a/sysdeps/ieee754/ldbl-96/gamma_productl.c b/sysdeps/ieee754/ldbl-96/gamma_productl.c new file mode 100644 index 0000000000..2c87649677 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/gamma_productl.c @@ -0,0 +1,75 @@ +/* Compute a product of X, X+1, ..., with an error estimate. + Copyright (C) 2013-2014 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + 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/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static inline void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N + - 1, in the form R * (1 + *EPS) where the return value R is an + approximation to the product and *EPS is set to indicate the + approximate error in the return value. X is such that all the + values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / + X is small enough that factors quadratic in it can be + neglected. */ + +long double +__gamma_productl (long double x, long double x_eps, int n, long double *eps) +{ + SET_RESTORE_ROUNDL (FE_TONEAREST); + long double ret = x; + *eps = x_eps / x; + for (int i = 1; i < n; i++) + { + *eps += x_eps / (x + i); + long double lo; + mul_split (&ret, &lo, ret, x + i); + *eps += lo / ret; + } + return ret; +} diff --git a/sysdeps/ieee754/ldbl-96/k_cosl.c b/sysdeps/ieee754/ldbl-96/k_cosl.c index 9e8f33a283..daf7c060d2 100644 --- a/sysdeps/ieee754/ldbl-96/k_cosl.c +++ b/sysdeps/ieee754/ldbl-96/k_cosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point cosine on <-pi/4,pi/4>. - Copyright (C) 1999-2012 Free Software Foundation, Inc. + Copyright (C) 1999-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision cosine by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c index feb24d9e79..f2d1c860e8 100644 --- a/sysdeps/ieee754/ldbl-96/k_sinl.c +++ b/sysdeps/ieee754/ldbl-96/k_sinl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine on <-pi/4,pi/4>. - Copyright (C) 1999-2012 Free Software Foundation, Inc. + Copyright (C) 1999-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision sine by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c index 7292b7463c..be8b93ffa5 100644 --- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995,1996,1997,1998,2002,2003 Free Software Foundation, Inc. +/* Copyright (C) 1995-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -62,7 +62,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, else { /* It is a denormal number, meaning it has no implicit leading - one bit, and its exponent is in fact the format minimum. */ + one bit, and its exponent is in fact the format minimum. */ int cnt; if (res_ptr[N - 1] != 0) diff --git a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c index 2d5e126021..8fdb335fe4 100644 --- a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995,1996,1997,1998,2002,2003 Free Software Foundation, Inc. +/* Copyright (C) 1995-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/printf_fphex.c b/sysdeps/ieee754/ldbl-96/printf_fphex.c index acb0508b87..798c788485 100644 --- a/sysdeps/ieee754/ldbl-96/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-96/printf_fphex.c @@ -1,5 +1,5 @@ /* Print floating point number in hexadecimal notation according to ISO C99. - Copyright (C) 1997-2012 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -25,11 +25,13 @@ do { \ /* The "strange" 80 bit format on ix86 and m68k has an explicit \ leading digit in the 64 bit mantissa. */ \ unsigned long long int num; \ + union ieee854_long_double u; \ + u.d = fpnum.ldbl; \ \ assert (sizeof (long double) == 12); \ \ - num = (((unsigned long long int) fpnum.ldbl.ieee.mantissa0) << 32 \ - | fpnum.ldbl.ieee.mantissa1); \ + num = (((unsigned long long int) u.ieee.mantissa0) << 32 \ + | u.ieee.mantissa1); \ \ zero_mantissa = num == 0; \ \ @@ -62,7 +64,7 @@ do { \ \ /* We have 3 bits from the mantissa in the leading nibble. \ Therefore we are here using `IEEE854_LONG_DOUBLE_BIAS + 3'. */ \ - exponent = fpnum.ldbl.ieee.exponent; \ + exponent = u.ieee.exponent; \ \ if (exponent == 0) \ { \ diff --git a/sysdeps/ieee754/ldbl-96/s_cbrtl.c b/sysdeps/ieee754/ldbl-96/s_cbrtl.c index 07236345b9..49a6891642 100644 --- a/sysdeps/ieee754/ldbl-96/s_cbrtl.c +++ b/sysdeps/ieee754/ldbl-96/s_cbrtl.c @@ -1,5 +1,5 @@ /* Compute cubic root of double value. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Dirk Alboth <dirka@uni-paderborn.de> and Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -45,7 +45,7 @@ __cbrtl (long double x) int xe; /* Reduce X. XM now is an range 1.0 to 0.5. */ - xm = __frexpl (fabs (x), &xe); + xm = __frexpl (fabsl (x), &xe); /* If X is not finite or is null return it (with raising exceptions if necessary. diff --git a/sysdeps/ieee754/ldbl-96/s_ceill.c b/sysdeps/ieee754/ldbl-96/s_ceill.c deleted file mode 100644 index aef8a32f63..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_ceill.c +++ /dev/null @@ -1,85 +0,0 @@ -/* s_ceill.c -- long double version of s_ceil.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * ceill(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). - */ - -#include <math.h> -#include <math_private.h> - -static const long double huge = 1.0e4930; - -long double __ceill(long double x) -{ - int32_t i1,j0; - u_int32_t i,j,se,i0,sx; - GET_LDOUBLE_WORDS(se,i0,i1,x); - sx = (se>>15)&1; - j0 = (se&0x7fff)-0x3fff; - if(j0<31) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(sx) {se=0x8000;i0=0;i1=0;} - else if((i0|i1)!=0) { se=0x3fff;i0=0;i1=0;} - } - } else { - i = (0x7fffffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(sx==0) { - if (j0>0 && (i0+(0x80000000>>j0))>i0) - i0+=0x80000000>>j0; - else - { - i = 0x7fffffff; - ++se; - } - } - i0 &= (~i); i1=0; - } - } - } else if (j0>62) { - if(j0==0x4000) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-31); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(sx==0) { - if(j0==31) i0+=1; - else { - j = i1 + (1<<(63-j0)); - if(j<i1) i0+=1; /* got a carry */ - i1 = j; - } - } - i1 &= (~i); - } - } - SET_LDOUBLE_WORDS(x,se,i0,i1); - return x; -} -weak_alias (__ceill, ceill) diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c index b49a49be98..47e4b9e909 100644 --- a/sysdeps/ieee754/ldbl-96/s_erfl.c +++ b/sysdeps/ieee754/ldbl-96/s_erfl.c @@ -11,9 +11,9 @@ /* Long double expansions are Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> - and are incorporated herein by permission of the author. The author + and are incorporated herein by permission of the author. The author reserves the right to distribute this material elsewhere under different - copying permissions. These modifications are distributed here under + copying permissions. These modifications are distributed here under the following terms: This library is free software; you can redistribute it and/or @@ -104,6 +104,7 @@ */ +#include <errno.h> #include <math.h> #include <math_private.h> @@ -422,14 +423,22 @@ __erfcl (long double x) r = __ieee754_expl (-z * z - 0.5625) * __ieee754_expl ((z - x) * (z + x) + R / S); if ((se & 0x8000) == 0) - return r / x; + { + long double ret = r / x; + if (ret == 0) + __set_errno (ERANGE); + return ret; + } else return two - r / x; } else { if ((se & 0x8000) == 0) - return tiny * tiny; + { + __set_errno (ERANGE); + return tiny * tiny; + } else return two - tiny; } diff --git a/sysdeps/ieee754/ldbl-96/s_fabsl.c b/sysdeps/ieee754/ldbl-96/s_fabsl.c deleted file mode 100644 index fdc70e0dcd..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_fabsl.c +++ /dev/null @@ -1,35 +0,0 @@ -/* s_fabsl.c -- long double version of s_fabs.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * fabsl(x) returns the absolute value of x. - */ - -#include <math.h> -#include <math_private.h> - -long double __fabsl(long double x) -{ - u_int32_t exp; - GET_LDOUBLE_EXP(exp,x); - SET_LDOUBLE_EXP(x,exp&0x7fff); - return x; -} -weak_alias (__fabsl, fabsl) diff --git a/sysdeps/ieee754/ldbl-96/s_finitel.c b/sysdeps/ieee754/ldbl-96/s_finitel.c deleted file mode 100644 index fbf4cc691c..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_finitel.c +++ /dev/null @@ -1,36 +0,0 @@ -/* s_finitel.c -- long double version of s_finite.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * finitel(x) returns 1 is x is finite, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -int __finitel(long double x) -{ - int32_t exp; - GET_LDOUBLE_EXP(exp,x); - return (int)((u_int32_t)((exp&0x7fff)-0x7fff)>>31); -} -hidden_def (__finitel) -weak_alias (__finitel, finitel) diff --git a/sysdeps/ieee754/ldbl-96/s_floorl.c b/sysdeps/ieee754/ldbl-96/s_floorl.c deleted file mode 100644 index cad7935b33..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_floorl.c +++ /dev/null @@ -1,86 +0,0 @@ -/* s_floorl.c -- long double version of s_floor.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * floorl(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). - */ - -#include <math.h> -#include <math_private.h> - -static const long double huge = 1.0e4930; - -long double __floorl(long double x) -{ - int32_t i1,j0; - u_int32_t i,j,se,i0,sx; - GET_LDOUBLE_WORDS(se,i0,i1,x); - sx = (se>>15)&1; - j0 = (se&0x7fff)-0x3fff; - if(j0<31) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(sx==0) {se=0;i0=i1=0;} - else if(((se&0x7fff)|i0|i1)!=0) - { se=0xbfff;i0=i1=0;} - } - } else { - i = (0x7fffffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(sx) { - if (j0>0 && (i0+(0x80000000>>j0))>i0) - i0 += (0x80000000)>>j0; - else - { - i = 0x7fffffff; - ++se; - } - } - i0 &= (~i); i1=0; - } - } - } else if (j0>62) { - if(j0==0x4000) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-31); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(sx) { - if(j0==31) i0+=1; - else { - j = i1+(1<<(63-j0)); - if(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } - } - i1 &= (~i); - } - } - SET_LDOUBLE_WORDS(x,se,i0,i1); - return x; -} -weak_alias (__floorl, floorl) diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c index 001d8063d4..fde2811040 100644 --- a/sysdeps/ieee754/ldbl-96/s_fma.c +++ b/sysdeps/ieee754/ldbl-96/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2012 Free Software Foundation, Inc. + Copyright (C) 2010-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -42,6 +42,10 @@ __fma (double x, double y, double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + feholdexcept (&env); + fesetround (FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1) long double x1 = (long double) x * C; @@ -60,9 +64,19 @@ __fma (double x, double y, double z) t1 = m1 - t1; t2 = z - t2; long double a2 = t1 + t2; + feclearexcept (FE_INEXACT); + + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } - fenv_t env; - feholdexcept (&env); fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ a2 = a2 + m2; diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c index b5fdcba6bc..0564321354 100644 --- a/sysdeps/ieee754/ldbl-96/s_fmal.c +++ b/sysdeps/ieee754/ldbl-96/s_fmal.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2012 Free Software Foundation, Inc. + Copyright (C) 2010-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -56,16 +56,18 @@ __fmal (long double x, long double y, long double z) underflows to 0. */ if (z == 0 && x != 0 && y != 0) return x * y; - /* If x or y or z is Inf/NaN, or if fma will certainly overflow, - or if x * y is zero, compute as x * y + z. */ + /* If x or y or z is Inf/NaN, or if x * y is zero, compute as + x * y + z. */ if (u.ieee.exponent == 0x7fff || v.ieee.exponent == 0x7fff || w.ieee.exponent == 0x7fff - || u.ieee.exponent + v.ieee.exponent - > 0x7fff + IEEE854_LONG_DOUBLE_BIAS || x == 0 || y == 0) return x * y + z; + /* If fma will certainly overflow, compute as x * y. */ + if (u.ieee.exponent + v.ieee.exponent + > 0x7fff + IEEE854_LONG_DOUBLE_BIAS) + return x * y; /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the result nor whether there is underflow depends on its exact value, only on its sign. */ @@ -114,8 +116,17 @@ __fmal (long double x, long double y, long double z) { /* Similarly. If z exponent is very large and x and y exponents are - very small, it doesn't matter if we don't adjust it. */ - if (u.ieee.exponent > v.ieee.exponent) + very small, adjust them up to avoid spurious underflows, + rather than down. */ + if (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) + { + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + else + v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + } + else if (u.ieee.exponent > v.ieee.exponent) { if (u.ieee.exponent > LDBL_MANT_DIG) u.ieee.exponent -= LDBL_MANT_DIG; @@ -145,15 +156,15 @@ __fmal (long double x, long double y, long double z) <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */ { if (u.ieee.exponent > v.ieee.exponent) - u.ieee.exponent += 2 * LDBL_MANT_DIG; + u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; else - v.ieee.exponent += 2 * LDBL_MANT_DIG; - if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4) + v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; + if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 6) { if (w.ieee.exponent) - w.ieee.exponent += 2 * LDBL_MANT_DIG; + w.ieee.exponent += 2 * LDBL_MANT_DIG + 2; else - w.d *= 0x1p128L; + w.d *= 0x1p130L; adjust = -1; } /* Otherwise x * y should just affect inexact @@ -168,6 +179,10 @@ __fmal (long double x, long double y, long double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + feholdexcept (&env); + fesetround (FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) long double x1 = x * C; @@ -186,9 +201,19 @@ __fmal (long double x, long double y, long double z) t1 = m1 - t1; t2 = z - t2; long double a2 = t1 + t2; + feclearexcept (FE_INEXACT); + + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } - fenv_t env; - feholdexcept (&env); fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; @@ -224,19 +249,19 @@ __fmal (long double x, long double y, long double z) /* If a1 + u.d is exact, the only rounding happens during scaling down. */ if (j == 0) - return v.d * 0x1p-128L; + return v.d * 0x1p-130L; /* If result rounded to zero is not subnormal, no double rounding will occur. */ - if (v.ieee.exponent > 128) - return (a1 + u.d) * 0x1p-128L; - /* If v.d * 0x1p-128L with round to zero is a subnormal above - or equal to LDBL_MIN / 2, then v.d * 0x1p-128L shifts mantissa + if (v.ieee.exponent > 130) + return (a1 + u.d) * 0x1p-130L; + /* If v.d * 0x1p-130L with round to zero is a subnormal above + or equal to LDBL_MIN / 2, then v.d * 0x1p-130L shifts mantissa down just by 1 bit, which means v.ieee.mantissa1 |= j would change the round bit, not sticky or guard bit. - v.d * 0x1p-128L never normalizes by shifting up, + v.d * 0x1p-130L never normalizes by shifting up, so round bit plus sticky bit should be already enough for proper rounding. */ - if (v.ieee.exponent == 128) + if (v.ieee.exponent == 130) { /* If the exponent would be in the normal range when rounding to normal precision with unbounded exponent @@ -246,8 +271,8 @@ __fmal (long double x, long double y, long double z) if (TININESS_AFTER_ROUNDING) { w.d = a1 + u.d; - if (w.ieee.exponent == 129) - return w.d * 0x1p-128L; + if (w.ieee.exponent == 131) + return w.d * 0x1p-130L; } /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding, v.ieee.mantissa1 & 1 is the round bit and j is our sticky @@ -256,12 +281,12 @@ __fmal (long double x, long double y, long double z) w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j; w.ieee.negative = v.ieee.negative; v.ieee.mantissa1 &= ~3U; - v.d *= 0x1p-128L; + v.d *= 0x1p-130L; w.d *= 0x1p-2L; return v.d + w.d; } v.ieee.mantissa1 |= j; - return v.d * 0x1p-128L; + return v.d * 0x1p-130L; } } weak_alias (__fmal, fmal) diff --git a/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c deleted file mode 100644 index 3df59c2239..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c +++ /dev/null @@ -1,44 +0,0 @@ -/* Return classification value corresponding to argument. - Copyright (C) 1997, 2002 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. - Fixed by Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - The GNU C Library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - 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/>. */ - -#include <math.h> - -#include <math_private.h> - - -int -__fpclassifyl (long double x) -{ - u_int32_t ex, hx, lx, m; - int retval = FP_NORMAL; - - GET_LDOUBLE_WORDS (ex, hx, lx, x); - m = (hx & 0x7fffffff) | lx; - ex &= 0x7fff; - if ((ex | m) == 0) - retval = FP_ZERO; - else if (ex == 0 && (hx & 0x80000000) == 0) - retval = FP_SUBNORMAL; - else if (ex == 0x7fff) - retval = m != 0 ? FP_NAN : FP_INFINITE; - - return retval; -} -libm_hidden_def (__fpclassifyl) diff --git a/sysdeps/ieee754/ldbl-96/s_isinfl.c b/sysdeps/ieee754/ldbl-96/s_isinfl.c deleted file mode 100644 index 94639f00f8..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_isinfl.c +++ /dev/null @@ -1,30 +0,0 @@ -/* - * Written by J.T. Conklin <jtc@netbsd.org>. - * Change for long double by Ulrich Drepper <drepper@cygnus.com>. - * Public domain. - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -int -__isinfl (long double x) -{ - int32_t se,hx,lx; - GET_LDOUBLE_WORDS(se,hx,lx,x); - lx |= (hx & 0x7fffffff) | ((se & 0x7fff) ^ 0x7fff); - lx |= -lx; - se &= 0x8000; - return ~(lx >> 31) & (1 - (se >> 14)); -} -hidden_def (__isinfl) -weak_alias (__isinfl, isinfl) diff --git a/sysdeps/ieee754/ldbl-96/s_isnanl.c b/sysdeps/ieee754/ldbl-96/s_isnanl.c deleted file mode 100644 index fd270fd849..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_isnanl.c +++ /dev/null @@ -1,40 +0,0 @@ -/* s_isnanl.c -- long double version of s_isnan.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * isnanl(x) returns 1 is x is nan, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -int __isnanl(long double x) -{ - int32_t se,hx,lx; - GET_LDOUBLE_WORDS(se,hx,lx,x); - se = (se & 0x7fff) << 1; - lx |= hx & 0x7fffffff; - se |= (u_int32_t)(lx|(-lx))>>31; - se = 0xfffe - se; - return (int)(((u_int32_t)(se))>>31); -} -hidden_def (__isnanl) -weak_alias (__isnanl, isnanl) diff --git a/sysdeps/ieee754/ldbl-96/s_issignalingl.c b/sysdeps/ieee754/ldbl-96/s_issignalingl.c new file mode 100644 index 0000000000..bbb98735f5 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_issignalingl.c @@ -0,0 +1,43 @@ +/* Test for signaling NaN. + Copyright (C) 2013-2014 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + 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/>. */ + +#include <math.h> +#include <math_private.h> + +int +__issignalingl (long double x) +{ + u_int32_t exi, hxi, lxi; + GET_LDOUBLE_WORDS (exi, hxi, lxi, x); +#ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN +# error not implemented +#else + /* To keep the following comparison simple, toggle the quiet/signaling bit, + so that it is set for sNaNs. This is inverse to IEEE 754-2008 (as well as + common practice for IEEE 754-1985). */ + hxi ^= 0x40000000; + /* If lxi != 0, then set any suitable bit of the significand in hxi. */ + hxi |= (lxi | -lxi) >> 31; + /* We do not recognize a pseudo NaN as sNaN; they're invalid on 80387 and + later. */ + /* We have to compare for greater (instead of greater or equal), because x's + significand being all-zero designates infinity not NaN. */ + return ((exi & 0x7fff) == 0x7fff) && (hxi > 0xc0000000); +#endif +} +libm_hidden_def (__issignalingl) diff --git a/sysdeps/ieee754/ldbl-96/s_llrintl.c b/sysdeps/ieee754/ldbl-96/s_llrintl.c index 32a614b66d..0102f18242 100644 --- a/sysdeps/ieee754/ldbl-96/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_llrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_llroundl.c b/sysdeps/ieee754/ldbl-96/s_llroundl.c index e75cd43301..1559540a7a 100644 --- a/sysdeps/ieee754/ldbl-96/s_llroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_llroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long long int. - Copyright (C) 1997, 2001 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_logbl.c b/sysdeps/ieee754/ldbl-96/s_logbl.c deleted file mode 100644 index 4289be1933..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_logbl.c +++ /dev/null @@ -1,51 +0,0 @@ -/* s_logbl.c -- long double version of s_logb.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * long double logbl(x) - * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. - * Use ilogb instead. - */ - -#include <math.h> -#include <math_private.h> - -long double -__logbl (long double x) -{ - int32_t es, lx, ix; - - GET_LDOUBLE_WORDS (es, ix, lx, x); - es &= 0x7fff; /* exponent */ - if ((es | ix | lx) == 0) - return -1.0 / fabs (x); - if (es == 0x7fff) - return x * x; - if (es == 0) /* IEEE 754 logb */ - { - /* POSIX specifies that denormal number is treated as - though it were normalized. */ - int ma; - if (ix == 0) - ma = __builtin_clz (lx) + 32; - else - ma = __builtin_clz (ix); - es -= ma - 1; - } - return (long double) (es - 16383); -} - -weak_alias (__logbl, logbl) diff --git a/sysdeps/ieee754/ldbl-96/s_lrintl.c b/sysdeps/ieee754/ldbl-96/s_lrintl.c index 3f61099dd0..a668934468 100644 --- a/sysdeps/ieee754/ldbl-96/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_lrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_lroundl.c b/sysdeps/ieee754/ldbl-96/s_lroundl.c index d41b7ca95b..a5c08d6391 100644 --- a/sysdeps/ieee754/ldbl-96/s_lroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_lroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long int. - Copyright (C) 1997, 2004 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_nearbyintl.c b/sysdeps/ieee754/ldbl-96/s_nearbyintl.c deleted file mode 100644 index ed9836c879..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_nearbyintl.c +++ /dev/null @@ -1,86 +0,0 @@ -/* s_rintl.c -- long double version of s_rint.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ -/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * rintl(x) - * Return x rounded to integral value according to the prevailing - * rounding mode. - * Method: - * Using floating addition. - * Exception: - * Inexact flag raised if x not equal to rintl(x). - */ - -#include <fenv.h> -#include <math.h> -#include <math_private.h> - -static const long double -TWO63[2]={ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - -long double __nearbyintl(long double x) -{ - fenv_t env; - int32_t se,j0,sx; - u_int32_t i,i0,i1; - long double w,t; - GET_LDOUBLE_WORDS(se,i0,i1,x); - sx = (se>>15)&1; - j0 = (se&0x7fff)-0x3fff; - if(j0<31) { - if(j0<0) { - if(((se&0x7fff)|i0|i1)==0) return x; - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1|-i1)&0x80000000; - SET_LDOUBLE_MSW(x,i0); - feholdexcept (&env); - w = TWO63[sx]+x; - t = w-TWO63[sx]; - fesetenv (&env); - GET_LDOUBLE_EXP(i0,t); - SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15)); - return t; - } else { - i = (0x7fffffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - i>>=1; - if(((i0&i)|i1)!=0) { - if (j0==30) i1 = 0x40000000; else - i0 = (i0&(~i))|((0x20000000)>>j0); - } - } - } else if (j0>62) { - if(j0==0x4000) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-31); - if((i1&i)==0) return x; /* x is integral */ - i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31)); - } - SET_LDOUBLE_WORDS(x,se,i0,i1); - feholdexcept (&env); - w = TWO63[sx]+x; - t = w-TWO63[sx]; - fesetenv (&env); - return t; -} -weak_alias (__nearbyintl, nearbyintl) diff --git a/sysdeps/ieee754/ldbl-96/s_nextafterl.c b/sysdeps/ieee754/ldbl-96/s_nextafterl.c deleted file mode 100644 index 6859349b7c..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_nextafterl.c +++ /dev/null @@ -1,96 +0,0 @@ -/* s_nextafterl.c -- long double version of s_nextafter.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* IEEE functions - * nextafterl(x,y) - * return the next machine floating-point number of x in the - * direction toward y. - * Special cases: - */ - -#include <math.h> -#include <math_private.h> - -long double __nextafterl(long double x, long double y) -{ - int32_t hx,hy,ix,iy; - u_int32_t lx,ly,esx,esy; - - GET_LDOUBLE_WORDS(esx,hx,lx,x); - GET_LDOUBLE_WORDS(esy,hy,ly,y); - ix = esx&0x7fff; /* |x| */ - iy = esy&0x7fff; /* |y| */ - - if (((ix==0x7fff)&&((hx|lx)!=0)) || /* x is nan */ - ((iy==0x7fff)&&((hy|ly)!=0))) /* y is nan */ - return x+y; - if(x==y) return y; /* x=y, return y */ - if((ix|hx|lx)==0) { /* x == 0 */ - long double u; - SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */ - u = math_opt_barrier (x); - u = u * u; - math_force_eval (u); /* raise underflow flag */ - return x; - } - if(esx<0x8000) { /* x > 0 */ - if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) { - /* x > y, x -= ulp */ - if(lx==0) { - if (hx==0) esx -= 1; - hx -= 1; - } - lx -= 1; - } else { /* x < y, x += ulp */ - lx += 1; - if(lx==0) { - hx += 1; - if (hx==0) - esx += 1; - } - } - } else { /* x < 0 */ - if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){ - /* x < y, x -= ulp */ - if(lx==0) { - if (hx==0) esx -= 1; - hx -= 1; - } - lx -= 1; - } else { /* x > y, x += ulp */ - lx += 1; - if(lx==0) { - hx += 1; - if (hx==0) esx += 1; - } - } - } - esy = esx&0x7fff; - if(esy==0x7fff) return x+x; /* overflow */ - if(esy==0) { - long double u = x*x; /* underflow */ - math_force_eval (u); /* raise underflow flag */ - } - SET_LDOUBLE_WORDS(x,esx,hx,lx); - return x; -} -weak_alias (__nextafterl, nextafterl) -strong_alias (__nextafterl, __nexttowardl) -weak_alias (__nextafterl, nexttowardl) diff --git a/sysdeps/ieee754/ldbl-96/s_remquol.c b/sysdeps/ieee754/ldbl-96/s_remquol.c index a0825c9521..866935c6e7 100644 --- a/sysdeps/ieee754/ldbl-96/s_remquol.c +++ b/sysdeps/ieee754/ldbl-96/s_remquol.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_rintl.c b/sysdeps/ieee754/ldbl-96/s_rintl.c deleted file mode 100644 index b6f899d4ef..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_rintl.c +++ /dev/null @@ -1,82 +0,0 @@ -/* s_rintl.c -- long double version of s_rint.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * rintl(x) - * Return x rounded to integral value according to the prevailing - * rounding mode. - * Method: - * Using floating addition. - * Exception: - * Inexact flag raised if x not equal to rintl(x). - */ - -#include <math.h> -#include <math_private.h> - -static const long double -TWO63[2]={ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - -long double __rintl(long double x) -{ - int32_t se,j0,sx; - u_int32_t i,i0,i1; - long double w,t; - GET_LDOUBLE_WORDS(se,i0,i1,x); - sx = (se>>15)&1; - j0 = (se&0x7fff)-0x3fff; - if(j0<31) { - if(j0<0) { - if(((se&0x7fff)|i0|i1)==0) return x; - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1|-i1)&0x80000000; - SET_LDOUBLE_MSW(x,i0); - w = TWO63[sx]+x; - t = w-TWO63[sx]; - GET_LDOUBLE_EXP(i0,t); - SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15)); - return t; - } else { - i = (0x7fffffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - i>>=1; - if(((i0&i)|i1)!=0) { - if(j0==30) i1 = 0x40000000; else - i0 = (i0&(~i))|((0x20000000)>>j0); - } - } - } else if (j0>62) { - if(j0==0x4000) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-31); - if((i1&i)==0) return x; /* x is integral */ - i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31)); - } - SET_LDOUBLE_WORDS(x,se,i0,i1); - w = TWO63[sx]+x; - return w-TWO63[sx]; -} -weak_alias (__rintl, rintl) diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c index b8626a509e..e16e9c6f65 100644 --- a/sysdeps/ieee754/ldbl-96/s_roundl.c +++ b/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -1,5 +1,5 @@ /* Round long double to integer away from zero. - Copyright (C) 1997, 2007, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_scalbnl.c b/sysdeps/ieee754/ldbl-96/s_scalbnl.c deleted file mode 100644 index 266a37b9c0..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_scalbnl.c +++ /dev/null @@ -1,61 +0,0 @@ -/* s_scalbnl.c -- long double version of s_scalbn.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * scalbnl (long double x, int n) - * scalbnl(x,n) returns x* 2**n computed by exponent - * manipulation rather than by actually performing an - * exponentiation or a multiplication. - */ - -#include <math.h> -#include <math_private.h> - -static const long double -two64 = 1.8446744073709551616e19L, -twom64 = 5.421010862427522170037e-20L, -huge = 1.0e+4900L, -tiny = 1.0e-4900L; - -long double -__scalbnl (long double x, int n) -{ - int32_t k,es,hx,lx; - GET_LDOUBLE_WORDS(es,hx,lx,x); - k = es&0x7fff; /* extract exponent */ - if (__builtin_expect(k==0, 0)) { /* 0 or subnormal x */ - if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ - x *= two64; - GET_LDOUBLE_EXP(hx,x); - k = (hx&0x7fff) - 64; - } - if (__builtin_expect(k==0x7fff, 0)) return x+x; /* NaN or Inf */ - if (__builtin_expect(n< -50000, 0)) - return tiny*__copysignl(tiny,x); - if (__builtin_expect(n> 50000 || k+n > 0x7ffe, 0)) - return huge*__copysignl(huge,x); /* overflow */ - /* Now k and n are bounded we know that k = k+n does not - overflow. */ - k = k+n; - if (__builtin_expect(k > 0, 1)) /* normal result */ - {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;} - if (k <= -64) - return tiny*__copysignl(tiny,x); /*underflow*/ - k += 64; /* subnormal result */ - SET_LDOUBLE_EXP(x,(es&0x8000)|k); - return x*twom64; -} -weak_alias (__scalbnl, scalbnl) diff --git a/sysdeps/ieee754/ldbl-96/s_signbitl.c b/sysdeps/ieee754/ldbl-96/s_signbitl.c index ffded99494..8256d9ddb0 100644 --- a/sysdeps/ieee754/ldbl-96/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-96/s_signbitl.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_sincosl.c b/sysdeps/ieee754/ldbl-96/s_sincosl.c index 2a819592c1..af8f899c24 100644 --- a/sysdeps/ieee754/ldbl-96/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/s_sincosl.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2012 Free Software Foundation, Inc. + Copyright (C) 1997-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_truncl.c b/sysdeps/ieee754/ldbl-96/s_truncl.c deleted file mode 100644 index 75a38261af..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_truncl.c +++ /dev/null @@ -1,56 +0,0 @@ -/* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997, 1999 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - The GNU C Library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - 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/>. */ - -#include <math.h> - -#include <math_private.h> - - -long double -__truncl (long double x) -{ - int32_t i0, j0; - u_int32_t se, i1; - int sx; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - sx = se & 0x8000; - j0 = (se & 0x7fff) - 0x3fff; - if (j0 < 31) - { - if (j0 < 0) - /* The magnitude of the number is < 1 so the result is +-0. */ - SET_LDOUBLE_WORDS (x, sx, 0, 0); - else - SET_LDOUBLE_WORDS (x, se, i0 & ~(0x7fffffff >> j0), 0); - } - else if (j0 > 63) - { - if (j0 == 0x4000) - /* x is inf or NaN. */ - return x + x; - } - else - { - SET_LDOUBLE_WORDS (x, se, i0, i1 & ~(0xffffffffu >> (j0 - 31))); - } - - return x; -} -weak_alias (__truncl, truncl) diff --git a/sysdeps/ieee754/ldbl-96/strtold_l.c b/sysdeps/ieee754/ldbl-96/strtold_l.c index bc5a52ea48..82c33e52d6 100644 --- a/sysdeps/ieee754/ldbl-96/strtold_l.c +++ b/sysdeps/ieee754/ldbl-96/strtold_l.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1999, 2004, 2005 Free Software Foundation, Inc. +/* Copyright (C) 1999-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -34,11 +34,10 @@ #define SET_MANTISSA(flt, mant) \ do { union ieee854_long_double u; \ u.d = (flt); \ - if ((mant & 0x7fffffffffffffffULL) == 0) \ - mant = 0x4000000000000000ULL; \ - u.ieee.mantissa0 = (((mant) >> 32) & 0x7fffffff) | 0x80000000; \ - u.ieee.mantissa1 = (mant) & 0xffffffff; \ - (flt) = u.d; \ + u.ieee_nan.mantissa0 = (mant) >> 32; \ + u.ieee_nan.mantissa1 = (mant); \ + if ((u.ieee.mantissa0 | u.ieee.mantissa1) != 0) \ + (flt) = u.d; \ } while (0) #include <stdlib/strtod_l.c> diff --git a/sysdeps/ieee754/ldbl-96/t_sincosl.c b/sysdeps/ieee754/ldbl-96/t_sincosl.c index 50dd84b948..da8fcfa112 100644 --- a/sysdeps/ieee754/ldbl-96/t_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/t_sincosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point sine and cosine tables. - Copyright (C) 1999-2012 Free Software Foundation, Inc. + Copyright (C) 1999-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision tables by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c index 79b10c5756..af221184d1 100644 --- a/sysdeps/ieee754/ldbl-96/w_expl.c +++ b/sysdeps/ieee754/ldbl-96/w_expl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011, 2012 Free Software Foundation, Inc. +/* Copyright (C) 2011-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1.c b/sysdeps/ieee754/ldbl-96/x2y2m1.c index 3086fa2c80..6ae2a63cbe 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1.c +++ b/sysdeps/ieee754/ldbl-96/x2y2m1.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012 Free Software Foundation, Inc. + Copyright (C) 2012-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1l.c b/sysdeps/ieee754/ldbl-96/x2y2m1l.c index a249cd3479..575bcd812d 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-96/x2y2m1l.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012 Free Software Foundation, Inc. + Copyright (C) 2012-2014 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or |