[libc][math][c23] Update for comments again.

Author: harrisonGPU

libc/src/math/generic/acoshf16.cpp

@@ -70,29 +70,37 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
     return r.value();
 
   float xf = x;
-  // High precision for inputs very close to 1.0
-  // For inputs close to 1 (1 <= x < 1.25), use polynomial approximation:
+  // High precision polynomial approximation for inputs very close to 1.0
+  // Specifically, for inputs within the range [1, 1.25), we employ the
+  // following step-by-step Taylor expansion derivation to maintain numerical
+  // accuracy:
   //
   // Step-by-step derivation:
-  // 1. Let y = acosh(x), thus x = cosh(y).
+  // 1. Define y = acosh(x), thus by definition x = cosh(y).
   //
-  // 2. Rewrite cosh identity using exponential form:
+  // 2. Expand cosh(y) using exponential identities:
   //      cosh(y) = (e^y + e^{-y}) / 2
-  //      For y close to 0, let y = sqrt(2 * delta), thus:
-  //      x = cosh(y) ≈ 1 + delta (since cosh(0) = 1, and delta is small)
-  //      thus delta = x - 1.
+  //      For small y, let us set y ≈ sqrt(2 * delta), thus:
+  //      x ≈ cosh(y) ≈ 1 + delta, for small delta
+  //      hence delta = x - 1.
   //
-  // 3: Express y in terms of delta (for small delta):
-  //      y ≈ sqrt(2 * delta)
+  // 3. Express y explicitly in terms of delta (for small delta):
+  //      y = acosh(1 + delta) ≈ sqrt(2 * delta) for very small delta.
   //
-  // 4: Expand acosh(1 + delta) using a Taylor expansion around delta = 0:
-  //      acosh(1 + delta) ≈ sqrt(2 * delta) * P(delta), where P(delta)
-  //      is a polynomial approximation obtained by fitting the function
-  //      precisely in the interval [0, 0.25].
+  // 4. Use Taylor expansion around delta = 0 to obtain a more accurate
+  // polynomial:
+  //      acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160 -
+  //      5*delta^3/896 + 35*delta^4/18432 + ...] For practical computation and
+  //      precision, truncate and fit the polynomial precisely in the range [0,
+  //      0.25].
   //
-  // Because delta = x - 1 and 0 <= delta < 0.25, the polynomial approximation
-  // remains numerically stable and accurate in this domain, ensuring high
-  // precision.
+  // 5. The implemented polynomial approximation (coefficients obtained from
+  // careful numerical fitting) is:
+  //      P(delta) ≈ 1 - 0x1.55551ap-4 * delta + 0x1.33160cp-6 * delta^2 -
+  //      0x1.6890f4p-8 * delta^3 + 0x1.8f3a62p-10 * delta^4
+  //
+  // Since delta = x - 1, and 0 <= delta < 0.25, this approximation achieves
+  // high precision and numerical stability.
   if (LIBC_UNLIKELY(xf < 1.25f)) {
     float delta = xf - 1.0f;
     float sqrt_2_delta = fputil::sqrt<float>(2.0 * delta);

libc/test/src/math/acoshf16_test.cpp

@@ -20,56 +20,9 @@ namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
 static constexpr uint16_t START = 0x3c00U;
 static constexpr uint16_t STOP = 0x7bffU;
 
-TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) {
-  LIBC_NAMESPACE::libc_errno = 0;
-
-  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(aNaN));
-  EXPECT_MATH_ERRNO(0);
-
-  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(sNaN), FE_INVALID);
-  EXPECT_MATH_ERRNO(0);
-
-  EXPECT_FP_EQ_ALL_ROUNDING(float16(0.0f),
-                            LIBC_NAMESPACE::acoshf16(float16(1.0f)));
-  EXPECT_MATH_ERRNO(0);
-
-  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f)));
-  EXPECT_MATH_ERRNO(EDOM);
-
-  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acoshf16(inf));
-  EXPECT_MATH_ERRNO(0);
-
-  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf));
-  EXPECT_MATH_ERRNO(EDOM);
-}
-
 TEST_F(LlvmLibcAcoshf16Test, PositiveRange) {
   for (uint16_t v = START; v <= STOP; ++v) {
     float16 x = FPBits(v).get_val();
-
-    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x,
-                                   LIBC_NAMESPACE::acoshf16(x), 0.5);
-  }
-}
-
-TEST_F(LlvmLibcAcoshf16Test, SpecificBitPatterns) {
-  constexpr int N = 12;
-  constexpr uint16_t INPUTS[N] = {
-      0x3C00, // x = 1.0
-      0x3C01, // x = just above 1.0 (minimally larger than 1)
-      0x3E00, // x = 1.5
-      0x4200, // x = 3.0
-      0x4500, // x = 5.0
-      0x4900, // x = 10.0
-      0x51FF, // x = ~47.94
-      0x5CB0, // x = ~300.0
-      0x643F, // x = ~1087.6
-      0x77FF, // x = just below next exponent interval (max for exponent 0x1D)
-      0x7801, // x = just above previous value (min for exponent 0x1E)
-      0x7BFF  // x = 65504.0 (max finite half)
-  };
-  for (int i = 0; i < N; ++i) {
-    float16 x = FPBits(INPUTS[i]).get_val();
     EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x,
                                    LIBC_NAMESPACE::acoshf16(x), 0.5);
   }