From 1319c4ddae89c41633a4cd894ba28896d3445e79 Mon Sep 17 00:00:00 2001 From: Danny Chan Date: Mon, 6 Feb 2012 21:24:10 +0100 Subject: add support for gsl_integration_qagi, gsl_integration_qagiu, gsl_integration_qagil --- lib/Numeric/GSL/Integration.hs | 104 ++++++++++++++++++++++++++++++++++++++++- lib/Numeric/GSL/gsl-aux.c | 42 +++++++++++++++++ 2 files changed, 145 insertions(+), 1 deletion(-) (limited to 'lib/Numeric') diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs index a0e922b..2330fc6 100644 --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs @@ -17,7 +17,10 @@ Numerical integration routines. module Numeric.GSL.Integration ( integrateQNG, - integrateQAGS + integrateQAGS, + integrateQAGI, + integrateQAGIU, + integrateQAGIL ) where import Foreign.C.Types @@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO $ do foreign import ccall "gsl-aux.h integrate_qng" c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). +For example: + +@\> let quad = integrateQAGI 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) +(2.5066282746310002,6.229215880648858e-11)@ + +-} + +integrateQAGI :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf) + -> (Double, Double) -- ^ result of the integration and error +integrateQAGI prec n f = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagi" + c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). +For example: + +@\> let quad = integrateQAGIU 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ a + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIU prec n f a = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagiu" + c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). +For example: + +@\> let quad = integrateQAGIL 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIL prec n f b = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagil" + c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 33d7dab..5c17836 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -761,6 +761,48 @@ int integrate_qags(double f(double,void*), double a, double b, double prec, int OK } +int integrate_qagi(double f(double,void*), double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagi"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagiu(double f(double,void*), double a, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagiu"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagil(double f(double,void*), double b, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagil"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + int polySolve(KRVEC(a), CVEC(z)) { DEBUGMSG("polySolve"); REQUIRES(an>1,BAD_SIZE); -- cgit v1.2.3