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