{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
-- | https://en.wikipedia.org/wiki/KBN_summation_algorithm
--
-- @math-functions@ has KBN-Babuška-Neumaier summation algorithm as well
-- in @Numeric.Sum@ module.
module Numeric.KBN (
    KBN (..),
    zeroKBN,
    getKBN,
    sumKBN,
    addKBN,
) where

import Control.DeepSeq (NFData (..))

import qualified Data.Foldable as F

-- $setup
-- >>> import qualified Data.List

-- | KBN summation accumulator.
data KBN = KBN !Double !Double
  deriving Int -> KBN -> ShowS
[KBN] -> ShowS
KBN -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [KBN] -> ShowS
$cshowList :: [KBN] -> ShowS
show :: KBN -> String
$cshow :: KBN -> String
showsPrec :: Int -> KBN -> ShowS
$cshowsPrec :: Int -> KBN -> ShowS
Show

instance NFData KBN where
    rnf :: KBN -> ()
rnf KBN {} = ()

getKBN :: KBN -> Double
getKBN :: KBN -> Double
getKBN (KBN Double
x Double
e) = Double
x forall a. Num a => a -> a -> a
+ Double
e

zeroKBN :: KBN
zeroKBN :: KBN
zeroKBN = Double -> Double -> KBN
KBN Double
0 Double
0

-- | KBN summation algorithm.
--
-- >>> sumKBN (replicate 10 0.1)
-- 1.0
--
-- >>> Data.List.foldl' (+) 0 (replicate 10 0.1) :: Double
-- 0.9999999999999999
--
-- >>> sumKBN [1, 1e100, 1, -1e100]
-- 2.0
--
-- >>> Data.List.foldl' (+) 0 [1, 1e100, 1, -1e100]
-- 0.0
--
sumKBN :: F.Foldable f => f Double -> Double
sumKBN :: forall (f :: * -> *). Foldable f => f Double -> Double
sumKBN = forall (f :: * -> *) a.
Foldable f =>
(a -> Double) -> f a -> Double
sumKBNWith forall a. a -> a
id

-- | Generalized version of 'sumKBN'.
sumKBNWith :: F.Foldable f => (a -> Double) -> f a -> Double
sumKBNWith :: forall (f :: * -> *) a.
Foldable f =>
(a -> Double) -> f a -> Double
sumKBNWith a -> Double
f f a
xs = KBN -> Double
getKBN (forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
F.foldl' (\KBN
k a
a -> KBN -> Double -> KBN
addKBN KBN
k (a -> Double
f a
a)) (Double -> Double -> KBN
KBN Double
0 Double
0) f a
xs)

-- | Add a 'Double' to 'KBN' accumulator.
addKBN :: KBN -> Double -> KBN
addKBN :: KBN -> Double -> KBN
addKBN (KBN Double
acc Double
c) Double
x = Double -> Double -> KBN
KBN Double
acc' Double
c' where
    acc' :: Double
acc' = Double
acc forall a. Num a => a -> a -> a
+ Double
x
    c' :: Double
c' | forall a. Num a => a -> a
abs Double
acc forall a. Ord a => a -> a -> Bool
>= forall a. Num a => a -> a
abs Double
x = Double
c forall a. Num a => a -> a -> a
+ ((Double
acc forall a. Num a => a -> a -> a
- Double
acc') forall a. Num a => a -> a -> a
+ Double
x)
       | Bool
otherwise        = Double
c forall a. Num a => a -> a -> a
+ ((Double
x forall a. Num a => a -> a -> a
- Double
acc') forall a. Num a => a -> a -> a
+ Double
acc)