logfloat-0.13.4: src/Data/Number/LogFloat.hs
-- 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 #-}
----------------------------------------------------------------
-- ~ 2017.06.18
-- |
-- Module : Data.Number.LogFloat
-- Copyright : Copyright (c) 2007--2017 wren gayle romano
-- License : BSD3
-- Maintainer : wren@community.haskell.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.List (foldl')
import Data.Number.Transfinite
import Data.Number.PartialOrd
-- 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
( Eq
, 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
, 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 = unsafeCoerce (bounds :: UArray i Double -> (i, i))
{-# INLINE numElements #-}
numElements :: forall i. Ix i => UArray i LogFloat -> Int
numElements = unsafeCoerce (numElements :: UArray i Double -> Int)
{-# INLINE unsafeArray #-}
unsafeArray :: forall i. Ix i => (i,i) -> [(Int,LogFloat)] -> UArray i LogFloat
unsafeArray = unsafeCoerce (unsafeArray :: (i,i) -> [(Int,Double)] -> UArray i Double)
{-# INLINE unsafeAt #-}
unsafeAt :: forall i. Ix i => UArray i LogFloat -> Int -> LogFloat
unsafeAt = unsafeCoerce (unsafeAt :: UArray i Double -> Int -> Double)
{-# INLINE unsafeReplace #-}
unsafeReplace :: forall i. Ix i => UArray i LogFloat -> [(Int,LogFloat)] -> UArray i LogFloat
unsafeReplace = unsafeCoerce (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 = unsafeCoerce (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 = unsafeCoerce (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)
-- 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 x) (LogFloat y)
| isNaN x || isNaN y = Nothing
| otherwise = Just $! x `compare` y
instance Read LogFloat where
readsPrec p s =
[(LogFloat (log x), r) | (x, r) <- readsPrec p s, not (isNaN x), x >= 0]
----------------------------------------------------------------
-- | Reduce the number of constant string literals we need to store.
errorOutOfRange :: String -> a
{-# NOINLINE errorOutOfRange #-}
errorOutOfRange fun =
error $! "Data.Number.LogFloat."++fun++ ": 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 fun x
| isNaN x || x < 0 = errorOutOfRange fun
| otherwise = x
guardIsANumber :: String -> Double -> Double
guardIsANumber fun x
| isNaN x = errorOutOfRange fun
| otherwise = 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 = LogFloat . log . guardNonNegative "logFloat"
-- | Constructor which assumes the argument is already in the
-- log-domain. Throws errors on @notANumber@ inputs.
logToLogFloat :: Double -> LogFloat
logToLogFloat = LogFloat . guardIsANumber "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 x) = exp x
-- | Return the log-domain value itself without conversion.
logFromLogFloat :: LogFloat -> Double
logFromLogFloat (LogFloat x) = 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 p (LogFloat x) =
let y = exp x in y `seq`
showParen (p > 9)
( showString "logFloat "
. showsPrec 11 y
)
----------------------------------------------------------------
-- Technically these should use 'Foreign.C.CDouble' however there's
-- no isomorphism provided to normal 'Double'. The former is
-- documented as being a newtype of the later, and so this should
-- be safe.
#ifdef __USE_FFI__
#define LOG1P_WHICH_VERSION FFI version.
#else
#define LOG1P_WHICH_VERSION naive version! \
Contact the maintainer with any FFI difficulties.
#endif
-- | Definition: @log1p == log . (1+)@. Standard C libraries provide
-- a special definition for 'log1p' which is more accurate than
-- doing the naive thing, especially for very small arguments. For
-- example, the naive version underflows around @2 ** -53@, whereas
-- the specialized version underflows around @2 ** -1074@. This
-- function is used by ('+') and ('-') on @LogFloat@.
--
-- N.B. The @statistics:Statistics.Math@ module provides a pure
-- Haskell implementation of @log1p@ for those who are interested.
-- We do not copy it here because it relies on the @vector@ package
-- which is non-portable. If there is sufficient interest, a portable
-- variant of that implementation could be made. Contact the
-- maintainer if the FFI and naive implementations are insufficient
-- for your needs.
--
-- /This installation was compiled to use the LOG1P_WHICH_VERSION/
#ifdef __USE_FFI__
foreign import ccall unsafe "math.h log1p"
log1p :: Double -> Double
#else
-- See statistics:Statistics.Math for a more accurate Haskell
-- implementation.
log1p :: Double -> Double
{-# INLINE [0] log1p #-}
log1p x = log (1 + x)
#endif
-- | Definition: @expm1 == subtract 1 . exp@. Standard C libraries
-- provide a special definition for 'expm1' which is more accurate
-- than doing the naive thing, especially for very small arguments.
-- This function isn't needed internally, but is provided for
-- symmetry with 'log1p'.
--
-- /This installation was compiled to use the LOG1P_WHICH_VERSION/
#ifdef __USE_FFI__
foreign import ccall unsafe "math.h expm1"
expm1 :: Double -> Double
#else
expm1 :: Double -> Double
{-# INLINE [0] expm1 #-}
expm1 x = exp x - 1
#endif
-- CPP guarded because they won't fire if we're using the FFI versions
#if !defined(__USE_FFI__)
{-# 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
#-}
#endif
----------------------------------------------------------------
-- 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
-- N.B. In Hugs (Sept2006) the (>=) always returns True if
-- either isNaN. This does not constitute a bug since we
-- maintain the invariant that values wrapped by 'LogFloat'
-- are not NaN.
(*) (LogFloat x) (LogFloat y)
| isInfinite x
&& isInfinite y
&& x == negate y = LogFloat negativeInfinity -- @0*infinity == 0@
| otherwise = LogFloat (x+y)
(+) (LogFloat x) (LogFloat y)
| x == y
&& isInfinite x
&& isInfinite y = LogFloat x -- @0+0 == 0@, @infinity+infinity == infinity@
| x >= y = LogFloat (x + log1p (exp (y - x)))
| otherwise = LogFloat (y + log1p (exp (x - y)))
(-) (LogFloat x) (LogFloat y)
| x == negativeInfinity
&& y == negativeInfinity = LogFloat negativeInfinity -- @0-0 == 0@
| otherwise =
-- BUG: Will throw error if x < y
-- TODO: flip @x@ and @y@ when @y > x@.
-- Also, will throw error if (x,y) is (infinity,infinity)
LogFloat (guardIsANumber "(-)" (x + log1p (negate (exp (y - x)))))
signum (LogFloat x)
| x == negativeInfinity = 0
| x > negativeInfinity = 1
| otherwise = errorOutOfRange "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.
negate _ = errorOutOfRange "negate"
abs = id
fromInteger = LogFloat . log
. guardNonNegative "fromInteger" . 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 x) (LogFloat y)
| x == y
&& isInfinite x
&& isInfinite y = errorOutOfRange "(/)"
| x == negativeInfinity = LogFloat negativeInfinity -- @0/infinity == 0@
| otherwise = LogFloat (x-y)
fromRational = LogFloat . log
. guardNonNegative "fromRational" . 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 x)
| isInfinite ex || isNaN ex = errorOutOfRange "toRational"
| otherwise = toRational ex
where
ex = exp 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 x) m
| isNaN mx = LogFloat 0
| otherwise = LogFloat 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 = m * 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 xs = LogFloat (theMax + log theSum)
where
LogFloat theMax = maximum xs
-- compute @\log \sum_{x \in xs} \exp(x - theMax)@
theSum = foldl' (\ acc (LogFloat x) -> acc + exp (x - theMax)) 0 xs
-- TODO: expose a single-pass version for the special case where
-- the first element of the list is (promised to be) the maximum
-- element?
-- | /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 = kahan 0 0
where
kahan t c _ | t `seq` c `seq` False = undefined
kahan t _ [] = LogFloat t
kahan t c (LogFloat x : xs)
-- Avoid NaN when there's a negInfty in the list. N.B.,
-- this causes zero to annihilate infinity.
| x == negativeInfinity = LogFloat negativeInfinity
| otherwise =
-- Beware this getting incorrectly optimized away by
-- constant folding!
let y = x - c
t' = t + y
c' = (t' - t) - y
in kahan t' c' xs
-- 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)
-}
----------------------------------------------------------------
----------------------------------------------------------- fin.