New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
sqrt(-0.) does not yield -0. #6020
Comments
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: |
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. |
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 _ = +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:
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 $ ocamlopt s.ml -o s.exe ; ./s.exe 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 i686-w64-mingw32-gcc is the native C compiler called by ocamlopt, given by ocamlopt -config. This compiler seems harmless ... |
Comment author: @xavierleroy
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:
|
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. |
Original bug ID: 6020
Reporter: dpariente
Status: closed (set by @xavierleroy on 2015-12-11T18:19:34Z)
Resolution: not a bug
Priority: normal
Severity: minor
Platform: Port MinGW
OS: Win7-64bits and XP-32bits
Version: 4.00.1
Category: ~DO NOT USE (was: OCaml general)
Related to: #5795
Monitored by: monate @yakobowski
Bug description
The following code:
Format.printf "sqrt(%.16e) = %.16e\n" (-0.0) (sqrt (-0.0)) ;;
gives the results below which does not follow IEEE754 (sqrt(-0.)==-0.).
1°) Under Win7 64 bits (OCaml 4.00.1 port MinGW):
$ ocamlopt s.ml -o s.exe ; ./s.exe
sqrt(-0.0000000000000000e+000) = nan
(note: bytecode is correct and returns -0.)
2°) Under WinXP 32 bits (OCaml 4.00.1 port MinGW):
$ ocamlopt s.ml -o s.exe ; ./s.exe
sqrt(0.0000000000000000e+000) = nan
$ ocamlc s.ml -o s.exe ; ./s.exe
sqrt(0.0000000000000000e+000) = 0.0000000000000000e+000
The text was updated successfully, but these errors were encountered: