[libc][math] Added sinhf function.

[orex]

Summary

Latest performance.

CORE_MATH_PERF_MODE=rdtsc PERF_ARGS='' ./perf.sh sinhf
GNU libc version: 2.31
GNU libc release: stable
17.897
38.125
13.712
CORE_MATH_PERF_MODE=rdtsc PERF_ARGS='--latency' ./perf.sh sinhf
GNU libc version: 2.31
GNU libc release: stable
48.845
113.762
51.198

[zimmermann6]

should I apply this patch also on top of revision 60d6be5 + D129005 + D129215 ?

[orex]

The easiest way, from my point of view, will be to apply it on top of working coshf branch, which you made and test.

[zimmermann6]

In D129278#3645231, @orex wrote:

The easiest way, from my point of view, will be to apply it on top of working coshf branch, which you made and test.

thank you, this worked perfectly. All exhaustive tests do pass, and the timings are as follows (same processor and compilers as D129275):

zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc.a CORE_MATH_LAUNCHER="/localdisk/zimmerma/glibc-2.35/install/lib/ld-linux-x86-64.so.2 --library-path /localdisk/zimmerma/glibc-2.35/install/lib" CORE_MATH_PERF_MODE=rdtsc ./perf.sh sinhf
GNU libc version: 2.35
GNU libc release: stable
17.417
50.935
32.774

thus 17 cycles for CORE-MATH, 51 cycles for glibc 2.35, and 33 cycles for llvm-libc.

[lntue]

libc/src/math/generic/sinhf.cpp

60

It should either be:

// Sollya: fpminimax(sinh(x),[|1, 3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);

or

// Sollya: fpminimax(sinh(x) - x,[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);

61–63

Will this change the accuracy / throughput / latency for this range of inputs?

double x3 = x2 * xdbl;
double pe = fputil::polyeval(x2, 0x1.5555555556583p-3, 0x1.111110d239f1fp-7, 0x1.a02b5a284013cp-13);
return fputil::multiply_add(x3, pe, xdbl);

68–70

Add comments about your expanded formula / computations:

exp(x) = ep_p.mult_exp * (ep_p.r + 1)
exp(-x) = ep_m.mult_exp * (ep_m.r + 1)
sinh(x) = (exp(x) - exp(-x)) / 2 = ...

[lntue]

libc/test/src/math/sinhf_test.cpp

53

0.5 for tolerance?

61

0.5 for tolerance?

66

0.5 for tolerance?

[orex]

Rebasing on main with small fixes.

[orex]

Review fixes. More verbose comments.

libc/src/math/generic/sinhf.cpp

61–63

Performance:

I've checked it with core-math, reducing the range to [-0.08, 0.08]. The performance was the same with ~5% precision.

Precision:

The changes did not improve precision (the exceptional values still exceptional).

[lntue]

libc/src/math/generic/CMakeLists.txt

1152

Add expxf.h to HDRS and nearest_integer to DEPENDS

[orex]

Review fixes.

[lntue]

Looks good to me. Let's wait for @zimmermann6 to confirm the accuracy and performance.

[lntue]

@zimmermann6 is on vacation. You can go ahead and land this patch. He can comment on it later if there is any problem.

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.remquol
     libc.src.math.rint
     libc.src.math.rintf
     libc.src.math.rintl
     libc.src.math.round
     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
     libc.src.math.sqrtl
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 )
 
 set(TARGET_LLVMLIBC_ENTRYPOINTS
   ${TARGET_LIBC_ENTRYPOINTS}
   ${TARGET_LIBM_ENTRYPOINTS}
 )

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.remquol
     libc.src.math.rint
     libc.src.math.rintf
     libc.src.math.rintl
     libc.src.math.round
     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
     libc.src.math.sqrtl
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.rint
     libc.src.math.rintf
     libc.src.math.rintl
     libc.src.math.round
     libc.src.math.roundf
     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
     libc.src.math.sqrtl
     libc.src.math.tan
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.rintf
     libc.src.math.rintl
     libc.src.math.round
     libc.src.math.roundf
     libc.src.math.roundl
     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
     libc.src.math.tan
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 )
 
 set(TARGET_LLVMLIBC_ENTRYPOINTS
   ${TARGET_LIBC_ENTRYPOINTS}
   ${TARGET_LIBM_ENTRYPOINTS}
 )
 

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"nearbyintf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"nearbyintl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"nextafterf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"nextafter", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"nextafterl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"coshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 
   HeaderSpec StdIO = HeaderSpec<
       "stdio.h",
       [
           Macro<"stderr">,
           Macro<"stdout">,

libc/src/math/CMakeLists.txt

 add_math_entrypoint_object(round)
 add_math_entrypoint_object(roundf)
 add_math_entrypoint_object(roundl)
 
 add_math_entrypoint_object(sincosf)
 
 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)
 add_math_entrypoint_object(sqrtl)
 
 add_math_entrypoint_object(tan)
 
 add_math_entrypoint_object(trunc)
 add_math_entrypoint_object(truncf)
 add_math_entrypoint_object(truncl)

libc/src/math/generic/CMakeLists.txt

DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
+add_entrypoint_object(
+  sinhf
+  SRCS
+    sinhf.cpp
+  HDRS
+    ../sinhf.h
+  DEPENDS
+    .common_constants
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.polyeval
+    libc.include.math
+  COMPILE_OPTIONS
+    -O3
+)
+

libc/src/math/generic/sinhf.cpp

+//===-- 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<float>;
+  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

libc/src/math/sinhf.h

+//===-- 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

libc/test/src/math/CMakeLists.txt

HDRS
     sdcomp26094.h
   DEPENDS
     libc.include.errno
     libc.src.math.coshf
     libc.src.__support.CPP.array
     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)

libc/test/src/math/exhaustive/CMakeLists.txt

DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.coshf
     libc.src.__support.FPUtil.fputil
   LINK_LIBRARIES
     -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
+)
+

libc/test/src/math/exhaustive/sinhf_test.cpp

+//===-- Exhaustive test for expf ------------------------------------------===//
+//
+// 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 <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcSinhfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+  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;
+  }
+};
+
+static const int NUM_THREADS = std::thread::hardware_concurrency();
+
+// 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);
+}

libc/test/src/math/sinhf_test.cpp

+//===-- 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 <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+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);
+}

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   Log2,
   Log10,
   Log1p,
   Mod2PI,
   ModPIOver2,
   ModPIOver4,
   Round,
   Sin,
+  Sinh,
   Sqrt,
   Tan,
   Trunc,
   EndUnaryOperationsSingleOutput,
 
   // Operations which take a single floating point nubmer as input
   // but produce two outputs. The first ouput is a floating point
   // number of the same type as the input. The second output is of type

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber sin() const {
     MPFRNumber result(*this);
     mpfr_sin(result.value, value, mpfr_rounding);
     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);
     return result;
   }
 
   MPFRNumber tan() const {
     MPFRNumber result(*this);
unary_operation(Operation op, InputType input, unsigned int precision,
   case Operation::ModPIOver2:
     return mpfrInput.mod_pi_over_2();
   case Operation::ModPIOver4:
     return mpfrInput.mod_pi_over_4();
   case Operation::Round:
     return mpfrInput.round();
   case Operation::Sin:
     return mpfrInput.sin();
+  case Operation::Sinh:
+    return mpfrInput.sinh();
   case Operation::Sqrt:
     return mpfrInput.sqrt();
   case Operation::Tan:
     return mpfrInput.tan();
   case Operation::Trunc:
     return mpfrInput.trunc();
   default:
     __builtin_unreachable();