[PATCH 2/5] math: Optimized generic exp10f with wrappers
Wilco Dijkstra
Wilco.Dijkstra@arm.com
Thu Jun 18 16:44:47 GMT 2020
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/flt32/e_exp10f.c  198 +++++++++++++++++++++++++++
sysdeps/ieee754/flt32/math_config.h  2 +
3 files changed, 199 insertions(+), 33 deletions()
delete mode 100644 math/e_exp10f.c
create mode 100644 sysdeps/ieee754/flt32/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) 19982020 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 <libmaliasfinite.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/flt32/e_exp10f.c b/sysdeps/ieee754/flt32/e_exp10f.c
new file mode 100644
index 0000000000..034a9e364a
 /dev/null
+++ b/sysdeps/ieee754/flt32/e_exp10f.c
@@ 0,0 +1,198 @@
+/* Singleprecision 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 <mathnarroweval.h>
+#include <stdint.h>
+#include <libmaliasfinite.h>
+#include <libmaliasfloat.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.
+ Nonnearest 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 degree3 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.694e6, C[0] * r < 8.47e7, and the absolute error on
+ C[0] * r is bounded by 1/2*ulp(8.47e7) = 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/flt32/math_config.h b/sysdeps/ieee754/flt32/math_config.h
index bf79274ce7..4817e500e1 100644
 a/sysdeps/ieee754/flt32/math_config.h
+++ b/sysdeps/ieee754/flt32/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
More information about the Libcalpha
mailing list