[libc] Update code.

Author: harrisonGPU

libc/src/math/generic/atanhf16.cpp

@@ -98,7 +98,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
   }
 
   float xf = x;
-  return fputil::cast<float16>(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f)));
+  return fputil::cast<float16>(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f)));
 }
 
 } // namespace LIBC_NAMESPACE_DECL

libc/src/math/generic/explogxf.h

@@ -133,29 +133,6 @@ 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:
@@ -321,12 +298,12 @@ LIBC_INLINE static double log2_eval(double x) {
 }
 
 // x should be positive, normal finite value
-template <typename T> LIBC_INLINE static T log_eval(T x) {
+LIBC_INLINE static float log_eval_f(float x) {
   // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx).
-  using FPBits = fputil::FPBits<T>;
+  using FPBits = fputil::FPBits<float>;
   FPBits xbits(x);
 
-  T ex = static_cast<T>(xbits.get_exponent());
+  float ex = static_cast<float>(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>(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7));
@@ -335,32 +312,63 @@ template <typename T> LIBC_INLINE static T log_eval(T x) {
   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)).
-  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 dx = (xbits.get_val() - 1.0f) * ONE_OVER_F_FLOAT[p1];
 
-  T dx2 = dx * dx;
+  // 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};
 
-  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;
-  }
+  float dx2 = dx * dx;
 
-  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]);
+  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 p = fputil::polyeval(dx2, dx, c1, c2, c3);
+  float p = fputil::polyeval(dx2, dx, c1, c2, c3);
 
-  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);
+  // Generated by Sollya with the following commands:
+  //   > display = hexadecimal;
+  //   > round(log(2), SG, RN);
+  static constexpr float LOGF_2 = 0x1.62e43p-1f;
+
+  float result = fputil::multiply_add(ex, LOGF_2, 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));
+
+  // 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];
+
+  // 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]);
+
+  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;
 }
 
 // Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We