diff options
Diffstat (limited to 'sysdeps/ieee754/flt-32')
87 files changed, 2122 insertions, 1243 deletions
diff --git a/sysdeps/ieee754/flt-32/e_acosf.c b/sysdeps/ieee754/flt-32/e_acosf.c index 6f792f6604..3b2e4f1bde 100644 --- a/sysdeps/ieee754/flt-32/e_acosf.c +++ b/sysdeps/ieee754/flt-32/e_acosf.c @@ -56,14 +56,14 @@ __ieee754_acosf(float x) z = (one+x)*(float)0.5; p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - s = __ieee754_sqrtf(z); + s = sqrtf(z); r = p/q; w = r*s-pio2_lo; return pi - (float)2.0*(s+w); } else { /* x > 0.5 */ int32_t idf; z = (one-x)*(float)0.5; - s = __ieee754_sqrtf(z); + s = sqrtf(z); df = s; GET_FLOAT_WORD(idf,df); SET_FLOAT_WORD(df,idf&0xfffff000); diff --git a/sysdeps/ieee754/flt-32/e_acoshf.c b/sysdeps/ieee754/flt-32/e_acoshf.c index aabfb85df7..49e64f3c43 100644 --- a/sysdeps/ieee754/flt-32/e_acoshf.c +++ b/sysdeps/ieee754/flt-32/e_acoshf.c @@ -40,10 +40,10 @@ float __ieee754_acoshf(float x) return 0.0; /* acosh(1) = 0 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ t=x*x; - return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one))); + return __ieee754_logf((float)2.0*x-one/(x+sqrtf(t-one))); } else { /* 1<x<2 */ t = x-one; - return __log1pf(t+__ieee754_sqrtf((float)2.0*t+t*t)); + return __log1pf(t+sqrtf((float)2.0*t+t*t)); } } strong_alias (__ieee754_acoshf, __acoshf_finite) diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c index 2ca2dbcb28..e03073b8bd 100644 --- a/sysdeps/ieee754/flt-32/e_asinf.c +++ b/sysdeps/ieee754/flt-32/e_asinf.c @@ -42,6 +42,7 @@ static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const float one = 1.0000000000e+00, /* 0x3F800000 */ @@ -85,7 +86,7 @@ float __ieee754_asinf(float x) w = one-fabsf(x); t = w*0.5f; p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); - s = __ieee754_sqrtf(t); + s = sqrtf(t); if(ix>=0x3F79999A) { /* if |x| > 0.975 */ t = pio2_hi-(2.0f*(s+s*p)-pio2_lo); } else { diff --git a/sysdeps/ieee754/flt-32/e_atan2f.c b/sysdeps/ieee754/flt-32/e_atan2f.c index 29eefc0dd6..ddc5873ade 100644 --- a/sysdeps/ieee754/flt-32/e_atan2f.c +++ b/sysdeps/ieee754/flt-32/e_atan2f.c @@ -81,7 +81,7 @@ __ieee754_atan2f (float y, float x) switch (m) { case 0: return z ; /* atan(+,+) */ case 1: { - u_int32_t zh; + uint32_t zh; GET_FLOAT_WORD(zh,z); SET_FLOAT_WORD(z,zh ^ 0x80000000); } diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c index 82de071ab8..207d759362 100644 --- a/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/sysdeps/ieee754/flt-32/e_atanhf.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2016 Free Software Foundation, Inc. +/* Copyright (C) 2011-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -38,7 +38,9 @@ #include <float.h> #include <inttypes.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> static const float huge = 1e30; diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c index 7b223758e1..a2aa83876d 100644 --- a/sysdeps/ieee754/flt-32/e_coshf.c +++ b/sysdeps/ieee754/flt-32/e_coshf.c @@ -15,6 +15,7 @@ */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> static const float huge = 1.0e30; diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 1723c482de..7218e5d254 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -1,7 +1,6 @@ -/* Single-precision floating point 2^x. - Copyright (C) 1997-2016 Free Software Foundation, Inc. +/* Single-precision 2^x function. + Copyright (C) 2017-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. - Contributed by Geoffrey Keating <geoffk@ozemail.com.au> The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -17,116 +16,81 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ -/* The basic design here is from - Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical - Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft., - 17 (1), March 1991, pp. 26-45. - It has been slightly modified to compute 2^x instead of e^x, and for - single-precision. - */ -#ifndef _GNU_SOURCE -# define _GNU_SOURCE -#endif -#include <stdlib.h> -#include <float.h> -#include <ieee754.h> #include <math.h> -#include <fenv.h> -#include <inttypes.h> -#include <math_private.h> - -#include "t_exp2f.h" - -static const float TWOM100 = 7.88860905e-31; -static const float TWO127 = 1.7014118346e+38; +#include <math-narrow-eval.h> +#include <stdint.h> +#include <shlib-compat.h> +#include <libm-alias-float.h> +#include "math_config.h" + +/* +EXP2F_TABLE_BITS = 5 +EXP2F_POLY_ORDER = 3 + +ULP error: 0.502 (nearest rounding.) +Relative error: 1.69 * 2^-34 in [-1/64, 1/64] (before rounding.) +Wrong count: 168353 (all nearest rounding wrong results with fma.) +Non-nearest ULP error: 1 (rounded ULP error) +*/ + +#define N (1 << EXP2F_TABLE_BITS) +#define T __exp2f_data.tab +#define C __exp2f_data.poly +#define SHIFT __exp2f_data.shift_scaled + +static inline uint32_t +top12 (float x) +{ + return asuint (x) >> 20; +} float -__ieee754_exp2f (float x) +__exp2f (float x) { - static const float himark = (float) FLT_MAX_EXP; - static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1); - - /* Check for usual case. */ - if (isless (x, himark) && isgreaterequal (x, lomark)) - { - static const float THREEp14 = 49152.0; - int tval, unsafe; - float rx, x22, result; - union ieee754_float ex2_u, scale_u; - - if (fabsf (x) < FLT_EPSILON / 4.0f) - return 1.0f + x; - - { - SET_RESTORE_ROUND_NOEXF (FE_TONEAREST); - - /* 1. Argument reduction. - Choose integers ex, -128 <= t < 128, and some real - -1/512 <= x1 <= 1/512 so that - x = ex + t/512 + x1. - - First, calculate rx = ex + t/256. */ - rx = x + THREEp14; - rx -= THREEp14; - x -= rx; /* Compute x=x1. */ - /* Compute tval = (ex*256 + t)+128. - Now, t = (tval mod 256)-128 and ex=tval/256 [that's mod, NOT %; - and /-round-to-nearest not the usual c integer /]. */ - tval = (int) (rx * 256.0f + 128.0f); - - /* 2. Adjust for accurate table entry. - Find e so that - x = ex + t/256 + e + x2 - where -7e-4 < e < 7e-4, and - (float)(2^(t/256+e)) - is accurate to one part in 2^-64. */ - - /* 'tval & 255' is the same as 'tval%256' except that it's always - positive. - Compute x = x2. */ - x -= __exp2f_deltatable[tval & 255]; - - /* 3. Compute ex2 = 2^(t/255+e+ex). */ - ex2_u.f = __exp2f_atable[tval & 255]; - tval >>= 8; - /* x2 is an integer multiple of 2^-30; avoid intermediate - underflow from the calculation of x22 * x. */ - unsafe = abs(tval) >= -FLT_MIN_EXP - 32; - ex2_u.ieee.exponent += tval >> unsafe; - scale_u.f = 1.0; - scale_u.ieee.exponent += tval - (tval >> unsafe); - - /* 4. Approximate 2^x2 - 1, using a second-degree polynomial, - with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14] - less than 1.3e-10. */ - - x22 = (.24022656679f * x + .69314736128f) * ex2_u.f; - } - - /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */ - result = x22 * x + ex2_u.f; - - if (!unsafe) - return result; - else - { - result *= scale_u.f; - math_check_force_underflow_nonneg (result); - return result; - } - } - /* Exceptional cases: */ - else if (isless (x, himark)) + uint32_t abstop; + uint64_t ki, t; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t kd, xd, z, r, r2, y, s; + + xd = (double_t) x; + abstop = top12 (x) & 0x7ff; + if (__glibc_unlikely (abstop >= top12 (128.0f))) { - if (isinf (x)) - /* e^-inf == 0, with no error. */ - return 0; - else - /* Underflow */ - return TWOM100 * TWOM100; + /* |x| >= 128 or x is nan. */ + if (asuint (x) == asuint (-INFINITY)) + return 0.0f; + if (abstop >= top12 (INFINITY)) + return x + x; + if (x > 0.0f) + return __math_oflowf (0); + if (x <= -150.0f) + return __math_uflowf (0); +#if WANT_ERRNO_UFLOW + if (x < -149.0f) + return __math_may_uflowf (0); +#endif } - else - /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ - return TWO127*x; + + /* x = k/N + r with r in [-1/(2N), 1/(2N)] and int k. */ + kd = math_narrow_eval ((double) (xd + SHIFT)); /* Needs to be double. */ + ki = asuint64 (kd); + kd -= SHIFT; /* k/N for int k. */ + r = xd - kd; + + /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return (float) y; } -strong_alias (__ieee754_exp2f, __exp2f_finite) +#ifndef __exp2f +strong_alias (__exp2f, __ieee754_exp2f) +strong_alias (__exp2f, __exp2f_finite) +versioned_symbol (libm, __exp2f, exp2f, GLIBC_2_27); +libm_alias_float_other (__exp2, exp2) +#endif diff --git a/sysdeps/ieee754/flt-32/e_exp2f_data.c b/sysdeps/ieee754/flt-32/e_exp2f_data.c new file mode 100644 index 0000000000..8d4e702c92 --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_exp2f_data.c @@ -0,0 +1,44 @@ +/* Shared data between expf, exp2f and powf. + Copyright (C) 2017-2018 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_config.h" + +#define N (1 << EXP2F_TABLE_BITS) + +const struct exp2f_data __exp2f_data = { + /* tab[i] = uint(2^(i/N)) - (i << 52-BITS) + used for computing 2^(k/N) for an int |k| < 150 N as + double(tab[k%N] + (k << 52-BITS)) */ + .tab = { +0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, 0x3fef9301d0125b51, +0x3fef72b83c7d517b, 0x3fef54873168b9aa, 0x3fef387a6e756238, 0x3fef1e9df51fdee1, +0x3fef06fe0a31b715, 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d, +0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429, 0x3feea47eb03a5585, +0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, 0x3feea11473eb0187, 0x3feea589994cce13, +0x3feeace5422aa0db, 0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d, +0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, 0x3fef3720dcef9069, +0x3fef5818dcfba487, 0x3fef7c97337b9b5f, 0x3fefa4afa2a490da, 0x3fefd0765b6e4540, + }, + .shift_scaled = 0x1.8p+52 / N, + .poly = { 0x1.c6af84b912394p-5, 0x1.ebfce50fac4f3p-3, 0x1.62e42ff0c52d6p-1 }, + .shift = 0x1.8p+52, + .invln2_scaled = 0x1.71547652b82fep+0 * N, + .poly_scaled = { +0x1.c6af84b912394p-5/N/N/N, 0x1.ebfce50fac4f3p-3/N/N, 0x1.62e42ff0c52d6p-1/N + }, +}; diff --git a/sysdeps/ieee754/flt-32/e_expf.c b/sysdeps/ieee754/flt-32/e_expf.c index 071f615ef4..f2238bfd74 100644 --- a/sysdeps/ieee754/flt-32/e_expf.c +++ b/sysdeps/ieee754/flt-32/e_expf.c @@ -1,7 +1,6 @@ -/* Single-precision floating point e^x. - Copyright (C) 1997-2016 Free Software Foundation, Inc. +/* Single-precision e^x function. + Copyright (C) 2017-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. - Contributed by Geoffrey Keating <geoffk@ozemail.com.au> The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -17,117 +16,102 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ -/* How this works: +#ifdef __expf +# undef libm_hidden_proto +# define libm_hidden_proto(ignored) +#endif - The input value, x, is written as - - x = n * ln(2) + t/512 + delta[t] + x; - - where: - - n is an integer, 127 >= n >= -150; - - t is an integer, 177 >= t >= -177 - - delta is based on a table entry, delta[t] < 2^-28 - - x is whatever is left, |x| < 2^-10 - - Then e^x is approximated as - - e^x = 2^n ( e^(t/512 + delta[t]) - + ( e^(t/512 + delta[t]) - * ( p(x + delta[t] + n * ln(2)) - delta ) ) ) - - where - - p(x) is a polynomial approximating e(x)-1; - - e^(t/512 + delta[t]) is obtained from a table. - - The table used is the same one as for the double precision version; - since we have the table, we might as well use it. - - It turns out to be faster to do calculations in double precision than - to perform an 'accurate table method' expf, because of the range reduction - overhead (compare exp2f). - */ -#include <float.h> -#include <ieee754.h> #include <math.h> -#include <fenv.h> -#include <inttypes.h> -#include <math_private.h> - -extern const float __exp_deltatable[178]; -extern const double __exp_atable[355] /* __attribute__((mode(DF))) */; - -static const float TWOM100 = 7.88860905e-31; -static const float TWO127 = 1.7014118346e+38; +#include <math-narrow-eval.h> +#include <stdint.h> +#include <shlib-compat.h> +#include <libm-alias-float.h> +#include "math_config.h" + +/* +EXP2F_TABLE_BITS = 5 +EXP2F_POLY_ORDER = 3 + +ULP error: 0.502 (nearest rounding.) +Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.) +Wrong count: 170635 (all nearest rounding wrong results with fma.) +Non-nearest ULP error: 1 (rounded ULP error) +*/ + +#define N (1 << EXP2F_TABLE_BITS) +#define InvLn2N __exp2f_data.invln2_scaled +#define T __exp2f_data.tab +#define C __exp2f_data.poly_scaled + +static inline uint32_t +top12 (float x) +{ + return asuint (x) >> 20; +} float -__ieee754_expf (float x) +__expf (float x) { - static const float himark = 88.72283935546875; - static const float lomark = -103.972084045410; - /* Check for usual case. */ - if (isless (x, himark) && isgreater (x, lomark)) - { - static const float THREEp42 = 13194139533312.0; - static const float THREEp22 = 12582912.0; - /* 1/ln(2). */ -#undef M_1_LN2 - static const float M_1_LN2 = 1.44269502163f; - /* ln(2) */ -#undef M_LN2 - static const double M_LN2 = .6931471805599452862; - - int tval; - double x22, t, result, dx; - float n, delta; - union ieee754_double ex2_u; - - { - SET_RESTORE_ROUND_NOEXF (FE_TONEAREST); - - /* Calculate n. */ - n = x * M_1_LN2 + THREEp22; - n -= THREEp22; - dx = x - n*M_LN2; - - /* Calculate t/512. */ - t = dx + THREEp42; - t -= THREEp42; - dx -= t; - - /* Compute tval = t. */ - tval = (int) (t * 512.0); - - if (t >= 0) - delta = - __exp_deltatable[tval]; - else - delta = __exp_deltatable[-tval]; - - /* Compute ex2 = 2^n e^(t/512+delta[t]). */ - ex2_u.d = __exp_atable[tval+177]; - ex2_u.ieee.exponent += (int) n; - - /* Approximate e^(dx+delta) - 1, using a second-degree polynomial, - with maximum error in [-2^-10-2^-28,2^-10+2^-28] - less than 5e-11. */ - x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta; - } - - /* Return result. */ - result = x22 * ex2_u.d + ex2_u.d; - return (float) result; - } - /* Exceptional cases: */ - else if (isless (x, himark)) + uint32_t abstop; + uint64_t ki, t; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t kd, xd, z, r, r2, y, s; + + xd = (double_t) x; + abstop = top12 (x) & 0x7ff; + if (__glibc_unlikely (abstop >= top12 (88.0f))) { - if (isinf (x)) - /* e^-inf == 0, with no error. */ - return 0; - else - /* Underflow */ - return TWOM100 * TWOM100; + /* |x| >= 88 or x is nan. */ + if (asuint (x) == asuint (-INFINITY)) + return 0.0f; + if (abstop >= top12 (INFINITY)) + return x + x; + if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ + return __math_oflowf (0); + if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ + return __math_uflowf (0); +#if WANT_ERRNO_UFLOW + if (x < -0x1.9d1d9ep6f) /* x < log(0x1p-149) ~= -103.28 */ + return __math_may_uflowf (0); +#endif } - else - /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ - return TWO127*x; + + /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ + z = InvLn2N * xd; + + /* Round and convert z to int, the result is in [-150*N, 128*N] and + ideally ties-to-even rule is used, otherwise the magnitude of r + can be bigger which gives larger approximation error. */ +#if TOINT_INTRINSICS + kd = roundtoint (z); + ki = converttoint (z); +#elif TOINT_RINT + kd = rint (z); + ki = (long) kd; +#elif TOINT_SHIFT +# define SHIFT __exp2f_data.shift + kd = math_narrow_eval ((double) (z + SHIFT)); /* Needs to be double. */ + ki = asuint64 (kd); + kd -= SHIFT; +#endif + r = z - kd; + + /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return (float) y; } -strong_alias (__ieee754_expf, __expf_finite) + +#ifndef __expf +hidden_def (__expf) +strong_alias (__expf, __ieee754_expf) +strong_alias (__expf, __expf_finite) +versioned_symbol (libm, __expf, expf, GLIBC_2_27); +libm_alias_float_other (__exp, exp) +#endif diff --git a/sysdeps/ieee754/flt-32/e_fmodf.c b/sysdeps/ieee754/flt-32/e_fmodf.c index 8d8fad4eb5..1a9407b517 100644 --- a/sysdeps/ieee754/flt-32/e_fmodf.c +++ b/sysdeps/ieee754/flt-32/e_fmodf.c @@ -41,7 +41,7 @@ __ieee754_fmodf (float x, float y) return (x*y)/(x*y); if(hx<hy) return x; /* |x|<|y| return x */ if(hx==hy) - return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ + return Zero[(uint32_t)sx>>31]; /* |x|=|y| return x*0*/ /* determine ix = ilogb(x) */ if(hx<0x00800000) { /* subnormal x */ @@ -74,7 +74,7 @@ __ieee754_fmodf (float x, float y) if(hz<0){hx = hx+hx;} else { if(hz==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; + return Zero[(uint32_t)sx>>31]; hx = hz+hz; } } @@ -83,7 +83,7 @@ __ieee754_fmodf (float x, float y) /* convert back to floating value and restore the sign */ if(hx==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; + return Zero[(uint32_t)sx>>31]; while(hx<0x00800000) { /* normalize x */ hx = hx+hx; iy -= 1; diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c index 19f51b0c8b..8b23add347 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,7 +18,9 @@ <http://www.gnu.org/licenses/>. */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> #include <float.h> /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's @@ -91,7 +93,7 @@ gammaf_positive (float x, int *exp2_adj) float ret = (__ieee754_powf (x_adj_mant, x_adj) * __ieee754_exp2f (x_adj_log2 * x_adj_frac) * __ieee754_expf (-x_adj) - * __ieee754_sqrtf (2 * (float) M_PI / x_adj) + * sqrtf (2 * (float) M_PI / x_adj) / prod); exp_adj += x_eps * __ieee754_logf (x_adj); float bsum = gamma_coeff[NCOEFF - 1]; @@ -118,7 +120,7 @@ __ieee754_gammaf_r (float x, int *signgamp) return 1.0 / x; } if (__builtin_expect (hx < 0, 0) - && (u_int32_t) hx < 0xff800000 && __rintf (x) == x) + && (uint32_t) hx < 0xff800000 && __rintf (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c index 717b82e42f..5336876cf4 100644 --- a/sysdeps/ieee754/flt-32/e_hypotf.c +++ b/sysdeps/ieee754/flt-32/e_hypotf.c @@ -26,9 +26,9 @@ __ieee754_hypotf(float x, float y) ha &= 0x7fffffff; GET_FLOAT_WORD(hb,y); hb &= 0x7fffffff; - if (ha == 0x7f800000) + if (ha == 0x7f800000 && !issignaling (y)) return fabsf(x); - else if (hb == 0x7f800000) + else if (hb == 0x7f800000 && !issignaling (x)) return fabsf(y); else if (ha > 0x7f800000 || hb > 0x7f800000) return fabsf(x) * fabsf(y); @@ -40,6 +40,6 @@ __ieee754_hypotf(float x, float y) d_x = (double) x; d_y = (double) y; - return (float) __ieee754_sqrt(d_x * d_x + d_y * d_y); + return (float) sqrt(d_x * d_x + d_y * d_y); } strong_alias (__ieee754_hypotf, __hypotf_finite) diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c index bd0b80fdb0..0efc646a12 100644 --- a/sysdeps/ieee754/flt-32/e_j0f.c +++ b/sysdeps/ieee754/flt-32/e_j0f.c @@ -14,6 +14,7 @@ */ #include <math.h> +#include <math-barriers.h> #include <math_private.h> static float pzerof(float), qzerof(float); @@ -58,10 +59,10 @@ __ieee754_j0f(float x) * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x); + if(ix>0x48000000) z = (invsqrtpi*cc)/sqrtf(x); else { u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x); + z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); } return z; } @@ -105,7 +106,7 @@ __ieee754_y0f(float x) ix = 0x7fffffff&hx; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */ if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ + if(ix==0) return -1/zero; /* -inf and divide by zero exception. */ if(hx<0) return zero/(zero*x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) @@ -131,10 +132,10 @@ __ieee754_y0f(float x) if ((s*c)<zero) cc = z/ss; else ss = z/cc; } - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); + if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x); else { u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); + z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c index f359a3d9ba..887b5f7d1e 100644 --- a/sysdeps/ieee754/flt-32/e_j1f.c +++ b/sysdeps/ieee754/flt-32/e_j1f.c @@ -16,7 +16,9 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static float ponef(float), qonef(float); @@ -61,10 +63,10 @@ __ieee754_j1f(float x) * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ - if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(y); + if(ix>0x48000000) z = (invsqrtpi*cc)/sqrtf(y); else { u = ponef(y); v = qonef(y); - z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(y); + z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); } if(hx<0) return -z; else return z; @@ -112,7 +114,7 @@ __ieee754_y1f(float x) /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ if(__builtin_expect(ix>=0x7f800000, 0)) return one/(x+x*x); if(__builtin_expect(ix==0, 0)) - return -HUGE_VALF+x; /* -inf and overflow exception. */ + return -1/zero; /* -inf and divide by zero exception. */ if(__builtin_expect(hx<0, 0)) return zero/(zero*x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ SET_RESTORE_ROUNDF (FE_TONEAREST); @@ -135,10 +137,10 @@ __ieee754_y1f(float x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); + if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x); else { u = ponef(x); v = qonef(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); + z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c index 4e634778d3..cd15ed7d4b 100644 --- a/sysdeps/ieee754/flt-32/e_jnf.c +++ b/sysdeps/ieee754/flt-32/e_jnf.c @@ -16,7 +16,9 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static const float two = 2.0000000000e+00, /* 0x40000000 */ @@ -186,7 +188,7 @@ __ieee754_ynf(int n, float x) float ret; { int32_t i,hx,ix; - u_int32_t ib; + uint32_t ib; int32_t sign; float a, b, temp; @@ -194,15 +196,15 @@ __ieee754_ynf(int n, float x) ix = 0x7fffffff&hx; /* if Y(n,NaN) is NaN */ if(__builtin_expect(ix>0x7f800000, 0)) return x+x; - if(__builtin_expect(ix==0, 0)) - return -HUGE_VALF+x; /* -inf and overflow exception. */ - if(__builtin_expect(hx<0, 0)) return zero/(zero*x); sign = 1; if(n<0){ n = -n; sign = 1 - ((n&1)<<1); } if(n==0) return(__ieee754_y0f(x)); + if(__builtin_expect(ix==0, 0)) + return -sign/zero; + if(__builtin_expect(hx<0, 0)) return zero/(zero*x); SET_RESTORE_ROUNDF (FE_TONEAREST); if(n==1) { ret = sign*__ieee754_y1f(x); diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c index 4d8a66bb39..8fdf9bb8bc 100644 --- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c +++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c @@ -13,9 +13,10 @@ * ==================================================== */ -#include <libc-internal.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libc-diag.h> static const float two23= 8.3886080000e+06, /* 0x4b000000 */ @@ -160,7 +161,7 @@ __ieee754_lgammaf_r(float x, int *signgamp) } if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ - return x/zero; + return fabsf (x)/zero; if (ix > 0x40000000 /* X < 2.0f. */ && ix < 0x41700000 /* X > -15.0f. */) return __lgamma_negf (x, signgamp); diff --git a/sysdeps/ieee754/flt-32/e_log10f.c b/sysdeps/ieee754/flt-32/e_log10f.c index 2cd01b4a50..7f1ffdad77 100644 --- a/sysdeps/ieee754/flt-32/e_log10f.c +++ b/sysdeps/ieee754/flt-32/e_log10f.c @@ -34,7 +34,7 @@ __ieee754_log10f(float x) k=0; if (hx < 0x00800000) { /* x < 2**-126 */ if (__builtin_expect((hx&0x7fffffff)==0, 0)) - return -two25/(x-x); /* log(+-0)=-inf */ + return -two25/fabsf (x); /* log(+-0)=-inf */ if (__builtin_expect(hx<0, 0)) return (x-x)/(x-x); /* log(-#) = NaN */ k -= 25; x *= two25; /* subnormal number, scale up x */ @@ -42,7 +42,7 @@ __ieee754_log10f(float x) } if (__builtin_expect(hx >= 0x7f800000, 0)) return x+x; k += (hx>>23)-127; - i = ((u_int32_t)k&0x80000000)>>31; + i = ((uint32_t)k&0x80000000)>>31; hx = (hx&0x007fffff)|((0x7f-i)<<23); y = (float)(k+i); if (FIX_INT_FP_CONVERT_ZERO && y == 0.0f) diff --git a/sysdeps/ieee754/flt-32/e_log2f.c b/sysdeps/ieee754/flt-32/e_log2f.c index 857d13fb9b..62edc59cde 100644 --- a/sysdeps/ieee754/flt-32/e_log2f.c +++ b/sysdeps/ieee754/flt-32/e_log2f.c @@ -1,86 +1,95 @@ -/* e_logf.c -- float version of e_log.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - * adapted for log2 by Ulrich Drepper <drepper@cygnus.com> - */ +/* Single-precision log2 function. + Copyright (C) 2017-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * ==================================================== - * 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. - * ==================================================== - */ + 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 <fix-int-fp-convert-zero.h> +#include <stdint.h> +#include <shlib-compat.h> +#include <libm-alias-float.h> +#include "math_config.h" + +/* +LOG2F_TABLE_BITS = 4 +LOG2F_POLY_ORDER = 4 -static const float -ln2 = 0.69314718055994530942, -two25 = 3.355443200e+07, /* 0x4c000000 */ -Lg1 = 6.6666668653e-01, /* 3F2AAAAB */ -Lg2 = 4.0000000596e-01, /* 3ECCCCCD */ -Lg3 = 2.8571429849e-01, /* 3E924925 */ -Lg4 = 2.2222198546e-01, /* 3E638E29 */ -Lg5 = 1.8183572590e-01, /* 3E3A3325 */ -Lg6 = 1.5313838422e-01, /* 3E1CD04F */ -Lg7 = 1.4798198640e-01; /* 3E178897 */ +ULP error: 0.752 (nearest rounding.) +Relative error: 1.9 * 2^-26 (before rounding.) +*/ -static const float zero = 0.0; +#define N (1 << LOG2F_TABLE_BITS) +#define T __log2f_data.tab +#define A __log2f_data.poly +#define OFF 0x3f330000 float -__ieee754_log2f(float x) +__log2f (float x) { - float hfsq,f,s,z,R,w,t1,t2,dk; - int32_t k,ix,i,j; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, p, y, y0, invc, logc; + uint32_t ix, iz, top, tmp; + int k, i; + + ix = asuint (x); +#if WANT_ROUNDING + /* Fix sign of zero with downward rounding when x==1. */ + if (__glibc_unlikely (ix == 0x3f800000)) + return 0; +#endif + if (__glibc_unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000)) + { + /* x < 0x1p-126 or inf or nan. */ + if (ix * 2 == 0) + return __math_divzerof (1); + if (ix == 0x7f800000) /* log2(inf) == inf. */ + return x; + if ((ix & 0x80000000) || ix * 2 >= 0xff000000) + return __math_invalidf (x); + /* x is subnormal, normalize it. */ + ix = asuint (x * 0x1p23f); + ix -= 23 << 23; + } + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N; + top = tmp & 0xff800000; + iz = ix - top; + k = (int32_t) tmp >> 23; /* arithmetic shift */ + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); - GET_FLOAT_WORD(ix,x); + /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ + r = z * invc - 1; + y0 = logc + (double_t) k; - k=0; - if (ix < 0x00800000) { /* x < 2**-126 */ - if (__builtin_expect((ix&0x7fffffff)==0, 0)) - return -two25/(x-x); /* log(+-0)=-inf */ - if (__builtin_expect(ix<0, 0)) - return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ - GET_FLOAT_WORD(ix,x); - } - if (__builtin_expect(ix >= 0x7f800000, 0)) return x+x; - k += (ix>>23)-127; - ix &= 0x007fffff; - i = (ix+(0x95f64<<3))&0x800000; - SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */ - k += (i>>23); - dk = (float)k; - f = x-(float)1.0; - if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */ - 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; - } - s = f/((float)2.0+f); - z = s*s; - i = ix-(0x6147a<<3); - w = z*z; - j = (0x6b851<<3)-ix; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - i |= j; - R = t2+t1; - if(i>0) { - hfsq=(float)0.5*f*f; - return dk-((hfsq-(s*(hfsq+R)))-f)/ln2; - } else { - return dk-((s*(f-R))-f)/ln2; - } + /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ + r2 = r * r; + y = A[1] * r + A[2]; + y = A[0] * r2 + y; + p = A[3] * r + y0; + y = y * r2 + p; + return (float) y; } -strong_alias (__ieee754_log2f, __log2f_finite) +#ifndef __log2f +strong_alias (__log2f, __ieee754_log2f) +strong_alias (__log2f, __log2f_finite) +versioned_symbol (libm, __log2f, log2f, GLIBC_2_27); +libm_alias_float_other (__log2, log2) +#endif diff --git a/sysdeps/ieee754/flt-32/e_log2f_data.c b/sysdeps/ieee754/flt-32/e_log2f_data.c new file mode 100644 index 0000000000..6390e9ce51 --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_log2f_data.c @@ -0,0 +1,44 @@ +/* Data definition for log2f. + Copyright (C) 2017-2018 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_config.h" + +const struct log2f_data __log2f_data = { + .tab = { + { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 }, + { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 }, + { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 }, + { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 }, + { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 }, + { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 }, + { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 }, + { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 }, + { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 }, + { 0x1p+0, 0x0p+0 }, + { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 }, + { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 }, + { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 }, + { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 }, + { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 }, + { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 }, + }, + .poly = { + -0x1.712b6f70a7e4dp-2, 0x1.ecabf496832ep-2, -0x1.715479ffae3dep-1, + 0x1.715475f35c8b8p0, + } +}; diff --git a/sysdeps/ieee754/flt-32/e_logf.c b/sysdeps/ieee754/flt-32/e_logf.c index cf75e11781..879a673e12 100644 --- a/sysdeps/ieee754/flt-32/e_logf.c +++ b/sysdeps/ieee754/flt-32/e_logf.c @@ -1,85 +1,94 @@ -/* e_logf.c -- float version of e_log.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +/* Single-precision log function. + Copyright (C) 2017-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * ==================================================== - * 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. - * ==================================================== - */ + 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 <stdint.h> +#include <shlib-compat.h> +#include <libm-alias-float.h> +#include "math_config.h" -static const float -ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ -ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ -two25 = 3.355443200e+07, /* 0x4c000000 */ -Lg1 = 6.6666668653e-01, /* 3F2AAAAB */ -Lg2 = 4.0000000596e-01, /* 3ECCCCCD */ -Lg3 = 2.8571429849e-01, /* 3E924925 */ -Lg4 = 2.2222198546e-01, /* 3E638E29 */ -Lg5 = 1.8183572590e-01, /* 3E3A3325 */ -Lg6 = 1.5313838422e-01, /* 3E1CD04F */ -Lg7 = 1.4798198640e-01; /* 3E178897 */ +/* +LOGF_TABLE_BITS = 4 +LOGF_POLY_ORDER = 4 -static const float zero = 0.0; +ULP error: 0.818 (nearest rounding.) +Relative error: 1.957 * 2^-26 (before rounding.) +*/ + +#define T __logf_data.tab +#define A __logf_data.poly +#define Ln2 __logf_data.ln2 +#define N (1 << LOGF_TABLE_BITS) +#define OFF 0x3f330000 float -__ieee754_logf(float x) +__logf (float x) { - float hfsq,f,s,z,R,w,t1,t2,dk; - int32_t k,ix,i,j; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, y, y0, invc, logc; + uint32_t ix, iz, tmp; + int k, i; + + ix = asuint (x); +#if WANT_ROUNDING + /* Fix sign of zero with downward rounding when x==1. */ + if (__glibc_unlikely (ix == 0x3f800000)) + return 0; +#endif + if (__glibc_unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000)) + { + /* x < 0x1p-126 or inf or nan. */ + if (ix * 2 == 0) + return __math_divzerof (1); + if (ix == 0x7f800000) /* log(inf) == inf. */ + return x; + if ((ix & 0x80000000) || ix * 2 >= 0xff000000) + return __math_invalidf (x); + /* x is subnormal, normalize it. */ + ix = asuint (x * 0x1p23f); + ix -= 23 << 23; + } + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; + k = (int32_t) tmp >> 23; /* arithmetic shift */ + iz = ix - (tmp & 0x1ff << 23); + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); - GET_FLOAT_WORD(ix,x); + /* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */ + r = z * invc - 1; + y0 = logc + (double_t) k * Ln2; - k=0; - if (ix < 0x00800000) { /* x < 2**-126 */ - if (__builtin_expect((ix&0x7fffffff)==0, 0)) - return -two25/zero; /* log(+-0)=-inf */ - if (__builtin_expect(ix<0, 0)) - return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ - GET_FLOAT_WORD(ix,x); - } - if (__builtin_expect(ix >= 0x7f800000, 0)) return x+x; - k += (ix>>23)-127; - ix &= 0x007fffff; - i = (ix+(0x95f64<<3))&0x800000; - SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */ - k += (i>>23); - f = x-(float)1.0; - if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */ - if(f==zero) { - if(k==0) return zero; else {dk=(float)k; - return dk*ln2_hi+dk*ln2_lo;} - } - R = f*f*((float)0.5-(float)0.33333333333333333*f); - if(k==0) return f-R; else {dk=(float)k; - return dk*ln2_hi-((R-dk*ln2_lo)-f);} - } - s = f/((float)2.0+f); - dk = (float)k; - z = s*s; - i = ix-(0x6147a<<3); - w = z*z; - j = (0x6b851<<3)-ix; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - i |= j; - R = t2+t1; - if(i>0) { - hfsq=(float)0.5*f*f; - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); - } else { - if(k==0) return f-s*(f-R); else - return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); - } + /* Pipelined polynomial evaluation to approximate log1p(r). */ + r2 = r * r; + y = A[1] * r + A[2]; + y = A[0] * r2 + y; + y = y * r2 + (y0 + r); + return (float) y; } -strong_alias (__ieee754_logf, __logf_finite) +#ifndef __logf +strong_alias (__logf, __ieee754_logf) +strong_alias (__logf, __logf_finite) +versioned_symbol (libm, __logf, logf, GLIBC_2_27); +libm_alias_float_other (__log, log) +#endif diff --git a/sysdeps/ieee754/flt-32/e_logf_data.c b/sysdeps/ieee754/flt-32/e_logf_data.c new file mode 100644 index 0000000000..0c71414a8a --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_logf_data.c @@ -0,0 +1,44 @@ +/* Data definition for logf. + Copyright (C) 2017-2018 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_config.h" + +const struct logf_data __logf_data = { + .tab = { + { 0x1.661ec79f8f3bep+0, -0x1.57bf7808caadep-2 }, + { 0x1.571ed4aaf883dp+0, -0x1.2bef0a7c06ddbp-2 }, + { 0x1.49539f0f010bp+0, -0x1.01eae7f513a67p-2 }, + { 0x1.3c995b0b80385p+0, -0x1.b31d8a68224e9p-3 }, + { 0x1.30d190c8864a5p+0, -0x1.6574f0ac07758p-3 }, + { 0x1.25e227b0b8eap+0, -0x1.1aa2bc79c81p-3 }, + { 0x1.1bb4a4a1a343fp+0, -0x1.a4e76ce8c0e5ep-4 }, + { 0x1.12358f08ae5bap+0, -0x1.1973c5a611cccp-4 }, + { 0x1.0953f419900a7p+0, -0x1.252f438e10c1ep-5 }, + { 0x1p+0, 0x0p+0 }, + { 0x1.e608cfd9a47acp-1, 0x1.aa5aa5df25984p-5 }, + { 0x1.ca4b31f026aap-1, 0x1.c5e53aa362eb4p-4 }, + { 0x1.b2036576afce6p-1, 0x1.526e57720db08p-3 }, + { 0x1.9c2d163a1aa2dp-1, 0x1.bc2860d22477p-3 }, + { 0x1.886e6037841edp-1, 0x1.1058bc8a07ee1p-2 }, + { 0x1.767dcf5534862p-1, 0x1.4043057b6ee09p-2 }, + }, + .ln2 = 0x1.62e42fefa39efp-1, + .poly = { + -0x1.00ea348b88334p-2, 0x1.5575b0be00b6ap-2, -0x1.ffffef20a4123p-2, + } +}; diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c index c72fe37d3b..ece83f0dd2 100644 --- a/sysdeps/ieee754/flt-32/e_powf.c +++ b/sysdeps/ieee754/flt-32/e_powf.c @@ -1,258 +1,224 @@ -/* e_powf.c -- float version of e_pow.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +/* Single-precision pow function. + Copyright (C) 2017-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * ==================================================== - * 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. - * ==================================================== - */ + 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. -#include <math.h> -#include <math_private.h> + 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. -static const float huge = 1.0e+30, tiny = 1.0e-30; + 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/>. */ -static const float -bp[] = {1.0, 1.5,}, -dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ -dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ -zero = 0.0, -one = 1.0, -two = 2.0, -two24 = 16777216.0, /* 0x4b800000 */ - /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ -L1 = 6.0000002384e-01, /* 0x3f19999a */ -L2 = 4.2857143283e-01, /* 0x3edb6db7 */ -L3 = 3.3333334327e-01, /* 0x3eaaaaab */ -L4 = 2.7272811532e-01, /* 0x3e8ba305 */ -L5 = 2.3066075146e-01, /* 0x3e6c3255 */ -L6 = 2.0697501302e-01, /* 0x3e53f142 */ -P1 = 1.6666667163e-01, /* 0x3e2aaaab */ -P2 = -2.7777778450e-03, /* 0xbb360b61 */ -P3 = 6.6137559770e-05, /* 0x388ab355 */ -P4 = -1.6533901999e-06, /* 0xb5ddea0e */ -P5 = 4.1381369442e-08, /* 0x3331bb4c */ -lg2 = 6.9314718246e-01, /* 0x3f317218 */ -lg2_h = 6.93145752e-01, /* 0x3f317200 */ -lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ -ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ -cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ -cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */ -cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */ -ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */ -ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ -ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ +#include <math.h> +#include <stdint.h> +#include <shlib-compat.h> +#include <libm-alias-float.h> +#include "math_config.h" -float -__ieee754_powf(float x, float y) +/* +POWF_LOG2_POLY_ORDER = 5 +EXP2F_TABLE_BITS = 5 + +ULP error: 0.82 (~ 0.5 + relerr*2^24) +relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2) +relerr_log2: 1.83 * 2^-33 (Relative error of logx.) +relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).) +*/ + +#define N (1 << POWF_LOG2_TABLE_BITS) +#define T __powf_log2_data.tab +#define A __powf_log2_data.poly +#define OFF 0x3f330000 + +/* Subnormal input is normalized so ix has negative biased exponent. + Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set. */ +static inline double_t +log2_inline (uint32_t ix) { - float z,ax,z_h,z_l,p_h,p_l; - float y1,t1,t2,r,s,t,u,v,w; - int32_t i,j,k,yisint,n; - int32_t hx,hy,ix,iy,is; - - GET_FLOAT_WORD(hx,x); - GET_FLOAT_WORD(hy,y); - ix = hx&0x7fffffff; iy = hy&0x7fffffff; - - /* y==zero: x**0 = 1 */ - if(iy==0) return one; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, r4, p, q, y, y0, invc, logc; + uint32_t iz, top, tmp; + int k, i; + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N; + top = tmp & 0xff800000; + iz = ix - top; + k = (int32_t) top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */ + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); + + /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ + r = z * invc - 1; + y0 = logc + (double_t) k; + + /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ + r2 = r * r; + y = A[0] * r + A[1]; + p = A[2] * r + A[3]; + r4 = r2 * r2; + q = A[4] * r + y0; + q = p * r2 + q; + y = y * r4 + q; + return y; +} - /* x==+-1 */ - if(x == 1.0) return one; - if(x == -1.0 && isinf(y)) return one; +#undef N +#undef T +#define N (1 << EXP2F_TABLE_BITS) +#define T __exp2f_data.tab +#define SIGN_BIAS (1 << (EXP2F_TABLE_BITS + 11)) + +/* The output of log2 and thus the input of exp2 is either scaled by N + (in case of fast toint intrinsics) or not. The unscaled xd must be + in [-1021,1023], sign_bias sets the sign of the result. */ +static inline double_t +exp2_inline (double_t xd, uint32_t sign_bias) +{ + uint64_t ki, ski, t; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t kd, z, r, r2, y, s; + +#if TOINT_INTRINSICS +# define C __exp2f_data.poly_scaled + /* N*x = k + r with r in [-1/2, 1/2] */ + kd = roundtoint (xd); /* k */ + ki = converttoint (xd); +#else +# define C __exp2f_data.poly +# define SHIFT __exp2f_data.shift_scaled + /* x = k/N + r with r in [-1/(2N), 1/(2N)] */ + kd = (double) (xd + SHIFT); /* Rounding to double precision is required. */ + ki = asuint64 (kd); + kd -= SHIFT; /* k/N */ +#endif + r = xd - kd; + + /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + ski = ki + sign_bias; + t += ski << (52 - EXP2F_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return y; +} - /* +-NaN return x+y */ - if(__builtin_expect(ix > 0x7f800000 || - iy > 0x7f800000, 0)) - return x+y; +/* Returns 0 if not int, 1 if odd int, 2 if even int. */ +static inline int +checkint (uint32_t iy) +{ + int e = iy >> 23 & 0xff; + if (e < 0x7f) + return 0; + if (e > 0x7f + 23) + return 2; + if (iy & ((1 << (0x7f + 23 - e)) - 1)) + return 0; + if (iy & (1 << (0x7f + 23 - e))) + return 1; + return 2; +} - /* determine if y is an odd int when x < 0 - * yisint = 0 ... y is not an integer - * yisint = 1 ... y is an odd int - * yisint = 2 ... y is an even int - */ - yisint = 0; - if(hx<0) { - if(iy>=0x4b800000) yisint = 2; /* even integer y */ - else if(iy>=0x3f800000) { - k = (iy>>23)-0x7f; /* exponent */ - j = iy>>(23-k); - if((j<<(23-k))==iy) yisint = 2-(j&1); - } - } +static inline int +zeroinfnan (uint32_t ix) +{ + return 2 * ix - 1 >= 2u * 0x7f800000 - 1; +} - /* special value of y */ - if (__builtin_expect(iy==0x7f800000, 0)) { /* y is +-inf */ - if (ix==0x3f800000) - return y - y; /* inf**+-1 is NaN */ - else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ - return (hy>=0)? y: zero; - else /* (|x|<1)**-,+inf = inf,0 */ - return (hy<0)?-y: zero; - } - if(iy==0x3f800000) { /* y is +-1 */ - if(hy<0) return one/x; else return x; - } - if(hy==0x40000000) return x*x; /* y is 2 */ - if(hy==0x3f000000) { /* y is 0.5 */ - if(__builtin_expect(hx>=0, 1)) /* x >= +0 */ - return __ieee754_sqrtf(x); +float +__powf (float x, float y) +{ + uint32_t sign_bias = 0; + uint32_t ix, iy; + + ix = asuint (x); + iy = asuint (y); + if (__glibc_unlikely (ix - 0x00800000 >= 0x7f800000 - 0x00800000 + || zeroinfnan (iy))) + { + /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + if (2 * iy == 0) + return issignalingf_inline (x) ? x + y : 1.0f; + if (ix == 0x3f800000) + return issignalingf_inline (y) ? x + y : 1.0f; + if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) + return x + y; + if (2 * ix == 2 * 0x3f800000) + return 1.0f; + if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) + return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ + return y * y; } - - ax = fabsf(x); - /* special value of x */ - if(__builtin_expect(ix==0x7f800000||ix==0||ix==0x3f800000, 0)){ - z = ax; /*x is +-0,+-inf,+-1*/ - if(hy<0) z = one/z; /* z = (1/|x|) */ - if(hx<0) { - if(((ix-0x3f800000)|yisint)==0) { - z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) - z = -z; /* (x<0)**odd = -(|x|**odd) */ + if (__glibc_unlikely (zeroinfnan (ix))) + { + float_t x2 = x * x; + if (ix & 0x80000000 && checkint (iy) == 1) + { + x2 = -x2; + sign_bias = 1; } - return z; - } - - /* (x<0)**(non-int) is NaN */ - if(__builtin_expect(((((u_int32_t)hx>>31)-1)|yisint)==0, 0)) - return (x-x)/(x-x); - - /* |y| is huge */ - if(__builtin_expect(iy>0x4d000000, 0)) { /* if |y| > 2**27 */ - /* over/underflow if x is not close to one */ - if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; - 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 = 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; - t1 = u+v; - GET_FLOAT_WORD(is,t1); - SET_FLOAT_WORD(t1,is&0xfffff000); - t2 = v-(t1-u); - } else { - float s2,s_h,s_l,t_h,t_l; - /* Avoid internal underflow for tiny y. The exact value - of y does not matter if |y| <= 2**-32. */ - if (iy < 0x2f800000) - SET_FLOAT_WORD (y, (hy & 0x80000000) | 0x2f800000); - n = 0; - /* take care subnormal number */ - if(ix<0x00800000) - {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } - n += ((ix)>>23)-0x7f; - j = ix&0x007fffff; - /* determine interval */ - ix = j|0x3f800000; /* normalize ix */ - if(j<=0x1cc471) k=0; /* |x|<sqrt(3/2) */ - else if(j<0x5db3d7) k=1; /* |x|<sqrt(3) */ - else {k=0;n+=1;ix -= 0x00800000;} - SET_FLOAT_WORD(ax,ix); - - /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ - u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ - v = one/(ax+bp[k]); - s = u*v; - s_h = s; - 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)+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) */ - s2 = s*s; - r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); - r += s_l*(s_h+s); - s2 = s_h*s_h; - t_h = (float)3.0+s2+r; - GET_FLOAT_WORD(is,t_h); - SET_FLOAT_WORD(t_h,is&0xfffff000); - t_l = r-((t_h-(float)3.0)-s2); - /* u+v = s*(1+...) */ - u = s_h*t_h; - v = s_l*t_h+t_l*s; - /* 2/(3log2)*(s+...) */ - p_h = u+v; - GET_FLOAT_WORD(is,p_h); - SET_FLOAT_WORD(p_h,is&0xfffff000); - p_l = v-(p_h-u); - z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ - z_l = cp_l*p_h+p_l*cp+dp_l[k]; - /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ - t = (float)n; - t1 = (((z_h+z_l)+dp_h[k])+t); - GET_FLOAT_WORD(is,t1); - SET_FLOAT_WORD(t1,is&0xfffff000); - t2 = z_l-(((t1-t)-dp_h[k])-z_h); - } - - s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ - if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0) - s = -one; /* (-ve)**(odd int) */ - - /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ - GET_FLOAT_WORD(is,y); - SET_FLOAT_WORD(y1,is&0xfffff000); - p_l = (y-y1)*t1+y*t2; - p_h = y1*t1; - z = p_l+p_h; - GET_FLOAT_WORD(j,z); - if (__builtin_expect(j>0x43000000, 0)) /* if z > 128 */ - return s*huge*huge; /* overflow */ - else if (__builtin_expect(j==0x43000000, 0)) { /* if z == 128 */ - if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ +#if WANT_ERRNO + if (2 * ix == 0 && iy & 0x80000000) + return __math_divzerof (sign_bias); +#endif + return iy & 0x80000000 ? 1 / x2 : x2; } - else if (__builtin_expect((j&0x7fffffff)>0x43160000, 0))/* z <= -150 */ - return s*tiny*tiny; /* underflow */ - else if (__builtin_expect((u_int32_t) j==0xc3160000, 0)){/* z == -150*/ - if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ + /* x and y are non-zero finite. */ + if (ix & 0x80000000) + { + /* Finite x < 0. */ + int yint = checkint (iy); + if (yint == 0) + return __math_invalidf (x); + if (yint == 1) + sign_bias = SIGN_BIAS; + ix &= 0x7fffffff; } - /* - * compute 2**(p_h+p_l) - */ - i = j&0x7fffffff; - k = (i>>23)-0x7f; - n = 0; - if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ - n = j+(0x00800000>>(k+1)); - k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */ - SET_FLOAT_WORD(t,n&~(0x007fffff>>k)); - n = ((n&0x007fffff)|0x00800000)>>(23-k); - if(j<0) n = -n; - p_h -= t; + if (ix < 0x00800000) + { + /* Normalize subnormal x so exponent becomes negative. */ + ix = asuint (x * 0x1p23f); + ix &= 0x7fffffff; + ix -= 23 << 23; } - t = p_l+p_h; - GET_FLOAT_WORD(is,t); - SET_FLOAT_WORD(t,is&0xfffff000); - u = t*lg2_h; - v = (p_l-(t-p_h))*lg2+t*lg2_l; - z = u+v; - w = v-(z-u); - t = z*z; - t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - r = (z*t1)/(t1-two)-(w+z*w); - z = one-(r-z); - GET_FLOAT_WORD(j,z); - j += (n<<23); - 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; + } + double_t logx = log2_inline (ix); + double_t ylogx = y * logx; /* Note: cannot overflow, y is single prec. */ + if (__glibc_unlikely ((asuint64 (ylogx) >> 47 & 0xffff) + >= asuint64 (126.0 * POWF_SCALE) >> 47)) + { + /* |y*log(x)| >= 126. */ + if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE) + return __math_oflowf (sign_bias); + if (ylogx <= -150.0 * POWF_SCALE) + return __math_uflowf (sign_bias); +#if WANT_ERRNO_UFLOW + if (ylogx < -149.0 * POWF_SCALE) + return __math_may_uflowf (sign_bias); +#endif + } + return (float) exp2_inline (ylogx, sign_bias); } -strong_alias (__ieee754_powf, __powf_finite) +#ifndef __powf +strong_alias (__powf, __ieee754_powf) +strong_alias (__powf, __powf_finite) +versioned_symbol (libm, __powf, powf, GLIBC_2_27); +libm_alias_float_other (__pow, pow) +#endif diff --git a/sysdeps/ieee754/flt-32/e_powf_log2_data.c b/sysdeps/ieee754/flt-32/e_powf_log2_data.c new file mode 100644 index 0000000000..27cd6b7a81 --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_powf_log2_data.c @@ -0,0 +1,45 @@ +/* Data definition for powf. + Copyright (C) 2017-2018 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_config.h" + +const struct powf_log2_data __powf_log2_data = { + .tab = { + { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 * POWF_SCALE }, + { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 * POWF_SCALE }, + { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 * POWF_SCALE }, + { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 * POWF_SCALE }, + { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 * POWF_SCALE }, + { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 * POWF_SCALE }, + { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 * POWF_SCALE }, + { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 * POWF_SCALE }, + { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 * POWF_SCALE }, + { 0x1p+0, 0x0p+0 * POWF_SCALE }, + { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 * POWF_SCALE }, + { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 * POWF_SCALE }, + { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 * POWF_SCALE }, + { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 * POWF_SCALE }, + { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 * POWF_SCALE }, + { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 * POWF_SCALE }, + }, + .poly = { + 0x1.27616c9496e0bp-2 * POWF_SCALE, -0x1.71969a075c67ap-2 * POWF_SCALE, + 0x1.ec70a6ca7baddp-2 * POWF_SCALE, -0x1.7154748bef6c8p-1 * POWF_SCALE, + 0x1.71547652ab82bp0 * POWF_SCALE, + } +}; diff --git a/sysdeps/ieee754/flt-32/e_rem_pio2f.c b/sysdeps/ieee754/flt-32/e_rem_pio2f.c index 0928373498..bd871a26c2 100644 --- a/sysdeps/ieee754/flt-32/e_rem_pio2f.c +++ b/sysdeps/ieee754/flt-32/e_rem_pio2f.c @@ -100,7 +100,7 @@ int32_t __ieee754_rem_pio2f(float x, float *y) if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */ if(hx>0) { z = x - pio2_1; - if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ + if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */ y[0] = z - pio2_1t; y[1] = (z-y[0])-pio2_1t; } else { /* near pi/2, use 24+24+24 bit pi */ @@ -111,7 +111,7 @@ int32_t __ieee754_rem_pio2f(float x, float *y) return 1; } else { /* negative x */ z = x + pio2_1; - if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ + if((ix&0xffffffc0)!=0x3fc90fc0) { /* 24+24 bit pi OK */ y[0] = z + pio2_1t; y[1] = (z-y[0])+pio2_1t; } else { /* near pi/2, use 24+24+24 bit pi */ @@ -131,7 +131,7 @@ int32_t __ieee754_rem_pio2f(float x, float *y) if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) { y[0] = r-w; /* quick check no cancellation */ } else { - u_int32_t high; + uint32_t high; j = ix>>23; y[0] = r-w; GET_FLOAT_WORD(high,y[0]); diff --git a/sysdeps/ieee754/flt-32/e_remainderf.c b/sysdeps/ieee754/flt-32/e_remainderf.c index cc0167862e..8e78784e0f 100644 --- a/sysdeps/ieee754/flt-32/e_remainderf.c +++ b/sysdeps/ieee754/flt-32/e_remainderf.c @@ -23,7 +23,7 @@ float __ieee754_remainderf(float x, float p) { int32_t hx,hp; - u_int32_t sx; + uint32_t sx; float p_half; GET_FLOAT_WORD(hx,x); diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c index 6100d95c55..20f7db81ea 100644 --- a/sysdeps/ieee754/flt-32/e_sinhf.c +++ b/sysdeps/ieee754/flt-32/e_sinhf.c @@ -15,7 +15,9 @@ #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static const float one = 1.0, shuge = 1.0e37; diff --git a/sysdeps/ieee754/flt-32/e_sqrtf.c b/sysdeps/ieee754/flt-32/e_sqrtf.c index c02206ac01..6025da19cf 100644 --- a/sysdeps/ieee754/flt-32/e_sqrtf.c +++ b/sysdeps/ieee754/flt-32/e_sqrtf.c @@ -24,7 +24,7 @@ __ieee754_sqrtf(float x) float z; int32_t sign = (int)0x80000000; int32_t ix,s,q,m,t,i; - u_int32_t r; + uint32_t r; GET_FLOAT_WORD(ix,x); diff --git a/sysdeps/ieee754/flt-32/k_rem_pio2f.c b/sysdeps/ieee754/flt-32/k_rem_pio2f.c index 392afdbb6c..ea4915b765 100644 --- a/sysdeps/ieee754/flt-32/k_rem_pio2f.c +++ b/sysdeps/ieee754/flt-32/k_rem_pio2f.c @@ -18,7 +18,9 @@ static char rcsid[] = "$NetBSD: k_rem_pio2f.c,v 1.4 1995/05/10 20:46:28 jtc Exp #endif #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libc-diag.h> /* In the float version, the input parameter x contains 8 bit integers, not 24 bit integers. 113 bit precision is not supported. */ @@ -122,7 +124,16 @@ recompute: j = 0; for (i=jz-1;i>=jk;i--) j |= iq[i]; if(j==0) { /* need recomputation */ + /* On s390x gcc 6.1 -O3 produces the warning "array subscript is + below array bounds [-Werror=array-bounds]". Only + __ieee754_rem_pio2f calls __kernel_rem_pio2f for normal + numbers and |x| ~> 2^7*(pi/2). Thus x can't be zero and + ipio2 is not zero, too. Thus not all iq[] values can't be + zero. */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_NEEDS_COMMENT (6.1, "-Warray-bounds"); for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ + DIAG_POP_NEEDS_COMMENT; for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ f[jx+i] = (float) ipio2[jv+i]; @@ -172,7 +183,17 @@ recompute: float fv = 0.0; for (i=jz;i>=0;i--) fv = math_narrow_eval (fv + fq[i]); y[0] = (ih==0)? fv: -fv; + /* GCC mainline (to be GCC 9), as of 2018-05-22 on + i686, warns that fq[0] may be used uninitialized. + This is not possible because jz is always + nonnegative when the above loop initializing fq is + executed, because the result is never zero to full + precision (this function is not called for zero + arguments). */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_NEEDS_COMMENT (9, "-Wmaybe-uninitialized"); fv = math_narrow_eval (fq[0]-fv); + DIAG_POP_NEEDS_COMMENT; for (i=1;i<=jz;i++) fv = math_narrow_eval (fv + fq[i]); y[1] = (ih==0)? fv: -fv; break; diff --git a/sysdeps/ieee754/flt-32/k_sinf.c b/sysdeps/ieee754/flt-32/k_sinf.c index a195d59466..dcf3c35358 100644 --- a/sysdeps/ieee754/flt-32/k_sinf.c +++ b/sysdeps/ieee754/flt-32/k_sinf.c @@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const float half = 5.0000000000e-01,/* 0x3f000000 */ diff --git a/sysdeps/ieee754/flt-32/k_tanf.c b/sysdeps/ieee754/flt-32/k_tanf.c index 9f0e55860f..228ece2075 100644 --- a/sysdeps/ieee754/flt-32/k_tanf.c +++ b/sysdeps/ieee754/flt-32/k_tanf.c @@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: k_tanf.c,v 1.4 1995/05/10 20:46:39 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const float one = 1.0000000000e+00, /* 0x3f800000 */ pio4 = 7.8539812565e-01, /* 0x3f490fda */ diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c index f15719b059..01edb0b8de 100644 --- a/sysdeps/ieee754/flt-32/lgamma_negf.c +++ b/sysdeps/ieee754/flt-32/lgamma_negf.c @@ -1,5 +1,5 @@ /* lgammaf expanding around zeros. - Copyright (C) 2015-2016 Free Software Foundation, Inc. + Copyright (C) 2015-2018 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 @@ -18,6 +18,7 @@ #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> static const float lgamma_zeros[][2] = diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h new file mode 100644 index 0000000000..9c4ef30173 --- /dev/null +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -0,0 +1,164 @@ +/* Configuration for math routines. + Copyright (C) 2017-2018 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/>. */ + +#ifndef _MATH_CONFIG_H +#define _MATH_CONFIG_H + +#include <math.h> +#include <math_private.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +#ifndef WANT_ROUNDING +/* Correct special case results in non-nearest rounding modes. */ +# define WANT_ROUNDING 1 +#endif +#ifndef WANT_ERRNO +/* Set errno according to ISO C with (math_errhandling & MATH_ERRNO) != 0. */ +# define WANT_ERRNO 1 +#endif +#ifndef WANT_ERRNO_UFLOW +/* Set errno to ERANGE if result underflows to 0 (in all rounding modes). */ +# define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO) +#endif + +#ifndef TOINT_INTRINSICS +# define TOINT_INTRINSICS 0 +#endif +#ifndef TOINT_RINT +# define TOINT_RINT 0 +#endif +#ifndef TOINT_SHIFT +# define TOINT_SHIFT 1 +#endif + +static inline uint32_t +asuint (float f) +{ + union + { + float f; + uint32_t i; + } u = {f}; + return u.i; +} + +static inline float +asfloat (uint32_t i) +{ + union + { + uint32_t i; + float f; + } u = {i}; + return u.f; +} + +static inline uint64_t +asuint64 (double f) +{ + union + { + double f; + uint64_t i; + } u = {f}; + return u.i; +} + +static inline double +asdouble (uint64_t i) +{ + union + { + uint64_t i; + double f; + } u = {i}; + return u.f; +} + +static inline int +issignalingf_inline (float x) +{ + uint32_t ix = asuint (x); + if (HIGH_ORDER_BIT_IS_SET_FOR_SNAN) + return (ix & 0x7fc00000) == 0x7fc00000; + return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; +} + +#define NOINLINE __attribute__ ((noinline)) + +attribute_hidden float __math_oflowf (uint32_t); +attribute_hidden float __math_uflowf (uint32_t); +attribute_hidden float __math_may_uflowf (uint32_t); +attribute_hidden float __math_divzerof (uint32_t); +attribute_hidden float __math_invalidf (float); + +/* Shared between expf, exp2f and powf. */ +#define EXP2F_TABLE_BITS 5 +#define EXP2F_POLY_ORDER 3 +extern const struct exp2f_data +{ + uint64_t tab[1 << EXP2F_TABLE_BITS]; + double shift_scaled; + double poly[EXP2F_POLY_ORDER]; + double shift; + double invln2_scaled; + double poly_scaled[EXP2F_POLY_ORDER]; +} __exp2f_data attribute_hidden; + +#define LOGF_TABLE_BITS 4 +#define LOGF_POLY_ORDER 4 +extern const struct logf_data +{ + struct + { + double invc, logc; + } tab[1 << LOGF_TABLE_BITS]; + double ln2; + double poly[LOGF_POLY_ORDER - 1]; /* First order coefficient is 1. */ +} __logf_data attribute_hidden; + +#define LOG2F_TABLE_BITS 4 +#define LOG2F_POLY_ORDER 4 +extern const struct log2f_data +{ + struct + { + double invc, logc; + } tab[1 << LOG2F_TABLE_BITS]; + double poly[LOG2F_POLY_ORDER]; +} __log2f_data attribute_hidden; + +#define POWF_LOG2_TABLE_BITS 4 +#define POWF_LOG2_POLY_ORDER 5 +#if TOINT_INTRINSICS +# define POWF_SCALE_BITS EXP2F_TABLE_BITS +#else +# define POWF_SCALE_BITS 0 +#endif +#define POWF_SCALE ((double) (1 << POWF_SCALE_BITS)) +extern const struct powf_log2_data +{ + struct + { + double invc, logc; + } tab[1 << POWF_LOG2_TABLE_BITS]; + double poly[POWF_LOG2_POLY_ORDER]; +} __powf_log2_data attribute_hidden; + +#endif diff --git a/sysdeps/ieee754/flt-32/math_errf.c b/sysdeps/ieee754/flt-32/math_errf.c new file mode 100644 index 0000000000..5bc7ac6ef5 --- /dev/null +++ b/sysdeps/ieee754/flt-32/math_errf.c @@ -0,0 +1,76 @@ +/* Single-precision math error handling. + Copyright (C) 2017-2018 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_config.h" + +#if WANT_ERRNO +# include <errno.h> +/* NOINLINE reduces code size. */ +NOINLINE static float +with_errnof (float y, int e) +{ + errno = e; + return y; +} +#else +# define with_errnof(x, e) (x) +#endif + +/* NOINLINE prevents fenv semantics breaking optimizations. */ +NOINLINE static float +xflowf (uint32_t sign, float y) +{ + y = (sign ? -y : y) * y; + return with_errnof (y, ERANGE); +} + +attribute_hidden float +__math_uflowf (uint32_t sign) +{ + return xflowf (sign, 0x1p-95f); +} + +#if WANT_ERRNO_UFLOW +/* Underflows to zero in some non-nearest rounding mode, setting errno + is valid even if the result is non-zero, but in the subnormal range. */ +attribute_hidden float +__math_may_uflowf (uint32_t sign) +{ + return xflowf (sign, 0x1.4p-75f); +} +#endif + +attribute_hidden float +__math_oflowf (uint32_t sign) +{ + return xflowf (sign, 0x1p97f); +} + +attribute_hidden float +__math_divzerof (uint32_t sign) +{ + float y = 0; + return with_errnof ((sign ? -1 : 1) / y, ERANGE); +} + +attribute_hidden float +__math_invalidf (float x) +{ + float y = (x - x) / (x - x); + return isnan (x) ? y : with_errnof (y, EDOM); +} diff --git a/sysdeps/ieee754/flt-32/mpn2flt.c b/sysdeps/ieee754/flt-32/mpn2flt.c index 7f40f56d33..accca55fde 100644 --- a/sysdeps/ieee754/flt-32/mpn2flt.c +++ b/sysdeps/ieee754/flt-32/mpn2flt.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2016 Free Software Foundation, Inc. +/* Copyright (C) 1995-2018 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 da9cafb600..0812b54dca 100644 --- a/sysdeps/ieee754/flt-32/s_asinhf.c +++ b/sysdeps/ieee754/flt-32/s_asinhf.c @@ -16,6 +16,8 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-float.h> static const float one = 1.0000000000e+00, /* 0x3F800000 */ @@ -39,12 +41,12 @@ __asinhf(float x) } else { float xa = fabsf(x); if (ix>0x40000000) { /* 2**14 > |x| > 2.0 */ - w = __ieee754_logf(2.0f*xa+one/(__ieee754_sqrtf(xa*xa+one)+xa)); + w = __ieee754_logf(2.0f*xa+one/(sqrtf(xa*xa+one)+xa)); } else { /* 2.0 > |x| > 2**-14 */ float t = xa*xa; - w =__log1pf(xa+t/(one+__ieee754_sqrtf(one+t))); + w =__log1pf(xa+t/(one+sqrtf(one+t))); } } return __copysignf(w, x); } -weak_alias (__asinhf, asinhf) +libm_alias_float (__asinh, asinh) diff --git a/sysdeps/ieee754/flt-32/s_atanf.c b/sysdeps/ieee754/flt-32/s_atanf.c index e322a1d41f..f61b0f5ef6 100644 --- a/sysdeps/ieee754/flt-32/s_atanf.c +++ b/sysdeps/ieee754/flt-32/s_atanf.c @@ -20,6 +20,8 @@ static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-float.h> static const float atanhi[] = { 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */ @@ -98,4 +100,4 @@ float __atanf(float x) return (hx<0)? -z:z; } } -weak_alias (__atanf, atanf) +libm_alias_float (__atan, atan) diff --git a/sysdeps/ieee754/flt-32/s_cbrtf.c b/sysdeps/ieee754/flt-32/s_cbrtf.c index b9e2b3738c..afa1454399 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 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. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> #define CBRT2 1.2599210498948731648 /* 2^(1/3) */ @@ -60,4 +61,4 @@ __cbrtf (float x) return __ldexpf (x > 0.0 ? ym : -ym, xe / 3); } -weak_alias (__cbrtf, cbrtf) +libm_alias_float (__cbrt, cbrt) diff --git a/sysdeps/ieee754/flt-32/s_ceilf.c b/sysdeps/ieee754/flt-32/s_ceilf.c index 37659ea2ae..f289ec2341 100644 --- a/sysdeps/ieee754/flt-32/s_ceilf.c +++ b/sysdeps/ieee754/flt-32/s_ceilf.c @@ -15,27 +15,25 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> -static const float huge = 1.0e30; - float __ceilf(float x) { int32_t i0,j0; - u_int32_t i; + uint32_t i; GET_FLOAT_WORD(i0,x); j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { - if(j0<0) { /* raise inexact if x != 0 */ - math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ if(i0<0) {i0=0x80000000;} else if(i0!=0) { i0=0x3f800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - math_force_eval(huge+x); /* raise inexact flag */ if(i0>0) i0 += (0x00800000)>>j0; i0 &= (~i); } @@ -47,5 +45,5 @@ __ceilf(float x) return x; } #ifndef __ceilf -weak_alias (__ceilf, ceilf) +libm_alias_float (__ceil, ceil) #endif diff --git a/sysdeps/ieee754/flt-32/s_copysignf.c b/sysdeps/ieee754/flt-32/s_copysignf.c index 1621836065..3c4ac7ce68 100644 --- a/sysdeps/ieee754/flt-32/s_copysignf.c +++ b/sysdeps/ieee754/flt-32/s_copysignf.c @@ -25,13 +25,14 @@ static char rcsid[] = "$NetBSD: s_copysignf.c,v 1.4 1995/05/10 20:46:59 jtc Exp #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> float __copysignf(float x, float y) { - u_int32_t ix,iy; + uint32_t ix,iy; GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(iy,y); SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000)); return x; } -weak_alias (__copysignf, copysignf) +libm_alias_float (__copysign, copysign) diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c index 0affd406bb..061264d259 100644 --- a/sysdeps/ieee754/flt-32/s_cosf.c +++ b/sysdeps/ieee754/flt-32/s_cosf.c @@ -1,25 +1,26 @@ -/* s_cosf.c -- float version of s_cos.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +/* Compute cosine of argument. + Copyright (C) 2017-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * ==================================================== - * 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. - * ==================================================== - */ + 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. -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_cosf.c,v 1.4 1995/05/10 20:47:03 jtc Exp $"; -#endif + 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 <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> +#include "s_sincosf.h" #ifndef COSF # define COSF_FUNC __cosf @@ -27,37 +28,123 @@ static char rcsid[] = "$NetBSD: s_cosf.c,v 1.4 1995/05/10 20:47:03 jtc Exp $"; # define COSF_FUNC COSF #endif -float COSF_FUNC(float x) +float +COSF_FUNC (float x) { - float y[2],z=0.0; - int32_t n,ix; - - GET_FLOAT_WORD(ix,x); - - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if(ix <= 0x3f490fd8) return __kernel_cosf(x,z); - - /* cos(Inf or NaN) is NaN */ - else if (ix>=0x7f800000) { - if (ix == 0x7f800000) - __set_errno (EDOM); - return x-x; + double theta = x; + double abstheta = fabs (theta); + if (isless (abstheta, M_PI_4)) + { + double cx; + if (abstheta >= 0x1p-5) + { + const double theta2 = theta * theta; + /* Chebyshev polynomial of the form for cos: + * 1 + x^2 (C0 + x^2 (C1 + x^2 (C2 + x^2 (C3 + x^2 * C4)))). */ + cx = C3 + theta2 * C4; + cx = C2 + theta2 * cx; + cx = C1 + theta2 * cx; + cx = C0 + theta2 * cx; + cx = 1. + theta2 * cx; + return cx; } - - /* argument reduction needed */ - else { - n = __ieee754_rem_pio2f(x,y); - switch(n&3) { - case 0: return __kernel_cosf(y[0],y[1]); - case 1: return -__kernel_sinf(y[0],y[1],1); - case 2: return -__kernel_cosf(y[0],y[1]); - default: - return __kernel_sinf(y[0],y[1],1); + else if (abstheta >= 0x1p-27) + { + /* A simpler Chebyshev approximation is close enough for this range: + * 1 + x^2 (CC0 + x^3 * CC1). */ + const double theta2 = theta * theta; + cx = CC0 + theta * theta2 * CC1; + cx = 1.0 + theta2 * cx; + return cx; + } + else + { + /* For small enough |theta|, this is close enough. */ + return 1.0 - abstheta; + } + } + else /* |theta| >= Pi/4. */ + { + if (isless (abstheta, 9 * M_PI_4)) + { + /* There are cases where FE_UPWARD rounding mode can + produce a result of abstheta * inv_PI_4 == 9, + where abstheta < 9pi/4, so the domain for + pio2_table must go to 5 (9 / 2 + 1). */ + unsigned int n = (abstheta * inv_PI_4) + 1; + theta = abstheta - pio2_table[n / 2]; + return reduced_cos (theta, n); + } + else if (isless (abstheta, INFINITY)) + { + if (abstheta < 0x1p+23) + { + unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; + double x = n / 2; + theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; + /* Argument reduction needed. */ + return reduced_cos (theta, n); } + else /* |theta| >= 2^23. */ + { + x = fabsf (x); + int exponent; + GET_FLOAT_WORD (exponent, x); + exponent = (exponent >> FLOAT_EXPONENT_SHIFT) + - FLOAT_EXPONENT_BIAS; + exponent += 3; + exponent /= 28; + double a = invpio4_table[exponent] * x; + double b = invpio4_table[exponent + 1] * x; + double c = invpio4_table[exponent + 2] * x; + double d = invpio4_table[exponent + 3] * x; + uint64_t l = a; + l &= ~0x7; + a -= l; + double e = a + b; + l = e; + e = a - l; + if (l & 1) + { + e -= 1.0; + e += b; + e += c; + e += d; + e *= M_PI_4; + return reduced_cos (e, l + 1); + } + else + { + e += b; + e += c; + e += d; + if (e <= 1.0) + { + e *= M_PI_4; + return reduced_cos (e, l + 1); + } + else + { + l++; + e -= 2.0; + e *= M_PI_4; + return reduced_cos (e, l + 1); + } + } + } + } + else + { + int32_t ix; + GET_FLOAT_WORD (ix, abstheta); + /* cos(Inf or NaN) is NaN. */ + if (ix == 0x7f800000) /* Inf. */ + __set_errno (EDOM); + return x - x; } + } } #ifndef COSF -weak_alias (__cosf, cosf) +libm_alias_float (__cos, cos) #endif diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c index c8b6287503..693934d59a 100644 --- a/sysdeps/ieee754/flt-32/s_erff.c +++ b/sysdeps/ieee754/flt-32/s_erff.c @@ -20,7 +20,10 @@ static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $"; #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-float.h> #include <fix-int-fp-convert-zero.h> static const float @@ -104,7 +107,7 @@ float __erff(float x) GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; if(ix>=0x7f800000) { /* erf(nan)=nan */ - i = ((u_int32_t)hx>>31)<<1; + i = ((uint32_t)hx>>31)<<1; return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ } @@ -152,7 +155,7 @@ float __erff(float x) r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); if(hx>=0) return one-r/x; else return r/x-one; } -weak_alias (__erff, erff) +libm_alias_float (__erf, erf) float __erfcf(float x) { @@ -162,7 +165,7 @@ float __erfcf(float x) ix = hx&0x7fffffff; if(ix>=0x7f800000) { /* erfc(nan)=nan */ /* erfc(+-inf)=0,2 */ - float ret = (float)(((u_int32_t)hx>>31)<<1)+one/x; + float ret = (float)(((uint32_t)hx>>31)<<1)+one/x; if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0f) return 0.0f; return ret; @@ -227,4 +230,4 @@ float __erfcf(float x) return two-tiny; } } -weak_alias (__erfcf, erfcf) +libm_alias_float (__erfc, erfc) diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c index c515d25e28..b72cc0d083 100644 --- a/sysdeps/ieee754/flt-32/s_expm1f.c +++ b/sysdeps/ieee754/flt-32/s_expm1f.c @@ -16,7 +16,10 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-float.h> static const float huge = 1.0e+30; static const float tiny = 1.0e-30; @@ -39,7 +42,7 @@ __expm1f(float x) { float y,hi,lo,c,t,e,hxs,hfx,r1; int32_t k,xsb; - u_int32_t hx; + uint32_t hx; GET_FLOAT_WORD(hx,x); xsb = hx&0x80000000; /* sign bit of x */ @@ -127,4 +130,4 @@ __expm1f(float x) } return y; } -weak_alias (__expm1f, expm1f) +libm_alias_float (__expm1, expm1) diff --git a/sysdeps/ieee754/flt-32/s_fabsf.c b/sysdeps/ieee754/flt-32/s_fabsf.c index 297abe64bd..d3f6eb830b 100644 --- a/sysdeps/ieee754/flt-32/s_fabsf.c +++ b/sysdeps/ieee754/flt-32/s_fabsf.c @@ -22,9 +22,10 @@ static char rcsid[] = "$NetBSD: s_fabsf.c,v 1.4 1995/05/10 20:47:15 jtc Exp $"; */ #include <math.h> +#include <libm-alias-float.h> float __fabsf(float x) { return __builtin_fabsf (x); } -weak_alias (__fabsf, fabsf) +libm_alias_float (__fabs, fabs) diff --git a/sysdeps/ieee754/flt-32/s_finitef.c b/sysdeps/ieee754/flt-32/s_finitef.c index 4c5b339235..b9ae048559 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&0x7f800000)-0x7f800000)>>31); + return (int)((uint32_t)((ix&0x7f800000)-0x7f800000)>>31); } hidden_def (__finitef) weak_alias (__finitef, finitef) diff --git a/sysdeps/ieee754/flt-32/s_floorf.c b/sysdeps/ieee754/flt-32/s_floorf.c index 99d6c01833..12aed343a0 100644 --- a/sysdeps/ieee754/flt-32/s_floorf.c +++ b/sysdeps/ieee754/flt-32/s_floorf.c @@ -18,32 +18,28 @@ * Return x rounded toward -inf to integral value * Method: * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorf(x). */ #include <math.h> #include <math_private.h> - -static const float huge = 1.0e30; +#include <libm-alias-float.h> float __floorf(float x) { int32_t i0,j0; - u_int32_t i; + uint32_t i; GET_FLOAT_WORD(i0,x); j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { - if(j0<0) { /* raise inexact if x != 0 */ - math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ if(i0>=0) {i0=0;} else if((i0&0x7fffffff)!=0) { i0=0xbf800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - math_force_eval(huge+x); /* raise inexact flag */ if(i0<0) i0 += (0x00800000)>>j0; i0 &= (~i); } @@ -55,5 +51,5 @@ __floorf(float x) return x; } #ifndef __floorf -weak_alias (__floorf, floorf) +libm_alias_float (__floor, floor) #endif diff --git a/sysdeps/ieee754/flt-32/s_fpclassifyf.c b/sysdeps/ieee754/flt-32/s_fpclassifyf.c index 480557e2ae..01df9f26f5 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -25,7 +25,7 @@ int __fpclassifyf (float x) { - u_int32_t wx; + uint32_t wx; int retval = FP_NORMAL; GET_FLOAT_WORD (wx, x); diff --git a/sysdeps/ieee754/flt-32/s_frexpf.c b/sysdeps/ieee754/flt-32/s_frexpf.c index 67a28d3603..b7403bf0d8 100644 --- a/sysdeps/ieee754/flt-32/s_frexpf.c +++ b/sysdeps/ieee754/flt-32/s_frexpf.c @@ -19,6 +19,7 @@ static char rcsid[] = "$NetBSD: s_frexpf.c,v 1.5 1995/05/10 20:47:26 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> static const float two25 = 3.3554432000e+07; /* 0x4c000000 */ @@ -29,7 +30,7 @@ float __frexpf(float x, int *eptr) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; *eptr = 0; - if(ix>=0x7f800000||(ix==0)) return x; /* 0,inf,nan */ + if(ix>=0x7f800000||(ix==0)) return x + x; /* 0,inf,nan */ if (ix<0x00800000) { /* subnormal */ x *= two25; GET_FLOAT_WORD(hx,x); @@ -41,4 +42,4 @@ float __frexpf(float x, int *eptr) SET_FLOAT_WORD(x,hx); return x; } -weak_alias (__frexpf, frexpf) +libm_alias_float (__frexp, frexp) diff --git a/sysdeps/ieee754/flt-32/s_fromfpf.c b/sysdeps/ieee754/flt-32/s_fromfpf.c new file mode 100644 index 0000000000..d0c83b8f5d --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_fromfpf.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 0 +#define FUNC __fromfpf +#include <s_fromfpf_main.c> +libm_alias_float (__fromfp, fromfp) diff --git a/sysdeps/ieee754/flt-32/s_fromfpf_main.c b/sysdeps/ieee754/flt-32/s_fromfpf_main.c new file mode 100644 index 0000000000..b220b7212c --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_fromfpf_main.c @@ -0,0 +1,83 @@ +/* Round to integer type. flt-32 version. + Copyright (C) 2016-2018 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 <errno.h> +#include <fenv.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-float.h> +#include <stdbool.h> +#include <stdint.h> + +#define BIAS 0x7f +#define MANT_DIG 24 + +#if UNSIGNED +# define RET_TYPE uintmax_t +#else +# define RET_TYPE intmax_t +#endif + +#include <fromfp.h> + +RET_TYPE +FUNC (float x, int round, unsigned int width) +{ + if (width > INTMAX_WIDTH) + width = INTMAX_WIDTH; + uint32_t ix; + GET_FLOAT_WORD (ix, x); + bool negative = (ix & 0x80000000) != 0; + if (width == 0) + return fromfp_domain_error (negative, width); + ix &= 0x7fffffff; + if (ix == 0) + return 0; + int exponent = ix >> (MANT_DIG - 1); + exponent -= BIAS; + int max_exponent = fromfp_max_exponent (negative, width); + if (exponent > max_exponent) + return fromfp_domain_error (negative, width); + + ix &= ((1U << (MANT_DIG - 1)) - 1); + ix |= 1U << (MANT_DIG - 1); + uintmax_t uret; + bool half_bit, more_bits; + if (exponent >= MANT_DIG - 1) + { + uret = ix; + uret <<= exponent - (MANT_DIG - 1); + half_bit = false; + more_bits = false; + } + else if (exponent >= -1) + { + uint32_t h = 1U << (MANT_DIG - 2 - exponent); + half_bit = (ix & h) != 0; + more_bits = (ix & (h - 1)) != 0; + uret = ix >> (MANT_DIG - 1 - exponent); + } + else + { + uret = 0; + half_bit = false; + more_bits = true; + } + return fromfp_round_and_return (negative, uret, half_bit, more_bits, round, + exponent, max_exponent, width); +} diff --git a/sysdeps/ieee754/flt-32/s_fromfpxf.c b/sysdeps/ieee754/flt-32/s_fromfpxf.c new file mode 100644 index 0000000000..01d9247c31 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_fromfpxf.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 1 +#define FUNC __fromfpxf +#include <s_fromfpf_main.c> +libm_alias_float (__fromfpx, fromfpx) diff --git a/sysdeps/ieee754/flt-32/s_getpayloadf.c b/sysdeps/ieee754/flt-32/s_getpayloadf.c new file mode 100644 index 0000000000..98a16f9347 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_getpayloadf.c @@ -0,0 +1,35 @@ +/* Get NaN payload. flt-32 version. + Copyright (C) 2016-2018 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 <fix-int-fp-convert-zero.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-float.h> +#include <stdint.h> + +float +__getpayloadf (const float *x) +{ + uint32_t ix; + GET_FLOAT_WORD (ix, *x); + ix &= 0x3fffff; + if (FIX_INT_FP_CONVERT_ZERO && ix == 0) + return 0.0f; + return (float) ix; +} +libm_alias_float (__getpayload, getpayload) diff --git a/sysdeps/ieee754/flt-32/s_isnanf.c b/sysdeps/ieee754/flt-32/s_isnanf.c index 820b31a2b4..52129735fe 100644 --- a/sysdeps/ieee754/flt-32/s_isnanf.c +++ b/sysdeps/ieee754/flt-32/s_isnanf.c @@ -32,7 +32,7 @@ int __isnanf(float x) GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; ix = 0x7f800000 - ix; - return (int)(((u_int32_t)(ix))>>31); + return (int)(((uint32_t)(ix))>>31); } hidden_def (__isnanf) weak_alias (__isnanf, isnanf) diff --git a/sysdeps/ieee754/flt-32/s_issignalingf.c b/sysdeps/ieee754/flt-32/s_issignalingf.c index 2409ff408c..a18ba3e586 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-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 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 @@ -18,13 +18,14 @@ #include <math.h> #include <math_private.h> +#include <nan-high-order-bit.h> int __issignalingf (float x) { - u_int32_t xi; + uint32_t xi; GET_FLOAT_WORD (xi, x); -#ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN /* We only have to care about the high-order bit of x's significand, because having it set (sNaN) already makes the significand different from that used to designate infinity. */ diff --git a/sysdeps/ieee754/flt-32/s_llrintf.c b/sysdeps/ieee754/flt-32/s_llrintf.c index 56415d34f8..7c64bb2db4 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,7 +22,9 @@ #include <limits.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libm-alias-float.h> #include <fix-fp-int-convert-overflow.h> static const float two23[2] = @@ -36,7 +38,7 @@ long long int __llrintf (float x) { int32_t j0; - u_int32_t i0; + uint32_t i0; float w; float t; long long int result; @@ -83,4 +85,4 @@ __llrintf (float x) return sx ? -result : result; } -weak_alias (__llrintf, llrintf) +libm_alias_float (__llrint, llrint) diff --git a/sysdeps/ieee754/flt-32/s_llroundf.c b/sysdeps/ieee754/flt-32/s_llroundf.c index 1d35f71321..5457f9fa88 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> #include <fix-fp-int-convert-overflow.h> @@ -29,7 +30,7 @@ long long int __llroundf (float x) { int32_t j0; - u_int32_t i; + uint32_t i; long long int result; int sign; @@ -70,4 +71,4 @@ __llroundf (float x) return sign * result; } -weak_alias (__llroundf, llroundf) +libm_alias_float (__llround, llround) diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index ade60a2e27..a06190a589 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -15,7 +15,10 @@ #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> +#include <libc-diag.h> static const float ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ @@ -97,6 +100,18 @@ __log1pf(float x) s = f/((float)2.0+f); z = s*s; R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); + if (k == 0) + return f - (hfsq - s * (hfsq + R)); + else + { + /* With GCC 7 when compiling with -Os the compiler warns + that c might be used uninitialized. This can't be true + because k must be 0 for c to be uninitialized and we + handled that computation earlier without using c. */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_Os_NEEDS_COMMENT (7, "-Wmaybe-uninitialized"); + return k * ln2_hi - ((hfsq - (s * (hfsq + R) + + (k * ln2_lo + c))) - f); + DIAG_POP_NEEDS_COMMENT; + } } diff --git a/sysdeps/ieee754/flt-32/s_logbf.c b/sysdeps/ieee754/flt-32/s_logbf.c index 9ae20e332a..e0d4f3fbab 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 <libm-alias-float.h> #include <fix-int-fp-convert-zero.h> float @@ -38,4 +39,4 @@ __logbf (float x) return 0.0f; return (float) (rix - 127); } -weak_alias (__logbf, logbf) +libm_alias_float (__logb, logb) diff --git a/sysdeps/ieee754/flt-32/s_lrintf.c b/sysdeps/ieee754/flt-32/s_lrintf.c index 3b480a2107..5171377a4e 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,7 +22,9 @@ #include <limits.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libm-alias-float.h> #include <fix-fp-int-convert-overflow.h> static const float two23[2] = @@ -36,7 +38,7 @@ long int __lrintf (float x) { int32_t j0; - u_int32_t i0; + uint32_t i0; float w; float t; long int result; @@ -83,4 +85,4 @@ __lrintf (float x) return sx ? -result : result; } -weak_alias (__lrintf, lrintf) +libm_alias_float (__lrint, lrint) diff --git a/sysdeps/ieee754/flt-32/s_lroundf.c b/sysdeps/ieee754/flt-32/s_lroundf.c index 116c9e0627..20e7216640 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> #include <fix-fp-int-convert-overflow.h> @@ -29,7 +30,7 @@ long int __lroundf (float x) { int32_t j0; - u_int32_t i; + uint32_t i; long int result; int sign; @@ -70,4 +71,4 @@ __lroundf (float x) return sign * result; } -weak_alias (__lroundf, lroundf) +libm_alias_float (__lround, lround) diff --git a/sysdeps/ieee754/flt-32/s_modff.c b/sysdeps/ieee754/flt-32/s_modff.c index 23f6a902b3..45a94594fe 100644 --- a/sysdeps/ieee754/flt-32/s_modff.c +++ b/sysdeps/ieee754/flt-32/s_modff.c @@ -15,6 +15,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> static const float one = 1.0; @@ -22,7 +23,7 @@ float __modff(float x, float *iptr) { int32_t i0,j0; - u_int32_t i; + uint32_t i; GET_FLOAT_WORD(i0,x); j0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */ if(__builtin_expect(j0<23, 1)) { /* integer part in x */ @@ -32,7 +33,7 @@ __modff(float x, float *iptr) } else { i = (0x007fffff)>>j0; if((i0&i)==0) { /* x is integral */ - u_int32_t ix; + uint32_t ix; *iptr = x; GET_FLOAT_WORD(ix,x); SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ @@ -51,4 +52,4 @@ __modff(float x, float *iptr) return x; } } -weak_alias (__modff, modff) +libm_alias_float (__modf, modf) diff --git a/sysdeps/ieee754/flt-32/s_nearbyintf.c b/sysdeps/ieee754/flt-32/s_nearbyintf.c index 5aebefafcf..4dfe491f27 100644 --- a/sysdeps/ieee754/flt-32/s_nearbyintf.c +++ b/sysdeps/ieee754/flt-32/s_nearbyintf.c @@ -17,7 +17,9 @@ #include <fenv.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-float.h> static const float TWO23[2]={ @@ -37,7 +39,7 @@ __nearbyintf(float x) if(j0<23) { if(j0<0) { libc_feholdexceptf (&env); - w = TWO23[sx]+x; + w = TWO23[sx] + math_opt_barrier (x); t = w-TWO23[sx]; math_force_eval (t); libc_fesetenvf (&env); @@ -50,10 +52,10 @@ __nearbyintf(float x) else return x; /* x is integral */ } libc_feholdexceptf (&env); - w = TWO23[sx]+x; + w = TWO23[sx] + math_opt_barrier (x); t = w-TWO23[sx]; math_force_eval (t); libc_fesetenvf (&env); return t; } -weak_alias (__nearbyintf, nearbyintf) +libm_alias_float (__nearbyint, nearbyint) diff --git a/sysdeps/ieee754/flt-32/s_nextafterf.c b/sysdeps/ieee754/flt-32/s_nextafterf.c index 625d54b768..aa49df5a9e 100644 --- a/sysdeps/ieee754/flt-32/s_nextafterf.c +++ b/sysdeps/ieee754/flt-32/s_nextafterf.c @@ -19,7 +19,9 @@ static char rcsid[] = "$NetBSD: s_nextafterf.c,v 1.4 1995/05/10 20:48:01 jtc Exp #include <errno.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-float.h> #include <float.h> float __nextafterf(float x, float y) @@ -70,4 +72,4 @@ float __nextafterf(float x, float y) SET_FLOAT_WORD(x,hx); return x; } -weak_alias (__nextafterf, nextafterf) +libm_alias_float (__nextafter, nextafter) diff --git a/sysdeps/ieee754/flt-32/s_nextupf.c b/sysdeps/ieee754/flt-32/s_nextupf.c new file mode 100644 index 0000000000..87ec7ba21f --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_nextupf.c @@ -0,0 +1,48 @@ +/* Return the least floating-point number greater than X. + Copyright (C) 2016-2018 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> +#include <libm-alias-float.h> + +/* Return the least floating-point number greater than X. */ +float +__nextupf (float x) +{ + int32_t hx, ix; + + GET_FLOAT_WORD (hx, x); + ix = hx & 0x7fffffff; + if (ix == 0) + return FLT_TRUE_MIN; + if (ix > 0x7f800000) /* x is nan. */ + return x + x; + if (hx >= 0) + { /* x > 0. */ + if (isinf (x)) + return x; + hx += 1; + } + else + hx -= 1; + SET_FLOAT_WORD (x, hx); + return x; +} + +libm_alias_float (__nextup, nextup) diff --git a/sysdeps/ieee754/flt-32/s_remquof.c b/sysdeps/ieee754/flt-32/s_remquof.c index ecf831deaf..66704d46c2 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -23,13 +23,14 @@ static const float zero = 0.0; +#include <libm-alias-float.h> float __remquof (float x, float y, int *quo) { int32_t hx,hy; - u_int32_t sx; + uint32_t sx; int cquo, qs; GET_FLOAT_WORD (hx, x); @@ -107,4 +108,4 @@ __remquof (float x, float y, int *quo) x = -x; return x; } -weak_alias (__remquof, remquof) +libm_alias_float (__remquo, remquo) diff --git a/sysdeps/ieee754/flt-32/s_rintf.c b/sysdeps/ieee754/flt-32/s_rintf.c index 8a907488f7..db6f260a0b 100644 --- a/sysdeps/ieee754/flt-32/s_rintf.c +++ b/sysdeps/ieee754/flt-32/s_rintf.c @@ -15,6 +15,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> static const float TWO23[2]={ @@ -46,5 +47,5 @@ __rintf(float x) return w-TWO23[sx]; } #ifndef __rintf -weak_alias (__rintf, rintf) +libm_alias_float (__rint, rint) #endif diff --git a/sysdeps/ieee754/flt-32/s_roundevenf.c b/sysdeps/ieee754/flt-32/s_roundevenf.c new file mode 100644 index 0000000000..90f991d5c4 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_roundevenf.c @@ -0,0 +1,70 @@ +/* Round to nearest integer value, rounding halfway cases to even. + flt-32 version. + Copyright (C) 2016-2018 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 <libm-alias-float.h> +#include <stdint.h> + +#define BIAS 0x7f +#define MANT_DIG 24 +#define MAX_EXP (2 * BIAS + 1) + +float +__roundevenf (float x) +{ + uint32_t ix, ux; + GET_FLOAT_WORD (ix, x); + ux = ix & 0x7fffffff; + int exponent = ux >> (MANT_DIG - 1); + if (exponent >= BIAS + MANT_DIG - 1) + { + /* Integer, infinity or NaN. */ + if (exponent == MAX_EXP) + /* Infinity or NaN; quiet signaling NaNs. */ + return x + x; + else + return x; + } + else if (exponent >= BIAS) + { + /* At least 1; not necessarily an integer. Locate the bits with + exponents 0 and -1 (when the unbiased exponent is 0, the bit + with exponent 0 is implicit, but as the bias is odd it is OK + to take it from the low bit of the exponent). */ + int int_pos = (BIAS + MANT_DIG - 1) - exponent; + int half_pos = int_pos - 1; + uint32_t half_bit = 1U << half_pos; + uint32_t int_bit = 1U << int_pos; + if ((ix & (int_bit | (half_bit - 1))) != 0) + /* Carry into the exponent works correctly. No need to test + whether HALF_BIT is set. */ + ix += half_bit; + ix &= ~(int_bit - 1); + } + else if (exponent == BIAS - 1 && ux > 0x3f000000) + /* Interval (0.5, 1). */ + ix = (ix & 0x80000000) | 0x3f800000; + else + /* Rounds to 0. */ + ix &= 0x80000000; + SET_FLOAT_WORD (x, ix); + return x; +} +libm_alias_float (__roundeven, roundeven) diff --git a/sysdeps/ieee754/flt-32/s_roundf.c b/sysdeps/ieee754/flt-32/s_roundf.c index a75d98f384..7c95125d9c 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,9 +20,7 @@ #include <math.h> #include <math_private.h> - - -static const float huge = 1.0e30; +#include <libm-alias-float.h> float @@ -36,21 +34,17 @@ __roundf (float x) { if (j0 < 0) { - math_force_eval (huge + x); - i0 &= 0x80000000; if (j0 == -1) i0 |= 0x3f800000; } else { - u_int32_t i = 0x007fffff >> j0; + uint32_t i = 0x007fffff >> j0; if ((i0 & i) == 0) /* X is integral. */ return x; - math_force_eval (huge + x); - /* Raise inexact if x != 0. */ i0 += 0x00400000 >> j0; i0 &= ~i; } @@ -67,4 +61,4 @@ __roundf (float x) SET_FLOAT_WORD (x, i0); return x; } -weak_alias (__roundf, roundf) +libm_alias_float (__round, round) diff --git a/sysdeps/ieee754/flt-32/s_setpayloadf.c b/sysdeps/ieee754/flt-32/s_setpayloadf.c new file mode 100644 index 0000000000..6faf26bdcc --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_setpayloadf.c @@ -0,0 +1,4 @@ +#define SIG 0 +#define FUNC __setpayloadf +#include <s_setpayloadf_main.c> +libm_alias_float (__setpayload, setpayload) diff --git a/sysdeps/ieee754/flt-32/s_setpayloadf_main.c b/sysdeps/ieee754/flt-32/s_setpayloadf_main.c new file mode 100644 index 0000000000..a01c8677a3 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_setpayloadf_main.c @@ -0,0 +1,54 @@ +/* Set NaN payload. flt-32 version. + Copyright (C) 2016-2018 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 <libm-alias-float.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +#define SET_HIGH_BIT (HIGH_ORDER_BIT_IS_SET_FOR_SNAN ? SIG : !SIG) +#define BIAS 0x7f +#define PAYLOAD_DIG 22 +#define EXPLICIT_MANT_DIG 23 + +int +FUNC (float *x, float payload) +{ + uint32_t ix; + GET_FLOAT_WORD (ix, payload); + int exponent = ix >> EXPLICIT_MANT_DIG; + /* Test if argument is (a) negative or too large; (b) too small, + except for 0 when allowed; (c) not an integer. */ + if (exponent >= BIAS + PAYLOAD_DIG + || (exponent < BIAS && !(SET_HIGH_BIT && ix == 0)) + || (ix & ((1U << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0) + { + SET_FLOAT_WORD (*x, 0); + return 1; + } + if (ix != 0) + { + ix &= (1U << EXPLICIT_MANT_DIG) - 1; + ix |= 1U << EXPLICIT_MANT_DIG; + ix >>= BIAS + EXPLICIT_MANT_DIG - exponent; + } + ix |= 0x7f800000 | (SET_HIGH_BIT ? 0x400000 : 0); + SET_FLOAT_WORD (*x, ix); + return 0; +} diff --git a/sysdeps/ieee754/flt-32/s_setpayloadsigf.c b/sysdeps/ieee754/flt-32/s_setpayloadsigf.c new file mode 100644 index 0000000000..f08c877dc6 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_setpayloadsigf.c @@ -0,0 +1,4 @@ +#define SIG 1 +#define FUNC __setpayloadsigf +#include <s_setpayloadf_main.c> +libm_alias_float (__setpayloadsig, setpayloadsig) diff --git a/sysdeps/ieee754/flt-32/s_signbitf.c b/sysdeps/ieee754/flt-32/s_signbitf.c index 3cdde778d5..e97c6261d2 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 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_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index e0737b5ce4..d4a5a1b22c 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,7 +1,6 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 2017-2018 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 @@ -19,8 +18,9 @@ #include <errno.h> #include <math.h> - #include <math_private.h> +#include <libm-alias-float.h> +#include "s_sincosf.h" #ifndef SINCOSF # define SINCOSF_FUNC __sincosf @@ -31,54 +31,141 @@ void SINCOSF_FUNC (float x, float *sinx, float *cosx) { - int32_t ix; - - /* High word of x. */ - GET_FLOAT_WORD (ix, x); - - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if (ix <= 0x3f490fd8) - { - *sinx = __kernel_sinf (x, 0.0, 0); - *cosx = __kernel_cosf (x, 0.0); - } - else if (ix>=0x7f800000) + double cx; + double theta = x; + double abstheta = fabs (theta); + /* If |x|< Pi/4. */ + if (isless (abstheta, M_PI_4)) { - /* sin(Inf or NaN) is NaN */ - *sinx = *cosx = x - x; - if (ix == 0x7f800000) - __set_errno (EDOM); + if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ + { + const double theta2 = theta * theta; + /* Chebyshev polynomial of the form for sin and cos. */ + cx = C3 + theta2 * C4; + cx = C2 + theta2 * cx; + cx = C1 + theta2 * cx; + cx = C0 + theta2 * cx; + cx = 1.0 + theta2 * cx; + *cosx = cx; + cx = S3 + theta2 * S4; + cx = S2 + theta2 * cx; + cx = S1 + theta2 * cx; + cx = S0 + theta2 * cx; + cx = theta + theta * theta2 * cx; + *sinx = cx; + } + else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ + { + /* A simpler Chebyshev approximation is close enough for this range: + for sin: x+x^3*(SS0+x^2*SS1) + for cos: 1.0+x^2*(CC0+x^3*CC1). */ + const double theta2 = theta * theta; + cx = CC0 + theta * theta2 * CC1; + cx = 1.0 + theta2 * cx; + *cosx = cx; + cx = SS0 + theta2 * SS1; + cx = theta + theta * theta2 * cx; + *sinx = cx; + } + else + { + /* Handle some special cases. */ + if (theta) + *sinx = theta - (theta * SMALL); + else + *sinx = theta; + *cosx = 1.0 - abstheta; + } } - else + else /* |x| >= Pi/4. */ { - /* Argument reduction needed. */ - float y[2]; - int n; - - n = __ieee754_rem_pio2f (x, y); - switch (n & 3) + unsigned int signbit = isless (x, 0); + if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ + { + /* There are cases where FE_UPWARD rounding mode can + produce a result of abstheta * inv_PI_4 == 9, + where abstheta < 9pi/4, so the domain for + pio2_table must go to 5 (9 / 2 + 1). */ + unsigned int n = (abstheta * inv_PI_4) + 1; + theta = abstheta - pio2_table[n / 2]; + *sinx = reduced_sin (theta, n, signbit); + *cosx = reduced_cos (theta, n); + } + else if (isless (abstheta, INFINITY)) + { + if (abstheta < 0x1p+23) /* |x| < 2^23. */ + { + unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; + double x = n / 2; + theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; + /* Argument reduction needed. */ + *sinx = reduced_sin (theta, n, signbit); + *cosx = reduced_cos (theta, n); + } + else /* |x| >= 2^23. */ + { + x = fabsf (x); + int exponent; + GET_FLOAT_WORD (exponent, x); + exponent + = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; + exponent += 3; + exponent /= 28; + double a = invpio4_table[exponent] * x; + double b = invpio4_table[exponent + 1] * x; + double c = invpio4_table[exponent + 2] * x; + double d = invpio4_table[exponent + 3] * x; + uint64_t l = a; + l &= ~0x7; + a -= l; + double e = a + b; + l = e; + e = a - l; + if (l & 1) + { + e -= 1.0; + e += b; + e += c; + e += d; + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + else + { + e += b; + e += c; + e += d; + if (e <= 1.0) + { + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + else + { + l++; + e -= 2.0; + e *= M_PI_4; + *sinx = reduced_sin (e, l + 1, signbit); + *cosx = reduced_cos (e, l + 1); + } + } + } + } + else { - case 0: - *sinx = __kernel_sinf (y[0], y[1], 1); - *cosx = __kernel_cosf (y[0], y[1]); - break; - case 1: - *sinx = __kernel_cosf (y[0], y[1]); - *cosx = -__kernel_sinf (y[0], y[1], 1); - break; - case 2: - *sinx = -__kernel_sinf (y[0], y[1], 1); - *cosx = -__kernel_cosf (y[0], y[1]); - break; - default: - *sinx = -__kernel_cosf (y[0], y[1]); - *cosx = __kernel_sinf (y[0], y[1], 1); - break; + int32_t ix; + /* High word of x. */ + GET_FLOAT_WORD (ix, abstheta); + /* sin/cos(Inf or NaN) is NaN. */ + *sinx = *cosx = x - x; + if (ix == 0x7f800000) + __set_errno (EDOM); } } } #ifndef SINCOSF -weak_alias (__sincosf, sincosf) +libm_alias_float (__sincos, sincos) #endif diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h new file mode 100644 index 0000000000..35b5eee536 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincosf.h @@ -0,0 +1,155 @@ +/* Used by sinf, cosf and sincosf functions. + Copyright (C) 2017-2018 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/>. */ + +/* Chebyshev constants for cos, range -PI/4 - PI/4. */ +static const double C0 = -0x1.ffffffffe98aep-2; +static const double C1 = 0x1.55555545c50c7p-5; +static const double C2 = -0x1.6c16b348b6874p-10; +static const double C3 = 0x1.a00eb9ac43ccp-16; +static const double C4 = -0x1.23c97dd8844d7p-22; + +/* Chebyshev constants for sin, range -PI/4 - PI/4. */ +static const double S0 = -0x1.5555555551cd9p-3; +static const double S1 = 0x1.1111110c2688bp-7; +static const double S2 = -0x1.a019f8b4bd1f9p-13; +static const double S3 = 0x1.71d7264e6b5b4p-19; +static const double S4 = -0x1.a947e1674b58ap-26; + +/* Chebyshev constants for sin, range 2^-27 - 2^-5. */ +static const double SS0 = -0x1.555555543d49dp-3; +static const double SS1 = 0x1.110f475cec8c5p-7; + +/* Chebyshev constants for cos, range 2^-27 - 2^-5. */ +static const double CC0 = -0x1.fffffff5cc6fdp-2; +static const double CC1 = 0x1.55514b178dac5p-5; + +/* PI/2 with 98 bits of accuracy. */ +static const double PI_2_hi = 0x1.921fb544p+0; +static const double PI_2_lo = 0x1.0b4611a626332p-34; + +static const double SMALL = 0x1p-50; /* 2^-50. */ +static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ + +#define FLOAT_EXPONENT_SHIFT 23 +#define FLOAT_EXPONENT_BIAS 127 + +static const double pio2_table[] = { + 0 * M_PI_2, + 1 * M_PI_2, + 2 * M_PI_2, + 3 * M_PI_2, + 4 * M_PI_2, + 5 * M_PI_2 +}; + +static const double invpio4_table[] = { + 0x0p+0, + 0x1.45f306cp+0, + 0x1.c9c882ap-28, + 0x1.4fe13a8p-58, + 0x1.f47d4dp-85, + 0x1.bb81b6cp-112, + 0x1.4acc9ep-142, + 0x1.0e4107cp-169 +}; + +static const double ones[] = { 1.0, -1.0 }; + +/* Compute the sine value using Chebyshev polynomials where + THETA is the range reduced absolute value of the input + and it is less than Pi/4, + N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide + whether a sine or cosine approximation is more accurate and + SIGNBIT is used to add the correct sign after the Chebyshev + polynomial is computed. */ +static inline float +reduced_sin (const double theta, const unsigned int n, + const unsigned int signbit) +{ + double sx; + const double theta2 = theta * theta; + /* We are operating on |x|, so we need to add back the original + signbit for sinf. */ + double sign; + /* Determine positive or negative primary interval. */ + sign = ones[((n >> 2) & 1) ^ signbit]; + /* Are we in the primary interval of sin or cos? */ + if ((n & 2) == 0) + { + /* Here sinf() is calculated using sin Chebyshev polynomial: + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + sx = S3 + theta2 * S4; /* S3+x^2*S4. */ + sx = S2 + theta2 * sx; /* S2+x^2*(S3+x^2*S4). */ + sx = S1 + theta2 * sx; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */ + sx = S0 + theta2 * sx; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */ + sx = theta + theta * theta2 * sx; + } + else + { + /* Here sinf() is calculated using cos Chebyshev polynomial: + 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ + sx = C3 + theta2 * C4; /* C3+x^2*C4. */ + sx = C2 + theta2 * sx; /* C2+x^2*(C3+x^2*C4). */ + sx = C1 + theta2 * sx; /* C1+x^2*(C2+x^2*(C3+x^2*C4)). */ + sx = C0 + theta2 * sx; /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))). */ + sx = 1.0 + theta2 * sx; + } + + /* Add in the signbit and assign the result. */ + return sign * sx; +} + +/* Compute the cosine value using Chebyshev polynomials where + THETA is the range reduced absolute value of the input + and it is less than Pi/4, + N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide + whether a sine or cosine approximation is more accurate and + the sign of the result. */ +static inline float +reduced_cos (double theta, unsigned int n) +{ + double sign, cx; + const double theta2 = theta * theta; + + /* Determine positive or negative primary interval. */ + n += 2; + sign = ones[(n >> 2) & 1]; + + /* Are we in the primary interval of sin or cos? */ + if ((n & 2) == 0) + { + /* Here cosf() is calculated using sin Chebyshev polynomial: + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + cx = S3 + theta2 * S4; + cx = S2 + theta2 * cx; + cx = S1 + theta2 * cx; + cx = S0 + theta2 * cx; + cx = theta + theta * theta2 * cx; + } + else + { + /* Here cosf() is calculated using cos Chebyshev polynomial: + 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */ + cx = C3 + theta2 * C4; + cx = C2 + theta2 * cx; + cx = C1 + theta2 * cx; + cx = C0 + theta2 * cx; + cx = 1. + theta2 * cx; + } + return sign * cx; +} diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c index 916e345571..138e318dcc 100644 --- a/sysdeps/ieee754/flt-32/s_sinf.c +++ b/sysdeps/ieee754/flt-32/s_sinf.c @@ -1,25 +1,26 @@ -/* s_sinf.c -- float version of s_sin.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ +/* Compute sine of argument. + Copyright (C) 2017-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. -/* - * ==================================================== - * 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. - * ==================================================== - */ + 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. -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_sinf.c,v 1.4 1995/05/10 20:48:16 jtc Exp $"; -#endif + 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 <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> +#include "s_sincosf.h" #ifndef SINF # define SINF_FUNC __sinf @@ -27,37 +28,129 @@ static char rcsid[] = "$NetBSD: s_sinf.c,v 1.4 1995/05/10 20:48:16 jtc Exp $"; # define SINF_FUNC SINF #endif -float SINF_FUNC(float x) +float +SINF_FUNC (float x) { - float y[2],z=0.0; - int32_t n, ix; - - GET_FLOAT_WORD(ix,x); - - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if(ix <= 0x3f490fd8) return __kernel_sinf(x,z,0); - - /* sin(Inf or NaN) is NaN */ - else if (ix>=0x7f800000) { - if (ix == 0x7f800000) - __set_errno (EDOM); - return x-x; + double cx; + double theta = x; + double abstheta = fabs (theta); + /* If |x|< Pi/4. */ + if (isless (abstheta, M_PI_4)) + { + if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ + { + const double theta2 = theta * theta; + /* Chebyshev polynomial of the form for sin + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + cx = S3 + theta2 * S4; + cx = S2 + theta2 * cx; + cx = S1 + theta2 * cx; + cx = S0 + theta2 * cx; + cx = theta + theta * theta2 * cx; + return cx; } - - /* argument reduction needed */ - else { - n = __ieee754_rem_pio2f(x,y); - switch(n&3) { - case 0: return __kernel_sinf(y[0],y[1],1); - case 1: return __kernel_cosf(y[0],y[1]); - case 2: return -__kernel_sinf(y[0],y[1],1); - default: - return -__kernel_cosf(y[0],y[1]); + else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ + { + /* A simpler Chebyshev approximation is close enough for this range: + for sin: x+x^3*(SS0+x^2*SS1). */ + const double theta2 = theta * theta; + cx = SS0 + theta2 * SS1; + cx = theta + theta * theta2 * cx; + return cx; + } + else + { + /* Handle some special cases. */ + if (theta) + return theta - (theta * SMALL); + else + return theta; + } + } + else /* |x| >= Pi/4. */ + { + unsigned int signbit = isless (x, 0); + if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ + { + /* There are cases where FE_UPWARD rounding mode can + produce a result of abstheta * inv_PI_4 == 9, + where abstheta < 9pi/4, so the domain for + pio2_table must go to 5 (9 / 2 + 1). */ + unsigned int n = (abstheta * inv_PI_4) + 1; + theta = abstheta - pio2_table[n / 2]; + return reduced_sin (theta, n, signbit); + } + else if (isless (abstheta, INFINITY)) + { + if (abstheta < 0x1p+23) /* |x| < 2^23. */ + { + unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; + double x = n / 2; + theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; + /* Argument reduction needed. */ + return reduced_sin (theta, n, signbit); + } + else /* |x| >= 2^23. */ + { + x = fabsf (x); + int exponent; + GET_FLOAT_WORD (exponent, x); + exponent + = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; + exponent += 3; + exponent /= 28; + double a = invpio4_table[exponent] * x; + double b = invpio4_table[exponent + 1] * x; + double c = invpio4_table[exponent + 2] * x; + double d = invpio4_table[exponent + 3] * x; + uint64_t l = a; + l &= ~0x7; + a -= l; + double e = a + b; + l = e; + e = a - l; + if (l & 1) + { + e -= 1.0; + e += b; + e += c; + e += d; + e *= M_PI_4; + return reduced_sin (e, l + 1, signbit); + } + else + { + e += b; + e += c; + e += d; + if (e <= 1.0) + { + e *= M_PI_4; + return reduced_sin (e, l + 1, signbit); + } + else + { + l++; + e -= 2.0; + e *= M_PI_4; + return reduced_sin (e, l + 1, signbit); + } + } } } + else + { + int32_t ix; + /* High word of x. */ + GET_FLOAT_WORD (ix, abstheta); + /* Sin(Inf or NaN) is NaN. */ + if (ix == 0x7f800000) + __set_errno (EDOM); + return x - x; + } + } } #ifndef SINF -weak_alias (__sinf, sinf) +libm_alias_float (__sin, sin) #endif diff --git a/sysdeps/ieee754/flt-32/s_tanf.c b/sysdeps/ieee754/flt-32/s_tanf.c index 685df8fa35..ba3af54913 100644 --- a/sysdeps/ieee754/flt-32/s_tanf.c +++ b/sysdeps/ieee754/flt-32/s_tanf.c @@ -20,6 +20,7 @@ static char rcsid[] = "$NetBSD: s_tanf.c,v 1.4 1995/05/10 20:48:20 jtc Exp $"; #include <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> float __tanf(float x) { @@ -46,4 +47,4 @@ float __tanf(float x) -1 -- n odd */ } } -weak_alias (__tanf, tanf) +libm_alias_float (__tan, tan) diff --git a/sysdeps/ieee754/flt-32/s_tanhf.c b/sysdeps/ieee754/flt-32/s_tanhf.c index f70702b29c..cc3d63d35d 100644 --- a/sysdeps/ieee754/flt-32/s_tanhf.c +++ b/sysdeps/ieee754/flt-32/s_tanhf.c @@ -20,6 +20,8 @@ static char rcsid[] = "$NetBSD: s_tanhf.c,v 1.4 1995/05/10 20:48:24 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-float.h> static const float one=1.0, two=2.0, tiny = 1.0e-30; @@ -59,4 +61,4 @@ float __tanhf(float x) } return (jx>=0)? z: -z; } -weak_alias (__tanhf, tanhf) +libm_alias_float (__tanh, tanh) diff --git a/sysdeps/ieee754/flt-32/s_totalorderf.c b/sysdeps/ieee754/flt-32/s_totalorderf.c new file mode 100644 index 0000000000..928a8d6d8f --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_totalorderf.c @@ -0,0 +1,48 @@ +/* Total order operation. flt-32 version. + Copyright (C) 2016-2018 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 <libm-alias-float.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalorderf (float x, float y) +{ + int32_t ix, iy; + GET_FLOAT_WORD (ix, x); + GET_FLOAT_WORD (iy, y); +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the arguments interpreted as + sign-magnitude integers. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if ((ix & 0x7fffffff) > 0x7f800000 && (iy & 0x7fffffff) > 0x7f800000) + { + ix ^= 0x00400000; + iy ^= 0x00400000; + } +#endif + uint32_t ix_sign = ix >> 31; + uint32_t iy_sign = iy >> 31; + ix ^= ix_sign >> 1; + iy ^= iy_sign >> 1; + return ix <= iy; +} +libm_alias_float (__totalorder, totalorder) diff --git a/sysdeps/ieee754/flt-32/s_totalordermagf.c b/sysdeps/ieee754/flt-32/s_totalordermagf.c new file mode 100644 index 0000000000..0e383653a5 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_totalordermagf.c @@ -0,0 +1,46 @@ +/* Total order operation on absolute values. flt-32 version. + Copyright (C) 2016-2018 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 <libm-alias-float.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalordermagf (float x, float y) +{ + uint32_t ix, iy; + GET_FLOAT_WORD (ix, x); + GET_FLOAT_WORD (iy, y); + ix &= 0x7fffffff; + iy &= 0x7fffffff; +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the absolute values of the + arguments. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if (ix > 0x7f800000 && iy > 0x7f800000) + { + ix ^= 0x00400000; + iy ^= 0x00400000; + } +#endif + return ix <= iy; +} +libm_alias_float (__totalordermag, totalordermag) diff --git a/sysdeps/ieee754/flt-32/s_truncf.c b/sysdeps/ieee754/flt-32/s_truncf.c index 43d35c7f6a..2e1464aeac 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-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-float.h> float @@ -48,4 +49,6 @@ __truncf (float x) return x; } -weak_alias (__truncf, truncf) +#ifndef __truncf +libm_alias_float (__trunc, trunc) +#endif diff --git a/sysdeps/ieee754/flt-32/s_ufromfpf.c b/sysdeps/ieee754/flt-32/s_ufromfpf.c new file mode 100644 index 0000000000..0d2b4493a8 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_ufromfpf.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 0 +#define FUNC __ufromfpf +#include <s_fromfpf_main.c> +libm_alias_float (__ufromfp, ufromfp) diff --git a/sysdeps/ieee754/flt-32/s_ufromfpxf.c b/sysdeps/ieee754/flt-32/s_ufromfpxf.c new file mode 100644 index 0000000000..81f56daf9d --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_ufromfpxf.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 1 +#define FUNC __ufromfpxf +#include <s_fromfpf_main.c> +libm_alias_float (__ufromfpx, ufromfpx) diff --git a/sysdeps/ieee754/flt-32/t_exp2f.h b/sysdeps/ieee754/flt-32/t_exp2f.h deleted file mode 100644 index 3045b82cd1..0000000000 --- a/sysdeps/ieee754/flt-32/t_exp2f.h +++ /dev/null @@ -1,351 +0,0 @@ -/* Accurate tables for exp2f(). - 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> - - 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/>. */ - -/* This table has the property that, for all integers -128 <= i <= 127, - exp(i/256.0 + __exp2f_deltatable[i-128]) == __exp2f_atable[i+128] + r - for some -2^-35 < r < 2^-35 (abs(r) < 2^-36 if i <= 0); and that - __exp2f_deltatable[i+128] == t * 2^-30 - for integer t so that abs(t) <= 43447 * 2^0. */ - -#define W30 (9.31322575e-10) -static const float __exp2f_deltatable[256] = { - -810*W30, 283*W30, -1514*W30, 1304*W30, - -1148*W30, -98*W30, -744*W30, -156*W30, - -419*W30, -155*W30, 474*W30, 167*W30, - -1984*W30, -826*W30, 692*W30, 781*W30, - -578*W30, -411*W30, -129*W30, -1500*W30, - 654*W30, -141*W30, -816*W30, -53*W30, - 148*W30, 493*W30, -2214*W30, 760*W30, - 260*W30, 750*W30, -1300*W30, 1424*W30, - -1445*W30, -339*W30, -680*W30, -349*W30, - -922*W30, 531*W30, 193*W30, -2892*W30, - 290*W30, -2145*W30, -276*W30, 485*W30, - -695*W30, 215*W30, -7093*W30, 412*W30, - -4596*W30, 367*W30, 592*W30, -615*W30, - -97*W30, -1066*W30, 972*W30, -226*W30, - -625*W30, -374*W30, -5647*W30, -180*W30, - 20349*W30, -447*W30, 111*W30, -4164*W30, - -87*W30, -21*W30, -251*W30, 66*W30, - -517*W30, 2093*W30, -263*W30, 182*W30, - -601*W30, 475*W30, -483*W30, -1251*W30, - -373*W30, 1471*W30, -92*W30, -215*W30, - -97*W30, -190*W30, 0*W30, -290*W30, - -2647*W30, 1940*W30, -582*W30, 28*W30, - 833*W30, 1493*W30, 34*W30, 321*W30, - 3327*W30, -35*W30, 177*W30, -135*W30, - -796*W30, -428*W30, 129*W30, 9332*W30, - -12*W30, -69*W30, -1743*W30, 6508*W30, - -60*W30, 359*W30, 43447*W30, 15*W30, - -23*W30, -305*W30, -375*W30, -652*W30, - 667*W30, 269*W30, -1575*W30, 185*W30, - -329*W30, 200*W30, 6002*W30, 163*W30, - -647*W30, 19*W30, -603*W30, -755*W30, - 742*W30, -438*W30, 3587*W30, 2560*W30, - 0*W30, -520*W30, -241*W30, -299*W30, - -1270*W30, -991*W30, -1138*W30, 255*W30, - -1192*W30, 1722*W30, 1023*W30, 3700*W30, - -1388*W30, -1551*W30, -2549*W30, 27*W30, - 282*W30, 673*W30, 113*W30, 1561*W30, - 72*W30, 873*W30, 87*W30, -395*W30, - -433*W30, 629*W30, 3440*W30, -284*W30, - -592*W30, -103*W30, -46*W30, -3844*W30, - 1712*W30, 303*W30, 1555*W30, -631*W30, - -1400*W30, -961*W30, -854*W30, -276*W30, - 407*W30, 833*W30, -345*W30, -1501*W30, - 121*W30, -1581*W30, 400*W30, 150*W30, - 1224*W30, -139*W30, -563*W30, 879*W30, - 933*W30, 2939*W30, 788*W30, 211*W30, - 530*W30, -192*W30, 706*W30, -13347*W30, - 1065*W30, 3*W30, 111*W30, -208*W30, - -360*W30, -532*W30, -291*W30, 483*W30, - 987*W30, -33*W30, -1373*W30, -166*W30, - -1174*W30, -3955*W30, 1601*W30, -280*W30, - 1405*W30, 600*W30, -1659*W30, -23*W30, - 390*W30, 449*W30, 570*W30, -13143*W30, - -9*W30, -1646*W30, 1201*W30, 294*W30, - 2181*W30, -1173*W30, 1388*W30, -4504*W30, - 190*W30, -2304*W30, 211*W30, 239*W30, - 48*W30, -817*W30, 1018*W30, 1828*W30, - -663*W30, 1408*W30, 408*W30, -36*W30, - 1295*W30, -230*W30, 1341*W30, 9*W30, - 40*W30, 705*W30, 186*W30, 376*W30, - 557*W30, 5866*W30, 363*W30, -1558*W30, - 718*W30, 669*W30, 1369*W30, -2972*W30, - -468*W30, -121*W30, -219*W30, 667*W30, - 29954*W30, 366*W30, 48*W30, -203*W30 -}; - -static const float __exp2f_atable[256] /* __attribute__((mode(SF))) */ = { - 0.707106411447, /* 0x0.b504ecfff */ - 0.709024071690, /* 0x0.b58299fff */ - 0.710945606239, /* 0x0.b60088000 */ - 0.712874472142, /* 0x0.b67ef1000 */ - 0.714806139464, /* 0x0.b6fd88fff */ - 0.716744661340, /* 0x0.b77c94000 */ - 0.718687653549, /* 0x0.b7fbea000 */ - 0.720636486992, /* 0x0.b87ba1fff */ - 0.722590208040, /* 0x0.b8fbabfff */ - 0.724549472323, /* 0x0.b97c12fff */ - 0.726514220228, /* 0x0.b9fcd5fff */ - 0.728483855735, /* 0x0.ba7deb000 */ - 0.730457961549, /* 0x0.baff4afff */ - 0.732438981522, /* 0x0.bb811efff */ - 0.734425544748, /* 0x0.bc0350000 */ - 0.736416816713, /* 0x0.bc85d0000 */ - 0.738412797450, /* 0x0.bd089efff */ - 0.740414917465, /* 0x0.bd8bd4fff */ - 0.742422521111, /* 0x0.be0f66fff */ - 0.744434773914, /* 0x0.be9346fff */ - 0.746454179287, /* 0x0.bf179f000 */ - 0.748477637755, /* 0x0.bf9c3afff */ - 0.750506639473, /* 0x0.c02133fff */ - 0.752541840064, /* 0x0.c0a694fff */ - 0.754582285889, /* 0x0.c12c4e000 */ - 0.756628334525, /* 0x0.c1b265000 */ - 0.758678436269, /* 0x0.c238bffff */ - 0.760736882681, /* 0x0.c2bfa6fff */ - 0.762799203401, /* 0x0.c346cf000 */ - 0.764867603790, /* 0x0.c3ce5d000 */ - 0.766940355298, /* 0x0.c45633fff */ - 0.769021093841, /* 0x0.c4de90fff */ - 0.771104693409, /* 0x0.c5671dfff */ - 0.773195922364, /* 0x0.c5f02afff */ - 0.775292098512, /* 0x0.c6798afff */ - 0.777394294745, /* 0x0.c70350000 */ - 0.779501736166, /* 0x0.c78d6d000 */ - 0.781615912910, /* 0x0.c817fafff */ - 0.783734917628, /* 0x0.c8a2d9fff */ - 0.785858273516, /* 0x0.c92e02000 */ - 0.787990570071, /* 0x0.c9b9c0000 */ - 0.790125787245, /* 0x0.ca45aefff */ - 0.792268991467, /* 0x0.cad223fff */ - 0.794417440881, /* 0x0.cb5ef0fff */ - 0.796570718287, /* 0x0.cbec0efff */ - 0.798730909811, /* 0x0.cc79a0fff */ - 0.800892710672, /* 0x0.cd074dfff */ - 0.803068041795, /* 0x0.cd95ddfff */ - 0.805242776881, /* 0x0.ce2464000 */ - 0.807428598393, /* 0x0.ceb3a3fff */ - 0.809617877002, /* 0x0.cf431dfff */ - 0.811812341211, /* 0x0.cfd2eefff */ - 0.814013659956, /* 0x0.d06333000 */ - 0.816220164311, /* 0x0.d0f3ce000 */ - 0.818434238424, /* 0x0.d184e7fff */ - 0.820652604094, /* 0x0.d21649fff */ - 0.822877407074, /* 0x0.d2a818000 */ - 0.825108587751, /* 0x0.d33a51000 */ - 0.827342867839, /* 0x0.d3ccbdfff */ - 0.829588949684, /* 0x0.d45ff1000 */ - 0.831849217401, /* 0x0.d4f411fff */ - 0.834093391880, /* 0x0.d58724fff */ - 0.836355149750, /* 0x0.d61b5f000 */ - 0.838620424257, /* 0x0.d6afd3fff */ - 0.840896368027, /* 0x0.d744fc000 */ - 0.843176305293, /* 0x0.d7da66fff */ - 0.845462262643, /* 0x0.d87037000 */ - 0.847754716864, /* 0x0.d90673fff */ - 0.850052893157, /* 0x0.d99d10fff */ - 0.852359056469, /* 0x0.da3433fff */ - 0.854668736446, /* 0x0.dacb91fff */ - 0.856986224651, /* 0x0.db6373000 */ - 0.859309315673, /* 0x0.dbfbb1fff */ - 0.861639738080, /* 0x0.dc946bfff */ - 0.863975346095, /* 0x0.dd2d7d000 */ - 0.866317391394, /* 0x0.ddc6f9fff */ - 0.868666708472, /* 0x0.de60f1000 */ - 0.871022939695, /* 0x0.defb5c000 */ - 0.873383641229, /* 0x0.df9611fff */ - 0.875751554968, /* 0x0.e03141000 */ - 0.878126025200, /* 0x0.e0ccde000 */ - 0.880506813521, /* 0x0.e168e4fff */ - 0.882894217966, /* 0x0.e2055afff */ - 0.885287821299, /* 0x0.e2a239000 */ - 0.887686729423, /* 0x0.e33f6ffff */ - 0.890096127973, /* 0x0.e3dd56fff */ - 0.892507970338, /* 0x0.e47b67000 */ - 0.894928157336, /* 0x0.e51a03000 */ - 0.897355020043, /* 0x0.e5b90efff */ - 0.899788379682, /* 0x0.e65888000 */ - 0.902227103705, /* 0x0.e6f85afff */ - 0.904673457151, /* 0x0.e798ae000 */ - 0.907128036008, /* 0x0.e8398afff */ - 0.909585535528, /* 0x0.e8da99000 */ - 0.912051796915, /* 0x0.e97c3a000 */ - 0.914524436003, /* 0x0.ea1e46000 */ - 0.917003571999, /* 0x0.eac0bf000 */ - 0.919490039339, /* 0x0.eb63b2fff */ - 0.921983361257, /* 0x0.ec071a000 */ - 0.924488604054, /* 0x0.ecab48fff */ - 0.926989555360, /* 0x0.ed4f30000 */ - 0.929502844812, /* 0x0.edf3e6000 */ - 0.932021975503, /* 0x0.ee98fdfff */ - 0.934553921208, /* 0x0.ef3eecfff */ - 0.937083780759, /* 0x0.efe4b8fff */ - 0.939624726786, /* 0x0.f08b3f000 */ - 0.942198514924, /* 0x0.f133ebfff */ - 0.944726586343, /* 0x0.f1d99a000 */ - 0.947287976728, /* 0x0.f28176fff */ - 0.949856162070, /* 0x0.f329c5fff */ - 0.952431440345, /* 0x0.f3d28bfff */ - 0.955013573175, /* 0x0.f47bc5000 */ - 0.957603693021, /* 0x0.f52584000 */ - 0.960199773321, /* 0x0.f5cfa7000 */ - 0.962801992906, /* 0x0.f67a31000 */ - 0.965413510788, /* 0x0.f72556fff */ - 0.968030691152, /* 0x0.f7d0dc000 */ - 0.970655620084, /* 0x0.f87ce2fff */ - 0.973290979849, /* 0x0.f92998fff */ - 0.975926160805, /* 0x0.f9d64bfff */ - 0.978571653370, /* 0x0.fa83ac000 */ - 0.981225252139, /* 0x0.fb3193fff */ - 0.983885228626, /* 0x0.fbdfe6fff */ - 0.986552715296, /* 0x0.fc8eb7fff */ - 0.989228487027, /* 0x0.fd3e14000 */ - 0.991909801964, /* 0x0.fdedcd000 */ - 0.994601726545, /* 0x0.fe9e38000 */ - 0.997297704209, /* 0x0.ff4ee6fff */ - 1.000000000000, /* 0x1.000000000 */ - 1.002710938457, /* 0x1.00b1aa000 */ - 1.005429744692, /* 0x1.0163d7ffe */ - 1.008155703526, /* 0x1.02167dffe */ - 1.010888457284, /* 0x1.02c995fff */ - 1.013629436498, /* 0x1.037d38000 */ - 1.016377568250, /* 0x1.043152000 */ - 1.019134163841, /* 0x1.04e5f9ffe */ - 1.021896362316, /* 0x1.059b00000 */ - 1.024668931945, /* 0x1.0650b3ffe */ - 1.027446627635, /* 0x1.0706be001 */ - 1.030234098408, /* 0x1.07bd6bffe */ - 1.033023953416, /* 0x1.087441ffe */ - 1.035824656494, /* 0x1.092bce000 */ - 1.038632392900, /* 0x1.09e3d0001 */ - 1.041450142840, /* 0x1.0a9c79ffe */ - 1.044273972530, /* 0x1.0b558a001 */ - 1.047105550795, /* 0x1.0c0f1c001 */ - 1.049944162390, /* 0x1.0cc924001 */ - 1.052791833895, /* 0x1.0d83c4001 */ - 1.055645227426, /* 0x1.0e3ec3fff */ - 1.058507919326, /* 0x1.0efa60001 */ - 1.061377286898, /* 0x1.0fb66bfff */ - 1.064254641510, /* 0x1.1072fdffe */ - 1.067140102389, /* 0x1.113018000 */ - 1.070034146304, /* 0x1.11edc1fff */ - 1.072937250162, /* 0x1.12ac04001 */ - 1.075843691823, /* 0x1.136a7dfff */ - 1.078760385496, /* 0x1.1429a3ffe */ - 1.081685543070, /* 0x1.14e958000 */ - 1.084618330005, /* 0x1.15a98c000 */ - 1.087556362176, /* 0x1.166a18001 */ - 1.090508937863, /* 0x1.172b98001 */ - 1.093464612954, /* 0x1.17ed4bfff */ - 1.096430182434, /* 0x1.18afa5ffe */ - 1.099401354802, /* 0x1.19725e000 */ - 1.102381587017, /* 0x1.1a35adfff */ - 1.105370759965, /* 0x1.1af994000 */ - 1.108367800686, /* 0x1.1bbdfdffe */ - 1.111373305331, /* 0x1.1c82f6000 */ - 1.114387035385, /* 0x1.1d4878001 */ - 1.117408752440, /* 0x1.1e0e7ffff */ - 1.120437502874, /* 0x1.1ed4fe000 */ - 1.123474478729, /* 0x1.1f9c06000 */ - 1.126521706601, /* 0x1.2063ba001 */ - 1.129574775716, /* 0x1.212bd0001 */ - 1.132638812065, /* 0x1.21f49e000 */ - 1.135709524130, /* 0x1.22bddbffe */ - 1.138789534565, /* 0x1.2387b5fff */ - 1.141876101508, /* 0x1.2451fe000 */ - 1.144971728301, /* 0x1.251cddffe */ - 1.148077130296, /* 0x1.25e861ffe */ - 1.151189923305, /* 0x1.26b462001 */ - 1.154312610610, /* 0x1.278107ffe */ - 1.157440662410, /* 0x1.284e08001 */ - 1.160578370109, /* 0x1.291baa001 */ - 1.163725256932, /* 0x1.29e9e6000 */ - 1.166879892324, /* 0x1.2ab8a3ffe */ - 1.170044302935, /* 0x1.2b8805fff */ - 1.173205971694, /* 0x1.2c5739ffe */ - 1.176397800428, /* 0x1.2d2867ffe */ - 1.179586529747, /* 0x1.2df962001 */ - 1.182784795737, /* 0x1.2ecafbffe */ - 1.185991406414, /* 0x1.2f9d21ffe */ - 1.189206838636, /* 0x1.306fdc001 */ - 1.192430973067, /* 0x1.314328000 */ - 1.195664167430, /* 0x1.32170c001 */ - 1.198906540890, /* 0x1.32eb8a001 */ - 1.202157497408, /* 0x1.33c098000 */ - 1.205416083326, /* 0x1.349625fff */ - 1.208683252332, /* 0x1.356c43fff */ - 1.211961269402, /* 0x1.364318001 */ - 1.215246438983, /* 0x1.371a64000 */ - 1.218539118740, /* 0x1.37f22dffe */ - 1.221847295770, /* 0x1.38cafc000 */ - 1.225158572187, /* 0x1.39a3fdfff */ - 1.228481650325, /* 0x1.3a7dc5ffe */ - 1.231811761846, /* 0x1.3b5803fff */ - 1.235149741144, /* 0x1.3c32c5ffe */ - 1.238499879811, /* 0x1.3d0e53ffe */ - 1.241858124726, /* 0x1.3dea69fff */ - 1.245225191102, /* 0x1.3ec713fff */ - 1.248601436624, /* 0x1.3fa458000 */ - 1.251975655584, /* 0x1.40817a001 */ - 1.255380749731, /* 0x1.4160a2001 */ - 1.258783102010, /* 0x1.423f9bffe */ - 1.262198328973, /* 0x1.431f6e000 */ - 1.265619754780, /* 0x1.43ffa7fff */ - 1.269052743928, /* 0x1.44e0a4001 */ - 1.272490739830, /* 0x1.45c1f4000 */ - 1.275942921659, /* 0x1.46a432001 */ - 1.279397487615, /* 0x1.478697ffe */ - 1.282870173427, /* 0x1.486a2dffe */ - 1.286346316319, /* 0x1.494dfdffe */ - 1.289836049094, /* 0x1.4a32b2001 */ - 1.293333172770, /* 0x1.4b17e1ffe */ - 1.296839594835, /* 0x1.4bfdadfff */ - 1.300354957560, /* 0x1.4ce40fffe */ - 1.303882122055, /* 0x1.4dcb38001 */ - 1.307417988757, /* 0x1.4eb2f1ffe */ - 1.310960650439, /* 0x1.4f9b1dfff */ - 1.314516782746, /* 0x1.50842bfff */ - 1.318079948424, /* 0x1.516daffff */ - 1.321653246888, /* 0x1.5257de000 */ - 1.325237751030, /* 0x1.5342c8001 */ - 1.328829526907, /* 0x1.542e2c000 */ - 1.332433700535, /* 0x1.551a5fffe */ - 1.336045145966, /* 0x1.56070dffe */ - 1.339667558645, /* 0x1.56f473ffe */ - 1.343300342533, /* 0x1.57e287ffe */ - 1.346941947961, /* 0x1.58d130001 */ - 1.350594043714, /* 0x1.59c087ffe */ - 1.354256033883, /* 0x1.5ab085fff */ - 1.357932448365, /* 0x1.5ba175ffe */ - 1.361609339707, /* 0x1.5c926dfff */ - 1.365299344044, /* 0x1.5d8441ffe */ - 1.369003057507, /* 0x1.5e76fc001 */ - 1.372714757920, /* 0x1.5f6a3c000 */ - 1.376437187179, /* 0x1.605e2fffe */ - 1.380165219333, /* 0x1.615282001 */ - 1.383909463864, /* 0x1.6247e3ffe */ - 1.387661933907, /* 0x1.633dd0000 */ - 1.391424179060, /* 0x1.64345fffe */ - 1.395197510706, /* 0x1.652ba9fff */ - 1.399006724329, /* 0x1.66254dffe */ - 1.402773022651, /* 0x1.671c22000 */ - 1.406576037403, /* 0x1.68155dfff */ - 1.410389423392, /* 0x1.690f48001 */ -}; diff --git a/sysdeps/ieee754/flt-32/w_exp2f.c b/sysdeps/ieee754/flt-32/w_exp2f.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_exp2f.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/w_expf.c b/sysdeps/ieee754/flt-32/w_expf.c index ed1550972f..1cc8931700 100644 --- a/sysdeps/ieee754/flt-32/w_expf.c +++ b/sysdeps/ieee754/flt-32/w_expf.c @@ -1,34 +1 @@ -/* 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. - - 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> - -/* wrapper expf */ -float -__expf (float x) -{ - float z = __ieee754_expf (x); - if (__builtin_expect (!isfinite (z) || z == 0, 0) - && isfinite (x) && _LIB_VERSION != _IEEE_) - return __kernel_standard_f (x, x, 106 + !!signbit (x)); - - return z; -} -hidden_def (__expf) -weak_alias (__expf, expf) +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/w_log2f.c b/sysdeps/ieee754/flt-32/w_log2f.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_log2f.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/w_logf.c b/sysdeps/ieee754/flt-32/w_logf.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_logf.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ieee754/flt-32/w_powf.c b/sysdeps/ieee754/flt-32/w_powf.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ieee754/flt-32/w_powf.c @@ -0,0 +1 @@ +/* Not needed. */ |