Comment author: @xavierleroy

I don't have access to a Windows machine right now to confirm, but I suspect the bug is in Microsoft's C run-time library. Indeed, in OCaml <= 4.00 and without special flags, ocamlopt directly translates calls to the "sqrt" Caml function into calls to the C library "sqrt()" function.

The correct result observed with MinGW64 in bytecode probably comes from GCC compiling the auxiliary routine caml_sqrt_float() to use the SSE2 "sqrtsd" instruction directly, bypassing the standard library.

[ Something is very wrong when a processor instruction implements a math function better / more exactly than a standard library function. In general, it's the other way around. ]

With OCaml <= 4.00 in 32 bits, you can use "ocamlopt -ffast-math", which will translate calls to "sqrt" into the equivalent x87 instruction, hopefully treating -0.0 correctly.

The forthcoming OCaml 4.01 version in 64 bit will also use a processor instruction (sqrtsd) to implement square root, again hopefully better than Microsoft's C library.

Comment author: dpariente

Thank you for these explanations.

However, for information, in 32 bits, with -ffast-math the result is still not the correct one:
$ ocamlopt -ffast-math s.ml -o s.exe ; ./s.exe
sqrt(0.0000000000000000e+000) = 0.0000000000000000e+000

Comment author: @yakobowski

@dillon: all your WinXP 32 bits tests seem to be missing the '-' in the argument of the calls to sqrt. Do you know the reason?

Comment author: dpariente

@boris: Right, the minus sign is systematically missing under XP32b.
I don't know why, and suspected a trouble in the pretty-printing, but maybe it's something else!

Comment author: dpariente

@note9338

Under toplevel:

printf "%.16e " (-2.**(-1022.));;

2.2250738585072014e-308 - : unit = ()

printf "%.16e " (0. -. 2.**(-1022.));;

-2.2250738585072014e-308 - : unit = ()

An issue when parsing the minus?

Comment author: @xavierleroy

Boris is right that you should not trust the printing of floating-point numbers either, since it goes through the C standard library as well. Better print the IEEE-754 binary64 encoding, like this:

open Printf

let print_hex_float oc f = fprintf oc "%Lx" (Int64.bits_of_float f)

let _ =
printf "sqrt(%a) = %a\n" print_hex_float (-0.0) print_hex_float (sqrt (-0.0))

+0.0 prints as "0" and -0.0 prints as "8000000000000000".

I confirm correct results under Linux, in 64-bit and 32-bit modes, including with "ocamlopt -ffast-math" in 32-bit mode. Moreover, the x87 "fsqrt" instruction and the SSE2 "sqrtsd" instructions produce the correct result also.

Concerning:

printf "%.16e " (-2.**(-1022.));;

2.2250738585072014e-308 - : unit = ()

Unary minus binds more tightly than "**", so you're raising -2.0 to the -1022.0 power, which is a positive float indeed.

Comment author: dpariente

(Ooops sorry for misused precedence!)

print_hex_float procedure yields:

$ ocamlc s.ml -o s.exe ; ./s.exe
sqrt(8000000000000000) = 8000000000000000
sqrt(0) = 0

$ ocamlopt s.ml -o s.exe ; ./s.exe
sqrt(8000000000000000) = fff8000000000000
sqrt(0) = 0

Trying to compile a C code equivalent embedding a call to C standard function sqrt(-0.), it returns:

$ i686-w64-mingw32-gcc -o s.exe s.c -O -mms-bitfields -Wall -Wno-unused ; ./s.exe
sqrt(0.000000e+000) = 0.000000e+000 = 0x8000000000000000000

i686-w64-mingw32-gcc is the native C compiler called by ocamlopt, given by ocamlopt -config. This compiler seems harmless ...

Comment author: @xavierleroy

i686-w64-mingw32-gcc is the native C compiler called by ocamlopt, given by ocamlopt -config. This compiler seems harmless ...

Maybe I wasn't clear in earlier messages, but when GCC sees a call to sqrt(), it generates inline code that first invokes the processor's sqrt instruction, then tests for errors, and only then calls the sqrt() function from the standard library to set errno correctly. Because of this, you're not observing the sqrt() function from Microsoft's CRT library, which is the real source of your problem. Try compiling your C test with a Microsoft compiler and see.

I move to "resolve/no change required" this PR, as:

• For x86-64, the issue goes away in the upcoming 4.01 release because of [patch] Generate sqrtsd opcode instead of external call to sqrt on amd64 #5795.
• For x86-32, -ffast-math also solves the problem
• OCaml isn't responsible for lack of standard-conformance in Microsoft's libraries. You use their stuff at your own risks.

Comment author: @yakobowski

Unfortunately, I'm not convinced that we can use -ffast-math within Frama-C. Am I correct in thinking ocamlopt may then use 80 bits float? If this is the case, I fear that some transfer functions for floating-point operations will have double rounding issues.

Comment author: @xavierleroy

@boris: on x86-32 bits, ocamlopt generates 80-bit x87 instructions for all float arithmetic, + - * / included. So, the double rounding problem is already there in a big way, and you should talk with Pascal to see how Frama-C value analysis is resilient against double rounding.

An example of non-IEEE754-compliant behavior with -ffast-math is that fsin and fcos produce a meaningless result if |argument| >= 2^64.