Mantis Bug Tracker

View Issue Details Jump to Notes ] Issue History ] Print ]
IDProjectCategoryView StatusDate SubmittedLast Update
0006020OCamlOCaml generalpublic2013-05-24 17:382013-06-10 13:39
Reporterdpariente 
Assigned To 
PrioritynormalSeverityminorReproducibilityalways
StatusresolvedResolutionno change required 
PlatformPort MinGWOSWin7-64bits and XP-32bitsOS Version
Product Version4.00.1 
Target VersionFixed in Version 
Summary0006020: sqrt(-0.) does not yield -0.
DescriptionThe 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
TagsNo tags attached.
Attached Files

- Relationships
related to 0005795resolvedlefessan [patch] Generate sqrtsd opcode instead of external call to sqrt on amd64 

-  Notes
(0009333)
xleroy (administrator)
2013-05-24 19:30

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.
(0009335)
dpariente (reporter)
2013-05-25 09:24

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
(0009336)
Boris Yakobowski (reporter)
2013-05-25 11:29

@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?
(0009338)
dpariente (reporter)
2013-05-25 12:22

@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!
(0009340)
dpariente (reporter)
2013-05-25 12:59

@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?
(0009342)
xleroy (administrator)
2013-05-25 14:15

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.
(0009343)
dpariente (reporter)
2013-05-25 18:52

(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 ...
(0009445)
xleroy (administrator)
2013-06-09 11:22

> 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 PR#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.
(0009446)
Boris Yakobowski (reporter)
2013-06-09 12:23

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.
(0009451)
xleroy (administrator)
2013-06-10 13:39

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

- Issue History
Date Modified Username Field Change
2013-05-24 17:38 dpariente New Issue
2013-05-24 19:30 xleroy Note Added: 0009333
2013-05-24 19:30 xleroy Description Updated View Revisions
2013-05-24 19:31 xleroy Status new => acknowledged
2013-05-25 09:24 dpariente Note Added: 0009335
2013-05-25 11:29 Boris Yakobowski Note Added: 0009336
2013-05-25 12:22 dpariente Note Added: 0009338
2013-05-25 12:59 dpariente Note Added: 0009340
2013-05-25 14:15 xleroy Note Added: 0009342
2013-05-25 18:52 dpariente Note Added: 0009343
2013-06-09 11:14 xleroy Relationship added related to 0005975
2013-06-09 11:14 xleroy Relationship deleted related to 0005975
2013-06-09 11:14 xleroy Relationship added related to 0005795
2013-06-09 11:22 xleroy Note Added: 0009445
2013-06-09 11:23 xleroy Status acknowledged => resolved
2013-06-09 11:23 xleroy Resolution open => no change required
2013-06-09 12:23 Boris Yakobowski Note Added: 0009446
2013-06-10 13:39 xleroy Note Added: 0009451


Copyright © 2000 - 2011 MantisBT Group
Powered by Mantis Bugtracker