[libc-commits] [libc] 4b41e7b - [libc][math] Universal exp function for cosh/sinh calculation.
Kirill Okhotnikov via libc-commits
libc-commits at lists.llvm.org
Thu Jul 28 02:24:31 PDT 2022
Author: Kirill Okhotnikov
Date: 2022-07-28T11:24:09+02:00
New Revision: 4b41e7b436917895a1d7033bfc4f005053c26643
URL: https://github.com/llvm/llvm-project/commit/4b41e7b436917895a1d7033bfc4f005053c26643
DIFF: https://github.com/llvm/llvm-project/commit/4b41e7b436917895a1d7033bfc4f005053c26643.diff
LOG: [libc][math] Universal exp function for cosh/sinh calculation.
Added a function and test, which can be used later for cosh/sinh
and possibly for expf/expm1f.
Differential Revision: https://reviews.llvm.org/D129215
Added:
libc/src/math/generic/expxf.h
libc/test/src/math/expxf_test.cpp
Modified:
libc/test/src/math/CMakeLists.txt
Removed:
################################################################################
diff --git a/libc/src/math/generic/expxf.h b/libc/src/math/generic/expxf.h
new file mode 100644
index 0000000000000..eefa9031770db
--- /dev/null
+++ b/libc/src/math/generic/expxf.h
@@ -0,0 +1,81 @@
+//===-- 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/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 4ac4e3c32acd4..266e5a41970f7 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1302,6 +1302,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fputil
)
+add_fp_unittest(
+ expxf_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ expxf_test.cpp
+ DEPENDS
+ libc.src.math.generic.common_constants
+ libc.src.__support.FPUtil.fputil
+)
+
add_subdirectory(generic)
add_subdirectory(exhaustive)
add_subdirectory(
diff erential_testing)
diff --git a/libc/test/src/math/expxf_test.cpp b/libc/test/src/math/expxf_test.cpp
new file mode 100644
index 0000000000000..b1d9a211022ab
--- /dev/null
+++ b/libc/test/src/math/expxf_test.cpp
@@ -0,0 +1,38 @@
+//===-- 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);
+ }
+}
More information about the libc-commits
mailing list