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