Bug 9937 – CTFE floats don't overflow correctly

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()); ---

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.

> 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.

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?

(In reply to comment #3) > 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!

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.

(In reply to comment #5) > 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.

I'd like to see an example of code that produces less accurate results if it had more precision.

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

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/

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.

(In reply to comment #10) > 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.

(In reply to comment #11) > 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.

(In reply to comment #12) > (In reply to comment #11) > > > 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.

(In reply to comment #13) > 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.

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!

(In reply to comment #8) > 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

(In reply to comment #15) > 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!

(In reply to comment #15) > 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....?

(In reply to comment #18) > (In reply to comment #15) > > 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.

The problem in issue 13474 is a case where extra precision causes tests to fail.

(In reply to Walter Bright from comment #1) > 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. So D is not system language? How this http://cage.ugent.be/~klein/papers/floating-point.pdf and msum from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090 can be implemented in D for float and double? For my implementation see PR https://github.com/D-Programming-Language/phobos/pull/2513 This PR works well with GDC, LDC2 and DMD64. But not with DMD32.

Just for reference, the exact cause of the trouble for the specific case in issue 13474 is that floating point expressions can be evaluated with different precision within the same code. eg this assertion could fail: double a; double b; double c; assert(a + b - c == a + b - c); Because it could be evaluated as (a + b - c) == cast(double)(cast(real)a + cast(real)b - cast(real)c) This is a little disturbing.

The real problem is that dmd (as only compiler) ignores `cast(float)myReal` when the float is converted to a higher precision after that. This optimization is semantically incorrect b/c rounding down can have side effects (setting exception flags). It's also unexpected b/c C (and gdc/ldc) correctly perform the rounding. Also see https://github.com/dlang/druntime/pull/1621#issuecomment-250426515

Also see the discussion in issue 13474 and "The semantics of float, double, and real" at http://dlang.org/d-floating-point.html, which gives some arguments why dmd uses -ffast-math behavior by default behavior. I think the spec argumentation is a bit weak (also see https://issues.dlang.org/show_bug.cgi?id=13474#c21). If someone asks to round down (cast(float) or assignment to float) it seems sensible to do it. If you're writing performance sensitive numeric code you will avoid to mix precisions. DMD is the only C/D/C++ compiler skipping such those roundings by default (w/o a -ffast-math switch). In any case we need to align the default behavior of dmd/gdc/ldc.

toPrec intrinsic implemented in compiler, and functions added to druntime. https://github.com/dlang/dmd/pull/10654 https://github.com/dlang/druntime/pull/2864

@ibuclaw created dlang/dmd pull request #10686 "Fix Issue 9937: CTFE floats don't overflow correctly" fixing this issue: - Fix Issue 9937: CTFE floats don't overflow correctly Adjusts original test to include new toPrec intrinsic https://github.com/dlang/dmd/pull/10686

dlang/dmd pull request #10686 "Fix Issue 9937: CTFE floats don't overflow correctly" was merged into master: - ea47332fb00829eb5d234daa9f2c618d544a596b by Iain Buclaw: Fix Issue 9937: CTFE floats don't overflow correctly Adjusts original test to include new toPrec intrinsic https://github.com/dlang/dmd/pull/10686

dlang/dmd pull request #11055 "[dmd-cxx] Fix Issue 9937: CTFE floats don't overflow correctly" was merged into dmd-cxx: - 843a7fe13a18fae66bcf2af45b0a9ef7b241da36 by Iain Buclaw: [dmd-cxx] Fix Issue 9937: CTFE floats don't overflow correctly Adjusts original test to include new toPrec intrinsic https://github.com/dlang/dmd/pull/11055