different behaviour or sqrt when compiled with 64 or 32 bits

Question

I'm using sqrt() function from math library, when I build for 64 bit using -m64 I'm getting correct result but when I build for 32 bit I have very inconsistent behaviour.

For example on 64bit

double dx = 0x1.fffffffffffffp+1023;
sqrt(dx); // => 0x1.fffffffffffffp+511
sqrt(0x1.fffffffffffffp+1023);// => 0x1.fffffffffffffp+511

(which I believe is the correctly rounded result, verified with mpfr)

But on 32 bit same input value it behaves differently.

double dx = 0x1.fffffffffffffp+1023;
sqrt(dx); // => 0x1.0p+512
sqrt(0x1.fffffffffffffp+1023); // => 0x1.fffffffffffffp+511

When the same value passed in a variable I'm getting wrong result. I checked rounding mode before and after each call and all are set to round to nearest. What the reason? I'm using gcc 4.6 on a 64bit machine, and options are -mfpmath=sse and -march=pentium for both x86 nad x64 cases.

Answer 1

You haven't said which compiler or architecure you're using, but assuming gcc on x86 / x86-64 then the difference is likely down to the fact that by default gcc uses 387 floating point instructions on 32 bit x86, whereas it uses SSE instructions on x86-64.

The 387 floating point registers are 80 bits wide, whereas double is 64 bits wide. This means that intermediate results can have higher precision using the 387 instructions, which can result in a slightly different answer after rounding. (The SSE2 instructions operate on packed 64 bit doubles).

There's a few ways to change the way the compiler operates, depending on what you want:

• If you use the -ffloat-store option on x86 builds, the compiler will discard extra precision whenever you store a value in a double variable;
• If you use the -mfpmath=sse options on x86 builds, along with -msse2 or an -march= switch that specifies an SSE2-supporting architecture, the compiler will use SSE instructions for floating point just as on x86-64. The code will only run on CPUs that support SSE2, though (Pentium-M / Pentium 4 and later).
• If you use the -mfpmath=387 option on x86-64 builds, the compiler will use 387 instructions for floating point just as on x86. This isn't recommended, though - the x86-64 ABI specifies that floating point values are passed in SSE registers, so the compiler has to do a lot of shuffling between 387 and SSE registers with this option.

Answer 2

Certain compilers, such as gcc, when they see certain math library functions performed on static literals, actually calculate the value during compile time, where-as with the variable, it's by necessity calculated at run-time. The compile-time value is typically calculated by the compiler using a math library like MPFR, GNU MP, etc., so the results will be more accurate, or at least as accurate as possible from platform-to-platform.