How can you easily calculate the square root of an unsigned long long in C?

Question

I was looking at another question (here) where someone was looking for a way to get the square root of a 64 bit integer in x86 assembly.

This turns out to be very simple. The solution is to convert to a floating point number, calculate the sqrt and then convert back.

I need to do something very similar in C however when I look into equivalents I'm getting a little stuck. I can only find a sqrt function which takes in doubles. Doubles do not have the precision to store large 64bit integers without introducing significant rounding error.

Is there a common math library that I can use which has a long double sqrt function?

Answer 1

There is no need for long double; the square root can be calculated with double (if it is IEEE-754 64-bit binary). The rounding error in converting a 64-bit integer to double is nearly irrelevant in this problem.

The rounding error is at most one part in 2⁵³. This causes an error in the square root of at most one part in 2⁵⁴. The sqrt itself has a rounding error of less than one part in 2⁵³, due to rounding the mathematical result to the double format. The sum of these errors is tiny; the largest possible square root of a 64-bit integer (rounded to 53 bits) is 2³², so an error of three parts in 2⁵⁴ is less than .00000072.

For a uint64_t x, consider sqrt(x). We know this value is within .00000072 of the exact square root of x, but we do not know its direction. If we adjust it to sqrt(x) - 0x1p-20, then we know we have a value that is less than, but very close to, the square root of x.

Then this code calculates the square root of x, truncated to an integer, provided the operations conform to IEEE 754:

uint64_t y = sqrt(x) - 0x1p-20;
if (2*y < x - y*y)
    ++y;

(2*y < x - y*y is equivalent to (y+1)*(y+1) <= x except that it avoids wrapping the 64-bit integer if y+1 is 2³².)

Answer 2

Function sqrtl(), taking a long double, is part of C99.

Note that your compilation platform does not have to implement long double as 80-bit extended-precision. It is only required to be as wide as double, and Visual Studio implements is as a plain double. GCC and Clang do compile long double to 80-bit extended-precision on Intel processors.

Answer 3

Yes, the standard library has sqrtl() (since C99).

Answer 4

If you only want to calculate sqrt for integers, using divide and conquer should find the result in max 32 iterations:

uint64_t mysqrt (uint64_t a)
{
  uint64_t min=0;
  //uint64_t max=1<<32;
  uint64_t max=((uint64_t) 1) << 32; //chux' bugfix
  while(1)
    {
       if (max <= 1 + min)
         return min;           

       uint64_t sqt = min + (max - min)/2;
       uint64_t sq = sqt*sqt;

       if (sq == a) 
         return sqt;

       if (sq > a)
         max = sqt;
       else
         min = sqt;
    }

Debugging is left as exercise for the reader.