module Numeric.Log
( Log(..)
, Precise(..)
, sum
) where
import Prelude hiding (maximum, sum)
import Control.Applicative
import Control.Comonad
import Control.DeepSeq
import Data.Binary as Binary
import Data.Complex
import Data.Data
import Data.Distributive
import Data.Foldable as Foldable hiding (sum)
import Data.Functor.Bind
import Data.Functor.Extend
import Data.Hashable
import Data.Int
import Data.List as List hiding (sum)
import Data.Monoid
import Data.SafeCopy
import Data.Semigroup.Foldable
import Data.Semigroup.Traversable
import Data.Serialize as Serialize
import Data.Traversable
import Foreign.Ptr
import Foreign.Storable
import Generics.Deriving
import Text.Read
newtype Log a = Log { runLog :: a } deriving (Eq,Ord,Data,Typeable,Generic)
deriveSafeCopy 1 'base ''Log
instance (Floating a, Show a) => Show (Log a) where
showsPrec d (Log a) = showsPrec d (exp a)
instance (Floating a, Read a) => Read (Log a) where
readPrec = Log . log <$> step readPrec
instance Binary a => Binary (Log a) where
put = Binary.put . runLog
get = Log <$> Binary.get
instance Serialize a => Serialize (Log a) where
put = Serialize.put . runLog
get = Log <$> Serialize.get
instance Functor Log where
fmap f (Log a) = Log (f a)
instance Hashable a => Hashable (Log a) where
hashWithSalt i (Log a) = hashWithSalt i a
instance Storable a => Storable (Log a) where
sizeOf = sizeOf . runLog
alignment = alignment . runLog
peek ptr = Log <$> peek (castPtr ptr)
poke ptr (Log a) = poke (castPtr ptr) a
instance NFData a => NFData (Log a) where
rnf (Log a) = rnf a
instance Foldable Log where
foldMap f (Log a) = f a
instance Foldable1 Log where
foldMap1 f (Log a) = f a
instance Traversable Log where
traverse f (Log a) = Log <$> f a
instance Traversable1 Log where
traverse1 f (Log a) = Log <$> f a
instance Distributive Log where
distribute = Log . fmap runLog
instance Extend Log where
extended f w@Log{} = Log (f w)
instance Comonad Log where
extract (Log a) = a
extend f w@Log{} = Log (f w)
instance Applicative Log where
pure = Log
Log f <*> Log a = Log (f a)
instance ComonadApply Log where
Log f <@> Log a = Log (f a)
instance Apply Log where
Log f <.> Log a = Log (f a)
instance Bind Log where
Log a >>- f = f a
instance Monad Log where
return = Log
Log a >>= f = f a
instance (RealFloat a, Precise a, Enum a) => Enum (Log a) where
succ a = a + 1
pred a = a 1
toEnum = fromIntegral
fromEnum = round . exp . runLog
enumFrom (Log a) = [ Log (log b) | b <- enumFrom (exp a) ]
enumFromThen (Log a) (Log b) = [ Log (log c) | c <- enumFromThen (exp a) (exp b) ]
enumFromTo (Log a) (Log b) = [ Log (log c) | c <- enumFromTo (exp a) (exp b) ]
enumFromThenTo (Log a) (Log b) (Log c) = [ Log (log d) | d <- enumFromThenTo (exp a) (exp b) (exp c) ]
negInf :: Fractional a => a
negInf = (1/0)
instance (Precise a, RealFloat a) => Num (Log a) where
Log a * Log b
| isInfinite a && isInfinite b && a == b = Log negInf
| otherwise = Log (a + b)
Log a + Log b
| a == b && isInfinite a && isInfinite b = Log a
| a >= b = Log (a + log1p (exp (b a)))
| otherwise = Log (b + log1p (exp (a b)))
Log a Log b
| a == negInf && b == negInf = Log negInf
| otherwise = Log (a + log1p (negate (exp (b a))))
signum (Log a)
| a == negInf = 0
| a > negInf = 1
| otherwise = negInf
negate _ = negInf
abs = id
fromInteger = Log . log . fromInteger
instance (Precise a, RealFloat a, Eq a) => Fractional (Log a) where
Log a / Log b
| a == b && isInfinite a && isInfinite b = Log negInf
| a == negInf = Log negInf
| otherwise = Log (ab)
fromRational = Log . log . fromRational
instance (Precise a, RealFloat a, Ord a) => Real (Log a) where
toRational (Log a) = toRational (exp a)
data Acc1 a = Acc1 !Int64 !a
instance (Precise a, RealFloat a) => Monoid (Log a) where
mempty = Log negInf
mappend = (+)
mconcat [] = 0
mconcat (Log z:zs) = Log $ case List.foldl' step1 (Acc1 0 z) zs of
Acc1 nm1 a
| isInfinite a -> a
| otherwise -> a + log1p (List.foldl' (step2 a) 0 zs + fromIntegral nm1)
where
step1 (Acc1 n y) (Log x) = Acc1 (n + 1) (max x y)
step2 a r (Log x) = r + expm1 (x a)
logMap :: Floating a => (a -> a) -> Log a -> Log a
logMap f = Log . log . f . exp . runLog
data Acc a = Acc !Int64 !a | None
sum :: (RealFloat a, Ord a, Precise a, Foldable f) => f (Log a) -> Log a
sum xs = Log $ case Foldable.foldl' step1 None xs of
None -> negInf
Acc nm1 a
| isInfinite a -> a
| otherwise -> a + log1p (Foldable.foldl' (step2 a) 0 xs + fromIntegral nm1)
where
step1 None (Log x) = Acc 0 x
step1 (Acc n y) (Log x) = Acc (n + 1) (max x y)
step2 a r (Log x) = r + expm1 (x a)
instance (RealFloat a, Precise a) => Floating (Log a) where
pi = Log (log pi)
exp (Log a) = Log (exp a)
log (Log a) = Log (log a)
sqrt (Log a) = Log (a / 2)
logBase (Log a) (Log b) = Log (log (logBase (exp a) (exp b)))
sin = logMap sin
cos = logMap cos
tan = logMap tan
asin = logMap asin
acos = logMap acos
atan = logMap atan
sinh = logMap sinh
cosh = logMap cosh
tanh = logMap tanh
asinh = logMap asinh
acosh = logMap acosh
atanh = logMap atanh
class Floating a => Precise a where
log1p :: a -> a
expm1 :: a -> a
instance Precise Double where
log1p = c_log1p
expm1 = c_expm1
instance Precise Float where
log1p = c_log1pf
expm1 = c_expm1f
instance (RealFloat a, Precise a) => Precise (Complex a) where
expm1 x@(a :+ b)
| a*a + b*b < 1, u <- expm1 a, v <- sin (b/2), w <- 2*v*v = (u*w+u+w) :+ (u+1)*sin b
| otherwise = exp x 1
log1p x@(a :+ b)
| abs a < 0.5 && abs b < 0.5, u <- 2*a+a*a+b*b = log1p (u/(1+sqrt (u+1))) :+ atan2 (1 + a) b
| otherwise = log (1 + x)
foreign import ccall unsafe "math.h log1p" c_log1p :: Double -> Double
foreign import ccall unsafe "math.h expm1" c_expm1 :: Double -> Double
foreign import ccall unsafe "math.h expm1f" c_expm1f :: Float -> Float
foreign import ccall unsafe "math.h log1pf" c_log1pf :: Float -> Float