www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Floating point rounding

reply Guillaume Chatelet <chatelet.guillaume gmail.com> writes:
I would expect that (1.0f + smallest float subnormal) > 1.0f when 
the Floating Point unit is set to Round Up.

Here is some C++ code:
#include <limits>
#include <cstdio>
#include <cfenv>

int main(int, char**) {
     std::fesetround(FE_UPWARD);
     printf("%.32g\n", std::numeric_limits<float>::denorm_min() + 
1.0f);
     return 0;
}

Execution on my machine yields:
clang++ --std=c++11 test_denormal.cc && ./a.out
1.00000011920928955078125

Here is the same code in D:
void main(string[] args)
{
     import std.math;
     FloatingPointControl fpctrl;
     fpctrl.rounding = FloatingPointControl.roundUp;
     writefln("%.32g", float.min_normal + 1.0f);
}

Execution on my machine yields:
dmd -run test_denormal.d
1

Did I miss something?
Mar 02 2017
parent reply Guillaume Chatelet <chatelet.guillaume gmail.com> writes:
On Thursday, 2 March 2017 at 20:30:47 UTC, Guillaume Chatelet 
wrote:

 Here is the same code in D:
 void main(string[] args)
 {
     import std.math;
     FloatingPointControl fpctrl;
     fpctrl.rounding = FloatingPointControl.roundUp;
     writefln("%.32g", float.min_normal + 1.0f);
 }

 Execution on my machine yields:
 dmd -run test_denormal.d
 1

 Did I miss something?


This example is closer to the C++ one:

void main(string[] args)
{
     import core.stdc.fenv;
     fesetround(FE_UPWARD);
     writefln("%.32g", float.min_normal + 1.0f);
}

It still yields "1"
Mar 02 2017
parent reply ag0aep6g <anonymous example.com> writes:
On 03/02/2017 10:10 PM, Guillaume Chatelet wrote:

 On Thursday, 2 March 2017 at 20:30:47 UTC, Guillaume Chatelet wrote:

 Here is the same code in D:
 void main(string[] args)
 {
     import std.math;
     FloatingPointControl fpctrl;
     fpctrl.rounding = FloatingPointControl.roundUp;
     writefln("%.32g", float.min_normal + 1.0f);
 }

 Execution on my machine yields:
 dmd -run test_denormal.d
 1

 Did I miss something?


 This example is closer to the C++ one:

 void main(string[] args)
 {
     import core.stdc.fenv;
     fesetround(FE_UPWARD);
     writefln("%.32g", float.min_normal + 1.0f);
 }

 It still yields "1"


This prints the same as the C++ version:

----
void main(string[] args)
{
     import std.stdio;
     import core.stdc.fenv;
     fesetround(FE_UPWARD);
     float x = 1.0f;
     x += float.min_normal;
     writefln("%.32g", x);
}
----

Soo, a bug/limitation of constant folding?

With FloatingPointControl it still prints "1". Does 
FloatingPointControl.rounding do something different than fesetround? 
The example in the docs [1] only shows how it changes rint's behavior.
Mar 02 2017
parent reply Guillaume Chatelet <chatelet.guillaume gmail.com> writes:
On Thursday, 2 March 2017 at 21:34:56 UTC, ag0aep6g wrote:

 On 03/02/2017 10:10 PM, Guillaume Chatelet wrote:

 On Thursday, 2 March 2017 at 20:30:47 UTC, Guillaume Chatelet 
 wrote:

 Here is the same code in D:
 void main(string[] args)
 {
     import std.math;
     FloatingPointControl fpctrl;
     fpctrl.rounding = FloatingPointControl.roundUp;
     writefln("%.32g", float.min_normal + 1.0f);
 }

 Execution on my machine yields:
 dmd -run test_denormal.d
 1

 Did I miss something?


 This example is closer to the C++ one:

 void main(string[] args)
 {
     import core.stdc.fenv;
     fesetround(FE_UPWARD);
     writefln("%.32g", float.min_normal + 1.0f);
 }

 It still yields "1"


 This prints the same as the C++ version:

 ----
 void main(string[] args)
 {
     import std.stdio;
     import core.stdc.fenv;
     fesetround(FE_UPWARD);
     float x = 1.0f;
     x += float.min_normal;
     writefln("%.32g", x);
 }
 ----

 Soo, a bug/limitation of constant folding?

 With FloatingPointControl it still prints "1". Does 
 FloatingPointControl.rounding do something different than 
 fesetround? The example in the docs [1] only shows how it 
 changes rint's behavior.





Thx for the investigation!
Here is the code for FloatingPointControl
https://github.com/dlang/phobos/blob/master/std/math.d#L4809

Other code (enableExceptions / disableExceptions) seems to have 
two code path depending on "version(X86_Any)", rounding doesn't.

Maybe that's the bug?
Mar 02 2017
parent reply ag0aep6g <anonymous example.com> writes:
On 03/02/2017 10:49 PM, Guillaume Chatelet wrote:

 Thx for the investigation!
 Here is the code for FloatingPointControl
 https://github.com/dlang/phobos/blob/master/std/math.d#L4809

 Other code (enableExceptions / disableExceptions) seems to have two code
 path depending on "version(X86_Any)", rounding doesn't.

 Maybe that's the bug?


I don't think that's it. But I think I've figured out what's wrong:

dmd generates SSE instructions for floating point math. 
FloatingPointControl only minds the control register for the FPU. But 
SSE instructions are not affected by that. SSE has a separate control 
register: MXCSR.

Too see that dmd generates SSE instructions:

     echo 'void main() { float x; x += 0.1f; }' > test.d && \
     dmd -c test.d && \
     objdump -d -Mintel -j.text._Dmain test.o

Note the instructions movss and addss, and the xmm registers.

To affect those instructions, FloatinPointControl.setsetControlState [1] 
must use ldmxcsr in addition to or instead of fldcw.

This should be a pretty obvious bug. Rounding (and floating point 
exceptions?) must not be covered in the phobos tests at all.


[1] 
https://github.com/dlang/phobos/blob/750593738bcbe9577c3f5951b5b639392871f90e/std/math.d#L4906-L4933
Mar 03 2017
parent ag0aep6g <anonymous example.com> writes:
On 03/03/2017 05:39 PM, ag0aep6g wrote:

 dmd generates SSE instructions for floating point math.
 FloatingPointControl only minds the control register for the FPU. But
 SSE instructions are not affected by that. SSE has a separate control
 register: MXCSR.


I've filed an issue:
     https://issues.dlang.org/show_bug.cgi?id=17243
and am attempting to fix it:
     https://github.com/dlang/phobos/pull/5240
Mar 04 2017