summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorRich Felker <dalias@aerifal.cx>2012-04-30 03:26:53 -0400
committerRich Felker <dalias@aerifal.cx>2012-04-30 03:26:53 -0400
commitf6819755779a084bf2f82cb90175a4d9a018de73 (patch)
tree94048c26dd47a1f82507421df1d8815e14675197
parent63374ee233c44891db80e6600d7a5a8c82e4ccca (diff)
downloadmusl-f6819755779a084bf2f82cb90175a4d9a018de73.tar.gz
first try at writing an efficient and "correct" exp10
this is a nonstandard function so it's not clear what conditions it should satisfy. my intent is that it be fast and exact for positive integral exponents when the result fits in the destination type, and fast and correctly rounded for small negative integral exponents. otherwise we aim for at most 1ulp error; it seems to differ from pow by at most 1ulp and it's often 2-5 times faster than pow.
-rw-r--r--include/math.h4
-rw-r--r--src/math/exp10.c19
-rw-r--r--src/math/exp10f.c17
-rw-r--r--src/math/exp10l.c19
4 files changed, 59 insertions, 0 deletions
diff --git a/include/math.h b/include/math.h
index bc0159e6..65b1e915 100644
--- a/include/math.h
+++ b/include/math.h
@@ -387,6 +387,10 @@ float y1f(float);
long double y1l(long double);
float ynf(int, float);
long double ynl(int, long double);
+
+double exp10(double);
+float exp10f(float);
+long double exp10l(long double);
#endif
#ifdef __cplusplus
diff --git a/src/math/exp10.c b/src/math/exp10.c
new file mode 100644
index 00000000..7fd86fba
--- /dev/null
+++ b/src/math/exp10.c
@@ -0,0 +1,19 @@
+#define _GNU_SOURCE
+#include <math.h>
+
+double exp10(double x)
+{
+ static const double p10[] = {
+ 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
+ 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+ 1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
+ 1e10, 1e11, 1e12, 1e13, 1e14, 1e15
+ };
+ double n, y = modf(x, &n);
+ if (fabs(n) < 16) {
+ if (!y) return p10[(int)n+15];
+ y = exp2(3.32192809488736234787031942948939 * y);
+ return y * p10[(int)n+15];
+ }
+ return pow(10.0, x);
+}
diff --git a/src/math/exp10f.c b/src/math/exp10f.c
new file mode 100644
index 00000000..c9521411
--- /dev/null
+++ b/src/math/exp10f.c
@@ -0,0 +1,17 @@
+#define _GNU_SOURCE
+#include <math.h>
+
+float exp10f(float x)
+{
+ static const float p10[] = {
+ 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+ 1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
+ };
+ float n, y = modff(x, &n);
+ if (fabsf(n) < 8) {
+ if (!y) return p10[(int)n+7];
+ y = exp2f(3.32192809488736234787031942948939f * y);
+ return y * p10[(int)n+7];
+ }
+ return exp2(3.32192809488736234787031942948939 * x);
+}
diff --git a/src/math/exp10l.c b/src/math/exp10l.c
new file mode 100644
index 00000000..4d0c5a01
--- /dev/null
+++ b/src/math/exp10l.c
@@ -0,0 +1,19 @@
+#define _GNU_SOURCE
+#include <math.h>
+
+long double exp10l(long double x)
+{
+ static const long double p10[] = {
+ 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
+ 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+ 1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
+ 1e10, 1e11, 1e12, 1e13, 1e14, 1e15
+ };
+ long double n, y = modfl(x, &n);
+ if (fabsl(n) < 16) {
+ if (!y) return p10[(int)n+15];
+ y = exp2l(3.32192809488736234787031942948939L * y);
+ return y * p10[(int)n+15];
+ }
+ return powl(10.0, x);
+}