Data/Floating/Environment.hs

   1 {-
   2  - Copyright (C) 2010 Nick Bowler.
   3  -
   4  - License BSD2:  2-clause BSD license.  See LICENSE for full terms.
   5  - This is free software: you are free to change and redistribute it.
   6  - There is NO WARRANTY, to the extent permitted by law.
   7  -}
   8
   9 -- | Access to the floating point environment.  Performing this access within
  10 -- a Haskell program turns out to be extremely problematic, because floating
  11 -- point operations are secretly impure.  For example, the innocent-looking
  12 -- function
  13 --
  14 -- @  (+) :: Double -> Double -> Double@
  15 --
  16 -- potentially both depends on and modifies the global floating point
  17 -- environment.
  18 --
  19 -- This module avoids the referential transparency problems that occur as a
  20 -- result of accessing the floating point environment by restricting when
  21 -- computations which access it are evaluated.  There is some minor discipline
  22 -- required of the programmer: she must arrange her code so that pure floating
  23 -- point expressions are not forced during a call to 'fenvEval'.
  24 -- See @fenv-impure.hs@ in the @examples/@ directory of the altfloat
  25 -- distribution for why this discipline is necessary.
  26 --
  27 -- FEnv instances the numeric classes, so it should be possible to use
  28 -- natural syntax.  Note that the operations done on FEnv are stored so that
  29 -- they can be performed later, thus one should be take care not to construct
  30 -- huge thunks when using this interface.
  31 --
  32 -- This interface has not been tested in multi-threaded programs.  It might
  33 -- work: more info is needed about GHC's threading support.
  34 {-# LANGUAGE CPP, ForeignFunctionInterface, ExistentialQuantification #-}
  35 {-# OPTIONS_GHC -I. #-}
  36 module Data.Floating.Environment (
  37     module Control.Applicative,
  38
  39     -- * Data types
  40     RoundingMode(..), FloatException(..), FEnvState, FEnv,
  41
  42     -- * Controlled access to the floating point environment
  43     -- | These functions can still break referential transparency, because it
  44     -- is possible to arrange for a pure floating point expression to be forced
  45     -- during the execution of 'fenvEval'.  The easiest way to ensure that this
  46     -- does not happen is to only use such expressions as the argument to
  47     -- 'pure'; never as the argument to 'fmap'.
  48     fenvEval, withRoundingMode, raiseExceptions, holdExceptions,
  49     holdExceptions_, fenvTrace,
  50
  51     -- * Direct access to the floating point environment
  52     -- | Special care must be taken when using these functions.  Modifying the
  53     -- floating point environment will affect all floating point computations
  54     -- that have not yet been evaluated.
  55     unsafeSaveEnvironment, unsafeRestoreEnvironment,
  56     unsafeTestExceptions,  unsafeRaiseExceptions,
  57     unsafeSetRoundingMode, getRoundingMode
  58 ) where
  59
  60 #include <config.h>
  61
  62 import Prelude hiding (Float, Double, Floating(..), RealFloat(..))
  63
  64 import Data.Floating.Types
  65 import Control.Exception
  66 import Control.Applicative
  67 import Control.Monad
  68 import Debug.Trace
  69
  70 import Foreign.C
  71 import Foreign
  72
  73 foreign import ccall unsafe "fenv_set_roundmode"
  74     fenv_set_roundmode :: CInt -> IO CInt
  75 foreign import ccall unsafe "fenv_get_roundmode"
  76     fenv_get_roundmode :: IO CInt
  77 foreign import ccall unsafe "fenv_test_excepts"
  78     fenv_test_excepts :: IO CUInt
  79 foreign import ccall unsafe "fenv_raise_excepts"
  80     fenv_raise_excepts :: CUInt -> IO CInt
  81
  82 foreign import ccall unsafe "fegetenv"
  83     c_fegetenv :: Ptr FEnvState -> IO CInt
  84 foreign import ccall unsafe "fesetenv"
  85     c_fesetenv :: Ptr FEnvState -> IO CInt
  86 foreign import ccall unsafe "feholdexcept"
  87     c_feholdexcept :: Ptr FEnvState -> IO CInt
  88
  89 data RoundingMode = ToNearest | Upward | Downward | TowardZero
  90     deriving (Show, Read, Enum, Bounded)
  91 data FloatException = DivByZero | Inexact | Invalid | Overflow | Underflow
  92     deriving (Show, Read, Enum, Bounded)
  93
  94 -- | Opaque type which stores the complete floating point environment.  It
  95 -- corresponds to the C type @fenv_t@.
  96 newtype FEnvState = FEnvState (ForeignPtr FEnvState)
  97
  98 instance Storable FEnvState where
  99     sizeOf    = const SIZEOF_FENV_T
 100     alignment = const ALIGNOF_FENV_T
 101
 102     peek ptr = do
 103         fp <- mallocForeignPtrBytes SIZEOF_FENV_T
 104         withForeignPtr fp (\p -> copyBytes p ptr SIZEOF_FENV_T)
 105         return (FEnvState fp)
 106     poke ptr (FEnvState fp) = do
 107         withForeignPtr fp (\p -> copyBytes ptr p SIZEOF_FENV_T)
 108
 109 -- | Container for computations which will be run in a modified floating point
 110 -- environment.  The FEnv container records all operations for later evaluation
 111 -- by 'fenvEval'.  Note that 'pure' is strict in order to force evaluation
 112 -- of floating point values stored in the container.
 113 --
 114 -- Do not use the 'Eq' or 'Show' instances, they are provided only because Num
 115 -- requires them.
 116 data FEnv a = forall b . FEnv (b -> a) !b
 117
 118 -- In the following instances, the two FEnv parts must be bashed together
 119 -- exactly once every time the contained value is extracted.  Care must be
 120 -- taken to avoid memoization of this result.  Interestingly, FEnv is not an
 121 -- instance of Monad: While join (FEnv f x) = f x has the right type, it does
 122 -- not satisfy this important property.
 123
 124 instance Functor FEnv where
 125     fmap f (FEnv g x) = FEnv (f . g) x
 126
 127 instance Applicative FEnv where
 128     pure = FEnv id
 129     (FEnv f x) <*> (FEnv g y) = FEnv (\(x',y') -> f x' . g $ y') (x, y)
 130
 131 -- For hysterical raisins, we need to instance Eq and Show since they are
 132 -- superclasses of Num.
 133 instance Eq a => Eq (FEnv a) where
 134     (==) = error "The Eq instance for FEnv is a lie."
 135 instance Show a => Show (FEnv a) where
 136     show = const "<<FEnv>>"
 137
 138 instance Num a => Num (FEnv a) where
 139     (+)         = liftA2 (+)
 140     (-)         = liftA2 (-)
 141     (*)         = liftA2 (*)
 142     negate      = liftA negate
 143     signum      = liftA signum
 144     abs         = liftA abs
 145     fromInteger = pure . fromInteger
 146
 147 instance Fractional a => Fractional (FEnv a) where
 148     (/)          = liftA2 (/)
 149     recip        = liftA recip
 150     fromRational = pure . fromRational
 151
 152 instance Floating a => Floating (FEnv a) where
 153     (**)  = liftA2 (**)
 154     sqrt  = liftA  sqrt
 155     acos  = liftA  acos
 156     asin  = liftA  asin
 157     atan  = liftA  atan
 158     cos   = liftA  cos
 159     sin   = liftA  sin
 160     tan   = liftA  tan
 161     cosh  = liftA  cosh
 162     sinh  = liftA  sinh
 163     tanh  = liftA  tanh
 164     exp   = liftA  exp
 165     log   = liftA  log
 166     acosh = liftA  acosh
 167     asinh = liftA  asinh
 168     atanh = liftA  atanh
 169
 170 instance RealFloat a => RealFloat (FEnv a) where
 171     fma       = liftA3 fma
 172     copysign  = liftA2 copysign
 173     nextafter = liftA2 nextafter
 174     fmod      = liftA2 fmod
 175     frem      = liftA2 frem
 176     atan2     = liftA2 atan2
 177     hypot     = liftA2 hypot
 178     cbrt      = liftA  cbrt
 179     exp2      = liftA  exp2
 180     expm1     = liftA  expm1
 181     log10     = liftA  log10
 182     log1p     = liftA  log1p
 183     log2      = liftA  log2
 184     erf       = liftA  erf
 185     erfc      = liftA  erfc
 186     gamma     = liftA  gamma
 187     lgamma    = liftA  lgamma
 188     nearbyint = liftA  nearbyint
 189     rint      = liftA  rint
 190
 191     infinity  = pure infinity
 192     nan       = pure nan
 193     pi        = pure pi
 194
 195 -- | Saves the current floating point environment and, optionally, clears all
 196 -- floating point exception flags and sets non-stop (continue on exceptions)
 197 -- mode.
 198 unsafeSaveEnvironment :: Bool -> IO FEnvState
 199 unsafeSaveEnvironment reset = alloca $ \env -> do
 200     rc <- saveEnv env
 201     unless (rc == 0) $ fail "Error saving floating point environment."
 202     peek env
 203     where
 204         saveEnv = if reset then c_feholdexcept else c_fegetenv
 205
 206 -- | Restores a previously-saved floating point environment and returns the
 207 -- list of floating point exceptions that occurred prior to restoring the
 208 -- environment.
 209 unsafeRestoreEnvironment :: FEnvState -> IO [FloatException]
 210 unsafeRestoreEnvironment (FEnvState fp) = do
 211     excepts <- unsafeTestExceptions
 212     rc <- withForeignPtr fp c_fesetenv
 213     unless (rc == 0) $ fail "Error restoring floating point environment."
 214     return excepts
 215
 216 -- | Return the currently raised floating point exceptions as a list.
 217 unsafeTestExceptions :: IO [FloatException]
 218 unsafeTestExceptions = do
 219     rawExcepts <- fenv_test_excepts
 220     return $! filter (testBit rawExcepts . fromEnum) [minBound..maxBound]
 221
 222 -- | Raises the given floating point exceptions.
 223 unsafeRaiseExceptions :: [FloatException] -> IO ()
 224 unsafeRaiseExceptions ex = do
 225     rc <- fenv_raise_excepts $ foldr (flip setBit . fromEnum) 0 ex
 226     unless (rc == 0) $ fail "Error raising floating point exceptions."
 227
 228 unsafeSetRoundingMode :: RoundingMode -> IO ()
 229 unsafeSetRoundingMode mode = do
 230     rc <- fenv_set_roundmode (fromIntegral (fromEnum mode))
 231     unless (rc == 0) $ fail "Error setting rounding mode"
 232
 233 getRoundingMode :: IO RoundingMode
 234 getRoundingMode = do
 235     rc <- fenv_get_roundmode
 236     unless (rc >= 0) $ fail "Error getting rounding mode"
 237     return . toEnum . fromIntegral $ rc
 238
 239 -- | Evaluate an FEnv using a specific rounding mode.  Rounding mode selections
 240 -- nest: subcomputations might use another mode.  The default rounding mode is
 241 -- unspecified.
 242 withRoundingMode :: RoundingMode -> FEnv a -> FEnv a
 243 withRoundingMode mode (FEnv f x) = FEnv unsafePerformIO $ do
 244     oldMode <- getRoundingMode
 245     unsafeSetRoundingMode mode
 246     rc <- evaluate $ f x
 247     unsafeSetRoundingMode oldMode
 248     return rc
 249
 250 -- | Raise floating point exceptions as part of an FEnv computation.
 251 raiseExceptions :: [FloatException] -> FEnv a -> FEnv a
 252 raiseExceptions ex = liftA2 seq $
 253     FEnv unsafePerformIO (unsafeRaiseExceptions ex)
 254
 255 -- | Save and restore the floating point state across an FEnv computation.
 256 -- The result includes the exceptions returned during this computation.
 257 holdExceptions :: FEnv a -> FEnv (a, [FloatException])
 258 holdExceptions = FEnv unsafePerformIO . fenvEval
 259
 260 -- | Same as 'holdExceptions', except that any exceptions raised during the
 261 -- computation are discarded.
 262 holdExceptions_ :: FEnv a -> FEnv a
 263 holdExceptions_ = fmap fst . holdExceptions
 264
 265 -- | This function is to help with debugging the floating point environment
 266 -- handling.  @fenvTrace msg x@ constructs an FEnv value containing @x@ that
 267 -- prints @msg@ (using 'Debug.Trace.trace') whenever the value is extracted.
 268 fenvTrace :: String -> a -> FEnv a
 269 fenvTrace s = fmap (trace s) . pure
 270
 271 -- | Runs all the computations which are recorded in an FEnv container.  The
 272 -- floating point environment is preserved across this call, and any floating
 273 -- point exceptions which were raised during the computation are returned.
 274 fenvEval :: FEnv a -> IO (a, [FloatException])
 275 fenvEval (FEnv f x) = do
 276     env <- unsafeSaveEnvironment True
 277     rc  <- evaluate $ f x
 278     ex  <- unsafeRestoreEnvironment env
 279     return (rc, ex)