{-# 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 Show
instance NFData KBN where
rnf KBN {} = ()
getKBN :: KBN -> Double
getKBN (KBN x e) = x + e
zeroKBN :: KBN
zeroKBN = KBN 0 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 = sumKBNWith id
-- | Generalized version of 'sumKBN'.
sumKBNWith :: F.Foldable f => (a -> Double) -> f a -> Double
sumKBNWith f xs = getKBN (F.foldl' (\k a -> addKBN k (f a)) (KBN 0 0) xs)
-- | Add a 'Double' to 'KBN' accumulator.
addKBN :: KBN -> Double -> KBN
addKBN (KBN acc c) x = KBN acc' c' where
acc' = acc + x
c' | abs acc >= abs x = c + ((acc - acc') + x)
| otherwise = c + ((x - acc') + acc)