www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.bugs - [Issue 9937] New: CTFE floats don't overflow correctly

reply d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937

           Summary: CTFE floats don't overflow correctly
           Product: D
           Version: D1 & D2
          Platform: All
        OS/Version: All
            Status: NEW
          Severity: normal
          Priority: P2
         Component: DMD
        AssignedTo: nobody puremagic.com
        ReportedBy: clugdbug yahoo.com.au



This should pass. The multiply by 2 must cause an infinity, at least when x is
stored.
---
int blah()
{
    float x = float.max;
    x *= 2;
    x /= 2;
    assert(x == float.infinity);
    return 1;
}

static assert(blah());
---

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 16 2013
next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937


Walter Bright <bugzilla digitalmars.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
                 CC|                            |bugzilla digitalmars.com
         Resolution|                            |INVALID



11:58:22 PDT ---
I'm going to disagree.

In D, any floating point algorithm that relies on a maximum precision is
broken. The compiler and runtime is allowed to do all such calculations at as
high a precision as they want to - the types only specify a minimum precision,
not a maximum.

This is unlike C, which requires a maximum precision upon assignment.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 16 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937


Don <clugdbug yahoo.com.au> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|RESOLVED                    |REOPENED
         Resolution|INVALID                     |




 In D, any floating point algorithm that relies on a maximum precision is


broken. The compiler and runtime is allowed to do all such calculations at as
high a precision as they want to - the types only specify a minimum precision,
not a maximum.

That's not exactly true. It's true that the intermediate results may be at
extended precision, for example, it would not be true that
    float x = float.max;
    assert( (x * 2) / 2 == float.infinity); // may fail

But float and double must be IEEE types.
Without this, it's impossible to write correct floating-point code.
The big problem here is that the CTFE behaviour is different to the runtime
behaviour.

I hit this with HalfFloat, I can't implement it correctly because of this
issue.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 18 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




13:42:16 PDT ---
You're correct for C, as those are the C rules. But this is not true for D.

Can you be more specific about what you're doing in CTFE that does not work?

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 18 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937






 You're correct for C, as those are the C rules. But this is not true for D.
 
 Can you be more specific about what you're doing in CTFE that does not work?


Define an IEEE single-precision floating-point constant.
Perhaps this is a better example:

void main()
{
         float f1 = float.max * 2;
   const float f2 = float.max * 2;
   assert(f1 == float.infinity);          // ok
   assert(f1 == f2);                      // ok
   assert(*&f2 == float.infinity);        // even this is OK
   assert(f2 == float.infinity);          // fails!!!!!!!!!!!!!!!!!
}

f2 apparently has two different values at the same time!

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 19 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




13:26:17 PDT ---
Yes, having a runtime computation be different from a compile time computation
is a potential result of doing floating point math at higher precision at
compile time.

However, this is as designed and is expected.

It is a programming error in D to rely on limited precision. The only time I've
seen legitimate code that relied on limited precision is as part of a test
suite - i.e. not a real application.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 19 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937






 Yes, having a runtime computation be different from a compile time computation
 is a potential result of doing floating point math at higher precision at
 compile time.
 
 However, this is as designed and is expected.
 
 It is a programming error in D to rely on limited precision. The only time I've
 seen legitimate code that relied on limited precision is as part of a test
 suite - i.e. not a real application.


Are you sure your recognized it? I've written very little floating point code
that didn't! Cephes relies on it, for example.

It's not difficult to cope with x87 using C's semantics, but this is not at all
the same as saying that you can add extra bits of precision *anywhere* without
impact.

It is impossible to generate correct results in the lowest bit, without knowing
the target precision. Otherwise you get double rounding.
I think this can be proved mathematically.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 22 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




13:00:12 PDT ---
I'd like to see an example of code that produces less accurate results if it
had more precision.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 22 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




There's a nice discussion here:

http://www.exploringbinary.com/double-rounding-errors-in-floating-point-conversions/

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 22 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




And a real-world example, of where it caused a vulnerability in both Java and
PHP.

http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 23 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




18:23:56 PDT ---
To quote from
http://docs.oracle.com/cd/E19060-01/stud8.compiler/817-0932/ncg_goldberg.html:

“Of course, this form of double-rounding is highly unlikely to affect any
practical program adversely.”

