Precise fsum

[9il]

Blocked

Issues [DMD]

1. [ICE] https://issues.dlang.org/show_bug.cgi?id=14158 (workaround)
2. https://issues.dlang.org/show_bug.cgi?id=13485 (fixed)
3. https://issues.dlang.org/show_bug.cgi?id=13474 (workaround)

Note issue 13469 while testing:
https://issues.dlang.org/show_bug.cgi?id=13469 (workaround)

Example 0:

import std.numeric.summation;
import std.math, std.algorithm, std.range;
auto r1 = iota(  1, 501 ).map!(a => (-1.0)^^a/a);
auto r2 = iota(501, 1001).map!(a => (-1.0)^^a/a);
Summator!double s1 = 0, s2 = 0;
foreach(e; r1) s1.put(e); 
foreach(e; r2) s2 -= e;  //put(-e)
s1 -= s2;
assert(s1.sum == -0.69264743055982025);

Example 1:

import std.algorithm : map;
import std.numeric : fsum;
auto ar = [1, 1e100, 1, -1e100].map!(a => a*10000);
double r = 20000;
assert(r != ar.fsum!(Summation.Naive));
assert(r != ar.fsum!(Summation.Pairwise));
assert(r != ar.fsum!(Summation.Kahan));
assert(r == ar.fsum!(Summation.KBN)); //Fails with dmd -O -m32: DMD Optimizer BUG.   :(
assert(r == ar.fsum!(Summation.KB2));
assert(r == ar.fsum);  //Summation.Precise

Example 2:

//Summation.Precise is default

// Works like Python's fsum, 
// but  current implementation re-establish special
// value semantics across iterations (i.e. handling -Inf + Inf).

import std.numeric : fsum;
import std.math, std.algorithm, std.range;
auto ar = iota(1000).map!(a => 1.7.pow(a+1) - 1.7.pow(a)).chain([-(1.7.pow(1000))]);
assert(ar.fsum  ==  -1.0);
assert(ar.fsum!real  ==  -1.0);

[quickfur]

ping

[andralex]

May I convince some of the dmd core team to take a look at https://issues.dlang.org/show_bug.cgi?id=13485, https://issues.dlang.org/show_bug.cgi?id=13474, and https://issues.dlang.org/show_bug.cgi?id=13469

Perhaps aim the fixes for 2.067?

cc @WalterBright @donc @9rnsr @yebblies

[yebblies]

I can probably have a look at the optimizer bugs when I'm done with va_copy. Do we have any sort of timeline on 2.067?

[andralex]

Let's aim for Sunday March 1st, which means we need to have a prerelease by Sun Feb 15.

[yebblies]

This needs rebasing.

[9il]

Whoohooo! I am starting to rewrite this code now! Thanks!

[9il]

Issue 13485 fixed (Thanks to @yebblies).
Issue 13474 still blocks this PR.

[yebblies]

Would doing all the internal calculations explicitly with real on x86 work around the failures?

[9il]

Yes, for real it works well.
I will try to create internal representation with real for x86.

[9il]

There are will be some minimal precision loss for float and double versions of fsum on x87:
nextDown(s) <= r && nextUp(s) >= r
instead of s == r. I think it is acceptable.

[9il]

Reopened in #2991.