[libc-commits] [PATCH] D118157: [libc] Improve hypotf performance with different algorithm correctly rounded to all rounding modes.

[Tue Ly via Phabricator via libc-commits]

lntue created this revision.
lntue added reviewers: sivachandra, michaelrj, zimmermann6, cqlauter.
Herald added subscribers: ecnelises, tschuett.
Herald added a project: libc-project.
lntue requested review of this revision.

Algorithm for hypotf: compute (a*a + b*b) in double precision, then use Dekker's algorithm to find the rounding error, and then correcting it after taking its square-root.


Repository:
  rG LLVM Github Monorepo

https://reviews.llvm.org/D118157

Files:
  libc/src/math/generic/hypotf.cpp


Index: libc/src/math/generic/hypotf.cpp
===================================================================
--- libc/src/math/generic/hypotf.cpp
+++ libc/src/math/generic/hypotf.cpp
@@ -6,13 +6,57 @@
 //
 //===----------------------------------------------------------------------===//
 #include "src/math/hypotf.h"
-#include "src/__support/FPUtil/Hypot.h"
+#include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/common.h"
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, hypotf, (float x, float y)) {
-  return __llvm_libc::fputil::hypot(x, y);
+  using DoubleBits = fputil::FPBits<double>;
+  using FPBits = fputil::FPBits<float>;
+
+  double xd = static_cast<double>(x);
+  double yd = static_cast<double>(y);
+
+  // These squares are exact.
+  double xSq = xd * xd;
+  double ySq = yd * yd;
+
+  // Compute the sum of squares.
+  double sumSq = xSq + ySq;
+
+  // Compute the rounding error with Dekker's algorithm.
+  double err = ((sumSq - xSq) - ySq) + ((sumSq - ySq) - xSq);
+
+  // Take sqrt in double precision.
+  DoubleBits result(__builtin_sqrt(sumSq));
+
+  if (!DoubleBits(sumSq).is_inf_or_nan()) {
+    // Correct rounding.
+    if (err != 0.0) {
+      double rSq = static_cast<double>(result) * static_cast<double>(result);
+      double diff = sumSq - rSq;
+      constexpr uint64_t mask = 0x0000'0000'3FFF'FFFFULL;
+      uint64_t lrs = result.uintval() & mask;
+
+      if (lrs == 0x0000'0000'1000'0000ULL && err < diff) {
+        result.bits |= 1ULL;
+      } else if (lrs == 0x0000'0000'3000'0000ULL && err > diff) {
+        result.bits -= 1ULL;
+      }
+    }
+  } else {
+    FPBits bits_x(x), bits_y(y);
+    if (bits_x.is_inf_or_nan() || bits_y.is_inf_or_nan()) {
+      if (bits_x.is_inf() || bits_y.is_inf())
+        return static_cast<float>(FPBits::inf());
+      if (bits_x.is_nan())
+        return x;
+      return y;
+    }
+  }
+
+  return static_cast<float>(static_cast<double>(result));
 }
 
 } // namespace __llvm_libc


-------------- next part --------------
A non-text attachment was scrubbed...
Name: D118157.402936.patch
Type: text/x-patch
Size: 1981 bytes
Desc: not available
URL: <http://lists.llvm.org/pipermail/libc-commits/attachments/20220125/9a378ded/attachment-0001.bin>