[libc-commits] [libc] [llvm] [libc][math] Implement double-precision acosh (PR #199953)

[via libc-commits]

================
@@ -112,50 +38,52 @@ LIBC_INLINE double acosh(double x) {
     return FPBits::quiet_nan().get_val();
   }
 
-  if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
-    if (xbits.is_signaling_nan()) {
-      fputil::raise_except_if_required(FE_INVALID);
-      return FPBits::quiet_nan().get_val();
-    }
-    return x;
-  }
-
-  // For x >= 2^52, the dropped term 1/(4x^2) is far below 0.5 ULP of
-  // acosh(x) = log(2x), and x^2 would overflow exact_mult for x > ~2^511.
-  // Redirect through math::log, which performs its own Ziv test.
+  // For x >= 2^52, the argument to log simplifies to 2x:
+  //   acosh(x) = log(x + sqrt(x^2 - 1))
+  //   sqrt(x^2 - 1) = x * sqrt(1 - 1/x^2)
+  //                 = x * (1 - 1/(2x^2) - 1/(8x^4) - ...)
+  //   x + sqrt(x^2 - 1) = 2x * (1 - 1/(4x^2) - 1/(16x^4) - ...)
+  //   acosh(x) = log(2x) + log(1 - 1/(4x^2) - ...) = log(2x) - 1/(4x^2) - ...
+  // At x = 2^52 the dropped term 1/(4x^2) < 2^-106, far below
+  // 0.5 * ULP(log(2^52)) ~ 2^-47, so log(2x) is correctly rounded for acosh.
+  // For x^2 > ~2^1022 exact_mult(x,x) would overflow, making this redirect
+  // necessary for all x >= 2^52. NaN and +inf also have uintval >= 0x4330...
   if (LIBC_UNLIKELY(xbits.uintval() >= 0x4330'0000'0000'0000ULL)) {
-    using namespace common_constants_internal;
-    // For x with biased exponent 2046 (x >= 2^1023), 2*x overflows; compute
-    // log(2x) = log(x/2) + 2*log(2) via compensated addition instead.
-    if (LIBC_UNLIKELY(xbits.uintval() >= 0x7FE0'0000'0000'0000ULL)) {
-      double log_xhalf = math::log(x * 0.5);
-      return (log_xhalf + 2.0 * LOG_2_HI) + 2.0 * LOG_2_LO;
+    if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
+      if (xbits.is_signaling_nan()) {
+        fputil::raise_except_if_required(FE_INVALID);
+        return FPBits::quiet_nan().get_val();
+      }
+      return x;
     }
-    return math::log(2.0 * x);
+    // Compute log(2x) = log(x) + log(2) via log_dd_core with e_adj = 1.
+    // Passing x rather than 2x avoids overflow for x >= 2^1023 and ensures
+    // correct rounding: the Ziv test inside log_dd_core runs on the full
+    // double-double sum (including the e_adj * log(2) term) before rounding.
+    return log_dd_core({0.0, x}, 1);
   }
 
-  // acosh(x) = log1p(u),  u = (x-1) + sqrt(x^2-1).
-  // Compute u as a double-double for a correctly-rounded result.
+  // acosh(x) = log(x + sqrt(x^2 - 1)) for x in (1, 2^52).
+  // Compute y = x + sqrt(x^2 - 1) as a double-double, then log(y).
 
-  // x^2 - 1 as double-double.
+  // x^2 - 1 as a double-double. x > 1 implies x^2 > 1 = |(-1)|, so
+  // Fast2Sum (exact_add) is valid.
   DoubleDouble x_sq = fputil::exact_mult(x, x);
-  DoubleDouble v_dd = fputil::exact_add<false>(x_sq.hi, -1.0);
+  DoubleDouble v_dd = fputil::exact_add(x_sq.hi, -1.0);
   v_dd.lo += x_sq.lo;
 
-  // sqrt(x^2-1) as double-double via one Newton correction.
-  // Use exact_mult for s_hi^2 so the correction is accurate on non-FMA
-  // targets too (fma(s,-s,v) without hardware FMA loses ~1/4 ULP).
+  // sqrt(x^2 - 1) as a double-double via one Newton step.
   double s_hi = fputil::sqrt<double>(v_dd.hi);
+  double s_inv = 0.5 / s_hi;
   DoubleDouble s_sq = fputil::exact_mult(s_hi, s_hi);
-  DoubleDouble r_dd = fputil::exact_add<false>(v_dd.hi, -s_sq.hi);
-  double s_lo = (r_dd.hi + (r_dd.lo - s_sq.lo) + v_dd.lo) / (2.0 * s_hi);
+  double s_lo = ((v_dd.hi - s_sq.hi) - s_sq.lo + v_dd.lo) * s_inv;
----------------
lntue wrote:

put `s_sq.lo` and `v_dd.lo` into parentheses to enforce the order of operations.

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