numerics - 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[2];
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 = Developer`FromPackedArray[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[0];
With[{x0 = RandomSample[x0]}, meps[{Plus @@ x0, Fold[Plus, x0]}, s0]]
(*  {0., -158.566}  *)

SeedRandom[1];
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.