www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Fun with floating point

reply "Kenny" <artemalive gmail.com> writes:
Hi, D community!

I have this program:

import std.stdio;
import std.conv;

int main(string[] argv)
{
       float eps = 1.0f;
       float f = 0.0f;
       while (f + eps != f)
           f += 1.0f;

       writeln("eps = " ~ to!string(eps) ~ ", max_f = " ~
to!string(f));
       return 0;
}

According to the languge specification what result would you
expect from its execution?

When running similar C++ program (VS 2013) the loop terminates
and I get t = 2^24 = 16777216.

Does D language specifies that loop will be terminated for this
program or  ?

I compiled with DMD and it hungs.
Details about assembly generated by DMD can be found here:
http://stackoverflow.com/questions/28380651/floating-point-maxing-out-loop-doesnt-terminate-in-d-works-in-c

Thanks.
Feb 07 2015
next sibling parent "Kenny" <artemalive gmail.com> writes:
And sory for the typos, cannot find edit functionality here..
Feb 07 2015
prev sibling next sibling parent reply "Kenny" <artemalive gmail.com> writes:
The above code snippet works correctly when I use LDC compiler 
(it finds expected 'f' value and prints it to console). I'm 
wondering is it a bug in DMD?

p.s. the final code used by both compilers:

import std.stdio;
import std.conv;

int main(string[] argv)
{
     const float eps = 1.0f;
     float f = 0.0f;
     while (f + eps != f)
         f += 1.0f;

     writeln("eps = ", eps, ", max_f = ", f);
     return 0;
}
Feb 07 2015
next sibling parent reply =?UTF-8?B?QWxpIMOHZWhyZWxp?= <acehreli yahoo.com> writes:
To answer your other question, there is no Edit because this is a 
newsgroup (see NNTP). The forum interface is supposed to be a 
convenience but it hides that fact.

On 02/07/2015 01:33 PM, Kenny wrote:


 The above code snippet works correctly


There is no right or wrong when you compare floating point values for 
equality (and inequality) unless those values can be represented exactly 
in the machine. 1.0 is famously not representable exactly. (It is 
similar to how 1/3 cannot be represented in the decimal system.)


 when I use LDC compiler (it finds
 expected 'f' value and prints it to console). I'm wondering is it a bug
 in DMD?


Not a bug.

Ali


 p.s. the final code used by both compilers:

 import std.stdio;
 import std.conv;

 int main(string[] argv)
 {
      const float eps = 1.0f;
      float f = 0.0f;
      while (f + eps != f)
          f += 1.0f;

      writeln("eps = ", eps, ", max_f = ", f);
      return 0;
 }
Feb 07 2015
next sibling parent reply "anonymous" <anonymous example.com> writes:
On Saturday, 7 February 2015 at 22:46:56 UTC, Ali Çehreli wrote:

 1.0 is famously not representable exactly.


1.0 is representable exactly, though.
Feb 07 2015
parent "Peter Alexander" <peter.alexander.au gmail.com> writes:
On Saturday, 7 February 2015 at 23:06:15 UTC, anonymous wrote:

 On Saturday, 7 February 2015 at 22:46:56 UTC, Ali Çehreli wrote:

 1.0 is famously not representable exactly.


 1.0 is representable exactly, though.


I think he meant 0.1 :-)
Feb 07 2015
prev sibling parent reply "Kenny" <artemalive gmail.com> writes:
 There is no right or wrong when you compare floating point 
 values for equality (and inequality) unless those values can be 
 represented exactly in the machine. 1.0 is famously not 
 representable exactly. (It is similar to how 1/3 cannot be 
 represented in the decimal system.)


Here I tested one aspect of IEEE-754 floating point behavior (not 
ordinary floating point calculations for daily job). All integers 
can be represented exactly until some maximum value.

I believe that first section from http://dlang.org/float.html 
explains why behavior of my program should not be considered a 
bug in the compiler. But does it mean that due to these rules the 
results of the computation are not according to IEEE-754 for some 
specified floating point format (32 bit single precision in my 
examples)?

For example, according to IEEE-754 specification if we work with 
32bit floating point numbers (floats in D) then the following 
pseudo code prints 'Works'.

F32 f = 16777216.0f;
F32 f2 = f + 0.1f;
if is_the_same_binary_presentation(f, f2)
   Print("Works");

