[libc][math][c23] Fix some test failures.

Author: harrisonGPU

libc/src/math/generic/acoshf16.cpp

@@ -23,10 +23,18 @@ namespace LIBC_NAMESPACE_DECL {
 LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
   using FPBits = fputil::FPBits<float16>;
   FPBits xbits(x);
-
   uint16_t x_u = xbits.uintval();
   uint16_t x_abs = x_u & 0x7fff;
 
+//   if (LIBC_UNLIKELY(x <= 1.0f)) {
+//     if (x == 1.0f)
+//       return 0.0f;
+//     // x < 1.
+//     fputil::set_errno_if_required(EDOM);
+//     fputil::raise_except_if_required(FE_INVALID);
+//     return FPBits::quiet_nan().get_val();
+//   }
+  
   // Check for NaN input first.
   if (LIBC_UNLIKELY(xbits.is_nan())) {
     if (xbits.is_signaling_nan()) {
@@ -62,11 +70,18 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
   // High precision for inputs very close to 1.0
   if (LIBC_UNLIKELY(xf32 < 1.25f)) {
     float delta = xf32 - 1.0f;
-    float sqrt_2_delta = fputil::sqrt<float>(2.0f * delta);
-    double correction = (double)delta / 12.0 - (3.0 * (double)delta * delta) / 160.0;
-    double precise_result = (double)sqrt_2_delta * (1.0 - correction);
-    return fputil::cast<float16>(static_cast<float>(precise_result));
+    float sqrt_2 = fputil::sqrt<float>(2.0f * delta);
+    float sqrt_2d = fputil::sqrt<float>(2.0f * delta);
+    float d32 = delta * fputil::sqrt<float>(delta);
+    float term2 = d32 / (6.0f * fputil::sqrt<float>(2.0f));
+    float d52 = d32 * delta;
+    float term3 = 3.0f * d52 / (80.0f * sqrt_2);
+    float d72 = d52 * delta;
+    float term4 = 5.0f * d72 / (1792.0f * sqrt_2);
+    float result = sqrt_2d - term2 + term3 - term4;
+    return fputil::cast<float16>(result);
   }
+
   // Special optimization for large input values.
   if (LIBC_UNLIKELY(xf32 >= 32.0f)) {
     float result = static_cast<float>(log_eval(2.0f * xf32));

libc/test/src/math/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 "src/__support/FPUtil/FPBits.h"
@@ -40,34 +41,38 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) {
 }
 
 TEST_F(LlvmLibcAcoshf16Test, InFloat16Range) {
-  constexpr uint16_t START = 0x3C00U; // 1.0
-  constexpr uint16_t STOP = 0x7BFFU;  // Largest finite float16 value
+  constexpr uint32_t COUNT = 100'000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
 
-  for (uint16_t bits = START; bits <= STOP; ++bits) {
-    float16 x = FPBits(bits).get_val();
-    if (FPBits(bits).is_nan() || FPBits(bits).is_inf())
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    LIBC_NAMESPACE::fputil::FPBits<float> bits(v);
+    float xf32 = bits.get_val();
+    if (bits.is_nan() || bits.is_inf())
       continue;
-    ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, float(x),
-                                   float(LIBC_NAMESPACE::acoshf16(x)),
-                                   5.0);
+    if (xf32 < 1.0f)
+      continue;
+    float16 xh = LIBC_NAMESPACE::fputil::cast<float16>(xf32);
+
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, xh,
+                                   LIBC_NAMESPACE::acoshf16(xh), 3.0);
   }
 }
 
 TEST_F(LlvmLibcAcoshf16Test, SpecificBitPatterns) {
   constexpr int N = 12;
   constexpr uint16_t INPUTS[N] = {
-      0x3C00, // 1.0
-      0x3C01, // just above 1.0 (minimally larger than 1)
-      0x3E00, // 1.5
-      0x4200, // 3.0
-      0x4500, // 5.0
-      0x4900, // 10.0
-      0x51FF, // ~47.94 (random mid-range value)
-      0x5CB0, // ~300.0 (random mid-range value)
-      0x643F, // ~1087.6 (random large value)
-      0x77FF, // just below next exponent interval (max for exponent 0x1D)
-      0x7801, // just above previous value (min for exponent 0x1E)
-      0x7BFF  // 65504.0 (max finite half)
+      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();