#ifndef MIN_VERSION_vector
#define MIN_VERSION_vector(x,y,z) 1
#endif
module Numeric.Compensated
( Compensable(..)
, _Compensated
, Overcompensated
, primal
, residual
, uncompensated
, fadd
, add, times, squared, divide, split
, kahan, (+^), (*^)
, square
) where
#if __GLASGOW_HASKELL__ < 710
import Control.Applicative
#endif
import Control.Lens as L
import Control.DeepSeq
import Control.Monad
import Data.Binary as Binary
import Data.Data
import Data.Foldable as Foldable
import Data.Function (on)
import Data.Hashable
import Data.Ratio
import Data.SafeCopy
import Data.Serialize as Serialize
import Data.Bytes.Serial as Bytes
import Data.Semigroup
import Data.Vector.Unboxed as U
import Data.Vector.Generic as G
import Data.Vector.Generic.Mutable as M
import Foreign.Ptr
import Foreign.Storable
import Numeric.Log
import Text.Read as T
import Text.Show as T
add :: Num a => a -> a -> (a -> a -> r) -> r
add a b k = k x y where
x = a + b
z = x a
y = (a (x z)) + (b z)
fadd :: Num a => a -> a -> (a -> a -> r) -> r
fadd a b k = k x (b (x a)) where
x = a + b
times :: Compensable a => a -> a -> (a -> a -> r) -> r
times a b k =
split a $ \a1 a2 ->
split b $ \b1 b2 ->
let x = a * b in k x (a2*b2 (((x a1*b1) a2*b1) a1*b2))
divide :: Compensable a => a -> a -> (a -> a -> r) -> r
divide a b = with (aX * ms) where
x0 = recip b
aX = times a x0 compensated
m = 1 +^ negate (times b x0 compensated)
mm = m*m
ms = 1+((m+mm)+m*mm)
renorm :: Compensable a => a -> a -> a -> Compensated a
renorm a b c =
fadd b c $ \x1 y1 ->
fadd a x1 $ \x2 y2 ->
fadd x2 (y1 + y2) compensated
squared :: Compensable a => a -> (a -> a -> r) -> r
squared a k =
split a $ \a1 a2 ->
let x = a * a in k x (a2*a2 ((x a1*a1) 2*(a2*a1)))
square :: Compensable a => Compensated a -> Compensated a
square m =
with m $ \a b ->
squared a $ \x1 y1 ->
times a b $ \x2 y2 ->
add y1 (x2*2) $ \x3 y3 ->
renorm x1 x3 (b*b + 2*y2 + y3)
split :: Compensable a => a -> (a -> a -> r) -> r
split a k = k x y where
c = magic*a
x = c (c a)
y = a x
(+^) :: Compensable a => a -> Compensated a -> Compensated a
a +^ m = with m $ \b c -> let y = a c; t = b + y in compensated t ((t b) y)
(*^) :: Compensable a => a -> Compensated a -> Compensated a
c *^ m =
with m $ \ a b ->
times c a $ \x1 y1 ->
times c b $ \x2 y2 ->
fadd x1 x2 $ \x3 y3 ->
add y1 y3 $ \x4 y4 ->
renorm x3 x4 (y4 + y2)
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
instance Compensable Double where
data Compensated Double = CD !Double !Double
with (CD a b) k = k a b
compensated = CD
magic = 134217729
instance Compensable Float where
data Compensated Float = CF !Float !Float
with (CF a b) k = k a b
compensated = CF
magic = 4097
instance Compensable a => Compensable (Compensated a) where
data Compensated (Compensated a) = CC !(Compensated a) !(Compensated a)
with (CC a b) k = k a b
compensated = CC
magic = times (magic 1) (magic 1) $ \ x y -> compensated x (y + 1)
#if __GLASGOW_HASKELL__ < 707
instance Typeable1 Compensated where
typeOf1 _ = mkTyConApp (mkTyCon3 "analytics" "Data.Analytics.Numeric.Compensated" "Compensated") []
#else
deriving instance Typeable Compensated
#endif
instance (Compensable a, Hashable a) => Hashable (Compensated a) where
hashWithSalt n m = with m $ \a b -> hashWithSalt n (a,b)
instance (Compensable a, Data a) => Data (Compensated a) where
gfoldl f z m = with m $ \a b -> z compensated `f` a `f` b
toConstr _ = compensatedConstr
gunfold k z c = case constrIndex c of
1 -> k (k (z compensated))
_ -> error "gunfold"
dataTypeOf _ = compensatedDataType
dataCast1 f = gcast1 f
compensatedConstr :: Constr
compensatedConstr = mkConstr compensatedDataType "compensated" [] Prefix
compensatedDataType :: DataType
compensatedDataType = mkDataType "Data.Analytics.Numeric.Compensated" [compensatedConstr]
instance (Compensable a, NFData a) => NFData (Compensated a) where
rnf m = with m $ \x y -> rnf x `seq` rnf y `seq` ()
instance (Compensable a, Show a) => Show (Compensated a) where
showsPrec d m = with m $ \a b -> showParen (d > 10) $
showString "compensated " . T.showsPrec 11 a . showChar ' ' . T.showsPrec 11 b
instance (Compensable a, Read a) => Read (Compensated a) where
readPrec = parens $ prec 10 $ do
Ident "compensated" <- lexP
a <- step T.readPrec
b <- step T.readPrec
return $ compensated a b
type Overcompensated a = Compensated (Compensated a)
_Compensated :: Compensable a => Iso' (Compensated a) (a, a)
_Compensated = iso (`with` (,)) (uncurry compensated)
primal :: Compensable a => Lens' (Compensated a) a
primal f c = with c $ \a b -> f a <&> \a' -> compensated a' b
residual :: Compensable a => Lens' (Compensated a) a
residual f c = with c $ \a b -> compensated a <$> f b
uncompensated :: Compensable a => Compensated a -> a
uncompensated c = with c const
type instance Index (Compensated a) = Int
instance (Compensable a, Compensable b) => Each (Compensated a) (Compensated b) a b where
each f m = with m $ \a b -> compensated <$> f a <*> f b
instance Compensable a => Eq (Compensated a) where
m == n = with m $ \a b -> with n $ \c d -> a == c && b == d
m /= n = with m $ \a b -> with n $ \c d -> a /= c || b /= d
instance Compensable a => Ord (Compensated a) where
compare m n = with m $ \a b -> with n $ \c d -> compare a c <> compare b d
m <= n = with m $ \a b -> with n $ \c d -> case compare a c of
LT -> True
EQ -> b <= d
GT -> False
m >= n = with m $ \a b -> with n $ \c d -> case compare a c of
LT -> False
EQ -> b >= d
GT -> a >= c
m > n = with m $ \a b -> with n $ \c d -> case compare a c of
LT -> False
EQ -> b > d
GT -> a > c
m < n = with m $ \a b -> with n $ \c d -> case compare a c of
LT -> True
EQ -> b < d
GT -> False
instance Compensable a => Semigroup (Compensated a) where
(<>) = (+)
instance Compensable a => Monoid (Compensated a) where
mempty = compensated 0 0
mappend = (+)
kahan :: (Foldable f, Compensable a) => f a -> Compensated a
kahan = Foldable.foldr (+^) mempty
instance Compensable a => Num (Compensated a) where
m + n =
with m $ \a b ->
with n $ \c d ->
add a c $ \x1 y1 ->
add y1 d $ \x2 y2 ->
add b x2 $ \x3 y3 ->
add x1 x3 $ \x4 y4 ->
add x4 (y2 + y3 + y4) compensated
m * n =
with m $ \a b ->
with n $ \c d ->
times a c $ \x1 y1 ->
times b c $ \x2 y2 ->
times a d $ \x3 y3 ->
add x1 x2 $ \x4 y4 ->
add x3 x4 $ \x5 y5 ->
add y1 y4 $ \x6 y6 ->
add y5 x6 $ \x7 y7 ->
add x5 x7 $ \x8 y8 ->
add x8 (b*d + y2 + y3 + y6 + y7 + y8) compensated
negate m = with m (on compensated negate)
x y = x + negate y
signum m = with m $ \a _ -> compensated (signum a) 0
abs m = with m $ \a b ->
if a < 0
then compensated (negate a) (negate b)
else compensated a b
fromInteger i = add x (fromInteger (i round x)) compensated where
x = fromInteger i
instance Compensable a => Enum (Compensated a) where
succ a = a + 1
pred a = a 1
toEnum i = add x (fromIntegral (i round x)) compensated where
x = fromIntegral i
fromEnum = round
enumFrom a = a : Prelude.enumFrom (a + 1)
enumFromThen a b = a : Prelude.enumFromThen b (b a + b)
enumFromTo a b
| a <= b = a : Prelude.enumFromTo (a + 1) b
| otherwise = []
enumFromThenTo a b c
| a <= b = up a
| otherwise = down a
where
delta = b a
up x | x <= c = x : up (x + delta)
| otherwise = []
down x | c <= x = x : down (x + delta)
| otherwise = []
instance Compensable a => Fractional (Compensated a) where
recip m = with m $ \a b -> add (recip a) (b / (a * a)) compensated
a / b = (a*x0) * (1+((m+mm)+m*mm)) where
x0 = recip b
m = 1 b*x0
mm = m*m
fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
instance Compensable a => Real (Compensated a) where
toRational m = with m (on (+) toRational)
instance Compensable a => RealFrac (Compensated a) where
properFraction m = with m $ \a b -> case properFraction a of
(w, p) -> add p b $ \ x y -> case properFraction x of
(w',q) -> (w + w', add q y compensated)
instance (Compensable a, Binary a) => Binary (Compensated a) where
get = compensated <$> Binary.get <*> Binary.get
put m = with m $ \a b -> do
Binary.put a
Binary.put b
instance (Compensable a, Serialize a) => Serialize (Compensated a) where
get = compensated <$> Serialize.get <*> Serialize.get
put m = with m $ \a b -> do
Serialize.put a
Serialize.put b
instance (Compensable a, Serial a) => Serial (Compensated a) where
deserialize = compensated <$> Bytes.deserialize <*> Bytes.deserialize
serialize m = with m $ \a b -> do
Bytes.serialize a
Bytes.serialize b
instance (Compensable a, Serialize a) => SafeCopy (Compensated a)
instance (Compensable a, Storable a) => Storable (Compensated a) where
sizeOf _ = sizeOf (undefined :: a) * 2
alignment _ = alignment (undefined :: a)
peekElemOff p o | q <- castPtr p, o2 <- o * 2 =
compensated <$> peekElemOff q o2 <*> peekElemOff q (o2+1)
pokeElemOff p o m | q <- castPtr p, o2 <- o * 2 = with m $ \a b -> do
pokeElemOff q o2 a
pokeElemOff q (o2+1) b
peekByteOff p o | q <- castPtr p =
compensated <$> peekByteOff q o <*> peekByteOff q (o + sizeOf (undefined :: a))
pokeByteOff p o m | q <- castPtr p = with m $ \a b -> do
pokeByteOff q o a
pokeByteOff q (o+sizeOf (undefined :: a)) b
peek p | q <- castPtr p = compensated <$> peek q <*> peekElemOff q 1
poke p m | q <- castPtr p = with m $ \a b -> do
poke q a
pokeElemOff q 1 b
newtype instance U.MVector s (Compensated a) = MV_Compensated (U.MVector s (a,a))
newtype instance U.Vector (Compensated a) = V_Compensated (U.Vector (a, a))
instance (Compensable a, Unbox a) => M.MVector U.MVector (Compensated a) where
basicLength (MV_Compensated v) = M.basicLength v
basicUnsafeSlice i n (MV_Compensated v) = MV_Compensated $ M.basicUnsafeSlice i n v
basicOverlaps (MV_Compensated v1) (MV_Compensated v2) = M.basicOverlaps v1 v2
basicUnsafeNew n = MV_Compensated `liftM` M.basicUnsafeNew n
basicUnsafeReplicate n m = with m $ \x y -> MV_Compensated `liftM` M.basicUnsafeReplicate n (x,y)
basicUnsafeRead (MV_Compensated v) i = uncurry compensated `liftM` M.basicUnsafeRead v i
basicUnsafeWrite (MV_Compensated v) i m = with m $ \ x y -> M.basicUnsafeWrite v i (x,y)
basicClear (MV_Compensated v) = M.basicClear v
basicSet (MV_Compensated v) m = with m $ \ x y -> M.basicSet v (x,y)
basicUnsafeCopy (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeCopy v1 v2
basicUnsafeMove (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeMove v1 v2
basicUnsafeGrow (MV_Compensated v) n = MV_Compensated `liftM` M.basicUnsafeGrow v n
#if MIN_VERSION_vector(0,11,0)
basicInitialize (MV_Compensated v) = M.basicInitialize v
#endif
instance (Compensable a, Unbox a) => G.Vector U.Vector (Compensated a) where
basicUnsafeFreeze (MV_Compensated v) = V_Compensated `liftM` G.basicUnsafeFreeze v
basicUnsafeThaw (V_Compensated v) = MV_Compensated `liftM` G.basicUnsafeThaw v
basicLength (V_Compensated v) = G.basicLength v
basicUnsafeSlice i n (V_Compensated v) = V_Compensated $ G.basicUnsafeSlice i n v
basicUnsafeIndexM (V_Compensated v) i
= uncurry compensated `liftM` G.basicUnsafeIndexM v i
basicUnsafeCopy (MV_Compensated mv) (V_Compensated v)
= G.basicUnsafeCopy mv v
elemseq _ m z = with m $ \x y -> G.elemseq (undefined :: U.Vector a) x
$ G.elemseq (undefined :: U.Vector a) y z
instance Compensable a => Floating (Compensated a) where
#ifdef SPECIALIZE_INSTANCES
#endif
exp m =
with m $ \a b ->
times (exp a) (exp b) compensated
sin m =
with m $ \a b ->
times (sin a) (cos b) $ \x1 y1 ->
times (sin b) (cos a) $ \x2 y2 ->
add x1 x2 $ \x3 y3 ->
add y1 y2 $ \x4 y4 ->
add x4 y3 $ \x5 y5 ->
add x5 x3 $ \x6 y6 ->
add (y4 + y5 + y6) x6 compensated
cos m =
with m $ \a b ->
times (cos a) (cos b) $ \x1 y1 ->
times (sin b) (sin a) $ \x2 y2 ->
add x1 x2 $ \x3 y3 ->
add y1 y2 $ \x4 y4 ->
add x4 y3 $ \x5 y5 ->
add x5 x3 $ \x6 y6 ->
add (y4 + y5 + y6) x6 compensated
tan m =
with m $ \a b ->
add (tan a) (tan b) compensated /
(1 +^ times (tan a) (tan b) compensated)
sinh m =
with m $ \a b ->
times (sinh a) (cosh b) $ \x1 y1 ->
times (cosh a) (sinh b) $ \x2 y2 ->
add x1 x2 $ \x3 y3 ->
add y1 y2 $ \x4 y4 ->
add x4 y3 $ \x5 y5 ->
add x5 x3 $ \x6 y6 ->
add (y4 + y5 + y6) x6 compensated
cosh m =
with m $ \a b ->
times (cosh a) (cosh b) $ \x1 y1 ->
times (sinh b) (sinh a) $ \x2 y2 ->
add x1 x2 $ \x3 y3 ->
add y1 y2 $ \x4 y4 ->
add x4 y3 $ \x5 y5 ->
add x5 x3 $ \x6 y6 ->
add (y4 + y5 + y6) x6 compensated
tanh m =
with m $ \a b ->
fadd (tanh a) (tanh b) compensated /
(1 +^ times (tanh a) (tanh b) compensated)
log m =
with m $ \ a b -> let
xy1 = add (log a) (b/a) compensated
xy2 = xy1 + m * exp (xy1) 1
in xy2 + m * exp (xy2) 1
sqrt m = with (z4 + m/z4) $ on compensated (/2) where
z0 = sqrt (m^.primal)
z1 = with (z0 +^ (m / compensated z0 0)) $ on compensated (/2)
z2 = with (z1 + m/z1) $ on compensated (/2)
z3 = with (z2 + m/z2) $ on compensated (/2)
z4 = with (z3 + m/z3) $ on compensated (/2)
pi = error "TODO"
asin = error "TODO"
atan = error "TODO"
acos = error "TODO"
asinh = error "TODO"
atanh = error "TODO"
acosh = error "TODO"
instance (Compensable a, Precise a) => Precise (Compensated a) where
log1p a = log (1 + a)
expm1 a = exp a 1
log1mexp a | a <= log 2 = log (negate (expm1 (negate a)))
| otherwise = log1p (negate (exp (negate a)))
log1pexp a
| a <= 18 = log1p (exp a)
| a <= 100 = a + exp (negate a)
| otherwise = a