[libc-commits] [libc] [libc][math] Switch log1pf to use the same log_eval from inverse hyperbolic functions. (PR #188388)

[Mohamed Emad via libc-commits]

================
@@ -35,46 +36,83 @@ LIBC_INLINE constexpr float acoshf(float x) {
     return FPBits_t::quiet_nan().get_val();
   }
 
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
   uint32_t x_u = xbits.uintval();
-  if (LIBC_UNLIKELY(x_u >= 0x4f8ffb03)) {
-    if (LIBC_UNLIKELY(xbits.is_inf_or_nan()))
+  double x_d = static_cast<double>(x);
+
+  if (LIBC_UNLIKELY(x_u >= 0x4580'0000U)) {
+    // x >= 2^12.
+    if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
+      if (xbits.is_signaling_nan()) {
+        fputil::raise_except_if_required(FE_INVALID);
+        return FPBits_t::quiet_nan().get_val();
+      }
       return x;
+    }
 
-    // Helper functions to set results for exceptional cases.
-    auto round_result_slightly_down = [](float r) -> float {
-      volatile float tmp = r;
-      tmp = tmp - 0x1.0p-25f;
-      return tmp;
-    };
-    auto round_result_slightly_up = [](float r) -> float {
-      volatile float tmp = r;
-      tmp = tmp + 0x1.0p-25f;
-      return tmp;
-    };
-
-    switch (x_u) {
-    case 0x4f8ffb03: // x = 0x1.1ff606p32f
-      return round_result_slightly_up(0x1.6fdd34p4f);
-    case 0x5c569e88: // x = 0x1.ad3d1p57f
-      return round_result_slightly_up(0x1.45c146p5f);
-    case 0x5e68984e: // x = 0x1.d1309cp61f
-      return round_result_slightly_up(0x1.5c9442p5f);
-    case 0x655890d3: // x = 0x1.b121a6p75f
-      return round_result_slightly_down(0x1.a9a3f2p5f);
-    case 0x6eb1a8ec: // x = 0x1.6351d8p94f
-      return round_result_slightly_down(0x1.08b512p6f);
-    case 0x7997f30a: // x = 0x1.2fe614p116f
-      return round_result_slightly_up(0x1.451436p6f);
+    // acosh(x) = log(x + sqrt(x^2 - 1))
+    // For large x:
+    //   log(x + sqrt(x^2 - 1)) = log(2x) + log((x + sqrt(x^2 - 1)) / (2x)).
+    // Let U = (x + sqrt(x^2 - 1))/(2x).
+    // Then U = 1 - (x - sqrt(x^2 - 1))/(2x)
+    //        = 1 - (1 - sqrt(1 - 1/x^2))/2
+    //        = 1 - (1/2) * (1/(2x^2) + 1/(8x^4) + ...)
+    //        = 1 - 1/(2x)^2 - 1/(2x)^4 - ...
+    // Hence log(U) = log(1 - 1/(2x)^2 - 1/(2x)^4 - ...)
+    //              = -(1/(2x)^2 - 1/(2x)^4 - ...) -
+    //                - (1/(2x)^2 - 1/(2x)^4 - ...)^2/2 - ...
+    //              ~ -1/(2x)^2 - 1/(2x^4) - ...
+    // For x >= 2^12:
+    //   acosh(x) ~ log(2x) - 1/(2x)^2.
+    // > g = log(2*x) + 1/(4 * x^2);
+    // > dirtyinfnorm((acosh(x) - g)/acosh(x), [2^12, 2^20]);
+    // 0x1.54eb81b0c0df3c9bf68c149748e507fa136e2294fp-55
+    //
+    // For x >= 2^25, 1/(2x)^2 <= 2^-54. So we just need log(2x).
+
+    double y = 2.0 * x_d;
+
+    if (x_u <= 0x4c80'0000U) {
+      // x <= 2^26
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+      if (LIBC_UNLIKELY(x_u == 0x45dc'6414U)) // x = 0x1.b8c828p12f
+        return fputil::round_result_slightly_up(0x1.31bcb6p3f);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+      double y_inv = 0.5 / x_d;
+      fputil::multiply_add(y_inv, -y_inv, log_eval(y));
+
+    } else {
+// x > 2^26
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+      switch (x_u) {
+      case 0x4c803f2c: // x = 0x1.007e58p26f
+        return fputil::round_result_slightly_down(0x1.2b786cp4f);
+      case 0x4f8ffb03: // x = 0x1.1ff606p32f
+        return fputil::round_result_slightly_up(0x1.6fdd34p4f);
+      case 0x5c569e88: // x = 0x1.ad3d1p57f
+        return fputil::round_result_slightly_up(0x1.45c146p5f);
+      case 0x5e68984e: // x = 0x1.d1309cp61f
+        return fputil::round_result_slightly_up(0x1.5c9442p5f);
+      case 0x655890d3: // x = 0x1.b121a6p75f
+        return fputil::round_result_slightly_down(0x1.a9a3f2p5f);
+      case 0x6eb1a8ec: // x = 0x1.6351d8p94f
+        return fputil::round_result_slightly_down(0x1.08b512p6f);
+      case 0x7997f30a: // x = 0x1.2fe614p116f
+        return fputil::round_result_slightly_up(0x1.451436p6f);
+#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+      case 0x65de7ca6: // x = 0x1.bcf94cp76f
+        return fputil::round_result_slightly_up(0x1.af66cp5f);
+      case 0x7967ec37: // x = 0x1.cfd86ep115f
+        return fputil::round_result_slightly_up(0x1.43ff6ep6f);
+#endif // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
----------------
hulxv wrote:

```suggestion
#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
      case 0x65de7ca6: // x = 0x1.bcf94cp76f
        return fputil::round_result_slightly_up(0x1.af66cp5f);
      case 0x7967ec37: // x = 0x1.cfd86ep115f
        return fputil::round_result_slightly_up(0x1.43ff6ep6f);
      case 0x4584ad64U" // x = 0x1.095ac8p12f
        return fputil::round_result_slightly_down(0x1.217f68p3f);
#endif // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
```

I think it fails on non-FMA machines because it fails on the CI. I am not sure if the CI machines are non-FMA or not, but the test is passing on my local machine.

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