@@ -115,7 +115,7 @@ type-double-routines := branred doasin dosincos halfulp mpa mpatan2 \
# float support
type-float-suffix := f
-type-float-routines := k_rem_pio2f
+type-float-routines := k_rem_pio2f math_errf e_exp2f_data
# _Float128 support
type-float128-suffix := f128
new file mode 100644
@@ -0,0 +1,35 @@
+/* Copyright (C) 2011-2017 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>
+#include <math-svid-compat.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)
@@ -319,6 +319,26 @@ libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx)
#define libc_feresetround_noexf_ctx libc_feresetround_noex_aarch64_ctx
#define libc_feresetround_noexl_ctx libc_feresetround_noex_aarch64_ctx
+/* Hack: only include the large arm_neon.h when needed. */
+#ifdef _MATH_CONFIG_H
+#include <arm_neon.h>
+
+/* ACLE intrinsics for frintn and fcvtns instructions. */
+#define TOINT_INTRINSICS 1
+
+static inline double_t
+roundtoint (double_t x)
+{
+ return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0);
+}
+
+static inline uint64_t
+converttoint (double_t x)
+{
+ return vcvtnd_s64_f64 (x);
+}
+#endif
+
#include_next <math_private.h>
#endif
@@ -1,7 +1,6 @@
-/* Single-precision floating point 2^x.
- Copyright (C) 1997-2017 Free Software Foundation, Inc.
+/* Single-precision 2^x function.
+ Copyright (C) 2017 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,75 @@
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 <stdint.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))
+ 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 (__builtin_expect (abstop >= top12 (128.0f), 0))
{
- 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))
- {
- 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 = (double) (xd + SHIFT); /* Rounding to double precision is required. */
+ 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)
+weak_alias (__exp2f, exp2f)
+strong_alias (__exp2f, __ieee754_exp2f)
+strong_alias (__exp2f, __exp2f_finite)
new file mode 100644
@@ -0,0 +1,44 @@
+/* Shared data between expf, exp2f and powf.
+ Copyright (C) 2017 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
+ },
+};
@@ -1,7 +1,6 @@
-/* Single-precision floating point e^x.
- Copyright (C) 1997-2017 Free Software Foundation, Inc.
+/* Single-precision e^x function.
+ Copyright (C) 2017 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,90 @@
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
-/* How this works:
-
- 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 <stdint.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))
+ 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 (__builtin_expect (abstop >= top12 (88.0f), 0))
{
- 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;
+ /* |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
}
- /* Exceptional cases: */
- else if (isless (x, himark))
- {
- if (isinf (x))
- /* e^-inf == 0, with no error. */
- return 0;
- else
- /* Underflow */
- return TWOM100 * TWOM100;
- }
- 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 = (double) (z + SHIFT); /* Rounding to double precision is required. */
+ 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)
+hidden_def (__expf)
+weak_alias (__expf, expf)
+strong_alias (__expf, __ieee754_expf)
+strong_alias (__expf, __expf_finite)
new file mode 100644
@@ -0,0 +1,128 @@
+/* Configuration for math routines.
+ Copyright (C) 2017 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 <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
+ieee_2008_issignaling (float x)
+{
+ uint32_t ix = asuint (x);
+ ix ^= 0x00400000; /* IEEE 754-2008 snan bit. */
+ return 2 * ix > 2u * 0x7fc00000;
+}
+
+#ifdef __GNUC__
+#define HIDDEN __attribute__ ((__visibility__ ("hidden")))
+#define NOINLINE __attribute__ ((noinline))
+#else
+#define HIDDEN
+#define NOINLINE
+#endif
+
+HIDDEN float __math_oflowf (unsigned long);
+HIDDEN float __math_uflowf (unsigned long);
+HIDDEN float __math_may_uflowf (unsigned long);
+HIDDEN float __math_divzerof (unsigned long);
+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;
+
+#endif
new file mode 100644
@@ -0,0 +1,77 @@
+/* Single-precision math error handling.
+ Copyright (C) 2017 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 and avoids making math functions non-leaf
+ when the error handling is inlined. */
+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 (unsigned long sign, float y)
+{
+ y = (sign ? -y : y) * y;
+ return with_errnof (y, ERANGE);
+}
+
+HIDDEN float
+__math_uflowf (unsigned long 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. */
+HIDDEN float
+__math_may_uflowf (unsigned long sign)
+{
+ return xflowf (sign, 0x1.4p-75f);
+}
+#endif
+
+HIDDEN float
+__math_oflowf (unsigned long sign)
+{
+ return xflowf (sign, 0x1p97f);
+}
+
+HIDDEN float
+__math_divzerof (unsigned long sign)
+{
+ float y = 0;
+ return with_errnof ((sign ? -1 : 1) / y, ERANGE);
+}
+
+HIDDEN float
+__math_invalidf (float x)
+{
+ float y = (x - x) / (x - x);
+ return isnan (x) ? y : with_errnof (y, EDOM);
+}
new file mode 100644
@@ -0,0 +1 @@
+/* Empty file to disable the exp2f wrapper. */
@@ -1,35 +1 @@
-/* Copyright (C) 2011-2017 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>
-#include <math-svid-compat.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)
+/* Empty file to disable the expf wrapper. */
new file mode 100644
@@ -0,0 +1 @@
+#include <math/w_expf_compat.c>