[libc][math][c23] Reduce MPFR to 1.

Author: harrisonGPU

libc/src/math/generic/acoshf16.cpp

@@ -7,34 +7,28 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/acoshf16.h"
-#include "hdr/errno_macros.h"
-#include "hdr/fenv_macros.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/generic/sqrt.h"
 #include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/FPUtil/sqrt.h"
 #include "src/__support/macros/optimization.h"
-#include "src/math/generic/common_constants.h"
 #include "src/math/generic/explogxf.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
+static constexpr size_t N_EXCEPTS = 1;
+static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSHF16_EXCEPTS{
+    {// (input, RZ output, RU offset, RD offset, RN offset)
+     {0x41B7, 0x3ED8, 0, 1, 0}}};
+
 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()) {
@@ -68,25 +62,20 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
   float xf32 = x;
 
   // High precision for inputs very close to 1.0
-  // if (LIBC_UNLIKELY(xf32 < 1.25f)) {
-  //   float delta = xf32 - 1.0f;
-  //   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);
-  // }
+  if (LIBC_UNLIKELY(xf32 < 1.25f)) {
+    float delta = xf32 - 1.0f;
+    float sqrt_2_delta = fputil::sqrt<float>(2.0 * delta);
+    float x2 = delta;
+    float pe = fputil::polyeval(x2, 0x1.0000000000000p+0f,
+                                -0x1.55551a83a9472p-4f, 0x1.331601c4b8ecfp-6f,
+                                -0x1.6890f49eb0acbp-8f, 0x1.8f3a617040a6ap-10f);
+    float approx = sqrt_2_delta * pe;
+    return fputil::cast<float16>(approx);
+  }
 
-  // Special optimization for large input values.
-  // if (LIBC_UNLIKELY(xf32 >= 32.0f)) {
-  //   float result = static_cast<float>(log_eval(2.0f * xf32));
-  //   return fputil::cast<float16>(result);
-  // }
+  if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval());
+      LIBC_UNLIKELY(r.has_value()))
+    return r.value();
 
   // Standard computation for general case.
   float sqrt_term =

libc/test/src/math/acoshf16_test.cpp

@@ -48,7 +48,7 @@ TEST_F(LlvmLibcAcoshf16Test, PositiveRange) {
     float16 x = FPBits(v).get_val();
 
     EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x,
-                                   LIBC_NAMESPACE::acoshf16(x), 3.0);
+                                   LIBC_NAMESPACE::acoshf16(x), 1.0);
   }
 }