Portability | non-portable |
---|---|

Stability | experimental |

Maintainer | Edward Kmett <ekmett@gmail.com> |

Safe Haskell | None |

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

- class (RealFrac a, Precise a, Floating a) => Compensable a where
- data Compensated a
- with :: Compensable a => Compensated a -> (a -> a -> r) -> r
- compensated :: Compensable a => a -> a -> Compensated a
- magic :: a

- _Compensated :: Compensable a => Iso' (Compensated a) (a, a)
- type Overcompensated a = Compensated (Compensated a)
- primal :: Compensable a => Lens' (Compensated a) a
- residual :: Compensable a => Lens' (Compensated a) a
- uncompensated :: Compensable a => Compensated a -> a
- fadd :: Num a => a -> a -> (a -> a -> r) -> r
- add :: Num a => a -> a -> (a -> a -> r) -> r
- times :: Compensable a => a -> a -> (a -> a -> r) -> r
- squared :: Compensable a => a -> (a -> a -> r) -> r
- divide :: Compensable a => a -> a -> (a -> a -> r) -> r
- split :: Compensable a => a -> (a -> a -> r) -> r
- kahan :: (Foldable f, Compensable a) => f a -> Compensated a
- (+^) :: Compensable a => a -> Compensated a -> Compensated a
- (*^) :: Compensable a => a -> Compensated a -> Compensated a
- square :: Compensable a => Compensated a -> Compensated a

# Documentation

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

data Compensated a Source

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.

`>>>`

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

`>>>`

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

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

This extracts both the `primal`

and `residual`

components of a `Compensated`

number.

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

Used internally to construct `compensated`

values that satisfy our residual contract.

When in doubt, use

instead of `add`

a b `compensated`

`compensated`

a b

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

type Overcompensated a = Compensated (Compensated a)Source

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

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

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

Extract the `primal`

component of a `compensated`

value, when and if compensation
is no longer required.

# lifting scalars

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

computes `add`

a b k`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) -> rSource

computes `times`

a b k`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) -> rSource

computes `squared`

a k`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) -> rSource

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

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 aSource

Perform Kahan summation over a list.

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

Calculate a scalar + compensated sum with Kahan summation.

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

Compute `a * `

`Compensated`

a

# compensated operators

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

Calculate a fast square of a compensated number.