math-functions-0.3.0.1: Collection of tools for numeric computations

Numeric.Sum

Description

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.

Synopsis

Summation type class

class Summation s where Source #

A class for summation of floating point numbers.

Minimal complete definition

Methods

zero :: s Source #

The identity for summation.

add :: s -> Double -> s Source #

Add a value to a sum.

sum :: Foldable f => (s -> Double) -> f Double -> Double Source #

Sum a collection of values.

Example: foo = sum kbn [1,2,3]

Instances
 Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (Double -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KahanSum -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KBNSum -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KB2Sum -> Double) -> f Double -> Double Source #

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

data KBNSum Source #

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

Constructors

 KBNSum !Double !Double
Instances
 Source # Instance detailsDefined in Numeric.Sum Methods(==) :: KBNSum -> KBNSum -> Bool #(/=) :: KBNSum -> KBNSum -> Bool # Source # Instance detailsDefined in Numeric.Sum Methodsgfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b) -> (forall g. g -> c g) -> KBNSum -> c KBNSum #gunfold :: (forall b r. Data b => c (b -> r) -> c r) -> (forall r. r -> c r) -> Constr -> c KBNSum #dataCast1 :: Typeable t => (forall d. Data d => c (t d)) -> Maybe (c KBNSum) #dataCast2 :: Typeable t => (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c KBNSum) #gmapT :: (forall b. Data b => b -> b) -> KBNSum -> KBNSum #gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> KBNSum -> r #gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> KBNSum -> r #gmapQ :: (forall d. Data d => d -> u) -> KBNSum -> [u] #gmapQi :: Int -> (forall d. Data d => d -> u) -> KBNSum -> u #gmapM :: Monad m => (forall d. Data d => d -> m d) -> KBNSum -> m KBNSum #gmapMp :: MonadPlus m => (forall d. Data d => d -> m d) -> KBNSum -> m KBNSum #gmapMo :: MonadPlus m => (forall d. Data d => d -> m d) -> KBNSum -> m KBNSum # Source # Instance detailsDefined in Numeric.Sum MethodsshowsPrec :: Int -> KBNSum -> ShowS #showList :: [KBNSum] -> ShowS # Source # Instance detailsDefined in Numeric.Sum Methods(<>) :: KBNSum -> KBNSum -> KBNSum #stimes :: Integral b => b -> KBNSum -> KBNSum # Source # Instance detailsDefined in Numeric.Sum Methodsmconcat :: [KBNSum] -> KBNSum # Source # Instance detailsDefined in Numeric.Sum Methodsrnf :: KBNSum -> () # Source # Instance detailsDefined in Numeric.Sum Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KBNSum -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeIndexM :: Monad m => Vector KBNSum -> Int -> m KBNSum #basicUnsafeCopy :: PrimMonad m => Mutable Vector (PrimState m) KBNSum -> Vector KBNSum -> m () #elemseq :: Vector KBNSum -> KBNSum -> b -> b # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeNew :: PrimMonad m => Int -> m (MVector (PrimState m) KBNSum) #basicInitialize :: PrimMonad m => MVector (PrimState m) KBNSum -> m () #basicUnsafeReplicate :: PrimMonad m => Int -> KBNSum -> m (MVector (PrimState m) KBNSum) #basicUnsafeRead :: PrimMonad m => MVector (PrimState m) KBNSum -> Int -> m KBNSum #basicUnsafeWrite :: PrimMonad m => MVector (PrimState m) KBNSum -> Int -> KBNSum -> m () #basicClear :: PrimMonad m => MVector (PrimState m) KBNSum -> m () #basicSet :: PrimMonad m => MVector (PrimState m) KBNSum -> KBNSum -> m () #basicUnsafeCopy :: PrimMonad m => MVector (PrimState m) KBNSum -> MVector (PrimState m) KBNSum -> m () #basicUnsafeMove :: PrimMonad m => MVector (PrimState m) KBNSum -> MVector (PrimState m) KBNSum -> m () #basicUnsafeGrow :: PrimMonad m => MVector (PrimState m) KBNSum -> Int -> m (MVector (PrimState m) KBNSum) # data Vector KBNSum Source # Instance detailsDefined in Numeric.Sum data Vector KBNSum = V_KBNSum (Vector (Double, Double)) data MVector s KBNSum Source # Instance detailsDefined in Numeric.Sum data MVector s KBNSum = MV_KBNSum (MVector s (Double, Double))