As I understand D does not guarantee this. Please confirm if it's 
correct.
Feb 08 2015
parent reply "Kenny" <artemalive gmail.com> writes:
 For example, according to IEEE-754 specification if we work 
 with 32bit floating point numbers (floats in D) then the 
 following pseudo code prints 'Works'.

 F32 f = 16777216.0f;
 F32 f2 = f + 0.1f;
 if is_the_same_binary_presentation(f, f2)
   Print("Works");

 As I understand D does not guarantee this. Please confirm if 
 it's correct.


One clarification:

In D the following should always print 'Works'
float f = 16777216.0f;
float f2 = f + 1.0f;
if (f == f2)
  writeln("Works");

I asked more about this case:
float f = 16777216.0f;
if (f == f + 1.0f)
  writeln("Works");

Although all operands are 32bit FP the result is not guaranteed 
to be equal to the result for FP32 computations as specified by 
IEEE-754.
Feb 08 2015
next sibling parent reply ketmar <ketmar ketmar.no-ip.org> writes:
On Sun, 08 Feb 2015 09:19:06 +0000, Kenny wrote:


 I asked more about this case:
 float f =3D 16777216.0f;
 if (f =3D=3D f + 1.0f)
   writeln("Works");
=20
 Although all operands are 32bit FP the result is not guaranteed to be
 equal to the result for FP32 computations as specified by IEEE-754.


