[libc] Add a new version log_eval.

Author: harrisonGPU

libc/src/math/generic/atanhf16.cpp

@@ -44,8 +44,13 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
 
   // |x| >= 1
   if (LIBC_UNLIKELY(x_abs >= 0x3c00U)) {
-    if (xbits.is_nan())
+    if (xbits.is_nan()) {
+      if (xbits.is_signaling_nan()) {
+        fputil::raise_except_if_required(FE_INVALID);
+        return FPBits::quiet_nan().get_val();
+      }
       return x;
+    }
 
     // |x| == 1.0
     if (x_abs == 0x3c00U) {
@@ -93,7 +98,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
   }
 
   float xf = x;
-  return fputil::cast<float16>(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f)));
+  return fputil::cast<float16>(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f)));
 }
 
 } // namespace LIBC_NAMESPACE_DECL

libc/src/math/generic/common_constants.cpp

@@ -15,6 +15,9 @@ namespace LIBC_NAMESPACE_DECL {
 
 // Lookup table for logf(f) = logf(1 + n*2^(-7)) where n = 0..127,
 // computed and stored as float precision constants.
+// Generated by Sollya with the following commands:
+//   display = hexadecimal;
+//   for n from 0 to 127 do { print(single(1 / (1 + n / 128.0))); };
 const float ONE_OVER_F_FLOAT[128] = {
     0x1p0f,         0x1.fc07fp-1f,  0x1.f81f82p-1f, 0x1.f4465ap-1f,
     0x1.f07c2p-1f,  0x1.ecc07cp-1f, 0x1.e9131ap-1f, 0x1.e573acp-1f,
@@ -97,6 +100,9 @@ const double ONE_OVER_F[128] = {
 
 // Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127,
 // computed and stored as float precision constants.
+// Generated by Sollya with the following commands:
+//   display = hexadecimal;
+//   for n from 0 to 127 do { print(single(log(1 + n / 128.0))); };
 const float LOG_F_FLOAT[128] = {
     0.0f,           0x1.fe02a6p-8f, 0x1.fc0a8cp-7f, 0x1.7b91bp-6f,
     0x1.f829bp-6f,  0x1.39e87cp-5f, 0x1.77459p-5f,  0x1.b42dd8p-5f,

libc/src/math/generic/explogxf.h

@@ -133,6 +133,29 @@ struct exp_b_reduc_t {
   double lo;
 };
 
+// Coefficients for double (6th-degree minimax polynomial on [0, 2^-7]).
+// Minimax polynomial of log(1 + dx) generated by Sollya with:
+// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
+static constexpr double LOG_COEFFS_DOUBLE[6] = {
+    -0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2,  -0x1.ffffffefe562dp-3,
+    0x1.9999817d3a50fp-3,  -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3};
+
+// Coefficients for float (6th-degree minimax polynomial on [0, 2^-7]).
+// Minimax polynomial of log(1 + dx) generated by Sollya with:
+// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
+static constexpr float LOG_COEFFS_FLOAT[6] = {-0x1.fffffep-2f, 0x1.555556p-2f,
+                                              -0x1.fffefep-3f, 0x1.99999ap-3f,
+                                              -0x1.554318p-3f, 0x1.1dc5c4p-3f};
+
+// log(2) in double precision.
+static constexpr double LOG2_DOUBLE = 0x1.62e42fefa39efp-1;
+
+// log(2) in float precision.
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log(2), SG, RN);
+static constexpr float LOG2_FLOAT = 0x1.62e43p-1f;
+
 // The function correctly calculates b^x value with at least float precision
 // in a limited range.
 // Range reduction:
@@ -297,73 +320,46 @@ LIBC_INLINE static double log2_eval(double x) {
   return result;
 }
 
-// x should be positive, normal finite value
-LIBC_INLINE static float log_eval_f(float x) {
+template <typename T> LIBC_INLINE static T log_eval(T x) {
   // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx).
-  using FPB = fputil::FPBits<float>;
-  FPB bs(x);
+  using FPBits = fputil::FPBits<T>;
+  FPBits xbits(x);
 
-  float ex = static_cast<float>(bs.get_exponent());
+  T ex = static_cast<T>(xbits.get_exponent());
   // p1 is the leading 7 bits of mx, i.e.
   // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
-  int p1 = static_cast<int>(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7));
+  int p1 = static_cast<int>(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7));
 
   // Set bs to (1 + (mx - p1*2^(-7))
-  bs.set_uintval(bs.uintval() & (FPB::FRACTION_MASK >> 7));
-  bs.set_biased_exponent(FPB::EXP_BIAS);
+  xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7));
+  xbits.set_biased_exponent(FPBits::EXP_BIAS);
   // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
-  float dx = (bs.get_val() - 1.0f) * ONE_OVER_F_FLOAT[p1];
-
-  // Minimax polynomial of log(1 + dx) generated by Sollya with:
-  // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
-  const float COEFFS[6] = {-0x1.fffffep-2f, 0x1.555556p-2f,  -0x1.fffefep-3f,
-                           0x1.99999ap-3f,  -0x1.554318p-3f, 0x1.1dc5c4p-3f};
-
-  float dx2 = dx * dx;
-
-  float c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
-  float c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
-  float c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+  T dx = static_cast<T>(xbits.get_val() - 1.0) *
+         (std::is_same<T, double>::value ? static_cast<T>(ONE_OVER_F[p1])
+                                         : ONE_OVER_F_FLOAT[p1]);
 
-  float p = fputil::polyeval(dx2, dx, c1, c2, c3);
+  T dx2 = dx * dx;
 
-  float result = fputil::multiply_add(ex, 0x1.62e42ep-1f, LOG_F_FLOAT[p1] + p);
-  return result;
-}
-
-// x should be positive, normal finite value
-LIBC_INLINE static double log_eval(double x) {
-  // For x = 2^ex * (1 + mx)
-  //   log(x) = ex * log(2) + log(1 + mx)
-  using FPB = fputil::FPBits<double>;
-  FPB bs(x);
-
-  double ex = static_cast<double>(bs.get_exponent());
-
-  // p1 is the leading 7 bits of mx, i.e.
-  // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
-  int p1 = static_cast<int>(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7));
+  const T *coeffs = nullptr;
+  T log2_val;
+  if constexpr (std::is_same<T, double>::value) {
+    coeffs = LOG_COEFFS_DOUBLE;
+    log2_val = LOG2_DOUBLE;
+  } else {
+    coeffs = LOG_COEFFS_FLOAT;
+    log2_val = LOG2_FLOAT;
+  }
 
