{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms       #-}
{-# LANGUAGE RebindableSyntax      #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}
-- |
-- Module      : Data.Array.Accelerate.Numeric.Sum
-- Copyright   : [2017..2020] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--
-- Functions for summing floating point numbers more accurately than the
-- straightforward 'Data.Array.Accelerate.sum' operation.
--
-- In the worst case, the 'Data.Array.Accelerate.sum' function accumulates error
-- at a rate proportional to the number of values being summed. The algorithms
-- in this module implement different methods of /compensated summation/, which
-- reduce the accumulation of numeric error so that it grows much more slowly
-- than the number of inputs (e.g. logarithmically), or remains constant.
--

-- TLM: The standard formulation of the algorithms implemented here are not
-- associative; e.g. they would have a type (KBN a -> a -> KBN a). I've
-- done what seems like the sensible conversion, but somebody versed in numeric
-- analysis should probably look...
--
-- See also: <https://hackage.haskell.org/package/math-functions>
--

module Data.Array.Accelerate.Numeric.Sum (

  -- * Summation type class
  Summation(..),
  sum,

  -- * Kahan-Babuška-Neumaier summation
  KBN(..), pattern KBN_,
  kbn,

  -- * Order-2 Kahan-Babuška summation
  KB2(..), pattern KB2_,
  kb2,

  -- * Kahan summation
  Kahan(..), pattern Kahan_,
  kahan,

) where

import Data.Array.Accelerate                                        as A hiding ( sum, fromInteger )
import Data.Array.Accelerate.Type                                   as A
import Data.Array.Accelerate.Numeric.Sum.Arithmetic                 as A

import Data.Proxy
import Prelude                                                      ( fromInteger )


-- | Sum an array using a particular compensation scheme.
--
-- >>> let xs = [1.0, 1.0e100, 1.0, -1.0e100] :: [Double]
-- >>> Prelude.sum xs
-- 0.0
--
-- >>> let ys = fromList (Z:.4) [1.0, 1.0e100, 1.0, -1.0e100] :: Vector Double
-- >>> sum kbn (use ys)
-- Scalar Z [2.0]
--
sum :: (Summation s a, Shape sh) => Proxy s -> Acc (Array (sh:.Int) a) -> Acc (Array sh a)
sum :: Proxy s -> Acc (Array (sh :. Int) a) -> Acc (Array sh a)
sum Proxy s
p = (Exp (s a) -> Exp a) -> Acc (Array sh (s a)) -> Acc (Array sh a)
forall sh a b.
(Shape sh, Elt a, Elt b) =>
(Exp a -> Exp b) -> Acc (Array sh a) -> Acc (Array sh b)
A.map (Proxy s -> Exp (s a) -> Exp a
forall (s :: * -> *) a.
Summation s a =>
Proxy s -> Exp (s a) -> Exp a
from Proxy s
p)
      (Acc (Array sh (s a)) -> Acc (Array sh a))
-> (Acc (Array (sh :. Int) a) -> Acc (Array sh (s a)))
-> Acc (Array (sh :. Int) a)
-> Acc (Array sh a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Exp (s a) -> Exp (s a) -> Exp (s a))
-> Exp (s a)
-> Acc (Array (sh :. Int) (s a))
-> Acc (Array sh (s a))
forall sh a.
(Shape sh, Elt a) =>
(Exp a -> Exp a -> Exp a)
-> Exp a -> Acc (Array (sh :. Int) a) -> Acc (Array sh a)
A.fold Exp (s a) -> Exp (s a) -> Exp (s a)
forall (s :: * -> *) a.
Summation s a =>
Exp (s a) -> Exp (s a) -> Exp (s a)
add Exp (s a)
forall (s :: * -> *) a. Summation s a => Exp (s a)
zero
      (Acc (Array (sh :. Int) (s a)) -> Acc (Array sh (s a)))
-> (Acc (Array (sh :. Int) a) -> Acc (Array (sh :. Int) (s a)))
-> Acc (Array (sh :. Int) a)
-> Acc (Array sh (s a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Exp a -> Exp (s a))
-> Acc (Array (sh :. Int) a) -> Acc (Array (sh :. Int) (s a))
forall sh a b.
(Shape sh, Elt a, Elt b) =>
(Exp a -> Exp b) -> Acc (Array sh a) -> Acc (Array sh b)
A.map (Proxy s -> Exp a -> Exp (s a)
forall (s :: * -> *) a.
Summation s a =>
Proxy s -> Exp a -> Exp (s a)
into Proxy s
p)


-- | A class for the summation of floating-point numbers
--
class (Elt a, Elt (s a)) => Summation s a where
  -- | Add a value to the sum
  add  :: Exp (s a) -> Exp (s a) -> Exp (s a)

  -- | The identity of the summation
  zero :: Exp (s a)

  -- | Insert a value into the summation
  into :: Proxy s -> Exp a -> Exp (s a)

  -- | Summarise the result of summation
  from :: Proxy s -> Exp (s a) -> Exp a


-- | Kahan-Babuška-Neumaier summation. This is a little more computationally
-- costly than plain Kahan summation, but is /always/ at least as accurate.
--
data KBN a = KBN a a
  deriving (Int -> KBN a -> ShowS
[KBN a] -> ShowS
KBN a -> String
(Int -> KBN a -> ShowS)
-> (KBN a -> String) -> ([KBN a] -> ShowS) -> Show (KBN a)
forall a. Show a => Int -> KBN a -> ShowS
forall a. Show a => [KBN a] -> ShowS
forall a. Show a => KBN a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [KBN a] -> ShowS
$cshowList :: forall a. Show a => [KBN a] -> ShowS
show :: KBN a -> String
$cshow :: forall a. Show a => KBN a -> String
showsPrec :: Int -> KBN a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> KBN a -> ShowS
Show, (forall x. KBN a -> Rep (KBN a) x)
-> (forall x. Rep (KBN a) x -> KBN a) -> Generic (KBN a)
forall x. Rep (KBN a) x -> KBN a
forall x. KBN a -> Rep (KBN a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (KBN a) x -> KBN a
forall a x. KBN a -> Rep (KBN a) x
$cto :: forall a x. Rep (KBN a) x -> KBN a
$cfrom :: forall a x. KBN a -> Rep (KBN a) x
Generic)

pattern KBN_ :: Elt a => Exp a -> Exp a -> Exp (KBN a)
pattern $bKBN_ :: Exp a -> Exp a -> Exp (KBN a)
$mKBN_ :: forall r a.
Elt a =>
Exp (KBN a) -> (Exp a -> Exp a -> r) -> (Void# -> r) -> r
KBN_ s c = Pattern (c, s)
{-# COMPLETE KBN_ #-}

instance Elt a => Elt (KBN a)

-- | Return the result of a Kahan-Babuška-Neumaier sum.
--
kbn :: Proxy KBN
kbn :: Proxy KBN
kbn = Proxy KBN
forall k (t :: k). Proxy t
Proxy

kbnAdd :: (Num a, Ord a, IsFloating a) => Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a)
kbnAdd :: Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a)
kbnAdd (KBN_ Exp a
s1 Exp a
c1) (KBN_ Exp a
s2 Exp a
c2) = Exp a -> Exp a -> Exp (KBN a)
forall a. Elt a => Exp a -> Exp a -> Exp (KBN a)
KBN_ Exp a
s' Exp a
c'
  where
    s' :: Exp a
s' = Exp a
s1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s2
    c' :: Exp a
c' = Exp a
c1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
c2 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` if Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
s1 Exp a -> Exp a -> Exp Bool
forall a. Ord a => Exp a -> Exp a -> Exp Bool
>= Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
s2
                               then (Exp a
s1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
s') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s2
                               else (Exp a
s2 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
s') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s1

-- instance (Num a, Ord a) => Summation KBN a where
--   zero      = lift $ KBN (0::Exp a) (0::Exp a)
--   add       = kbnAdd
--   into _ x  = lift (KBN x 0)
--   from _ x  = let KBN s c = unlift x in s + c

instance Summation KBN Float where
  zero :: Exp (KBN Float)
zero      = KBN Float -> Exp (KBN Float)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Float -> Float -> KBN Float
forall a. a -> a -> KBN a
KBN Float
0 Float
0)
  add :: Exp (KBN Float) -> Exp (KBN Float) -> Exp (KBN Float)
add       = Exp (KBN Float) -> Exp (KBN Float) -> Exp (KBN Float)
forall a.
(Num a, Ord a, IsFloating a) =>
Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a)
kbnAdd
  into :: Proxy KBN -> Exp Float -> Exp (KBN Float)
into Proxy KBN
_ Exp Float
x  = KBN (Exp Float) -> Exp (Plain (KBN (Exp Float)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Float -> Exp Float -> KBN (Exp Float)
forall a. a -> a -> KBN a
KBN Exp Float
x Exp Float
0)
  from :: Proxy KBN -> Exp (KBN Float) -> Exp Float
from Proxy KBN
_ Exp (KBN Float)
x  = let KBN Exp Float
s Exp Float
c = Exp (Plain (KBN (Exp Float))) -> KBN (Exp Float)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (KBN (Exp Float)))
Exp (KBN Float)
x in Exp Float
s Exp Float -> Exp Float -> Exp Float
forall a. Num a => a -> a -> a
+ Exp Float
c

instance Summation KBN Double where
  zero :: Exp (KBN Double)
zero      = KBN Double -> Exp (KBN Double)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Double -> Double -> KBN Double
forall a. a -> a -> KBN a
KBN Double
0 Double
0)
  add :: Exp (KBN Double) -> Exp (KBN Double) -> Exp (KBN Double)
add       = Exp (KBN Double) -> Exp (KBN Double) -> Exp (KBN Double)
forall a.
(Num a, Ord a, IsFloating a) =>
Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a)
kbnAdd
  into :: Proxy KBN -> Exp Double -> Exp (KBN Double)
into Proxy KBN
_ Exp Double
x  = KBN (Exp Double) -> Exp (Plain (KBN (Exp Double)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Double -> Exp Double -> KBN (Exp Double)
forall a. a -> a -> KBN a
KBN Exp Double
x Exp Double
0)
  from :: Proxy KBN -> Exp (KBN Double) -> Exp Double
from Proxy KBN
_ Exp (KBN Double)
x  = let KBN Exp Double
s Exp Double
c = Exp (Plain (KBN (Exp Double))) -> KBN (Exp Double)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (KBN (Exp Double)))
Exp (KBN Double)
x in Exp Double
s Exp Double -> Exp Double -> Exp Double
forall a. Num a => a -> a -> a
+ Exp Double
c

instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KBN a) where
  type Plain (KBN a) = KBN (Plain a)
  lift :: KBN a -> Exp (Plain (KBN a))
lift (KBN a
a a
b)     = Exp (Plain a) -> Exp (Plain a) -> Exp (KBN (Plain a))
forall a. Elt a => Exp a -> Exp a -> Exp (KBN a)
KBN_ (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
a) (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
b)

instance Elt a => Unlift Exp (KBN (Exp a)) where
  unlift :: Exp (Plain (KBN (Exp a))) -> KBN (Exp a)
unlift (KBN_ a b) = Exp a -> Exp a -> KBN (Exp a)
forall a. a -> a -> KBN a
KBN Exp a
a Exp a
b


-- | Second-order Kahan-Babuška summation.  This is more computationally costly
-- than Kahan-Babuška-Neumaier summation. Its advantage is that it can lose less
-- precision (in admittedly obscure cases).
--
-- This method compensates for error in both the sum and the first-order
-- compensation term, hence the use of \"second order\" in the name.
--
data KB2 a = KB2 a a a
  deriving (Int -> KB2 a -> ShowS
[KB2 a] -> ShowS
KB2 a -> String
(Int -> KB2 a -> ShowS)
-> (KB2 a -> String) -> ([KB2 a] -> ShowS) -> Show (KB2 a)
forall a. Show a => Int -> KB2 a -> ShowS
forall a. Show a => [KB2 a] -> ShowS
forall a. Show a => KB2 a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [KB2 a] -> ShowS
$cshowList :: forall a. Show a => [KB2 a] -> ShowS
show :: KB2 a -> String
$cshow :: forall a. Show a => KB2 a -> String
showsPrec :: Int -> KB2 a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> KB2 a -> ShowS
Show, (forall x. KB2 a -> Rep (KB2 a) x)
-> (forall x. Rep (KB2 a) x -> KB2 a) -> Generic (KB2 a)
forall x. Rep (KB2 a) x -> KB2 a
forall x. KB2 a -> Rep (KB2 a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (KB2 a) x -> KB2 a
forall a x. KB2 a -> Rep (KB2 a) x
$cto :: forall a x. Rep (KB2 a) x -> KB2 a
$cfrom :: forall a x. KB2 a -> Rep (KB2 a) x
Generic)

pattern KB2_ :: Elt a => Exp a -> Exp a -> Exp a -> Exp (KB2 a)
pattern $bKB2_ :: Exp a -> Exp a -> Exp a -> Exp (KB2 a)
$mKB2_ :: forall r a.
Elt a =>
Exp (KB2 a) -> (Exp a -> Exp a -> Exp a -> r) -> (Void# -> r) -> r
KB2_ s c cc = Pattern (s, c, cc)
{-# COMPLETE KB2_ #-}

instance Elt a => Elt (KB2 a)

-- | Return the result of a second-order Kahan-Babuška sum.
--
kb2 :: Proxy KB2
kb2 :: Proxy KB2
kb2 = Proxy KB2
forall k (t :: k). Proxy t
Proxy

kb2Add :: (Num a, Ord a, IsFloating a) => Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a)
kb2Add :: Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a)
kb2Add (Exp (KB2 a) -> KB2 (Exp a)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift -> KB2 Exp a
s1 Exp a
c1 Exp a
cc1) (Exp (KB2 a) -> KB2 (Exp a)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift -> KB2 Exp a
s2 Exp a
c2 Exp a
cc2) = KB2 (Exp a) -> Exp (Plain (KB2 (Exp a)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp a -> Exp a -> Exp a -> KB2 (Exp a)
forall a. a -> a -> a -> KB2 a
KB2 Exp a
sum' Exp a
c' Exp a
cc')
  where
    sum' :: Exp a
sum'  = Exp a
s1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s2
    c' :: Exp a
c'    = Exp a
t  Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
k
    cc' :: Exp a
cc'   = Exp a
cc1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
cc2 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` if Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
t Exp a -> Exp a -> Exp Bool
forall a. Ord a => Exp a -> Exp a -> Exp Bool
>= Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
k
                                    then (Exp a
t Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
c') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
k
                                    else (Exp a
k Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
c') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
t
    t :: Exp a
t     = Exp a
c1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
c2
    k :: Exp a
k     = if Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
s1 Exp a -> Exp a -> Exp Bool
forall a. Ord a => Exp a -> Exp a -> Exp Bool
>= Exp a -> Exp a
forall a. Num a => a -> a
abs Exp a
s2
              then (Exp a
s1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
sum') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s2
              else (Exp a
s2 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
sum') Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
s1

-- instance (Num a, Ord a) => Summation KB2 a where
--   zero      = lift $ KB2 (0::Exp a) (0::Exp a) (0::Exp a)
--   add       = kb2Add
--   into _ x  = lift (KB2 x 0 0)
--   from _ x  = let KB2 s c cc = unlift x in s + c + cc

instance Summation KB2 Float where
  zero :: Exp (KB2 Float)
zero      = KB2 Float -> Exp (KB2 Float)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Float -> Float -> Float -> KB2 Float
forall a. a -> a -> a -> KB2 a
KB2 Float
0 Float
0 Float
0)
  add :: Exp (KB2 Float) -> Exp (KB2 Float) -> Exp (KB2 Float)
add       = Exp (KB2 Float) -> Exp (KB2 Float) -> Exp (KB2 Float)
forall a.
(Num a, Ord a, IsFloating a) =>
Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a)
kb2Add
  into :: Proxy KB2 -> Exp Float -> Exp (KB2 Float)
into Proxy KB2
_ Exp Float
x  = KB2 (Exp Float) -> Exp (Plain (KB2 (Exp Float)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Float -> Exp Float -> Exp Float -> KB2 (Exp Float)
forall a. a -> a -> a -> KB2 a
KB2 Exp Float
x Exp Float
0 Exp Float
0)
  from :: Proxy KB2 -> Exp (KB2 Float) -> Exp Float
from Proxy KB2
_ Exp (KB2 Float)
x  = let KB2 Exp Float
s Exp Float
c Exp Float
cc = Exp (Plain (KB2 (Exp Float))) -> KB2 (Exp Float)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (KB2 (Exp Float)))
Exp (KB2 Float)
x in Exp Float
s Exp Float -> Exp Float -> Exp Float
forall a. Num a => a -> a -> a
+ Exp Float
c Exp Float -> Exp Float -> Exp Float
forall a. Num a => a -> a -> a
+ Exp Float
cc

instance Summation KB2 Double where
  zero :: Exp (KB2 Double)
zero      = KB2 Double -> Exp (KB2 Double)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Double -> Double -> Double -> KB2 Double
forall a. a -> a -> a -> KB2 a
KB2 Double
0 Double
0 Double
0)
  add :: Exp (KB2 Double) -> Exp (KB2 Double) -> Exp (KB2 Double)
add       = Exp (KB2 Double) -> Exp (KB2 Double) -> Exp (KB2 Double)
forall a.
(Num a, Ord a, IsFloating a) =>
Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a)
kb2Add
  into :: Proxy KB2 -> Exp Double -> Exp (KB2 Double)
into Proxy KB2
_ Exp Double
x  = KB2 (Exp Double) -> Exp (Plain (KB2 (Exp Double)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Double -> Exp Double -> Exp Double -> KB2 (Exp Double)
forall a. a -> a -> a -> KB2 a
KB2 Exp Double
x Exp Double
0 Exp Double
0)
  from :: Proxy KB2 -> Exp (KB2 Double) -> Exp Double
from Proxy KB2
_ Exp (KB2 Double)
x  = let KB2 Exp Double
s Exp Double
c Exp Double
cc = Exp (Plain (KB2 (Exp Double))) -> KB2 (Exp Double)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (KB2 (Exp Double)))
Exp (KB2 Double)
x in Exp Double
s Exp Double -> Exp Double -> Exp Double
forall a. Num a => a -> a -> a
+ Exp Double
c Exp Double -> Exp Double -> Exp Double
forall a. Num a => a -> a -> a
+ Exp Double
cc

instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KB2 a) where
  type Plain (KB2 a) = KB2 (Plain a)
  lift :: KB2 a -> Exp (Plain (KB2 a))
lift (KB2 a
a a
b a
c)   = Exp (Plain a)
-> Exp (Plain a) -> Exp (Plain a) -> Exp (KB2 (Plain a))
forall a. Elt a => Exp a -> Exp a -> Exp a -> Exp (KB2 a)
KB2_ (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
a) (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
b) (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
c)

instance Elt a => Unlift Exp (KB2 (Exp a)) where
  unlift :: Exp (Plain (KB2 (Exp a))) -> KB2 (Exp a)
unlift (KB2_ a b c) = Exp a -> Exp a -> Exp a -> KB2 (Exp a)
forall a. a -> a -> a -> KB2 a
KB2 Exp a
a Exp a
b Exp a
c


-- | Kahan summation. This is the least accurate of the compensated summation
-- methods. This summation method is included only for completeness.
--
data Kahan a = Kahan a a
  deriving (Int -> Kahan a -> ShowS
[Kahan a] -> ShowS
Kahan a -> String
(Int -> Kahan a -> ShowS)
-> (Kahan a -> String) -> ([Kahan a] -> ShowS) -> Show (Kahan a)
forall a. Show a => Int -> Kahan a -> ShowS
forall a. Show a => [Kahan a] -> ShowS
forall a. Show a => Kahan a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Kahan a] -> ShowS
$cshowList :: forall a. Show a => [Kahan a] -> ShowS
show :: Kahan a -> String
$cshow :: forall a. Show a => Kahan a -> String
showsPrec :: Int -> Kahan a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Kahan a -> ShowS
Show, (forall x. Kahan a -> Rep (Kahan a) x)
-> (forall x. Rep (Kahan a) x -> Kahan a) -> Generic (Kahan a)
forall x. Rep (Kahan a) x -> Kahan a
forall x. Kahan a -> Rep (Kahan a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (Kahan a) x -> Kahan a
forall a x. Kahan a -> Rep (Kahan a) x
$cto :: forall a x. Rep (Kahan a) x -> Kahan a
$cfrom :: forall a x. Kahan a -> Rep (Kahan a) x
Generic)

pattern Kahan_ :: Elt a => Exp a -> Exp a -> Exp (Kahan a)
pattern $bKahan_ :: Exp a -> Exp a -> Exp (Kahan a)
$mKahan_ :: forall r a.
Elt a =>
Exp (Kahan a) -> (Exp a -> Exp a -> r) -> (Void# -> r) -> r
Kahan_ s c = Pattern (s, c)
{-# COMPLETE Kahan_ #-}

instance Elt a => Elt (Kahan a)

-- | Return the result of a Kahan sum.
--
kahan :: Proxy Kahan
kahan :: Proxy Kahan
kahan = Proxy Kahan
forall k (t :: k). Proxy t
Proxy

kahanAdd :: (Num a, IsFloating a) => Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a)
kahanAdd :: Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a)
kahanAdd (Kahan_ Exp a
s1 Exp a
c1) (Kahan_ Exp a
s2 Exp a
c2) = Exp a -> Exp a -> Exp (Kahan a)
forall a. Elt a => Exp a -> Exp a -> Exp (Kahan a)
Kahan_ Exp a
s' Exp a
c'
  where
    s' :: Exp a
s'  = Exp a
s1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fadd` Exp a
y
    c' :: Exp a
c'  = (Exp a
s' Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
s1) Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
y
    y :: Exp a
y   = Exp a
s2 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
c1 Exp a -> Exp a -> Exp a
forall a. (Num a, IsFloating a) => Exp a -> Exp a -> Exp a
`fsub` Exp a
c2

-- instance (Num a, Ord a) => Summation Kahan a where
--   zero      = lift $ Kahan (0::Exp a) (0::Exp a)
--   add       = kahanAdd
--   into _ x  = lift (Kahan x 0)
--   from _ x  = let Kahan s _ = unlift x in s

instance Summation Kahan Float where
  zero :: Exp (Kahan Float)
zero      = Kahan Float -> Exp (Kahan Float)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Float -> Float -> Kahan Float
forall a. a -> a -> Kahan a
Kahan Float
0 Float
0)
  add :: Exp (Kahan Float) -> Exp (Kahan Float) -> Exp (Kahan Float)
add       = Exp (Kahan Float) -> Exp (Kahan Float) -> Exp (Kahan Float)
forall a.
(Num a, IsFloating a) =>
Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a)
kahanAdd
  into :: Proxy Kahan -> Exp Float -> Exp (Kahan Float)
into Proxy Kahan
_ Exp Float
x  = Kahan (Exp Float) -> Exp (Plain (Kahan (Exp Float)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Float -> Exp Float -> Kahan (Exp Float)
forall a. a -> a -> Kahan a
Kahan Exp Float
x Exp Float
0)
  from :: Proxy Kahan -> Exp (Kahan Float) -> Exp Float
from Proxy Kahan
_ Exp (Kahan Float)
x  = let Kahan Exp Float
s Exp Float
_ = Exp (Plain (Kahan (Exp Float))) -> Kahan (Exp Float)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (Kahan (Exp Float)))
Exp (Kahan Float)
x in Exp Float
s

instance Summation Kahan Double where
  zero :: Exp (Kahan Double)
zero      = Kahan Double -> Exp (Kahan Double)
forall e. (HasCallStack, Elt e) => e -> Exp e
constant (Double -> Double -> Kahan Double
forall a. a -> a -> Kahan a
Kahan Double
0 Double
0)
  add :: Exp (Kahan Double) -> Exp (Kahan Double) -> Exp (Kahan Double)
add       = Exp (Kahan Double) -> Exp (Kahan Double) -> Exp (Kahan Double)
forall a.
(Num a, IsFloating a) =>
Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a)
kahanAdd
  into :: Proxy Kahan -> Exp Double -> Exp (Kahan Double)
into Proxy Kahan
_ Exp Double
x  = Kahan (Exp Double) -> Exp (Plain (Kahan (Exp Double)))
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift (Exp Double -> Exp Double -> Kahan (Exp Double)
forall a. a -> a -> Kahan a
Kahan Exp Double
x Exp Double
0)
  from :: Proxy Kahan -> Exp (Kahan Double) -> Exp Double
from Proxy Kahan
_ Exp (Kahan Double)
x  = let Kahan Exp Double
s Exp Double
_ = Exp (Plain (Kahan (Exp Double))) -> Kahan (Exp Double)
forall (c :: * -> *) e. Unlift c e => c (Plain e) -> e
unlift Exp (Plain (Kahan (Exp Double)))
Exp (Kahan Double)
x in Exp Double
s

instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Kahan a) where
  type Plain (Kahan a) = Kahan (Plain a)
  lift :: Kahan a -> Exp (Plain (Kahan a))
lift (Kahan a
a a
b)     = Exp (Plain a) -> Exp (Plain a) -> Exp (Kahan (Plain a))
forall a. Elt a => Exp a -> Exp a -> Exp (Kahan a)
Kahan_ (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
a) (a -> Exp (Plain a)
forall (c :: * -> *) e. Lift c e => e -> c (Plain e)
lift a
b)

instance Elt a => Unlift Exp (Kahan (Exp a)) where
  unlift :: Exp (Plain (Kahan (Exp a))) -> Kahan (Exp a)
unlift (Kahan_ a b) = Exp a -> Exp a -> Kahan (Exp a)
forall a. a -> a -> Kahan a
Kahan Exp a
a Exp a
b