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     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.Core
  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 "set_roundmode"
  74     set_roundmode :: CInt -> IO CInt
  75 foreign import ccall unsafe "get_roundmode"
  76     get_roundmode :: IO CInt
  77
  78 foreign import ccall unsafe "fegetenv"
  79     c_fegetenv :: Ptr FEnvState -> IO CInt
  80 foreign import ccall unsafe "feholdexcept"
  81     c_feholdexcept :: Ptr FEnvState -> IO CInt
  82 foreign import ccall unsafe "fenv_restore"
  83     fenv_restore :: Ptr FEnvState -> Ptr CUInt -> IO CInt
  84 foreign import ccall unsafe "fenv_raise_excepts"
  85     fenv_raise_excepts :: CUInt -> IO CInt
  86
  87 data RoundingMode = ToNearest | Upward | Downward | TowardZero
  88     deriving (Show, Read, Enum, Bounded)
  89 data FloatException = DivByZero | Inexact | Invalid | Overflow | Underflow
  90     deriving (Show, Read, Enum, Bounded)
  91
  92 -- | Opaque type which stores the complete floating point environment.  It
  93 -- corresponds to the C type @fenv_t@.
  94 newtype FEnvState = FEnvState (ForeignPtr FEnvState)
  95
  96 instance Storable FEnvState where
  97     sizeOf    = const SIZEOF_FENV_T
  98     alignment = const ALIGNOF_FENV_T
  99
 100     peek ptr = do
 101         fp <- mallocForeignPtrBytes SIZEOF_FENV_T
 102         withForeignPtr fp (\p -> copyBytes p ptr SIZEOF_FENV_T)
 103         return (FEnvState fp)
 104     poke ptr (FEnvState fp) = do
 105         withForeignPtr fp (\p -> copyBytes ptr p SIZEOF_FENV_T)
 106
 107 -- | Container for computations which will be run in a modified floating point
 108 -- environment.  The FEnv container records all operations for later evaluation
 109 -- by 'fenvEval'.  Note that 'pure' is strict in order to force evaluation
 110 -- of floating point values stored in the container.
 111 --
 112 -- Do not use the 'Eq' or 'Show' instances, they are provided only because Num
 113 -- requires them.
 114 data FEnv a = forall b . FEnv (b -> a) !b
 115
 116 -- In the following instances, the two FEnv parts must be bashed together
 117 -- exactly once every time the contained value is extracted.  Care must be
 118 -- taken to avoid memoization of this result.  Interestingly, FEnv is not an
 119 -- instance of Monad: While join (FEnv f x) = f x has the right type, it does
 120 -- not satisfy this important property.
 121
 122 instance Functor FEnv where
 123     fmap f (FEnv g x) = FEnv (f . g) x
 124
 125 instance Applicative FEnv where
 126     pure = FEnv id
 127     (FEnv f x) <*> (FEnv g y) = FEnv (\(x',y') -> f x' . g $ y') (x, y)
 128
 129 -- For hysterical raisins, we need to instance Eq and Show since they are
 130 -- superclasses of Num.
 131 instance Eq a => Eq (FEnv a) where
 132     (==) = error "The Eq instance for FEnv is a lie."
 133 instance Show a => Show (FEnv a) where
 134     show = const "<<FEnv>>"
 135
 136 instance Num a => Num (FEnv a) where
 137     (+)         = liftA2 (+)
 138     (-)         = liftA2 (-)
 139     (*)         = liftA2 (*)
 140     negate      = liftA negate
 141     signum      = liftA signum
 142     abs         = liftA abs
 143     fromInteger = pure . fromInteger
 144
 145 instance Fractional a => Fractional (FEnv a) where
 146     (/)          = liftA2 (/)
 147     recip        = liftA recip
 148     fromRational = pure . fromRational
 149
 150 instance Floating a => Floating (FEnv a) where
 151     (**)  = liftA2 (**)
 152     sqrt  = liftA  sqrt
 153     acos  = liftA  acos
 154     asin  = liftA  asin
 155     atan  = liftA  atan
 156     cos   = liftA  cos
 157     sin   = liftA  sin
 158     tan   = liftA  tan
 159     cosh  = liftA  cosh
 160     sinh  = liftA  sinh
 161     tanh  = liftA  tanh
 162     exp   = liftA  exp
 163     log   = liftA  log
 164     acosh = liftA  acosh
 165     asinh = liftA  asinh
 166     atanh = liftA  atanh
 167
 168 instance RealFloat a => RealFloat (FEnv a) where
 169     fma       = liftA3 fma
 170     copysign  = liftA2 copysign
 171     nextafter = liftA2 nextafter
 172     fmod      = liftA2 fmod
 173     frem      = liftA2 frem
 174     atan2     = liftA2 atan2
 175     hypot     = liftA2 hypot
 176     cbrt      = liftA  cbrt
 177     exp2      = liftA  exp2
 178     expm1     = liftA  expm1
 179     log10     = liftA  log10
 180     log1p     = liftA  log1p
 181     log2      = liftA  log2
 182     erf       = liftA  erf
 183     erfc      = liftA  erfc
 184     gamma     = liftA  gamma
 185     lgamma    = liftA  lgamma
 186     nearbyint = liftA  nearbyint
 187     rint      = liftA  rint
 188
 189     infinity  = pure infinity
 190     nan       = pure nan
 191     pi        = pure pi
 192
 193 -- | Saves the current floating point environment and, optionally, clears all
 194 -- floating point exception flags and sets non-stop (continue on exceptions)
 195 -- mode.
 196 unsafeSaveEnvironment :: Bool -> IO FEnvState
 197 unsafeSaveEnvironment reset = alloca $ \env -> do
 198     rc <- saveEnv env
 199     unless (rc == 0) $ fail "Error saving floating point environment."
 200     peek env
 201     where
 202         saveEnv = if reset then c_feholdexcept else c_fegetenv
 203
 204 -- | Restores a previously-saved floating point environment and returns the
 205 -- list of floating point exceptions that occurred prior to restoring the
 206 -- environment.
 207 unsafeRestoreEnvironment :: FEnvState -> IO [FloatException]
 208 unsafeRestoreEnvironment (FEnvState fp) = alloca $ \pe -> do
 209     rc <- withForeignPtr fp (flip fenv_restore pe)
 210     unless (rc == 0) $ fail "Error restoring floating point environment."
 211     rawExcepts <- peek pe
 212     return $! filter (testBit rawExcepts . fromEnum) [minBound..maxBound]
 213
 214 -- | Raises the given floating point exceptions.
 215 unsafeRaiseExceptions :: [FloatException] -> IO ()
 216 unsafeRaiseExceptions ex = do
 217     rc <- fenv_raise_excepts $ foldr (flip setBit . fromEnum) 0 ex
 218     unless (rc == 0) $ fail "Error raising floating point exceptions."
 219
 220 unsafeSetRoundingMode :: RoundingMode -> IO ()
 221 unsafeSetRoundingMode mode = do
 222     rc <- set_roundmode (fromIntegral (fromEnum mode))
 223     unless (rc == 0) $ fail "Error setting rounding mode"
 224
 225 getRoundingMode :: IO RoundingMode
 226 getRoundingMode = do
 227     rc <- get_roundmode
 228     unless (rc >= 0) $ fail "Error getting rounding mode"
 229     return . toEnum . fromIntegral $ rc
 230
 231 -- | Evaluate an FEnv using a specific rounding mode.  Rounding mode selections
 232 -- nest: subcomputations might use another mode.  The default rounding mode is
 233 -- unspecified.
 234 withRoundingMode :: RoundingMode -> FEnv a -> FEnv a
 235 withRoundingMode mode (FEnv f x) = FEnv unsafePerformIO $ do
 236     oldMode <- getRoundingMode
 237     unsafeSetRoundingMode mode
 238     rc <- evaluate $ f x
 239     unsafeSetRoundingMode oldMode
 240     return rc
 241
 242 -- | Raise floating point exceptions as part of an FEnv computation.
 243 raiseExceptions :: [FloatException] -> FEnv a -> FEnv a
 244 raiseExceptions ex = liftA2 seq $
 245     FEnv unsafePerformIO (unsafeRaiseExceptions ex)
 246
 247 -- | Save and restore the floating point state across an FEnv computation.
 248 -- The result includes the exceptions returned during this computation.
 249 holdExceptions :: FEnv a -> FEnv (a, [FloatException])
 250 holdExceptions = FEnv unsafePerformIO . fenvEval
 251
 252 -- | Same as 'holdExceptions', except that any exceptions raised during the
 253 -- computation are discarded.
 254 holdExceptions_ :: FEnv a -> FEnv a
 255 holdExceptions_ = fmap fst . holdExceptions
 256
 257 -- | This function is to help with debugging the floating point environment
 258 -- handling.  @fenvTrace msg x@ constructs an FEnv value containing @x@ that
 259 -- prints @msg@ (using 'Debug.Trace.trace') whenever the value is extracted.
 260 fenvTrace :: String -> a -> FEnv a
 261 fenvTrace s = fmap (trace s) . pure
 262
 263 -- | Runs all the computations which are recorded in an FEnv container.  The
 264 -- floating point environment is preserved across this call, and any floating
 265 -- point exceptions which were raised during the computation are returned.
 266 fenvEval :: FEnv a -> IO (a, [FloatException])
 267 fenvEval (FEnv f x) = do
 268     env <- unsafeSaveEnvironment True
 269     rc  <- evaluate $ f x
 270     ex  <- unsafeRestoreEnvironment env
 271     return (rc, ex)