Last@Accumulate not giving same result as Total

Question

If I run the code:

arr = Table[0.979, 249];
Last@Accumulate[arr] - Total[arr]

on Mathematica 11.0 (Windows), I get the result 1.53477*10^-12. I should instead be getting a result of 0, since Accumulate's last element should be the total of the array.

Also, this only happens when the array's length length is less than 250. The difference between the Total and Acumulate also seems to vary depending on the value of the elements inside the array, and for some values, Total and Accumulate will give the same result.

Something such as ListPlot@Table[(Accumulate[Table[a, 249]] - Total[Table[a, 249]])[[-1]], {a, 0, 1, 0.1}] shows an intresting graph of the difference between the two for various values:

My questions are:

• Why is this happening? Is this a bug?
• Why doesn't this happen for arrays larger than 250?
• I need Accumulate to give me the same result as Total, so I'm currently zero-padding my arrays to 250 as a workaround. Is there a better soultion instead?

---

Edit: Even though Last@Accumulate[arr] - Total[arr] is 1.53477*10^-12, Table[Last@Accumulate[arr] - Total[arr], 250] gives an array of 0s (Table[1.53477*10^-12, 250] does not).

Answer 1

Let's turn the comments into an answer. The issue seemingly arises because Accumulate and Total use different summation methods by default. As ciao pointed out,

arr = Table[0.979, 249];
Last@Accumulate[arr, Method -> "CompensatedSummation"] - 
 Total[arr, Method -> "CompensatedSummation"]

0.

I would also note that the difference that OP found also appears in other places, for example:

{
 Last@Accumulate[arr] - Total[arr],
 Plus @@ arr - Total[arr],
 Fold[Plus, arr] - Total[arr]
 }

{1.53477*10^-12, 0., 1.53477*10^-12}

Hence, the workaround seems to be to specify the summation method in Accumulate/Total. The dependence on the array size seems to have disappeared in V12 according to MichaelE2's comment.

Answer 2

^{(An addition to the answer by C. E.)}

It is interesting that behavior of Total in Mathematica 12 indeed depends on the length of the list according to Trace. Let us compare:

arr = Table[0.979, 249];
Trace[Total[arr], TraceInternal -> True] /. {arr -> "arr", 
  p : HoldForm[Plus[__]] :> Shallow[p, 3]}

{{arr, "arr"}, Total["arr"], {$MessageList={}, {}}, {Plus@@"arr", Plus[<<249>>], 243.771}, {$MessageList, {}}, 243.771}

arr = Table[0.979, 250];
Trace[Total[arr], TraceInternal -> True] /. {arr -> "arr", 
  p : HoldForm[Plus[__]] :> Shallow[p, 3]}    

{{arr, "arr"}, Total["arr"], 244.75}

Using Developer`PackedArrayQ we find that in the first case the array generated by Table isn't packed, while in the second case it is a packed array. Let us compare the results for the same unpacked and packed array:

arr = Table[0.979, 249];
arrPacked = Developer`ToPackedArray@arr;
{Plus @@ arr - Total[arr], Plus @@ arrPacked - Total[arrPacked], 
 Plus @@ arr - Total[arrPacked]}

{0., 1.13687*10^-13, 1.13687*10^-13}

We see that the result of such a simple summation of identical numbers by Total strongly depends on whether the array is packed or not.

Comparison with other summation methods reveals that we have at least 4 built-in methods giving different results:

arr = Table[0.979, 249];
sum = 0; Do[sum += 0.979, 249];
{Fold[Plus, arr],
 Plus @@ arr,
 Total[arr, Method -> "CompensatedSummation"],
 Total[Developer`ToPackedArray@arr]} - sum

{0., -1.53477*10^-12, -1.50635*10^-12, -1.64846*10^-12}

It is interesting to compare results of summations with multiplication:

{Fold[Plus, arr],
 Plus @@ arr,
 Total[arr, Method -> "CompensatedSummation"],
 Total[Developer`ToPackedArray@arr]} - 0.979*249

 {1.50635*10^-12, -2.84217*10^-14, 0., -1.42109*10^-13}

Surprisingly, Method -> "CompensatedSummation" produces the same result as multiplication!

Obviously the order of summation cannot matter here because all the numbers are identical. What can matter, for example, is the size of a chunk when the summation is performed in chunks:

sums = Table[N@ArrayReshape[arr, Table[n, Ceiling@N@Log[n, 249]]]
    //. {x__Real} :> Plus[x], {n, 2, 50}];
sums - Plus @@ arr
% // Abs // Max

{5.68434*10^-14, 5.68434*10^-14, 5.68434*10^-14, 2.84217*10^-14, 5.68434*10^-14, 
 2.84217*10^-14, 5.68434*10^-14, 5.68434*10^-14, 2.84217*10^-14, 0., 5.68434*10^-14, 
 2.84217*10^-14, 2.84217*10^-14, 2.84217*10^-14, 5.68434*10^-14, 5.68434*10^-14, 
 5.68434*10^-14, 2.84217*10^-14, 2.84217*10^-14, 0., 0., 0., 5.68434*10^-14, 
 5.68434*10^-14, 5.68434*10^-14, 5.68434*10^-14, 2.84217*10^-14, 2.84217*10^-14, 
 2.84217*10^-14, 2.84217*10^-14, 5.68434*10^-14, 5.68434*10^-14, 5.68434*10^-14, 
 2.84217*10^-14, 5.68434*10^-14, 5.68434*10^-14, 2.84217*10^-14, 0., 2.84217*10^-14, 
 2.84217*10^-14, 0., 0., 0., 0., 0., 0., 2.84217*10^-14, 5.68434*10^-14, 5.68434*10^-14}

5.68434*10^-14

These differences doesn't explain however the discrepancy of size 10^-12 which we get with built-in methods.

Let us compare with one-by-one summation in chunks:

sums = Table[N@ArrayReshape[arr, Table[n, Ceiling@N@Log[n, 249]]]
    //. l : {__Real} :> Fold[Plus, l], {n, 2, 249}];
sums - Plus @@ arr // Abs // Max
ListPlot[sums - Plus @@ arr]

1.53477*10^-12

The linear growth of the discrepancy starts at n = 131. It isn't clear for me from where comes this digit: "FoldCompileLength" is by default set to 100 and changing this option doesn't affect the plot.