[libc-commits] [libc] 89ed5b7 - [libc][math] Added auxiliary function log2_eval for asinhf/acoshf/atanhf.
Kirill Okhotnikov via libc-commits
libc-commits at lists.llvm.org
Tue Aug 30 13:40:44 PDT 2022
Author: Kirill Okhotnikov
Date: 2022-08-30T22:39:54+02:00
New Revision: 89ed5b7c50cd37454cddd433a49e19194ae353ce
URL: https://github.com/llvm/llvm-project/commit/89ed5b7c50cd37454cddd433a49e19194ae353ce
DIFF: https://github.com/llvm/llvm-project/commit/89ed5b7c50cd37454cddd433a49e19194ae353ce.diff
LOG: [libc][math] Added auxiliary function log2_eval for asinhf/acoshf/atanhf.
1) `double log2_eval(double)` function added with better than float precision is added.
2) Some refactoring done to put all auxiliary functions and corresponding data
to one place to reuse the code.
3) Added tests for new functions.
4) Performance and precision tests of the function shows, that it more precise than exiting log2,
(no exceptional cases), but timing is ~5% higer that on current one.
Differential Revision: https://reviews.llvm.org/D132809
Added:
libc/src/math/generic/explogxf.cpp
libc/src/math/generic/explogxf.h
libc/test/src/math/explogxf_test.cpp
libc/test/src/math/in_float_range_test_helper.h
Modified:
libc/src/math/generic/CMakeLists.txt
libc/src/math/generic/common_constants.cpp
libc/src/math/generic/common_constants.h
libc/src/math/generic/coshf.cpp
libc/src/math/generic/exp2f.cpp
libc/src/math/generic/sinhf.cpp
libc/src/math/generic/tanhf.cpp
libc/test/src/math/CMakeLists.txt
Removed:
libc/src/math/generic/expxf.h
libc/test/src/math/expxf_test.cpp
################################################################################
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index cf7facb6b5770..ac7a0e6f0a3d3 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -536,7 +536,7 @@ add_entrypoint_object(
HDRS
../exp2f.h
DEPENDS
- .common_constants
+ .explogxf
libc.src.__support.FPUtil.fputil
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
@@ -1166,15 +1166,22 @@ add_entrypoint_object(
-O3
)
+add_object_library(
+ explogxf
+ HDRS
+ explogxf.h
+ SRCS
+ explogxf.cpp
+)
+
add_entrypoint_object(
coshf
SRCS
coshf.cpp
HDRS
../coshf.h
- expxf.h
DEPENDS
- .common_constants
+ .explogxf
libc.src.__support.FPUtil.fputil
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
@@ -1190,9 +1197,8 @@ add_entrypoint_object(
sinhf.cpp
HDRS
../sinhf.h
- expxf.h
DEPENDS
- .common_constants
+ .explogxf
libc.src.__support.FPUtil.fputil
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
@@ -1208,9 +1214,8 @@ add_entrypoint_object(
tanhf.cpp
HDRS
../tanhf.h
- expxf.h
DEPENDS
- .common_constants
+ .explogxf
libc.src.__support.FPUtil.fputil
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp
index 5f7fe56a4d541..d56f620ad3706 100644
--- a/libc/src/math/generic/common_constants.cpp
+++ b/libc/src/math/generic/common_constants.cpp
@@ -225,19 +225,5 @@ const double EXP_M2[128] = {
0x1.4e9c56c731f5dp1, 0x1.513c2e6c731d7p1, 0x1.53e14b042f9cap1,
0x1.568bb722dd593p1, 0x1.593b7d72305bbp1,
};
-// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
-// printf("%.13a,\n", d[i]);
-const double EXP_2_POW[EXP_num_p] = {
- -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2,
- -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3,
- -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3,
- -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4,
- -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
- -0x1.5f134923757f3p-6, 0x0.0000000000000p+0, 0x1.66c34c5615d0fp-6,
- 0x1.6ab0d9f3121ecp-5, 0x1.1301d0125b50ap-4, 0x1.72b83c7d517aep-4,
- 0x1.d4873168b9aa8p-4, 0x1.1c3d373ab11c3p-3, 0x1.4f4efa8fef709p-3,
- 0x1.837f0518db8a9p-3, 0x1.b8d39b9d54e55p-3, 0x1.ef5326091a112p-3,
- 0x1.13821818624b4p-2, 0x1.2ff6b54d8a89cp-2, 0x1.4d0ad5a753e07p-2,
- 0x1.6ac1f752150a5p-2, 0x1.891fac0e95613p-2};
} // namespace __llvm_libc
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index 998d1eb2cd068..4aecff9429df5 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -31,13 +31,6 @@ extern const double EXP_M1[195];
// > for i from 0 to 127 do { D(exp(i / 128)); };
extern const double EXP_M2[128];
-static constexpr int EXP_bits_p = 5;
-static constexpr int EXP_num_p = 1 << EXP_bits_p;
-
-// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
-// printf("%.13a,\n", d[i]);
-extern const double EXP_2_POW[EXP_num_p];
-
} // namespace __llvm_libc
#endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H
diff --git a/libc/src/math/generic/coshf.cpp b/libc/src/math/generic/coshf.cpp
index ce1f8085adeb9..95f7f7775cf62 100644
--- a/libc/src/math/generic/coshf.cpp
+++ b/libc/src/math/generic/coshf.cpp
@@ -9,7 +9,7 @@
#include "src/math/coshf.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/multiply_add.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
namespace __llvm_libc {
diff --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp
index 9756b1ac51901..5d6b6bea1f22b 100644
--- a/libc/src/math/generic/exp2f.cpp
+++ b/libc/src/math/generic/exp2f.cpp
@@ -8,6 +8,7 @@
#include "src/math/exp2f.h"
#include "common_constants.h"
+#include "explogxf.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -19,9 +20,6 @@
namespace __llvm_libc {
-constexpr float mlp = EXP_num_p;
-constexpr float mmld = -1.0 / mlp;
-
constexpr uint32_t exval1 = 0x3b42'9d37U;
constexpr uint32_t exval2 = 0xbcf3'a937U;
constexpr uint32_t exval_mask = exval1 & exval2;
@@ -78,24 +76,7 @@ LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
}
}
- float kf = fputil::nearest_integer(x * mlp);
- double dx = fputil::multiply_add(mmld, kf, x);
- double mult_f, ml;
- {
- uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
- (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
- fputil::FPBits<double> bs;
- bs.set_unbiased_exponent(ps >> EXP_bits_p);
- ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
- mult_f = bs.get_val();
- }
-
- // N[Table[Ln[2]^n/n!,{n,1,6}],30]
- double pe = fputil::polyeval(
- dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
- 0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
-
- return mult_f * ml * pe;
+ return exp2_eval(x);
}
} // namespace __llvm_libc
diff --git a/libc/src/math/generic/explogxf.cpp b/libc/src/math/generic/explogxf.cpp
new file mode 100644
index 0000000000000..47d54e70b04e4
--- /dev/null
+++ b/libc/src/math/generic/explogxf.cpp
@@ -0,0 +1,89 @@
+//===-- Single-precision general exp/log functions ------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "explogxf.h"
+
+namespace __llvm_libc {
+
+// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
+// printf("%.13a,\n", d[i]);
+alignas(64) const double EXP_2_POW[EXP_num_p] = {
+ -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2,
+ -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3,
+ -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3,
+ -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4,
+ -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
+ -0x1.5f134923757f3p-6, 0x0.0000000000000p+0, 0x1.66c34c5615d0fp-6,
+ 0x1.6ab0d9f3121ecp-5, 0x1.1301d0125b50ap-4, 0x1.72b83c7d517aep-4,
+ 0x1.d4873168b9aa8p-4, 0x1.1c3d373ab11c3p-3, 0x1.4f4efa8fef709p-3,
+ 0x1.837f0518db8a9p-3, 0x1.b8d39b9d54e55p-3, 0x1.ef5326091a112p-3,
+ 0x1.13821818624b4p-2, 0x1.2ff6b54d8a89cp-2, 0x1.4d0ad5a753e07p-2,
+ 0x1.6ac1f752150a5p-2, 0x1.891fac0e95613p-2};
+
+// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
+alignas(64) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
+ 0x0.0000000000000p+0, 0x1.6e79685c2d22ap-6, 0x1.6bad3758efd87p-5,
+ 0x1.0eb389fa29f9bp-4, 0x1.663f6fac91316p-4, 0x1.bc84240adabbap-4,
+ 0x1.08c588cda79e4p-3, 0x1.32ae9e278ae1ap-3, 0x1.5c01a39fbd688p-3,
+ 0x1.84c2bd02f03b3p-3, 0x1.acf5e2db4ec94p-3, 0x1.d49ee4c325970p-3,
+ 0x1.fbc16b902680ap-3, 0x1.11307dad30b76p-2, 0x1.24407ab0e073ap-2,
+ 0x1.37124cea4cdedp-2, 0x1.49a784bcd1b8bp-2, 0x1.5c01a39fbd688p-2,
+ 0x1.6e221cd9d0cdep-2, 0x1.800a563161c54p-2, 0x1.91bba891f1709p-2,
+ 0x1.a33760a7f6051p-2, 0x1.b47ebf73882a1p-2, 0x1.c592fad295b56p-2,
+ 0x1.d6753e032ea0fp-2, 0x1.e726aa1e754d2p-2, 0x1.f7a8568cb06cfp-2,
+ 0x1.03fda8b97997fp-1, 0x1.0c10500d63aa6p-1, 0x1.140c9faa1e544p-1,
+ 0x1.1bf311e95d00ep-1, 0x1.23c41d42727c8p-1, 0x1.2b803473f7ad1p-1,
+ 0x1.3327c6ab49ca7p-1, 0x1.3abb3faa02167p-1, 0x1.423b07e986aa9p-1,
+ 0x1.49a784bcd1b8bp-1, 0x1.510118708a8f9p-1, 0x1.5848226989d34p-1,
+ 0x1.5f7cff41e09afp-1, 0x1.66a008e4788ccp-1, 0x1.6db196a76194ap-1,
+ 0x1.74b1fd64e0754p-1, 0x1.7ba18f93502e4p-1, 0x1.82809d5be7073p-1,
+ 0x1.894f74b06ef8bp-1, 0x1.900e6160002cdp-1, 0x1.96bdad2acb5f6p-1,
+ 0x1.9d5d9fd5010b3p-1, 0x1.a3ee7f38e181fp-1, 0x1.aa708f58014d3p-1,
+ 0x1.b0e4126bcc86cp-1, 0x1.b74948f5532dap-1, 0x1.bda071cc67e6ep-1,
+ 0x1.c3e9ca2e1a055p-1, 0x1.ca258dca93316p-1, 0x1.d053f6d260896p-1,
+ 0x1.d6753e032ea0fp-1, 0x1.dc899ab3ff56cp-1, 0x1.e29142e0e0140p-1,
+ 0x1.e88c6b3626a73p-1, 0x1.ee7b471b3a950p-1, 0x1.f45e08bcf0655p-1,
+ 0x1.fa34e1177c233p-1,
+};
+
+// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
+alignas(64) const double LOG_P1_1_OVER[LOG_P1_SIZE] = {
+ 0x1.0000000000000p+0, 0x1.f81f81f81f820p-1, 0x1.f07c1f07c1f08p-1,
+ 0x1.e9131abf0b767p-1, 0x1.e1e1e1e1e1e1ep-1, 0x1.dae6076b981dbp-1,
+ 0x1.d41d41d41d41dp-1, 0x1.cd85689039b0bp-1, 0x1.c71c71c71c71cp-1,
+ 0x1.c0e070381c0e0p-1, 0x1.bacf914c1bad0p-1, 0x1.b4e81b4e81b4fp-1,
+ 0x1.af286bca1af28p-1, 0x1.a98ef606a63bep-1, 0x1.a41a41a41a41ap-1,
+ 0x1.9ec8e951033d9p-1, 0x1.999999999999ap-1, 0x1.948b0fcd6e9e0p-1,
+ 0x1.8f9c18f9c18fap-1, 0x1.8acb90f6bf3aap-1, 0x1.8618618618618p-1,
+ 0x1.8181818181818p-1, 0x1.7d05f417d05f4p-1, 0x1.78a4c8178a4c8p-1,
+ 0x1.745d1745d1746p-1, 0x1.702e05c0b8170p-1, 0x1.6c16c16c16c17p-1,
+ 0x1.6816816816817p-1, 0x1.642c8590b2164p-1, 0x1.6058160581606p-1,
+ 0x1.5c9882b931057p-1, 0x1.58ed2308158edp-1, 0x1.5555555555555p-1,
+ 0x1.51d07eae2f815p-1, 0x1.4e5e0a72f0539p-1, 0x1.4afd6a052bf5bp-1,
+ 0x1.47ae147ae147bp-1, 0x1.446f86562d9fbp-1, 0x1.4141414141414p-1,
+ 0x1.3e22cbce4a902p-1, 0x1.3b13b13b13b14p-1, 0x1.3813813813814p-1,
+ 0x1.3521cfb2b78c1p-1, 0x1.323e34a2b10bfp-1, 0x1.2f684bda12f68p-1,
+ 0x1.2c9fb4d812ca0p-1, 0x1.29e4129e4129ep-1, 0x1.27350b8812735p-1,
+ 0x1.2492492492492p-1, 0x1.21fb78121fb78p-1, 0x1.1f7047dc11f70p-1,
+ 0x1.1cf06ada2811dp-1, 0x1.1a7b9611a7b96p-1, 0x1.1811811811812p-1,
+ 0x1.15b1e5f75270dp-1, 0x1.135c81135c811p-1, 0x1.1111111111111p-1,
+ 0x1.0ecf56be69c90p-1, 0x1.0c9714fbcda3bp-1, 0x1.0a6810a6810a7p-1,
+ 0x1.0842108421084p-1, 0x1.0624dd2f1a9fcp-1, 0x1.0410410410410p-1,
+ 0x1.0204081020408p-1};
+
+// Taylos series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
+// K_LOG2_ODD starts from x^3
+alignas(64) const
+ double K_LOG2_ODD[4] = {0x1.ec709dc3a03fdp-2, 0x1.2776c50ef9bfep-2,
+ 0x1.a61762a7aded9p-3, 0x1.484b13d7c02a9p-3};
+
+alignas(64) const
+ double K_LOG2_EVEN[4] = {-0x1.71547652b82fep-1, -0x1.71547652b82fep-2,
+ -0x1.ec709dc3a03fdp-3, -0x1.2776c50ef9bfep-3};
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h
new file mode 100644
index 0000000000000..5f79f57fdfa7d
--- /dev/null
+++ b/libc/src/math/generic/explogxf.h
@@ -0,0 +1,159 @@
+//===-- Single-precision general exp/log functions ------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
+
+#include "common_constants.h" // Lookup tables EXP_M
+#include "math_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include <src/__support/FPUtil/NearestIntegerOperations.h>
+
+#include <errno.h>
+
+namespace __llvm_libc {
+
+static constexpr int EXP_bits_p = 5;
+static constexpr int EXP_num_p = 1 << EXP_bits_p;
+constexpr double mlp = EXP_num_p;
+constexpr double mmld = -1.0 / mlp;
+
+// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
+// printf("%.13a,\n", d[i]);
+extern const double EXP_2_POW[EXP_num_p];
+
+constexpr int LOG_P1_BITS = 6;
+constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS;
+
+// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_LOG2[LOG_P1_SIZE];
+
+// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_1_OVER[LOG_P1_SIZE];
+
+// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
+// K_LOG2_ODD starts from x^3
+extern const double K_LOG2_ODD[4];
+extern const double K_LOG2_EVEN[4];
+
+// The algorithm represents exp(x) as
+// exp(x) = 2^(ln(2) * i) * 2^(ln(2) * j / NUM_P )) * exp(dx)
+// where i integer value, j integer in range [-NUM_P/2, NUM_P/2).
+// 2^(ln(2) * j / NUM_P )) is a table values: 1.0 + EXP_M
+// exp(dx) calculates by taylor expansion.
+
+// Inversion of ln(2). Multiplication by EXP_num_p due to sampling by 1 /
+// EXP_num_p Precise value of the constant is not needed.
+static constexpr double LN2_INV = 0x1.71547652b82fep+0 * EXP_num_p;
+
+// LN2_HIGH + LN2_LOW = ln(2) with precision higher than double(ln(2))
+// Minus sign is to use FMA directly.
+static constexpr double LN2_HIGH = -0x1.62e42fefa0000p-1 / EXP_num_p;
+static constexpr double LN2_LOW = -0x1.cf79abc9e3b3ap-40 / EXP_num_p;
+
+struct exe_eval_result_t {
+ // exp(x) = 2^MULT_POWER2 * mult_exp * (r + 1.0)
+ // where
+ // MULT_POWER2 template parameter;
+ // mult_exp = 2^e;
+ // r in range [~-0.3, ~0.41]
+ double mult_exp;
+ double r;
+};
+
+// The function correctly calculates exp value with at least float precision
+// in range not narrow than [-log(2^-150), 90]
+template <int MULT_POWER2 = 0>
+inline static exe_eval_result_t exp_eval(double x) {
+ double ps_dbl = fputil::nearest_integer(LN2_INV * x);
+ // Negative sign due to multiply_add optimization
+ double mult_e1, ml;
+ {
+ int ps =
+ static_cast<int>(ps_dbl) + (1 << (EXP_bits_p - 1)) +
+ ((fputil::FPBits<double>::EXPONENT_BIAS + MULT_POWER2) << EXP_bits_p);
+ int table_index = ps & (EXP_num_p - 1);
+ fputil::FPBits<double> bs;
+ bs.set_unbiased_exponent(ps >> EXP_bits_p);
+ ml = EXP_2_POW[table_index];
+ mult_e1 = bs.get_val();
+ }
+ double dx = fputil::multiply_add(ps_dbl, LN2_LOW,
+ fputil::multiply_add(ps_dbl, LN2_HIGH, x));
+
+ // Taylor series coefficients
+ double pe = dx * fputil::polyeval(dx, 1.0, 0x1.0p-1, 0x1.5555555555555p-3,
+ 0x1.5555555555555p-5, 0x1.1111111111111p-7,
+ 0x1.6c16c16c16c17p-10);
+
+ double r = fputil::multiply_add(ml, pe, pe) + ml;
+ return {mult_e1, r};
+}
+
+inline static double exp2_eval(double x) {
+ double kf = fputil::nearest_integer(x * mlp);
+ double dx = fputil::multiply_add(mmld, kf, x);
+ double mult_f, ml;
+ {
+ uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
+ (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
+ fputil::FPBits<double> bs;
+ bs.set_unbiased_exponent(ps >> EXP_bits_p);
+ ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
+ mult_f = bs.get_val();
+ }
+
+ // Taylor series coefficients for 2^x
+ double pe = fputil::polyeval(
+ dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
+ 0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
+
+ return mult_f * ml * pe;
+}
+
+// x should be positive, normal finite value
+inline static double log2_eval(double x) {
+ using FPB = fputil::FPBits<double>;
+ FPB bs(x);
+
+ double result = 0;
+ result += bs.get_exponent();
+
+ int p1 =
+ (bs.get_mantissa() >> (FPB::FloatProp::MANTISSA_WIDTH - LOG_P1_BITS)) &
+ (LOG_P1_SIZE - 1);
+
+ bs.bits &= FPB::FloatProp::MANTISSA_MASK >> LOG_P1_BITS;
+ bs.set_unbiased_exponent(FPB::FloatProp::EXPONENT_BIAS);
+ double dx = (bs.get_val() - 1.0) * LOG_P1_1_OVER[p1];
+
+ // Taylor series for log(2,1+x)
+ double c1 = fputil::multiply_add(dx, K_LOG2_ODD[0], K_LOG2_EVEN[0]);
+ double c2 = fputil::multiply_add(dx, K_LOG2_ODD[1], K_LOG2_EVEN[1]);
+ double c3 = fputil::multiply_add(dx, K_LOG2_ODD[2], K_LOG2_EVEN[2]);
+ double c4 = fputil::multiply_add(dx, K_LOG2_ODD[3], K_LOG2_EVEN[3]);
+
+ // c0 = dx * (1.0 / ln(2)) + LOG_P1_LOG2[p1]
+ double c0 = fputil::multiply_add(dx, 0x1.71547652b82fep+0, LOG_P1_LOG2[p1]);
+ result += __llvm_libc::fputil::polyeval(dx * dx, c0, c1, c2, c3, c4);
+ return result;
+}
+
+// x should be positive, normal finite value
+inline static double log_eval(double x) {
+ // ln(x) = log[2,x] * ln(2)
+ return log2_eval(x) * 0x1.62e42fefa39efp-1;
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
diff --git a/libc/src/math/generic/expxf.h b/libc/src/math/generic/expxf.h
deleted file mode 100644
index eefa9031770db..0000000000000
--- a/libc/src/math/generic/expxf.h
+++ /dev/null
@@ -1,81 +0,0 @@
-//===-- Single-precision general exp function -----------------------------===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-//
-//===----------------------------------------------------------------------===//
-
-#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
-#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
-
-#include "common_constants.h" // Lookup tables EXP_M
-#include "math_utils.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/nearest_integer.h"
-#include "src/__support/common.h"
-#include <src/__support/FPUtil/NearestIntegerOperations.h>
-
-#include <errno.h>
-
-namespace __llvm_libc {
-
-// The algorithm represents exp(x) as
-// exp(x) = 2^(ln(2) * i) * 2^(ln(2) * j / NUM_P )) * exp(dx)
-// where i integer value, j integer in range [-NUM_P/2, NUM_P/2).
-// 2^(ln(2) * j / NUM_P )) is a table values: 1.0 + EXP_M
-// exp(dx) calculates by taylor expansion.
-
-// Inversion of ln(2). Multiplication by EXP_num_p due to sampling by 1 /
-// EXP_num_p Precise value of the constant is not needed.
-static constexpr double LN2_INV = 0x1.71547652b82fep+0 * EXP_num_p;
-
-// LN2_HIGH + LN2_LOW = ln(2) with precision higher than double(ln(2))
-// Minus sign is to use FMA directly.
-static constexpr double LN2_HIGH = -0x1.62e42fefa0000p-1 / EXP_num_p;
-static constexpr double LN2_LOW = -0x1.cf79abc9e3b3ap-40 / EXP_num_p;
-
-struct exe_eval_result_t {
- // exp(x) = 2^MULT_POWER2 * mult_exp * (r + 1.0)
- // where
- // MULT_POWER2 template parameter;
- // mult_exp = 2^e;
- // r in range [~-0.3, ~0.41]
- double mult_exp;
- double r;
-};
-
-// The function correctly calculates exp value with at least float precision
-// in range not narrow than [-log(2^-150), 90]
-template <int MULT_POWER2 = 0>
-inline static exe_eval_result_t exp_eval(double x) {
- double ps_dbl = fputil::nearest_integer(LN2_INV * x);
- // Negative sign due to multiply_add optimization
- double mult_e1, ml;
- {
- int ps =
- static_cast<int>(ps_dbl) + (1 << (EXP_bits_p - 1)) +
- ((fputil::FPBits<double>::EXPONENT_BIAS + MULT_POWER2) << EXP_bits_p);
- int table_index = ps & (EXP_num_p - 1);
- fputil::FPBits<double> bs;
- bs.set_unbiased_exponent(ps >> EXP_bits_p);
- ml = EXP_2_POW[table_index];
- mult_e1 = bs.get_val();
- }
- double dx = fputil::multiply_add(ps_dbl, LN2_LOW,
- fputil::multiply_add(ps_dbl, LN2_HIGH, x));
-
- // Taylor series coefficients
- double pe = dx * fputil::polyeval(dx, 1.0, 0x1.0p-1, 0x1.5555555555555p-3,
- 0x1.5555555555555p-5, 0x1.1111111111111p-7,
- 0x1.6c16c16c16c17p-10);
-
- double r = fputil::multiply_add(ml, pe, pe) + ml;
- return {mult_e1, r};
-}
-
-} // namespace __llvm_libc
-
-#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp
index 7c7e6ae3e9cc2..7b5ac60454a3b 100644
--- a/libc/src/math/generic/sinhf.cpp
+++ b/libc/src/math/generic/sinhf.cpp
@@ -8,7 +8,7 @@
#include "src/math/sinhf.h"
#include "src/__support/FPUtil/FPBits.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
namespace __llvm_libc {
diff --git a/libc/src/math/generic/tanhf.cpp b/libc/src/math/generic/tanhf.cpp
index fd374b4c016d9..e1f753c12161a 100644
--- a/libc/src/math/generic/tanhf.cpp
+++ b/libc/src/math/generic/tanhf.cpp
@@ -8,7 +8,7 @@
#include "src/math/tanhf.h"
#include "src/__support/FPUtil/FPBits.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
namespace __llvm_libc {
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index e51ae24428ae7..8d43f16513eae 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1319,14 +1319,16 @@ add_fp_unittest(
)
add_fp_unittest(
- expxf_test
+ explogxf_test
NEED_MPFR
SUITE
libc_math_unittests
+ HDRS
+ in_float_range_test_helper.h
SRCS
- expxf_test.cpp
+ explogxf_test.cpp
DEPENDS
- libc.src.math.generic.common_constants
+ libc.src.math.generic.explogxf
libc.src.__support.FPUtil.fputil
)
diff --git a/libc/test/src/math/explogxf_test.cpp b/libc/test/src/math/explogxf_test.cpp
new file mode 100644
index 0000000000000..acaaacacc98c3
--- /dev/null
+++ b/libc/test/src/math/explogxf_test.cpp
@@ -0,0 +1,55 @@
+//===-- Unittests for explogxf --------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "in_float_range_test_helper.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/explogxf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+constexpr int def_count = 100003;
+constexpr float def_prec = 0.500001f;
+
+TEST(LlvmLibcExpxfTest, InFloatRange) {
+ auto fx = [](float x) -> float {
+ auto result = __llvm_libc::exp_eval<-1>(x);
+ return static_cast<float>(2 * result.mult_exp * result.r +
+ 2 * result.mult_exp);
+ };
+ auto f_check = [](float x) -> bool {
+ return !(
+ (isnan(x) || isinf(x) || x < -70 || x > 70 || fabsf(x) < 0x1.0p-10));
+ };
+ CHECK_DATA(0.0f, neg_inf, mpfr::Operation::Exp, fx, f_check, def_count,
+ def_prec);
+}
+
+TEST(LlvmLibcExp2xfTest, InFloatRange) {
+ auto f_check = [](float x) -> bool {
+ return !(
+ (isnan(x) || isinf(x) || x < -130 || x > 130 || fabsf(x) < 0x1.0p-10));
+ };
+ CHECK_DATA(0.0f, neg_inf, mpfr::Operation::Exp2, __llvm_libc::exp2_eval,
+ f_check, def_count, def_prec);
+}
+
+TEST(LlvmLibcLog2xfTest, InFloatRange) {
+ CHECK_DATA(0.0f, inf, mpfr::Operation::Log2, __llvm_libc::log2_eval, isnormal,
+ def_count, def_prec);
+}
+
+TEST(LlvmLibcLogxfTest, InFloatRange) {
+ CHECK_DATA(0.0f, inf, mpfr::Operation::Log, __llvm_libc::log_eval, isnormal,
+ def_count, def_prec);
+}
diff --git a/libc/test/src/math/expxf_test.cpp b/libc/test/src/math/expxf_test.cpp
deleted file mode 100644
index b1d9a211022ab..0000000000000
--- a/libc/test/src/math/expxf_test.cpp
+++ /dev/null
@@ -1,38 +0,0 @@
-//===-- Unittests for expxf -----------------------------------------------===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-//
-//===----------------------------------------------------------------------===//
-
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/math/generic/expxf.h"
-#include "utils/MPFRWrapper/MPFRUtils.h"
-#include "utils/UnitTest/FPMatcher.h"
-#include "utils/UnitTest/Test.h"
-#include <math.h>
-
-#include <errno.h>
-#include <stdint.h>
-
-namespace mpfr = __llvm_libc::testing::mpfr;
-
-DECLARE_SPECIAL_CONSTANTS(float)
-
-TEST(LlvmLibcExpxfTest, InFloatRange) {
- constexpr uint32_t COUNT = 1000000;
- constexpr uint32_t STEP = UINT32_MAX / COUNT;
- auto fx = [](float x) -> float {
- auto result = __llvm_libc::exp_eval<-1>(x);
- return static_cast<float>(2 * result.mult_exp * result.r +
- 2 * result.mult_exp);
- };
- for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
- float x = float(FPBits(v));
- if (isnan(x) || isinf(x) || x < -70 || x > 70 || fabsf(x) < 0x1.0p-10)
- continue;
-
- EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, fx(x), 0.5);
- }
-}
diff --git a/libc/test/src/math/in_float_range_test_helper.h b/libc/test/src/math/in_float_range_test_helper.h
new file mode 100644
index 0000000000000..5f345c0cf17af
--- /dev/null
+++ b/libc/test/src/math/in_float_range_test_helper.h
@@ -0,0 +1,26 @@
+//
+// Created by kirill on 8/30/22.
+//
+
+#ifndef LLVM_IN_FLOAT_RANGE_TEST_HELPER_H
+#define LLVM_IN_FLOAT_RANGE_TEST_HELPER_H
+
+#include <stdint.h>
+
+#define CHECK_DATA(start, stop, mfp_op, f, f_check, count, prec) \
+ { \
+ uint64_t ustart = FPBits(start).uintval(); \
+ uint64_t ustop = FPBits(stop).uintval(); \
+ for (uint64_t i = 0;; ++i) { \
+ uint64_t v = ustart + (ustop - ustart) * i / count; \
+ if (v > ustop) \
+ break; \
+ float x = FPBits(uint32_t(v)).get_val(); \
+ if ((f_check)(x)) { \
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mfp_op, x, static_cast<float>((f)(x)), \
+ (prec)); \
+ } \
+ } \
+ }
+
+#endif // LLVM_IN_FLOAT_RANGE_TEST_HELPER_H
More information about the libc-commits
mailing list