| Title: | Accurate Floating Point Sums and Products |
|---|---|
| Description: | Most of the time floating point arithmetic does approximately the right thing. When adding sums or having products of numbers that greatly differ in magnitude, the floating point arithmetic may be incorrect. This package implements the Kahan (1965) sum <doi:10.1145/363707.363723>, Neumaier (1974) sum <doi:10.1002/zamm.19740540106>, pairwise-sum (adapted from 'NumPy', See Castaldo (2008) <doi:10.1137/070679946> for a discussion of accuracy), and arbitrary precision sum (adapted from the fsum in 'Python' ; Shewchuk (1997) <https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf>). In addition, products are changed to long double precision for accuracy, or changed into a log-sum for accuracy. |
| Authors: | Matthew Fidler [aut, cre, cph], Raymond Hettinger [cph, aut], Jonathan Shewchuk [cph, aut], Julian Taylor [cph, aut], Nathaniel Smith [cph, aut], NumPy Team [cph], Python Team [cph] |
| Maintainer: | Matthew Fidler <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.7 |
| Built: | 2024-08-25 16:20:55 UTC |
| Source: | https://github.com/nlmixr2/PreciseSums |
This method avoids loss of precision by tracking multiple intermediate partial sums. Based on python's math.fsum
fsum(numbers)fsum(numbers)
numbers |
A vector of numbers to sum. |
Sum of numbers without loss of precision
The algorithm's accuracy depends on IEEE-754 arithmetic guarantees and the typical case where the rounding mode is half-even. On some non-Windows builds, the underlying C library uses extended precision addition and may occasionally double-round an intermediate sum causing it to be off in its least significant bit.
Matthew Fidler (R implementation), Raymond Hettinger, Jonathan Shewchuk, Python Team
https://docs.python.org/2/library/math.html
https://code.activestate.com/recipes/393090/
https://github.com/python/cpython/blob/a0ce375e10b50f7606cb86b072fed7d8cd574fe7/Modules/mathmodule.c
Shewchuk, JR. (1996) Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates. http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 fsum(c(1,1e100,1,-1e100)) ## Gives 2.sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 fsum(c(1,1e100,1,-1e100)) ## Gives 2.
Using the Kahan method, take a more accurate sum
kahanSum(numbers)kahanSum(numbers)
numbers |
A vector of numbers to sum. |
Sum of numbers
https://en.wikipedia.org/wiki/Kahan_summation_algorithm
sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 kahanSum(c(1,1e100,1,-1e100)) ## Not accurate enough for the correct result. (still = 0)sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 kahanSum(c(1,1e100,1,-1e100)) ## Not accurate enough for the correct result. (still = 0)
Using the Neumaier method, take a more accurate sum
neumaierSum(numbers)neumaierSum(numbers)
numbers |
A vector of numbers to sum. |
Sum of numbers, a bit more accurate than kahanSum
https://en.wikipedia.org/wiki/Kahan_summation_algorithm
sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 neumaierSum(c(1,1e100,1,-1e100)) ## Gives 2sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 neumaierSum(c(1,1e100,1,-1e100)) ## Gives 2
This was taken by NumPy and adapted for use here. It is more accurate than a standard sum, but still has numerical issues. Its main benefit is that it is about the same amount of time as a standard time with the added accuracy.
pairwiseSum(numbers)pairwiseSum(numbers)
numbers |
A vector of numbers to sum. |
A sum of numbers with a rounding error of O(lg n) instead of O(n).
Matthew Fidler (R implementation), Julian Taylor, Nathaniel J Smith, and NumPy team. #' @examples sum(c(1,1e100,1,-1e100)) ## Should be 2, gives 0 pairwiseSum(c(1,1e100,1,-1e100)) ## Should be 2, still 0
https://github.com/numpy/numpy/pull/3685
Using PreciceSums's default method, take a product
psProd(numbers)psProd(numbers)
numbers |
A vector of numbers to sum. |
Product of numbers
prod blocksChoose the type of product to use in PreciceSums. These are used in the
PreciceSums prod blocks
psSetProd(type = c("long double", "double", "logify"))psSetProd(type = c("long double", "double", "logify"))
type |
Product to use for
|
nothing
Matthew L. Fidler
Choose the types of sums to use in PreciceSums. These are used in the
PreciceSums sum blocks and the psSum function
psSetSum(type = c("pairwise", "fsum", "kahan", "neumaier", "klein", "c"))psSetSum(type = c("pairwise", "fsum", "kahan", "neumaier", "klein", "c"))
type |
Sum type to use for
|
nothing
Matthew L. Fidler
Using PreciceSums's default method, take a sum
psSum(numbers)psSum(numbers)
numbers |
A vector of numbers to sum. |
Sum of numbers