Skip to content
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

Closed
vicuna opened this issue May 24, 2013 · 10 comments

Comments

Projects
None yet
1 participant
@vicuna
Copy link

commented May 24, 2013

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

@vicuna

This comment has been minimized.

Copy link
Author

commented May 24, 2013

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.

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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?

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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!

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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?

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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.

@vicuna

This comment has been minimized.

Copy link
Author

commented May 25, 2013

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 ...

@vicuna

This comment has been minimized.

Copy link
Author

commented Jun 9, 2013

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 #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.
@vicuna

This comment has been minimized.

Copy link
Author

commented Jun 9, 2013

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.

@vicuna

This comment has been minimized.

Copy link
Author

commented Jun 10, 2013

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.