Wrong double initialization on Windows/MSYS2 compared to Linux

Question

This simple program gives a wrong result on Windows/MSYS2. On Ubuntu 18.04, it works properly:

#include <math.h>
#include <stdio.h>

int main(int argc, char *argv[]) {
  double a = -1.5708;
  double c = cos(-1.5708);
  double d = cos(a);
  unsigned long long *lc = (unsigned long long *)&c;
  unsigned long long *ld = (unsigned long long *)&d;
  printf("c = %.17e %llu\n", c, *lc);
  printf("d = %.17e %llu\n", d, *ld);
  if (c == d)
    printf("correct result\n");
  else
    printf("wrong result\n");
  return 0;
}

Windows/MSYS2:

$ gcc cos.c -o cos -lm ; ./cos
c = -3.67320510334657393e-06 13749155573439758092
d = -3.67320510334657563e-06 13749155573439758096
wrong result

Ubuntu 18.04:

olivier@ubuntu:~$ gcc cos.c -o cos -lm ; ./cos
c = -3.67320510334657393e-06 13749155573439758092
d = -3.67320510334657393e-06 13749155573439758092
correct result

I am doing something wrong, is it a bug of gcc on MSYS2 or something else?

Maybe related: https://stackoverflow.com/questions/16395615/is-there-any-way-to-make-sure-the-floating-point-arithmetic-result-the-same-in-b https://stackoverflow.com/questions/15062728/gcc-and-sin-cos-transcendental-functions-precision-like-in-windows https://stackoverflow.com/questions/21183477/windows-vs-linux-math-result-difference — KamilCuk, Jun 03 '21 at 07:43
Confirmed on Win7 with MinGW using `gcc -msse2 -mfpmath=sse ...` provides the correct output. — David C. Rankin, Jun 03 '21 at 09:22
Using `gcc cos.c -msse2 -mfpmath=sse -o cos -lm` did work for me on MSYS2 with gcc 64 bits. — Olivier Michel, Jun 04 '21 at 10:08

Andreas H. · Accepted Answer · 2021-06-03T07:51:19.880

It cannot be expected to get bitwise identical results for floating point functions other than sqrt across systems.

More specifically cos is not guarantueed to be correctly rounded so the result of the run-time library (your second result d) is not necessarily identical to the result of the constant folding happending during compilation (your result c).

For instance the page on GCC FP Math states that during constant folding all operations are correctly rounded (even the transcendental functions), but this is not generally the case for the run-time library due to efficiency reasons. (Also not all compilers would give you this correctly rounded guarantuee for constant folding, so for cross compilers or compilers which use different libm than the generated executable you may see this issue.)

GCC internally uses the MPFR library which implements correctly rounded operations with intermediate arbitrary precision arithmetic, but in the run-time library this would be too slow.

Even for linux you probably can find aguments to the cos function so that there is disagreement between your c and d (unless they are using crlibm under the hood). So the identical of cand d results for linux are coincidence in my opionion.

If you look for more predictable FP math results during run-time you can check out crlibm, which is an efficient implementation of many floating point functions that give correctly rounded results. It is much more efficient than MPFR but not as efficient as the standard libm functions.

Wrong double initialization on Windows/MSYS2 compared to Linux

1 Answers1