Search code examples
c++assemblyx86standard-libraryx87

is this a bug in sqrt function

I create an application to compute the primes numbers in 64-bit range so when i tried to compute the square root of 64-bit number using sqrt function from math.h i found the answer is not accurate for example when the input is ~0ull the answer should be ~0u but the one I get is 0x100000000 which is not right, so i decided to create my own version using assembly x86 language to see if this is a bug, here is my function:

inline unsigned prime_isqrt(unsigned long long value)
{
    const unsigned one = 1;
    const unsigned two = 2;

    __asm
    {
        test dword ptr [value+4], 0x80000000
        jz ZERO
        mov eax, dword ptr [value]
        mov ecx, dword ptr [value + 4]

        shrd eax, ecx, 1
        shr  ecx, 1
        mov  dword ptr [value],eax 
        mov  dword ptr [value+4],ecx

        fild value
        fimul two
        fiadd one
        jmp REST
ZERO: 
        fild value
REST: 
        fsqrt
        fisttp value
        mov eax, dword ptr [value]
    }
}

the input is an odd number to get its square root. When i test my function with same input the result was the same.

What i don't get is why those functions round the result or to be specific why sqrt instruction round the result?

Solution

  • sqrt doesn't round anything - you do when you convert your integer into a double. A double can't represent all the numbers that a 64-bit integer can without loss of precision. Specifically starting at 253, there are multiple integers that will be represented as the same double value.

    So if you convert an integer above 253 to double, you lose some of the least significant bits, which is why (double)(~0ull) is 18446744073709552000.0, not 18446744073709551615.0 (or to be more precise the latter is actually equal to the former because they represent the same double number).