Added Newton-Raphson iterations

-The accuracy improved drastically, but it still fails

Author: amemov

libc/src/math/generic/CMakeLists.txt

@@ -967,7 +967,7 @@ add_entrypoint_object(
     libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
-    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.fma
     libc.src.__support.FPUtil.manipulation_functions
     libc.src.__support.FPUtil.polyeval
     libc.src.__support.macros.optimization

libc/src/math/generic/rsqrtf16.cpp

@@ -10,11 +10,11 @@
 #include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FMA.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" // to remove
 #include "src/__support/macros/optimization.h"
 
 namespace LIBC_NAMESPACE_DECL {
@@ -55,11 +55,6 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
     return fputil::cast<float16>(0.0f);
   }
 
-  // x = 1 => rsqrt(x) = 1
-  if (LIBC_UNLIKELY(x_u == 0x1)) {
-    return fputil::cast<float16>(1.0f);
-  }
-
   // x is valid, estimate the result
   // Range reduction:
   // x can be expressed as m*2^e, where e - int exponent and m - mantissa
@@ -72,23 +67,41 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
   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]);
+  // bigger polynomial doesn't generate better results-> the current one
+  // produces the least number of errors but still errors are presents 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);
 
+  // Apply one Newton-Raphson iteration to refine the approximation of
+  // 1/sqrt(mantissa) y_new = y_old * (1.5 - 0.5 * mantissa * y_old^2) Using
+  // fputil::fma for potential precision benefits in the factor calculation
+  float interm_sq = interm * interm;
+  float factor = fputil::fma<float>(-0.5f * mantissa, interm_sq, 1.5f);
+  float interm_refined = interm * factor; // Final multiplication
+
+  // Apply a second Newton-Raphson iteration
+  // y_new = y_old * (1.5 - 0.5 * mantissa * y_old^2)
+  // y_old is now interm_refined
+  float interm_refined_sq = interm_refined * interm_refined;
+  float factor2 = fputil::fma<float>(-0.5f * mantissa, interm_refined_sq, 1.5f);
+  float interm_refined2 = interm_refined * factor2;
+
   // 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);
+  // float result = fputil::ldexp(interm, exp_floored);
+  float result = fputil::ldexp(interm_refined2, exp_floored);
 
   // Handle the case where exponent is odd
   if (exponent & 1) {
-    const float ONE_OVER_SQRT2 =
-        0x1.6a09e667f3bcc908b2fb1366ea957d3e3adec1751p-1f;
-    result *= ONE_OVER_SQRT2;
+    const float ONE_OVER_SQRT2 = 0x1.6a09e6p-1f;
+    // result *= ONE_OVER_SQRT2;
+    result = fputil::fma<float>(result, ONE_OVER_SQRT2,
+                                0.0f); // Use FMA for multiplication
   }
 
   return fputil::cast<float16>(result);