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