[libc-commits] [libc] [libc][math][c23] Add atanhf16 C23 math function. (PR #132612)

[via libc-commits]

================
@@ -297,6 +297,40 @@ 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) {
+  // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx).
+  using FPB = fputil::FPBits<float>;
+  FPB bs(x);
+
+  float ex = static_cast<float>(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)).
+  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]);
+
+  float p = fputil::polyeval(dx2, dx, c1, c2, c3);
+
+  float result = fputil::multiply_add(ex, 0x1.62e42ep-1f, LOG_F_FLOAT[p1] + p);
----------------
overmighty wrote:

Using a named constant for `0x1.62e42ep-1f` would increase readability. libc/src/math/generic/expxf16.h defines `LOGF_2`:

https://github.com/llvm/llvm-project/blob/595cc960b6da7f55d64d2025abb67a9aabb01a85/libc/src/math/generic/expxf16.h#L183-L186

But explogxf.h probably shouldn't include expxf16.h. I'm not sure which file we could move `LOGF_2`'s definition to.

https://github.com/llvm/llvm-project/pull/132612