Basic Sum Operations without Round-off Errors

Functions for performing basic sum operations without round-off errors

sumexact(..., na.rm = FALSE)
  cumsumexact(x)

Arguments

x

numeric vector

...

numeric vector(s), numbers or other objects to be summed

na.rm

logical. Should missing values be removed?

Details

All three functions use full precision summation using multiple doubles for intermediate values. The sum of numbers x & y is a=x+y with error term b=error(a+b). That way a+b is equal exactly x+y, so sum of 2 numbers is stored as 2 or fewer values, which when added would under-flow. By extension sum of n numbers is calculated with intermediate results stored as array of numbers that can not be added without introducing an error. Only final result is converted to a single number

Value

Function sumexact returns single number. Function cumsumexact returns vector of the same length as x.

References

Round-off error correction is based on: Shewchuk, Jonathan, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates

McCullough, D.B., (1998) Assessing the Reliability of Statistical Software, Part I, The American Statistician, Vol. 52 No.

McCullough, D.B., (1999) Assessing the Reliability of Statistical Software, Part II, The American Statistician, Vol. 53 No 2

NIST Statistical Reference Datasets (StRD) website

Author

Jarek Tuszynski (SAIC) jaroslaw.w.tuszynski@saic.com based on code by Vadim Ogranovich, which is based on algorithms described in references, pointed out by Gabor Grothendieck.

See also

• sum is faster but not error-save version of sumexact

• cumsum is equivalent to cumsumexact

Examples

  x = c(1, 1e20, 1e40, -1e40, -1e20, -1)
  a = sum(x);         print(a)
#> [1] -1e+20
  b = sumexact(x);    print(b)
#> [1] 0
  stopifnot(b==0)
  a = cumsum(x);      print(a)
#> [1]  1e+00  1e+20  1e+40  0e+00 -1e+20 -1e+20
  b = cumsumexact(x); print(b)
#> [1] 1e+00 1e+20 1e+40 1e+20 1e+00 0e+00
  stopifnot(b[6]==0)