Therefore, rather than retard the precision of all programs, I suggest creating
two standard functions,

   double roundToDouble(real x);
   float roundToFloat(real x);

which will guarantee a result rounded to the return type precision.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 23 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937






 To quote from
 http://docs.oracle.com/cd/E19060-01/stud8.compiler/817-0932/ncg_goldberg.html:
 
 “Of course, this form of double-rounding is highly unlikely to affect any
 practical program adversely.”
 
 Therefore, rather than retard the precision of all programs,


This is not about retarding precision of all programs (the way Java does). This
is about making compile-time behave the same as runtime.

Specifically:
1. A manifest constant should not be carrying extra precision. That's simply
incorrect.
2. JIT compilation onto the target machine should be valid for CTFE.
3. The results of CTFE should be as reproducible as runtime results are.

Right now, there is lots of stuff which appears to work, but relies on the
current DMD implementation.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 24 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




14:24:50 PDT ---



 This is about making compile-time behave the same as runtime.


Floating point answers can (and do) vary depending on optimization settings.
This is true for C and C++ compilers, too, even though the arithmetic is
Standard compliant.


 Specifically:
 1. A manifest constant should not be carrying extra precision. That's simply
 incorrect.


This is an assertion that earlier I would have accepted without question. I
question it now.


 2. JIT compilation onto the target machine should be valid for CTFE.


I agree.


 3. The results of CTFE should be as reproducible as runtime results are.


The definition should provide a minimum precision, not a maximum. As long as
the runtime results comply with a minimum precision, they are D standard
compliant.

I feel kinda strongly about this.

(Note that Go uses infinite precision for all integer constant folding. I think
that is superior to what D does, which limits the precision to the target
integer size.)

My earlier proposal to provide functions roundToFloat and roundToDouble at
least provide solid documentation at points in an algorithm where precision
ceilings are required.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 24 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937







 

 This is about making compile-time behave the same as runtime.


 
 Floating point answers can (and do) vary depending on optimization settings.
 This is true for C and C++ compilers, too, even though the arithmetic is
 Standard compliant.


Sure, but in practice the semantics are predictable. It's not "all bets are
off". The reality is, x87, 68K and Itanium may do intermediate operations at 80
bit precision, some embedded systems use alias float = double
and all other systems use the precision of the largest operands.

There are only these three possibilities to worry about. And there will never
be any more. All new systems will use the precision of the largest operands.

That's completely different from "the semantics are not even consistent during
compilation of a single module".


 Specifically:
 1. A manifest constant should not be carrying extra precision. That's simply
 incorrect.


 
 This is an assertion that earlier I would have accepted without question. I
 question it now.


If you can have a floating point constant that's between float.max and
float.infinity, it's not an IEEE value. That sucks.


 2. JIT compilation onto the target machine should be valid for CTFE.


 
 I agree.
 

 3. The results of CTFE should be as reproducible as runtime results are.


 
 The definition should provide a minimum precision, not a maximum. As long as
 the runtime results comply with a minimum precision, they are D standard
 compliant.
 
 I feel kinda strongly about this.


So do I. I think it is completely wrong.
The D spec also says that you can assume IEEE behaviour. If extra arbitrary
precision is possible, it is not IEEE. It's back to the bad old days.
Basically it makes ALL floating point operations into implementation-defined
behaviour!


 (Note that Go uses infinite precision for all integer constant folding. I think
 that is superior to what D does, which limits the precision to the target
 integer size.)


I think infinite precision is a big mistake. When CTFE becomes capable of
running BigInt, it can be used in those rare cases when it's actually required.


 My earlier proposal to provide functions roundToFloat and roundToDouble at
 least provide solid documentation at points in an algorithm where precision
 ceilings are required.


That would be a useful guarantee. That would be enough to deal with my original
test case.
It does not deal with the second one.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 25 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937


bearophile_hugs eml.cc changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bearophile_hugs eml.cc






 I think infinite precision is a big mistake. When CTFE becomes capable of
 running BigInt, it can be used in those rare cases when it's actually required.


I think multi-precision for integral constant folding is a different discussion
topic. I think it's better to keep this discussion focused on floating point
values only.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 25 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937




