Anonymous  Login  Signup for a new account  20160731 07:35 CEST 
Main  My View  View Issues  Change Log  Roadmap 
View Issue Details [ Jump to Notes ]  [ Issue History ] [ Print ]  
ID  Project  Category  View Status  Date Submitted  Last Update  
0006020  OCaml  OCaml general  public  20130524 17:38  20151211 19:19  
Reporter  dpariente  
Assigned To  
Priority  normal  Severity  minor  Reproducibility  always  
Status  closed  Resolution  no change required  
Platform  Port MinGW  OS  Win764bits and XP32bits  OS Version  
Product Version  4.00.1  
Target Version  Fixed in Version  
Summary  0006020: sqrt(0.) does not yield 0.  
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  
Tags  No tags attached.  
Attached Files  
Relationships  

Notes  
(0009333) xleroy (administrator) 20130524 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 runtime 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 ffastmath", 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) 20130525 09:24 
Thank you for these explanations. However, for information, in 32 bits, with ffastmath the result is still not the correct one: $ ocamlopt ffastmath s.ml o s.exe ; ./s.exe sqrt(0.0000000000000000e+000) = 0.0000000000000000e+000 
(0009336) Boris Yakobowski (reporter) 20130525 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) 20130525 12:22 
@Boris: Right, the minus sign is systematically missing under XP32b. I don't know why, and suspected a trouble in the prettyprinting, but maybe it's something else! 
(0009340) dpariente (reporter) 20130525 12:59 
@note9338 Under toplevel: # printf "%.16e " (2.**(1022.));; 2.2250738585072014e308  : unit = () # printf "%.16e " (0. . 2.**(1022.));; 2.2250738585072014e308  : unit = () An issue when parsing the minus? 
(0009342) xleroy (administrator) 20130525 14:15 
Boris is right that you should not trust the printing of floatingpoint numbers either, since it goes through the C standard library as well. Better print the IEEE754 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 64bit and 32bit modes, including with "ocamlopt ffastmath" in 32bit mode. Moreover, the x87 "fsqrt" instruction and the SSE2 "sqrtsd" instructions produce the correct result also. Concerning: > # printf "%.16e " (2.**(1022.));; > 2.2250738585072014e308  : 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) 20130525 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: $ i686w64mingw32gcc o s.exe s.c O mmsbitfields Wall Wnounused ; ./s.exe sqrt(0.000000e+000) = 0.000000e+000 = 0x8000000000000000000 i686w64mingw32gcc is the native C compiler called by ocamlopt, given by ocamlopt config. This compiler seems harmless ... 
(0009445) xleroy (administrator) 20130609 11:22 
> i686w64mingw32gcc 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 x8664, the issue goes away in the upcoming 4.01 release because of PR#5795.  For x8632, ffastmath also solves the problem  OCaml isn't responsible for lack of standardconformance in Microsoft's libraries. You use their stuff at your own risks. 
(0009446) Boris Yakobowski (reporter) 20130609 12:23 
Unfortunately, I'm not convinced that we can use ffastmath within FramaC. Am I correct in thinking ocamlopt may then use 80 bits float? If this is the case, I fear that some transfer functions for floatingpoint operations will have double rounding issues. 
(0009451) xleroy (administrator) 20130610 13:39 
@Boris: on x8632 bits, ocamlopt generates 80bit 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 FramaC value analysis is resilient against double rounding. An example of nonIEEE754compliant behavior with ffastmath is that fsin and fcos produce a meaningless result if argument >= 2^64. 
Issue History  
Date Modified  Username  Field  Change 
20130524 17:38  dpariente  New Issue  
20130524 19:30  xleroy  Note Added: 0009333  
20130524 19:30  xleroy  Description Updated  View Revisions 
20130524 19:31  xleroy  Status  new => acknowledged 
20130525 09:24  dpariente  Note Added: 0009335  
20130525 11:29  Boris Yakobowski  Note Added: 0009336  
20130525 12:22  dpariente  Note Added: 0009338  
20130525 12:59  dpariente  Note Added: 0009340  
20130525 14:15  xleroy  Note Added: 0009342  
20130525 18:52  dpariente  Note Added: 0009343  
20130609 11:14  xleroy  Relationship added  related to 0005975 
20130609 11:14  xleroy  Relationship deleted  related to 0005975 
20130609 11:14  xleroy  Relationship added  related to 0005795 
20130609 11:22  xleroy  Note Added: 0009445  
20130609 11:23  xleroy  Status  acknowledged => resolved 
20130609 11:23  xleroy  Resolution  open => no change required 
20130609 12:23  Boris Yakobowski  Note Added: 0009446  
20130610 13:39  xleroy  Note Added: 0009451  
20151211 19:19  xleroy  Status  resolved => closed 
Copyright © 2000  2011 MantisBT Group 