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

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 > 264 - 2*232 + 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.

• I think that's a really good idea, at least if the tests have reasonably low frequency. Trade a sloooow sqrt for a bunch of faster MULs. However, in the current case (sieve) the outer loop accounts for a significant part of the runtime, since most primes step less than a single time on average in a typical sieve window, and the number of prime factors is the same order of magnitude as the sieve window size (203,280,221 primes vs. sieve window between 1e6 and 1e9). In this case hoisting + sqrt is easier on impatient nerves... ;-) – DarthGizka Nov 6 '14 at 17:25
• However, if the loop contains the product anyway then the test is essentially free (since the compiler computes it only once). This doesn't answer the floaty question but it's a valuable optimisation tip that can often avoid the floaty question, up to a certain point. Unfortunately I can give you only one uptick... It's worth at least three. – DarthGizka Nov 6 '14 at 17:33
• Also, there may be a way to make your idea work even for high-frequency loops, by using strength reduction (n += (r << 1) + 1), especially since the term ((r << 1) + 1) occurs naturally in a sieve loop. New ideas, new avenues... Thanks! – DarthGizka Nov 6 '14 at 17:41
• The strength reduction idea is neat! On modern processors a branch is generally more expensive than a multiply, even a highly predictable branch like this, so I would expect the multiply to be effectively free anyway. – mattnewport Nov 6 '14 at 17:46

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

• Very nice indeed! I took a long hard look at (lo + hi + 1) since expressions like that are very delicate, but it is perfectly all right for the parameters it gets from ct_sqrt(). I wish I had a compiler that accepts -std=c++14... – DarthGizka Dec 31 '14 at 14:21
• Have a look at my post on stackoverflow, it contains a C++11 version: stackoverflow.com/a/27709195/363778 – Linoliumz Dec 31 '14 at 14:34

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.

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.

• With many compilers, long double is just the same 64-bit float as double; the standard doesn't mandate that it actually be bigger (i.e. Intel's 80-bit float as in BC++, or a 128-bit float as on some workstation compilers). The topic you mentioned didn't have any solutions that weren't either 'use the sqrt(), Luke' or some slow iterative computation... I'd really like to crack the correctness thing for this particular example, so that I can use what I learnt in other, similar cases. – DarthGizka Nov 6 '14 at 11:56
• Even if a long double has sufficient bits to represent a 64-bit integer exactly, the rounding issue could still occur. Unless I'm mistaken, even a hypothetical sqrtlllll(UINT64_MAX) would still return 2^32 instead of 2^32-1. – DarthGizka Nov 6 '14 at 14:32