Data/Floating/Environment.hs

   1 {-
   2  - Copyright (C) 2009 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.  Floating point exceptions are not
  29 -- handled yet, so anything in this documentation which refers to \"the
  30 -- floating point environment\" actually means \"the current rounding mode\".
  31 --
  32 -- Multi-threaded programs are almost certainly /not/ supported at this time.
  33 {-# LANGUAGE ForeignFunctionInterface, ExistentialQuantification #-}
  34 {-# INCLUDE "cfloat.h" #-}
  35 module Data.Floating.Environment (
  36     RoundingMode(..), FEnv,
  37     fenvEval, withRoundingMode, fenvTrace,
  38     unsafeSetRoundingMode, getRoundingMode
  39 ) where
  40
  41 import Prelude hiding (Float, Double, Floating(..), RealFloat(..))
  42
  43 import Data.Floating.Classes
  44 import Control.Exception
  45 import Control.Applicative
  46 import Control.Monad
  47
  48 import System.IO.Unsafe
  49 import Debug.Trace
  50
  51 import Foreign.C
  52
  53 foreign import ccall unsafe "set_roundmode"
  54     set_roundmode :: CInt -> IO CInt
  55 foreign import ccall unsafe "get_roundmode"
  56     get_roundmode :: IO CInt
  57
  58 data RoundingMode = ToNearest | Upward | Downward | TowardZero
  59     deriving (Show, Read, Enum, Bounded)
  60
  61 -- | Container for computations which will be run in a modified floating point
  62 -- environment.  The FEnv container records all operations for later evaluation
  63 -- by 'fenvEval'.  Do not use the 'Eq' or 'Show' instances, they are provided
  64 -- only because Num requires them.
  65 data FEnv a = forall b . FEnv (b -> a) b
  66
  67 -- In the following instances, the two FEnv parts must be bashed together
  68 -- exactly once every time the contained value is extracted.  Care must be
  69 -- taken to avoid memoization of this result.  Interestingly, FEnv is not an
  70 -- instance of Monad: While join (FEnv f x) = f x has the right type, it does
  71 -- not satisfy this important property.
  72
  73 instance Functor FEnv where
  74     fmap f (FEnv g x) = FEnv (f . g) x
  75
  76 instance Applicative FEnv where
  77     pure x = seq x (FEnv id x)
  78     (FEnv f x) <*> (FEnv g y) = FEnv (\(x',y') -> f x' . g $ y') (x, y)
  79
  80 -- For hysterical raisins, we need to instance Eq and Show since they are
  81 -- superclasses of Num.
  82 instance Eq a => Eq (FEnv a) where
  83     (==) = error "The Eq instance for FEnv is a lie."
  84 instance Show a => Show (FEnv a) where
  85     show = const "<<FEnv>>"
  86
  87 instance Num a => Num (FEnv a) where
  88     (+)         = liftA2 (+)
  89     (-)         = liftA2 (-)
  90     (*)         = liftA2 (*)
  91     negate      = liftA negate
  92     signum      = liftA signum
  93     abs         = liftA abs
  94     fromInteger = pure . fromInteger
  95
  96 instance Fractional a => Fractional (FEnv a) where
  97     (/)          = liftA2 (+)
  98     recip        = liftA recip
  99     fromRational = pure . fromRational
 100
 101 instance Floating a => Floating (FEnv a) where
 102     (**)  = liftA2 (**)
 103     sqrt  = liftA  sqrt
 104     acos  = liftA  acos
 105     asin  = liftA  asin
 106     atan  = liftA  atan
 107     cos   = liftA  cos
 108     sin   = liftA  sin
 109     tan   = liftA  tan
 110     cosh  = liftA  cosh
 111     sinh  = liftA  sinh
 112     tanh  = liftA  tanh
 113     exp   = liftA  exp
 114     log   = liftA  log
 115     acosh = liftA  acosh
 116     asinh = liftA  asinh
 117     atanh = liftA  atanh
 118
 119 instance RealFloat a => RealFloat (FEnv a) where
 120     fma       = liftA3 fma
 121     copysign  = liftA2 copysign
 122     nextafter = liftA2 nextafter
 123     fmod      = liftA2 fmod
 124     frem      = liftA2 frem
 125     atan2     = liftA2 atan2
 126     hypot     = liftA2 hypot
 127     cbrt      = liftA  cbrt
 128     exp2      = liftA  exp2
 129     expm1     = liftA  expm1
 130     log10     = liftA  log10
 131     log1p     = liftA  log1p
 132     log2      = liftA  log2
 133     logb      = liftA  logb
 134     erf       = liftA  erf
 135     erfc      = liftA  erfc
 136     lgamma    = liftA  lgamma
 137     tgamma    = liftA  tgamma
 138
 139     classify  = error "classify is not supported on FEnv."
 140     fquotRem  = error "fquotRem is not supported on FEnv."
 141
 142 -- | Sets the floating point rounding mode.  This function is considered unsafe
 143 -- because it affects unevaluated thunks, breaking referential transparency.
 144 unsafeSetRoundingMode :: RoundingMode -> IO ()
 145 unsafeSetRoundingMode mode = do
 146     rc <- set_roundmode (fromIntegral (fromEnum mode))
 147     unless (rc == 0) $ fail "Error setting rounding mode"
 148
 149 -- | Gets the current floating point rounding mode.
 150 getRoundingMode :: IO RoundingMode
 151 getRoundingMode = do
 152     rc <- get_roundmode
 153     unless (rc >= 0) $ fail "Error getting rounding mode"
 154     return . toEnum . fromIntegral $ rc
 155
 156 -- | Evaluate an FEnv using a specific rounding mode.  Rounding mode selections
 157 -- nest: subcomputations might use another mode.  The default rounding mode is
 158 -- unspecified.
 159 withRoundingMode :: RoundingMode -> FEnv a -> FEnv a
 160 withRoundingMode mode (FEnv f x) = FEnv unsafePerformIO $ do
 161     oldMode <- getRoundingMode
 162     unsafeSetRoundingMode mode
 163     rc <- evaluate $ f x
 164     unsafeSetRoundingMode oldMode
 165     return rc
 166
 167 -- | This function is to help with debugging the floating point environment
 168 -- handling.  @fenvTrace msg x@ constructs an FEnv value containing @x@ that
 169 -- prints @msg@ (using 'Debug.Trace.trace') whenever the value is extracted.
 170 fenvTrace :: String -> a -> FEnv a
 171 fenvTrace s = fmap (trace s) . pure
 172
 173 -- | Runs all the computations which are recorded in an FEnv container.  The
 174 -- floating point environment is preserved across this call.
 175 fenvEval :: FEnv a -> IO a
 176 fenvEval (FEnv f x) = evaluate $ f x