-  // Set bs to (1 + (mx - p1*2^(-7))
-  bs.set_uintval(bs.uintval() & (FPB::FRACTION_MASK >> 7));
-  bs.set_biased_exponent(FPB::EXP_BIAS);
-  // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
-  double dx = (bs.get_val() - 1.0) * ONE_OVER_F[p1];
+  T c1 = fputil::multiply_add(dx, coeffs[1], coeffs[0]);
+  T c2 = fputil::multiply_add(dx, coeffs[3], coeffs[2]);
+  T c3 = fputil::multiply_add(dx, coeffs[5], coeffs[4]);
 
-  // Minimax polynomial of log(1 + dx) generated by Sollya with:
-  // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
-  const double COEFFS[6] = {-0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2,
-                            -0x1.ffffffefe562dp-3, 0x1.9999817d3a50fp-3,
-                            -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3};
-  double dx2 = dx * dx;
-  double c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
-  double c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
-  double c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+  T p = fputil::polyeval(dx2, dx, c1, c2, c3);
 
-  double p = fputil::polyeval(dx2, dx, c1, c2, c3);
-  double result =
-      fputil::multiply_add(ex, /*log(2)*/ 0x1.62e42fefa39efp-1, LOG_F[p1] + p);
-  return result;
+  if constexpr (std::is_same<T, double>::value)
+    return fputil::multiply_add(ex, log2_val, LOG_F[p1] + p);
+  else
+    return fputil::multiply_add(ex, log2_val, LOG_F_FLOAT[p1] + p);
 }
 
 // Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We

libc/test/src/math/smoke/atanhf16_test.cpp

@@ -16,6 +16,10 @@ using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
 
 TEST_F(LlvmLibcAtanhf16Test, SpecialNumbers) {
   LIBC_NAMESPACE::libc_errno = 0;
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atanhf16(sNaN),
+                                           FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
   EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atanhf16(aNaN));
   EXPECT_MATH_ERRNO(0);