logfloat 0.12.1 → 0.13.0.1
raw patch · 9 files changed
+493/−168 lines, 9 filesdep ~basePVP ok
version bump matches the API change (PVP)
Dependency ranges changed: base
API changes (from Hackage documentation)
+ Data.Number.LogFloat: pow :: LogFloat -> Double -> LogFloat
+ Data.Number.LogFloat: product :: [LogFloat] -> LogFloat
+ Data.Number.LogFloat: sum :: [LogFloat] -> LogFloat
- Data.Number.LogFloat: fromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a) => LogFloat -> a
+ Data.Number.LogFloat: fromLogFloat :: LogFloat -> Double
- Data.Number.LogFloat: logFloat :: (Real a, RealToFrac a Double) => a -> LogFloat
+ Data.Number.LogFloat: logFloat :: Double -> LogFloat
- Data.Number.LogFloat: logFromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a) => LogFloat -> a
+ Data.Number.LogFloat: logFromLogFloat :: LogFloat -> Double
- Data.Number.LogFloat: logToLogFloat :: (Real a, RealToFrac a Double) => a -> LogFloat
+ Data.Number.LogFloat: logToLogFloat :: Double -> LogFloat
- Data.Number.PartialOrd: class PartialOrd a
+ Data.Number.PartialOrd: class PartialOrd a where gt x y = case x `cmp` y of { Just GT -> Just True Just _ -> Just False Nothing -> Nothing } ge x y = case x `cmp` y of { Just LT -> Just False Just _ -> Just True Nothing -> Nothing } eq x y = case x `cmp` y of { Just EQ -> Just True Just _ -> Just False Nothing -> Nothing } ne x y = case x `cmp` y of { Just EQ -> Just False Just _ -> Just True Nothing -> Nothing } le x y = case x `cmp` y of { Just GT -> Just False Just _ -> Just True Nothing -> Nothing } lt x y = case x `cmp` y of { Just LT -> Just True Just _ -> Just False Nothing -> Nothing } maxPO x y = do { o <- x `cmp` y; case o of { GT -> Just x EQ -> Just x LT -> Just y } } minPO x y = do { o <- x `cmp` y; case o of { GT -> Just y EQ -> Just x LT -> Just x } }
Files
- INSTALL +125/−0
- LICENSE +1/−1
- VERSION +72/−0
- logfloat.cabal +27/−23
- src/Data/Number/LogFloat.hs +248/−114
- src/Data/Number/PartialOrd.hs +2/−2
- src/Data/Number/RealToFrac.hs +4/−3
- src/Data/Number/Transfinite.hs +6/−14
- src/Hugs/RealFloat.hs +8/−11
+ INSTALL view
@@ -0,0 +1,125 @@+logfloat-0.12+=============++In general, logfloat 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. This document will+walk you through the ins and outs for all supported configurations.++If you run into issues about missing <math.h> or similar, then see+the FFI Problems section. If you want to compile Haddock documentation,+then see that section.+++Building for GHC (6.8 and above)+================================++Installation under GHC is simple using the standard Cabal build+loop:++ runhaskell Setup.hs configure+ runhaskell Setup.hs build+ runhaskell Setup.hs copy+ runhaskell Setup.hs register++The logfloat package is no longer supported for GHC 6.6, though it+could probably be made to work with it by changing LANGUAGE pragma+into the appropriate OPTIONS_GHC equivalents. Since it's not supported+anymore, more significant changes may have become necessary as well.+++Building for Hugs (September 2006)+==================================++Installation under Hugs is unfortunately more complicated due to+some option-passing bugs in Cabal. You should be able to build+things with these modifications to the usual Cabal build loop:++ runhaskell Setup.hs configure --hugs \+ --with-cpphs="`which cpphs-hugs`" \+ --ffihugs-options="-98 +o" \+ --ffihugs-option=-F'cpp -P -traditional -D__HUGS__=200609 -D__USE_FFI__'+ runhaskell Setup.hs build+ runhaskell Setup.hs copy+ runhaskell Setup.hs register++Hugs installs cpphs under the name "cpphs-hugs" by default. The+`which` command will try to resolve the location, assuming it's on+your $PATH. If it isn't, then change the --with-cpphs= flag to point+to where cpphs is installed on your system.++Some additional details about difficulties with building+for Hugs can be found in this blog post and the bug reports:++ http://winterkoninkje.livejournal.com/60707.html+ http://hackage.haskell.org/trac/hackage/ticket/526+ http://hackage.haskell.org/trac/hackage/ticket/527++And while previous versions of Cabal could build this package (with+the above commandline to work around those bugs), Hugs support is+broken in Cabal 1.8 (including at least 1.8.0.2 through 1.8.0.6):++ http://hackage.haskell.org/trac/hackage/ticket/633+++FFI Problems+============++The logfloat package uses the FFI to access functions in libm to+improve accuracy. These functions are part of the ISO C 1999 standard+and are supported on most systems, however they're not part of the+ANSI C standard nor the System V standard and so they might be+unavailable on some systems. If you encounter errors during build+about not being able to find <math.h>, not having prototypes for+functions, or linking errors, then you are probably on such a system.++In order to use pure Haskell and disable the FFI, pass -f-useFFI+during the configure phase to disable the useFFI flag. (If compiling+for Hugs, also be sure *not* to pass -D__USE_FFI__ to the cpp filter+for ffihugs.) This will degrade the accuracy of certain operations,+but should still compile cleanly.++If you run into errors like "error: 'log1p' redeclared as different+kind of symbol" this is due to a bug in GHC 6.10.1 where you can't+use both -fvia-C and the FFI. This bug should be circumvented for+the moment, but resolving it is still a work in progress. See the+bug report (which is resolved in GHC 6.10.2):++ http://hackage.haskell.org/trac/ghc/ticket/3117+++Windows FFI+===========++The logfloat package builds and installs cleanly on Windows XP using+GHC 6.10.1 (without needing Cygwin nor Mingw/Msys). The C functions+reside in libmingwex.a which is bundled with GHC on WinXP.++However, it does not work from the GHCi debugger and gives errors+like "can't load .so/.DLL for: m (addDLL: could not load DLL)". If+you know how to fix this issue or have further details on Windows+compatability, please contact the maintainer.+++Haddock Documentation+=====================++In Cabal 1.2 there is a bug in the handling of building Haddock+documentation when CPP is involved. These issues have been fixed+in Cabal 1.6, but here are the instructions if you're on older+systems.++In Cabal 1.2 the cpp-options field is not passed to Haddock, and+therefore any macros defined there are not seen, which can cause+Haddock to fail. The old workaround was to define CPP macros in the+ghc-options field which *does* get passed to Haddock. This is now+considered bad style and is forbidden by Hackage. It appears that+passing --haddock-option=... or --ghc-option=... flags during+configure does not have the same effect as defining the field.++Therefore, in order to properly compile Haddock documentation on+Cabal 1.2, you should go into logfloat.cabal and uncomment the+ghc-options fields which declare CPP macros -D__USE_FFI__ and+-D__HUGS__=200609. If you know of a better workaround for this+configuration, contact the maintainer.
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2007--2010, wren ng thornton.+Copyright (c) 2007--2015, wren gayle romano. ALL RIGHTS RESERVED. Redistribution and use in source and binary forms, with or without
+ VERSION view
@@ -0,0 +1,72 @@+0.13 (2015-02-17):+ - monomorphized logFloat, logToLogFloat, fromLogFloat, logFromLogFloat+ - added sum, product, pow+ - cleaned up some extraneous rewrite rules+ - removed support for older versions of GHC, in order to clean things up++0.12.1 (2010-03-19):+ - Fixed some NaN injection bugs in the Num instance. These would only have been triggered by non-probabilistic uses of LogFloat (i.e., if you used @logFloat infinity@)+0.12.0.3 (2009-04-03):+ - Added notes to INSTALL about Windows compatibility+0.12.0.2 (2009-04-03):+ - Updated TODO and VERSION+0.12.0.1 (2009-03-25):+ - Disabled -fvia-C so the FFI can be used in GHC 6.10+0.12 (2009-03-10):+ - Added log1p.+ - Added Storable instance (for GHC only)+ - removed orphaned toRational/fromRational rules. The RealToFrac module obviates those optimizations.+ - Adjusted Real LogFloat instance to throw errors on transfinite values++0.11.2 (2009-03-09):+ - Moved log/exp rules from LogFloat into Transfinite+0.11.1.1 (2009-03-09):+ - Documentation updates+0.11.1 (2009-03-08):+ - Added IArray UArray LogFloat instance (thanks to Felipe Lessa).+0.11 (2009-01-29):+ - Moved the RealToFrac class from Data.Number.Transfinite to Data.Number.RealToFrac. This breaks backwards compatibility if you depended directly on that module for the class. If you imported the class from Data.Number.LogFloat then you're still fine.+ - Cabal build fail due to missing new module (Fixed in 0.11.0.1).++0.10 (2009-01-29):+ - Hugs' Prelude definitions of isInfinite and isNaN for Float and Double are buggy. The new Hugs.RealFloat module provides correct implementations (defaulting to the Prelude definition for non-Hugs compilers).+ - The PartialOrd and Transfinite instances for Double and Float have been updated to use these corrected functions.+ - Added maxPO and minPO to PartialOrd and added comparingPO to Data.Number.PartialOrd.+ - Minor changes to circumvent bugs in Hackage's new version of Haddock. You should be able to see all the documentation now.++0.9.1 (2008-08-30):+ - Fixed some PartialOrd stuff and sanitized documentation+ - Build fail on GHC 6.10 due to GHC.Prim (Fixed in 0.9.1.1).+0.9.0 (2008-08-30):+ - s/toFractional/realToFrac/g.+ - Also moved realToFrac and log to Transfinite+ - Build fail on GHC 6.10 due to GHC.Prim++0.8.6 (2008-08-17):+ - Removed buggy RULES+ - Build fail on GHC 6.10 due to -Werror and removal of -fno-warn-orphans+0.8.5 (2008-08-17):+ - Gave up and converted from lhs to hs so Hackage docs work+0.8.4 (2008-08-17):+ - Broke out Transfinite+0.8.3 (2008-08-16):+ - Documentation updates+0.8.2 (2008-08-14):+ - Announced release+0.8:+ - Did a bunch of tweaking. Things should be decent now++0.7:+ - Haddockified+0.6:+ - Fixed monomorphism of LogFloat constructors/destructors.+0.5:+ - Added optimization rules.+0.4:+ - Translated from Perl to Haskell at revision 2007.12.20.+0.3:+ - Converted extensive comments to POD format.+0.2:+ - Did a bunch of profiling, optimizing, and debugging.+0.1:+ - Initial version (in Perl) created for hw5 for NLP with Jason Eisner.
logfloat.cabal view
@@ -1,18 +1,22 @@ ------------------------------------------------------------------- wren ng thornton <wren@community.haskell.org> ~ 2010.03.19+-- wren gayle romano <wren@community.haskell.org> ~ 2015.02.17 ---------------------------------------------------------------- -Name: logfloat-Version: 0.12.1-Cabal-Version: >= 1.2+-- By and large Cabal >=1.2 is fine; but >= 1.6 gives tested-with:+-- and source-repository:.+Cabal-Version: >= 1.6 Build-Type: Simple++Name: logfloat+Version: 0.13.0.1 Stability: experimental-Copyright: Copyright (c) 2007--2010 wren ng thornton+Homepage: http://code.haskell.org/~wren/+Author: wren gayle romano+Maintainer: wren@community.haskell.org+Copyright: Copyright (c) 2007--2015 wren gayle romano License: BSD3 License-File: LICENSE-Author: wren ng thornton-Maintainer: wren@community.haskell.org-Homepage: http://code.haskell.org/~wren/+ Category: Data, Math, Natural Language Processing, Statistics Synopsis: Log-domain floating point numbers Description: This module presents a type for storing numbers@@ -22,16 +26,23 @@ It is also helpful for preventing overflow. -Flag splitBase- Description: base-3.0 broke out array and other packages- Default: False+-- TODO: does GHC 6.12.1 and 7.6.1 still work? do we care?+Tested-With:+ GHC == 7.8.3+Extra-source-files:+ INSTALL, VERSION+Source-Repository head+ Type: darcs+ Location: http://community.haskell.org/~wren/logfloat +---------------------------------------------------------------- Flag useFFI+ Default: True Description: Use FFI to link C's log1p. Improves accuracy and same performance. Disable only if your compiler doesn't support FFI.- Default: True +---------------------------------------------------------------- Library Hs-Source-Dirs: src Exposed-Modules: Data.Number.LogFloat@@ -39,24 +50,17 @@ , Data.Number.Transfinite , Data.Number.PartialOrd , Hugs.RealFloat- if flag(splitBase)- Build-depends: base >= 3.0, array- else- Build-depends: base < 3.0 + Build-Depends: base == 4.*+ -- TODO: do we need any version constraints on array?+ , array+ if flag(useFFI) -- BUG (Cabal 1.2 + Haddock): See the INSTALL file. --GHC-Options: -D__USE_FFI__ CPP-Options: -D__USE_FFI__ includes: math.h extra-libraries: m- - -- BUG (GHC 6.10.1): Can't mix FFI and -fvia-C- -- <http://hackage.haskell.org/trac/ghc/ticket/3117>- -- This is a momentary hack to keep it for 6.8- -- TODO: see if -fasm suffices for the benefits- if impl(ghc < 6.10)- GHC-Options: -fvia-C -optc-O3 -- BUG (Cabal <= 1.6 + Hugs): See the INSTALL file. Hugs-Options: -98 +o
src/Data/Number/LogFloat.hs view
@@ -1,30 +1,18 @@--- FlexibleContexts needed by our RealToFrac contexts--- CPP needed for IArray UArray instance+-- CPP and GeneralizedNewtypeDeriving are needed for IArray UArray instance -- FFI is for log1p--- N.B. can't mix FFI and -fvia-C under ghc==6.10.1--- <http://hackage.haskell.org/trac/ghc/ticket/3117>--- TODO: figure out how to conditionally use "-fvia-C -optc-O3"--- for ghc6.8 when not using Cabal to compile.-{-# LANGUAGE FlexibleContexts- , CPP- , ForeignFunctionInterface- #-}+{-# LANGUAGE CPP, ForeignFunctionInterface #-}+-- We don't put these in LANGUAGE, because it's CPP guarded for GHC only+{-# OPTIONS_GHC -XGeneralizedNewtypeDeriving #-} --- Removed -Wall because -fno-warn-orphans was removed in GHC 6.10 {-# OPTIONS_GHC -Wall -fwarn-tabs #-} --- Unfortunately we need -fglasgow-exts in order to actually pick--- up on the rules (see -ddump-rules). The -frewrite-rules flag--- doesn't do what you want.--- cf <http://hackage.haskell.org/trac/ghc/ticket/2213>--- cf <http://www.mail-archive.com/glasgow-haskell-users@haskell.org/msg14313.html>-{-# OPTIONS_GHC -O2 -fexcess-precision -fglasgow-exts #-}+{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-} ------------------------------------------------------------------- ~ 2010.03.19+-- ~ 2015.02.17 -- | -- Module : Data.Number.LogFloat--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2015 wren gayle romano -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable@@ -51,7 +39,6 @@ ( -- * Exceptional numeric values module Data.Number.Transfinite- , module Data.Number.RealToFrac -- * @LogFloat@ data type , LogFloat@@ -61,14 +48,17 @@ -- ** Isomorphism to log-domain , logToLogFloat , logFromLogFloat+ -- ** Additional operations+ , sum, product+ , pow -- * Accurate versions of logarithm\/exponentiation , log1p, expm1 ) where -import Prelude hiding (log, realToFrac, isInfinite, isNaN)+import Prelude hiding (log, sum, product, isInfinite, isNaN)+import Data.List (foldl') -import Data.Number.RealToFrac import Data.Number.Transfinite import Data.Number.PartialOrd @@ -95,9 +85,24 @@ -- | 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 so will+-- instance will convert back to the normal-domain, and hence will -- underflow at that point. 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@@ -105,17 +110,22 @@ -- 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.+-- 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 it's necessary at times--- and the proper implementation is not immediately transparent.--- However, between two @LogFloat@s addition requires crossing the--- exp\/log boundary twice; with a @LogFloat@ and a regular number--- 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 plain operations first, do it!+-- 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@@ -141,9 +151,9 @@ infixr 2 ~> f ~> g = (. f) . (g .) -{-# INLINE ($.) #-}-infixl 1 $.-($.) = flip ($)+{-# INLINE ($::) #-}+infixl 1 $::+($::) = flip ($) {-# INLINE logFromLFAssocs #-}@@ -162,9 +172,9 @@ -- Named unsafe because it could allow injecting NaN if misused {-# INLINE unsafeLogToLFFunc #-} unsafeLogToLFFunc :: (LogFloat -> a -> LogFloat) -> (Double -> a -> Double)-unsafeLogToLFFunc = ($. unsafeLogToLogFloat ~> id ~> logFromLogFloat)+unsafeLogToLFFunc = ($:: unsafeLogToLogFloat ~> id ~> logFromLogFloat) --- | Remove the extranious 'isNaN' test of 'logToLogFloat', when+-- | Remove the extraneous 'isNaN' test of 'logToLogFloat', when -- we know it's safe. {-# INLINE unsafeLogToLogFloat #-} unsafeLogToLogFloat :: Double -> LogFloat@@ -184,30 +194,30 @@ {-# INLINE unsafeArray #-} unsafeArray =- unsafeArray $. id ~> logFromLFAssocs ~> unsafeLogToLFUArray+ unsafeArray $:: id ~> logFromLFAssocs ~> unsafeLogToLFUArray {-# INLINE unsafeAt #-} unsafeAt =- unsafeAt $. logFromLFUArray ~> id ~> unsafeLogToLogFloat+ unsafeAt $:: logFromLFUArray ~> id ~> unsafeLogToLogFloat {-# INLINE unsafeReplace #-} unsafeReplace =- unsafeReplace $. logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray+ unsafeReplace $:: logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray {-# INLINE unsafeAccum #-} unsafeAccum =- unsafeAccum $. unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray+ unsafeAccum $:: unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray {-# INLINE unsafeAccumArray #-} unsafeAccumArray =- unsafeAccumArray $. unsafeLogToLFFunc ~> logFromLogFloat ~> id ~> id ~> unsafeLogToLFUArray+ 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) + cmp (LogFloat x) (LogFloat y) | isNaN x || isNaN y = Nothing | otherwise = Just $! x `compare` y @@ -216,82 +226,75 @@ -- | 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"+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+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+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 input.-logFloat :: (Real a, RealToFrac a Double) => a -> LogFloat-{-# SPECIALIZE logFloat :: Double -> LogFloat #-}-logFloat = LogFloat . log . guardNonNegative "logFloat" . realToFrac+-- 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" --- This is simply a polymorphic version of the 'LogFloat' data--- constructor. We present it mainly because we hide the constructor--- in order to make the type a bit more opaque. If the polymorphism--- turns out to be a performance liability because the rewrite rules--- can't remove it, then we need to rethink all four--- constructors\/destructors.--- -- | Constructor which assumes the argument is already in the--- log-domain. Throws errors on @notANumber@ input.-logToLogFloat :: (Real a, RealToFrac a Double) => a -> LogFloat-{-# SPECIALIZE logToLogFloat :: Double -> LogFloat #-}-logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" . realToFrac+-- log-domain. Throws errors on @notANumber@ inputs.+logToLogFloat :: Double -> LogFloat+logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" --- | Return our log-domain value back into normal-domain. Beware--- of overflow\/underflow.-fromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a)- => LogFloat -> a-{-# SPECIALIZE fromLogFloat :: LogFloat -> Double #-}-fromLogFloat (LogFloat x) = realToFrac (exp x)+-- | 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 :: (Fractional a, Transfinite a, RealToFrac Double a)- => LogFloat -> a-{-# SPECIALIZE logFromLogFloat :: LogFloat -> Double #-}-logFromLogFloat (LogFloat x) = realToFrac x+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 may need to delay inlining--- of the four con-\/destructors, like we do for 'realToFrac'.--- Unfortunately, I'm not entirely sure whether they will be inlined--- already or not (and whether they are may be fragile) and I don't--- want to inline them excessively and lead to code bloat in the--- off chance that we could prune some of it away.--- TODO: thoroughly investigate this.+-- 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-"log.fromLogFloat" log . fromLogFloat = logFromLogFloat- "logFloat/fromLogFloat" forall x. logFloat (fromLogFloat x) = x-"logFloat.fromLogFloat" logFloat . fromLogFloat = id -- Into log-domain and back out "fromLogFloat/logFloat" forall x. fromLogFloat (logFloat x) = x-"fromLogFloat.logFloat" fromLogFloat . logFloat = id #-} ----------------------------------------------------------------@@ -305,8 +308,9 @@ -- log-domain value instead. instance Show LogFloat where- show (LogFloat x) = let y = exp x- in y `seq` "LogFloat "++show y+ show (LogFloat x) =+ let y = exp x+ in y `seq` "LogFloat "++show y ----------------------------------------------------------------@@ -316,34 +320,44 @@ -- be safe. #ifdef __USE_FFI__-#define LOG1P_WHICH_VERSION specialized version.+#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+)@. The C language provides 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@.+-- | 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 log1p #-}+{-# INLINE [0] log1p #-} log1p x = log (1 + x) #endif --- | Definition: @expm1 == (subtract 1) . exp@. The C language--- provides a special definition for 'expm1' which is more accurate+-- | 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'.@@ -355,10 +369,19 @@ expm1 :: Double -> Double #else expm1 :: Double -> Double-{-# INLINE expm1 #-}+{-# INLINE [0] expm1 #-} expm1 x = exp x - 1 #endif ++{-# 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+ #-}+ ---------------------------------------------------------------- -- These all work without causing underflow. However, do note that -- they tend to induce more of the floating-point fuzz than using@@ -368,7 +391,7 @@ -- multiplying many large numbers) so we're not too worried about -- +\/- 4e-16. -instance Num LogFloat where +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'@@ -391,7 +414,9 @@ | x == negativeInfinity && y == negativeInfinity = LogFloat negativeInfinity -- @0-0 == 0@ | otherwise =- -- Will throw error if x < y or if (x,y) is (infinity,infinity)+ -- 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)@@ -412,7 +437,8 @@ 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@@ -424,7 +450,7 @@ . guardNonNegative "fromRational" . fromRational --- Just for fun. The more coersion functions the better. Though+-- 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)@@ -434,17 +460,125 @@ ex = exp x -{- -- Commented out because I'm not sure about requiring MPTCs. Of course, those are already required by "Data.Number.Transfinite" so it's pretty moot...+----------------------------------------------------------------+-- | /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 --- LogFloat->LogFloat is already given via generic (a->a)--- No LogFloat->Rational since LogFloat can have 'infinity'--- Can't have LogFloat->a using fromLogFloat because Hugs dislikes incoherence. Adding an explicit LogFloat->LogFloat instance doesn't help like it does for GHC. -instance RealToFrac LogFloat Double where- realToFrac = fromLogFloat+-- 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 + theSum)+ where+ theMax = logFromLogFloat (maximum xs) -instance RealToFrac LogFloat Float where- realToFrac = fromLogFloat+ -- compute @\log \sum_{x \in xs} \exp(x - theMax)@+ theSum = foldl' (\ acc x -> acc + exp (logFromLogFloat 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) =+ -- 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) -} ----------------------------------------------------------------
src/Data/Number/PartialOrd.hs view
@@ -9,7 +9,7 @@ -- ~ 2009.01.29 -- | -- Module : Data.Number.PartialOrd--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2015 wren gayle romano -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable@@ -18,7 +18,7 @@ -- 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--- are moreso. This module presents a class for partially ordered+-- are more so. This module presents a class for partially ordered -- types. ---------------------------------------------------------------- module Data.Number.PartialOrd
src/Data/Number/RealToFrac.hs view
@@ -5,16 +5,17 @@ , OverlappingInstances , FlexibleInstances , CPP- , MagicHash #-}+-- We don't put these in LANGUAGE, because it's CPP guarded for GHC only+{-# OPTIONS_GHC -XMagicHash #-} {-# OPTIONS_GHC -Wall -fwarn-tabs #-} ------------------------------------------------------------------- ~ 2009.01.29+-- ~ 2013.05.11 -- | -- Module : Data.Number.RealToFrac--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2015 wren gayle romano -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable
src/Data/Number/Transfinite.hs view
@@ -1,17 +1,12 @@ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} --- Unfortunately we need -fglasgow-exts in order to actually pick--- up on the rules (see -ddump-rules). The -frewrite-rules flag--- doesn't do what you want.--- <http://hackage.haskell.org/trac/ghc/ticket/2213>--- <http://www.mail-archive.com/glasgow-haskell-users@haskell.org/msg14313.html>-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# OPTIONS_GHC -O2 -fenable-rewrite-rules #-} ------------------------------------------------------------------- ~ 2009.03.09+-- ~ 2013.05.11 -- | -- Module : Data.Number.Transfinite--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2015 wren gayle romano -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable@@ -112,11 +107,11 @@ -- -- * @inf \/ inf@ --- -- * @inf `div` inf@+ -- * @inf \`div\` inf@ -- -- * @0 \/ 0@ --- -- * @0 `div` 0@+ -- * @0 \`div\` 0@ -- -- Additionally, any mathematical operations on @notANumber@ -- must also return @notANumber@, and any equality or ordering@@ -163,7 +158,7 @@ -- This function will raise an error when taking the log of negative -- numbers, rather than returning 'notANumber' as the newer GHC -- implementation does. The reason being that typically this is a--- logical error, and @notANumber@ allows the error to propegate+-- logical error, and @notANumber@ allows the error to propagate -- silently. -- -- In order to improve portability, the 'Transfinite' class is@@ -203,10 +198,7 @@ -- These rules moved here from "LogFloat" in v0.11.2 {-# RULES "log/exp" forall x. log (exp x) = x-"log.exp" log . exp = id- "exp/log" forall x. exp (log x) = x-"exp.log" exp . log = id #-} -- We'd like to be able to take advantage of general rule versions
src/Hugs/RealFloat.hs view
@@ -11,10 +11,10 @@ #define REALFLOAT_VERSION normal Prelude version. This could be buggy. #endif ------------------------------------------------------------------- ~ 2010.03.19+-- ~ 2013.05.11 -- | -- Module : Hugs.RealFloat--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2015 wren gayle romano -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable@@ -44,9 +44,7 @@ import qualified Prelude ---------------------------------------------------------------- -isInfinite :: (RealFloat a) => a -> Bool-{-# SPECIALIZE isInfinite :: Double -> Bool #-}-{-# SPECIALIZE isInfinite :: Float -> Bool #-}+isInfinite :: (RealFloat a) => a -> Bool {-# INLINE isInfinite #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isInfinite x = (1/0) == abs x@@ -56,8 +54,6 @@ isNaN :: (RealFloat a) => a -> Bool-{-# SPECIALIZE isNaN :: Double -> Bool #-}-{-# SPECIALIZE isNaN :: Float -> Bool #-} {-# INLINE isNaN #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isNaN x = compareEQ x 0 && compareEQ x 1@@ -65,10 +61,11 @@ -- | In Hugs (September 2006), 'compare' always returns @EQ@ if one -- of the arguments is not a number. Thus, if a number is @compareEQ@ -- against multiple different numbers, then it must be @isNaN@.-compareEQ :: (Ord a) => a -> a -> Bool-compareEQ x y = case compare x y of- EQ -> True- _ -> False+compareEQ :: (Ord a) => a -> a -> Bool+compareEQ x y =+ case compare x y of+ EQ -> True+ _ -> False #else isNaN = Prelude.isNaN #endif