diff --git a/logfloat.cabal b/logfloat.cabal --- a/logfloat.cabal +++ b/logfloat.cabal @@ -1,9 +1,9 @@ ---------------------------------------------------------------- --- wren ng thornton ~ 2009.01.29 +-- wren ng thornton ~ 2009.03.17 ---------------------------------------------------------------- Name: logfloat -Version: 0.11.1 +Version: 0.12 Cabal-Version: >= 1.2 Build-Type: Simple Stability: experimental @@ -26,6 +26,11 @@ Description: base-3.0 broke out array and other packages Default: False +Flag useFFI + 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 @@ -39,9 +44,19 @@ else Build-depends: base < 3.0 - Hugs-Options: -98 +o -F'cpp -P -traditional -D__HUGS__=200609' - if impl(ghc < 6.10) - GHC-Options: -fno-warn-orphans + 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: (Cabal <= 1.6 + Hugs) See the INSTALL file. + Hugs-Options: -98 +o + if impl(hugs) + -- BUG: (Cabal 1.2 + Haddock) See the INSTALL file. + --GHC-Options: -D__HUGS__=200609 + CPP-Options: -D__HUGS__=200609 ---------------------------------------------------------------- ----------------------------------------------------------- fin. 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 @@ -1,11 +1,14 @@ -- FlexibleContexts needed by our RealToFrac contexts -- CPP needed for IArray UArray instance +-- FFI is for log1p {-# LANGUAGE FlexibleContexts - , CPP #-} + , CPP + , ForeignFunctionInterface + #-} -- Removed -Wall because -fno-warn-orphans was removed in GHC 6.10 -{-# OPTIONS_GHC -fwarn-tabs #-} +{-# 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 @@ -14,37 +17,15 @@ -- cf {-# OPTIONS_GHC -O2 -fvia-C -optc-O3 -fexcess-precision -fglasgow-exts #-} --- Version History --- (v0.11.1) Added IArray UArray instance --- (v0.11) Broke Data.Number.RealToFrac out --- (v0.10) Fixed bugs in Hugs for PartialOrd and Transfinite. --- Also added maxPO, minPO, comparingPO --- (v0.9.1) Fixed some PartialOrd stuff and sanitized documentation --- (v0.9.0) s/toFractional/realToFrac/g. --- Also moved realToFrac and log to Transfinite --- (v0.8.6) Removed buggy RULES --- (v0.8.5) Gave up and converted from lhs to hs so Hackage docs work --- (v0.8.4) Broke out Transfinite --- (v0.8.3) Documentation updates --- (v0.8.2) Announced release --- (v0.8) Did a bunch of tweaking. Things should be decent now --- (v0.7) Haddockified --- (v0.6) Fixed monomorphism. --- (v0.5) Added optimization rules. --- (v0.4) Translated to Haskell at revision 2007.12.20. --- (v0.3) Converted extensive comments to POD format. --- (v0.2) Did a bunch of profiling, optimizing, and debugging. --- (v0.1) Initial version created for hw5 for NLP with Jason Eisner. --- ---------------------------------------------------------------- --- ~ 2009.03.07 +-- ~ 2009.03.10 -- | -- Module : Data.Number.LogFloat -- Copyright : Copyright (c) 2007--2009 wren ng thornton -- License : BSD3 -- Maintainer : wren@community.haskell.org -- Stability : stable --- Portability : portable (with CPP) +-- Portability : portable (with CPP, FFI) -- -- This module presents a type for storing numbers in the log-domain. -- The main reason for doing this is to prevent underflow when @@ -61,7 +42,7 @@ -- The 'LogFloat' of this module is restricted to non-negative -- numbers for efficiency's sake, see the forthcoming -- "Data.Number.LogFloat.Signed" for doing signed log-domain --- calculations. +-- calculations. (Or harass the maintainer to write it already.) ---------------------------------------------------------------- module Data.Number.LogFloat @@ -70,10 +51,17 @@ module Data.Number.Transfinite , module Data.Number.RealToFrac - -- * @LogFloat@ data type and conversion functions + -- * @LogFloat@ data type , LogFloat - , logFloat, logToLogFloat - , fromLogFloat, logFromLogFloat + -- ** Isomorphism to normal-domain + , logFloat + , fromLogFloat + -- ** Isomorphism to log-domain + , logToLogFloat + , logFromLogFloat + + -- * Accurate versions of logarithm\/exponentiation + , log1p, expm1 ) where import Prelude hiding (log, realToFrac, isInfinite, isNaN) @@ -96,73 +84,9 @@ import NonStdUnsafeCoerce (unsafeCoerce) #endif ----------------------------------------------------------------- --- --- Try to add in some optimizations. Why these need to be --- down here and localized to the module, I don't know. We don't --- do anything foolish like this, but our clients might, or they --- might be generated by other code transformations. Note that due --- to the fuzz, these equations are not strictly true, even though --- they are mathematically correct. - -{-# 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 --- of our operators for 'LogFloat', with rules like @log x + log y --- = log (x * y)@ and @log x - log y = log (x / y)@. However the --- problem is that those equations could be driven in either direction --- depending on whether we think time performance or non-underflow --- performance is more important, and the answers may be different --- at every call site. --- --- Since we implore users to do normal-domain computations whenever --- it would not degenerate accuracy, we should not rewrite their --- decisions in any way. The log\/exp fusion strictly improves both --- time and accuracy, so those are safe. But the buck stops with --- them. - - --- These should only fire when it's type-safe --- This should already happen, but... --- TODO: Check the logs to see if it ever fires --- N.B. these are orphaned -{-# RULES -"toRational/fromRational" forall x. toRational (fromRational x) = x -"toRational.fromRational" toRational . fromRational = id - #-} - - ----------------------------------------------------------------- - --- | Reduce the number of constant string literals we need to store. -errorOutOfRange :: String -> a -errorOutOfRange fun = error $! "Data.Number.LogFloat."++fun - ++ ": argument out of range" - - --- | We need these guards in order to ensure some invariants. -guardNonNegative :: String -> Double -> Double -guardNonNegative fun x | x >= 0 = x - | otherwise = errorOutOfRange fun - - --- TODO: since we're using Hugs.RealFloat instead of Prelude now, --- is it still non-portable? --- --- | It's unfortunate that 'notANumber' is not equal to itself, but --- we can hack around that. GHC gives NaN for the log of negatives --- and so we could ideally take advantage of @log . guardNonNegative --- fun = guardIsANumber fun . log@ to simplify things, but Hugs --- raises an error so that's non-portable. -guardIsANumber :: String -> Double -> Double -guardIsANumber fun x | isNaN x = errorOutOfRange fun - | otherwise = x +#ifdef __GLASGOW_HASKELL__ +import Foreign.Storable (Storable) +#endif ---------------------------------------------------------------- -- @@ -196,15 +120,17 @@ ( Eq , Ord -- Should we really perpetuate the Ord lie? #ifdef __GLASGOW_HASKELL__ - , IArray UArray - -- At least GHC 6.8.2 can derive IArray UArray (without + -- At least GHC 6.8.2 can derive these (without -- GeneralizedNewtypeDeriving). The H98 Report doesn't include - -- that among the options for automatic derivation though. + -- them among the options for automatic derivation though. + , IArray UArray + , Storable #endif ) #if __HUGS__ || __NHC__ +-- TODO: Storable instance. Though Foreign.Storable isn't in Hugs(Sept06) -- These two operators make it much easier to read the instance. -- Hopefully inlining everything will get rid of the eta overhead. @@ -235,7 +161,7 @@ logToLFFunc = ($. unsafeLogToLogFloat ~> id ~> logFromLogFloat) -- | Remove the extranious 'isNaN' test of 'logToLogFloat', when --- we know we can. +-- we know it's safe. {-# INLINE unsafeLogToLogFloat #-} unsafeLogToLogFloat :: Double -> LogFloat unsafeLogToLogFloat = LogFloat @@ -281,8 +207,31 @@ ---------------------------------------------------------------- --- | A constructor which does semantic conversion from normal-domain --- to log-domain. +-- | 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" + + +-- | We need these guards in order to ensure some invariants. +guardNonNegative :: String -> Double -> Double +guardNonNegative fun x | x >= 0 = x + | otherwise = errorOutOfRange fun + + +-- | In general @log . guardNonNegative fun == guardIsANumber fun . log@ +-- This even holds on Hugs now. However, the latter issues an error +-- from 'Data.Number.Transfinite.log' whereas the former issues an +-- error from the calling function and is therefore more helpful +-- for debugging. +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 @@ -296,7 +245,7 @@ -- constructors\/destructors. -- -- | Constructor which assumes the argument is already in the --- log-domain. +-- log-domain. Throws errors on @notANumber@ input. logToLogFloat :: (Real a, RealToFrac a Double) => a -> LogFloat {-# SPECIALIZE logToLogFloat :: Double -> LogFloat #-} logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" . realToFrac @@ -310,7 +259,7 @@ fromLogFloat (LogFloat x) = realToFrac (exp x) --- | Return the log-domain value itself without costly conversion +-- | Return the log-domain value itself without conversion. logFromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a) => LogFloat -> a {-# SPECIALIZE logFromLogFloat :: LogFloat -> Double #-} @@ -358,6 +307,56 @@ ---------------------------------------------------------------- +#ifdef __USE_FFI__ +#define LOG1P_WHICH_VERSION specialized 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@. +-- +-- /This installation was compiled to use the LOG1P_WHICH_VERSION/ + +#ifdef __USE_FFI__ +foreign import ccall unsafe "math.h log1p" + log1p :: Double -> Double + +-- Technically we 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. + +#else + +log1p :: Double -> Double +{-# INLINE 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 +-- than doing the naive thing, especially for very small arguments. +-- This function isn't needed internally, but is provided for +-- symmetry with 'log1p'. +-- +-- /This installation was compiled to use the LOG1P_WHICH_VERSION/ + +#ifdef __USE_FFI__ +foreign import ccall unsafe "math.h expm1" + expm1 :: Double -> Double +#else +expm1 :: Double -> Double +{-# INLINE expm1 #-} +expm1 x = exp x - 1 +#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 @@ -376,16 +375,16 @@ (*) (LogFloat x) (LogFloat y) = LogFloat (x+y) - + (+) (LogFloat x) (LogFloat y) - | x >= y = LogFloat (x + log (1 + exp (y - x))) - | otherwise = LogFloat (y + log (1 + exp (x - y))) - + | x >= y = LogFloat (x + log1p (exp (y - x))) + | otherwise = LogFloat (y + log1p (exp (x - y))) + -- Without the guard this would return NaN instead of error (-) (LogFloat x) (LogFloat y) - | x >= y = LogFloat (x + log (1 - exp (y - x))) + | x >= y = LogFloat (x + log1p (negate (exp (y - x)))) | otherwise = errorOutOfRange "(-)" - + signum (LogFloat x) | x == negativeInfinity = 0 | x > negativeInfinity = 1 @@ -417,7 +416,11 @@ -- Just for fun. The more coersion functions the better. Though -- Rationals are very buggy when it comes to transfinite values instance Real LogFloat where - toRational (LogFloat x) = toRational (exp x) + toRational (LogFloat x) + | isInfinite ex || isNaN ex = errorOutOfRange "toRational" + | otherwise = toRational ex + where + 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... 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 @@ -97,10 +97,12 @@ instance RealToFrac Integer Float where -- TODO: is there a more primitive way? + {-# INLINE realToFrac #-} realToFrac j = Prelude.realToFrac j instance RealToFrac Integer Double where -- TODO: is there a more primitive way? + {-# INLINE realToFrac #-} realToFrac j = Prelude.realToFrac j 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,7 +1,14 @@ {-# 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. +-- +-- +{-# OPTIONS_GHC -fglasgow-exts #-} + ---------------------------------------------------------------- --- ~ 2009.01.29 +-- ~ 2009.03.09 -- | -- Module : Data.Number.Transfinite -- Copyright : Copyright (c) 2007--2009 wren ng thornton @@ -91,10 +98,14 @@ -- operations must return @notANumber@, where @inf@ is any value -- which @isInfinite@: -- - -- * @inf + inf@ + -- * @infinity + negativeInfinity@ -- - -- * @inf - inf@ + -- * @negativeInfinity + infinity@ -- + -- * @infinity - infinity@ + -- + -- * @negativeInfinity - negativeInfinity@ + -- -- * @inf * 0@ -- -- * @0 * inf@ @@ -163,6 +174,14 @@ -- for your @0@ and @negativeInfinity@. If it doesn't, then you -- should avoid importing our @log@ and will probably want converters -- to handle the discrepancy. +-- +-- For GHC, this version of @log@ has rules for fusion with @exp@. +-- These can give different behavior by preventing overflow to +-- @infinity@ and preventing errors for taking the logarithm of +-- negative values. For 'Double' and 'Float' they can also give +-- different answers due to eliminating floating point fuzz. The +-- rules strictly improve mathematical accuracy, however they should +-- be noted in case your code depends on the implementation details. log :: (Floating a, Transfinite a) => a -> a {-# SPECIALIZE log :: Double -> Double #-} @@ -180,5 +199,28 @@ -- hack our way around without using PartialOrd by using isNaN, (== -- 0), ((>0).signum) but that would be less efficient. +---------------------------------------------------------------- +-- 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 +-- of our operators for 'LogFloat', with rules like @log x + log y +-- = log (x * y)@ and @log x - log y = log (x / y)@. However the +-- problem is that those equations could be driven in either direction +-- depending on whether we think time performance or non-underflow +-- performance is more important, and the answers may be different +-- at every call site. +-- +-- Since we implore users to do normal-domain computations whenever +-- it would not degenerate accuracy, we should not rewrite their +-- decisions in any way. The log\/exp fusion strictly improves both +-- time and accuracy, so those are safe. But the buck stops with +-- them. ---------------------------------------------------------------- ----------------------------------------------------------- fin. diff --git a/src/Hugs/RealFloat.hs b/src/Hugs/RealFloat.hs --- a/src/Hugs/RealFloat.hs +++ b/src/Hugs/RealFloat.hs @@ -3,8 +3,15 @@ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} +#if defined(__HUGS__) && (__HUGS__ <= 200609) +#define REALFLOAT_VERSION corrected Hugs version. +#elif defined(__GLASGOW_HASKELL__) || defined(__NHC__) +#define REALFLOAT_VERSION normal Prelude version. This should be correct. +#else +#define REALFLOAT_VERSION normal Prelude version. This could be buggy. +#endif ---------------------------------------------------------------- --- ~ 2009.01.29 +-- ~ 2009.03.10 -- | -- Module : Hugs.RealFloat -- Copyright : Copyright (c) 2007--2009 wren ng thornton @@ -25,6 +32,8 @@ -- N.B. The corrected definitions have only been tested to work for -- 'Float' and 'Double'. These definitions should probably not be -- used for other 'RealFloat' types. +-- +-- /This installation was compiled with the REALFLOAT_VERSION/ ---------------------------------------------------------------- module Hugs.RealFloat ( isInfinite @@ -36,6 +45,7 @@ ---------------------------------------------------------------- isInfinite :: (RealFloat a) => a -> Bool +{-# SPECIALIZE isInfinite :: Double -> Bool #-} {-# INLINE isInfinite #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isInfinite x = (1/0) == abs x @@ -45,6 +55,7 @@ isNaN :: (RealFloat a) => a -> Bool +{-# SPECIALIZE isNaN :: Double -> Bool #-} {-# INLINE isNaN #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isNaN x = compareEQ x 0 && compareEQ x 1