diff options
Diffstat (limited to 'sysdeps/ieee754/flt-32')
47 files changed, 457 insertions, 208 deletions
diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c index 00bad4239b..2ca2dbcb28 100644 --- a/sysdeps/ieee754/flt-32/e_asinf.c +++ b/sysdeps/ieee754/flt-32/e_asinf.c @@ -73,11 +73,7 @@ float __ieee754_asinf(float x) return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix<0x3f000000) { /* |x|<0.5 */ if(ix<0x32000000) { /* if |x| < 2**-27 */ - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if(huge+x>one) return x;/* return x with inexact if x!=0*/ } else { t = x*x; diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c index bc74960e16..82de071ab8 100644 --- a/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/sysdeps/ieee754/flt-32/e_atanhf.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2015 Free Software Foundation, Inc. +/* Copyright (C) 2011-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -52,11 +52,7 @@ __ieee754_atanhf (float x) if (__glibc_unlikely (xa < 0x1.0p-28f)) { math_force_eval (huge + x); - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c index dedda47c09..7b223758e1 100644 --- a/sysdeps/ieee754/flt-32/e_coshf.c +++ b/sysdeps/ieee754/flt-32/e_coshf.c @@ -58,6 +58,6 @@ __ieee754_coshf (float x) if(ix>=0x7f800000) return x*x; /* |x| > overflowthresold, cosh(x) overflow */ - return huge*huge; + return math_narrow_eval (huge*huge); } strong_alias (__ieee754_coshf, __coshf_finite) diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 0b75a7ea2a..1723c482de 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -1,5 +1,5 @@ /* Single-precision floating point 2^x. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> @@ -37,8 +37,8 @@ #include "t_exp2f.h" -static const volatile float TWOM100 = 7.88860905e-31; -static const volatile float TWO127 = 1.7014118346e+38; +static const float TWOM100 = 7.88860905e-31; +static const float TWO127 = 1.7014118346e+38; float __ieee754_exp2f (float x) @@ -109,12 +109,16 @@ __ieee754_exp2f (float x) if (!unsafe) return result; else - return result * scale_u.f; + { + result *= scale_u.f; + math_check_force_underflow_nonneg (result); + return result; + } } /* Exceptional cases: */ else if (isless (x, himark)) { - if (__isinf_nsf (x)) + if (isinf (x)) /* e^-inf == 0, with no error. */ return 0; else diff --git a/sysdeps/ieee754/flt-32/e_expf.c b/sysdeps/ieee754/flt-32/e_expf.c index abf9111a74..071f615ef4 100644 --- a/sysdeps/ieee754/flt-32/e_expf.c +++ b/sysdeps/ieee754/flt-32/e_expf.c @@ -1,5 +1,5 @@ /* Single-precision floating point e^x. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c index 29fe8b46c2..19f51b0c8b 100644 --- a/sysdeps/ieee754/flt-32/e_gammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -69,11 +69,7 @@ gammaf_positive (float x, int *exp2_adj) /* Adjust into the range for applying Stirling's approximation. */ float n = __ceilf (4.0f - x); -#if FLT_EVAL_METHOD != 0 - volatile -#endif - float x_tmp = x + n; - x_adj = x_tmp; + x_adj = math_narrow_eval (x + n); x_eps = (x - (x_adj - n)); prod = __gamma_productf (x_adj - n, x_eps, n, &eps); } @@ -111,9 +107,6 @@ float __ieee754_gammaf_r (float x, int *signgamp) { int32_t hx; -#if FLT_EVAL_METHOD != 0 - volatile -#endif float ret; GET_FLOAT_WORD (hx, x); @@ -149,7 +142,7 @@ __ieee754_gammaf_r (float x, int *signgamp) { /* Overflow. */ *signgamp = 0; - ret = FLT_MAX * FLT_MAX; + ret = math_narrow_eval (FLT_MAX * FLT_MAX); return ret; } else @@ -186,29 +179,31 @@ __ieee754_gammaf_r (float x, int *signgamp) float tret = (float) M_PI / (-x * sinpix * gammaf_positive (-x, &exp2_adj)); ret = __scalbnf (tret, -exp2_adj); + math_check_force_underflow_nonneg (ret); } } + ret = math_narrow_eval (ret); } if (isinf (ret) && x != 0) { if (*signgamp < 0) { - ret = -__copysignf (FLT_MAX, ret) * FLT_MAX; + ret = math_narrow_eval (-__copysignf (FLT_MAX, ret) * FLT_MAX); ret = -ret; } else - ret = __copysignf (FLT_MAX, ret) * FLT_MAX; + ret = math_narrow_eval (__copysignf (FLT_MAX, ret) * FLT_MAX); return ret; } else if (ret == 0) { if (*signgamp < 0) { - ret = -__copysignf (FLT_MIN, ret) * FLT_MIN; + ret = math_narrow_eval (-__copysignf (FLT_MIN, ret) * FLT_MIN); ret = -ret; } else - ret = __copysignf (FLT_MIN, ret) * FLT_MIN; + ret = math_narrow_eval (__copysignf (FLT_MIN, ret) * FLT_MIN); return ret; } else diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c index 8f499b5fde..717b82e42f 100644 --- a/sysdeps/ieee754/flt-32/e_hypotf.c +++ b/sysdeps/ieee754/flt-32/e_hypotf.c @@ -27,17 +27,9 @@ __ieee754_hypotf(float x, float y) GET_FLOAT_WORD(hb,y); hb &= 0x7fffffff; if (ha == 0x7f800000) - { - if (x == y) - return fabsf(y); - return fabsf(x); - } + return fabsf(x); else if (hb == 0x7f800000) - { - if (x == y) - return fabsf(x); - return fabsf(y); - } + return fabsf(y); else if (ha > 0x7f800000 || hb > 0x7f800000) return fabsf(x) * fabsf(y); else if (ha == 0) diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c index 63de21f609..f359a3d9ba 100644 --- a/sysdeps/ieee754/flt-32/e_j1f.c +++ b/sysdeps/ieee754/flt-32/e_j1f.c @@ -71,11 +71,10 @@ __ieee754_j1f(float x) } if(__builtin_expect(ix<0x32000000, 0)) { /* |x|<2**-27 */ if(huge+x>one) { /* inexact if x!=0 necessary */ - float ret = (float) 0.5 * x; - if (fabsf (ret) < FLT_MIN) { - float force_underflow = ret * ret; - math_force_eval (force_underflow); - } + float ret = math_narrow_eval ((float) 0.5 * x); + math_check_force_underflow (ret); + if (ret == 0 && x != 0) + __set_errno (ERANGE); return ret; } } diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c index 44a3761adb..4e634778d3 100644 --- a/sysdeps/ieee754/flt-32/e_jnf.c +++ b/sysdeps/ieee754/flt-32/e_jnf.c @@ -167,13 +167,15 @@ __ieee754_jnf(int n, float x) } } if(sgn==1) ret = -b; else ret = b; + ret = math_narrow_eval (ret); } if (ret == 0) - ret = __copysignf (FLT_MIN, ret) * FLT_MIN; - else if (fabsf (ret) < FLT_MIN) { - float force_underflow = ret * ret; - math_force_eval (force_underflow); - } + { + ret = math_narrow_eval (__copysignf (FLT_MIN, ret) * FLT_MIN); + __set_errno (ERANGE); + } + else + math_check_force_underflow (ret); return ret; } strong_alias (__ieee754_jnf, __jnf_finite) diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c index c0bf4156ff..4d8a66bb39 100644 --- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c @@ -161,6 +161,9 @@ __ieee754_lgammaf_r(float x, int *signgamp) if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ return x/zero; + if (ix > 0x40000000 /* X < 2.0f. */ + && ix < 0x41700000 /* X > -15.0f. */) + return __lgamma_negf (x, signgamp); t = sin_pif(x); if(t==zero) return one/fabsf(t); /* -integer */ nadj = __ieee754_logf(pi/fabsf(t*x)); @@ -229,17 +232,13 @@ __ieee754_lgammaf_r(float x, int *signgamp) r = (x-half)*(t-one)+w; } else /* 2**26 <= x <= inf */ - r = x*(__ieee754_logf(x)-one); + r = math_narrow_eval (x*(__ieee754_logf(x)-one)); /* NADJ is set for negative arguments but not otherwise, resulting in warnings that it may be used uninitialized although in the cases where it is used it has always been set. */ DIAG_PUSH_NEEDS_COMMENT; -#if __GNUC_PREREQ (4, 7) DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized"); -#else - DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized"); -#endif if(hx<0) r = nadj - r; DIAG_POP_NEEDS_COMMENT; return r; diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c index 96f0e81c52..2cd01b4a50 100644 --- a/sysdeps/ieee754/flt-32/e_log10f.c +++ b/sysdeps/ieee754/flt-32/e_log10f.c @@ -15,6 +15,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const float two25 = 3.3554432000e+07, /* 0x4c000000 */ @@ -22,8 +23,6 @@ ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ log10_2lo = 7.9034151668e-07; /* 0x355427db */ -static const float zero = 0.0; - float __ieee754_log10f(float x) { @@ -46,6 +45,8 @@ __ieee754_log10f(float x) i = ((u_int32_t)k&0x80000000)>>31; hx = (hx&0x007fffff)|((0x7f-i)<<23); y = (float)(k+i); + if (FIX_INT_FP_CONVERT_ZERO && y == 0.0f) + y = 0.0f; SET_FLOAT_WORD(x,hx); z = y*log10_2lo + ivln10*__ieee754_logf(x); return z+y*log10_2hi; diff --git a/sysdeps/ieee754/flt-32/e_log2f.c b/sysdeps/ieee754/flt-32/e_log2f.c index 245be4e6f7..857d13fb9b 100644 --- a/sysdeps/ieee754/flt-32/e_log2f.c +++ b/sysdeps/ieee754/flt-32/e_log2f.c @@ -17,6 +17,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const float ln2 = 0.69314718055994530942, @@ -57,7 +58,12 @@ __ieee754_log2f(float x) dk = (float)k; f = x-(float)1.0; if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */ - if(f==zero) return dk; + if(f==zero) + { + if (FIX_INT_FP_CONVERT_ZERO && dk == 0.0f) + dk = 0.0f; + return dk; + } R = f*f*((float)0.5-(float)0.33333333333333333*f); return dk-(R-f)/ln2; } diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c index 12c408f93c..c72fe37d3b 100644 --- a/sysdeps/ieee754/flt-32/e_powf.c +++ b/sysdeps/ieee754/flt-32/e_powf.c @@ -131,7 +131,7 @@ __ieee754_powf(float x, float y) if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ - t = x-1; /* t has 20 trailing zeros */ + t = ax-1; /* t has 20 trailing zeros */ w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25)); u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ v = t*ivln2_l-w*ivln2; @@ -166,7 +166,9 @@ __ieee754_powf(float x, float y) GET_FLOAT_WORD(is,s_h); SET_FLOAT_WORD(s_h,is&0xfffff000); /* t_h=ax+bp[k] High */ - SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21)); + SET_FLOAT_WORD (t_h, + ((((ix>>1)|0x20000000)+0x00400000+(k<<21)) + & 0xfffff000)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ @@ -244,7 +246,12 @@ __ieee754_powf(float x, float y) z = one-(r-z); GET_FLOAT_WORD(j,z); j += (n<<23); - if((j>>23)<=0) z = __scalbnf(z,n); /* subnormal output */ + if((j>>23)<=0) /* subnormal output */ + { + z = __scalbnf (z, n); + float force_underflow = z * z; + math_force_eval (force_underflow); + } else SET_FLOAT_WORD(z,j); return s*z; } diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c index 0a3cc90cd1..6100d95c55 100644 --- a/sysdeps/ieee754/flt-32/e_sinhf.c +++ b/sysdeps/ieee754/flt-32/e_sinhf.c @@ -13,6 +13,7 @@ * ==================================================== */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -34,8 +35,10 @@ __ieee754_sinhf(float x) if (jx<0) h = -h; /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x41b00000) { /* |x|<22 */ - if (__builtin_expect(ix<0x31800000, 0)) /* |x|<2**-28 */ + if (__builtin_expect(ix<0x31800000, 0)) { /* |x|<2**-28 */ + math_check_force_underflow (x); if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + } t = __expm1f(fabsf(x)); if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one)); return h*(t+t/(t+one)); @@ -52,6 +55,6 @@ __ieee754_sinhf(float x) } /* |x| > overflowthresold, sinh(x) overflow */ - return x*shuge; + return math_narrow_eval (x*shuge); } strong_alias (__ieee754_sinhf, __sinhf_finite) diff --git a/sysdeps/ieee754/flt-32/k_rem_pio2f.c b/sysdeps/ieee754/flt-32/k_rem_pio2f.c index 6f14d5bac7..392afdbb6c 100644 --- a/sysdeps/ieee754/flt-32/k_rem_pio2f.c +++ b/sysdeps/ieee754/flt-32/k_rem_pio2f.c @@ -65,7 +65,9 @@ int __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const int32 /* compute q[0],q[1],...q[jk] */ for (i=0;i<=jk;i++) { - for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; + for(j=0,fw=0.0;j<=jx;j++) + fw += x[j]*f[jx+i-j]; + q[i] = fw; } jz = jk; @@ -167,30 +169,21 @@ recompute: break; case 1: case 2:; -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif float fv = 0.0; - for (i=jz;i>=0;i--) fv += fq[i]; + for (i=jz;i>=0;i--) fv = math_narrow_eval (fv + fq[i]); y[0] = (ih==0)? fv: -fv; - fv = fq[0]-fv; - for (i=1;i<=jz;i++) fv += fq[i]; + fv = math_narrow_eval (fq[0]-fv); + for (i=1;i<=jz;i++) fv = math_narrow_eval (fv + fq[i]); y[1] = (ih==0)? fv: -fv; break; case 3: /* painful */ for (i=jz;i>0;i--) { -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif - float fv = fq[i-1]+fq[i]; + float fv = math_narrow_eval (fq[i-1]+fq[i]); fq[i] += fq[i-1]-fv; fq[i-1] = fv; } for (i=jz;i>1;i--) { -#if __FLT_EVAL_METHOD__ != 0 - volatile -#endif - float fv = fq[i-1]+fq[i]; + float fv = math_narrow_eval (fq[i-1]+fq[i]); fq[i] += fq[i-1]-fv; fq[i-1] = fv; } diff --git a/sysdeps/ieee754/flt-32/k_sinf.c b/sysdeps/ieee754/flt-32/k_sinf.c index 0c98a2ae91..a195d59466 100644 --- a/sysdeps/ieee754/flt-32/k_sinf.c +++ b/sysdeps/ieee754/flt-32/k_sinf.c @@ -38,11 +38,7 @@ float __kernel_sinf(float x, float y, int iy) ix &= 0x7fffffff; /* high word of x */ if(ix<0x32000000) /* |x| < 2**-27 */ { - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if ((int) x == 0) return x; /* generate inexact */ } diff --git a/sysdeps/ieee754/flt-32/k_tanf.c b/sysdeps/ieee754/flt-32/k_tanf.c index a67f36e283..9f0e55860f 100644 --- a/sysdeps/ieee754/flt-32/k_tanf.c +++ b/sysdeps/ieee754/flt-32/k_tanf.c @@ -17,6 +17,7 @@ static char rcsid[] = "$NetBSD: k_tanf.c,v 1.4 1995/05/10 20:46:39 jtc Exp $"; #endif +#include <float.h> #include <math.h> #include <math_private.h> static const float @@ -48,7 +49,13 @@ float __kernel_tanf(float x, float y, int iy) if(ix<0x39000000) /* x < 2**-13 */ {if((int)x==0) { /* generate inexact */ if((ix|(iy+1))==0) return one/fabsf(x); - else return (iy==1)? x: -one/x; + else if (iy == 1) + { + math_check_force_underflow (x); + return x; + } + else + return -one / x; } } if(ix>=0x3f2ca140) { /* |x|>=0.6744 */ diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c new file mode 100644 index 0000000000..f15719b059 --- /dev/null +++ b/sysdeps/ieee754/flt-32/lgamma_negf.c @@ -0,0 +1,280 @@ +/* lgammaf expanding around zeros. + Copyright (C) 2015-2016 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 <float.h> +#include <math.h> +#include <math_private.h> + +static const float lgamma_zeros[][2] = + { + { -0x2.74ff94p+0f, 0x1.3fe0f2p-24f }, + { -0x2.bf682p+0f, -0x1.437b2p-24f }, + { -0x3.24c1b8p+0f, 0x6.c34cap-28f }, + { -0x3.f48e2cp+0f, 0x1.707a04p-24f }, + { -0x4.0a13ap+0f, 0x1.e99aap-24f }, + { -0x4.fdd5ep+0f, 0x1.64454p-24f }, + { -0x5.021a98p+0f, 0x2.03d248p-24f }, + { -0x5.ffa4cp+0f, 0x2.9b82fcp-24f }, + { -0x6.005ac8p+0f, -0x1.625f24p-24f }, + { -0x6.fff3p+0f, 0x2.251e44p-24f }, + { -0x7.000dp+0f, 0x8.48078p-28f }, + { -0x7.fffe6p+0f, 0x1.fa98c4p-28f }, + { -0x8.0001ap+0f, -0x1.459fcap-28f }, + { -0x8.ffffdp+0f, -0x1.c425e8p-24f }, + { -0x9.00003p+0f, 0x1.c44b82p-24f }, + { -0xap+0f, 0x4.9f942p-24f }, + { -0xap+0f, -0x4.9f93b8p-24f }, + { -0xbp+0f, 0x6.b9916p-28f }, + { -0xbp+0f, -0x6.b9915p-28f }, + { -0xcp+0f, 0x8.f76c8p-32f }, + { -0xcp+0f, -0x8.f76c7p-32f }, + { -0xdp+0f, 0xb.09231p-36f }, + { -0xdp+0f, -0xb.09231p-36f }, + { -0xep+0f, 0xc.9cba5p-40f }, + { -0xep+0f, -0xc.9cba5p-40f }, + { -0xfp+0f, 0xd.73f9fp-44f }, + }; + +static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f; + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const float lgamma_coeff[] = + { + 0x1.555556p-4f, + -0xb.60b61p-12f, + 0x3.403404p-12f, + }; + +#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0])) + +/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is + the integer end-point of the half-integer interval containing x and + x0 is the zero of lgamma in that half-integer interval. Each + polynomial is expressed in terms of x-xm, where xm is the midpoint + of the interval for which the polynomial applies. */ + +static const float poly_coeff[] = + { + /* Interval [-2.125, -2] (polynomial degree 5). */ + -0x1.0b71c6p+0f, + -0xc.73a1ep-4f, + -0x1.ec8462p-4f, + -0xe.37b93p-4f, + -0x1.02ed36p-4f, + -0xe.cbe26p-4f, + /* Interval [-2.25, -2.125] (polynomial degree 5). */ + -0xf.29309p-4f, + -0xc.a5cfep-4f, + 0x3.9c93fcp-4f, + -0x1.02a2fp+0f, + 0x9.896bep-4f, + -0x1.519704p+0f, + /* Interval [-2.375, -2.25] (polynomial degree 5). */ + -0xd.7d28dp-4f, + -0xe.6964cp-4f, + 0xb.0d4f1p-4f, + -0x1.9240aep+0f, + 0x1.dadabap+0f, + -0x3.1778c4p+0f, + /* Interval [-2.5, -2.375] (polynomial degree 6). */ + -0xb.74ea2p-4f, + -0x1.2a82cp+0f, + 0x1.880234p+0f, + -0x3.320c4p+0f, + 0x5.572a38p+0f, + -0x9.f92bap+0f, + 0x1.1c347ep+4f, + /* Interval [-2.625, -2.5] (polynomial degree 6). */ + -0x3.d10108p-4f, + 0x1.cd5584p+0f, + 0x3.819c24p+0f, + 0x6.84cbb8p+0f, + 0xb.bf269p+0f, + 0x1.57fb12p+4f, + 0x2.7b9854p+4f, + /* Interval [-2.75, -2.625] (polynomial degree 6). */ + -0x6.b5d25p-4f, + 0x1.28d604p+0f, + 0x1.db6526p+0f, + 0x2.e20b38p+0f, + 0x4.44c378p+0f, + 0x6.62a08p+0f, + 0x9.6db3ap+0f, + /* Interval [-2.875, -2.75] (polynomial degree 5). */ + -0x8.a41b2p-4f, + 0xc.da87fp-4f, + 0x1.147312p+0f, + 0x1.7617dap+0f, + 0x1.d6c13p+0f, + 0x2.57a358p+0f, + /* Interval [-3, -2.875] (polynomial degree 5). */ + -0xa.046d6p-4f, + 0x9.70b89p-4f, + 0xa.a89a6p-4f, + 0xd.2f2d8p-4f, + 0xd.e32b4p-4f, + 0xf.fb741p-4f, + }; + +static const size_t poly_deg[] = + { + 5, + 5, + 5, + 6, + 6, + 6, + 5, + 5, + }; + +static const size_t poly_end[] = + { + 5, + 11, + 17, + 24, + 31, + 38, + 44, + 50, + }; + +/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_sinpi (float x) +{ + if (x <= 0.25f) + return __sinf ((float) M_PI * x); + else + return __cosf ((float) M_PI * (0.5f - x)); +} + +/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_cospi (float x) +{ + if (x <= 0.25f) + return __cosf ((float) M_PI * x); + else + return __sinf ((float) M_PI * (0.5f - x)); +} + +/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ + +static float +lg_cotpi (float x) +{ + return lg_cospi (x) / lg_sinpi (x); +} + +/* Compute lgamma of a negative argument -15 < X < -2, setting + *SIGNGAMP accordingly. */ + +float +__lgamma_negf (float x, int *signgamp) +{ + /* Determine the half-integer region X lies in, handle exact + integers and determine the sign of the result. */ + int i = __floorf (-2 * x); + if ((i & 1) == 0 && i == -2 * x) + return 1.0f / 0.0f; + float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); + i -= 4; + *signgamp = ((i & 2) == 0 ? -1 : 1); + + SET_RESTORE_ROUNDF (FE_TONEAREST); + + /* Expand around the zero X0 = X0_HI + X0_LO. */ + float x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; + float xdiff = x - x0_hi - x0_lo; + + /* For arguments in the range -3 to -2, use polynomial + approximations to an adjusted version of the gamma function. */ + if (i < 2) + { + int j = __floorf (-8 * x) - 16; + float xm = (-33 - 2 * j) * 0.0625f; + float x_adj = x - xm; + size_t deg = poly_deg[j]; + size_t end = poly_end[j]; + float g = poly_coeff[end]; + for (size_t j = 1; j <= deg; j++) + g = g * x_adj + poly_coeff[end - j]; + return __log1pf (g * xdiff / (x - xn)); + } + + /* The result we want is log (sinpi (X0) / sinpi (X)) + + log (gamma (1 - X0) / gamma (1 - X)). */ + float x_idiff = fabsf (xn - x), x0_idiff = fabsf (xn - x0_hi - x0_lo); + float log_sinpi_ratio; + if (x0_idiff < x_idiff * 0.5f) + /* Use log not log1p to avoid inaccuracy from log1p of arguments + close to -1. */ + log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff) + / lg_sinpi (x_idiff)); + else + { + /* Use log1p not log to avoid inaccuracy from log of arguments + close to 1. X0DIFF2 has positive sign if X0 is further from + XN than X is from XN, negative sign otherwise. */ + float x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5f; + float sx0d2 = lg_sinpi (x0diff2); + float cx0d2 = lg_cospi (x0diff2); + log_sinpi_ratio = __log1pf (2 * sx0d2 + * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); + } + + float log_gamma_ratio; + float y0 = math_narrow_eval (1 - x0_hi); + float y0_eps = -x0_hi + (1 - y0) - x0_lo; + float y = math_narrow_eval (1 - x); + float y_eps = -x + (1 - y); + /* We now wish to compute LOG_GAMMA_RATIO + = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF + accurately approximates the difference Y0 + Y0_EPS - Y - + Y_EPS. Use Stirling's approximation. */ + float log_gamma_high + = (xdiff * __log1pf ((y0 - e_hi - e_lo + y0_eps) / e_hi) + + (y - 0.5f + y_eps) * __log1pf (xdiff / y)); + /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ + float y0r = 1 / y0, yr = 1 / y; + float y0r2 = y0r * y0r, yr2 = yr * yr; + float rdiff = -xdiff / (y * y0); + float bterm[NCOEFF]; + float dlast = rdiff, elast = rdiff * yr * (yr + y0r); + bterm[0] = dlast * lgamma_coeff[0]; + for (size_t j = 1; j < NCOEFF; j++) + { + float dnext = dlast * y0r2 + elast; + float enext = elast * yr2; + bterm[j] = dnext * lgamma_coeff[j]; + dlast = dnext; + elast = enext; + } + float log_gamma_low = 0; + for (size_t j = 0; j < NCOEFF; j++) + log_gamma_low += bterm[NCOEFF - 1 - j]; + log_gamma_ratio = log_gamma_high + log_gamma_low; + + return log_sinpi_ratio + log_gamma_ratio; +} diff --git a/sysdeps/ieee754/flt-32/lgamma_productf.c b/sysdeps/ieee754/flt-32/lgamma_productf.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/lgamma_productf.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/math_private.h b/sysdeps/ieee754/flt-32/math_private.h deleted file mode 100644 index e33db02b4c..0000000000 --- a/sysdeps/ieee754/flt-32/math_private.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef _MATH_PRIVATE_H_ - -#include_next <math_private.h> - -#ifndef __isnanf -extern __always_inline int -__isnanf (float d) -{ - u_int32_t di; - GET_FLOAT_WORD (di, d); - return (di & 0x7fffffff) > 0x7f800000; -} -#endif - -#ifndef __isinf_nsf -extern __always_inline int -__isinf_nsf (float d) -{ - u_int32_t di; - GET_FLOAT_WORD (di, d); - return (di & 0x7fffffff) == 0x7f800000; -} -#endif - -#ifndef __finitef -extern __always_inline int -__finitef (float d) -{ - u_int32_t di; - GET_FLOAT_WORD (di, d); - return (di & 0x7fffffff) < 0x7f800000; -} -#endif - -#endif /* _MATH_PRIVATE_H_ */ diff --git a/sysdeps/ieee754/flt-32/mpn2flt.c b/sysdeps/ieee754/flt-32/mpn2flt.c index 36a6dc3ad3..7f40f56d33 100644 --- a/sysdeps/ieee754/flt-32/mpn2flt.c +++ b/sysdeps/ieee754/flt-32/mpn2flt.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2015 Free Software Foundation, Inc. +/* Copyright (C) 1995-2016 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/flt-32/s_asinhf.c b/sysdeps/ieee754/flt-32/s_asinhf.c index 697faa843c..da9cafb600 100644 --- a/sysdeps/ieee754/flt-32/s_asinhf.c +++ b/sysdeps/ieee754/flt-32/s_asinhf.c @@ -30,11 +30,7 @@ __asinhf(float x) GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; if(__builtin_expect(ix< 0x38000000, 0)) { /* |x|<2**-14 */ - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if(huge+x>one) return x; /* return x inexact except 0 */ } if(__builtin_expect(ix>0x47000000, 0)) { /* |x| > 2**14 */ diff --git a/sysdeps/ieee754/flt-32/s_atanf.c b/sysdeps/ieee754/flt-32/s_atanf.c index be2addbdff..e322a1d41f 100644 --- a/sysdeps/ieee754/flt-32/s_atanf.c +++ b/sysdeps/ieee754/flt-32/s_atanf.c @@ -67,11 +67,7 @@ float __atanf(float x) else return -atanhi[3]-atanlo[3]; } if (ix < 0x3ee00000) { /* |x| < 0.4375 */ if (ix < 0x31000000) { /* |x| < 2^-29 */ - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if(huge+x>one) return x; /* raise inexact */ } id = -1; diff --git a/sysdeps/ieee754/flt-32/s_cbrtf.c b/sysdeps/ieee754/flt-32/s_cbrtf.c index 6d6c366c6f..b9e2b3738c 100644 --- a/sysdeps/ieee754/flt-32/s_cbrtf.c +++ b/sysdeps/ieee754/flt-32/s_cbrtf.c @@ -1,5 +1,5 @@ /* Compute cubic root of float value. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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. diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c index 864ab27a4f..0affd406bb 100644 --- a/sysdeps/ieee754/flt-32/s_cosf.c +++ b/sysdeps/ieee754/flt-32/s_cosf.c @@ -21,8 +21,6 @@ static char rcsid[] = "$NetBSD: s_cosf.c,v 1.4 1995/05/10 20:47:03 jtc Exp $"; #include <math.h> #include <math_private.h> -static const float one=1.0; - #ifndef COSF # define COSF_FUNC __cosf #else diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c index 2be44cc40c..c8b6287503 100644 --- a/sysdeps/ieee754/flt-32/s_erff.c +++ b/sysdeps/ieee754/flt-32/s_erff.c @@ -21,6 +21,7 @@ static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> static const float tiny = 1e-30, @@ -113,11 +114,7 @@ float __erff(float x) { /* Avoid spurious underflow. */ float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x); - if (fabsf (ret) < FLT_MIN) - { - float force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); return ret; } return x + efx*x; @@ -165,7 +162,10 @@ float __erfcf(float x) ix = hx&0x7fffffff; if(ix>=0x7f800000) { /* erfc(nan)=nan */ /* erfc(+-inf)=0,2 */ - return (float)(((u_int32_t)hx>>31)<<1)+one/x; + float ret = (float)(((u_int32_t)hx>>31)<<1)+one/x; + if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0f) + return 0.0f; + return ret; } if(ix < 0x3f580000) { /* |x|<0.84375 */ @@ -213,10 +213,7 @@ float __erfcf(float x) r = __ieee754_expf(-z*z-(float)0.5625)* __ieee754_expf((z-x)*(z+x)+R/S); if(hx>0) { -#if FLT_EVAL_METHOD != 0 - volatile -#endif - float ret = r/x; + float ret = math_narrow_eval (r/x); if (ret == 0) __set_errno (ERANGE); return ret; diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c index c81b057f24..c515d25e28 100644 --- a/sysdeps/ieee754/flt-32/s_expm1f.c +++ b/sysdeps/ieee754/flt-32/s_expm1f.c @@ -81,11 +81,7 @@ __expm1f(float x) c = (hi-x)-lo; } else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); t = huge+x; /* return x with inexact flags when x!=0 */ return x - (t-(huge+x)); } diff --git a/sysdeps/ieee754/flt-32/s_finitef.c b/sysdeps/ieee754/flt-32/s_finitef.c index 4ea270ae07..4c5b339235 100644 --- a/sysdeps/ieee754/flt-32/s_finitef.c +++ b/sysdeps/ieee754/flt-32/s_finitef.c @@ -35,7 +35,7 @@ int FINITEF(float x) { int32_t ix; GET_FLOAT_WORD(ix,x); - return (int)((u_int32_t)((ix&0x7fffffff)-0x7f800000)>>31); + return (int)((u_int32_t)((ix&0x7f800000)-0x7f800000)>>31); } hidden_def (__finitef) weak_alias (__finitef, finitef) diff --git a/sysdeps/ieee754/flt-32/s_fpclassifyf.c b/sysdeps/ieee754/flt-32/s_fpclassifyf.c index 8053184ae7..480557e2ae 100644 --- a/sysdeps/ieee754/flt-32/s_fpclassifyf.c +++ b/sysdeps/ieee754/flt-32/s_fpclassifyf.c @@ -1,5 +1,5 @@ /* Return classification value corresponding to argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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/flt-32/s_isinf_nsf.c b/sysdeps/ieee754/flt-32/s_isinf_nsf.c deleted file mode 100644 index 284d61926a..0000000000 --- a/sysdeps/ieee754/flt-32/s_isinf_nsf.c +++ /dev/null @@ -1,20 +0,0 @@ -/* - * Written by Ulrich Drepper <drepper@gmail.com>. - */ - -/* - * __isinf_nsf(x) returns != 0 if x is ±inf, else 0; - * no branching! - */ - -#include <math.h> -#include <math_private.h> - -#undef __isinf_nsf -int -__isinf_nsf (float x) -{ - int32_t ix; - GET_FLOAT_WORD(ix,x); - return (ix & 0x7fffffff) == 0x7f800000; -} diff --git a/sysdeps/ieee754/flt-32/s_issignalingf.c b/sysdeps/ieee754/flt-32/s_issignalingf.c index 093630c12e..2409ff408c 100644 --- a/sysdeps/ieee754/flt-32/s_issignalingf.c +++ b/sysdeps/ieee754/flt-32/s_issignalingf.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 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/flt-32/s_llrintf.c b/sysdeps/ieee754/flt-32/s_llrintf.c index ac8c1a2216..56415d34f8 100644 --- a/sysdeps/ieee754/flt-32/s_llrintf.c +++ b/sysdeps/ieee754/flt-32/s_llrintf.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,9 +18,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> static const float two23[2] = { @@ -34,7 +37,7 @@ __llrintf (float x) { int32_t j0; u_int32_t i0; - volatile float w; + float w; float t; long long int result; int sx; @@ -52,7 +55,7 @@ __llrintf (float x) result = (long long int) i0 << (j0 - 23); else { - w = two23[sx] + x; + w = math_narrow_eval (two23[sx] + x); t = w - two23[sx]; GET_FLOAT_WORD (i0, t); j0 = ((i0 >> 23) & 0xff) - 0x7f; @@ -64,8 +67,16 @@ __llrintf (float x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_FLT_LLONG_CONVERT_OVERFLOW && x != (float) LLONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LLONG_MAX : LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/flt-32/s_llroundf.c b/sysdeps/ieee754/flt-32/s_llroundf.c index f695c63710..1d35f71321 100644 --- a/sysdeps/ieee754/flt-32/s_llroundf.c +++ b/sysdeps/ieee754/flt-32/s_llroundf.c @@ -1,5 +1,5 @@ /* Round float value to long long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -17,9 +17,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> long long int @@ -51,8 +54,16 @@ __llroundf (float x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LLONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_FLT_LLONG_CONVERT_OVERFLOW && x != (float) LLONG_MIN) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LLONG_MAX : LLONG_MIN; + } +#endif return (long long int) x; } diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index 83a09f1414..ade60a2e27 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -50,11 +50,7 @@ __log1pf(float x) math_force_eval(two25+x); /* raise inexact */ if (ax<0x24800000) /* |x| < 2**-54 */ { - if (fabsf (x) < FLT_MIN) - { - float force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; } else diff --git a/sysdeps/ieee754/flt-32/s_logbf.c b/sysdeps/ieee754/flt-32/s_logbf.c index ba0267ebcb..9ae20e332a 100644 --- a/sysdeps/ieee754/flt-32/s_logbf.c +++ b/sysdeps/ieee754/flt-32/s_logbf.c @@ -15,6 +15,7 @@ #include <math.h> #include <math_private.h> +#include <fix-int-fp-convert-zero.h> float __logbf (float x) @@ -33,6 +34,8 @@ __logbf (float x) though it were normalized. */ rix -= __builtin_clz (ix) - 9; } + if (FIX_INT_FP_CONVERT_ZERO && rix == 127) + return 0.0f; return (float) (rix - 127); } weak_alias (__logbf, logbf) diff --git a/sysdeps/ieee754/flt-32/s_lrintf.c b/sysdeps/ieee754/flt-32/s_lrintf.c index 7581a8d286..3b480a2107 100644 --- a/sysdeps/ieee754/flt-32/s_lrintf.c +++ b/sysdeps/ieee754/flt-32/s_lrintf.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,9 +18,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> static const float two23[2] = { @@ -34,7 +37,7 @@ __lrintf (float x) { int32_t j0; u_int32_t i0; - volatile float w; + float w; float t; long int result; int sx; @@ -52,7 +55,7 @@ __lrintf (float x) result = (long int) i0 << (j0 - 23); else { - w = two23[sx] + x; + w = math_narrow_eval (two23[sx] + x); t = w - two23[sx]; GET_FLOAT_WORD (i0, t); j0 = ((i0 >> 23) & 0xff) - 0x7f; @@ -64,8 +67,16 @@ __lrintf (float x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_FLT_LONG_CONVERT_OVERFLOW && x != (float) LONG_MIN) + { + feraiseexcept (FE_INVALID); + return sx == 0 ? LONG_MAX : LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/flt-32/s_lroundf.c b/sysdeps/ieee754/flt-32/s_lroundf.c index f20352a3c8..116c9e0627 100644 --- a/sysdeps/ieee754/flt-32/s_lroundf.c +++ b/sysdeps/ieee754/flt-32/s_lroundf.c @@ -1,5 +1,5 @@ /* Round float value to long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -17,9 +17,12 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> +#include <fix-fp-int-convert-overflow.h> long int @@ -51,8 +54,16 @@ __lroundf (float x) } else { - /* The number is too large. It is left implementation defined - what happens. */ +#ifdef FE_INVALID + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ + if (FIX_FLT_LONG_CONVERT_OVERFLOW && x != (float) LONG_MIN) + { + feraiseexcept (FE_INVALID); + return sign == 1 ? LONG_MAX : LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/flt-32/s_nextafterf.c b/sysdeps/ieee754/flt-32/s_nextafterf.c index 22e0b3d4ed..625d54b768 100644 --- a/sysdeps/ieee754/flt-32/s_nextafterf.c +++ b/sysdeps/ieee754/flt-32/s_nextafterf.c @@ -17,6 +17,7 @@ static char rcsid[] = "$NetBSD: s_nextafterf.c,v 1.4 1995/05/10 20:48:01 jtc Exp $"; #endif +#include <errno.h> #include <math.h> #include <math_private.h> #include <float.h> @@ -59,10 +60,12 @@ float __nextafterf(float x, float y) if(hy>=0x7f800000) { float u = x+x; /* overflow */ math_force_eval (u); + __set_errno (ERANGE); } if(hy<0x00800000) { float u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } SET_FLOAT_WORD(x,hx); return x; diff --git a/sysdeps/ieee754/flt-32/s_remquof.c b/sysdeps/ieee754/flt-32/s_remquof.c index 36cf359b9d..ecf831deaf 100644 --- a/sysdeps/ieee754/flt-32/s_remquof.c +++ b/sysdeps/ieee754/flt-32/s_remquof.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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/flt-32/s_roundf.c b/sysdeps/ieee754/flt-32/s_roundf.c index 2c2a9e198d..a75d98f384 100644 --- a/sysdeps/ieee754/flt-32/s_roundf.c +++ b/sysdeps/ieee754/flt-32/s_roundf.c @@ -1,5 +1,5 @@ /* Round float to integer away from zero. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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/flt-32/s_scalbnf.c b/sysdeps/ieee754/flt-32/s_scalbnf.c index 5d9bdd647b..f36ae241b2 100644 --- a/sysdeps/ieee754/flt-32/s_scalbnf.c +++ b/sysdeps/ieee754/flt-32/s_scalbnf.c @@ -50,4 +50,3 @@ __scalbnf (float x, int n) SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x*twom25; } -weak_alias (__scalbnf, scalbnf) diff --git a/sysdeps/ieee754/flt-32/s_signbitf.c b/sysdeps/ieee754/flt-32/s_signbitf.c index 169820e9df..3cdde778d5 100644 --- a/sysdeps/ieee754/flt-32/s_signbitf.c +++ b/sysdeps/ieee754/flt-32/s_signbitf.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -19,13 +19,8 @@ #include <math.h> -#include <math_private.h> - int __signbitf (float x) { - int32_t hx; - - GET_FLOAT_WORD (hx, x); - return hx & 0x80000000; + return __builtin_signbitf (x); } diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index 596e076609..e0737b5ce4 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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/flt-32/s_tanhf.c b/sysdeps/ieee754/flt-32/s_tanhf.c index dc96da9a5b..f70702b29c 100644 --- a/sysdeps/ieee754/flt-32/s_tanhf.c +++ b/sysdeps/ieee754/flt-32/s_tanhf.c @@ -17,6 +17,7 @@ static char rcsid[] = "$NetBSD: s_tanhf.c,v 1.4 1995/05/10 20:48:24 jtc Exp $"; #endif +#include <float.h> #include <math.h> #include <math_private.h> @@ -41,7 +42,10 @@ float __tanhf(float x) if (ix == 0) return x; /* x == +-0 */ if (ix<0x24000000) /* |x|<2**-55 */ + { + math_check_force_underflow (x); return x*(one+x); /* tanh(small) = small */ + } if (ix>=0x3f800000) { /* |x|>=1 */ t = __expm1f(two*fabsf(x)); z = one - two/(t+two); diff --git a/sysdeps/ieee754/flt-32/s_truncf.c b/sysdeps/ieee754/flt-32/s_truncf.c index 1c4f9c94ed..43d35c7f6a 100644 --- a/sysdeps/ieee754/flt-32/s_truncf.c +++ b/sysdeps/ieee754/flt-32/s_truncf.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 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/flt-32/t_exp2f.h b/sysdeps/ieee754/flt-32/t_exp2f.h index db528038f1..3045b82cd1 100644 --- a/sysdeps/ieee754/flt-32/t_exp2f.h +++ b/sysdeps/ieee754/flt-32/t_exp2f.h @@ -1,5 +1,5 @@ /* Accurate tables for exp2f(). - Copyright (C) 1998-2015 Free Software Foundation, Inc. + Copyright (C) 1998-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> diff --git a/sysdeps/ieee754/flt-32/w_expf.c b/sysdeps/ieee754/flt-32/w_expf.c index cc5ff76421..ed1550972f 100644 --- a/sysdeps/ieee754/flt-32/w_expf.c +++ b/sysdeps/ieee754/flt-32/w_expf.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2015 Free Software Foundation, Inc. +/* Copyright (C) 2011-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. |