diff --git a/numhask.cabal b/numhask.cabal --- a/numhask.cabal +++ b/numhask.cabal @@ -1,5 +1,5 @@ name: numhask -version: 0.2.2.0 +version: 0.2.3.0 synopsis: numeric classes description: A numeric class heirarchy. category: mathematics @@ -50,6 +50,8 @@ NumHask.Algebra.Module NumHask.Algebra.Multiplicative NumHask.Algebra.Singleton + NumHask.Data + NumHask.Data.Complex + NumHask.Data.LogField other-modules: - Paths_numhask default-language: Haskell2010 diff --git a/src/NumHask/Algebra.hs b/src/NumHask/Algebra.hs --- a/src/NumHask/Algebra.hs +++ b/src/NumHask/Algebra.hs @@ -20,10 +20,10 @@ , module NumHask.Algebra.Multiplicative , module NumHask.Algebra.Rational , module NumHask.Algebra.Ring - , Complex(..) + , module NumHask.Data.Complex ) where -import Data.Complex (Complex(..)) +import NumHask.Data.Complex (Complex(..)) import NumHask.Algebra.Additive import NumHask.Algebra.Basis import NumHask.Algebra.Distribution diff --git a/src/NumHask/Algebra/Field.hs b/src/NumHask/Algebra/Field.hs --- a/src/NumHask/Algebra/Field.hs +++ b/src/NumHask/Algebra/Field.hs @@ -2,6 +2,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE DefaultSignatures #-} {-# OPTIONS_GHC -Wall #-} -- | Field classes @@ -115,11 +116,12 @@ -- > round a == floor (a + one/(one+one)) -- -- fixme: had to redefine Signed operators here because of the Field import in Metric, itself due to Complex being defined there -class (P.Ord a, Field a, P.Eq b, Integral b, AdditiveGroup b, MultiplicativeUnital b) => +class (Field a, Integral b, AdditiveGroup b, MultiplicativeUnital b) => QuotientField a b where properFraction :: a -> (b, a) round :: a -> b + default round :: (P.Ord a, P.Eq b) => a -> b round x = case properFraction x of (n,r) -> let m = bool (n+one) (n-one) (r P.< zero) @@ -134,10 +136,12 @@ P.GT -> m ceiling :: a -> b + default ceiling :: (P.Ord a) => a -> b ceiling x = bool n (n+one) (r P.> zero) where (n,r) = properFraction x floor :: a -> b + default floor :: (P.Ord a) => a -> b floor x = bool n (n-one) (r P.< zero) where (n,r) = properFraction x diff --git a/src/NumHask/Algebra/Metric.hs b/src/NumHask/Algebra/Metric.hs --- a/src/NumHask/Algebra/Metric.hs +++ b/src/NumHask/Algebra/Metric.hs @@ -1,5 +1,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE DefaultSignatures #-} {-# OPTIONS_GHC -Wall #-} -- | Metric classes @@ -306,12 +307,13 @@ distanceLp p a b = fromInteger (normLp (toInteger p) (toInteger a - toInteger b)) -- | todo: This should probably be split off into some sort of alternative Equality logic, but to what end? -class (Eq a, AdditiveGroup a) => +class (Eq a, AdditiveUnital a) => Epsilon a where nearZero :: a -> Bool nearZero a = a == zero aboutEqual :: a -> a -> Bool + default aboutEqual :: AdditiveGroup a => a -> a -> Bool aboutEqual a b = nearZero $ a - b positive :: (Signed a) => a -> Bool @@ -337,7 +339,7 @@ instance Epsilon Integer -instance (Epsilon a) => Epsilon (Complex a) where +instance (Epsilon a, AdditiveGroup a) => Epsilon (Complex a) where nearZero (rx :+ ix) = nearZero rx && nearZero ix aboutEqual a b = nearZero $ a - b diff --git a/src/NumHask/Algebra/Rational.hs b/src/NumHask/Algebra/Rational.hs --- a/src/NumHask/Algebra/Rational.hs +++ b/src/NumHask/Algebra/Rational.hs @@ -30,8 +30,21 @@ import NumHask.Algebra.Ring import NumHask.Algebra.Field -data Ratio a = !a :% !a deriving (P.Eq, P.Show) +data Ratio a = !a :% !a deriving (P.Show) +instance (P.Eq a, AdditiveUnital a) => P.Eq (Ratio a) where + a == b + | (isRNaN a P.|| isRNaN b) = P.False + | P.otherwise = (x P.== x') P.&& (y P.== y') + where + (x:%y) = a + (x':%y') = b + +isRNaN :: (P.Eq a, AdditiveUnital a) => Ratio a -> P.Bool +isRNaN (x :% y) | (x P.== zero P.&& y P.== zero) = P.True + | P.otherwise = P.False + + type Rational = Ratio Integer instance (P.Ord a, Multiplicative a, Integral a) => P.Ord (Ratio a) where @@ -39,8 +52,11 @@ (x:%y) < (x':%y') = x * y' P.< x' * y instance (P.Ord a, Integral a, Signed a, AdditiveInvertible a) => AdditiveMagma (Ratio a) where - (x:%y) `plus` (x':%y') = - reduce ((x `times` y') `plus` (x' `times` y)) (y `times` y') + (x :% y) `plus` (x' :% y') + | (y P.== zero P.&& y' P.== zero) = sign (x `plus` x') :% zero + | (y P.== zero) = x :% y + | (y' P.== zero) = x' :% y' + | P.otherwise = reduce ((x `times` y') `plus` (x' `times` y)) (y `times` y') instance (P.Ord a, Integral a, Signed a, AdditiveInvertible a) => AdditiveUnital (Ratio a) where zero = zero :% one diff --git a/src/NumHask/Data.hs b/src/NumHask/Data.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Data.hs @@ -0,0 +1,94 @@ +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveFunctor, GeneralizedNewtypeDeriving #-} +module NumHask.Data where + +import GHC.Generics +import Data.Coerce (coerce) + +import NumHask.Algebra + +import Prelude hiding (Num(..), sum, recip) + +-- | Monoid under addition. +-- +-- >>> getSum (Sum 1 <> Sum 2 <> mempty) +-- 3 +newtype Sum a = Sum { getSum :: a } + deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Functor) + +-- | @since 4.8.0.0 +instance Applicative Sum where + pure = Sum + (<*>) = coerce + +-- | @since 4.8.0.0 +instance Monad Sum where + m >>= k = k (getSum m) + +instance AdditiveMagma a => AdditiveMagma (Sum a) where + (Sum x) `plus` (Sum y) = Sum (x `plus` y) + +instance AdditiveUnital a => AdditiveUnital (Sum a) where + zero = Sum zero + +instance AdditiveMagma a => AdditiveAssociative (Sum a) + +instance AdditiveInvertible a => AdditiveInvertible (Sum a) where + negate (Sum x) = Sum (negate x) + +instance AdditiveMagma a => AdditiveCommutative (Sum a) where + +instance (AdditiveUnital a, AdditiveMagma a) => Additive (Sum a) where + +instance (AdditiveInvertible a, AdditiveUnital a) => AdditiveGroup (Sum a) where + + +instance AdditiveMagma a => Semigroup (Sum a) where + (Sum x) <> (Sum y) = Sum $ x `plus` y + +instance AdditiveUnital a => Monoid (Sum a) where + mempty = Sum zero + + + + +-- | Monoid under multiplication. +-- +-- >>> getProduct (Product 3 <> Product 4 <> mempty) +-- 12 +newtype Product a = Product { getProduct :: a } + deriving (Eq, Ord, Read, Show, Bounded, Generic, Generic1, Functor) + +-- | @since 4.8.0.0 +instance Applicative Product where + pure = Product + (<*>) = coerce + +-- | @since 4.8.0.0 +instance Monad Product where + m >>= k = k (getProduct m) + + +instance MultiplicativeMagma a => MultiplicativeMagma (Product a) where + (Product x) `times` (Product y) = Product (x `times` y) + +instance MultiplicativeUnital a => MultiplicativeUnital (Product a) where + one = Product one + +instance MultiplicativeMagma a => MultiplicativeAssociative (Product a) + +instance MultiplicativeInvertible a => MultiplicativeInvertible (Product a) where + recip (Product x) = Product (recip x) + +instance MultiplicativeMagma a => MultiplicativeCommutative (Product a) + +instance MultiplicativeUnital a => Multiplicative (Product a) where + +instance (MultiplicativeUnital a, MultiplicativeInvertible a) => MultiplicativeGroup (Product a) where + + +instance MultiplicativeMagma a => Semigroup (Product a) where + (Product x) <> (Product y) = Product $ x `times` y + +instance MultiplicativeUnital a => Monoid (Product a) where + mempty = Product one diff --git a/src/NumHask/Data/Complex.hs b/src/NumHask/Data/Complex.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Data/Complex.hs @@ -0,0 +1,174 @@ +{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, DeriveFunctor, GeneralizedNewtypeDeriving, DeriveFoldable, DeriveTraversable #-} +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} +module NumHask.Data.Complex where + +import GHC.Generics (Generic, Generic1) +import Data.Data (Data) + +import NumHask.Algebra.Additive +import NumHask.Algebra.Multiplicative +import NumHask.Algebra.Ring +import NumHask.Algebra.Distribution +import NumHask.Algebra.Field +import NumHask.Algebra.Metric + +import Prelude hiding (Num(..), negate, sin, cos, sqrt, (/), atan, pi, exp, log, recip, (**)) +import qualified Prelude as P ( (&&), (>), (<=), (<), (==), otherwise, Ord(..) ) + +-- ----------------------------------------------------------------------------- +-- The Complex type + +infix 6 :+ + +-- | Complex numbers are an algebraic type. +-- +-- For a complex number @z@, @'abs' z@ is a number with the magnitude of @z@, +-- but oriented in the positive real direction, whereas @'sign' z@ +-- has the phase of @z@, but unit magnitude. +-- +-- The 'Foldable' and 'Traversable' instances traverse the real part first. +data Complex a + = !a :+ !a -- ^ forms a complex number from its real and imaginary + -- rectangular components. + deriving (Eq, Show, Read, Data, Generic, Generic1 + , Functor, Foldable, Traversable) + +-- | Extracts the real part of a complex number. +realPart :: Complex a -> a +realPart (x :+ _) = x + +-- | Extracts the imaginary part of a complex number. +imagPart :: Complex a -> a +imagPart (_ :+ y) = y + + + + +instance (AdditiveMagma a) => AdditiveMagma (Complex a) where + (rx :+ ix) `plus` (ry :+ iy) = (rx `plus` ry) :+ (ix `plus` iy) + +instance (AdditiveUnital a) => AdditiveUnital (Complex a) where + zero = zero :+ zero + +instance (AdditiveAssociative a) => AdditiveAssociative (Complex a) + +instance (AdditiveCommutative a) => AdditiveCommutative (Complex a) + +instance (Additive a) => Additive (Complex a) + +instance (AdditiveInvertible a) => AdditiveInvertible (Complex a) where + negate (rx :+ ix) = negate rx :+ negate ix + +instance (AdditiveGroup a) => AdditiveGroup (Complex a) + + +instance (Distribution a, AdditiveGroup a) => Distribution (Complex a) + + +instance (AdditiveUnital a, AdditiveGroup a, MultiplicativeUnital a) => MultiplicativeUnital (Complex a) where + one = one :+ zero + +instance (MultiplicativeMagma a, AdditiveGroup a) => MultiplicativeMagma (Complex a) where + (rx :+ ix) `times` (ry :+ iy) = + (rx `times` ry - ix `times` iy) :+ (ix `times` ry + iy `times` rx) + +instance (MultiplicativeMagma a, AdditiveGroup a) => MultiplicativeCommutative (Complex a) + +instance (MultiplicativeUnital a, MultiplicativeAssociative a, AdditiveGroup a) => Multiplicative (Complex a) + + +instance (AdditiveGroup a, MultiplicativeInvertible a) => MultiplicativeInvertible (Complex a) where + recip (rx :+ ix) = (rx `times` d) :+ (negate ix `times` d) + where + d = recip ((rx `times` rx) `plus` (ix `times` ix)) + + + +instance (MultiplicativeUnital a, MultiplicativeAssociative a, MultiplicativeInvertible a, AdditiveGroup a) => MultiplicativeGroup (Complex a) + + + + +instance (AdditiveGroup a, MultiplicativeAssociative a) => + MultiplicativeAssociative (Complex a) + + +instance (Semiring a, AdditiveGroup a) => Semiring (Complex a) + +instance (Semiring a, AdditiveGroup a) => Ring (Complex a) + +instance (Semiring a, AdditiveGroup a) => InvolutiveRing (Complex a) + +instance (MultiplicativeAssociative a, MultiplicativeUnital a, AdditiveGroup a, Semiring a) => + CRing (Complex a) + +instance (MultiplicativeGroup a, AdditiveGroup a, Semiring a) => Field (Complex a) + + +instance (Multiplicative a, ExpField a, Normed a a) => + Normed (Complex a) a where + normL1 (rx :+ ix) = normL1 rx + normL1 ix + normL2 (rx :+ ix) = sqrt (rx * rx + ix * ix) + normLp p (rx :+ ix) = (normL1 rx ** p + normL1 ix ** p) ** (one / p) + +instance (Multiplicative a, ExpField a, Normed a a) => Metric (Complex a) a where + distanceL1 a b = normL1 (a - b) + distanceL2 a b = normL2 (a - b) + distanceLp p a b = normLp p (a - b) + + + +-- | todo: bottom is here somewhere??? +instance (P.Ord a, TrigField a, ExpField a) => ExpField (Complex a) where + exp (rx :+ ix) = exp rx * cos ix :+ exp rx * sin ix + log (rx :+ ix) = log (sqrt (rx * rx + ix * ix)) :+ atan2' ix rx + where + atan2' y x + | x P.> zero = atan (y / x) + | x P.== zero P.&& y P.> zero = pi / (one + one) + | x P.< one P.&& y P.> one = pi + atan (y / x) + | (x P.<= zero P.&& y P.< zero) || (x P.< zero) = + negate (atan2' (negate y) x) + | y P.== zero = pi -- must be after the previous test on zero y + | x P.== zero P.&& y P.== zero = y -- must be after the other double zero tests + | P.otherwise = x + y -- x or y is a NaN, return a NaN (via +) + + + + + + +-- * Helpers from Data.Complex + +mkPolar :: TrigField a => a -> a -> Complex a +mkPolar r theta = r * cos theta :+ r * sin theta + + +-- | @'cis' t@ is a complex value with magnitude @1@ +-- and phase @t@ (modulo @2*'pi'@). +{-# SPECIALISE cis :: Double -> Complex Double #-} +cis :: TrigField a => a -> Complex a +cis theta = cos theta :+ sin theta + +-- | The function 'polar' takes a complex number and +-- returns a (magnitude, phase) pair in canonical form: +-- the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@; +-- if the magnitude is zero, then so is the phase. +{-# SPECIALISE polar :: Complex Double -> (Double,Double) #-} +polar :: (RealFloat a, ExpField a) => Complex a -> (a,a) +polar z = (magnitude z, phase z) + +-- | The nonnegative magnitude of a complex number. +{-# SPECIALISE magnitude :: Complex Double -> Double #-} +magnitude :: (ExpField a, RealFloat a) => Complex a -> a +magnitude (x :+ y) = scaleFloat k (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y))) + where k = max (exponent x) (exponent y) + mk = - k + sqr z = z * z + +-- | The phase of a complex number, in the range @(-'pi', 'pi']@. +-- If the magnitude is zero, then so is the phase. +{-# SPECIALISE phase :: Complex Double -> Double #-} +phase :: (RealFloat a) => Complex a -> a +phase (0 :+ 0) = 0 -- SLPJ July 97 from John Peterson +phase (x :+ y) = atan2 y x diff --git a/src/NumHask/Data/LogField.hs b/src/NumHask/Data/LogField.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Data/LogField.hs @@ -0,0 +1,360 @@ +{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, DeriveFunctor, GeneralizedNewtypeDeriving, DeriveFoldable, DeriveTraversable, GADTs #-} +{-# LANGUAGE FlexibleInstances, FlexibleContexts, UndecidableInstances, MultiParamTypeClasses #-} +module NumHask.Data.LogField + ( + -- * @LogField@ + LogField() + -- ** Isomorphism to normal-domain + , logField + , fromLogField + -- ** Isomorphism to log-domain + , logToLogField + , logFromLogField + -- ** Additional operations + , accurateSum, accurateProduct + , pow + )where + +import GHC.Generics ( Generic + , Generic1 + ) +import Data.Data ( Data ) + +import NumHask.Algebra.Additive +import NumHask.Algebra.Multiplicative +import NumHask.Algebra.Distribution +import NumHask.Algebra.Field +import NumHask.Algebra.Integral +import NumHask.Algebra.Rational +import NumHask.Algebra.Metric + +import Prelude hiding ( Num(..) + , negate + , sin + , cos + , sqrt + , (/) + , atan + , pi + , exp + , log + , recip + , (**) + , toInteger + ) +import qualified Data.Foldable as F + +-- LogField is adapted from LogFloat +---------------------------------------------------------------- +-- ~ 2015.08.06 +-- | +-- Module : Data.Number.LogFloat +-- Copyright : Copyright (c) 2007--2015 wren gayle romano +-- License : BSD3 +-- Maintainer : wren@community.haskell.org +-- Stability : stable +-- Portability : portable (with CPP, FFI) +-- Link : https://hackage.haskell.org/package/logfloat +---------------------------------------------------------------- + +---------------------------------------------------------------- +-- +-- | A @LogField@ is just a 'Field' with a special interpretation. +-- The 'LogField' function is presented instead of the constructor, +-- in order to ensure semantic conversion. At present the 'Show' +-- instance will convert back to the normal-domain, and hence will +-- underflow at that point. This behavior may change in the future. +-- +-- Because 'logField' performs the semantic conversion, we can use +-- operators which say what we *mean* rather than saying what we're +-- actually doing to the underlying representation. That is, +-- equivalences like the following are true[1] thanks to type-class +-- overloading: +-- +-- > logField (p + q) == logField p + logField q +-- > logField (p * q) == logField p * logField q +-- +-- +-- Performing operations in the log-domain is cheap, prevents +-- underflow, and is otherwise very nice for dealing with miniscule +-- probabilities. However, crossing into and out of the log-domain +-- is expensive and should be avoided as much as possible. In +-- particular, if you're doing a series of multiplications as in +-- @lp * LogField q * LogField r@ it's faster to do @lp * LogField +-- (q * r)@ if you're reasonably sure the normal-domain multiplication +-- won't underflow; because that way you enter the log-domain only +-- once, instead of twice. Also note that, for precision, if you're +-- doing more than a few multiplications in the log-domain, you +-- should use 'product' rather than using '(*)' repeatedly. +-- +-- Even more particularly, you should /avoid addition/ whenever +-- possible. Addition is provided because sometimes we need it, and +-- the proper implementation is not immediately apparent. However, +-- between two @LogField@s addition requires crossing the exp\/log +-- boundary twice; with a @LogField@ and a 'Double' it's three +-- times, since the regular number needs to enter the log-domain +-- first. This makes addition incredibly slow. Again, if you can +-- parenthesize to do normal-domain operations first, do it! +-- +-- [1] That is, true up-to underflow and floating point fuzziness. +-- Which is, of course, the whole point of this module. +newtype LogField a = LogField a + deriving (Eq, Ord, Read, Data, Generic, Generic1, Functor, Foldable, Traversable) + +---------------------------------------------------------------- +-- To show it, we want to show the normal-domain value rather than +-- the log-domain value. Also, if someone managed to break our +-- invariants (e.g. by passing in a negative and noone's pulled on +-- the thunk yet) then we want to crash before printing the +-- constructor, rather than after. N.B. This means the show will +-- underflow\/overflow in the same places as normal doubles since +-- we underflow at the @exp@. Perhaps this means we should show the +-- log-domain value instead. + +instance (ExpField a, Show a) => Show (LogField a) where + showsPrec p (LogField x) = + let y = exp x in y `seq` + showParen (p > 9) + ( showString "LogField " + . showsPrec 11 y + ) + +---------------------------------------------------------------- +-- | Constructor which does semantic conversion from normal-domain +-- to log-domain. Throws errors on negative and NaN inputs. If @p@ +-- is non-negative, then following equivalence holds: +-- +-- > logField p == logToLogField (log p) +logField :: (ExpField a) => a -> LogField a +{-# INLINE [0] logField #-} +logField = LogField . log + + +-- TODO: figure out what to do here, removed guards +-- | Constructor which assumes the argument is already in the +-- log-domain. +logToLogField :: a -> LogField a +logToLogField = LogField + + +-- | Semantically convert our log-domain value back into the +-- normal-domain. Beware of overflow\/underflow. The following +-- equivalence holds (without qualification): +-- +-- > fromLogField == exp . logFromLogField +-- +fromLogField :: ExpField a => LogField a -> a +{-# INLINE [0] fromLogField #-} +fromLogField (LogField x) = exp x + + +-- | Return the log-domain value itself without conversion. +logFromLogField :: LogField a -> a +logFromLogField (LogField x) = x + + +-- These are our module-specific versions of "log\/exp" and "exp\/log"; +-- They do the same things but also have a @LogField@ in between +-- the logarithm and exponentiation. In order to ensure these rules +-- fire, we have to delay the inlining on two of the four +-- con-\/destructors. + +{-# RULES +-- Out of log-domain and back in +"log/fromLogField" forall x. log (fromLogField x) = logFromLogField x +-- TODO: Rewrite-rule too complicated +"LogField/fromLogField" forall x. LogField (fromLogField x) = x + +-- Into log-domain and back out +"fromLogField/LogField" forall x. fromLogField (LogField x) = x + #-} + + +log1p :: ExpField a => a -> a +{-# INLINE [0] log1p #-} +log1p x = log (one + x) + +expm1 :: ExpField a => a -> a +{-# INLINE [0] expm1 #-} +expm1 x = exp x - one + +{-# RULES +-- Into log-domain and back out +"expm1/log1p" forall x. expm1 (log1p x) = x + +-- Out of log-domain and back in +"log1p/expm1" forall x. log1p (expm1 x) = x + #-} + +instance (ExpField a, LowerBoundedField a, Ord a) => AdditiveMagma (LogField a) where + x@(LogField x') `plus` y@(LogField y') + | x == zero && y == zero = zero + | x == zero = y + | y == zero = x + | x >= y = LogField (x' + log1p (exp (y' - x'))) + | otherwise = LogField (y' + log1p (exp (x' - y'))) + +instance (LowerBoundedField a, ExpField a, Ord a) => AdditiveUnital (LogField a) where + zero = LogField negInfinity + +instance (LowerBoundedField a, ExpField a, Ord a) => AdditiveAssociative (LogField a) + +instance (LowerBoundedField a,ExpField a, Ord a) => AdditiveCommutative (LogField a) + +instance (LowerBoundedField a, ExpField a, Ord a) => Additive (LogField a) + +instance (AdditiveMagma a, LowerBoundedField a, Eq a) => MultiplicativeMagma (LogField a) where + (LogField x) `times ` (LogField y) + | x == negInfinity || y == negInfinity = LogField negInfinity + | otherwise = LogField (x `plus` y) + +instance (AdditiveUnital a, LowerBoundedField a, Eq a) => MultiplicativeUnital (LogField a) where + one = LogField zero + +instance (AdditiveAssociative a, LowerBoundedField a, Eq a) => MultiplicativeAssociative (LogField a) + +instance (AdditiveCommutative a, LowerBoundedField a, Eq a) => MultiplicativeCommutative (LogField a) + +instance (AdditiveInvertible a, LowerBoundedField a, Eq a) => MultiplicativeInvertible (LogField a) where + recip (LogField x) = LogField $ negate x + +instance (AdditiveUnital a + , AdditiveAssociative a + , AdditiveCommutative a + , Additive a + , LowerBoundedField a + , Eq a) => Multiplicative (LogField a) + +instance (AdditiveUnital a + , AdditiveAssociative a + , AdditiveInvertible a + , AdditiveLeftCancellative a + , LowerBoundedField a + , Eq a) => MultiplicativeLeftCancellative (LogField a) + +instance (AdditiveUnital a + , AdditiveAssociative a + , AdditiveInvertible a + , AdditiveRightCancellative a + , LowerBoundedField a + , Eq a) => MultiplicativeRightCancellative (LogField a) + +instance (Multiplicative (LogField a), AdditiveInvertible a, AdditiveGroup a, LowerBoundedField a, Eq a) => MultiplicativeGroup (LogField a) + +instance (LowerBoundedField a, ExpField a, Ord a, AdditiveMagma a) => Distribution (LogField a) + +-- unable to provide this instance because there is no Field (LogField a) instance +-- instance (Field (LogField a), ExpField a, LowerBoundedField a, Ord a) => ExpField (LogField a) where +-- exp (LogField x) = (LogField $ exp x) +-- log (LogField x) = (LogField $ log x) +-- (**) x (LogField y) = pow x $ exp y + +instance (FromInteger a, ExpField a) => FromInteger (LogField a) where + fromInteger = logField . fromInteger + +instance (ToInteger a, ExpField a) => ToInteger (LogField a) where + toInteger = toInteger . fromLogField + +instance (FromRatio a, ExpField a) => FromRatio (LogField a) where + fromRatio = logField . fromRatio + +instance (ToRatio a, ExpField a) => ToRatio (LogField a) where + toRatio = toRatio . fromLogField + +instance (Epsilon a, ExpField a, LowerBoundedField a, Ord a) => Epsilon (LogField a) where + nearZero (LogField x) = nearZero $ exp x + aboutEqual (LogField x) (LogField y) = aboutEqual (exp x) (exp y) + + +---------------------------------------------------------------- +-- | /O(1)/. Compute powers in the log-domain; that is, the following +-- equivalence holds (modulo underflow and all that): +-- +-- > LogField (p ** m) == LogField p `pow` m +-- +-- /Since: 0.13/ +pow :: (ExpField a, LowerBoundedField a, Ord a) => LogField a -> a -> LogField a +{-# INLINE pow #-} +infixr 8 `pow` +pow x@(LogField x') m + | x == zero && m == zero = LogField zero + | x == zero = x + | otherwise = LogField $ m * x' + + +-- Some good test cases: +-- for @logsumexp == log . accurateSum . map exp@: +-- logsumexp[0,1,0] should be about 1.55 +-- for correctness of avoiding underflow: +-- logsumexp[1000,1001,1000] ~~ 1001.55 == 1000 + 1.55 +-- logsumexp[-1000,-999,-1000] ~~ -998.45 == -1000 + 1.55 +-- +-- | /O(n)/. Compute the sum of a finite list of 'LogField's, being +-- careful to avoid underflow issues. That is, the following +-- equivalence holds (modulo underflow and all that): +-- +-- > LogField . accurateSum == accurateSum . map LogField +-- +-- /N.B./, this function requires two passes over the input. Thus, +-- it is not amenable to list fusion, and hence will use a lot of +-- memory when summing long lists. +{-# INLINE accurateSum #-} +accurateSum :: (ExpField a, Foldable f, Ord a) => f (LogField a) -> LogField a +accurateSum xs = LogField (theMax + log theSum) + where + LogField theMax = maximum xs + + -- compute @\log \sum_{x \in xs} \exp(x - theMax)@ + theSum = F.foldl' (\acc (LogField x) -> acc + exp (x - theMax)) zero xs + +-- | /O(n)/. Compute the product of a finite list of 'LogField's, +-- being careful to avoid numerical error due to loss of precision. +-- That is, the following equivalence holds (modulo underflow and +-- all that): +-- +-- > LogField . accurateProduct == accurateProduct . map LogField +{-# INLINE accurateProduct #-} +accurateProduct :: (ExpField a, Foldable f) => f (LogField a) -> LogField a +accurateProduct = LogField . fst . F.foldr kahanPlus (zero, zero) + where + kahanPlus (LogField x) (t, c) = + let y = x - c + t' = t + y + c' = (t' - t) - y + in (t', c') + +-- This version *completely* eliminates rounding errors and loss +-- of significance due to catastrophic cancellation during summation. +-- Also see the other +-- implementations given there. For Python's actual C implementation, +-- see math_fsum in +-- +-- +-- For merely *mitigating* errors rather than completely eliminating +-- them, see . +-- +-- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@ +{- +-- For proof of correctness, see +-- +def msum(xs): + partials = [] # sorted, non-overlapping partial sums + # N.B., the actual C implementation uses a 32 array, doubling size as needed + for x in xs: + i = 0 + for y in partials: # for(i = j = 0; j < n; j++) + if abs(x) < abs(y): + x, y = y, x + hi = x + y + lo = y - (hi - x) + if lo != 0.0: + partials[i] = lo + i += 1 + x = hi + # does an append of x while dropping all the partials after + # i. The C version does n=i; and leaves the garbage in place + partials[i:] = [x] + # BUG: this last step isn't entirely correct and can lose + # precision + return sum(partials, 0.0) +-}