mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-26 23:10:06 +00:00
43cfdf8f48
This patch currently only affects aarch64. The roundtoint and converttoint internal functions are only called with small values, so 32 bit result is enough for converttoint and it is a signed int conversion so the return type is changed to int32_t. The original idea was to help the compiler keeping the result in uint64_t, then it's clear that no sign extension is needed and there is no accidental undefined or implementation defined signed int arithmetics. But it turns out gcc does a good job with inlining so changing the type has no overhead and the semantics of the conversion is less surprising this way. Since we want to allow the asuint64 (x + 0x1.8p52) style conversion, the top bits were never usable and the existing code ensures that only the bottom 32 bits of the conversion result are used. On aarch64 the neon intrinsics (which round ties to even) are changed to round and lround (which round ties away from zero) this does not affect the results in a significant way, but more portable (relies on round and lround being inlined which works with -fno-math-errno). The TOINT_SHIFT and TOINT_RINT macros were removed, only keep separate code paths for TOINT_INTRINSICS and !TOINT_INTRINSICS. * sysdeps/aarch64/fpu/math_private.h (roundtoint): Use round. (converttoint): Use lround. * sysdeps/ieee754/flt-32/math_config.h (roundtoint): Declare and document the semantics when TOINT_INTRINSICS is set. (converttoint): Likewise. (TOINT_RINT): Remove. (TOINT_SHIFT): Remove. * sysdeps/ieee754/flt-32/e_expf.c (__expf): Remove the TOINT_RINT code path.
115 lines
3.2 KiB
C
115 lines
3.2 KiB
C
/* Single-precision e^x function.
|
|
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/>. */
|
|
|
|
#ifdef __expf
|
|
# undef libm_hidden_proto
|
|
# define libm_hidden_proto(ignored)
|
|
#endif
|
|
|
|
#include <math.h>
|
|
#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
|
|
__expf (float x)
|
|
{
|
|
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)))
|
|
{
|
|
/* |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
|
|
}
|
|
|
|
/* 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);
|
|
#else
|
|
# 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;
|
|
}
|
|
|
|
#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
|