diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -68,6 +68,7 @@ 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 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -68,6 +68,7 @@ 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 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -393,6 +393,7 @@ FunctionSpec<"expf", RetValSpec, [ArgSpec]>, FunctionSpec<"exp2f", RetValSpec, [ArgSpec]>, + FunctionSpec<"expm1f", RetValSpec, [ArgSpec]>, FunctionSpec<"remainderf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"remainder", RetValSpec, [ArgSpec, ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -70,6 +70,8 @@ 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 --- /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 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -489,6 +489,18 @@ 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 --- /dev/null +++ b/libc/src/math/generic/expm1f.cpp @@ -0,0 +1,72 @@ +//===-- 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 { + +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) { + // We divide [-Ln2; Ln2] into 3 subintervals [-Ln2; -1/8], [-1/8; 1/8], + // [1/8; Ln2]. And we use a degree-5 polynomial to approximate exp(x) - 1 in + // each interval. The coefficients were generated by Sollya's fpminmax. + if (abs_x < 0.25f) { + // From Sollya command: + // > fpminimax(f, [|1, ..., 6|], [|1, 24...|], [-1/4, 1/4], relative); + // x * (1.0f + x * (0.5f + x * (0.166666328907012939453125f + x * + // (4.1666679084300994873046875e-2f + x * (8.350677788257598876953125e-3f + // + x * 1.39005668461322784423828125e-3f))))); + return x * __llvm_libc::fputil::polyeval( + x, 1.0f, 0.5f, 0.166666328907012939453125f, + 4.1666679084300994873046875e-2f, + 8.350677788257598876953125e-3f, + 1.39005668461322784423828125e-3f); + } + if (x >= 0.25f) { + // From Sollya command: + // > fpminimax(f, [|0, ..., 5|], [|24...|], [1/8, log(2)], relative); + // 8.6965019363560713827610015869140625e-7f + x * + // (0.999985158443450927734375f + x * (0.500105917453765869140625f + x * + // (0.16625757515430450439453125f + x * (4.25875484943389892578125e-2f + x + // * (7.1335560642182826995849609375e-3f + x + // * 2.208019606769084930419921875e-3f))))); + return __llvm_libc::fputil::polyeval( + x, 8.6965019363560713827610015869140625e-7f, + 0.999985158443450927734375f, 0.500105917453765869140625f, + 0.16625757515430450439453125f, 4.25875484943389892578125e-2f, + 7.1335560642182826995849609375e-3f, + 2.208019606769084930419921875e-3f); + } + // From Sollya command: + // > fpminimax(f, [|0, ..., 5|], [|24...|], [-log(2), -1/4], relative); + // -4.05176479034707881510257720947265625e-7f + x * + // (0.999993026256561279296875f + x * (0.4999495446681976318359375f + x * + // (0.166467368602752685546875f + x * (4.12006340920925140380859375e-2f + x + // * (7.68242776393890380859375e-3f + x + // * 8.735431474633514881134033203125e-4f))))); + return __llvm_libc::fputil::polyeval( + x, -4.05176479034707881510257720947265625e-7f, + 0.999993026256561279296875f, 0.4999495446681976318359375f, + 0.166467368602752685546875f, 4.12006340920925140380859375e-2f, + 7.68242776393890380859375e-3f, 8.735431474633514881134033203125e-4f); + } + // When |x| > Ln2, catastrophic cancellation does not occur with the + // subtruction expf(x) - 1.0f. + return expf(x) - 1.0f; +} + +} // namespace __llvm_libc diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1132,5 +1132,19 @@ 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) diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -12,3 +12,17 @@ libc.src.math.sqrtf 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 +) \ No newline at end of file diff --git a/libc/test/src/math/exhaustive/expm1f_test.cpp b/libc/test/src/math/exhaustive/expm1f_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exhaustive/expm1f_test.cpp @@ -0,0 +1,29 @@ +//===-- 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 + +using FPBits = __llvm_libc::fputil::FPBits; + +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 --- /dev/null +++ b/libc/test/src/math/expm1f_test.cpp @@ -0,0 +1,142 @@ +//===-- 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 + +#include +#include + +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; + +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), 2.0); + } +} diff --git a/libc/utils/FPUtil/BitPatterns.h b/libc/utils/FPUtil/BitPatterns.h --- a/libc/utils/FPUtil/BitPatterns.h +++ b/libc/utils/FPUtil/BitPatterns.h @@ -32,6 +32,7 @@ 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 --- a/libc/utils/FPUtil/CMakeLists.txt +++ b/libc/utils/FPUtil/CMakeLists.txt @@ -27,6 +27,7 @@ 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 --- /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 static inline T polyeval(T x, T a0) { return a0; } + +template +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 static inline T polyeval(T x, T a0) { return a0; } + +template +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 --- a/libc/utils/FPUtil/generic/FMA.h +++ b/libc/utils/FPUtil/generic/FMA.h @@ -69,6 +69,4 @@ } // namespace fputil } // namespace __llvm_libc -#endif // Generic fma implementations - #endif // LLVM_LIBC_UTILS_FPUTIL_GENERIC_FMA_H diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -28,6 +28,7 @@ Cos, Exp, Exp2, + Expm1, Floor, Round, Sin, diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -139,6 +139,12 @@ 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); @@ -303,6 +309,8 @@ return mpfrInput.exp(); case Operation::Exp2: return mpfrInput.exp2(); + case Operation::Expm1: + return mpfrInput.expm1(); case Operation::Floor: return mpfrInput.floor(); case Operation::Round: