diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -187,6 +187,7 @@ libc.src.math.roundf libc.src.math.roundl libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf 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 @@ -206,6 +206,7 @@ libc.src.math.roundf libc.src.math.roundl libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf 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 @@ -213,6 +213,7 @@ libc.src.math.roundl libc.src.math.sin libc.src.math.sincosf + libc.src.math.sinhf libc.src.math.sinf libc.src.math.sqrt libc.src.math.sqrtf diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -191,6 +191,7 @@ libc.src.math.sin libc.src.math.sincosf libc.src.math.sinf + libc.src.math.sinhf libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -472,6 +472,7 @@ FunctionSpec<"nextafterl", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"coshf", RetValSpec, [ArgSpec]>, + FunctionSpec<"sinhf", RetValSpec, [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 @@ -182,6 +182,7 @@ add_math_entrypoint_object(sin) add_math_entrypoint_object(sinf) +add_math_entrypoint_object(sinhf) add_math_entrypoint_object(sqrt) add_math_entrypoint_object(sqrtf) 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 @@ -1147,3 +1147,21 @@ -O3 ) +add_entrypoint_object( + sinhf + SRCS + sinhf.cpp + HDRS + ../sinhf.h + expxf.h + DEPENDS + .common_constants + libc.src.__support.FPUtil.fputil + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.FPUtil.polyeval + libc.include.math + COMPILE_OPTIONS + -O3 +) + diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/sinhf.cpp @@ -0,0 +1,84 @@ +//===-- Single-precision sinh 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/sinhf.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/generic/expxf.h" + +namespace __llvm_libc { + +LLVM_LIBC_FUNCTION(float, sinhf, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + bool sign = xbits.get_sign(); + uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK; + + // |x| <= 2^-26 + if (unlikely(x_abs <= 0x3280'0000U)) { + return unlikely(x_abs == 0) ? x : (x + 0.25 * x * x * x); + } + + // When |x| >= 90, or x is inf or nan + if (unlikely(x_abs >= 0x42b4'0000U)) { + if (xbits.is_nan()) + return x + 1.0f; // sNaN to qNaN + signal + + if (xbits.is_inf()) + return x; + + int rounding = fputil::get_round(); + if (sign) { + if (unlikely(rounding == FE_UPWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK) + .get_val(); + } else { + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); + } + + errno = ERANGE; + + return x + FPBits::inf(sign).get_val(); + } + + // |x| <= 0.078125 + if (unlikely(x_abs <= 0x3da0'0000U)) { + // |x| = 0.0005589424981735646724700927734375 + if (unlikely(x_abs == 0x3a12'85ffU)) { + if (fputil::get_round() == FE_TONEAREST) + return x; + } + + double xdbl = x; + double x2 = xdbl * xdbl; + // Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x); + // Sollya output: x * (0x1p0 + x^0x1p1 * (0x1.5555555556583p-3 + x^0x1p1 + // * (0x1.111110d239f1fp-7 + // + x^0x1p1 * 0x1.a02b5a284013cp-13))) + // Therefore, output of Sollya = x * pe; + double pe = fputil::polyeval(x2, 0.0, 0x1.5555555556583p-3, + 0x1.111110d239f1fp-7, 0x1.a02b5a284013cp-13); + return fputil::multiply_add(xdbl, pe, xdbl); + } + + // MULT_POWER2 = -1 + auto ep_p = exp_eval<-1>(x); // 0.5 * exp(x) + auto ep_m = exp_eval<-1>(-x); // 0.5 * exp(-x) + + // 0.5 * expm1(x) = ep_p.mult_exp * (ep_p.r + 1) - 0.5 + // = ep_p.mult_exp * ep_p.r + ep_p.mult_exp - 0.5 + // 0.5 * expm1(-x) = ep_m.mult_exp * (ep_m.r + 1) - 0.5 + // = ep_m.mult_exp * ep_m.r + ep_m.mult_exp - 0.5 + // sinh(x) = 0.5 * expm1(x) - 0.5 * expm1(-x) + // Using expm1 instead of exp improved precision around zero. + double ep = fputil::multiply_add(ep_p.mult_exp, ep_p.r, ep_p.mult_exp - 0.5) - + fputil::multiply_add(ep_m.mult_exp, ep_m.r, ep_m.mult_exp - 0.5); + return ep; +} + +} // namespace __llvm_libc diff --git a/libc/src/math/sinhf.h b/libc/src/math/sinhf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/sinhf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for sinhf -------------------------*- 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_SINHF_H +#define LLVM_LIBC_SRC_MATH_SINHF_H + +namespace __llvm_libc { + +float sinhf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINHF_H 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 @@ -1330,6 +1330,22 @@ libc.src.__support.FPUtil.fputil ) +add_fp_unittest( + sinhf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + sinhf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + libc.include.errno + libc.src.math.sinhf + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fputil +) + add_subdirectory(generic) add_subdirectory(exhaustive) add_subdirectory(differential_testing) 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 @@ -219,3 +219,20 @@ -lpthread ) +add_fp_unittest( + sinhf_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + sinhf_test.cpp + DEPENDS + .exhaustive_test + libc.include.math + libc.src.math.sinhf + libc.src.__support.FPUtil.fputil + LINK_LIBRARIES + -lpthread +) + diff --git a/libc/test/src/math/exhaustive/sinhf_test.cpp b/libc/test/src/math/exhaustive/sinhf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exhaustive/sinhf_test.cpp @@ -0,0 +1,75 @@ +//===-- Exhaustive test for sinhf -----------------------------------------===// +// +// 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 "exhaustive_test.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/sinhf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" + +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcSinhfExhaustiveTest : public LlvmLibcExhaustiveTest { + bool check(uint32_t start, uint32_t stop, + mpfr::RoundingMode rounding) override { + mpfr::ForceRoundingMode r(rounding); + uint32_t bits = start; + bool result = true; + do { + FPBits xbits(bits); + float x = float(xbits); + result &= EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, + __llvm_libc::sinhf(x), 0.5, rounding); + } while (bits++ < stop); + return result; + } +}; + +// Range: [0, 90]; +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x42b4'0000U; + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero); +} + +// Range: [-90, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xc2b4'0000U; + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundUp) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundDown) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero); +} diff --git a/libc/test/src/math/sinhf_test.cpp b/libc/test/src/math/sinhf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/sinhf_test.cpp @@ -0,0 +1,89 @@ +//===-- Unittests for sinhf -----------------------------------------------===// +// +// 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/CPP/Array.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/sinhf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" +#include + +#include +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +DECLARE_SPECIAL_CONSTANTS(float) + +TEST(LlvmLibcSinhfTest, SpecialNumbers) { + errno = 0; + + EXPECT_FP_EQ(aNaN, __llvm_libc::sinhf(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(0.0f, __llvm_libc::sinhf(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(-0.0f, __llvm_libc::sinhf(-0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(neg_inf, __llvm_libc::sinhf(neg_inf)); + EXPECT_MATH_ERRNO(0); +} + +TEST(LlvmLibcSinhfTest, 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 = float(FPBits(v)); + if (isnan(x) || isinf(x)) + continue; + ASSERT_MPFR_MATCH(mpfr::Operation::Sinh, x, __llvm_libc::sinhf(x), 0.5); + } +} + +// For small values, sinh(x) is x. +TEST(LlvmLibcSinhfTest, SmallValues) { + float x = float(FPBits(uint32_t(0x17800000))); + float result = __llvm_libc::sinhf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 0.5); + EXPECT_FP_EQ(x, result); + + x = float(FPBits(uint32_t(0x00400000))); + result = __llvm_libc::sinhf(x); + EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 0.5); + EXPECT_FP_EQ(x, result); +} + +TEST(LlvmLibcSinhfTest, Overflow) { + errno = 0; + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x7f7fffffU)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42cffff8U)))); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42d00008U)))); + EXPECT_MATH_ERRNO(ERANGE); +} + +TEST(LlvmLibcSinhfTest, ExceptionalValues) { + float x = float(FPBits(uint32_t(0x3a12'85ffU))); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x, + __llvm_libc::sinhf(x), 0.5); + + x = -float(FPBits(uint32_t(0x3a12'85ffU))); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x, + __llvm_libc::sinhf(x), 0.5); +} 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 @@ -41,6 +41,7 @@ ModPIOver4, Round, Sin, + Sinh, Sqrt, Tan, Trunc, 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 @@ -337,6 +337,12 @@ return result; } + MPFRNumber sinh() const { + MPFRNumber result(*this); + mpfr_sinh(result.value, value, mpfr_rounding); + return result; + } + MPFRNumber sqrt() const { MPFRNumber result(*this); mpfr_sqrt(result.value, value, mpfr_rounding); @@ -515,6 +521,8 @@ return mpfrInput.round(); case Operation::Sin: return mpfrInput.sin(); + case Operation::Sinh: + return mpfrInput.sinh(); case Operation::Sqrt: return mpfrInput.sqrt(); case Operation::Tan: