[libc-commits] [libc] 4e5f8b4 - [libc] Add implementation of expm1f.
Tue Ly via libc-commits
libc-commits at lists.llvm.org
Thu Jun 10 12:04:58 PDT 2021
Author: Tue Ly
Date: 2021-06-10T14:58:34-04:00
New Revision: 4e5f8b4d8d9d7a6039e10b9507dac896eed92040
URL: https://github.com/llvm/llvm-project/commit/4e5f8b4d8d9d7a6039e10b9507dac896eed92040
DIFF: https://github.com/llvm/llvm-project/commit/4e5f8b4d8d9d7a6039e10b9507dac896eed92040.diff
LOG: [libc] Add implementation of expm1f.
Use expm1f(x) = exp(x) - 1 for |x| > ln(2).
For |x| <= ln(2), divide it into 3 subintervals: [-ln2, -1/8], [-1/8, 1/8], [1/8, ln2]
and use a degree-6 polynomial approximation generated by Sollya's fpminmax for each interval.
Errors < 1.5 ULPs when we use fma to evaluate the polynomials.
Differential Revision: https://reviews.llvm.org/D101134
Added:
libc/src/math/expm1f.h
libc/src/math/generic/expm1f.cpp
libc/test/src/math/differential_testing/expm1f_diff.cpp
libc/test/src/math/differential_testing/expm1f_perf.cpp
libc/test/src/math/exhaustive/expm1f_test.cpp
libc/test/src/math/expm1f_test.cpp
libc/utils/FPUtil/PolyEval.h
libc/utils/mathtools/expm1f.sollya
Modified:
libc/config/linux/aarch64/entrypoints.txt
libc/config/linux/x86_64/entrypoints.txt
libc/spec/stdc.td
libc/src/math/CMakeLists.txt
libc/src/math/generic/CMakeLists.txt
libc/test/src/math/CMakeLists.txt
libc/test/src/math/differential_testing/CMakeLists.txt
libc/test/src/math/exhaustive/CMakeLists.txt
libc/utils/FPUtil/BitPatterns.h
libc/utils/FPUtil/CMakeLists.txt
libc/utils/FPUtil/generic/FMA.h
libc/utils/MPFRWrapper/MPFRUtils.cpp
libc/utils/MPFRWrapper/MPFRUtils.h
Removed:
################################################################################
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 211d918bc54bd..95e8b8c637578 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -68,6 +68,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.cosf
libc.src.math.expf
libc.src.math.exp2f
+ libc.src.math.expm1f
libc.src.math.fabs
libc.src.math.fabsf
libc.src.math.fabsl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 331ea16299d89..5143864e013bf 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -69,6 +69,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.cosf
libc.src.math.expf
libc.src.math.exp2f
+ libc.src.math.expm1f
libc.src.math.fabs
libc.src.math.fabsf
libc.src.math.fabsl
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index c2a90450dbd2a..c6621812258cc 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -396,6 +396,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+ FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"remainderf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
FunctionSpec<"remainder", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index e2abe21a213bf..ddee3ca6c3104 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -79,6 +79,8 @@ add_math_entrypoint_object(expf)
add_math_entrypoint_object(exp2f)
+add_math_entrypoint_object(expm1f)
+
add_math_entrypoint_object(fabs)
add_math_entrypoint_object(fabsf)
add_math_entrypoint_object(fabsl)
diff --git a/libc/src/math/expm1f.h b/libc/src/math/expm1f.h
new file mode 100644
index 0000000000000..2577e4a492980
--- /dev/null
+++ b/libc/src/math/expm1f.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for expm1f ------------------------*- C++ -*-===//
+//
+// 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_EXPM1F_H
+#define LLVM_LIBC_SRC_MATH_EXPM1F_H
+
+namespace __llvm_libc {
+
+float expm1f(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_EXPM1F_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 45dcc5ae27b31..686d2757842c4 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -489,6 +489,18 @@ add_entrypoint_object(
libc.include.math
)
+add_entrypoint_object(
+ expm1f
+ SRCS
+ expm1f.cpp
+ HDRS
+ ../expm1f.h
+ DEPENDS
+ libc.include.math
+ libc.src.math.expf
+ libc.src.math.fabsf
+)
+
add_entrypoint_object(
copysign
SRCS
diff --git a/libc/src/math/generic/expm1f.cpp b/libc/src/math/generic/expm1f.cpp
new file mode 100644
index 0000000000000..61ebbf288c53a
--- /dev/null
+++ b/libc/src/math/generic/expm1f.cpp
@@ -0,0 +1,57 @@
+//===-- Implementation of expm1f 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/expm1f.h"
+#include "src/__support/common.h"
+#include "src/math/expf.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/PolyEval.h"
+
+namespace __llvm_libc {
+
+// When |x| > Ln2, catastrophic cancellation does not occur with the
+// subtraction expf(x) - 1.0f, so we use it to compute expm1f(x).
+//
+// We divide [-Ln2; Ln2] into 3 subintervals [-Ln2; -1/8], [-1/8; 1/8],
+// [1/8; Ln2]. And we use a degree-6 polynomial to approximate exp(x) - 1 in
+// each interval. The coefficients were generated by Sollya's fpminmax.
+//
+// See libc/utils/mathtools/expm1f.sollya for more detail.
+LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
+ const float Ln2 =
+ 0.69314718055994530941723212145817656807550013436025f; // For C++17:
+ // 0x1.62e'42ffp-1
+ float abs_x = __llvm_libc::fputil::abs(x);
+
+ if (abs_x <= Ln2) {
+ if (abs_x <= 0.125f) {
+ return x * __llvm_libc::fputil::polyeval(
+ x, 1.0f, 0.5f, 0.16666664183139801025390625f,
+ 4.1666664183139801025390625e-2f,
+ 8.3379410207271575927734375e-3f,
+ 1.3894210569560527801513671875e-3f);
+ }
+ if (x > 0.125f) {
+ return __llvm_libc::fputil::polyeval(
+ x, 1.23142086749794543720781803131103515625e-7f,
+ 0.9999969005584716796875f, 0.500031292438507080078125f,
+ 0.16650259494781494140625f, 4.21491153538227081298828125e-2f,
+ 7.53940828144550323486328125e-3f,
+ 2.05591344274580478668212890625e-3f);
+ }
+ return __llvm_libc::fputil::polyeval(
+ x, -6.899231408397099585272371768951416015625e-8f,
+ 0.999998271465301513671875f, 0.4999825656414031982421875f,
+ 0.16657467186450958251953125f, 4.1390590369701385498046875e-2f,
+ 7.856394164264202117919921875e-3f,
+ 9.380675037391483783721923828125e-4f);
+ }
+ return expf(x) - 1.0f;
+}
+
+} // namespace __llvm_libc
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 49b529048c31e..be0683932520a 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1168,6 +1168,20 @@ add_fp_unittest(
libc.utils.FPUtil.fputil
)
+add_fp_unittest(
+ expm1f_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ expm1f_test.cpp
+ DEPENDS
+ libc.include.errno
+ libc.include.math
+ libc.src.math.expm1f
+ libc.utils.FPUtil.fputil
+)
+
add_subdirectory(generic)
add_subdirectory(exhaustive)
add_subdirectory(
diff erential_testing)
diff --git a/libc/test/src/math/
diff erential_testing/CMakeLists.txt b/libc/test/src/math/
diff erential_testing/CMakeLists.txt
index 94a0aa6a73903..82f912cc7315d 100644
--- a/libc/test/src/math/
diff erential_testing/CMakeLists.txt
+++ b/libc/test/src/math/
diff erential_testing/CMakeLists.txt
@@ -106,3 +106,23 @@ add_
diff _binary(
COMPILE_OPTIONS
-fno-builtin
)
+
+add_
diff _binary(
+ expm1f_
diff
+ SRCS
+ expm1f_
diff .cpp
+ DEPENDS
+ .single_input_single_output_
diff
+ libc.src.math.expm1f
+)
+
+add_
diff _binary(
+ expm1f_perf
+ SRCS
+ expm1f_perf.cpp
+ DEPENDS
+ .single_input_single_output_
diff
+ libc.src.math.expm1f
+ COMPILE_OPTIONS
+ -fno-builtin
+)
diff --git a/libc/test/src/math/
diff erential_testing/expm1f_
diff .cpp b/libc/test/src/math/
diff erential_testing/expm1f_
diff .cpp
new file mode 100644
index 0000000000000..e12072efec57d
--- /dev/null
+++ b/libc/test/src/math/
diff erential_testing/expm1f_
diff .cpp
@@ -0,0 +1,16 @@
+//===-- Differential test for expm1f --------------------------------------===//
+//
+// 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 "SingleInputSingleOutputDiff.h"
+
+#include "src/math/expm1f.h"
+
+#include <math.h>
+
+SINGLE_INPUT_SINGLE_OUTPUT_DIFF(float, __llvm_libc::expm1f, ::expm1f,
+ "expm1f_
diff .log")
diff --git a/libc/test/src/math/
diff erential_testing/expm1f_perf.cpp b/libc/test/src/math/
diff erential_testing/expm1f_perf.cpp
new file mode 100644
index 0000000000000..9e8e2dfb8cd92
--- /dev/null
+++ b/libc/test/src/math/
diff erential_testing/expm1f_perf.cpp
@@ -0,0 +1,16 @@
+//===-- Differential test for expm1f --------------------------------------===//
+//
+// 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 "SingleInputSingleOutputDiff.h"
+
+#include "src/math/expm1f.h"
+
+#include <math.h>
+
+SINGLE_INPUT_SINGLE_OUTPUT_PERF(float, __llvm_libc::expm1f, ::expm1f,
+ "expm1f_perf.log")
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index e7b5b308f53f1..f6580da5966d8 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -38,3 +38,16 @@ add_fp_unittest(
libc.src.math.cosf
libc.utils.FPUtil.fputil
)
+
+add_fp_unittest(
+ expm1f_test
+ NEED_MPFR
+ SUITE
+ libc_math_exhaustive_tests
+ SRCS
+ expm1f_test.cpp
+ DEPENDS
+ libc.include.math
+ libc.src.math.expm1f
+ libc.utils.FPUtil.fputil
+)
diff --git a/libc/test/src/math/exhaustive/expm1f_test.cpp b/libc/test/src/math/exhaustive/expm1f_test.cpp
new file mode 100644
index 0000000000000..30d588f83dbdc
--- /dev/null
+++ b/libc/test/src/math/exhaustive/expm1f_test.cpp
@@ -0,0 +1,28 @@
+//===-- Exhaustive test for expm1f-----------------------------------------===//
+//
+// 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/math/expm1f.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include <math.h>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+TEST(LlvmLibcExpm1fExhaustiveTest, AllValues) {
+ uint32_t bits = 0;
+ do {
+ FPBits x(bits);
+ if (!x.isInfOrNaN() && float(x) < 88.70f) {
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, float(x),
+ __llvm_libc::expm1f(float(x)), 1.5);
+ }
+ } while (bits++ < 0xffff'ffffU);
+}
diff --git a/libc/test/src/math/expm1f_test.cpp b/libc/test/src/math/expm1f_test.cpp
new file mode 100644
index 0000000000000..d4b06c29c4bbb
--- /dev/null
+++ b/libc/test/src/math/expm1f_test.cpp
@@ -0,0 +1,141 @@
+//===-- Unittests for expm1f-----------------------------------------------===//
+//
+// 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/math/expm1f.h"
+#include "utils/FPUtil/BitPatterns.h"
+#include "utils/FPUtil/ClassificationFunctions.h"
+#include "utils/FPUtil/FloatOperations.h"
+#include "utils/FPUtil/FloatProperties.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+using __llvm_libc::fputil::isNegativeQuietNaN;
+using __llvm_libc::fputil::isQuietNaN;
+using __llvm_libc::fputil::valueAsBits;
+using __llvm_libc::fputil::valueFromBits;
+
+using BitPatterns = __llvm_libc::fputil::BitPatterns<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+TEST(LlvmLibcExpm1fTest, SpecialNumbers) {
+ errno = 0;
+
+ EXPECT_TRUE(
+ isQuietNaN(__llvm_libc::expm1f(valueFromBits(BitPatterns::aQuietNaN))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_TRUE(isNegativeQuietNaN(
+ __llvm_libc::expm1f(valueFromBits(BitPatterns::aNegativeQuietNaN))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_TRUE(isQuietNaN(
+ __llvm_libc::expm1f(valueFromBits(BitPatterns::aSignallingNaN))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_TRUE(isNegativeQuietNaN(
+ __llvm_libc::expm1f(valueFromBits(BitPatterns::aNegativeSignallingNaN))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_EQ(BitPatterns::inf,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(BitPatterns::inf))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_EQ(BitPatterns::negOne, valueAsBits(__llvm_libc::expm1f(
+ valueFromBits(BitPatterns::negInf))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_EQ(BitPatterns::zero,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(BitPatterns::zero))));
+ EXPECT_EQ(errno, 0);
+
+ EXPECT_EQ(BitPatterns::negZero, valueAsBits(__llvm_libc::expm1f(
+ valueFromBits(BitPatterns::negZero))));
+ EXPECT_EQ(errno, 0);
+}
+
+TEST(LlvmLibcExpm1fTest, Overflow) {
+ errno = 0;
+ EXPECT_EQ(BitPatterns::inf,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0x7f7fffffU))));
+ EXPECT_EQ(errno, ERANGE);
+
+ errno = 0;
+ EXPECT_EQ(BitPatterns::inf,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0x42cffff8U))));
+ EXPECT_EQ(errno, ERANGE);
+
+ errno = 0;
+ EXPECT_EQ(BitPatterns::inf,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0x42d00008U))));
+ EXPECT_EQ(errno, ERANGE);
+}
+
+TEST(LlvmLibcExpm1fTest, Underflow) {
+ errno = 0;
+ EXPECT_EQ(BitPatterns::negOne,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0xff7fffffU))));
+ EXPECT_EQ(errno, ERANGE);
+
+ errno = 0;
+ EXPECT_EQ(BitPatterns::negOne,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0xc2cffff8U))));
+ EXPECT_EQ(errno, ERANGE);
+
+ errno = 0;
+ EXPECT_EQ(BitPatterns::negOne,
+ valueAsBits(__llvm_libc::expm1f(valueFromBits(0xc2d00008U))));
+ EXPECT_EQ(errno, ERANGE);
+}
+
+// Test with inputs which are the borders of underflow/overflow but still
+// produce valid results without setting errno.
+TEST(LlvmLibcExpm1fTest, Borderline) {
+ float x;
+
+ errno = 0;
+ x = valueFromBits(0x42affff8U);
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+ EXPECT_EQ(errno, 0);
+
+ x = valueFromBits(0x42b00008U);
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+ EXPECT_EQ(errno, 0);
+
+ x = valueFromBits(0xc2affff8U);
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+ EXPECT_EQ(errno, 0);
+
+ x = valueFromBits(0xc2b00008U);
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+ EXPECT_EQ(errno, 0);
+}
+
+TEST(LlvmLibcExpm1fTest, InFloatRange) {
+ constexpr uint32_t count = 1000000;
+ constexpr uint32_t step = UINT32_MAX / count;
+ for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) {
+ float x = valueFromBits(v);
+ if (isnan(x) || isinf(x))
+ continue;
+ errno = 0;
+ float result = __llvm_libc::expm1f(x);
+
+ // If the computation resulted in an error or did not produce valid result
+ // in the single-precision floating point range, then ignore comparing with
+ // MPFR result as MPFR can still produce valid results because of its
+ // wider precision.
+ if (isnan(result) || isinf(result) || errno != 0)
+ continue;
+ ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.5);
+ }
+}
diff --git a/libc/utils/FPUtil/BitPatterns.h b/libc/utils/FPUtil/BitPatterns.h
index 439ccbcdcb92b..a26d403e1bf69 100644
--- a/libc/utils/FPUtil/BitPatterns.h
+++ b/libc/utils/FPUtil/BitPatterns.h
@@ -32,6 +32,7 @@ template <> struct BitPatterns<float> {
static constexpr BitsType negZero = 0x80000000U;
static constexpr BitsType one = 0x3f800000U;
+ static constexpr BitsType negOne = 0xbf800000U;
// Examples of quiet NAN.
static constexpr BitsType aQuietNaN = 0x7fc00000U;
diff --git a/libc/utils/FPUtil/CMakeLists.txt b/libc/utils/FPUtil/CMakeLists.txt
index c4600fd7f2bf8..75fa89d36b1d6 100644
--- a/libc/utils/FPUtil/CMakeLists.txt
+++ b/libc/utils/FPUtil/CMakeLists.txt
@@ -27,6 +27,7 @@ add_header_library(
ManipulationFunctions.h
NearestIntegerOperations.h
NormalFloat.h
+ PolyEval.h
DEPENDS
libc.include.math
libc.include.errno
diff --git a/libc/utils/FPUtil/PolyEval.h b/libc/utils/FPUtil/PolyEval.h
new file mode 100644
index 0000000000000..ce8377fe895a1
--- /dev/null
+++ b/libc/utils/FPUtil/PolyEval.h
@@ -0,0 +1,54 @@
+//===-- Common header for PolyEval implementations --------------*- C++ -*-===//
+//
+// 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_UTILS_FPUTIL_POLYEVAL_H
+#define LLVM_LIBC_UTILS_FPUTIL_POLYEVAL_H
+
+#include "utils/CPP/TypeTraits.h"
+
+// Evaluate polynomial using Horner's Scheme:
+// With polyeval(x, a_0, a_1, ..., a_n) = a_n * x^n + ... + a_1 * x + a_0, we
+// evaluated it as: a_0 + x * (a_1 + x * ( ... (a_(n-1) + x * a_n) ... ) ) ).
+// We will use fma instructions if available.
+// Example: to evaluate x^3 + 2*x^2 + 3*x + 4, call
+// polyeval( x, 4.0, 3.0, 2.0, 1.0 )
+
+#if defined(__x86_64__) || defined(__aarch64__)
+#include "FMA.h"
+
+namespace __llvm_libc {
+namespace fputil {
+
+template <typename T> static inline T polyeval(T x, T a0) { return a0; }
+
+template <typename T, typename... Ts>
+static inline T polyeval(T x, T a0, Ts... a) {
+ return fma(x, polyeval(x, a...), a0);
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#else
+
+namespace __llvm_libc {
+namespace fputil {
+
+template <typename T> static inline T polyeval(T x, T a0) { return a0; }
+
+template <typename T, typename... Ts>
+static inline T polyeval(T x, T a0, Ts... a) {
+ return x * polyeval(x, a...) + a0;
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif
+
+#endif // LLVM_LIBC_UTILS_FPUTIL_FMA_H
diff --git a/libc/utils/FPUtil/generic/FMA.h b/libc/utils/FPUtil/generic/FMA.h
index f88e28932932a..9471830af2c7d 100644
--- a/libc/utils/FPUtil/generic/FMA.h
+++ b/libc/utils/FPUtil/generic/FMA.h
@@ -69,6 +69,4 @@ static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y,
} // namespace fputil
} // namespace __llvm_libc
-#endif // Generic fma implementations
-
#endif // LLVM_LIBC_UTILS_FPUTIL_GENERIC_FMA_H
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 7eb2e4f5243f8..e040393bc26ce 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -139,6 +139,12 @@ class MPFRNumber {
return result;
}
+ MPFRNumber expm1() const {
+ MPFRNumber result;
+ mpfr_expm1(result.value, value, MPFR_RNDN);
+ return result;
+ }
+
MPFRNumber floor() const {
MPFRNumber result;
mpfr_floor(result.value, value);
@@ -309,6 +315,8 @@ unaryOperation(Operation op, InputType input) {
return mpfrInput.exp();
case Operation::Exp2:
return mpfrInput.exp2();
+ case Operation::Expm1:
+ return mpfrInput.expm1();
case Operation::Floor:
return mpfrInput.floor();
case Operation::Round:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 604a20d64c99a..ee98c8f7ecc2d 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -28,6 +28,7 @@ enum class Operation : int {
Cos,
Exp,
Exp2,
+ Expm1,
Floor,
Round,
Sin,
diff --git a/libc/utils/mathtools/expm1f.sollya b/libc/utils/mathtools/expm1f.sollya
new file mode 100644
index 0000000000000..d9dff4d413398
--- /dev/null
+++ b/libc/utils/mathtools/expm1f.sollya
@@ -0,0 +1,41 @@
+# Scripts to generate polynomial approximations for expm1f function using Sollya.
+#
+# To compute expm1f(x), for |x| > Ln(2), using expf(x) - 1.0f is accurate enough, since catastrophic
+# cancellation does not occur with the subtraction.
+#
+# For |x| <= Ln(2), we divide [-Ln2; Ln2] into 3 subintervals: [-Ln2; -1/8], [-1/8, 1/8], [1/8, Ln2],
+# and use a degree-6 polynomial to approximate expm1f in each interval.
+
+> f := expm1(x);
+
+# Polynomial approximation for e^(x) - 1 on [-Ln2, -1/8].
+> P1 := fpminmax(f, [|0, ..., 6|], [|24...], [-log(2), -1/8], relative);
+
+> log2(supnorm(P1, f, [-log(2), -1/8], relative, 2^(-50)));
+[-29.718757839645220560605567049447893449270454705067;-29.7187578396452193192777968211678241631166415833034]
+
+> P1;
+-6.899231408397099585272371768951416015625e-8 + x * (0.999998271465301513671875 + x * (0.4999825656414031982421875
++ x * (0.16657467186450958251953125 + x * (4.1390590369701385498046875e-2 + x * (7.856394164264202117919921875e-3
++ x * 9.380675037391483783721923828125e-4)))))
+
+# Polynomial approximation for e^(x) - 1 on [-1/8, 1/8].
+> P2 := fpminimax(f, [|1,...,6|], [|24...|], [-1/8, 1/8], relative);
+
+> log2(supnorm(P2, f, [-1/8, 1/8], relative, 2^(-50)));
+[-34.542864999883718873453825391741639571826398336605;-34.542864999883717632126055163461570285672585214842]
+
+> P2;
+x * (1 + x * (0.5 + x * (0.16666664183139801025390625 + x * (4.1666664183139801025390625e-2
++ x * (8.3379410207271575927734375e-3 + x * 1.3894210569560527801513671875e-3)))))
+
+# Polynomial approximation for e^(x) - 1 on [1/8, Ln2].
+> P3 := fpminimax(f, [|0,...,6|], [|24...|], [1/8, log(2)], relative);
+
+> log2(supnorm(P3, f, [1/8, log(2)], relative, 2^(-50)));
+[-29.189438260653683379922869677995123967174571911561;-29.1894382606536821385950994497150546810207587897976]
+
+> P3;
+1.23142086749794543720781803131103515625e-7 + x * (0.9999969005584716796875 + x * (0.500031292438507080078125
++ x * (0.16650259494781494140625 + x * (4.21491153538227081298828125e-2 + x * (7.53940828144550323486328125e-3
++ x * 2.05591344274580478668212890625e-3)))))
More information about the libc-commits
mailing list