[libc][math] Added tanhf function.

[orex]

Summary

Correct rounding function. Performance ~2x faster than glibc analog.

Performance (llvm 12 intel):

CORE_MATH_PERF_MODE=rdtsc PERF_ARGS='' ./perf.sh tanhf
GNU libc version: 2.31
GNU libc release: stable
13.279
37.492
18.145
CORE_MATH_PERF_MODE=rdtsc PERF_ARGS='--latency' ./perf.sh tanhf
GNU libc version: 2.31
GNU libc release: stable
40.658
109.582
66.568

[lntue]

The implementation passed exhaustive tests for all 3 modes: SSE2, SSE4.2, and FMA.

And here is the performance number on Ryzen 1700:

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh tanhf
CORE-MATH reciprocal throughput   : 13.379
System LIBC reciprocal throughput : 55.843
LIBC reciprocal throughput        : 25.130

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh tanhf --latency
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH latency   : 43.814
System LIBC latency : 125.128
LIBC latency        : 94.906

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     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.tanhf
     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.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.tanhf
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 
     # sys/mman.h entrypoints
     libc.src.sys.mman.mmap
     libc.src.sys.mman.munmap
 )

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     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.tanhf
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 )
 
 if(LLVM_LIBC_FULL_BUILD)
   list(APPEND TARGET_LIBC_ENTRYPOINTS
     # assert.h entrypoints

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     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.tanhf
     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<"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>]>,
+          FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 
   HeaderSpec StdIO = HeaderSpec<
       "stdio.h",
       [
           Macro<"stderr">,
           Macro<"stdout">,

libc/src/math/CMakeLists.txt

 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(tanhf)
 
 add_math_entrypoint_object(trunc)
 add_math_entrypoint_object(truncf)
 add_math_entrypoint_object(truncl)

libc/src/math/generic/CMakeLists.txt

DEPENDS
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
+add_entrypoint_object(
+  tanhf
+  SRCS
+    tanhf.cpp
+  HDRS
+    ../tanhf.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
+)
+

libc/src/math/generic/tanhf.cpp

+//===-- Single-precision tanh 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/tanhf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/expxf.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, tanhf, (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 - 0x1.5555555555555p-2 * x * x * x);
+  }
+
+  // When |x| >= 15, or x is inf or nan
+  if (unlikely(x_abs >= 0x4170'0000U)) {
+    if (xbits.is_nan())
+      return x + 1.0f; // sNaN to qNaN + signal
+
+    if (xbits.is_inf())
+      return sign ? -1.0f : 1.0f;
+
+    if (sign) {
+      return -1.0f + opt_barrier(FPBits(FPBits::MIN_NORMAL).get_val());
+    } else
+      return 1.0f - opt_barrier(FPBits(FPBits::MIN_NORMAL).get_val());
+  }
+
+  // |x| <= 0.078125
+  if (unlikely(x_abs <= 0x3da0'0000U)) {
+    double xdbl = x;
+    double x2 = xdbl * xdbl;
+    // Pure Taylor series.
+    double pe = fputil::polyeval(x2, 0.0, -0x1.5555555555555p-2,
+                                 0x1.1111111111111p-3, -0x1.ba1ba1ba1ba1cp-5,
+                                 0x1.664f4882c10fap-6, -0x1.226e355e6c23dp-7);
+    return fputil::multiply_add(xdbl, pe, xdbl);
+  }
+
+  if (unlikely(xbits.bits == 0x4058'e0a3U)) {
+    if (fputil::get_round() == FE_DOWNWARD)
+      return FPBits(0x3f7f'6ad9U).get_val();
+  }
+
+  auto ep = exp_eval(2.0f * (sign ? x : -x)); // exp(-2 * x)
+  double result = fputil::multiply_add(ep.mult_exp, ep.r, ep.mult_exp - 1.0) /
+                  (fputil::multiply_add(ep.mult_exp, ep.r, ep.mult_exp + 1.0));
+  return sign ? result : -result;
+}
+
+} // namespace __llvm_libc

libc/src/math/tanhf.h

+//===-- Implementation header for tanhf -------------------------*- 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_TANHF_H
+#define LLVM_LIBC_SRC_MATH_TANHF_H
+
+namespace __llvm_libc {
+
+float tanhf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_TANHF_H

libc/test/src/math/CMakeLists.txt

HDRS
     sdcomp26094.h
   DEPENDS
     libc.include.errno
     libc.src.math.sinhf
     libc.src.__support.CPP.array
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  tanhf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    tanhf_test.cpp
+  DEPENDS
+    libc.src.math.tanhf
+    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.sinhf
     libc.src.__support.FPUtil.fputil
   LINK_LIBRARIES
     -lpthread
 )
 
+add_fp_unittest(
+  tanhf_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    tanhf_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.tanhf
+    libc.src.__support.FPUtil.fputil
+  LINK_LIBRARIES
+    -lpthread
+)
+

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

+//===-- Exhaustive test for tanhf -----------------------------------------===//
+//
+// 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/tanhf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcTanhfExhaustiveTest : 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::Tanh, x,
+                                  __llvm_libc::tanhf(x), 0.5, rounding);
+    } while (bits++ < stop);
+    return result;
+  }
+};
+
+static const int NUM_THREADS = std::thread::hardware_concurrency();
+
+// Range: [0, INF];
+static const uint32_t POS_START = 0x0000'0000U;
+static const uint32_t POS_STOP = FPBits::inf(false).uintval();
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, PostiveRangeRoundUp) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, PostiveRangeRoundDown) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, PostiveRangeRoundTowardZero) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
+}
+
+// Range: [-INF, 0];
+static const uint32_t NEG_START = 0x8000'0000U;
+static const uint32_t NEG_STOP = FPBits::inf(true).uintval();
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, NegativeRangeRoundUp) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, NegativeRangeRoundDown) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcTanhfExhaustiveTest, NegativeRangeRoundTowardZero) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/tanhf_test.cpp

