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.Classes
  65 import Control.Exception
  66 import Control.Applicative
  67 import Control.Monad
  68
  69 import System.IO.Unsafe
  70 import Debug.Trace
  71
  72 import Foreign.C
  73 import Foreign
  74
  75 foreign import ccall unsafe "set_roundmode"
  76     set_roundmode :: CInt -> IO CInt
  77 foreign import ccall unsafe "get_roundmode"
  78     get_roundmode :: IO CInt
  79
  80 foreign import ccall unsafe "fegetenv"
  81     c_fegetenv :: Ptr FEnvState -> IO CInt
  82 foreign import ccall unsafe "feholdexcept"
  83     c_feholdexcept :: Ptr FEnvState -> IO CInt
  84 foreign import ccall unsafe "fenv_restore"
  85     fenv_restore :: Ptr FEnvState -> Ptr CUInt -> IO CInt
  86 foreign import ccall unsafe "fenv_raise_excepts"
  87     fenv_raise_excepts :: CUInt -> 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) = alloca $ \pe -> do
 211     rc <- withForeignPtr fp (flip fenv_restore pe)
 212     unless (rc == 0) $ fail "Error restoring floating point environment."
 213     rawExcepts <- peek pe
 214     return $! filter (testBit rawExcepts . fromEnum) [minBound..maxBound]
 215
 216 -- | Raises the given floating point exceptions.
 217 unsafeRaiseExceptions :: [FloatException] -> IO ()
 218 unsafeRaiseExceptions ex = do
 219     rc <- fenv_raise_excepts $ foldr (flip setBit . fromEnum) 0 ex
 220     unless (rc == 0) $ fail "Error raising floating point exceptions."
 221
 222 unsafeSetRoundingMode :: RoundingMode -> IO ()
 223 unsafeSetRoundingMode mode = do
 224     rc <- set_roundmode (fromIntegral (fromEnum mode))
 225     unless (rc == 0) $ fail "Error setting rounding mode"
 226
 227 getRoundingMode :: IO RoundingMode
 228 getRoundingMode = do
 229     rc <- get_roundmode
 230     unless (rc >= 0) $ fail "Error getting rounding mode"
 231     return . toEnum . fromIntegral $ rc
 232
 233 -- | Evaluate an FEnv using a specific rounding mode.  Rounding mode selections
 234 -- nest: subcomputations might use another mode.  The default rounding mode is
 235 -- unspecified.
 236 withRoundingMode :: RoundingMode -> FEnv a -> FEnv a
 237 withRoundingMode mode (FEnv f x) = FEnv unsafePerformIO $ do
 238     oldMode <- getRoundingMode
 239     unsafeSetRoundingMode mode
 240     rc <- evaluate $ f x
 241     unsafeSetRoundingMode oldMode
 242     return rc
 243
 244 -- | Raise floating point exceptions as part of an FEnv computation.
 245 raiseExceptions :: [FloatException] -> FEnv a -> FEnv a
 246 raiseExceptions ex = liftA2 seq $
 247     FEnv unsafePerformIO (unsafeRaiseExceptions ex)
 248
 249 -- | Save and restore the floating point state across an FEnv computation.
 250 -- The result includes the exceptions returned during this computation.
 251 holdExceptions :: FEnv a -> FEnv (a, [FloatException])
 252 holdExceptions = FEnv unsafePerformIO . fenvEval
 253
 254 -- | Same as 'holdExceptions', except that any exceptions raised during the
 255 -- computation are discarded.
 256 holdExceptions_ :: FEnv a -> FEnv a
 257 holdExceptions_ = fmap fst . holdExceptions
 258
 259 -- | This function is to help with debugging the floating point environment
 260 -- handling.  @fenvTrace msg x@ constructs an FEnv value containing @x@ that
 261 -- prints @msg@ (using 'Debug.Trace.trace') whenever the value is extracted.
 262 fenvTrace :: String -> a -> FEnv a
 263 fenvTrace s = fmap (trace s) . pure
 264
 265 -- | Runs all the computations which are recorded in an FEnv container.  The
 266 -- floating point environment is preserved across this call, and any floating
 267 -- point exceptions which were raised during the computation are returned.
 268 fenvEval :: FEnv a -> IO (a, [FloatException])
 269 fenvEval (FEnv f x) = do
 270     env <- unsafeSaveEnvironment True
 271     rc  <- evaluate $ f x
 272     ex  <- unsafeRestoreEnvironment env
 273     return (rc, ex)