Computing the square root of a 64-bit integer

Question

When sieving primes up to some 64-bit number n I need to determine an upper bound for the greatest potential factor of n, which is the greatest integer that is not greater than the (true) square root of n.

Computing this as std::sqrt(double(n)) is fraught with pitfalls. The standard IEEE double has only 53 bits in its significand, which means that some rounding is often inevitable. One consequence is that the result cannot safely be converted to uint32_t even though the true result could never exceed 32 bits. For example, sqrt(double(uint64_t(0) - 1)) usually rounds up to \$2^{32}\$ and if that is cast to uint32_t then the result becomes 0.

Currently I'm using a function like this:

uint32_t max_factor32 (uint64_t n)
{
   double r = std::sqrt(double(n));

   if (r < UINT32_MAX)
   {
      uint32_t r32 = uint32_t(r);

      return r32 - (uint64_t(r32) * r32 > n);
   }

   return UINT32_MAX;
}

Obviously, this works only if the value returned by sqrt() is reasonably close to the true result, so that the cast to integer yields either the desired result or a number that is exactly one greater.

I think this code should be reasonably safe, even in the face of floating-point hardware/software operating with unknown rounding modes and at an unknown level of strictness. But is that truly so?

Different rounding modes should not affect the correctness of the function, but what about different compiler settings for optimisation and floating point strictness? It seems to me that normally a compiler should only err on the side of excess precision by not rounding intermediate results to IEEE double, in which case my code should work fine as it is. Can different compiler switches throw a monkey wrench into my scheme or is the code safely portable to worlds and compilers unknown?

P.S.: what I want to avoid is the 'tail wagging the dog' style caveats like "this code works only with a strictly-conforming compiler and with strict floating-point precision enabled". On the other hand, if the function max_factor32() were to mess up quietly then the result of the containing program would likely be wrong, and even crashes are possible in cases like the above example where 0 is returned when the true result would be UINT32_MAX.

Answer 1

Your code is OK for prime sieving.

If the rounding mode is nearest-ties-to-even and std::sqrt is correctly rounded (as required by IEEE 754), your code returns

• \$\lfloor\sqrt{n}\rfloor + 1\$ for inputs in the range 0xfffffffdfffffc00 .. 0xfffffffe00000000, and
• \$\lfloor\sqrt{n}\rfloor\$ for all other inputs.

If you make the if-condition r <= UINT32_MAX instead of <, the code will always return \$\lfloor\sqrt{n}\rfloor\$.
Your version doesn't consider the few cases where the correct result is UINT32_MAX - 1 but sqrt(double(n)) rounds to UINT32_MAX from below.

Here's the proof:

• We want to compute \$r = \lfloor\sqrt{n}\rfloor\$, which we can rewrite to r*r <= n < (r+1)*(r+1).
• Casting uint64_t → double and std::sqrt are lossy but non-decreasing operations, so we have
  sqrt(r*r) <= sqrt(n) <= sqrt((r+1)*(r+1)). Both inequalities are now <=.
  We have to special-case UINT32_MAX because (UINT32_MAX+1)*(UINT32_MAX+1) doesn't fit in uint64_t.
• If the rounding mode is nearest-ties-to-even and std::sqrt is correctly rounded, then
  x = sqrt((uint64_t)x * x) for any uint32_t x. (You can test it exhaustively.)
• So we can rewrite the inequalities to r <= sqrt(n) <= r+1 and check only these two cases.

Some considerations:

• Directed rounding modes are messier.
  The proof is the same for rounding mode toward +∞. But when rounding toward −∞ or toward 0, we only get sqrt(double((uint64_t)x * x)) ∈ {x-1, x} for x >= 94906267. The final inequality r-1 <= sqrt(double(n)) <= r+1 means that you have to look for errors in both directions.
• If std::sqrt isn't correctly rounded, the final step breaks down and you're better off with an exact integer square root.

Here's my version (only for the nearest-ties-to-even rounding mode). I like to handle the special case first.

uint32_t max_factor32(uint64_t n)
{
    // handle special case before doing the sqrt
    if (n >= uint64_t(UINT32_MAX) * UINT32_MAX) return UINT32_MAX;

    // now it fits in 32 bits
    uint32_t r = std::sqrt(n);
    return r - (uint64_t(r) * r > n);
}

Answer 2

Usually you can just test n < r*r instead which also has the benefit of being much cheaper than a sqrt. You do need to handle potential integer overflow for n > 2⁶⁴ - 2*2³² + 1.

You can't hoist r*r out of the loop the way you can hoist sqrt(n) out but integer multiplication is cheap on any modern processor.

Answer 3

In order to get a fix on the issue I dug into scripture (C/C++ standards, math textbook) and ran extensive experiments.

The standards mandate that the result of sqrt(n) be within half an ulp of the true result. Since the subsequent rounding to uint32_t sheds 'excess' precision, the final result must be within one ulp of the true value (depending on rounding mode).

However, the square root is exact in this sense only for the double \$n' = n + h\$, which is the result of the initial conversion of the uint64_t n to a double n' with only 53 bits of precision.

