[libc][math] Added atanhf function.

[orex]

Summary

Performance by core-math (core-math/glibc 2.31/current llvm-14):
10.845/43.174/13.467

The review is done on top of D132809.

[zimmermann6]

ok for me. I get slightly better figures for llvm-libc:

zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc.a CORE_MATH_PERF_MODE=rdtsc ./perf.sh atanhf
GNU libc version: 2.34
GNU libc release: stable
23.547
70.432
20.065

This is on a AMD EPYC 7282.

[orex]

Thank you, Paul! I'll land the revision after Tue review.

[lntue]

libc/src/math/generic/atanhf.cpp

22

|x| >= 1.0 ?

[orex]

Fixed comments.

[orex]

Fixed entrypoint.

[orex]

Review fixes.

[lntue]

libc/src/math/generic/CMakeLists.txt

1235

Update this dependency.

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.coshf
     libc.src.math.cosf

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanhf    
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.coshf
     libc.src.math.cosf

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanhf    
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     libc.src.math.coshf

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     libc.src.math.cosf

libc/spec/stdc.td

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

libc/src/math/CMakeLists.txt

HDRS
     fma.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.fma
   COMPILE_OPTIONS
     -O3
 )
 
+add_math_entrypoint_object(atanhf)
+
 add_math_entrypoint_object(ceil)
 add_math_entrypoint_object(ceilf)
 add_math_entrypoint_object(ceill)
 
 add_math_entrypoint_object(copysign)
 add_math_entrypoint_object(copysignf)
 add_math_entrypoint_object(copysignl)
 

libc/src/math/atanhf.h

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

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(
+  atanhf
+  SRCS
+    atanhf.cpp
+  HDRS
+    ../atanhf.h
+  DEPENDS
+    .supfuncf
+    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/atanhf.cpp

+//===-- Single-precision atanh 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/atanhf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/supfuncf.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, atanhf, (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| >= 1.0
+  if (unlikely(x_abs >= 0x3F80'0000U)) {
+    if (xbits.is_nan()) {
+      return x + 1.0f;
+    }
+    // |x| == 0
+    if (x_abs == 0x3F80'0000U) {
+      fputil::set_except(FE_DIVBYZERO);
+      return with_errno(FPBits::inf(sign).get_val(), ERANGE);
+    } else {
+      fputil::set_except(FE_INVALID);
+      return with_errno(
+          FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1)),
+          EDOM);
+    }
+  }
+
+  // |x| < ~0.10
+  if (unlikely(x_abs <= 0x3dcc'0000U)) {
+    // |x| <= 2^-26
+    if (unlikely(x_abs <= 0x3280'0000U)) {
+      return unlikely(x_abs == 0) ? x : (x + 0x1.5555555555555p-2 * x * x * x);
+    }
+
+    double xdbl = x;
+    double x2 = xdbl * xdbl;
+    // Pure Taylor series.
+    double pe = fputil::polyeval(x2, 0.0, 0x1.5555555555555p-2,
+                                 0x1.999999999999ap-3, 0x1.2492492492492p-3,
+                                 0x1.c71c71c71c71cp-4, 0x1.745d1745d1746p-4);
+    return fputil::multiply_add(xdbl, pe, xdbl);
+  }
+  double xdbl = x;
+  return 0.5 * log_eval((xdbl + 1.0) / (xdbl - 1.0));
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

SUITE
     libc_math_unittests
   SRCS
     tanhf_test.cpp
   DEPENDS
     libc.src.math.tanhf
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  atanhf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    atanhf_test.cpp
+  DEPENDS
+    libc.src.math.atanhf
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/atanhf_test.cpp

+//===-- Unittests for atanhf ----------------------------------------------===//
+//
+// 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/FPUtil/FPBits.h"
+#include "src/math/atanhf.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(LlvmLibcAtanhfTest, SpecialNumbers) {
+  errno = 0;
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(aNaN));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(0.0f, __llvm_libc::atanhf(0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::atanhf(-0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(inf, __llvm_libc::atanhf(1.0f));
+  EXPECT_FP_EXCEPTION(FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(neg_inf, __llvm_libc::atanhf(-1.0f));
+  EXPECT_FP_EXCEPTION(FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  auto bt = FPBits(1.0f);
+  bt.bits += 1;
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(bt.get_val()));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  bt.set_sign(true);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(bt.get_val()));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(2.0f));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(-2.0f));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(inf));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  bt.set_sign(true);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(neg_inf));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST(LlvmLibcAtanhfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  const uint32_t STEP = FPBits(1.0f).uintval() / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    ASSERT_MPFR_MATCH(mpfr::Operation::Atanh, x, __llvm_libc::atanhf(x), 0.5);
+    ASSERT_MPFR_MATCH(mpfr::Operation::Atanh, -x, __llvm_libc::atanhf(-x), 0.5);
+  }
+}
+
+// For small values, atanh(x) is x.
+TEST(LlvmLibcAtanhfTest, SmallValues) {
+  float x = float(FPBits(uint32_t(0x17800000)));
+  float result = __llvm_libc::atanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Atanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+
+  x = float(FPBits(uint32_t(0x00400000)));
+  result = __llvm_libc::atanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Atanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+}

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

DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.tanhf
     libc.src.__support.FPUtil.fputil
   LINK_LIBRARIES
     -lpthread
 )
 
+add_fp_unittest(
+  atanhf_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    atanhf_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.atanhf
+    libc.src.__support.FPUtil.fputil
+  LINK_LIBRARIES
+    -lpthread
+)

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

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

libc/utils/MPFRWrapper/MPFRUtils.h

 namespace mpfr {
 
 enum class Operation : int {
   // Operations with take a single floating point number as input
   // and produce a single floating point number as output. The input
   // and output floating point numbers are of the same kind.
   BeginUnaryOperationsSingleOutput,
   Abs,
+  Atanh,
   Ceil,
   Cos,
   Cosh,
   Exp,
   Exp2,
   Expm1,
   Floor,
   Log,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber abs() const {
     MPFRNumber result(*this);
     mpfr_abs(result.value, value, mpfr_rounding);
     return result;
   }
 
+  MPFRNumber atanh() const {
+    MPFRNumber result(*this);
+    mpfr_atanh(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber ceil() const {
     MPFRNumber result(*this);
     mpfr_ceil(result.value, value);
     return result;
   }
 
   MPFRNumber cos() const {
     MPFRNumber result(*this);

 template <typename InputType>
 cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
 unary_operation(Operation op, InputType input, unsigned int precision,
                 RoundingMode rounding) {
   MPFRNumber mpfrInput(input, precision, rounding);
   switch (op) {
   case Operation::Abs:
     return mpfrInput.abs();
+  case Operation::Atanh:
+    return mpfrInput.atanh();
   case Operation::Ceil:
     return mpfrInput.ceil();
   case Operation::Cos:
     return mpfrInput.cos();
   case Operation::Cosh:
     return mpfrInput.cosh();
   case Operation::Exp:
     return mpfrInput.exp();