logfloat 0.13.4 → 0.14.0
raw patch · 9 files changed
+585/−265 lines, 9 files
Files
- CHANGELOG +8/−0
- README.md +6/−31
- logfloat.cabal +21/−11
- src/Data/Number/LogFloat.hs +89/−184
- src/Data/Number/LogFloat/Raw.hs +403/−0
- src/Data/Number/PartialOrd.hs +12/−12
- src/Data/Number/RealToFrac.hs +6/−6
- src/Data/Number/Transfinite.hs +36/−17
- src/Hugs/RealFloat.hs +4/−4
CHANGELOG view
@@ -1,6 +1,14 @@+0.14.0 (2021-xx-xx):+ - Separating out "Data.Number.LogFloat.Raw"+ - Added `Tested-With: GHC == 9.2.4, 9.4.4, 9.6.1` (didn't actually need to+ nudge the upper bound on 'base', because it's already lenient)+ - Silencing warnings about rewrite rules, as a hack for working around:+ <https://gitlab.haskell.org/ghc/ghc/-/issues/10595> 0.13.4 (2017-06-18): - Added Read LogFloat instance. (h/t Rob Zinkov) - Corrected a bug where `product` returns NaN on `logFloat 0`.+0.13.3.3 (2015-10-02):+ - ??? 0.13.3.2 (2015-08-06): - Fixed the buggy Show LogFloat instance 0.13.3.1 (2015-05-30):
README.md view
@@ -1,9 +1,8 @@ logfloat ======== [](https://hackage.haskell.org/package/logfloat)-[](http://packdeps.haskellers.com/specific?package=logfloat)-[](https://travis-ci.org/wrengr/logfloat)-[](https://circleci.com/gh/wrengr/logfloat)+[](https://github.com/wrengr/logfloat/actions?query=workflow%3Aci)+[](http://packdeps.haskellers.com/specific?package=logfloat) This package provides a type for storing numbers in the log-domain, primarily useful for preventing underflow when multiplying many@@ -19,38 +18,17 @@ point between backwards compatability and adding new features (see below). -Note that the GitHub repository is just a clone of [the Darcs-repo](http://code.haskell.org/~wren/logfloat/). I'm testing out-whether to switch things over to GitHub in order to use TravisCI,-and an official ticket tracker, etc.-- ## Install In general, this is a simple package and should be easy to install. The specifics are a bit murky however, since we use CPP and the FFI and want to retain portability, and due to the rapid evolution of-Cabal and other development tools. You should be able to use one-of the following standard installation methods:+Cabal and other development tools. However, you should be able to+use the standard: - -- With cabal-install and without the source: $> cabal install logfloat - -- With cabal-install and with the source already:- $> cd logfloat- $> cabal install - -- Without cabal-install, but with the source already:- $> cd logfloat- $> runhaskell Setup.hs configure --user- $> runhaskell Setup.hs build- $> runhaskell Setup.hs haddock --hyperlink-source- $> runhaskell Setup.hs copy- $> runhaskell Setup.hs register--The Haddock step is optional.-- ### FFI Problems The logfloat package uses the FFI to access functions in libm to@@ -252,11 +230,8 @@ ## Links -* [Website](http://cl.indiana.edu/~wren/)+* [Website](http://wrengr.org/) * [Blog](http://winterkoninkje.dreamwidth.org/) * [Twitter](https://twitter.com/wrengr) * [Hackage](http://hackage.haskell.org/package/logfloat)-* [Darcs](http://code.haskell.org/~wren/logfloat)-* [GitHub (clone)](https://github.com/wrengr/logfloat)-* [Haddock (Darcs version)- ](http://code.haskell.org/~wren/logfloat/dist/doc/html/logfloat)+* [GitHub](https://github.com/wrengr/logfloat)
logfloat.cabal view
@@ -1,19 +1,25 @@+Cabal-Version: 2.2+-- Cabal >=2.2 is required for:+-- <https://cabal.readthedocs.io/en/latest/cabal-package.html#common-stanzas>+-- Since 2.1, the Cabal-Version must be the absolutely first thing+-- in the file, even before comments. Also, no longer uses ">=".+-- <https://github.com/haskell/cabal/issues/4899>+ ------------------------------------------------------------------- wren gayle romano <wren@community.haskell.org> ~ 2021.10.16+-- wren gayle romano <wren@cpan.org> ~ 2023.03.19 ---------------------------------------------------------------- --- Cabal >=1.10 is required by Hackage.-Cabal-Version: >= 1.10-Build-Type: Simple- Name: logfloat-Version: 0.13.4+Version: 0.14.0+Build-Type: Simple Stability: experimental-Homepage: http://wrengr.org+Homepage: https://wrengr.org/software/hackage.html+Bug-Reports: https://github.com/wrengr/logfloat/issues Author: wren gayle romano Maintainer: wren@cpan.org-Copyright: Copyright (c) 2007--2021 wren gayle romano-License: BSD3+Copyright: 2007–2023 wren romano+-- Cabal-2.2 requires us to say "BSD-3-Clause" not "BSD3"+License: BSD-3-Clause License-File: LICENSE Category: Data, Math, Natural Language Processing, Statistics@@ -37,12 +43,15 @@ GHC ==8.6.5, GHC ==8.8.4, GHC ==8.10.3,- GHC ==9.0.1+ GHC ==9.0.1,+ GHC ==9.2.4,+ GHC ==9.4.4,+ GHC ==9.6.1 ---------------------------------------------------------------- Source-Repository head Type: git- Location: git://github.com/wrengr/logfloat.git+ Location: https://github.com/wrengr/logfloat.git ---------------------------------------------------------------- Flag useFFI@@ -57,6 +66,7 @@ Default-Language: Haskell2010 Hs-Source-Dirs: src Exposed-Modules: Data.Number.LogFloat+ , Data.Number.LogFloat.Raw , Data.Number.RealToFrac , Data.Number.Transfinite , Data.Number.PartialOrd
src/Data/Number/LogFloat.hs view
@@ -14,12 +14,12 @@ {-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-} ------------------------------------------------------------------- ~ 2017.06.18+-- ~ 2021.10.17 -- | -- Module : Data.Number.LogFloat--- Copyright : Copyright (c) 2007--2017 wren gayle romano+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable (with CPP, FFI) --@@ -62,10 +62,10 @@ ) where import Prelude hiding (log, sum, product, isInfinite, isNaN)-import Data.List (foldl') import Data.Number.Transfinite import Data.Number.PartialOrd+import Data.Number.LogFloat.Raw -- GHC can derive (IArray UArray LogFloat), but Hugs needs to coerce@@ -103,7 +103,7 @@ -- 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+-- 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:@@ -179,28 +179,44 @@ 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)+ 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)+ 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)+ 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)+ 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) +-- 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>@@ -258,28 +274,26 @@ #endif {-# INLINE unsafeArray #-}- unsafeArray =- unsafeArray $:: id ~> logFromLFAssocs ~> unsafeLogToLFUArray+ unsafeArray = unsafeArray $:: id ~> logFromLFAssocs ~> unsafeLogToLFUArray {-# INLINE unsafeAt #-}- unsafeAt =- unsafeAt $:: logFromLFUArray ~> id ~> unsafeLogToLogFloat+ unsafeAt = unsafeAt $:: logFromLFUArray ~> id ~> unsafeLogToLogFloat {-# INLINE unsafeReplace #-}- unsafeReplace =- unsafeReplace $:: logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray+ unsafeReplace = unsafeReplace+ $:: logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray {-# INLINE unsafeAccum #-}- unsafeAccum =- unsafeAccum $:: unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray+ unsafeAccum = unsafeAccum+ $:: unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray {-# INLINE unsafeAccumArray #-}- unsafeAccumArray =- unsafeAccumArray $:: unsafeLogToLFFunc ~> logFromLogFloat ~> id ~> id ~> unsafeLogToLFUArray+ 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,+-- 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)@@ -385,77 +399,34 @@ ------------------------------------------------------------------- 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.+-- | A curried function for converting arbitrary pairs into ordered+-- pairs. The continuation recieves the minimum first and the maximum+-- second. ----- /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+-- 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 x y k+ | x <= y = k x y+ | otherwise = k y 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 #-} --- | 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'.+-- TODO: Do we need to add explicit INLINE pragmas here? Or will+-- GHC automatically see that they're small enough to want inlining? ----- /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@@ -463,34 +434,32 @@ -- 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)+ | isInfinite x && isInfinite y && x == negate y =+ LogFloat negativeInfinity -- @0 * infinity == 0@+ | otherwise =+ -- This includes the @0 * 0 == 0@ and @infty * infty == infty@+ -- cases, since @(+)@ treats them appropriately.+ 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)))+ | isInfinite x && isInfinite y && x == y =+ LogFloat x -- @0 + 0 == 0@ and @infty + infty == infty@+ | otherwise =+ -- This includes the @0 + infinity == infinity@ case,+ -- since 'log1pexp' (and 'ordered') treats them appropriately.+ ordered x y $ \n m ->+ LogFloat (m + log1pexp (n - m)) + -- TODO: give a better error message in the (infinity,infinity) case.+ -- TODO: does 'log1mexp' handle the (+infty,-infty) cases correctly? (-) (LogFloat x) (LogFloat y)- | x == negativeInfinity- && y == negativeInfinity = LogFloat negativeInfinity -- @0-0 == 0@+ | 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)))))+ ordered x y $ \n m ->+ LogFloat (guardIsANumber "(-)" (m + log1mexp (n - m))) signum (LogFloat x) | x == negativeInfinity = 0@@ -500,24 +469,20 @@ -- 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 _ = errorOutOfRange "negate"- abs = id-- fromInteger = LogFloat . log- . guardNonNegative "fromInteger" . fromInteger+ 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.+ -- @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 "(/)"+ | isInfinite x && isInfinite y && x == y = errorOutOfRange "(/)" | x == negativeInfinity = LogFloat negativeInfinity -- @0/infinity == 0@- | otherwise = LogFloat (x-y)+ | otherwise = LogFloat (x - y) fromRational = LogFloat . log . guardNonNegative "fromRational" . fromRational@@ -585,17 +550,7 @@ -- -- /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?-+sum = LogFloat . logSumExp . fmap logFromLogFloat -- | /O(n)/. Compute the product of a finite list of 'LogFloat's,@@ -607,57 +562,7 @@ -- -- /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)--}+product = LogFloat . kahanSum . fmap logFromLogFloat ---------------------------------------------------------------- ----------------------------------------------------------- fin.
+ src/Data/Number/LogFloat/Raw.hs view
@@ -0,0 +1,403 @@+{-# LANGUAGE CPP, ForeignFunctionInterface, BangPatterns #-}+{-# OPTIONS_GHC -Wall -fwarn-tabs #-}+{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-}++----------------------------------------------------------------+-- ~ 2022.03.14+-- |+-- Module : Data.Number.LogFloat.Raw+-- Copyright : 2007--2022 wren romano+-- License : BSD-3-Clause+-- Maintainer : wren@cpan.org+-- Stability : provisional+-- Portability : portable (with CPP, FFI)+--+-- This module provides implementations for computing various+-- logarithmic and exponential functions without losing precision+-- (as the naive implementations do). These are the \"raw\"+-- implementations; i.e., sans newtypes and other conveniences.+-- Since the lack of newtypes means we can't rely on types to clarify+-- things, we use the traditional baroque names for things. The+-- design considerations behind (most of) these implementations are+-- documented at:+-- <https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf>+--+-- In base-4.9.0.0 GHC added some of these to the 'Floating' class+-- exported from "Numeric". Alas, they provide default definitions+-- using the naive implementations, so one can't really rely on the+-- 'Floating' methods being precision preserving. Overall, the+-- specific instance for 'Double' looks fine (though they use+-- different cutoffs for 'log1pexp' for some reason); but it's easy+-- enough to reimplement here, to make absolutely sure we're getting+-- the right thing.+--+-- @since: 0.14.0+----------------------------------------------------------------+module Data.Number.LogFloat.Raw+ (+ -- * Logarithmic\/exponential basics+ expm1+ , log1p+ , log1mexp+ , log1pexp+ -- * Summation+ , logSumExp+ , kahanSum+ -- , neumaierSum+ -- * Softmax+ , logSoftmax+ , softmax+ -- * Sigmoid and related functions+ , sigmoid+ , logit+ , logitExp+ ) where++import Data.List (foldl')+import Data.Number.Transfinite (negativeInfinity)++----------------------------------------------------------------+-- 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+++-- | Compute @log (1 + x)@ without losing precision.+--+-- Standard C libraries provide a special definition for this+-- function, 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@.+--+-- 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__+-- TODO: verify that the Haddock comes out as intended...+foreign import ccall unsafe "math.h log1p"+ log1p+ :: Double -- ^ N.B., only defined on the @[-1,infty]@ interval.+ -> Double+#else+-- See @statistics@:"Statistics.Math" for a more accurate Haskell+-- implementation.+log1p+ :: Double -- ^ N.B., only defined on the @[-1,infty]@ interval.+ -> Double+{-# INLINE [0] log1p #-}+log1p x = log (1 + x)+#endif+++-- | Compute @exp x - 1@ without losing precision.+--+-- Standard C libraries provide a special definition for 'expm1'+-- which is more accurate than doing the naive thing, especially+-- for very small arguments.+--+-- /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.+-- TODO: can we get them to fire if we to the standard thing about+-- naming the FFI version @c_foo@ and then defining a Haskell+-- function @foo = c_foo@?+#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+++-- | Compute @log (1 - exp x)@ without losing precision.+log1mexp+ :: Double -- ^ N.B., only defined on the @[-infty,0]@ interval.+ -> Double+log1mexp x+ | x <= log 2 = (log . negate . expm1) x+ | otherwise = (log1p . negate . exp) x+{-# INLINE log1mexp #-}+++-- | Compute @log (1 + exp x)@ without losing precision. Algebraically+-- this is @0 ⊔ x@, which is the log-domain's analogue of @1 + x@.+log1pexp :: Double -> Double+log1pexp x+ | x <= -37 = exp x+ | x <= 18 = log1p (exp x)+ | x <= 33.3 = x + exp (negate x)+ | otherwise = x+{-# INLINE log1pexp #-}+++-- TODO: bring back 'expm1c' and 'log1pc'+++----------------------------------------------------------------+-- | The logistic function; aka, the inverse of 'logit'.+-- > sigmoid x = 1 / (1 + exp (-x))+-- > sigmoid x = exp x / (exp x + 1)+-- > sigmoid x = (1 + tanh (x/2)) / 2+sigmoid :: Double -> Double+sigmoid x = (1 + tanh (x/2)) / 2+{-# INLINE sigmoid #-}+-- We prefer the 'tanh'-based definition because it's (exactly!)+-- symmetric about zero, whereas the naive version isn't (due to+-- floating-point fuzz).+-- TODO(b/68203642): Properly analyze the accuracy and precision+-- of the 'tanh' version.+++-- | The quantile function; aka, the inverse of 'sigmoid'.+-- > logit x = log (x / (1 - x))+-- > logit x = 2 * atanh (2*x - 1)+logit+ :: Double -- ^ N.B., only defined on the @[0,1]@ interval.+ -> Double+logit x = 2 * atanh (2*x - 1)+{-# INLINE logit #-}+-- TODO(b/68203642): properly analyze the precision of the 'atanh' version.+++-- | A variant of 'logit' for when the argument is already in the+-- log-domain; hence, @logitExp = logit . exp@+logitExp+ :: Double -- ^ N.B., only defined on the @[-infty,0]@ interval.+ -> Double+logitExp x = x - log1mexp x+{-# INLINE logitExp #-}+-- TODO(b/68203642): properly analyze the precision of this+-- implementation with respect to the 'logit' implementation.+++----------------------------------------------------------------+-- TODO: double check that everything inlines away, so this data+-- type doesn't introduce any slowdown.+--+-- | A helper type for 'logSumExp'. As a semigroup, this is isomorphic to:+-- @(WrappedMonoid (Sum Int), Max Double)@; however, we strictify and+-- flatten everything to improve performance.+data LSE = LSE+ {-# UNPACK #-}!Int -- The length, minus one.+ {-# UNPACK #-}!Double -- The maximum.++-- | Compute the length and maximum of a list. This is a semigroup+-- reduction. However we roll it ourselves rather than using the+-- semigroup class: since that would incur an otherwise unnecessary+-- dependency on @base >= 4.9.0.0@.+foldLSE :: Double -> [Double] -> LSE+foldLSE = foldl' step . LSE 0+ where+ step (LSE lm1 m) x = LSE (lm1 + 1) (m `max` x)+++-- 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)/. Log-domain summation, aka: @(log . sum . fmap exp)@.+-- Algebraically this is @⨆ xs@, which is the log-domain equivalent+-- of @∑ xs@.+--+-- /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.+logSumExp :: [Double] -> Double+logSumExp [] = (-1)/0+logSumExp xs0@(x:xs) =+ case foldLSE x xs of+ LSE lm1 m+ | isInfinite m -> m+ | otherwise ->+ -- TODO: push the addition of @lm1@ into the 'kahanSum',+ -- but making sure to add it in only at the very end.+ -- TODO: would using 'neumaierSum' be better? Should+ -- we factor the summation function out as an argument?+ -- TODO: is using 'log1p' here /really/ any better than+ -- just using 'log'?+ -- TODO: does that 'fmap' properly fuse into the+ -- 'kahanSum', or need we inline it ourselves?+ m + log1p (fromIntegral lm1 + kahanSum (fmap (expm1 . subtract m) xs0))++{-+-- TODO(wrengr): Compare precision of the following implementations.+-- We need to make sure to structure it in such a way that the @m@+-- doesn't obliterate the whole purpose of using @exp (x - m)@ in+-- the first place; but supposing we can do that, then it might+-- could help++sumExp = exp . logSumExp++sumExp [] = 0+sumExp xs0@(x:xs) =+ case foldLSE x xs of+ LSE lm1 m+ | isInfinite m -> m+ | otherwise ->+ exp m * kahanSum (fromIntegral lm1 : fmap (expm1 . subtract m) xs0)+-}+++----------------------------------------------------------------+-- | /O(n)/. Log-domain softmax, aka: @(fmap log . softmax)@.+--+-- /N.B./, this requires three passes over the data: two for the+-- 'logSumExp', and a third for the normalization itself. Thus,+-- it is not amenable to list fusion, and hence will use a lot of+-- memory when summing long lists.+logSoftmax :: [Double] -> [Double]+logSoftmax xs = let z = logSumExp xs in z `seq` fmap (subtract z) xs+-- TODO(wrengr): alternatively we could use a variant of 'logSumExp'+-- which doesn't add the maximum back in, and do the final rescaling+-- by subtracting both the maximum and the summation; that is, a more+-- efficient\/straightforward variant of:+-- > logSoftmax xs =+-- > subtract z <$> xs' -- aka @subtract (m + z) <$> xs@+-- > where+-- > m = maximum xs+-- > xs' = subtract m <$> xs+-- > z = logSumExp xs'+-- This works because for any constant @c@, @softmax xs == softmax ((+c)+-- <$> xs)@. Of course, I don't know that doing that would really help+-- precision by much (given the improved performance of using 'logSumExp'+-- in the first place), and saving a single add won't really matter+-- performance-wise. Perhaps if instead of the thing just proposed+-- about avoiding adding the max back in, what if instead we did things+-- exactly as written above: so we subtract off the maximum, but then+-- also do 'logSumExp' such that it subtracts off the maximum of those+-- differences. We could get the top-2 maxima in a single pass without+-- much extra work; but again, unclear whether it'd really help...+++-- | /O(n)/. Normal-domain softmax:+-- > softmax xs = [ exp x / sum [ exp y | y <- xs] | x <- xs ]+--+-- /N.B./, this requires three passes over the data: same as 'logSoftmax'.+softmax :: [Double] -> [Double]+softmax = fmap exp . logSoftmax+{-# INLINE softmax #-}+++----------------------------------------------------------------+-- TODO: double check that everything inlines away, so this data+-- type doesn't introduce any slowdown.+--+-- | A helper type for 'kahanSum'. As a data type, this is really+-- just so we can phrase things as using 'foldl''.+data Kahan = Kahan+ {-# UNPACK #-}!Double -- The total.+ {-# UNPACK #-}!Double -- The error correction.++kahanZero :: Kahan+kahanZero = Kahan 0 0+{-# INLINE kahanZero #-}++-- DONOTSUBMIT: if @x == negativeInfinity@ then our use case demands we return negativeInfinity (so that @0 * infinity == 0@ as desired). But moreover, we really want to short-circuit things to avoid even scanning the rest of the list. To do that, we need to re-inline everything and use recursion directly instead of using 'foldl''.+kahanPlus :: Kahan -> Double -> Kahan+kahanPlus (Kahan t c) x = Kahan t' c'+ where+ -- Beware this getting incorrectly optimized away by constant folding!+ x' = x - c+ t' = t + x'+ c' = (t' - t) - x'+{-# INLINE kahanPlus #-}++fromKahan :: Kahan -> Double+fromKahan (Kahan t _) = t+{-# INLINE fromKahan #-}++-- | /O(n)/. Floating-point summation, via Kahan's algorithm. This+-- is nominally equivalent to 'sum', but greatly mitigates the+-- problem of losing precision.+--+-- /N.B./, this only requires a single pass over the data; but we+-- use a strict left fold for performance, so it's still not amenable+-- to list fusion.+kahanSum :: [Double] -> Double+{-+-- Alas, this implementation loses the optimization below where we+-- avoid NaN and short-circuit to return @LogFloat -infty@ aka 0.+kahanSum = fromKahan . foldl' kahanPlus kahanZero+-}+kahanSum = go kahanZero+ where+ go tc _ | tc `seq` False = undefined+ go tc [] = fromKahan tc+ go tc (x:xs)+ -- Avoid NaN when there's a negInfty in the list. N.B.,+ -- this causes zero to annihilate infinity.+ | x == negativeInfinity = negativeInfinity+ | otherwise = go (kahanPlus tc x) xs+++-- TODO: bring back the 'neumaierSum'+++-- 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.
src/Data/Number/PartialOrd.hs view
@@ -10,15 +10,15 @@ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} ------------------------------------------------------------------- ~ 2015.03.29+-- ~ 2021.10.17 -- | -- Module : Data.Number.PartialOrd--- Copyright : Copyright (c) 2007--2015 wren gayle romano+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : semi-portable (OverlappingInstances,...)--- +-- -- The Prelude's 'Ord' class for dealing with ordered types is often -- onerous to use because it requires 'Eq' as well as a total -- ordering. While such total orderings are common, partial orderings@@ -49,49 +49,49 @@ class PartialOrd a where -- | like 'compare' cmp :: a -> a -> Maybe Ordering- + -- | like ('>') gt :: a -> a -> Maybe Bool gt x y = case x `cmp` y of Just GT -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like ('>=') ge :: a -> a -> Maybe Bool ge x y = case x `cmp` y of Just LT -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('==') eq :: a -> a -> Maybe Bool eq x y = case x `cmp` y of Just EQ -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like ('/=') ne :: a -> a -> Maybe Bool ne x y = case x `cmp` y of Just EQ -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('<=') le :: a -> a -> Maybe Bool le x y = case x `cmp` y of Just GT -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('<') lt :: a -> a -> Maybe Bool lt x y = case x `cmp` y of Just LT -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like 'max'. The default instance returns the left argument -- when they're equal. maxPO :: a -> a -> Maybe a@@ -100,7 +100,7 @@ GT -> Just x EQ -> Just x LT -> Just y- + -- | like 'min'. The default instance returns the left argument -- when they're equal. minPO :: a -> a -> Maybe a
src/Data/Number/RealToFrac.hs view
@@ -9,22 +9,22 @@ #if __GLASGOW_HASKELL__ < 710 {-# LANGUAGE OverlappingInstances #-} #endif- + -- We don't put these in LANGUAGE, because it's CPP guarded for GHC only {-# OPTIONS_GHC -XMagicHash #-} {-# OPTIONS_GHC -Wall -fwarn-tabs #-} ------------------------------------------------------------------- ~ 2013.05.29+-- ~ 2021.10.17 -- | -- Module : Data.Number.RealToFrac--- Copyright : Copyright (c) 2007--2015 wren gayle romano+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : semi-portable (CPP, MPTC, OverlappingInstances)--- +-- -- This module presents a type class for generic conversion between -- numeric types, generalizing @realToFrac@ in order to overcome -- problems with pivoting through 'Rational'@@ -138,7 +138,7 @@ RealToFrac Float Double where {-# INLINE realToFrac #-} realToFrac (F# f) = D# (float2Double# f)- + instance #if __GLASGOW_HASKELL__ >= 710 {-# OVERLAPPING #-}
src/Data/Number/Transfinite.hs view
@@ -1,22 +1,33 @@ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} {-# OPTIONS_GHC -O2 -fenable-rewrite-rules #-}++-- FIXME(2023-03-19): Since recent versions of GHC the rewrite rules+-- in this file generate warnings that they may never fire because+-- the rule "Class op exp" may fire first. Although the warning+-- suggests adding a phase limit, that doesn't actually help because+-- the "Class op exp" rule is a built-in. So it's not actually+-- clear how to silence this warning:+-- <https://gitlab.haskell.org/ghc/ghc/-/issues/10595>+--+-- So for now we just silence the warnings.+{-# OPTIONS_GHC -fno-warn-inline-rule-shadowing #-} ------------------------------------------------------------------- ~ 2015.03.29+-- ~ 2021.10.17 -- | -- Module : Data.Number.Transfinite--- Copyright : Copyright (c) 2007--2015 wren gayle romano+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable--- +-- -- This module presents a type class for numbers which have -- representations for transfinite values. The idea originated from -- the IEEE-754 floating-point special values, used by -- "Data.Number.LogFloat". However not all 'Fractional' types -- necessarily support transfinite values. In particular, @Ratio@ -- types including 'Rational' do not have portable representations.--- +-- -- For the Glasgow compiler (GHC 6.8.2), "GHC.Real" defines @1%0@ -- and @0%0@ as representations for 'infinity' and 'notANumber', -- but most operations on them will raise exceptions. If 'toRational'@@ -29,7 +40,7 @@ -- * <http://www.haskell.org/pipermail/haskell-prime/2006-February/000791.html> -- -- * <http://www.haskell.org/pipermail/haskell-prime/2006-February/000821.html>--- +-- -- Hugs (September 2006) stays closer to the haskell98 spec and -- offers no way of constructing those values, raising arithmetic -- overflow errors if attempted.@@ -68,24 +79,24 @@ -- is compiled correctly. class (PartialOrd a) => Transfinite a where- + -- | A transfinite value which is greater than all finite values. -- Adding or subtracting any finite value is a no-op. As is -- multiplying by any non-zero positive value (including -- @infinity@), and dividing by any positive finite value. Also -- obeys the law @negate infinity = negativeInfinity@ with all -- appropriate ramifications.- + infinity :: a- - ++ -- | A transfinite value which is less than all finite values. -- Obeys all the same laws as @infinity@ with the appropriate -- changes for the sign difference.- + negativeInfinity :: a- - ++ -- | An exceptional transfinite value for dealing with undefined -- results when manipulating infinite values. The following -- operations must return @notANumber@, where @inf@ is any value@@ -118,14 +129,14 @@ -- for 'Eq'; thus, 'eq' and 'ne' are preferred over ('==') and -- ('/=')). Since it returns false for equality, there may be -- more than one machine representation of this `value'.- + notANumber :: a- - ++ -- | Return true for both @infinity@ and @negativeInfinity@, -- false for all other values. isInfinite :: a -> Bool- + -- | Return true only for @notANumber@. isNaN :: a -> Bool @@ -194,6 +205,14 @@ ---------------------------------------------------------------- -- These rules moved here from "LogFloat" in v0.11.2+--+-- FIXME(2023-03-19): Since recent versions of GHC these rules+-- generate warnings that they may never fire because the rule+-- "Class op exp" may fire first. Although the warning suggests+-- adding a phase limit, that doesn't actually help because the+-- "Class op exp" rule is a built-in. So it's not actually clear+-- how to silence this warning:+-- <https://gitlab.haskell.org/ghc/ghc/-/issues/10595> {-# RULES "log/exp" forall x. log (exp x) = x "exp/log" forall x. exp (log x) = x
src/Hugs/RealFloat.hs view
@@ -11,15 +11,15 @@ #define REALFLOAT_VERSION normal Prelude version. This could be buggy. #endif ------------------------------------------------------------------- ~ 2013.05.11+-- ~ 2021.10.17 -- | -- Module : Hugs.RealFloat--- Copyright : Copyright (c) 2007--2015 wren gayle romano+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable (with CPP)--- +-- -- Hugs (September 2006) has buggy definitions for 'Prelude.isNaN' -- and 'Prelude.isInfinite' on Float and Double. If this module is -- run through CPP with the macro @__HUGS__@ set to a value no