+//===-- Unittests for tanhf -----------------------------------------------===//
+//
+// 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/tanhf.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(LlvmLibcTanhfTest, SpecialNumbers) {
+  errno = 0;
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::tanhf(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(0.0f, __llvm_libc::tanhf(0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::tanhf(-0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(1.0f, __llvm_libc::tanhf(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(-1.0f, __llvm_libc::tanhf(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST(LlvmLibcTanhfTest, 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::Tanh, x, __llvm_libc::tanhf(x), 0.5);
+  }
+}
+
+// For small values, tanh(x) is x.
+TEST(LlvmLibcTanhfTest, SmallValues) {
+  float x = float(FPBits(uint32_t(0x17800000)));
+  float result = __llvm_libc::tanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Tanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+
+  x = float(FPBits(uint32_t(0x00400000)));
+  result = __llvm_libc::tanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Tanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+}
+
+TEST(LlvmLibcTanhfTest, ExceptionalValues) {
+  float x = float(FPBits(uint32_t(0x3a12'85ffU)));
+  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, x,
+                                 __llvm_libc::tanhf(x), 0.5);
+
+  x = -float(FPBits(uint32_t(0x3a12'85ffU)));
+  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, x,
+                                 __llvm_libc::tanhf(x), 0.5);
+}

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   Mod2PI,
   ModPIOver2,
   ModPIOver4,
   Round,
   Sin,
   Sinh,
   Sqrt,
   Tan,
+  Tanh,
   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
   // 'int'.
   BeginUnaryOperationsTwoOutputs,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber tan() const {
     MPFRNumber result(*this);
     mpfr_tan(result.value, value, mpfr_rounding);
     return result;
   }
 
+  MPFRNumber tanh() const {
+    MPFRNumber result(*this);
+    mpfr_tanh(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber trunc() const {
     MPFRNumber result(*this);
     mpfr_trunc(result.value, value);
     return result;
   }
 
   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
     MPFRNumber result(*this);
unary_operation(Operation op, InputType input, unsigned int precision,
   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::Tanh:
+    return mpfrInput.tanh();
   case Operation::Trunc:
     return mpfrInput.trunc();
   default:
     __builtin_unreachable();
   }
 }
 
 template <typename InputType>