diff --git a/CHANGELOG b/CHANGELOG --- a/CHANGELOG +++ b/CHANGELOG @@ -1,6 +1,21 @@ +0.14.0.2 (2026-02-28): + - Nudged version bounds for GHC 9.14 (didn't actually need to + change the upper bounds on 'base' and 'array', because they + were already lenient) + - Silenced GHC >= 9.10 warnings about importing 'Data.List' +0.14.0.1 (2024-08-29): + - Nudged version bounds for ghc 9.10 +0.14.0 (2021-03-19): + - 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: + 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): diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ logfloat ======== +[![CI Status](https://github.com/wrengr/logfloat/actions/workflows/ci.yml/badge.svg)](https://github.com/wrengr/logfloat/actions?query=workflow%3Aci+-event%3Apull_request) [![Hackage version](https://img.shields.io/hackage/v/logfloat.svg?style=flat)](https://hackage.haskell.org/package/logfloat) -[![Hackage-Deps](https://img.shields.io/hackage-deps/v/logfloat.svg?style=flat)](http://packdeps.haskellers.com/specific?package=logfloat) -[![TravisCI Build Status](https://img.shields.io/travis/wrengr/logfloat.svg?style=flat)](https://travis-ci.org/wrengr/logfloat) -[![CircleCI Build Status](https://circleci.com/gh/wrengr/logfloat.svg?style=shield&circle-token=b57517657c556be6fd8fca92b843f9e4cffaf8d1)](https://circleci.com/gh/wrengr/logfloat) +[![Stackage LTS version](https://stackage.org/package/logfloat/badge/lts)](https://stackage.org/lts/package/logfloat) +[![Stackage Nightly version](https://stackage.org/package/logfloat/badge/nightly)](https://stackage.org/nightly/package/logfloat) This package provides a type for storing numbers in the log-domain, primarily useful for preventing underflow when multiplying many @@ -19,38 +19,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 +231,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) diff --git a/logfloat.cabal b/logfloat.cabal --- a/logfloat.cabal +++ b/logfloat.cabal @@ -1,19 +1,25 @@ +Cabal-Version: 2.2 +-- Cabal >=2.2 is required for: +-- +-- Since 2.1, the Cabal-Version must be the absolutely first thing +-- in the file, even before comments. Also, no longer uses ">=". +-- + ---------------------------------------------------------------- --- wren gayle romano ~ 2021.10.16 +-- wren gayle romano ~ 2026-02-26 ---------------------------------------------------------------- --- Cabal >=1.10 is required by Hackage. -Cabal-Version: >= 1.10 -Build-Type: Simple - Name: logfloat -Version: 0.13.4 -Stability: experimental -Homepage: http://wrengr.org +Version: 0.14.0.2 +Build-Type: Simple +Stability: provisional +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–2026 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,19 @@ 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.8, + GHC ==9.6.5, + GHC ==9.8.2, + GHC ==9.10.1, + GHC ==9.12.1, + GHC ==9.14.1 ---------------------------------------------------------------- Source-Repository head Type: git - Location: git://github.com/wrengr/logfloat.git + Location: https://github.com/wrengr/logfloat.git ---------------------------------------------------------------- Flag useFFI @@ -57,6 +70,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 diff --git a/src/Data/Number/LogFloat.hs b/src/Data/Number/LogFloat.hs --- a/src/Data/Number/LogFloat.hs +++ b/src/Data/Number/LogFloat.hs @@ -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. -- @@ -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. --- Also see the other --- implementations given there. For Python's actual C implementation, --- see math_fsum in --- --- --- For merely *mitigating* errors rather than completely eliminating --- them, see . --- --- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@ -{- --- For proof of correctness, see --- -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 - return sum(partials, 0.0) --} +product = LogFloat . kahanSum . fmap logFromLogFloat ---------------------------------------------------------------- ----------------------------------------------------------- fin. diff --git a/src/Data/Number/LogFloat/Raw.hs b/src/Data/Number/LogFloat/Raw.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/LogFloat/Raw.hs @@ -0,0 +1,405 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface, BangPatterns #-} +{-# OPTIONS_GHC -Wall -fwarn-tabs #-} +{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-} + +---------------------------------------------------------------- +-- ~ 2026-02-28 +-- | +-- Module : Data.Number.LogFloat.Raw +-- Copyright : 2007--2026 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: +-- +-- +-- 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 + +#if __GLASGOW_HASKELL__ < 910 +import Data.List (foldl') +#endif +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. +-- Also see the other +-- implementations given there. For Python's actual C implementation, +-- see math_fsum in +-- +-- +-- For merely *mitigating* errors rather than completely eliminating +-- them, see . +-- +-- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@ +{- +-- For proof of correctness, see +-- +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 + return sum(partials, 0.0) +-} + + +---------------------------------------------------------------- +----------------------------------------------------------- fin. diff --git a/src/Data/Number/PartialOrd.hs b/src/Data/Number/PartialOrd.hs --- a/src/Data/Number/PartialOrd.hs +++ b/src/Data/Number/PartialOrd.hs @@ -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 diff --git a/src/Data/Number/RealToFrac.hs b/src/Data/Number/RealToFrac.hs --- a/src/Data/Number/RealToFrac.hs +++ b/src/Data/Number/RealToFrac.hs @@ -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 #-} diff --git a/src/Data/Number/Transfinite.hs b/src/Data/Number/Transfinite.hs --- a/src/Data/Number/Transfinite.hs +++ b/src/Data/Number/Transfinite.hs @@ -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: +-- +-- +-- 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 @@ -- * -- -- * --- +-- -- 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: +-- {-# RULES "log/exp" forall x. log (exp x) = x "exp/log" forall x. exp (log x) = x diff --git a/src/Hugs/RealFloat.hs b/src/Hugs/RealFloat.hs --- a/src/Hugs/RealFloat.hs +++ b/src/Hugs/RealFloat.hs @@ -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