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