5 files changed, 287 insertions, 1 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index dddd269..acefe92 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -82,3 +82,7 @@ isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
 
 scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
 scast = fromJust . cast
+
+{- | conversion of Haskell functions into function pointers that can be used in the C side
+-}
+foreign import ccall "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
diff --git a/lib/GSL/Differentiation.hs b/lib/GSL/Differentiation.hs
new file mode 100644
index 0000000..6fd7565
--- /dev/null
+++ b/lib/GSL/Differentiation.hs
@@ -0,0 +1,79 @@
+{-# OPTIONS  #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  GSL.Differentiation
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Numerical differentiation.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>
+
+From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
+-}
+-----------------------------------------------------------------------------
+module GSL.Differentiation (
+    derivCentral,
+    derivForward,
+    derivBackward
+) where
+
+import Foreign
+import Data.Packed.Internal.Common(mkfun,check,(//))
+
+derivGen ::
+    Int                   -- ^ type: 0 central, 1 forward, 2 backward
+    -> Double             -- ^ initial step size
+    -> (Double -> Double) -- ^ function
+    -> Double             -- ^ point where the derivative is taken
+    -> (Double, Double)   -- ^ result and error
+derivGen c h f x = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x) 
+    c_deriv c fp x h r e // check "deriv" []
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall "gsl-aux.h deriv" 
+ c_deriv :: Int -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double 
+                    -> Ptr Double -> Ptr Double -> IO Int
+
+
+{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
+
+> > let deriv = derivCentral 0.01 
+> > deriv sin (pi/4)
+>(0.7071067812000676,1.0600063101654055e-10)
+> > cos (pi/4)
+>0.7071067811865476 
+
+-}
+derivCentral :: Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivCentral = derivGen 0
+
+-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
+derivForward :: Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivForward = derivGen 1
+
+-- | Adaptive backward difference algorithm, /gsl_deriv_backward/. 
+derivBackward ::Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivBackward = derivGen 2
diff --git a/lib/GSL/Integration.hs b/lib/GSL/Integration.hs
new file mode 100644
index 0000000..dd76be8
--- /dev/null
+++ b/lib/GSL/Integration.hs
@@ -0,0 +1,90 @@
+{-# OPTIONS #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  GSL.Integration
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Numerical integration routines.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
+-}
+-----------------------------------------------------------------------------
+
+module GSL.Integration (
+    integrateQNG,
+    integrateQAGS
+) where
+
+import Foreign
+import Data.Packed.Internal.Common(mkfun,check,(//))
+
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
+
+@\> let quad = integrateQAGS 1E-9 1000 
+\> let f a x = x**(-0.5) * log (a*x)
+\> quad (f 1) 0 1
+(-3.999999999999974,4.871658632055187e-13)@
+ 
+-}
+
+integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
+                 -> Double           -- ^ a
+                 -> Double           -- ^ b
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQAGS prec n f a b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x) 
+    c_integrate_qags fp a b prec n r e // check "integrate_qags" []
+    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_qags" 
+ c_integrate_qags :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> Int 
+                     -> Ptr Double -> Ptr Double -> IO Int
+
+-----------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
+
+@\> let quad = integrateQNG 1E-6 
+\> quad (\\x -> 4\/(1+x*x)) 0 1 
+(3.141592653589793,3.487868498008632e-14)@
+ 
+-}
+
+integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
+                 -> Double           -- ^ a
+                 -> Double           -- ^ b
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQNG prec f a b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x) 
+    c_integrate_qng fp a b prec r e  // check "integrate_qng" []
+    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_qng" 
+ c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double 
+                    -> Ptr Double -> Ptr Double -> IO Int
+
diff --git a/lib/GSL/Special.hs b/lib/GSL/Special.hs
new file mode 100644
index 0000000..643c499
--- /dev/null
+++ b/lib/GSL/Special.hs
@@ -0,0 +1,101 @@
+{-# OPTIONS #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  GSL.Special
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Special functions.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Special-Functions.html#Special-Functions>
+-}
+-----------------------------------------------------------------------------
+
+module GSL.Special (
+    erf,
+    erf_Z,
+    bessel_J0_e,
+    exp_e10_e
+)
+where
+
+import Foreign
+import Data.Packed.Internal.Common(check,(//))
+
+----------------------------------------------------------------
+-------------------- simple functions --------------------------
+
+{- | The error function (/gsl_sf_erf/), defined as 2\/ \\sqrt \\pi * \int\_0\^t \\exp -t\^2 dt.
+
+@> erf 1.5
+0.9661051464753108@
+
+-}
+foreign import ccall "gsl-aux.h gsl_sf_erf" erf :: Double -> Double
+
+{- | The Gaussian probability density function (/gsl_sf_erf_Z/), defined as (1\/\\sqrt\{2\\pi\}) \\exp(-x\^2\/2).
+
+>> erf_Z 1.5
+>0.12951759566589172
+
+-}
+foreign import ccall "gsl-aux.h gsl_sf_erf_Z" erf_Z :: Double -> Double
+
+----------------------------------------------------------------
+-- the sf_result struct is equivalent to an array of two doubles
+
+createSFR s f = unsafePerformIO $ do
+    p <- mallocArray 2
+    f p // check "createSFR" []
+    [val,err] <- peekArray 2 p
+    free p
+    return (val,err)
+
+-------------------- functions returning sf_result -------------
+
+{- | The regular cylindrical Bessel function of zeroth order, J_0(x). This is
+   the example in the GSL manual, returning the value of the function and
+   the error term: 
+
+@\> bessel_J0_e 5.0
+(-0.1775967713143383,1.9302109579684196e-16)@
+
+-}
+bessel_J0_e :: Double -> (Double,Double)
+bessel_J0_e x = createSFR "bessel_J0_e" (gsl_sf_bessel_J0_e x)
+foreign import ccall "gsl-aux.h gsl_sf_bessel_J0_e" gsl_sf_bessel_J0_e :: Double -> Ptr Double -> IO Int
+
+---------------------------------------------------------------------
+-- the sf_result_e10 contains two doubles and the exponent
+
+createSFR_E10 s f = unsafePerformIO $ do
+    let sd = sizeOf (0::Double)
+    let si = sizeOf (0::Int)
+    p <- mallocBytes (2*sd + si)
+    f p // check "createSFR_E10" []
+    val <- peekByteOff p 0
+    err <- peekByteOff p sd
+    expo <- peekByteOff p (2*sd) 
+    free p
+    return (val,expo,err)
+
+-------------------- functions returning sf_result_e10 -------------
+
+{- | (From the GSL manual) \"This function computes the exponential \exp(x) using the @gsl_sf_result_e10@ type to return a result with extended range. This function may be useful if the value of \exp(x) would overflow the numeric range of double\". 
+
+For example:
+
+@\> exp_e10_e 30.0
+(1.0686474581524432,1.4711818964275088e-14,13)@
+
+@\> exp 30.0
+1.0686474581524463e13@
+
+-}
+exp_e10_e :: Double -> (Double,Int,Double)
+exp_e10_e x = createSFR_E10 "exp_e10_e" (gsl_sf_exp_e10_e x)
+foreign import ccall "gsl-aux.h gsl_sf_exp_e10_e" gsl_sf_exp_e10_e :: Double -> Ptr Double -> IO Int
diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs
index e538b7e..3225139 100644
--- a/lib/GSL/Vector.hs
+++ b/lib/GSL/Vector.hs
@@ -18,7 +18,7 @@ module GSL.Vector (
     FunCodeV(..), vectorMapR, vectorMapC,
     FunCodeSV(..), vectorMapValR, vectorMapValC,
     FunCodeVV(..), vectorZipR, vectorZipC,
-    scale, addConstant
+    scale, addConstant, add, mul,
 ) where
 
 import Data.Packed.Internal
@@ -79,6 +79,18 @@ addConstant x v | isReal baseOf v = scast $ vectorMapValR AddConstant (scast x)
                 | isComp baseOf v = scast $ vectorMapValC AddConstant (scast x) (scast v)
                 | otherwise   = fromList $ map (*x) $ toList v
 
+add :: (Num a, Field a) => Vector a -> Vector a -> Vector a
+add u v | isReal baseOf v = scast $ vectorZipR Add (scast u) (scast v)
+                | isComp baseOf v = scast $ vectorZipC Add (scast u) (scast v)
+                | otherwise   = fromList $ zipWith (+) (toList u) (toList v)
+
+mul :: (Num a, Field a) => Vector a -> Vector a -> Vector a
+mul u v | isReal baseOf v = scast $ vectorZipR Mul (scast u) (scast v)
+                | isComp baseOf v = scast $ vectorZipC Mul (scast u) (scast v)
+                | otherwise   = fromList $ zipWith (*) (toList u) (toList v)
+
+
+
 ------------------------------------------------------------------
 
 toScalarAux fun code v = unsafePerformIO $ do