add vector Product to Vectors

2 files changed, 84 insertions, 0 deletions

diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 14ba0ff..e8df27c 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -15,6 +15,7 @@
 
 module Numeric.GSL.Vector (
     sumF, sumR, sumQ, sumC,
+    prodF, prodR, prodQ, prodC,
     dotF, dotR, dotQ, dotC,
     FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
     FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
@@ -111,6 +112,39 @@ foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV
 foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV
 foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV
 
+-- | product of elements
+prodF :: Vector Float -> Float
+prodF x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodF vec x vec r "prodF"
+           return $ r @> 0
+
+-- | product of elements
+prodR :: Vector Double -> Double
+prodR x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodR vec x vec r "prodR"
+           return $ r @> 0
+
+-- | product of elements
+prodQ :: Vector (Complex Float) -> Complex Float
+prodQ x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodQ vec x vec r "prodQ"
+           return $ r @> 0
+
+-- | product of elements
+prodC :: Vector (Complex Double) -> Complex Double
+prodC x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodC vec x vec r "prodC"
+           return $ r @> 0
+
+foreign import ccall safe "gsl-aux.h prodF" c_prodF :: TFF
+foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV
+foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV
+foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV
+
 -- | dot product
 dotF :: Vector Float -> Vector Float -> Float
 dotF x y = unsafePerformIO $ do
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 689989d..3e9e2cb 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -153,6 +153,56 @@ int sumC(KCVEC(x),CVEC(r)) {
     OK
 }
 
+int prodF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("prodF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 0;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int prodR(KRVEC(x),RVEC(r)) {
+    DEBUGMSG("prodR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 0;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int prodQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("prodQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex_float res;
+    res.dat[0] = 0;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++) {
+      res.dat[0] = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
+    }
+    rp[0] = res;
+    OK
+}
+    
+int prodC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("prodC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex res;
+    res.dat[0] = 0;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++)  {
+      res.dat[0] = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
+    }
+    rp[0] = res;
+    OK
+}
+
 int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
     DEBUGMSG("dotF");
     REQUIRES(xn==yn,BAD_SIZE);