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
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
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. [1] http://dlang.org/phobos/std_math.html#.FloatingPointControl
Mar 02
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. [1] http://dlang.org/phobos/std_math.html#.FloatingPointControl
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
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
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