[libc-commits] [libc] [llvm] [libc][math] Refactor sqrtf128 to header only (PR #177760)

[Muhammad Bassiouni via libc-commits]

================
@@ -0,0 +1,446 @@
+//===-- Implementation header of sqrtf128 function
+//-------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SQRTF128_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SQRTF128_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/uint128.h"
+#include "src/math/sqrtf128.h"
+
+// Compute sqrtf128 with correct rounding for all rounding modes using integer
+// arithmetic by Alexei Sibidanov (sibid at uvic.ca):
+//   https://github.com/sibidanov/llvm-project/tree/as_sqrt_v2
+//   https://github.com/sibidanov/llvm-project/tree/as_sqrt_v3
+// TODO: Update the reference once Alexei's implementation is in the CORE-MATH
+// project. https://github.com/llvm/llvm-project/issues/126794
+
+// Let the input be expressed as x = 2^e * m_x,
+// - Step 1: Range reduction
+//   Let x_reduced = 2^(e % 2) * m_x,
+//   Then sqrt(x) = 2^(e / 2) * sqrt(x_reduced), with
+//     1 <= x_reduced < 4.
+// - Step 2: Polynomial approximation
+//   Approximate 1/sqrt(x_reduced) using polynomial approximation with the
+//   result errors bounded by:
+//     |r0 - 1/sqrt(x_reduced)| < 2^-32.
+//   The computations are done in uint64_t.
+// - Step 3: First Newton iteration
+//   Let the scaled error defined by:
+//     h0 = r0^2 * x_reduced - 1.
+//   Then we compute the first Newton iteration:
+//     r1 = r0 - r0 * h0 / 2.
+//   The result is then bounded by:
+//     |r1 - 1 / sqrt(x_reduced)| < 2^-62.
+// - Step 4: Second Newton iteration
+//   We calculate the scaled error from Step 3:
+//     h1 = r1^2 * x_reduced - 1.
+//   Then the second Newton iteration is computed by:
+//     r2 = x_reduced * (r1 - r1 * h0 / 2)
+//        ~ x_reduced * (1/sqrt(x_reduced)) = sqrt(x_reduced)
+// - Step 5: Perform rounding test and correction if needed.
+//     Rounding correction is done by computing the exact rounding errors:
+//       x_reduced - r2^2.
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+using FPBits = fputil::FPBits<float128>;
----------------
bassiounix wrote:

```suggestion
```

https://github.com/llvm/llvm-project/pull/177760