[libc][math][c23] Add rsqrtf16() function

[amemov]

Addresses #132818
Part of #95250

[github-actions]

⚠️ C/C++ code formatter, clang-format found issues in your code. ⚠️

You can test this locally with the following command:

git-clang-format --diff HEAD~1 HEAD --extensions h,cpp -- libc/src/math/generic/rsqrtf16.cpp libc/src/math/rsqrtf16.h libc/test/src/math/rsqrtf16_test.cpp libc/test/src/math/smoke/rsqrtf16_test.cpp libc/utils/MPFRWrapper/MPCommon.cpp libc/utils/MPFRWrapper/MPCommon.h libc/utils/MPFRWrapper/MPFRUtils.cpp libc/utils/MPFRWrapper/MPFRUtils.h

View the diff from clang-format here.

diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
index 6ad9f5f96..b53331a2d 100644
--- a/libc/src/math/generic/rsqrtf16.cpp
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -70,9 +70,9 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
   int exp_floored = -(exponent >> 1);
 
   if (mantissa == 0.5f) {
-    // When mantissa is 0.5f, x was a power of 2 (or subnormal that normalizes this way).
-    // 1/sqrt(0.5f) = sqrt(2.0f) = 0x1.6a09e6p0f
-    // If exponent is odd (exponent = 2k + 1):
+    // When mantissa is 0.5f, x was a power of 2 (or subnormal that normalizes
+    // this way). 1/sqrt(0.5f) = sqrt(2.0f) = 0x1.6a09e6p0f If exponent is odd
+    // (exponent = 2k + 1):
     //   rsqrt(x) = (1/sqrt(0.5)) * 2^(-(2k+1)/2) = sqrt(2) * 2^(-k-0.5)
     //            = sqrt(2) * 2^(-k) * (1/sqrt(2)) = 2^(-k)
     //   exp_floored = -((2k+1)>>1) = -(k) = -k
@@ -90,13 +90,13 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
   } else {
     // 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);
-    
+    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  
+    // 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;

[amemov]

Trying to figure out what would be the best option to compute the result.
I found that the current polynomial produces the least errors ( bigger ones yield negligible results )
P = fpminimax(1/sqrt(x), [|0,1,2,3,4,5|], [|SG...|], [0.5, 1]);
And has ULP Error = 1.0

Also found this already existing implementation:

llvm-project/libc/src/__support/fixed_point/sqrt.h

Line 39 in ae6b4b2

// P = fpminimax(sqrt(x), 1, [|8, 8|], [i * 2^-4, (i + 1)*2^-4],

It has some other interesting points that I found when I was doing my research: specifically, Newton's method.

Upd: Tried adding 2 iterations of Newton's method. Each significantly reduced number of errors, but there are still some