From f09aef7698bbae79a67287a353c63c1ca31055b0 Mon Sep 17 00:00:00 2001 From: Reid Spencer Date: Fri, 2 Mar 2007 04:21:55 +0000 Subject: [PATCH] Use a better algorithm for rounding sqrt results. Change the FIXME about this to a NOTE: because pari/gp results start to get rounded incorrectly after 192 bits of precision. APInt and pari/gp never differ by more than 1, but APInt is more accurate because it does not lose precision after 192 bits as does pari/gp. git-svn-id: https://llvm.org/svn/llvm-project/llvm/trunk@34834 91177308-0d34-0410-b5e6-96231b3b80d8 --- lib/Support/APInt.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/lib/Support/APInt.cpp b/lib/Support/APInt.cpp index 5ab16448b77..5a6bb66247c 100644 --- a/lib/Support/APInt.cpp +++ b/lib/Support/APInt.cpp @@ -1239,19 +1239,23 @@ APInt APInt::sqrt() const { } // Make sure we return the closest approximation - // FIXME: This still has an off-by-one error in it. Test case: - // 190 bits: sqrt(694114394047834196220892040454508646882614255319893124270) = - // 26346050824513229049493703285 (not 26346050824513229049493703284) + // NOTE: The rounding calculation below is correct. It will produce an + // off-by-one discrepancy with results from pari/gp. That discrepancy has been + // determined to be a rounding issue with pari/gp as it begins to use a + // floating point representation after 192 bits. There are no discrepancies + // between this algorithm and pari/gp for bit widths < 192 bits. APInt square(x_old * x_old); APInt nextSquare((x_old + 1) * (x_old +1)); if (this->ult(square)) return x_old; - else if (this->ule(nextSquare)) - if ((nextSquare - *this).ult(*this - square)) - return x_old + 1; - else + else if (this->ule(nextSquare)) { + APInt midpoint((nextSquare - square).udiv(two)); + APInt offset(*this - square); + if (offset.ult(midpoint)) return x_old; - else + else + return x_old + 1; + } else assert(0 && "Error in APInt::sqrt computation"); return x_old + 1; }