Always false float comparisons

[Timon Gehr via Digitalmars-d]

On 20.05.2016 11:14, Joakim wrote:
> On Thursday, 19 May 2016 at 18:22:48 UTC, Timon Gehr wrote:
>> On 19.05.2016 08:04, Joakim wrote:
>>> On Wednesday, 18 May 2016 at 17:10:25 UTC, Timon Gehr wrote:
>>>> It's not just slightly worse, it can cut the number of useful bits in
>>>> half or more! It is not unusual, I have actually run into those
>>>> problems in the past, and it can break an algorithm that is in Phobos
>>>> today!
>>>
>>> I wouldn't call that broken.  Looking at the hex output by replacing %f
>>> with %A in writefln, it appears the only differences in all those
>>> results is the last byte in the significand.
>>
>> Argh...
>>
>> // ...
>>
>> void main(){
>>     //double[] data=[1e16,1,-9e15];
>>     import std.range;
>>     double[] data=1e16~repeat(1.0,100000000).array~(-9e15);
>>     import std.stdio;
>>     writefln("%f",sum(data)); // baseline
>>     writefln("%f",kahan(data)); // kahan
>>     writefln("%f",kahanBroken(data)); // broken kahan
>> }
>>
>>
>> dmd -run kahanDemo.d
>> 1000000000000000.000000
>> 1000000100000000.000000
>> 1000000000000000.000000
>>
>> dmd -m32 -O -run kahanDemo.d
>> 1000000000000000.000000
>> 1000000000000000.000000
>> 1000000000000000.000000
>>
>>
>> Better?
>>
>> Obviously there is more structure in the data that I invent manually
>> than in a real test case where it would go wrong. The problems carry
>> over though.
>
> I looked over your code a bit.  If I define sum and c as reals in
> "kahanBroken" at runtime, this problem goes away.

Yes. That's absolutely obvious, and I have noted it before, but thanks. 
Maybe try to understand why this problem occurs in the first place.

> Since that's what the
> CTFE rule is actually doing, ie extending all floating-point to reals at
> compile-time, I don't see what you're complaining about.  Try it, run
> even your original naive summation algorithm through CTFE and it will
> produce the result you want:
>
> enum double[] ctData=[1e16,1,-9e15];
> enum ctSum = sum(ctData);
> writefln("%f", ctSum);
> ...

This example wasn't specifically about CTFE, but just imagine that only 
part of the computation is done at CTFE, all local variables are 
transferred to runtime and the computation is completed there.


>>> As Don's talk pointed out,
>>> all floating-point calculations will see loss of precision starting
>>> there.
>>> ...
>>
>>
>> This is implicitly assuming a development model where the programmer
>> first writes down the computation as it would be correct in the real
>> number system and then naively replaces every operation by the
>> rounding equivalent and hopes for the best.
>
> No, it is intrinsic to any floating-point calculation.
> ...

How do you even define accuracy if you don't specify an infinitely 
precise reference result?

>> It is a useful rule if that is what you're doing. One might be doing
>> something else. Consider the following paper for an example where the
>> last bit in the significant actually carries useful information for
>> many of the values used in the program.
>>
>> http://www.jaist.ac.jp/~s1410018/papers/qd.pdf
>
> Did you link to the wrong paper? ;)

No. That paper uses multiple doubles per approximated real value to 
implement arithmetic that is more precise than using just plain doubles. 
If any bit in the first double is off, this is no better than using a 
single double.

> I skimmed it and that paper
> explicitly talks about error bounds all over the place.

It is enough to read the abstract to figure out what the problem is. 
This demonstrates a non-contrived case where CTFE using enhanced 
precision throughout can break your program. Compute something as a 
double-double at compile-time, and when it is transferred to runtime you 
lose all the nice extra precision, because bits in the middle of the 
(conceptual) mantissa are lost.


> The only mention of "the last bit" is

This part is actually funny. Thanks for the laugh. :-)
I was going to say that your text search was too naive, but then I 
double-checked your claim and there are actually two mentions of "the 
last bit", and close by to the other mention, the paper says that "the 
first double a_0 is a double-precision approximation to the number a, 
accurate to almost half an ulp."

> when they say they calculated their
> constants in arbitrary precision before rounding them for runtime use,
> which is ironically similar to what Walter suggested doing for D's CTFE
> also.
> ...

Nothing "ironic" about that. It is sometimes a good idea and I can do 
this explicitly and make sure the rounding is done correctly, just like 
they did. Also, it is a lot more flexible if I can specify the exact way 
the computation is done and the result is rounded. 80 bits might not be 
enough anyway. There is no reason for the language to apply potentially 
incorrect "hand holding" here.

Again, please understand that my point is not that lower precision is 
better. My point is that doing the same thing in every context and 
allowing the programmer to specify what happens is better.

>>> In this case, not increasing precision gets the more accurate result,
>>> but other examples could be constructed that _heavily_ favor increasing
>>> precision.
>>
>> Sure. In such cases, you should use higher precision. What is the
>> problem? This is already supported (the compiler is not allowed to use
>> lower precision than requested).
>
> I'm not the one with the problem, you're the one complaining.
> ...

So you see no problem with my requested semantics for the built-in 
floating point types?

>>> In fact, almost any real-world, non-toy calculation would
>>> favor it.
>>>
>>> In any case, nobody should depend on the precision out that far being
>>> accurate or "reliable."
>>
>>
>> IEEE floating point has well-defined behaviour and there is absolutely
>> nothing wrong with code that delivers more accurate results just
>> because it is actually aware of the actual semantics of the operations
>> being carried out.
>
> You just made the case for Walter doing what he did. :)

No.