numhask 0.2.2.0 → 0.2.3.0
raw patch · 8 files changed
+662/−10 lines, 8 files
Files
- numhask.cabal +4/−2
- src/NumHask/Algebra.hs +2/−2
- src/NumHask/Algebra/Field.hs +5/−1
- src/NumHask/Algebra/Metric.hs +4/−2
- src/NumHask/Algebra/Rational.hs +19/−3
- src/NumHask/Data.hs +94/−0
- src/NumHask/Data/Complex.hs +174/−0
- src/NumHask/Data/LogField.hs +360/−0
numhask.cabal view
@@ -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
src/NumHask/Algebra.hs view
@@ -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
src/NumHask/Algebra/Field.hs view
@@ -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
src/NumHask/Algebra/Metric.hs view
@@ -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
src/NumHask/Algebra/Rational.hs view
@@ -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
+ src/NumHask/Data.hs view
@@ -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
+ src/NumHask/Data/Complex.hs view
@@ -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
+ src/NumHask/Data/LogField.hs view
@@ -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.+-- <http://code.activestate.com/recipes/393090/> Also see the other+-- implementations given there. For Python's actual C implementation,+-- see math_fsum in+-- <http://svn.python.org/view/python/trunk/Modules/mathmodule.c?view=markup>+--+-- For merely *mitigating* errors rather than completely eliminating+-- them, see <http://code.activestate.com/recipes/298339/>.+--+-- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@+{-+-- For proof of correctness, see+-- <www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps>+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 <http://stackoverflow.com/a/2704565/358069>+ return sum(partials, 0.0)+-}