3

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?

George Koehler
  • 1,560
  • 17
  • 23
Muhammad
  • 1,598
  • 3
  • 19
  • 31

2 Answers2

8

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

sepp2k
  • 363,768
  • 54
  • 674
  • 675
  • thank you, so i can take the lower 32bit calculate the square root for it - call it (a) - then subtract the lower 32bit from the 64bit number and get the square root for that number - call it (b) - now all i need to do is multiply (a) by (b) and we're done. – Muhammad May 30 '12 at 12:51
  • 1
    @Muhammadalaa: Eh, no. The 64 bits can be written mathematically as (H * 2^32 + L) where H and L are 32 bits each. You state that `sqrt(H * 2^32 + L)` = `sqrt(L) * sqrt(H*2^32)`. That's not the case. However, you _can_ calculate it as `sqrt(H) * sqrt(2^32+L/H)`. – MSalters May 30 '12 at 21:00
  • Rounding would give you an error of about 2^-53, but the error here is a factor of 2: `~0u == UINT_MAX` versus `0x80000000 = UINT_MAX/2`. I don't think this explains it. – MSalters May 30 '12 at 21:04
  • @MSalters What? The error here is 1. He expected ~0u = 0xFFFFFFFF. He got 0x100000000 = 0xFFFFFFFF + 1. – sepp2k May 31 '12 at 05:03
  • @sepp2k: Sorry, I got confused again (not that familiar with x86 assembly) so the correct result should in fact be `0x00000000FFFFFFFF` ? As I originally read it, the result was taken from `EAX` and couldn't have been `0x100000000`. – MSalters May 31 '12 at 08:08
  • 1
    @MSalters Now that you mention it, his code can't return 0x100000000 because that doesn't fit in an unsigned it... That's what he said he got though (and it's what he would get if he returned an unsigned long long and set the registers accordingly). Probably he didn't look at the return value (which would be 0), but he used the debugger to find out that the value of `value` is 0x100000000 before the function returns. Or he did look at the return value and then used the debugger to find out *why* it was 0 and then found that `value` was 0x100000000. – sepp2k May 31 '12 at 08:23
  • @sepp2k: you're right, i know that the square root for the largest unsigned 64bit can fit into unsigned 32bit but when i saw `EAX` with zero i debugged the code and i find out the problem that you already explained. – Muhammad May 31 '12 at 11:28
0

You're not very clear about the C++ function that you're calling. sqrt is an overloaded name. You probably wanted sqrt(double(~0ull)). There's no sqrt overload which takes an unsigned long long.

MSalters
  • 173,980
  • 10
  • 155
  • 350