ieee-utils 0.3 → 0.4.0

raw patch · 2 files changed

+109/−2 lines, 2 files

Files

• ieee-utils.cabal +3/−2
• src/Numeric/IEEE/Monad.hs +106/−0

@@ -1,5 +1,5 @@ name:               ieee-utils-version:            0.3+version:            0.4.0 cabal-version:      >= 1.3 build-type:         Simple license:            BSD3@@ -22,4 +22,5 @@   c-sources:        cbits/ieee_utils.c   extra-libraries:   exposed-modules:  Numeric.IEEE.RoundMode,-                    Numeric.IEEE.FloatExceptions+                    Numeric.IEEE.FloatExceptions,+                    Numeric.IEEE.Monad

@@ -0,0 +1,106 @@+{-# LANGUAGE ImplicitParams, RankNTypes, ExistentialQuantification, ScopedTypeVariables, GeneralizedNewtypeDeriving #-}++{- |+  Module      :  Numeric.IEEE.Monad+  Copyright   :  (c) Sterling Clover 2008 <s.clover@gmail.com>+  License     :  BSD3+  Maintainer  :  Matt Morrow <mjm2002@gmail.com>+  Stability   :  provisional+  Portability :  portable (FFI)++  The IEEE monad provides tools for enforcing sequencing of calculations such that+  fine grained control is provided over triggering exceptions, evaluations within+  particular rounding modes, etc. The perturb family of functions is built using this,+  allowing pure computations paramaterized over an arbitrary Floating type to be+  tested for numeric stability.++-}++module Numeric.IEEE.Monad where++import Numeric.IEEE.RoundMode (RoundMode(..))+import qualified Numeric.IEEE.RoundMode as RM+import qualified Numeric.IEEE.FloatExceptions as FE+import Control.Applicative+import Control.Concurrent+import Control.Exception+import Control.Monad++-- | All uses of the IEEE monad need to be wrapped inside a top level call to+-- withIeeeDo. This ensures that access to floating point internals+-- is serialized properly, even between multiple threads.+withIeeeDo :: ((?ieeeMutex :: MVar ()) => IO a) -> IO a+withIeeeDo f = newMVar () >>= \x -> let ?ieeeMutex = x in f++newtype IEEE a = IEEE {unIEEE :: IO a} deriving (Monad, Functor, Applicative)++runIEEE :: (?ieeeMutex :: MVar()) => IEEE a -> IO a+runIEEE f = do+  takeMVar ?ieeeMutex+  ret <- unIEEE f+  putMVar ?ieeeMutex ()+  return ret++getRound :: IEEE RoundMode+getRound = IEEE RM.getRound++setRound :: RoundMode -> IEEE Bool+setRound m = IEEE $ RM.setRound m++clearFloatExcepts :: [ArithException] -> IEEE Bool+clearFloatExcepts xs = IEEE $ FE.clearFloatExcepts xs++getFloatExcepts :: IEEE [ArithException]+getFloatExcepts = IEEE FE.getFloatExcepts++-- | Forces strict evaluation of the enclosed numeric argument.+calculate :: a -> IEEE a+calculate = IEEE . evaluate++-- | Calculate, but also returns any floating exceptions triggered.+calculate' :: a -> IEEE (a,[FE.ArithException])+calculate' f = do+  getFloatExcepts >>= clearFloatExcepts+  ret <- IEEE . evaluate $ f+  exs <- getFloatExcepts+  return (ret,exs)++-- | Executes the specified IEEE action within a specific round mode.+withRoundMode :: RoundMode -> IEEE a -> IEEE a+withRoundMode r f = do+  x <- getRound+  if x == r then f else setRound r >> (f >>= calculate) >>= \ret -> setRound x >> return ret++-- | Given something of (forall a. Floating a => IEEE a) produces a four-tuple of+-- the value as calculated rounding up, down, to nearest, and towards zero.+perturb' :: (?ieeeMutex :: MVar (), Floating b) => (forall a. Floating a => IEEE a) -> IO (b, b, b, b)+perturb' f = runIEEE $ do+  x <- getRound+  setRound Upward+  u <- calculate =<< foo f+  setRound Downward+  d <- calculate =<< foo f+  setRound ToNearest+  tn <- calculate =<< foo f+  setRound TowardZero+  tz <- calculate =<< foo f+  setRound x+  return (u,d,tn,tz)+  where foo :: Floating b => (forall c. Floating c => IEEE c) -> IEEE b+        foo = id++-- | Given something that produces a Floating, returns a representation of the+-- absolute difference between the results as calculated rounding upwards and downwards.+perturb :: (?ieeeMutex :: MVar (), Floating b) => (forall a. Floating a => IEEE a) -> IO b+perturb f = do+  (u,d,_,_) <- perturb' f+  return $ abs (u - d)++-- | Given something that produces a Floating, returns the magnitude of instability+-- introduced by perturbing the equation by rounding upwards and then downwards. This is+-- the absolute difference between the results as calculated rounding upwards and downwards,+-- and then divided by the averaged result.+perturbedMag :: (?ieeeMutex :: MVar (), Floating b) => (forall a. Floating a => IEEE a) -> IO b+perturbedMag f = do+  (u,d,_,_) <- perturb' f+  return $ abs ((u - d) / (u + d))