XSum

Fast Sum with error compensation

XSum - SUM with error compensation
The accuracy of the sum of floating point numbers is limited by the truncation error. E.g. SUM([1e16, 1, -1e16]) replies 0 instead of 1 and the error of SUM(RANDN(N, 1)) is about EPS*(N / 10).
Kahan, Knuth, Dekker, Ogita and Rump (and others) have derived some methods to reduce the influence of rounding errors, which are implemented here as fast C-Mex: XSum(RANDN(N, 1), 'Knuth') is exact to all 15 digits.

Y = XSum(X, N, Method)
INPUT:
X: Double array of any size.
N: Dimension to operate on.
Method: String: 'Double', 'Long', 'Kahan', 'Knuth', 'KnuthLong', 'Knuth2'.

OUTPUT:
Y: Double array, equivalent to SUM, but with compensated error depending
on the Method. The high-precision result is rounded to double precision.

METHODS: (speed and accuracy compared to SUM)
- Double: A single threaded implementation of Matlab's SUM. At least in Matlab 2008a to 2009b the results of the multi-threaded SUM can differ slightly from call to call. Equivalent accuracy. 1.1 to 2 times slower than SUM.
- Long: Accumulated in a 80 bit long double, if the compiler support this (e.g. LCC v3.8). 3.5 more valid digits, 2 times slower.
- Kahan: The local error is subtracted from the next element. 1 to 3 more valid digits, 2 to 9 times slower.
- Knuth: As if the sum is accumulated in a 128 bit float: about 15 more valid digits. 1.4 to 4 times slower. This is suitable for the most real world problems.
- Knuth2: 30 more valid digits as if it is accumulated in a 196 bit float. 2 to 8 times slower.
- KnuthLong: As Knuth, but using long doubles to get about 21 more valid digits, if supported by the compiler. 2.5 times slower.

COMPILATION: mex -O XSum.c
Tested: Matlab 6.5, 7.7, 7.8, WinXP, BCC5.5, LCC2.4/3.8, Open Watcom 1.8, MSVC++ 2008
TestXSum checks the validity, speed and accuracy after compiling (see screen shot).
Pre-compiled MEX files: http://www.n-simon.de/mex

References: Takeshi Ogita and Siegfried M. Rump and Shin'ichi Oishi: "Accurate Sum and Dot Product with Applications"
See also: INTLAB, S. M. Rump: http://www.ti3.tu-harburg.de/rump/intlab