i think you are mixing two things here. IEEE doesn't specify which=20
internal representation compilers should use, it only specifies the=20
results for chosen representation. so if D specs states that `float`=20
calculations are always performing with `float` precision (and specs=20
aren't), your sample should work.

but the specs says that compiler is free to promote such expression to=20
any type it wants, it just should not loose precision. so the actual type=20
of `f` in `f =3D=3D f + 1.0f` can be freely promoted to `double`, `real` or=
=20
even to some GMP representation.=
Feb 08 2015
parent "Kenny" <artemalive gmail.com> writes:
 i think you are mixing two things here. IEEE doesn't specify 
 which
 internal representation compilers should use, it only specifies 
 the
 results for chosen representation. so if D specs states that 
 `float`
 calculations are always performing with `float` precision (and 
 specs
 aren't), your sample should work.

 but the specs says that compiler is free to promote such 
 expression to
 any type it wants, it just should not loose precision. so the 
 actual type
 of `f` in `f == f + 1.0f` can be freely promoted to `double`, 
 `real` or
 even to some GMP representation.


It's clear now, thanks!. Also thanks to everyone for the answers, 
it was very helpful.
Feb 08 2015
prev sibling parent reply "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Sunday, 8 February 2015 at 09:19:08 UTC, Kenny wrote:

 For example, according to IEEE-754 specification if we work 
 with 32bit floating point numbers (floats in D) then the 
 following pseudo code prints 'Works'.

 F32 f = 16777216.0f;
 F32 f2 = f + 0.1f;
 if is_the_same_binary_presentation(f, f2)
  Print("Works");

 As I understand D does not guarantee this. Please confirm if 
 it's correct.


 One clarification:

 In D the following should always print 'Works'
 float f = 16777216.0f;
 float f2 = f + 1.0f;
 if (f == f2)
  writeln("Works");

 I asked more about this case:
 float f = 16777216.0f;
 if (f == f + 1.0f)
  writeln("Works");

 Although all operands are 32bit FP the result is not guaranteed 
 to be equal to the result for FP32 computations as specified by 
 IEEE-754.


«Algorithms should be written to work based on the minimum 
precision of the calculation. They should not degrade or fail if 
the actual precision is greater.»

http://dlang.org/float.html

:-/
Feb 08 2015
parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
Also note that denormal numbers is an issue, e.g. D assumes that 
they are always supported, I think. Which frequently is not the 
case...
Feb 08 2015
prev sibling next sibling parent "Peter Alexander" <peter.alexander.au gmail.com> writes:
On Saturday, 7 February 2015 at 21:33:51 UTC, Kenny wrote:

 The above code snippet works correctly when I use LDC compiler 
 (it finds expected 'f' value and prints it to console). I'm 
 wondering is it a bug in DMD?

 p.s. the final code used by both compilers:

 import std.stdio;
 import std.conv;

 int main(string[] argv)
 {
     const float eps = 1.0f;
     float f = 0.0f;
     while (f + eps != f)
         f += 1.0f;

     writeln("eps = ", eps, ", max_f = ", f);
     return 0;
 }


Intermediate calculations may be performed at higher precision 
than the precision of the values themselves. In particular, the f 
+ eps may be performed with 80 bits of precision, even though 
both values are 32-bit. The comparison will then fail.

The reason for the difference between DMD and LDC is that DMD 
tends to use the FPU more with 80 bits of precision, whereas LDC 
and GDC will use the SSE2 instructions, which only support 32-bit 
and 64-bit precision.
Feb 07 2015
prev sibling next sibling parent =?UTF-8?B?QWxpIMOHZWhyZWxp?= <acehreli yahoo.com> writes:
On 02/07/2015 01:33 PM, Kenny wrote:

 The above code snippet works correctly when I use LDC compiler (it finds
 expected 'f' value and prints it to console). I'm wondering is it a bug
 in DMD?

 p.s. the final code used by both compilers:

 import std.stdio;
 import std.conv;

 int main(string[] argv)
 {
      const float eps = 1.0f;
      float f = 0.0f;
      while (f + eps != f)
          f += 1.0f;

      writeln("eps = ", eps, ", max_f = ", f);
      return 0;
 }


OK, ignore some of my earlier response. :)

The code above works with dmd git head 64-bit compilation and prints the 
following:

eps = 1, max_f = 1.67772e+07

You can use the %a format specifier when debugging this issue. It allows 
you see the bits of the floating point value:

       writefln("eps: %a", 0.1);

double 0.1 on my system:

eps: 0x1.999999999999ap-4

Ali
Feb 07 2015
prev sibling parent reply ketmar <ketmar ketmar.no-ip.org> writes:
On Sat, 07 Feb 2015 21:33:46 +0000, Kenny wrote:


 The above code snippet works correctly when I use LDC compiler (it finds
 expected 'f' value and prints it to console). I'm wondering is it a bug
 in DMD?


nope, this is a bug in your code. compiler (by the specs) is free to=20
perform intermediate calculations with any precision that is not lower=20
than a highest used type (i.e. not lower that `float`'s one for `while`=20
condition (`f + eps !=3D f`). it may be even infinite precision, so your=20
code may not exit the loop at all.=
Feb 08 2015
parent reply "Kenny" <artemalive gmail.com> writes:
 nope, this is a bug in your code. compiler (by the specs) is 
 free to
 perform intermediate calculations with any precision that is 
 not lower
 than a highest used type (i.e. not lower that `float`'s one for 
 `while`
 condition (`f + eps != f`). it may be even infinite precision, 
 so your
 code may not exit the loop at all.


Thanks, it's clear now. I still have one question in the above 
post, I would appreciate if you check it too.
Feb 08 2015
parent ketmar <ketmar ketmar.no-ip.org> writes:
On Sun, 08 Feb 2015 09:05:30 +0000, Kenny wrote:


 Thanks, it's clear now. I still have one question in the above post, I
 would appreciate if you check it too.


i've seen that, but i don't know the answer, sorry. here we have to=20
summon Walter to explain what his intentions was, how it should work and=20
why.=
Feb 08 2015
prev sibling parent "Meta" <jared771 gmail.com> writes:
On Saturday, 7 February 2015 at 16:06:14 UTC, Kenny wrote:

 Hi, D community!

 I have this program:

 import std.stdio;
 import std.conv;

 int main(string[] argv)
 {
       float eps = 1.0f;
       float f = 0.0f;
       while (f + eps != f)
           f += 1.0f;

       writeln("eps = " ~ to!string(eps) ~ ", max_f = " ~
 to!string(f));
       return 0;
 }

 According to the languge specification what result would you
 expect from its execution?

 When running similar C++ program (VS 2013) the loop terminates
 and I get t = 2^24 = 16777216.

 Does D language specifies that loop will be terminated for this
 program or  ?

 I compiled with DMD and it hungs.
 Details about assembly generated by DMD can be found here:
 http://stackoverflow.com/questions/28380651/floating-point-maxing-out-loop-doesnt-terminate-in-d-works-in-c

 Thanks.


A point of advice that Walter gives for situations like these is 
to ensure that your algorithm has a minimum required precision to 
work correctly but not a maximum. As Peter Alexander mentioned, 
DMD performs intermediate calculations at higher precision than 
GDC/LDC and thus your algorithm breaks.
Feb 07 2015