[PATCH 2/5] math: Optimized generic exp10f with wrappers
Adhemerval Zanella
adhemerval.zanella@linaro.org
Wed May 20 20:18:56 GMT 2020
Ping (x2).
On 29/04/2020 14:11, Adhemerval Zanella wrote:
> Ping.
>
> On 09/04/2020 16:59, Adhemerval Zanella wrote:
>> From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
>>
>> It is inspired by expf and reuses its tables and internal functions.
>> The error checks are inlined and errno setting is in separate tail
>> called functions, but the wrappers are kept in this patch to handle
>> the _LIB_VERSION==_SVID_ case.
>>
>> Double precision arithmetics is used which is expected to be faster on
>> most targets (including soft-float) than using single precision and it
>> is easier to get good precision result with it.
>>
>> Result for x86_64 (i7-4790K CPU @ 4.00GHz) are:
>>
>> Before new code:
>> "exp10f": {
>> "": {
>> "duration": 1.00923e+09,
>> "iterations": 5.4832e+07,
>> "max": 249.361,
>> "min": 12.5,
>> "mean": 18.4059
>> }
>>
>> With new code:
>> "exp10f": {
>> "": {
>> "duration": 9.7639e+08,
>> "iterations": 8.1056e+07,
>> "max": 563.323,
>> "min": 11.54,
>> "mean": 12.0459
>> }
>>
>> Result for aarch64 (A72 @ 2GHz) are:
>>
>> Before new code:
>> "exp10f": {
>> "": {
>> "duration": 1.00923e+09,
>> "iterations": 5.4832e+07,
>> "max": 249.361,
>> "min": 12.5,
>> "mean": 18.4059
>> }
>>
>> With new code:
>> "exp10f": {
>> "": {
>> "duration": 9.7639e+08,
>> "iterations": 8.1056e+07,
>> "max": 563.323,
>> "min": 11.54,
>> "mean": 12.0459
>> }
>>
>> Checked on x86_64-linux-gnu, powerpc64le-linux-gnu, aarch64-linux-gnu,
>> and sparc64-linux-gnu.
>> ---
>> math/e_exp10f.c | 32 -----
>> sysdeps/ieee754/flt-32/e_exp10f.c | 198 +++++++++++++++++++++++++++
>> sysdeps/ieee754/flt-32/math_config.h | 2 +-
>> 3 files changed, 199 insertions(+), 33 deletions(-)
>> delete mode 100644 math/e_exp10f.c
>> create mode 100644 sysdeps/ieee754/flt-32/e_exp10f.c
>>
>> diff --git a/math/e_exp10f.c b/math/e_exp10f.c
>> deleted file mode 100644
>> index 93c41d00a6..0000000000
>> --- a/math/e_exp10f.c
>> +++ /dev/null
>> @@ -1,32 +0,0 @@
>> -/* Copyright (C) 1998-2020 Free Software Foundation, Inc.
>> - This file is part of the GNU C Library.
>> - Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998.
>> -
>> - 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
>> - <https://www.gnu.org/licenses/>. */
>> -
>> -#include <math.h>
>> -#include <math_private.h>
>> -#include <libm-alias-finite.h>
>> -
>> -float
>> -__ieee754_exp10f (float arg)
>> -{
>> - /* The argument to exp needs to be calculated in enough precision
>> - that the fractional part has as much precision as float, in
>> - addition to the bits in the integer part; using double ensures
>> - this. */
>> - return __ieee754_exp (M_LN10 * arg);
>> -}
>> -libm_alias_finite (__ieee754_exp10f, __exp10f)
>> diff --git a/sysdeps/ieee754/flt-32/e_exp10f.c b/sysdeps/ieee754/flt-32/e_exp10f.c
>> new file mode 100644
>> index 0000000000..034a9e364a
>> --- /dev/null
>> +++ b/sysdeps/ieee754/flt-32/e_exp10f.c
>> @@ -0,0 +1,198 @@
>> +/* Single-precision 10^x function.
>> + Copyright (C) 2020 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
>> + <https://www.gnu.org/licenses/>. */
>> +
>> +#include <math.h>
>> +#include <math-narrow-eval.h>
>> +#include <stdint.h>
>> +#include <libm-alias-finite.h>
>> +#include <libm-alias-float.h>
>> +#include "math_config.h"
>> +
>> +/*
>> + EXP2F_TABLE_BITS 5
>> + EXP2F_POLY_ORDER 3
>> +
>> + Max. ULP error: 0.502 (normal range, nearest rounding).
>> + Max. relative error: 2^-33.240 (before rounding to float).
>> + Wrong count: 169839.
>> + Non-nearest ULP error: 1 (rounded ULP error).
>> +
>> + Detailed error analysis (for EXP2F_TABLE_BITS=3 thus N=32):
>> +
>> + - We first compute z = RN(InvLn10N * x) where
>> +
>> + InvLn10N = N*log(10)/log(2) * (1 + theta1) with |theta1| < 2^-54.150
>> +
>> + since z is rounded to nearest:
>> +
>> + z = InvLn10N * x * (1 + theta2) with |theta2| < 2^-53
>> +
>> + thus z = N*log(10)/log(2) * x * (1 + theta3) with |theta3| < 2^-52.463
>> +
>> + - Since |x| < 38 in the main branch, we deduce:
>> +
>> + z = N*log(10)/log(2) * x + theta4 with |theta4| < 2^-40.483
>> +
>> + - We then write z = k + r where k is an integer and |r| <= 0.5 (exact)
>> +
>> + - We thus have
>> +
>> + x = z*log(2)/(N*log(10)) - theta4*log(2)/(N*log(10))
>> + = z*log(2)/(N*log(10)) + theta5 with |theta5| < 2^-47.215
>> +
>> + 10^x = 2^(k/N) * 2^(r/N) * 10^theta5
>> + = 2^(k/N) * 2^(r/N) * (1 + theta6) with |theta6| < 2^-46.011
>> +
>> + - Then 2^(k/N) is approximated by table lookup, the maximal relative error
>> + is for (k%N) = 5, with
>> +
>> + s = 2^(i/N) * (1 + theta7) with |theta7| < 2^-53.249
>> +
>> + - 2^(r/N) is approximated by a degree-3 polynomial, where the maximal
>> + mathematical relative error is 2^-33.243.
>> +
>> + - For C[0] * r + C[1], assuming no FMA is used, since |r| <= 0.5 and
>> + |C[0]| < 1.694e-6, |C[0] * r| < 8.47e-7, and the absolute error on
>> + C[0] * r is bounded by 1/2*ulp(8.47e-7) = 2^-74. Then we add C[1] with
>> + |C[1]| < 0.000235, thus the absolute error on C[0] * r + C[1] is bounded
>> + by 2^-65.994 (z is bounded by 0.000236).
>> +
>> + - For r2 = r * r, since |r| <= 0.5, we have |r2| <= 0.25 and the absolute
>> + error is bounded by 1/4*ulp(0.25) = 2^-56 (the factor 1/4 is because |r2|
>> + cannot exceed 1/4, and on the left of 1/4 the distance between two
>> + consecutive numbers is 1/4*ulp(1/4)).
>> +
>> + - For y = C[2] * r + 1, assuming no FMA is used, since |r| <= 0.5 and
>> + |C[2]| < 0.0217, the absolute error on C[2] * r is bounded by 2^-60,
>> + and thus the absolute error on C[2] * r + 1 is bounded by 1/2*ulp(1)+2^60
>> + < 2^-52.988, and |y| < 1.01085 (the error bound is better if a FMA is
>> + used).
>> +
>> + - for z * r2 + y: the absolute error on z is bounded by 2^-65.994, with
>> + |z| < 0.000236, and the absolute error on r2 is bounded by 2^-56, with
>> + r2 < 0.25, thus |z*r2| < 0.000059, and the absolute error on z * r2
>> + (including the rounding error) is bounded by:
>> +
>> + 2^-65.994 * 0.25 + 0.000236 * 2^-56 + 1/2*ulp(0.000059) < 2^-66.429.
>> +
>> + Now we add y, with |y| < 1.01085, and error on y bounded by 2^-52.988,
>> + thus the absolute error is bounded by:
>> +
>> + 2^-66.429 + 2^-52.988 + 1/2*ulp(1.01085) < 2^-51.993
>> +
>> + - Now we convert the error on y into relative error. Recall that y
>> + approximates 2^(r/N), for |r| <= 0.5 and N=32. Thus 2^(-0.5/32) <= y,
>> + and the relative error on y is bounded by
>> +
>> + 2^-51.993/2^(-0.5/32) < 2^-51.977
>> +
>> + - Taking into account the mathematical relative error of 2^-33.243 we have:
>> +
>> + y = 2^(r/N) * (1 + theta8) with |theta8| < 2^-33.242
>> +
>> + - Since we had s = 2^(k/N) * (1 + theta7) with |theta7| < 2^-53.249,
>> + after y = y * s we get y = 2^(k/N) * 2^(r/N) * (1 + theta9)
>> + with |theta9| < 2^-33.241
>> +
>> + - Finally, taking into account the error theta6 due to the rounding error on
>> + InvLn10N, and when multiplying InvLn10N * x, we get:
>> +
>> + y = 10^x * (1 + theta10) with |theta10| < 2^-33.240
>> +
>> + - Converting into binary64 ulps, since y < 2^53*ulp(y), the error is at most
>> + 2^19.76 ulp(y)
>> +
>> + - If the result is a binary32 value in the normal range (i.e., >= 2^-126),
>> + then the error is at most 2^-9.24 ulps, i.e., less than 0.00166 (in the
>> + subnormal range, the error in ulps might be larger).
>> +
>> + Note that this bound is tight, since for x=-0x25.54ac0p0 the final value of
>> + y (before conversion to float) differs from 879853 ulps from the correctly
>> + rounded value, and 879853 ~ 2^19.7469. */
>> +
>> +#define N (1 << EXP2F_TABLE_BITS)
>> +
>> +#define InvLn10N (0x3.5269e12f346e2p0 * N) /* log(10)/log(2) to nearest */
>> +#define T __exp2f_data.tab
>> +#define C __exp2f_data.poly_scaled
>> +#define SHIFT __exp2f_data.shift
>> +
>> +static inline uint32_t
>> +top13 (float x)
>> +{
>> + return asuint (x) >> 19;
>> +}
>> +
>> +float
>> +__ieee754_exp10f (float x)
>> +{
>> + uint32_t abstop;
>> + uint64_t ki, t;
>> + double kd, xd, z, r, r2, y, s;
>> +
>> + xd = (double) x;
>> + abstop = top13 (x) & 0xfff; /* Ignore sign. */
>> + if (__glibc_unlikely (abstop >= top13 (38.0f)))
>> + {
>> + /* |x| >= 38 or x is nan. */
>> + if (asuint (x) == asuint (-INFINITY))
>> + return 0.0f;
>> + if (abstop >= top13 (INFINITY))
>> + return x + x;
>> + /* 0x26.8826ap0 is the largest value such that 10^x < 2^128. */
>> + if (x > 0x26.8826ap0f)
>> + return __math_oflowf (0);
>> + /* -0x2d.278d4p0 is the smallest value such that 10^x > 2^-150. */
>> + if (x < -0x2d.278d4p0f)
>> + return __math_uflowf (0);
>> +#if WANT_ERRNO_UFLOW
>> + if (x < -0x2c.da7cfp0)
>> + return __math_may_uflowf (0);
>> +#endif
>> + /* the smallest value such that 10^x >= 2^-126 (normal range)
>> + is x = -0x25.ee060p0 */
>> + /* we go through here for 2014929 values out of 2060451840
>> + (not counting NaN and infinities, i.e., about 0.1% */
>> + }
>> +
>> + /* x*N*Ln10/Ln2 = k + r with r in [-1/2, 1/2] and int k. */
>> + z = InvLn10N * xd;
>> + /* |xd| < 38 thus |z| < 1216 */
>> +#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;
>> +
>> + /* 10^x = 10^(k/N) * 10^(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;
>> +}
>> +libm_alias_finite (__ieee754_exp10f, __exp10f)
>> diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
>> index bf79274ce7..4817e500e1 100644
>> --- a/sysdeps/ieee754/flt-32/math_config.h
>> +++ b/sysdeps/ieee754/flt-32/math_config.h
>> @@ -109,7 +109,7 @@ 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. */
>> +/* Shared between expf, exp2f, exp10f, and powf. */
>> #define EXP2F_TABLE_BITS 5
>> #define EXP2F_POLY_ORDER 3
>> extern const struct exp2f_data
>>
More information about the Libc-alpha
mailing list