Copyright | (c) Edward Kmett 2013-2015 |
---|---|
License | BSD3 |
Maintainer | Edward Kmett <ekmett@gmail.com> |
Stability | experimental |
Portability | non-portable |
Safe Haskell | Trustworthy |
Language | Haskell98 |
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
Synopsis
- class (RealFrac 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, Floating a) => Compensable a where Source #
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.
>>>
round (Prelude.product [2..100] :: Compensated (Compensated (Compensated Double)))
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
>>>
Prelude.product [2..100]
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
with :: Compensable a => Compensated a -> (a -> a -> r) -> r Source #
This extracts both the primal
and residual
components of a Compensated
number.
compensated :: Compensable a => a -> a -> Compensated a Source #
Used internally to construct compensated
values that satisfy our residual contract.
When in doubt, use
instead of add
a b compensated
compensated
a b
Instances
Compensable Double Source # | |
Defined in Numeric.Compensated data Compensated Double :: Type Source # with :: Compensable Double => Compensated Double -> (Double -> Double -> r) -> r Source # compensated :: Double -> Double -> Compensated Double Source # | |
Compensable Float Source # | |
Defined in Numeric.Compensated data Compensated Float :: Type Source # with :: Compensable Float => Compensated Float -> (Float -> Float -> r) -> r Source # compensated :: Float -> Float -> Compensated Float Source # | |
Compensable a => Compensable (Compensated a) Source # | |
Defined in Numeric.Compensated data Compensated (Compensated a) :: Type Source # with :: Compensable (Compensated a) => Compensated (Compensated a) -> (Compensated a -> Compensated a -> r) -> r Source # compensated :: Compensated a -> Compensated a -> Compensated (Compensated a) Source # magic :: Compensated a Source # |
_Compensated :: Compensable a => Iso' (Compensated a) (a, a) Source #
type Overcompensated a = Compensated (Compensated a) Source #
primal :: Compensable a => Lens' (Compensated a) a Source #
residual :: Compensable a => Lens' (Compensated a) a Source #
uncompensated :: Compensable a => Compensated a -> a Source #
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) -> r Source #
computes add
a b kk 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 Source #
computes times
a b kk 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 Source #
computes squared
a kk 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 Source #
split :: Compensable a => a -> (a -> a -> r) -> r Source #
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 Source #
Perform Kahan summation over a list.
(+^) :: Compensable a => a -> Compensated a -> Compensated a Source #
Calculate a scalar + compensated sum with Kahan summation.
(*^) :: Compensable a => a -> Compensated a -> Compensated a Source #
Compute a *
Compensated
a
compensated operators
square :: Compensable a => Compensated a -> Compensated a Source #
Calculate a fast square of a compensated number.