[libc][math][c23] Do some change.

Author: harrisonGPU

libc/src/math/generic/acoshf16.cpp

@@ -1,4 +1,4 @@
-//===-- Half-precision acoshf16 function ----------------------------------===//
+//===-- Half-precision acosh(x) function ----------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/acoshf16.h"
+#include "explogxf.h"
 #include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
@@ -16,16 +17,15 @@
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/sqrt.h"
 #include "src/__support/macros/optimization.h"
-#include "src/math/generic/explogxf.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
 static constexpr size_t N_EXCEPTS = 2;
 static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSHF16_EXCEPTS{{
     // (input, RZ output, RU offset, RD offset, RN offset)
-    // x = 0x1.6ep+1, acoshf16(x) = 0x1.dbp+0 (RZ)
+    // x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ)
     {0x41B7, 0x3ED8, 1, 0, 0},
-    // x = 0x1.c8p+0, acoshf16(x) = 0x1.27cp-1 (RZ)
+    // x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ)
     {0x3CE4, 0x393E, 1, 0, 1},
 }};
 
@@ -70,42 +70,24 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
     return r.value();
 
   float xf = x;
-  // 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:
+  // High-precision polynomial approximation for inputs close to 1.0
+  // ([1, 1.25)).
   //
-  // Step-by-step derivation:
-  // 1. Define y = acosh(x), thus by definition x = cosh(y).
-  //
-  // 2. Expand cosh(y) using exponential identities:
-  //      cosh(y) = (e^y + e^{-y}) / 2
-  //      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 explicitly in terms of delta (for small delta):
-  //      y = acosh(1 + delta) ≈ sqrt(2 * delta) for very small delta.
-  //
-  // 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].
-  //
-  // 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.
+  // Brief derivation:
+  // 1. Start from the definition: acosh(x) = y means x = cosh(y).
+  // 2. Set x = 1 + delta. For small delta, y ≈ sqrt(2 * delta).
+  // 3. Expand acosh(1 + delta) using Taylor series around delta=0:
+  //    acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160
+  //                     - 5*delta^3/896 + 35*delta^4/18432 + ...]
+  // 4. Truncate the series to fit accurately for delta in [0, 0.25].
+  // 5. Polynomial coefficients (from sollya) used here are:
+  //    P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2
+  //               - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4
   if (LIBC_UNLIKELY(xf < 1.25f)) {
     float delta = xf - 1.0f;
     float sqrt_2_delta = fputil::sqrt<float>(2.0 * delta);
-    float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.55551ap-4f, 0x1.33160cp-6f,
-                                -0x1.6890f4p-8f, 0x1.8f3a62p-10f);
+    float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f,
+                                -0x1.6db6dcp-8f, 0x1.f1c71cp-10f);
     float approx = sqrt_2_delta * pe;
     return fputil::cast<float16>(approx);
   }

libc/test/src/math/acoshf16_test.cpp

@@ -6,11 +6,10 @@
 //
 //===----------------------------------------------------------------------===//
 
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/acoshf16.h"
 #include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
 #include "utils/MPFRWrapper/MPFRUtils.h"
 #include <stdint.h>
 

libc/test/src/math/smoke/CMakeLists.txt

@@ -3953,7 +3953,8 @@ add_fp_unittest(
     acoshf16_test.cpp
   DEPENDS
     libc.src.errno.errno
-    libc.src.math.acoshf16  
+    libc.src.math.acoshf16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(

libc/test/src/math/smoke/acoshf16_test.cpp

@@ -6,6 +6,7 @@
 //
 //===----------------------------------------------------------------------===//
 
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/acoshf16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -33,9 +34,11 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) {
   EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf));
   EXPECT_MATH_ERRNO(EDOM);
 
-  EXPECT_FP_EQ(float16(0.0f), LIBC_NAMESPACE::acoshf16(float16(1.0f)));
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::acoshf16(
+                         LIBC_NAMESPACE::fputil::cast<float16>(1.0)));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f)));
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(
+                         LIBC_NAMESPACE::fputil::cast<float16>(0.5)));
   EXPECT_MATH_ERRNO(EDOM);
 }