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

[via libc-commits]

================
@@ -297,6 +297,45 @@ 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 FPBits = fputil::FPBits<float>;
+  FPBits xbits(x);
+
+  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));
+
+  // Set bs to (1 + (mx - p1*2^(-7))
+  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 = (xbits.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]);
----------------
lntue wrote:

```
> dirtyinfnorm(log(1 + x)/x - P, [0, 2^-7]);
0x1.13796f285210e7f1dc74fb3d17a812b992a9369b5p-65
```
The coefficients should be in single, replacing `[|D...|]` with `[|SG...|]`.
Also the errors of order `2^-64` is way overkill, and will be outside of `float`'s precision (aka wasteful).  For half precision inputs, the lower bound of `|x|` that needs this part of computations would be ~ 2^-6, so the outputs would be at least `2^-6`, and hence the precision needed would be `~ 2^(-6 - 2*(precision of float16) + some more)`.  So if you use the entire precision of float, errors ~ 2^-30 should be fine.  You can now play around with the range reduction (how big is the table) and the degree of the polynomials.

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