C language

In floating-point arithmetic, the order in which numbers are added matters. There is a famous algorithm originally due to Kahan (1965) called compensated summation that estimates the error produced when summing two floating-point numbers and accumulates this error estimate during the computation. It is based on the observation that the round-off error in the sum of two floating-point numbers can be estimated by subtracting one of the numbers from the sum.

Let a and b be two numbers such that |a| |b|. Their exact sum is s = a + b = s + e.

Let a = a₁ + a₂ and b = b₁ + b₂, where a₂ and b₁ are determined by the significance overlap of the numbers. Then s = a₁ + (a₂ + b₁).

Now s a = b₁,

and thus we can get an estimate of the error from

(s a) b = b₂ = e.

Generate an array of one million random numbers in single precision from a uniform distribution on [0, 1].

Cast the numbers as double precision and sum them in the order in which they were generated. We will take this as the reference solution, i.e., a highly accurate estimate of the exact solution.

Sum the numbers in the order in which they were generated and compute the error in the sum.

Sort the numbers in ascending order, sum them, and compute the error in the sum. Sum the numbers using compensated summation and compute the error in the sum.

Comment on the performance of the different methods for computing the sum.