The delta h is bounded by the ulp of n, which corresponds to \$(63 - 53 + 1) = 11\$ bits for values between \$2^{63}\$ and \$2^{64}-1\$. The estimate for the absolute error of \$\sqrt{n + h}\$ is \$\frac h {2 \sqrt n}\$; the ulp is constant over that range and equal to \$2^{11}\$. This maximises the absolute error for values close to the beginning of the range, e.g. \$n = 2^{63}\$. The value there is \$\frac{2^{11}}{(2 \sqrt{2^{63}})}\$ or \$\frac{2^{-21}}{\sqrt 2}\$. With other words, the error is negligible.

Since there are two precision-change roundings (casts) and one rounding computation involved, the absolute error for the final result can be greater than one even for a standards-conforming compiler operating in a fully standards-conforming mode.

When testing uncovered the existence of cases where my original max_factor32() diverged from the ideal result, I switched to the following 'luxury' version that should be exact under any and all circumstances, regardless of floating-strictness, optimisation settings or whatnot:

template<typename word_t>
word_t square (word_t n)  {  return word_t(n * n);  }

uint32_t max_factor32_b (uint64_t n)
{
   uint64_t r = uint64_t(std::sqrt(double(n)));

   while (r < UINT32_MAX && square(r) < n)  
      ++r;            
   while (r > UINT32_MAX || square(r) > n)  
      --r;

   return uint32_t(r);
}

That code can be optimised by hoisting a test out of a loop and similar operations, but I left it as is, to preserve the strange antisymmetry between the two cases.

Tests with trillions of computations in near the upper end of the range uncovered only 1025 cases where the two versions of the function returned different results, and in all these cases the simpler code reported a value that was higher by one than the ideal result (\$2^{32}-1\$ instead of \$2^{31}-2\$). For the Sieve of Eratosthenes that makes virtually no difference but other code might not be so forgiving, and other compilers might not be as well-behaved as gcc 4.8.1 and VC++ 2013...

The outlier results were for values of n between 0xFFFFFFFDFFFFFC00 and 0xFFFFFFFE00000000.

Answer 4

I like DarthGizka's integer square root function with error correction, so I made a generic version of it. The code requires a C++14 compiler so It currently only works with clang.

template <typename T>
constexpr T sqrt_helper(T x, T lo, T hi)
{
  if (lo == hi)
    return lo;

  const T mid = (lo + hi + 1) / 2;
  if (x / mid < mid)
    return sqrt_helper<T>(x, lo, mid - 1);
  else
    return sqrt_helper(x, mid, hi);
}

template <typename T>
constexpr T ct_sqrt(T x)
{
  return sqrt_helper<T>(x, 0, x / 2 + 1);
}

template <typename T>
T isqrt(T x)
{
  T r = (T) std::sqrt((double) x);
  T sqrt_max = ct_sqrt(std::numeric_limits<T>::max());

  while (r < sqrt_max && r * r < x)
    r++;
  while (r > sqrt_max || r * r > x)
    r--;

  return r;
}

The above version uses two correction iterations for most numbers. As DarthGizka mentioned it is possible to optimize this, below is a version that uses no correction iteration for most numbers. Depending on your compiler and CPU this version can be faster.

template <typename T>
T isqrt(T x)
{
  T r = (T) std::sqrt((double) x);
  r = std::min(r, ct_sqrt(std::numeric_limits<T>::max()));

  while (r * r > x)
    r--;
  while (x - r * r > r * 2)
    r++;

  return r;
}

Answer 5

If you have C99, I believe you can use sqrtl(), which will return a long double with the answer you need.

However, given your case, where all you really want is the integer you should run your sieve till, I'd recommend taking a look at the solutions provided here.

Answer 6

An alterative to risking floating point issues is to perform an integer based square root:

C based code

uint32_t usqrt64(uint64_t n) {
  if (n < 2) {
    return (uint32_t)n;
  }
  uint32_t s = usqrt64(n >> 2) << 1;
  uint32_t l = s + 1;
  if ((uint64_t)l * l > n) {
    return s;
  }
  return l;
}

Or perhaps a non-recursive version

uint32_t usqrt64(uint64_t t) {
  uint64_t s, b;
  for (b = 0, s = t; s;  b++, s >>= 1) {
    ;
  }

  s = (uint64_t)1 << (b >> 1);
  if (b & 1) {
    s += s >> 1;
  }

  do {
    b = t / s;
    s = (s + b) >> 1;
  } while (b < s);

  return (uint32_t)s;
}

or (adapted from Integer_square_root)

// Binary search
uint32_t usqrt(uint64_t y) {
    if (y == UINT64_MAX) {
      return UINT32_MAX;
    }
    uint64_t L = 0;
    uint64_t R = y + 1;

    while (L + 1 != R) {
      uint64_t M = (L + R) / 2;
      if (M * M <= y) {
        L = M;
      } else {
        R = M;
      }
    }

    return (uint32_t) L;
}