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
.