You should be okay with a sum but interval arithmetic can blow up in a way that is not correlated with the numerical accuracy of the algorithm. Newton's method, for example, is strongly convergent but produces diverging intervals.

Kahan summation is better than naive summation, but not perfect. The stable_mean algorithm in the article is not perfect either. It just has more precision than the naive one, but it more or less has the same precision as twice-as-precise floats.

>>> data = (1_000_000.1, 1.1, -1_000_000.1, -1.1) * 25_000

>>> sum(data) / len(data) == -2.3283153183228934e-16

>>> sum(sorted(data)) / len(data) == 2.7122441679239272e-11

>>> stable_mean(data) == -2.571225228716577e-13

>>> float(sum(numpy.array(data, dtype=numpy.float128)) / len(data)) == 2.2648549702353195e-19

Note that in python there is a math.fsum function which return a perfectly rounded summation of floats. There are multiple different techniques for that, but they are more complex than Kahan summation.

>>> math.fsum(data) / len(data) == 0.0

Interestingly enough I was given this problem (well, numerically stable online variance, but close enough) a while ago as a take-home job interview and I ended up researching this quite a bit.

It has been a while, but if I remember correctly pairwise sum was in practice as good or better as Kahan summation (and both more stable than Walford's) but faster.

But by far the fastest and still very stable solution was to use the legacy 80 bit x87 fpu as accumulator. Which is not surprising as Kahan explicitly added 80 bit precision to x87 for intermediate computations.

Or can we sum more accurately by taking the stable_mean and multiplying by n?

Ruby seems to handle this out of the box:

  >> RUBY_VERSION
  ⇒  "2.7.2"
  >> data = [1_000_000_000.1, 1.1] * 50_000 ; nil
  ⇒  nil
  >> data.sum
  ⇒  50000000060000.0
  >> data.sum / data.length
  ⇒  500000000.6
  >> data.sort.sum / data.length
  ⇒  500000000.6
  >> data.sort.reverse.sum / data.length
  ⇒  500000000.6

I guess the question is does it still use floats internally or use something like BigNum and handle the precision issue (and thus more memory/cpu)

Not sure what this could possibly have to do with floating point numbers, nor am I aware of any function that does what you're describing. Is this some kind of JavaScript quirky behavior because nothing about floating point numbers has anything to do with coercing strings.

I've gone back and forth over the years. Having implemented some of the lowest-level floating point routines, taking care to handle all the corner cases, it definitely is a chore. But is there an alternative?

After all of this, I still think that IEEE 754 is pretty darn good. The fact that it is available and fast in pretty much all hardware is a major selling point.

That said, I absolutely hate -0 (minus zero). It's a stupid design wart that I've never seen any good use case for. But it is observable and therefore cannot be ignored, leading to even more icky special cases.

Don't blame the floating point numbers.

Just use a programming language that doesn't coerce strings into numbers. Or use a comparison function/operator that explicitely compares arguments as strings and not as numbers.

That isn't floating points fault, but whoever wrote the buggy parser..

Floating point is overall pretty good, sometimes I wish there was a version without some many ways to store NAN(for 16 bit floats it is so wasteful), and maybe versions of sqrt/div/rcp/rsqrt that would clamp invalid inputs to 0, but other than that it is hard to see how you make any major improvements given the fixed # of bits.

Actually one other possible improvement would be a way to request the hardware to dither when rounding down to the fixed bit count, although I'm not sure how it could do anything besides pure white noise, but better than nothing..

This is exactly what happens when you sort the number ascending and compute the sum starting from the smallest. Like in the article.

Except this doesn't work if you must handle negative numbers. So what matters is not the smallest numbers but the number nearer from zero.

> the only real question is how much error you are willing to tolerate

So, you are, in fact, gonna need an algorithm for summation that controls that error, if you're not willing to tolerate very much.

Well, his final example suggests otherwise, but I kind of suspect that there must be something else going on there (that I'm too lazy to try to track down). Still, while fixing most floating point errors seem like rearranging the deck chairs on the titanic, if there's one algorithm that's 99.999% accurate and another that's only 99.9% accurate, might as well use and be aware of the more accurate one.