Artimosha - 1 month ago 13

C++ Question

`acos(double)`

`printf("%.30g\n", double(acosl(0.49990774364240564)));`

printf("%.30g\n", acos(0.49990774364240564));

on x64:

`1.0473040763868076`

on x32:

`1.0473040763868078`

on linux4.4 x32 and x64 with sse enabled:

`1.0473040763868078`

is there a way to make VSx64

`acos()`

`1.0473040763868078`

Answer Source

**TL:DR: this is normal and you can't reasonably change it.**

The 32-bit library may be using 80-bit FP values in x87 registers for its temporaries, avoiding rounding off to 64-bit `double`

after every operation. (Unless there's a whole separate library, compiling your own code to use SSE doesn't change what's inside the library, or even the calling convention for passing data to the library. But since 32-bit passes `double`

and `float`

in memory on the stack, a library is free to load it with SSE2 or with x87. Still, you don't get the performance advantage of passing FP values in xmm registers unless it's impossible for non-SSE code to use the library.)

It's also possible that they're different simply because they use a different order of operations, producing different temporaries along the way. That's less plausible, unless they're separately hand-written in asm. If they're built from the same C source (without "unsafe" FP optimizations), then the compiler isn't allowed to reorder things, because of this non-associative behaviour of FP math.

glibc's libm (used on Linux) typically favours precision over speed, so its giving you the correctly-rounded result out to the last bit of the mantissa for both 32 and 64-bit. The IEEE FP standard only requires the basic operations (+ - * / FMA and FP remainder) to be "correctly rounded" out to the last bit of the mantissa. (i.e. rounding error of at most 0.5 ulp). (The exact result, according to `calc`

, is `1.047304076386807714...`

. Keep in mind that `double`

(on x86 with normal compilers) is IEEE754 binary64, so internally the mantissa and exponent are in base2. If you print enough extra decimal digits, though, you can tell that `...7714`

should round up to `...78`

, although really you should print more digits in case they're not zero beyond that. I'm just assuming it's `...78000`

.)

**So Microsoft's 64-bit library implementation produces 1.0473040763868076 and there's pretty much nothing you can do about it, other than not use it**. (e.g. find your own

`acos()`

implementation and use it.) But FP determinism is `acos()`

.You might be able to get the 32-bit library version to produce the same value as the 64-bit version, if it uses x87 and changing the x87 precision setting affects it. But the other way around is not possible: SSE2 has separate instructions for 64-bit `double`

and 32-bit `float`

, and always rounds after every instruction, so you can't change any setting that will increase precision result. (You could change the SSE rounding mode, and that will change the result, but not in a good way!)

See also:

Intermediate Floating-Point Precision and the rest of Bruce Dawson's

**excellent**series of articles about floating point. (table of contents.The linked article describes how some versions of VC++'s CRT runtime startup set the x87 FP register precision to 53-bit mantissa instead of 80-bit full precision. Also that D3D9 will set it to 24, so even

`double`

only has the precision of`float`

if done with x87.- https://en.wikipedia.org/wiki/Rounding#Table-maker.27s_dilemma
- What Every Computer Scientist Should Know About Floating-Point Arithmetic