Copyright | (c) 2014 Bryan O'Sullivan |
---|---|

License | BSD3 |

Maintainer | bos@serpentine.com |

Stability | experimental |

Portability | portable |

Safe Haskell | None |

Language | Haskell98 |

Functions for summing floating point numbers more accurately than
the naive `sum`

function and its counterparts in the
`vector`

package and elsewhere.

When used with floating point numbers, in the worst case, the
`sum`

function accumulates numeric error at a rate
proportional to the number of values being summed. The algorithms
in this module implement different methods of /compensated
summation/, which reduce the accumulation of numeric error so that
it either grows much more slowly than the number of inputs
(e.g. logarithmically), or remains constant.

- class Summation s where
- sumVector :: (Vector v Double, Summation s) => (s -> Double) -> v Double -> Double
- data KBNSum = KBNSum !Double !Double
- kbn :: KBNSum -> Double
- data KB2Sum = KB2Sum !Double !Double !Double
- kb2 :: KB2Sum -> Double
- data KahanSum = KahanSum !Double !Double
- kahan :: KahanSum -> Double
- pairwiseSum :: Vector v Double => v Double -> Double

# Summation type class

class Summation s where Source #

A class for summation of floating point numbers.

sumVector :: (Vector v Double, Summation s) => (s -> Double) -> v Double -> Double Source #

*O(n)* Sum a vector of values.

## Usage

Most of these summation algorithms are intended to be used via the
`Summation`

typeclass interface. Explicit type annotations should
not be necessary, as the use of a function such as `kbn`

or `kb2`

to extract the final sum out of a `Summation`

instance gives the
compiler enough information to determine the precise type of
summation algorithm to use.

As an example, here is a (somewhat silly) function that manually computes the sum of elements in a list.

sillySumList :: [Double] -> Double sillySumList = loop`zero`

where loop s [] =`kbn`

s loop s (x:xs) =`seq`

s' loop s' xs where s' =`add`

s x

In most instances, you can simply use the much more general `sum`

function instead of writing a summation function by hand.

-- Avoid ambiguity around which sum function we are using. import Prelude hiding (sum) -- betterSumList :: [Double] -> Double betterSumList xs =`sum`

`kbn`

xs

# Kahan-Babuška-Neumaier summation

Kahan-Babuška-Neumaier summation. This is a little more
computationally costly than plain Kahan summation, but is *always*
at least as accurate.

# Order-2 Kahan-Babuška summation

Second-order Kahan-Babuška summation. This is more computationally costly than Kahan-Babuška-Neumaier summation, running at about a third the speed. Its advantage is that it can lose less precision (in admittedly obscure cases).

This method compensates for error in both the sum and the first-order compensation term, hence the use of "second order" in the name.

# Less desirable approaches

## Kahan summation

Kahan summation. This is the least accurate of the compensated
summation methods. In practice, it only beats naive summation for
inputs with large magnitude. Kahan summation can be *less*
accurate than naive summation for small-magnitude inputs.

This summation method is included for completeness. Its use is not
recommended. In practice, `KBNSum`

is both 30% faster and more
accurate.

## Pairwise summation

pairwiseSum :: Vector v Double => v Double -> Double Source #

*O(n)* Sum a vector of values using pairwise summation.

This approach is perhaps 10% faster than `KBNSum`

, but has poorer
bounds on its error growth. Instead of having roughly constant
error regardless of the size of the input vector, in the worst case
its accumulated error grows with *O(log n)*.

# References

- Kahan, W. (1965), Further remarks on reducing truncation
errors.
*Communications of the ACM*8(1):40. - Neumaier, A. (1974), Rundungsfehleranalyse einiger Verfahren zur
Summation endlicher Summen.
*Zeitschrift für Angewandte Mathematik und Mechanik*54:39–51. - Klein, A. (2006), A Generalized
Kahan-Babuška-Summation-Algorithm.
*Computing*76(3):279-293. - Higham, N.J. (1993), The accuracy of floating point
summation.
*SIAM Journal on Scientific Computing*14(4):783–799.