Return the result of a Kahan-Babuška-Neumaier sum.

Order-2 Kahan-Babuška summation

data KB2Sum Source #

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.

Constructors

 KB2Sum !Double !Double !Double
Instances
 Source # Instance detailsDefined in Numeric.Sum Methods(==) :: KB2Sum -> KB2Sum -> Bool #(/=) :: KB2Sum -> KB2Sum -> Bool # Source # Instance detailsDefined in Numeric.Sum Methodsgfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b) -> (forall g. g -> c g) -> KB2Sum -> c KB2Sum #gunfold :: (forall b r. Data b => c (b -> r) -> c r) -> (forall r. r -> c r) -> Constr -> c KB2Sum #dataCast1 :: Typeable t => (forall d. Data d => c (t d)) -> Maybe (c KB2Sum) #dataCast2 :: Typeable t => (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c KB2Sum) #gmapT :: (forall b. Data b => b -> b) -> KB2Sum -> KB2Sum #gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> KB2Sum -> r #gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> KB2Sum -> r #gmapQ :: (forall d. Data d => d -> u) -> KB2Sum -> [u] #gmapQi :: Int -> (forall d. Data d => d -> u) -> KB2Sum -> u #gmapM :: Monad m => (forall d. Data d => d -> m d) -> KB2Sum -> m KB2Sum #gmapMp :: MonadPlus m => (forall d. Data d => d -> m d) -> KB2Sum -> m KB2Sum #gmapMo :: MonadPlus m => (forall d. Data d => d -> m d) -> KB2Sum -> m KB2Sum # Source # Instance detailsDefined in Numeric.Sum MethodsshowsPrec :: Int -> KB2Sum -> ShowS #showList :: [KB2Sum] -> ShowS # Source # Instance detailsDefined in Numeric.Sum Methods(<>) :: KB2Sum -> KB2Sum -> KB2Sum #stimes :: Integral b => b -> KB2Sum -> KB2Sum # Source # Instance detailsDefined in Numeric.Sum Methodsmconcat :: [KB2Sum] -> KB2Sum # Source # Instance detailsDefined in Numeric.Sum Methodsrnf :: KB2Sum -> () # Source # Instance detailsDefined in Numeric.Sum Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KB2Sum -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeIndexM :: Monad m => Vector KB2Sum -> Int -> m KB2Sum #basicUnsafeCopy :: PrimMonad m => Mutable Vector (PrimState m) KB2Sum -> Vector KB2Sum -> m () #elemseq :: Vector KB2Sum -> KB2Sum -> b -> b # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeNew :: PrimMonad m => Int -> m (MVector (PrimState m) KB2Sum) #basicInitialize :: PrimMonad m => MVector (PrimState m) KB2Sum -> m () #basicUnsafeReplicate :: PrimMonad m => Int -> KB2Sum -> m (MVector (PrimState m) KB2Sum) #basicUnsafeRead :: PrimMonad m => MVector (PrimState m) KB2Sum -> Int -> m KB2Sum #basicUnsafeWrite :: PrimMonad m => MVector (PrimState m) KB2Sum -> Int -> KB2Sum -> m () #basicClear :: PrimMonad m => MVector (PrimState m) KB2Sum -> m () #basicSet :: PrimMonad m => MVector (PrimState m) KB2Sum -> KB2Sum -> m () #basicUnsafeCopy :: PrimMonad m => MVector (PrimState m) KB2Sum -> MVector (PrimState m) KB2Sum -> m () #basicUnsafeMove :: PrimMonad m => MVector (PrimState m) KB2Sum -> MVector (PrimState m) KB2Sum -> m () #basicUnsafeGrow :: PrimMonad m => MVector (PrimState m) KB2Sum -> Int -> m (MVector (PrimState m) KB2Sum) # data Vector KB2Sum Source # Instance detailsDefined in Numeric.Sum data Vector KB2Sum = V_KB2Sum (Vector (Double, Double, Double)) data MVector s KB2Sum Source # Instance detailsDefined in Numeric.Sum data MVector s KB2Sum = MV_KB2Sum (MVector s (Double, Double, Double))

