summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/GSL/Integration.hs104
-rw-r--r--lib/Numeric/GSL/gsl-aux.c42
2 files changed, 145 insertions, 1 deletions
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.
17 17
18module Numeric.GSL.Integration ( 18module Numeric.GSL.Integration (
19 integrateQNG, 19 integrateQNG,
20 integrateQAGS 20 integrateQAGS,
21 integrateQAGI,
22 integrateQAGIU,
23 integrateQAGIL
21) where 24) where
22 25
23import Foreign.C.Types 26import Foreign.C.Types
@@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO $ do
94foreign import ccall "gsl-aux.h integrate_qng" 97foreign import ccall "gsl-aux.h integrate_qng"
95 c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double 98 c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double
96 -> Ptr Double -> Ptr Double -> IO CInt 99 -> Ptr Double -> Ptr Double -> IO CInt
100
101--------------------------------------------------------------------
102{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
103For example:
104
105@\> let quad = integrateQAGI 1E-9 1000
106\> let f a x = exp(-a * x^2)
107\> quad (f 0.5)
108(2.5066282746310002,6.229215880648858e-11)@
109
110-}
111
112integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
113 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
114 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
115 -> (Double, Double) -- ^ result of the integration and error
116integrateQAGI prec n f = unsafePerformIO $ do
117 r <- malloc
118 e <- malloc
119 fp <- mkfun (\x _ -> f x)
120 c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi"
121 vr <- peek r
122 ve <- peek e
123 let result = (vr,ve)
124 free r
125 free e
126 freeHaskellFunPtr fp
127 return result
128
129foreign import ccall "gsl-aux.h integrate_qagi"
130 c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt
131 -> Ptr Double -> Ptr Double -> IO CInt
132
133--------------------------------------------------------------------
134{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
135For example:
136
137@\> let quad = integrateQAGIU 1E-9 1000
138\> let f a x = exp(-a * x^2)
139\> quad (f 0.5) 0
140(1.2533141373155001,3.114607940324429e-11)@
141
142-}
143
144integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
145 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
146 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
147 -> Double -- ^ a
148 -> (Double, Double) -- ^ result of the integration and error
149integrateQAGIU prec n f a = unsafePerformIO $ do
150 r <- malloc
151 e <- malloc
152 fp <- mkfun (\x _ -> f x)
153 c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu"
154 vr <- peek r
155 ve <- peek e
156 let result = (vr,ve)
157 free r
158 free e
159 freeHaskellFunPtr fp
160 return result
161
162foreign import ccall "gsl-aux.h integrate_qagiu"
163 c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt
164 -> Ptr Double -> Ptr Double -> IO CInt
165
166--------------------------------------------------------------------
167{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
168For example:
169
170@\> let quad = integrateQAGIL 1E-9 1000
171\> let f a x = exp(-a * x^2)
172\> quad (f 0.5) 0
173(1.2533141373155001,3.114607940324429e-11)@
174
175-}
176
177integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
178 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
179 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
180 -> Double -- ^ b
181 -> (Double, Double) -- ^ result of the integration and error
182integrateQAGIL prec n f b = unsafePerformIO $ do
183 r <- malloc
184 e <- malloc
185 fp <- mkfun (\x _ -> f x)
186 c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil"
187 vr <- peek r
188 ve <- peek e
189 let result = (vr,ve)
190 free r
191 free e
192 freeHaskellFunPtr fp
193 return result
194
195foreign import ccall "gsl-aux.h integrate_qagil"
196 c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt
197 -> Ptr Double -> Ptr Double -> IO CInt
198
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
761 OK 761 OK
762} 762}
763 763
764int integrate_qagi(double f(double,void*), double prec, int w,
765 double *result, double* error) {
766 DEBUGMSG("integrate_qagi");
767 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
768 gsl_function F;
769 F.function = f;
770 F.params = NULL;
771 int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error);
772 CHECK(res,res);
773 gsl_integration_workspace_free (wk);
774 OK
775}
776
777
778int integrate_qagiu(double f(double,void*), double a, double prec, int w,
779 double *result, double* error) {
780 DEBUGMSG("integrate_qagiu");
781 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
782 gsl_function F;
783 F.function = f;
784 F.params = NULL;
785 int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error);
786 CHECK(res,res);
787 gsl_integration_workspace_free (wk);
788 OK
789}
790
791
792int integrate_qagil(double f(double,void*), double b, double prec, int w,
793 double *result, double* error) {
794 DEBUGMSG("integrate_qagil");
795 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
796 gsl_function F;
797 F.function = f;
798 F.params = NULL;
799 int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error);
800 CHECK(res,res);
801 gsl_integration_workspace_free (wk);
802 OK
803}
804
805
764int polySolve(KRVEC(a), CVEC(z)) { 806int polySolve(KRVEC(a), CVEC(z)) {
765 DEBUGMSG("polySolve"); 807 DEBUGMSG("polySolve");
766 REQUIRES(an>1,BAD_SIZE); 808 REQUIRES(an>1,BAD_SIZE);