{-# LANGUAGE ForeignFunctionInterface #-}
module Numerical.Integration
(integration, IntegralResult(..))
where
import Foreign.Marshal.Alloc (free, mallocBytes)
import Foreign.Ptr (FunPtr, Ptr, freeHaskellFunPtr)
import Foreign.Storable (peek, sizeOf)
import Foreign.C (CDouble(..), CInt(..))
zeroDouble :: Double
zeroDouble = 0.0
-- zeroCDouble :: CDouble
-- zeroCDouble = 0.0
nanDouble :: Double
nanDouble = zeroDouble / zeroDouble
-- nanCDouble :: CDouble
-- nanCDouble = zeroCDouble / zeroCDouble
double2Cdouble :: Double -> CDouble
double2Cdouble = realToFrac
-- cdouble2double :: CDouble -> Double
-- cdouble2double = realToFrac
data IntegralResult = IntegralResult {
_value :: Double,
_error :: Double,
_code :: Int
} deriving Show
foreign import ccall safe "wrapper" funPtr
:: (CDouble -> CDouble) -> IO (FunPtr (CDouble -> CDouble))
foreign import ccall safe "integration" c_integration
:: FunPtr (CDouble -> CDouble) -> CDouble -> CDouble -> CDouble -> CInt
-> Ptr CDouble -> Ptr CInt -> IO CDouble
-- | Numerical integration.
integration :: (CDouble -> CDouble) -- ^ integrand
-> Double -- ^ lower bound
-> Double -- ^ upper bound
-> Double -- ^ desired relative error
-> Int -- ^ number of subdivisions
-> IO IntegralResult -- ^ value, error estimate, error code
integration f lower upper relError subdiv = do
errorEstimatePtr <- mallocBytes (sizeOf (0 :: CDouble))
errorCodePtr <- mallocBytes (sizeOf (0 :: CInt))
fPtr <- funPtr f
let lower' = double2Cdouble lower
upper' = double2Cdouble upper
relError' = double2Cdouble relError
subdiv' = fromIntegral subdiv
result <-
c_integration fPtr lower' upper' relError' subdiv' errorEstimatePtr errorCodePtr
let result' = if isNaN result
then nanDouble
else realToFrac result
errorEstimate <- peek errorEstimatePtr
let errorEstimate' = if isNaN errorEstimate
then nanDouble
else realToFrac errorEstimate
errorCode <- fromIntegral <$> peek errorCodePtr
let out = IntegralResult {_value = result', _error = errorEstimate', _code = errorCode}
free errorEstimatePtr
free errorCodePtr
freeHaskellFunPtr fPtr
return out