Exactly Rounded Sums, Means and Dot Products

Sums of double-precision numbers computed exactly and rounded once, so the result is the correctly rounded true sum and does not depend on the order of the elements, the platform, or whether R was built with long doubles. Built on the small and large superaccumulators of Neal (2015) . Provides exactly rounded sums, means, cumulative sums, dot products and norms, row, column and grouped sums, accumulators that can be built in pieces and merged across processes, and a registered 'C' interface for other packages.

zuexact

zuexact computes sums of doubles exactly and rounds once, at the end. The result is the correctly rounded true sum. It does not depend on the order of the elements, on the platform, or on how R was built. It is built on Radford Neal's xsum superaccumulators.

Base R's sum() accumulates in long double. That is 80-bit on x86-64 Linux and Windows, but the same as double on Apple Silicon and on R built without long doubles. So the same call gives different answers on different machines:

sum(c(1e16, 1, -1e16))
#> 1 on x86-64 Linux and Windows
#> 0 on Apple Silicon

The exact answer is 1, everywhere, and xsum() returns it on every platform:

library(zuexact)
xsum(c(1e16, 1, -1e16))
#> [1] 1
xsum(c(1e20, 1, -1e20))      # wrong on every platform with sum()
#> [1] 1
xsum_mean(c(1e15, -1e15, 0.1))
#> [1] 0.03333333

Sums in pieces, and across processes

An accumulator holds an exact sum in progress. It is an ordinary R value, so an accumulator built on a worker comes back to the caller and merges there, in any order, with the same bits as summing everything at once:

chunks <- split(x, rep(1:8, length.out = length(x)))
accs <- parallel::mclapply(chunks, xsum_acc)
xsum_value(Reduce(xsum_merge, accs)) == xsum(x)
#> [1] TRUE

The same works for Arrow record batches, DuckDB result chunks or files read in pieces: acc <- xsum_add(acc, chunk) for each one, then xsum_value(acc). xsum_merge(a, xsum_negate(b)) is the exact difference of two sums.

What else is in the package

• xsum_mean(), xsum_weighted_mean(): exact means.
• xsum_dot(), xsum_sqnorm(), xsum_norm(): dot products and norms with each product's rounding error included.
• xsum_cumsum(), xsum_rowsums(), xsum_colsums(), xsum_by(): exact versions of cumsum(), rowSums(), colSums() and rowsum().
• A registered C interface, <zuexact.h>, for other packages' compiled code.

vignette("exact-summation") explains the details. For compensated summation (Kahan, Neumaier, pairwise), which is faster but not exact, see PreciseSums.

Installation

Install the released version from CRAN with:

install.packages("zuexact")

or the development version from GitHub with:

# install.packages("pak")
pak::pak("pedrobtz/zuexact")