12:11:21 PDT ---
Don, I think we've reached an impasse here. How about we duke it out at Dcon? I
think some beers will get us to a solution!

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 25 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937






 There's a nice discussion here:
 
 http://www.exploringbinary.com/double-rounding-errors-in-floating-point-conversions/


That's a nice article. I have converted his two examples to a single D program:


import std.stdio;

void main() {
    real   r = 0.5000000894069671353303618843710864894092082977294921875L;
    double d = 0.5000000894069671353303618843710864894092082977294921875;
    float  f = 0.5000000894069671353303618843710864894092082977294921875f;
    float fd = 0.5000000894069671353303618843710864894092082977294921875;

    writefln("r =  %a", r);
    writefln("d =  %a", d);
    writefln("f =  %a", f);
    writefln("fd = %a", fd);

    double r2 = 0.5000000298023224154508881156289135105907917022705078125L;
    double d2 = 0.5000000298023224154508881156289135105907917022705078125;
    float  f2 = 0.5000000298023224154508881156289135105907917022705078125f;
    float fd2 = 0.5000000298023224154508881156289135105907917022705078125;

    writefln("r2  = %a", r2);
    writefln("d2  = %a", d2);
    writefln("f2  = %a", f2);
    writefln("fd2 = %a", fd2);
}


It seems dmd on Windows lacks both rounding problems discussed there. Output
(the gdc and ldc tests are done on dpaste):

dmd 32 bit Windows:
r =  0x1.000002fffffffbfep-1
d =  0x1.000003p-1
f =  0x1.000002p-1
fd = 0x1.000002p-1
r2  = 0x1.000001p-1
d2  = 0x1.000001p-1
f2  = 0x1.000002p-1
fd2 = 0x1.000002p-1

GDC 2.060 64 bit Linux:
r = 0x8.000017ffffffep-4
d = 0x1.000003p-1
f = 0x1.000002p-1
fd = 0x1.000002p-1
r2 = 0x1.000001p-1
d2 = 0x1.000001p-1
f2 = 0x1.000002p-1
fd2 = 0x1.000002p-1

LDC 2.060 64 bit Linux:
r = 0x8.000017ffffffep-4
d = 0x1.000003p-1
f = 0x1.000004p-1
fd = 0x1.000004p-1
r2 = 0x1.000001p-1
d2 = 0x1.000001p-1
f2 = 0x1p-1
fd2 = 0x1p-1

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 25 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937






 Don, I think we've reached an impasse here. How about we duke it out at Dcon? I
 think some beers will get us to a solution!


I was going to suggest the same thing!

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Apr 26 2013
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937


yebblies <yebblies gmail.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |yebblies gmail.com





 Don, I think we've reached an impasse here. How about we duke it out at Dcon? I
 think some beers will get us to a solution!


And the result is....?

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Aug 01 2013
prev sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=9937








 Don, I think we've reached an impasse here. How about we duke it out at Dcon? I
 think some beers will get us to a solution!


 
 And the result is....?


Well, we didn't spend much time on it. I've realized a few things since then:

(1) My statement about x87 being the only case where you can get extra
precision, is wrong. Machines with FMA can effectively have twice double
precision on some operations. The 80-precision that DMD currently uses, is
actually not enough, even for x86-64.

(2) We have a few things conflated here.
The most severe bug in the current situation is the example in comment 4: the
compiler acts as though extra bits of precision can be transferred from compile
time to run time. This has a simple solution, see below.

(3) The solution to my initial test case is less clear. There needs to be a way
to reliably dump extra precision from a calculation. Walter's example in
comment 10 is one possibility.

----
Now, to fix comment 4:
All that needs to happen to fix this, is that floating point values should be
rounded to fit the target type *at the end of CTFE*, as part of the existing
scrubbing process.
They should not be carrying extra precision that could not possibly exist at
runtime. This rounding does not apply to const-folding (where things such as
inlining could happen which might preserve more precision).

The philosophy then would be, const-folding never reduces precision below what
the machine can support. But in cases where a CTFE value must be created, such
as in a template value parameter, the result has a single value throughout the
program. Any extra precision it used in generating that value is entirely
CTFE's business and should not leak into the rest of the compiler.

Now that I've split CTFE off completely from const-folding, it's possible to
make this distinction. This would fix the issue in comment 4, but doesn't
affect anything else.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Aug 22 2013