Added range reduction to the approximation

Author: amemov

libc/src/math/generic/CMakeLists.txt

@@ -968,8 +968,8 @@ add_entrypoint_object(
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.manipulation_functions
     libc.src.__support.FPUtil.polyeval
-    libc.src.__support.FPUtil.sqrt
     libc.src.__support.macros.optimization
     libc.src.__support.macros.properties.types  
 )

libc/src/math/generic/rsqrtf16.cpp

@@ -11,10 +11,10 @@
 #include "hdr/fenv_macros.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/ManipulationFunctions.h"
 #include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/FPUtil/cast.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/FPUtil/multiply_add.h" // to remove
 #include "src/__support/macros/optimization.h"
 
 namespace LIBC_NAMESPACE_DECL {
@@ -40,7 +40,7 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
   if (LIBC_UNLIKELY(x_abs == 0x0)) {
     fputil::raise_except_if_required(FE_DIVBYZERO);
     fputil::set_errno_if_required(ERANGE);
-    return FPBits::quiet_nan().get_val();
+    return FPBits::inf(Sign::POS).get_val();
   }
 
   // -inf <= x < 0
@@ -61,11 +61,36 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
   }
 
   // x is valid, estimate the result
-  // 3-degree polynomial generated using Sollya
-  // P = fpminimax(1/sqrt(x), [|1, 2, 3|], [|SG...|], [0.5, 1]);
+  // Range reduction:
+  // x can be expressed as m*2^e, where e - int exponent and m - mantissa
+  // rsqrtf16(x) = rsqrtf16(m*2^e)
+  // rsqrtf16(m*2^e) = 1/sqrt(m) * 1/sqrt(2^e) = 1/sqrt(m) * 1/2^(e/2)
+  // 1/sqrt(m) * 1/2^(e/2) = 1/sqrt(m) * 2^(-e/2)
+
   float xf = x;
-  float result =
-      fputil::polyeval(xf, 0x1.d42408p2f, -0x1.7cc4fep3f, 0x1.66cb6ap2f);
+  int exponent;
+  float mantissa = fputil::frexp(xf, exponent);
+
+  // 6-degree polynomial generated using Sollya
+  // P = fpminimax(1/sqrt(x), [|0,1,2,3,4,5|], [|SG...|], [0.5, 1]);
+  float interm =
+      fputil::polyeval(mantissa, 0x1.9c81c4p1f, -0x1.e2c57cp2f, 0x1.91e8bp3f,
+                       -0x1.899954p3f, 0x1.9edcp2f, -0x1.6bd93cp0f);
+
+  // Round (-e/2)
+  int exp_floored = -(exponent >> 1);
+
+  // rsqrt(x) = 1/sqrt(mantissa) * 2^(-e/2)
+  // rsqrt(x) = P(mantissa) * 2*(exp_floored)
+  float result = fputil::ldexp(interm, exp_floored);
+
+  // Handle the case where exponent is odd
+  if (exponent & 1) {
+    const float ONE_OVER_SQRT2 =
+        0x1.6a09e667f3bcc908b2fb1366ea957d3e3adec1751p-1f;
+    result *= ONE_OVER_SQRT2;
+  }
+
   return fputil::cast<float16>(result);
 }
 } // namespace LIBC_NAMESPACE_DECL

libc/test/src/math/rsqrtf16_test.cpp

@@ -39,4 +39,4 @@ TEST_F(LlvmLibcRsqrtf16Test, NegativeRange) {
     EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Rsqrt, x,
                                    LIBC_NAMESPACE::rsqrtf16(x), 0.5);
   }
-}
+}

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

@@ -20,7 +20,7 @@ TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
   EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::rsqrtf16(sNaN), FE_INVALID);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(0.0f));
+  EXPECT_FP_EQ(inf, LIBC_NAMESPACE::rsqrtf16(0.0f));
   EXPECT_MATH_ERRNO(ERANGE);
 
   EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::rsqrtf16(1.0f));
@@ -34,4 +34,4 @@ TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
 
   EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(-2.0f));
   EXPECT_MATH_ERRNO(EDOM);
-}
+}