# Last@Accumulate not giving same result as Total

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).

• The 250 limit has to do with whether Table returns a packed array. Try arr = DeveloperToPackedArray[arr]. I would have guessed that the different result is simply due to the summation being done in a different order, but I don't understand how this relates to the array being packed. Sep 21 '16 at 4:44
• ... and use Method->"CompensatedSummation" for both and compare...
– ciao
Sep 21 '16 at 4:47
• Jul 8 '19 at 3:38
• I get a similar plot with tables over 250 in length in V12. Jul 8 '19 at 3:45
• Strongly related: "Why is Plus so much slower than Total?" Sep 24 '19 at 11:49

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.

• Please see my addition to your answer. What I really don't understand is why method or order of summation can matter when all the numbers are identical? Sep 23 '19 at 13:38
• @AlexeyPopkov It would be strange if the order of summation mattered when the numbers are all identical. But I don't think it's strange that the method of summation matters. Real numbers cannot be represented exactly in computers. As a result, the error builds up when you add many numbers together. My understanding is that there exists different methods for decreasing the error in the end result. Since the methods approach this error fixing differently, with different degrees of success, I imagine, it seems natural that the end outcome should be different for different methods. Sep 24 '19 at 6:57
• @AlexeyPopkov I would like to ask this on CS Stack Exchange since it is outside of my expertise. It sounds like it's too much though, I agree. Sep 24 '19 at 9:00
• I updated my answer with additional observations which potentially can help to figure out what happens. Sep 24 '19 at 13:04

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 DeveloperPackedArrayQ 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 = DeveloperToPackedArray@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[DeveloperToPackedArray@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[DeveloperToPackedArray@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.