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

[Muhammad Bassiouni via libc-commits]

================
@@ -0,0 +1,456 @@
+//===-- Implementation header of sqrtf128 ---------------===//
+//
+// 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 "include/llvm-libc-types/float128.h"
+
+#ifdef LIBC_TYPES_HAS_FLOAT128
+
+#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"
+
+// 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 {
+
+namespace sqrtf128_internal {
+
+using FPBits = fputil::FPBits<float128>;
+
+template <typename T, typename U = T>
+LIBC_INLINE static constexpr T prod_hi(T, U);
+
+// Get high part of integer multiplications.
+// Use template to prevent implicit conversion.
+template <>
+LIBC_INLINE constexpr uint64_t prod_hi<uint64_t>(uint64_t x, uint64_t y) {
+  return static_cast<uint64_t>(
+      (static_cast<UInt128>(x) * static_cast<UInt128>(y)) >> 64);
+}
+
+// Get high part of unsigned 128x64 bit multiplication.
+template <>
+LIBC_INLINE constexpr UInt128 prod_hi<UInt128, uint64_t>(UInt128 x,
+                                                         uint64_t y) {
+  uint64_t x_lo = static_cast<uint64_t>(x);
+  uint64_t x_hi = static_cast<uint64_t>(x >> 64);
+  UInt128 xyl = static_cast<UInt128>(x_lo) * static_cast<UInt128>(y);
+  UInt128 xyh = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y);
+  return xyh + (xyl >> 64);
+}
+
+// Get high part of signed 64x64 bit multiplication.
+template <>
+LIBC_INLINE constexpr int64_t prod_hi<int64_t>(int64_t x, int64_t y) {
+  return static_cast<int64_t>(
+      (static_cast<Int128>(x) * static_cast<Int128>(y)) >> 64);
+}
+
+// Get high 128-bit part of unsigned 128x128 bit multiplication.
+template <>
+LIBC_INLINE constexpr UInt128 prod_hi<UInt128>(UInt128 x, UInt128 y) {
+  uint64_t x_lo = static_cast<uint64_t>(x);
+  uint64_t x_hi = static_cast<uint64_t>(x >> 64);
+  uint64_t y_lo = static_cast<uint64_t>(y);
+  uint64_t y_hi = static_cast<uint64_t>(y >> 64);
+
+  UInt128 xh_yh = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y_hi);
+  UInt128 xh_yl = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y_lo);
+  UInt128 xl_yh = static_cast<UInt128>(x_lo) * static_cast<UInt128>(y_hi);
+
+  xh_yh += xh_yl >> 64;
+
+  return xh_yh + (xl_yh >> 64);
+}
+
+// Get high 128-bit part of mixed sign 128x128 bit multiplication.
+template <>
+LIBC_INLINE constexpr Int128 prod_hi<Int128, UInt128>(Int128 x, UInt128 y) {
----------------
bassiounix wrote:

mark as `static`

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