Imprecise summation for Float32

[hsgg]

This package can have less precision than the standard sum() on julia-1.7.3. For example:

julia> using KahanSummation

julia> x = rand(Float32, 10_000_000) .+ 1f5;

julia> sum_kbn(x)
9.9479867f11

julia> sum(x)
1.00000504f12

Is this expected?

[ararslan]

How do you mean less precise?

If you mean f11 vs. f12, note that 9.9479867f11 === 0.99479867f12. (Typically in scientific notation you change the exponent rather than have a leading zero.)

[hsgg]

How do you mean less precise?

Sorry, what I mean is that the correct result is >1f12, and sum() gets it pretty well.

Actually, a better test case is without randomness. For example,

julia> x = fill(1f5, 10_000_000);

julia> sum_kbn(x)
9.9477704f11

julia> sum(x)
1.0f12

julia> 10_000_000 * 1f5
1.0f12

The error for sum_kbn() is ~0.5%.

[ararslan]

Interestingly, the original Kahan-Babuška algorithm gets it right but Neumaier's extension (i.e. sum_kbn) gets it wrong: the running compensation term ends up overcompensating.

julia> function sum_kb(x)
           s = c = zero(eltype(x))
           for xi in x
               y = xi - c
               t = s + y
               c = (t - s) - y
               s = t
           end
           return s
       end
sum_kb (generic function with 1 method)

julia> sum_kb(x)
1.0f12