{-# 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 g = g

-- | 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))