Return the result of an order-2 Kahan-Babuška sum.

Less desirable approaches

Kahan summation

data KahanSum Source #

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.

Constructors

 KahanSum !Double !Double
Instances
 Source # Instance detailsDefined in Numeric.Sum Methods Source # Instance detailsDefined in Numeric.Sum Methodsgfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b) -> (forall g. g -> c g) -> KahanSum -> c KahanSum #gunfold :: (forall b r. Data b => c (b -> r) -> c r) -> (forall r. r -> c r) -> Constr -> c KahanSum #dataCast1 :: Typeable t => (forall d. Data d => c (t d)) -> Maybe (c KahanSum) #dataCast2 :: Typeable t => (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c KahanSum) #gmapT :: (forall b. Data b => b -> b) -> KahanSum -> KahanSum #gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> KahanSum -> r #gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> KahanSum -> r #gmapQ :: (forall d. Data d => d -> u) -> KahanSum -> [u] #gmapQi :: Int -> (forall d. Data d => d -> u) -> KahanSum -> u #gmapM :: Monad m => (forall d. Data d => d -> m d) -> KahanSum -> m KahanSum #gmapMp :: MonadPlus m => (forall d. Data d => d -> m d) -> KahanSum -> m KahanSum #gmapMo :: MonadPlus m => (forall d. Data d => d -> m d) -> KahanSum -> m KahanSum # Source # Instance detailsDefined in Numeric.Sum MethodsshowList :: [KahanSum] -> ShowS # Source # Instance detailsDefined in Numeric.Sum Methodsstimes :: Integral b => b -> KahanSum -> KahanSum # Source # Instance detailsDefined in Numeric.Sum Methodsmconcat :: [KahanSum] -> KahanSum # Source # Instance detailsDefined in Numeric.Sum Methodsrnf :: KahanSum -> () # Source # Instance detailsDefined in Numeric.Sum Source # Instance detailsDefined in Numeric.Sum Methodssum :: Foldable f => (KahanSum -> Double) -> f Double -> Double Source # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeIndexM :: Monad m => Vector KahanSum -> Int -> m KahanSum #basicUnsafeCopy :: PrimMonad m => Mutable Vector (PrimState m) KahanSum -> Vector KahanSum -> m () #elemseq :: Vector KahanSum -> KahanSum -> b -> b # Source # Instance detailsDefined in Numeric.Sum MethodsbasicUnsafeNew :: PrimMonad m => Int -> m (MVector (PrimState m) KahanSum) #basicInitialize :: PrimMonad m => MVector (PrimState m) KahanSum -> m () #basicUnsafeReplicate :: PrimMonad m => Int -> KahanSum -> m (MVector (PrimState m) KahanSum) #basicUnsafeRead :: PrimMonad m => MVector (PrimState m) KahanSum -> Int -> m KahanSum #basicUnsafeWrite :: PrimMonad m => MVector (PrimState m) KahanSum -> Int -> KahanSum -> m () #basicClear :: PrimMonad m => MVector (PrimState m) KahanSum -> m () #basicSet :: PrimMonad m => MVector (PrimState m) KahanSum -> KahanSum -> m () #basicUnsafeCopy :: PrimMonad m => MVector (PrimState m) KahanSum -> MVector (PrimState m) KahanSum -> m () #basicUnsafeMove :: PrimMonad m => MVector (PrimState m) KahanSum -> MVector (PrimState m) KahanSum -> m () #basicUnsafeGrow :: PrimMonad m => MVector (PrimState m) KahanSum -> Int -> m (MVector (PrimState m) KahanSum) # data Vector KahanSum Source # Instance detailsDefined in Numeric.Sum data MVector s KahanSum Source # Instance detailsDefined in Numeric.Sum data MVector s KahanSum = MV_KahanSum (MVector s (Double, Double))

Return the result of a Kahan sum.

Pairwise summation

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.