-- CPP and GeneralizedNewtypeDeriving are needed for IArray UArray instance
-- FFI is for log1p
{-# LANGUAGE CPP, ForeignFunctionInterface, MultiParamTypeClasses #-}
-- We don't put these in LANGUAGE, because it's CPP guarded for GHC only
-- HACK: ScopedTypeVariables and InstanceSigs are for GHC 7.10 only...
{-# OPTIONS_GHC
    -XGeneralizedNewtypeDeriving
    -XScopedTypeVariables
    -XInstanceSigs
    #-}

{-# OPTIONS_GHC -Wall -fwarn-tabs #-}

{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-}

----------------------------------------------------------------
--                                                  ~ 2021.10.17
-- |
-- Module      :  Data.Number.LogFloat
-- Copyright   :  Copyright (c) 2007--2021 wren gayle romano
-- License     :  BSD3
-- Maintainer  :  wren@cpan.org
-- Stability   :  stable
-- Portability :  portable (with CPP, FFI)
--
-- This module presents a type for storing numbers in the log-domain.
-- The main reason for doing this is to prevent underflow when
-- multiplying many small probabilities as is done in Hidden Markov
-- Models and other statistical models often used for natural
-- language processing. The log-domain also helps prevent overflow
-- when multiplying many large numbers. In rare cases it can speed
-- up numerical computation (since addition is faster than
-- multiplication, though logarithms are exceptionally slow), but
-- the primary goal is to improve accuracy of results. A secondary
-- goal has been to maximize efficiency since these computations
-- are frequently done within a /O(n^3)/ loop.
--
-- The 'LogFloat' of this module is restricted to non-negative
-- numbers for efficiency's sake, see "Data.Number.LogFloat.Signed"
-- for doing signed log-domain calculations.
----------------------------------------------------------------

module Data.Number.LogFloat
    (
    -- * Exceptional numeric values
      module Data.Number.Transfinite

    -- * @LogFloat@ data type
    , LogFloat()
    -- ** Isomorphism to normal-domain
    , logFloat
    , fromLogFloat
    -- ** Isomorphism to log-domain
    , logToLogFloat
    , logFromLogFloat
    -- ** Additional operations
    , sum, product
    , pow

    -- * Accurate versions of logarithm\/exponentiation
    , log1p, expm1
    ) where

import Prelude hiding (log, sum, product, isInfinite, isNaN)

import Data.Number.Transfinite
import Data.Number.PartialOrd
import Data.Number.LogFloat.Raw


-- GHC can derive (IArray UArray LogFloat), but Hugs needs to coerce
-- TODO: see about nhc98/yhc, jhc/lhc
import Data.Array.Base    (IArray(..))
import Data.Array.Unboxed (UArray)

-- Hugs (Sept 2006) doesn't use the generic wrapper in base:Unsafe.Coerce
-- so we'll just have to go back to the original source.
#ifdef __HUGS__
import Hugs.IOExts        (unsafeCoerce)
#elif __NHC__
import NonStdUnsafeCoerce (unsafeCoerce)
#elif __GLASGOW_HASKELL__ >= 710
-- For when the *heap* representations are the same
--import Data.Coerce        (coerce)
-- For when the *unboxed array* storage representations are the same
import Unsafe.Coerce      (unsafeCoerce)
import Data.Ix            (Ix)
#endif

#ifdef __GLASGOW_HASKELL__
import Foreign.Storable (Storable)
#endif

----------------------------------------------------------------
--
-- | A @LogFloat@ is just a 'Double' with a special interpretation.
-- The 'logFloat' 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.
-- At present, the 'Read' instance parses things in the normal-domain
-- and then converts them to the log-domain. Again, this behavior
-- may change in the future.
--
-- Because 'logFloat' 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:
--
-- > logFloat (p + q) == logFloat p + logFloat q
-- > logFloat (p * q) == logFloat p * logFloat q
--
-- (Do note, however, that subtraction can and negation will throw
-- errors: since @LogFloat@ can only represent the positive half of
-- 'Double'. 'Num' is the wrong abstraction to put at the bottom
-- of the numeric type-class hierarchy; but alas, we're stuck with
-- it.)
--
-- 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 * logFloat q * logFloat r@ it's faster to do @lp * logFloat
-- (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 @LogFloat@s addition requires crossing the exp\/log
-- boundary twice; with a @LogFloat@ 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 LogFloat = LogFloat Double
    deriving
    ( LogFloat -> LogFloat -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: LogFloat -> LogFloat -> Bool
$c/= :: LogFloat -> LogFloat -> Bool
== :: LogFloat -> LogFloat -> Bool
$c== :: LogFloat -> LogFloat -> Bool
Eq
    , Eq LogFloat
LogFloat -> LogFloat -> Bool
LogFloat -> LogFloat -> Ordering
LogFloat -> LogFloat -> LogFloat
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: LogFloat -> LogFloat -> LogFloat
$cmin :: LogFloat -> LogFloat -> LogFloat
max :: LogFloat -> LogFloat -> LogFloat
$cmax :: LogFloat -> LogFloat -> LogFloat
>= :: LogFloat -> LogFloat -> Bool
$c>= :: LogFloat -> LogFloat -> Bool
> :: LogFloat -> LogFloat -> Bool
$c> :: LogFloat -> LogFloat -> Bool
<= :: LogFloat -> LogFloat -> Bool
$c<= :: LogFloat -> LogFloat -> Bool
< :: LogFloat -> LogFloat -> Bool
$c< :: LogFloat -> LogFloat -> Bool
compare :: LogFloat -> LogFloat -> Ordering
$ccompare :: LogFloat -> LogFloat -> Ordering
Ord -- Should we really perpetuate the Ord lie?
#ifdef __GLASGOW_HASKELL__
    -- The H98 Report doesn't include these among the options for
    -- automatic derivation. But GHC >= 6.8.2 (at least) can derive
    -- them (even without GeneralizedNewtypeDeriving). However,
    -- GHC >= 7.10 can't derive @IArray UArray@ thanks to the new
    -- type-role stuff! since 'UArray' is declared to be nominal
    -- in both arguments... and that seems to be necessary:
    -- cf., <https://ghc.haskell.org/trac/ghc/ticket/9220>
#if __GLASGOW_HASKELL__ < 710
    , IArray UArray
#endif
    , Ptr LogFloat -> IO LogFloat
Ptr LogFloat -> Int -> IO LogFloat
Ptr LogFloat -> Int -> LogFloat -> IO ()
Ptr LogFloat -> LogFloat -> IO ()
LogFloat -> Int
forall b. Ptr b -> Int -> IO LogFloat
forall b. Ptr b -> Int -> LogFloat -> IO ()
forall a.
(a -> Int)
-> (a -> Int)
-> (Ptr a -> Int -> IO a)
-> (Ptr a -> Int -> a -> IO ())
-> (forall b. Ptr b -> Int -> IO a)
-> (forall b. Ptr b -> Int -> a -> IO ())
-> (Ptr a -> IO a)
-> (Ptr a -> a -> IO ())
-> Storable a
poke :: Ptr LogFloat -> LogFloat -> IO ()
$cpoke :: Ptr LogFloat -> LogFloat -> IO ()
peek :: Ptr LogFloat -> IO LogFloat
$cpeek :: Ptr LogFloat -> IO LogFloat
pokeByteOff :: forall b. Ptr b -> Int -> LogFloat -> IO ()
$cpokeByteOff :: forall b. Ptr b -> Int -> LogFloat -> IO ()
peekByteOff :: forall b. Ptr b -> Int -> IO LogFloat
$cpeekByteOff :: forall b. Ptr b -> Int -> IO LogFloat
pokeElemOff :: Ptr LogFloat -> Int -> LogFloat -> IO ()
$cpokeElemOff :: Ptr LogFloat -> Int -> LogFloat -> IO ()
peekElemOff :: Ptr LogFloat -> Int -> IO LogFloat
$cpeekElemOff :: Ptr LogFloat -> Int -> IO LogFloat
alignment :: LogFloat -> Int
$calignment :: LogFloat -> Int
sizeOf :: LogFloat -> Int
$csizeOf :: LogFloat -> Int
Storable
#endif
    )

#if __GLASGOW_HASKELL__ >= 710
-- TODO: this version should also work for NHC and Hugs, I think...
-- HACK: we should be able to just unsafeCoerce the functions
-- themselves, instead of coercing the inputs and the outputs; but,
-- GHC 7.10 seems to get confused about trying to coerce the index
-- types too... To fix this we give explicit signatures, as below,
-- but this requires both ScopedTypeVariables and InstanceSigs; and
-- I'm not sure when InstanceSigs was introduced.

instance IArray UArray LogFloat where
    {-# INLINE bounds #-}
    bounds :: forall i. Ix i => UArray i LogFloat -> (i, i)
    bounds :: forall i. Ix i => UArray i LogFloat -> (i, i)
bounds = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
bounds :: UArray i Double -> (i, i))

    {-# INLINE numElements #-}
    numElements :: forall i. Ix i => UArray i LogFloat -> Int
    numElements :: forall i. Ix i => UArray i LogFloat -> Int
numElements = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i. (IArray a e, Ix i) => a i e -> Int
numElements :: UArray i Double -> Int)

    {-# INLINE unsafeArray #-}
    unsafeArray
        :: forall i. Ix i => (i,i) -> [(Int,LogFloat)] -> UArray i LogFloat
    unsafeArray :: forall i. Ix i => (i, i) -> [(Int, LogFloat)] -> UArray i LogFloat
unsafeArray = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [(Int, e)] -> a i e
unsafeArray
        :: (i,i) -> [(Int,Double)] -> UArray i Double)

    {-# INLINE unsafeAt #-}
    unsafeAt :: forall i. Ix i => UArray i LogFloat -> Int -> LogFloat
    unsafeAt :: forall i. Ix i => UArray i LogFloat -> Int -> LogFloat
unsafeAt = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt :: UArray i Double -> Int -> Double)

    {-# INLINE unsafeReplace #-}
    unsafeReplace
        :: forall i. Ix i
        => UArray i LogFloat -> [(Int,LogFloat)] -> UArray i LogFloat
    unsafeReplace :: forall i.
Ix i =>
UArray i LogFloat -> [(Int, LogFloat)] -> UArray i LogFloat
unsafeReplace = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> [(Int, e)] -> a i e
unsafeReplace
        :: UArray i Double -> [(Int,Double)] -> UArray i Double)

    {-# INLINE unsafeAccum #-}
    unsafeAccum
        :: forall i e. Ix i
        => (LogFloat -> e -> LogFloat)
        -> UArray i LogFloat -> [(Int,e)] -> UArray i LogFloat
    unsafeAccum :: forall i e'.
Ix i =>
(LogFloat -> e' -> LogFloat)
-> UArray i LogFloat -> [(Int, e')] -> UArray i LogFloat
unsafeAccum = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i e'.
(IArray a e, Ix i) =>
(e -> e' -> e) -> a i e -> [(Int, e')] -> a i e
unsafeAccum
        :: (Double -> e -> Double)
        -> UArray i Double -> [(Int,e)] -> UArray i Double)

    {-# INLINE unsafeAccumArray #-}
    unsafeAccumArray
        :: forall i e. Ix i
        => (LogFloat -> e -> LogFloat)
        -> LogFloat -> (i,i) -> [(Int,e)] -> UArray i LogFloat
    unsafeAccumArray :: forall i e'.
Ix i =>
(LogFloat -> e' -> LogFloat)
-> LogFloat -> (i, i) -> [(Int, e')] -> UArray i LogFloat
unsafeAccumArray = forall a b. a -> b
unsafeCoerce (forall (a :: * -> * -> *) e i e'.
(IArray a e, Ix i) =>
(e -> e' -> e) -> e -> (i, i) -> [(Int, e')] -> a i e
unsafeAccumArray
        :: (Double -> e -> Double)
        -> Double -> (i,i) -> [(Int,e)] -> UArray i Double)

#elif __HUGS__ || __NHC__
-- TODO: Storable instance. Though Foreign.Storable isn't in Hugs(Sept06)

-- TODO: depend on my @pointless-fun@ package rather than repeating things here...
-- These two operators make it much easier to read the instance.
-- Hopefully inlining everything will get rid of the eta overhead.
-- <http://matt.immute.net/content/pointless-fun>
(~>) :: (a -> b) -> (d -> c) -> (b -> d) -> a -> c
{-# INLINE (~>) #-}
infixr 2 ~>
f ~> g = (. f) . (g .)

($::) :: a -> (a -> b) -> b
{-# INLINE ($::) #-}
infixl 1 $::
($::) = flip ($)


{-# INLINE logFromLFAssocs #-}
logFromLFAssocs :: [(Int, LogFloat)] -> [(Int, Double)]
#if __GLASGOW_HASKELL__ >= 710
logFromLFAssocs = coerce
#else
logFromLFAssocs = unsafeCoerce
#endif

-- HACK: can't coerce, cf: <https://ghc.haskell.org/trac/ghc/ticket/9220>
{-# INLINE logFromLFUArray #-}
logFromLFUArray :: UArray a LogFloat -> UArray a Double
logFromLFUArray = unsafeCoerce

-- Named unsafe because it could allow injecting NaN if misused
-- HACK: can't coerce, cf: <https://ghc.haskell.org/trac/ghc/ticket/9220>
{-# INLINE unsafeLogToLFUArray #-}
unsafeLogToLFUArray :: UArray a Double -> UArray a LogFloat
unsafeLogToLFUArray = unsafeCoerce

-- Named unsafe because it could allow injecting NaN if misused
{-# INLINE unsafeLogToLFFunc #-}
unsafeLogToLFFunc :: (LogFloat -> a -> LogFloat) -> (Double -> a -> Double)
unsafeLogToLFFunc = ($:: unsafeLogToLogFloat ~> id ~> logFromLogFloat)

-- | Remove the extraneous 'isNaN' test of 'logToLogFloat', when
-- we know it's safe.
{-# INLINE unsafeLogToLogFloat #-}
unsafeLogToLogFloat :: Double -> LogFloat
unsafeLogToLogFloat = LogFloat


instance IArray UArray LogFloat where
    {-# INLINE bounds #-}
    bounds = bounds . logFromLFUArray

-- Apparently this method was added in base-2.0/GHC-6.6 but Hugs
-- (Sept 2006) doesn't have it. Not sure about NHC's base
#if (!(defined(__HUGS__))) || (__HUGS__ > 200609)
    {-# INLINE numElements #-}
    numElements = numElements . logFromLFUArray
#endif

    {-# INLINE unsafeArray #-}
    unsafeArray = unsafeArray $:: id ~> logFromLFAssocs ~> unsafeLogToLFUArray

    {-# INLINE unsafeAt #-}
    unsafeAt = unsafeAt $:: logFromLFUArray ~> id ~> unsafeLogToLogFloat

    {-# INLINE unsafeReplace #-}
    unsafeReplace = unsafeReplace
        $:: logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray

    {-# INLINE unsafeAccum #-}
    unsafeAccum = unsafeAccum
        $:: unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray

    {-# INLINE unsafeAccumArray #-}
    unsafeAccumArray = unsafeAccumArray
        $:: unsafeLogToLFFunc ~> logFromLogFloat ~> id ~> id ~> unsafeLogToLFUArray
#endif

-- TODO: the Nothing branch should never be reachable. Once we get
-- a test suite up and going to /verify/ the never-NaN invariant,
-- we should be able to eliminate the branch and the isNaN checks.
instance PartialOrd LogFloat where
    cmp :: LogFloat -> LogFloat -> Maybe Ordering
cmp (LogFloat Double
x) (LogFloat Double
y)
        | forall a. Transfinite a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| forall a. Transfinite a => a -> Bool
isNaN Double
y = forall a. Maybe a
Nothing
        | Bool
otherwise          = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$! Double
x forall a. Ord a => a -> a -> Ordering
`compare` Double
y

instance Read LogFloat where
    readsPrec :: Int -> ReadS LogFloat
readsPrec Int
p String
s =
        [(Double -> LogFloat
LogFloat (forall a. (Floating a, Transfinite a) => a -> a
log Double
x), String
r) | (Double
x, String
r) <- forall a. Read a => Int -> ReadS a
readsPrec Int
p String
s, Bool -> Bool
not (forall a. Transfinite a => a -> Bool
isNaN Double
x), Double
x forall a. Ord a => a -> a -> Bool
>= Double
0]

----------------------------------------------------------------
-- | Reduce the number of constant string literals we need to store.
errorOutOfRange :: String -> a
{-# NOINLINE errorOutOfRange #-}
errorOutOfRange :: forall a. String -> a
errorOutOfRange String
fun =
    forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$! String
"Data.Number.LogFloat."forall a. [a] -> [a] -> [a]
++String
funforall a. [a] -> [a] -> [a]
++ String
": argument out of range"

-- Both guards are redundant due to the subsequent call to
-- 'Data.Number.Transfinite.log' at all use sites. However we use
-- this function to give local error messages. Perhaps we should
-- catch the exception and throw the new message instead? Portability?

guardNonNegative :: String -> Double -> Double
guardNonNegative :: String -> Double -> Double
guardNonNegative String
fun Double
x
    | forall a. Transfinite a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double
x forall a. Ord a => a -> a -> Bool
< Double
0 = forall a. String -> a
errorOutOfRange String
fun
    | Bool
otherwise        = Double
x

guardIsANumber :: String -> Double -> Double
guardIsANumber :: String -> Double -> Double
guardIsANumber String
fun Double
x
    | forall a. Transfinite a => a -> Bool
isNaN Double
x   = forall a. String -> a
errorOutOfRange String
fun
    | Bool
otherwise = Double
x

----------------------------------------------------------------
-- | 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:
--
-- > logFloat p == logToLogFloat (log p)
--
-- If @p@ is NaN or negative, then the two sides differ only in
-- which error is thrown.
logFloat :: Double -> LogFloat
{-# INLINE [0] logFloat #-}
-- TODO: should we use NOINLINE or [~0] to avoid the possibility of code bloat?
logFloat :: Double -> LogFloat
logFloat = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (Floating a, Transfinite a) => a -> a
log forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> Double -> Double
guardNonNegative String
"logFloat"


-- | Constructor which assumes the argument is already in the
-- log-domain. Throws errors on @notANumber@ inputs.
logToLogFloat :: Double -> LogFloat
logToLogFloat :: Double -> LogFloat
logToLogFloat = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> Double -> Double
guardIsANumber String
"logToLogFloat"


-- | Semantically convert our log-domain value back into the
-- normal-domain. Beware of overflow\/underflow. The following
-- equivalence holds (without qualification):
--
-- > fromLogFloat == exp . logFromLogFloat
--
fromLogFloat :: LogFloat -> Double
{-# INLINE [0] fromLogFloat #-}
-- TODO: should we use NOINLINE or [~0] to avoid the possibility of code bloat?
fromLogFloat :: LogFloat -> Double
fromLogFloat (LogFloat Double
x) = forall a. Floating a => a -> a
exp Double
x


-- | Return the log-domain value itself without conversion.
logFromLogFloat :: LogFloat -> Double
logFromLogFloat :: LogFloat -> Double
logFromLogFloat (LogFloat Double
x) = Double
x


-- These are our module-specific versions of "log\/exp" and "exp\/log";
-- They do the same things but also have a @LogFloat@ 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/fromLogFloat"       forall x. log (fromLogFloat x) = logFromLogFloat x
"logFloat/fromLogFloat"  forall x. logFloat (fromLogFloat x) = x

-- Into log-domain and back out
"fromLogFloat/logFloat"  forall x. fromLogFloat (logFloat x) = x
    #-}

----------------------------------------------------------------
-- 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 Show LogFloat where
    showsPrec :: Int -> LogFloat -> ShowS
showsPrec Int
p (LogFloat Double
x) =
        let y :: Double
y = forall a. Floating a => a -> a
exp Double
x in Double
y seq :: forall a b. a -> b -> b
`seq`
        Bool -> ShowS -> ShowS
showParen (Int
p forall a. Ord a => a -> a -> Bool
> Int
9)
            ( String -> ShowS
showString String
"logFloat "
            forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 Double
y
            )


----------------------------------------------------------------
-- | A curried function for converting arbitrary pairs into ordered
-- pairs. The continuation recieves the minimum first and the maximum
-- second.
--
-- This combinator is primarily intended to reduce repetition in
-- the source code; but hopefully it should also help reduce bloat
-- in the compiled code, by sharing the continuation and just
-- swapping the variables in place. Of course, if the continuation
-- is very small, then requiring a join point after the conditional
-- swap may end up being more expensive than simply duplicating the
-- continuation. Also, given as we're inlining it, I'm not sure
-- whether GHC will decide to keep the sharing we introduced or
-- whether it'll end up duplicating the continuation into the two
-- call sites.
ordered :: Ord a => a -> a -> (a -> a -> b) -> b
ordered :: forall a b. Ord a => a -> a -> (a -> a -> b) -> b
ordered a
x a
y a -> a -> b
k
    | a
x forall a. Ord a => a -> a -> Bool
<= a
y    = a -> a -> b
k a
x a
y
    | Bool
otherwise = a -> a -> b
k a
y a
x
    -- N.B., the implementation of @(>=)@ in Hugs (Sept2006) will
    -- always returns True if either argument isNaN. This does not
    -- constitute a bug for us, since we maintain the invariant that
    -- values wrapped by 'LogFloat' are not NaN.
{-# INLINE ordered #-}


-- TODO: Do we need to add explicit INLINE pragmas here? Or will
-- GHC automatically see that they're small enough to want inlining?
--
-- These all work without causing underflow. However, do note that
-- they tend to induce more of the floating-point fuzz than using
-- regular floating numbers because @exp . log@ doesn't really equal
-- @id@. In any case, our main aim is for preventing underflow when
-- multiplying many small numbers (and preventing overflow for
-- multiplying many large numbers) so we're not too worried about
-- +\/- 4e-16.
instance Num LogFloat where
    * :: LogFloat -> LogFloat -> LogFloat
(*) (LogFloat Double
x) (LogFloat Double
y)
        | forall a. Transfinite a => a -> Bool
isInfinite Double
x Bool -> Bool -> Bool
&& forall a. Transfinite a => a -> Bool
isInfinite Double
y Bool -> Bool -> Bool
&& Double
x forall a. Eq a => a -> a -> Bool
== forall a. Num a => a -> a
negate Double
y =
            Double -> LogFloat
LogFloat forall a. Transfinite a => a
negativeInfinity -- @0 * infinity == 0@
        | Bool
otherwise =
            -- This includes the @0 * 0 == 0@ and @infty * infty == infty@
            -- cases, since @(+)@ treats them appropriately.
            Double -> LogFloat
LogFloat (Double
x forall a. Num a => a -> a -> a
+ Double
y)

    + :: LogFloat -> LogFloat -> LogFloat
(+) (LogFloat Double
x) (LogFloat Double
y)
        | forall a. Transfinite a => a -> Bool
isInfinite Double
x Bool -> Bool -> Bool
&& forall a. Transfinite a => a -> Bool
isInfinite Double
y Bool -> Bool -> Bool
&& Double
x forall a. Eq a => a -> a -> Bool
== Double
y =
            Double -> LogFloat
LogFloat Double
x -- @0 + 0 == 0@ and @infty + infty == infty@
        | Bool
otherwise =
            -- This includes the @0 + infinity == infinity@ case,
            -- since 'log1pexp' (and 'ordered') treats them appropriately.
            forall a b. Ord a => a -> a -> (a -> a -> b) -> b
ordered Double
x Double
y forall a b. (a -> b) -> a -> b
$ \Double
n Double
m ->
            Double -> LogFloat
LogFloat (Double
m forall a. Num a => a -> a -> a
+ Double -> Double
log1pexp (Double
n forall a. Num a => a -> a -> a
- Double
m))

    -- TODO: give a better error message in the (infinity,infinity) case.
    -- TODO: does 'log1mexp' handle the (+infty,-infty) cases correctly?
    (-) (LogFloat Double
x) (LogFloat Double
y)
        | Double
x forall a. Eq a => a -> a -> Bool
== forall a. Transfinite a => a
negativeInfinity Bool -> Bool -> Bool
&& Double
y forall a. Eq a => a -> a -> Bool
== forall a. Transfinite a => a
negativeInfinity =
            Double -> LogFloat
LogFloat forall a. Transfinite a => a
negativeInfinity -- @0 - 0 == 0@
        | Bool
otherwise =
            forall a b. Ord a => a -> a -> (a -> a -> b) -> b
ordered Double
x Double
y forall a b. (a -> b) -> a -> b
$ \Double
n Double
m ->
            Double -> LogFloat
LogFloat (String -> Double -> Double
guardIsANumber String
"(-)" (Double
m forall a. Num a => a -> a -> a
+ Double -> Double
log1mexp (Double
n forall a. Num a => a -> a -> a
- Double
m)))

    signum :: LogFloat -> LogFloat
signum (LogFloat Double
x)
        | Double
x forall a. Eq a => a -> a -> Bool
== forall a. Transfinite a => a
negativeInfinity = LogFloat
0
        | Double
x forall a. Ord a => a -> a -> Bool
>  forall a. Transfinite a => a
negativeInfinity = LogFloat
1
        | Bool
otherwise             = forall a. String -> a
errorOutOfRange String
"signum"
        -- The extra guard protects against NaN, in case someone
        -- broke the invariant. That shouldn't be possible and
        -- so noone else bothers to check, but we check here just
        -- in case.
        -- TODO: wouldn't @not (isNaN x)@ be a better guard to use?

    negate :: LogFloat -> LogFloat
negate LogFloat
_    = forall a. String -> a
errorOutOfRange String
"negate"
    abs :: LogFloat -> LogFloat
abs         = forall a. a -> a
id
    fromInteger :: Integer -> LogFloat
fromInteger = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (Floating a, Transfinite a) => a -> a
log forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> Double -> Double
guardNonNegative String
"fromInteger" forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => Integer -> a
fromInteger


instance Fractional LogFloat where
    -- @n / 0 == infinity@ is handled seamlessly for us. We must catch
    -- @0 / 0@ and @infinity / infinity@ NaNs, and handle @0 / infinity@.
    / :: LogFloat -> LogFloat -> LogFloat
(/) (LogFloat Double
x) (LogFloat Double
y)
        | forall a. Transfinite a => a -> Bool
isInfinite Double
x Bool -> Bool -> Bool
&& forall a. Transfinite a => a -> Bool
isInfinite Double
y Bool -> Bool -> Bool
&& Double
x forall a. Eq a => a -> a -> Bool
== Double
y = forall a. String -> a
errorOutOfRange String
"(/)"
        | Double
x forall a. Eq a => a -> a -> Bool
== forall a. Transfinite a => a
negativeInfinity = Double -> LogFloat
LogFloat forall a. Transfinite a => a
negativeInfinity -- @0/infinity == 0@
        | Bool
otherwise             = Double -> LogFloat
LogFloat (Double
x forall a. Num a => a -> a -> a
- Double
y)

    fromRational :: Rational -> LogFloat
fromRational = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (Floating a, Transfinite a) => a -> a
log
                 forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> Double -> Double
guardNonNegative String
"fromRational" forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Fractional a => Rational -> a
fromRational


-- Just for fun. The more coercion functions the better. Though
-- Rationals are very buggy when it comes to transfinite values
instance Real LogFloat where
    toRational :: LogFloat -> Rational
toRational (LogFloat Double
x)
        | forall a. Transfinite a => a -> Bool
isInfinite Double
ex Bool -> Bool -> Bool
|| forall a. Transfinite a => a -> Bool
isNaN Double
ex = forall a. String -> a
errorOutOfRange String
"toRational"
        | Bool
otherwise                 = forall a. Real a => a -> Rational
toRational Double
ex
        where
        ex :: Double
ex = forall a. Floating a => a -> a
exp Double
x


----------------------------------------------------------------
-- | /O(1)/. Compute powers in the log-domain; that is, the following
-- equivalence holds (modulo underflow and all that):
--
-- > logFloat (p ** m) == logFloat p `pow` m
--
-- /Since: 0.13/
pow :: LogFloat -> Double -> LogFloat
{-# INLINE pow #-}
infixr 8 `pow`
pow :: LogFloat -> Double -> LogFloat
pow (LogFloat Double
x) Double
m
    | forall a. Transfinite a => a -> Bool
isNaN Double
mx  = Double -> LogFloat
LogFloat Double
0
    | Bool
otherwise = Double -> LogFloat
LogFloat Double
mx
    where
    -- N.B., will be NaN when @x == negativeInfinity && m == 0@
    -- (which is true when m is -0 as well as +0). We check for NaN
    -- after multiplying, rather than checking this precondition
    -- before multiplying, in an attempt to simplify/optimize the
    -- generated code.
    -- TODO: benchmark.
    mx :: Double
mx = Double
m forall a. Num a => a -> a -> a
* Double
x


-- N.B., the default implementation of (**) for Complex is wrong.
-- It can be fixed by using the definition:
-- > x ** y = if x == 0 then 0 else exp (log x * y)
-- cf., <https://ghc.haskell.org/trac/ghc/ticket/8539>
-- TODO: Is this relevant to us at all?


-- TODO: check out ekmett's compensated library.


-- Some good test cases:
-- for @logsumexp == log . sum . 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 'LogFloat's, being
-- careful to avoid underflow issues. That is, the following
-- equivalence holds (modulo underflow and all that):
--
-- > logFloat . sum == sum . map logFloat
--
-- /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.
--
-- /Since: 0.13/
sum :: [LogFloat] -> LogFloat
sum :: [LogFloat] -> LogFloat
sum = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Double
logSumExp forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap LogFloat -> Double
logFromLogFloat


-- | /O(n)/. Compute the product of a finite list of 'LogFloat's,
-- being careful to avoid numerical error due to loss of precision.
-- That is, the following equivalence holds (modulo underflow and
-- all that):
--
-- > logFloat . product == product . map logFloat
--
-- /Since: 0.13/
product :: [LogFloat] -> LogFloat
product :: [LogFloat] -> LogFloat
product = Double -> LogFloat
LogFloat forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Double
kahanSum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap LogFloat -> Double
logFromLogFloat

----------------------------------------------------------------
----------------------------------------------------------- fin.