module Numeric.Compensated
( Compensable(..)
, _Compensated
, Overcompensated
, primal
, residual
, uncompensated
, fadd
, add, times, squared, divide, split
, kahan, (+^), (*^)
, square
) where
import Control.Applicative
import Control.Lens as L
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.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)
instance Typeable1 Compensated where
typeOf1 _ = mkTyConApp (mkTyCon3 "analytics" "Data.Analytics.Numeric.Compensated" "Compensated") []
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, 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 (Applicative f, Compensable a, Compensable b) => Each f (Compensated a) (Compensated b) a b where
each f m = with m $ \a b -> compensated <$> L.indexed f (0 :: Int) a <*> L.indexed f (1 :: Int) 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 (Reviewable p, Functor f, Compensable a, a ~ b) => Cons p f (Compensated a) (Compensated b) a b where
_Cons = unto $ \(a, e) -> with e $ \b c -> let y = a c; t = b + y in compensated t ((t b) y)
instance (Reviewable p, Functor f, Compensable a, a ~ b) => Snoc p f (Compensated a) (Compensated b) a b where
_Snoc = unto $ \(e, a) -> with e $ \b c -> let y = a c; t = b + y in compensated t ((t b) y)
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, 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
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