> Are you saying that "fsum(seq) / len(seq)" is incorrect because on some
> older builds there is a rare possibility of an error of 1 unit in the last
> place?

Just to be clear, it's not that rare a possibility, and it's not restricted to older builds. I think Raymond is referring to the bug where on some machines that are still using the x87 for arithmetic, double rounding can lead to the fsum result being out (by potentially an arbitrary amount, not just 1 ulp). But even without that consideration, we'll still often be out by 1 ulp or more on typical new systems.

However, if I got my sums right (and assuming that fsum _is_ correctly rounded, and IEEE 754 arithmetic is in use, and the rounding mode hasn't been changed from its default of round-ties-to-even, and we're excluding corner cases like overflow, underflow, etc., etc.), the result of fsum(seq) / len(seq) can never be out by more than 1.5 ulps. That's better than NumPy can promise (even with its use of pairwise summation in some -- but not all -- cases), and should be good enough for almost any practical purpose.

I think I'd rather see the regular statistics.mean sacrifice the good-to-the-last-bit accuracy (and it's really not that much of a sacrifice: accurate to 1.5 ulps is a _very_ respectable compromise) than have a new function. I don't know how feasible that is, though, given all the type-handling in statistics.mean.

Bikeshedding time: the proposed fmean name seems a bit unfortunate in that it's reversing the sense of sum and fsum, where fsum is the correctly-rounded, slightly slower variant. My first assumption on seeing the name was that fmean was supposed to be to mean as fsum is to sum.