Computing the square root of a 64-bit integer

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 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

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



• Casting uint64_t → double and std::sqrt are lossy but non-decreasing operations, so we have




sqrt(r*r)

• 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

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

• 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);
}

StackExchange Code Review Q#69069, answer score: 6