How does Plus work on machine precision Real arguments?

I thought Kahan's summation method would make a nice example for students to use to think about round-off error [W. Kahan, Pracniques: Further Remarks on Reducing Truncation Errors, Commun. ACM 8  (1965), 40]. The method is available (I surmise) through the option Method -> "CompensatedSummation" of Total, but it's an easy program to write that is nicely laid out in the half-page article.

Below the vector x0 is a list whose elements are to be added. Note that s0 is the "exact" sum, that is, the sum of the exact numbers represented by the floating-point numbers in x0, rounded to machine precision. I consider it the target of the summation methods below.

SeedRandom;
x0 = 2 + RandomReal[1, 1000000];
s0 = N@Total[SetPrecision[x0, Infinity]];

Here is a norm (with sign) that measures the relative error in x - x0 in units of $MachineEpsilon. meps[x_, x0_, eps_:$MachineEpsilon] := ((x - x0)/x0)/eps;

Then here are some ways to add up the vector. The timings are not particularly important, but they're somewhat interesting. I unpack x0 before apply Plus, since it would be unpacked anyway by Apply (@@), in order that the unpacking not be included in the timing; the sum is the same we would get with Plus @@ x0.

{With[{x0 = DeveloperFromPackedArray[x0]},
sP = Plus @@ x0; // AbsoluteTiming],
sF = Fold[Plus, x0]; // AbsoluteTiming,
sT = Total[x0]; // AbsoluteTiming,
sC = Total[x0, Method -> "CompensatedSummation"]; // AbsoluteTiming
}[[All, 1]]
(*  {0.131215, 0.120635, 0.000571, 0.01399}  *)

The example was chosen to show that Plus does something nearly like compensated summation but not exactly the same.

meps[{sP, sF, sT, sC}, s0]
(*  Plus       Fold     Total     Compensated  Version
{-0.838972, 149.337, -56.2111, 0.}          10.4.1
{-0.838972, 149.337, -2.51692, 0.}          11.1.1  *)

Aside from being surprised that Plus is such a slow method, I got to wondering what Plus is doing. With some effort, one can find that the order of the elements in x0 can affect the result of Plus. The only results I can find have errors of either 1 or 0 in the last bit. They are compared to the results of iterative adding using Fold.

meps[{Plus @@ Sort@x0, Fold[Plus, Sort@x0]}, s0]
(*  {-0.838972, -33.5589}  *)

meps[{Plus @@ -Sort[-x0], Fold[Plus, -Sort[-x0]]}, s0]
(*  {-0.838972, 182.057}  *)

SeedRandom;
With[{x0 = RandomSample[x0]}, meps[{Plus @@ x0, Fold[Plus, x0]}, s0]]
(*  {0., -158.566}  *)

SeedRandom;
With[{x0 = RandomSample[x0]}, meps[{Plus @@ x0, Fold[Plus, x0]}, s0]]
(*  {-0.838972, -93.9648}  *)

Does anyone know how Plus works? It would be nice to be able to explain it.

• Might want to check on progressively larger examples to see if it scales as n log(n), that might indicate sorting under the hood. I remark that I do not know the answer offhand and probably will not have time to peek into the code to find out. – Daniel Lichtblau Jan 10 '16 at 23:31
• @george2079 Perhaps you mean extended-precision registers on the CPU/FPU are used to accumulate the result? That may be possible if the summand are chunked and the intermediate sums rounded to MP from time to time. If the whole computation is done with as little as 6 extra bits, the result rounds to the "exact" 53-bit result. – Michael E2 Jan 11 '16 at 17:19
• If it is a matter of extended precision registers, this is likely to vary based on operating system, whether Compile was used somewhere, and probably a few other incidentals (like, phase of the moon). – Daniel Lichtblau Jan 12 '16 at 15:56
• @AlexeyPopkov I did check out the Orderless attribute: Plus can give different results (of at most one ulp) for different orderings of the arguments (via RandomSample). If the Orderless attribute affects evaluation, it is not by the arguments being sorted into canonical order before evaluation. It makes sense to me that orderless numerical functions like Plus would go directly to evaluation on numeric input. Simple example: Plus[1., -1., $MachineEpsilon/2] vs. Plus[1.,$MachineEpsilon/2, -1.]. Orderless is more important for symbolic computations. – Michael E2 May 6 '17 at 20:07
• Interesting question. Incidentally you forgot to include Tr, my favorite method to sum a list. Including it in v10.1 I find it twice as fast as (plain) Total with the same error. – Mr.Wizard Aug 1 '17 at 9:14