Numeric.Compensated

This module provides a fairly extensive API for compensated floating point arithmetic based on Knuth's error free transformation, various algorithms by Ogita, Rump and Oishi, Hida, Li and Bailey, Kahan summation, etc. with custom compensated arithmetic circuits to do multiplication, division, etc. of compensated numbers.

In general if a has x bits of significand, Compensated a gives you twice that. You can iterate this construction for arbitrary precision.

References:

• http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf
• http://www.ti3.tuhh.de/paper/rump/OgRuOi05.pdf
• Donald Knuth's "The Art of Computer Programming, Volume 2: Seminumerical Algorithms"
• http://en.wikipedia.org/wiki/Kahan_summation_algorithm

Documentation

class (RealFrac a, Precise a, Floating a) => Compensable a where

Associated Types

data Compensated a

This provides a numeric data type with effectively doubled precision by using Knuth's error free transform and a number of custom compensated arithmetic circuits.

This construction can be iterated, doubling precision each time.

>>> round (Prelude.product [2..100] :: Compensated (Compensated (Compensated Double)))
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

>>> Prelude.product [2..100]
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Methods

with :: Compensable a => Compensated a -> (a -> a -> r) -> r

This extracts both the primal and residual components of a Compensated number.

compensated :: Compensable a => a -> a -> Compensated a

Used internally to construct compensated values that satisfy our residual contract.

When in doubt, use add a b compensated instead of compensated a b

magic :: a

This magic number is used to split the significand in half, so we can multiply them separately without losing precision in times.

Instances

Compensable Double

Compensable Float

Compensable a => Compensable (Compensated a)

_Compensated :: Compensable a => Iso' (Compensated a) (a, a)

This provides the isomorphism between the compact representation we store these in internally and the naive pair of the primal and residual components.

type Overcompensated a = Compensated (Compensated a)

primal :: Compensable a => Lens' (Compensated a) a

This Lens lets us edit the primal directly, leaving the residual untouched.

residual :: Compensable a => Lens' (Compensated a) a

This Lens lets us edit the residual directly, leaving the primal untouched.

uncompensated :: Compensable a => Compensated a -> a

Extract the primal component of a compensated value, when and if compensation is no longer required.

fadd :: Num a => a -> a -> (a -> a -> r) -> r

fadd a b k computes k x y such that

x + y = a + b
x = fl(a + b)

but only under the assumption that abs a >= abs b. If you aren't sure, use add.

Which is to say that x is the floating point image of (a + b) and y stores the residual error term.

lifting scalars

add :: Num a => a -> a -> (a -> a -> r) -> r

add a b k computes k x y such that

x + y = a + b
x = fl(a + b)

Which is to say that x is the floating point image of (a + b) and y stores the residual error term.

times :: Compensable a => a -> a -> (a -> a -> r) -> r

times a b k computes k x y such that

x + y = a * b
x = fl(a * b)

Which is to say that x is the floating point image of (a * b) and y stores the residual error term.

This could be nicer if we had access to a hardware fused multiply-add.

squared :: Compensable a => a -> (a -> a -> r) -> r

squared a k computes k x y such that

x + y = a * a
x = fl(a * a)

Which is to say that x is the floating point image of (a * a) and y stores the residual error term.

divide :: Compensable a => a -> a -> (a -> a -> r) -> r

split :: Compensable a => a -> (a -> a -> r) -> r

error-free split of a floating point number into two parts.

Note: these parts do not satisfy the compensated contract

kahan :: (Foldable f, Compensable a) => f a -> Compensated a

Perform Kahan summation over a list.

(+^) :: Compensable a => a -> Compensated a -> Compensated a

Calculate a scalar + compensated sum with Kahan summation.

(*^) :: Compensable a => a -> Compensated a -> Compensated a

Compute a * Compensated a

compensated operators

square :: Compensable a => Compensated a -> Compensated a

Calculate a fast square of a compensated number.