* [PATCH 2/5] math: Optimized generic exp10f with wrappers
@ 2020-06-18 16:44 Wilco Dijkstra
0 siblings, 0 replies; 6+ messages in thread
From: Wilco Dijkstra @ 2020-06-18 16:44 UTC (permalink / raw)
To: Adhemerval Zanella; +Cc: 'GNU C Library', paul zimmermann
Hi,
This looks good (I like the extra comments with the bounds), OK with adding
the missing 'f' on the float constant below.
Wilco
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)
OK
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)
This is missing an 'f' on the float constant.
+ 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)
OK
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
OK
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 2/5] math: Optimized generic exp10f with wrappers
2020-04-29 17:11 ` Adhemerval Zanella
@ 2020-05-20 20:18 ` Adhemerval Zanella
0 siblings, 0 replies; 6+ messages in thread
From: Adhemerval Zanella @ 2020-05-20 20:18 UTC (permalink / raw)
To: libc-alpha
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
>>
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 2/5] math: Optimized generic exp10f with wrappers
2020-04-09 19:59 ` [PATCH 2/5] math: Optimized generic exp10f with wrappers Adhemerval Zanella
2020-04-10 6:32 ` paul zimmermann
@ 2020-04-29 17:11 ` Adhemerval Zanella
2020-05-20 20:18 ` Adhemerval Zanella
1 sibling, 1 reply; 6+ messages in thread
From: Adhemerval Zanella @ 2020-04-29 17:11 UTC (permalink / raw)
To: libc-alpha
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
>
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 2/5] math: Optimized generic exp10f with wrappers
2020-04-10 6:32 ` paul zimmermann
@ 2020-04-16 20:42 ` Adhemerval Zanella
0 siblings, 0 replies; 6+ messages in thread
From: Adhemerval Zanella @ 2020-04-16 20:42 UTC (permalink / raw)
To: paul zimmermann; +Cc: libc-alpha
On 10/04/2020 03:32, paul zimmermann wrote:
> Dear Adhemerval,
>
>> 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
>> }
>
> the timings for x86_64 and aarch64 are identical.
> A copy/paste error?
>
Indeed, I have updated all timing based on new benchmark version as well.
The new results are:
Result for x86_64 (i7-4790K CPU @ 4.00GHz) are:
Before new code:
"exp10f": {
"workload-spec2017.wrf (adapted)": {
"duration": 4.0414e+09,
"iterations": 1.00128e+08,
"reciprocal-throughput": 26.6818,
"latency": 54.043,
"max-throughput": 3.74787e+07,
"min-throughput": 1.85038e+07
}
With new code:
"exp10f": {
"workload-spec2017.wrf (adapted)": {
"duration": 4.11951e+09,
"iterations": 1.23968e+08,
"reciprocal-throughput": 21.0581,
"latency": 45.4028,
"max-throughput": 4.74876e+07,
"min-throughput": 2.20251e+07
}
Result for aarch64 (A72 @ 2GHz) are:
Before new code:
"exp10f": {
"workload-spec2017.wrf (adapted)": {
"duration": 4.62362e+09,
"iterations": 3.3376e+07,
"reciprocal-throughput": 127.698,
"latency": 149.365,
"max-throughput": 7.831e+06,
"min-throughput": 6.69501e+06
}
With new code:
"exp10f": {
"workload-spec2017.wrf (adapted)": {
"duration": 4.29108e+09,
"iterations": 6.6752e+07,
"reciprocal-throughput": 51.2111,
"latency": 77.3568,
"max-throughput": 1.9527e+07,
"min-throughput": 1.29271e+07
}
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 2/5] math: Optimized generic exp10f with wrappers
2020-04-09 19:59 ` [PATCH 2/5] math: Optimized generic exp10f with wrappers Adhemerval Zanella
@ 2020-04-10 6:32 ` paul zimmermann
2020-04-16 20:42 ` Adhemerval Zanella
2020-04-29 17:11 ` Adhemerval Zanella
1 sibling, 1 reply; 6+ messages in thread
From: paul zimmermann @ 2020-04-10 6:32 UTC (permalink / raw)
To: Adhemerval Zanella; +Cc: libc-alpha
Dear Adhemerval,
> 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
> }
the timings for x86_64 and aarch64 are identical.
A copy/paste error?
Paul
^ permalink raw reply [flat|nested] 6+ messages in thread
* [PATCH 2/5] math: Optimized generic exp10f with wrappers
2020-04-09 19:59 [PATCH 1/5] benchtests: Add exp10f benchmark Adhemerval Zanella
@ 2020-04-09 19:59 ` Adhemerval Zanella
2020-04-10 6:32 ` paul zimmermann
2020-04-29 17:11 ` Adhemerval Zanella
0 siblings, 2 replies; 6+ messages in thread
From: Adhemerval Zanella @ 2020-04-09 19:59 UTC (permalink / raw)
To: libc-alpha
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
--
2.17.1
^ permalink raw reply [flat|nested] 6+ messages in thread
end of thread, other threads:[~2020-06-18 16:44 UTC | newest]
Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2020-06-18 16:44 [PATCH 2/5] math: Optimized generic exp10f with wrappers Wilco Dijkstra
-- strict thread matches above, loose matches on Subject: below --
2020-04-09 19:59 [PATCH 1/5] benchtests: Add exp10f benchmark Adhemerval Zanella
2020-04-09 19:59 ` [PATCH 2/5] math: Optimized generic exp10f with wrappers Adhemerval Zanella
2020-04-10 6:32 ` paul zimmermann
2020-04-16 20:42 ` Adhemerval Zanella
2020-04-29 17:11 ` Adhemerval Zanella
2020-05-20 20:18 ` Adhemerval Zanella
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).