The idea for my blog topic this week came from a question from Ron Kneusel about the loss of precision using IDL’s TOTAL function. IDL functions are often implemented in a way that makes them run fast, since processing speed is important. This is certainly true for the TOTAL function. If you are not using any of the keywords and additional arguments, then it is basically adding to the sum in the input array order. The exception is if you have enough elements in the array to trigger the multi-threaded code path, which I will ignore for now.

The issue with loss of precision is therefore best seen if you start with the largest term and keep adding smaller and smaller terms. At some point, the small term gets so small compared to the large sum that the additional terms get lost in round off.

Here is an example of a sum that we know should approach 128 (with a deviation in the 52 decimal place):

So, we are getting an error in the 13^th decimal place in this case. This array starts with the largest term and ends with the smallest term, so that is the worst case for the TOTAL implementation. If we reverse the order of the terms, we get:

The error gets 10 times smaller in this case. I decided to compare the algorithm that Ron suggested, which is called Kahan sum, with a divide-and-conquer scheme that I have previously used for execution on a massively parallel GPU (which runs fast with 10000’s of independent threads). The Kahan sum algorithm can be found on Wikipedia and the IDL code is pretty short, (but very slow):

The massively parallel algorithm will run much faster than the Kahan Sum, and the IDL code is listed here:

Here are the result with the decreasing magnitude terms:

So, in this case the error is similar to the best case sorted input to TOTAL. As an independent test, I also tried randomly ordering the terms, and both KahanSum and TOTAL_MP, are still consistent on the order of 5e-14:

My conclusion is that the TOTAL_MP example is just as accurate as the Kahan sum, and has the added advantage of allowing for massively parallel execution, as opposed to the highly sequential execution needed for the Kahan algorithm.