Linear and Floating (Complex Float)

2 files changed, 147 insertions, 3 deletions

diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 0148c4f..14ba0ff 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -17,9 +17,9 @@ module Numeric.GSL.Vector (
     sumF, sumR, sumQ, sumC,
     dotF, dotR, dotQ, dotC,
     FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
-    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF,
-    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF,
-    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF,
+    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
+    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
+    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
     RandDist(..), randomVector
 ) where
 
@@ -214,6 +214,12 @@ vectorMapF = vectorMapAux c_vectorMapF
 
 foreign import ccall safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF
 
+-- | map of real vectors with given function
+vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float)
+vectorMapQ = vectorMapAux c_vectorMapQ
+
+foreign import ccall safe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV
+
 -------------------------------------------------------------------
 
 -- | map of real vectors with given function
@@ -234,6 +240,12 @@ vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)
 
 foreign import ccall safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF
 
+-- | map of complex vectors with given function
+vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)
+vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)
+
+foreign import ccall safe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV
+
 -------------------------------------------------------------------
 
 -- | elementwise operation on real vectors
@@ -254,6 +266,12 @@ vectorZipF = vectorZipAux c_vectorZipF
 
 foreign import ccall safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF
 
+-- | elementwise operation on complex vectors
+vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)
+vectorZipQ = vectorZipAux c_vectorZipQ
+
+foreign import ccall safe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV
+
 -----------------------------------------------------------------------
 
 data RandDist = Uniform  -- ^ uniform distribution in [0,1)
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 3f9eeba..689989d 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -393,6 +393,83 @@ int mapC(int code, KCVEC(x), CVEC(r)) {
 }
 
 
+gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
+{
+  gsl_complex c;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c.dat[0] = a.dat[0];
+  c.dat[1] = a.dat[1];
+
+  r = (*cf)(c);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), 
+                                        gsl_complex_float a,gsl_complex_float b)
+{
+  gsl_complex c1;
+  gsl_complex c2;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c1.dat[0] = a.dat[0];
+  c1.dat[1] = a.dat[1];
+
+  c2.dat[0] = b.dat[0];
+  c2.dat[1] = b.dat[1];
+
+  r = (*cf)(c1,c2);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
+#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
+int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapQ");
+    switch (code) {
+        OPC(0,gsl_complex_sin)
+        OPC(1,gsl_complex_cos)
+        OPC(2,gsl_complex_tan)
+        OPC(3,complex_abs)
+        OPC(4,gsl_complex_arcsin)
+        OPC(5,gsl_complex_arccos)
+        OPC(6,gsl_complex_arctan)
+        OPC(7,gsl_complex_sinh)
+        OPC(8,gsl_complex_cosh)
+        OPC(9,gsl_complex_tanh)
+        OPC(10,gsl_complex_arcsinh)
+        OPC(11,gsl_complex_arccosh)
+        OPC(12,gsl_complex_arctanh)
+        OPC(13,gsl_complex_exp)
+        OPC(14,gsl_complex_log)
+        OPC(15,complex_signum)
+        OPC(16,gsl_complex_sqrt)
+
+        // gsl_complex_arg
+        // gsl_complex_abs
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapQ(int code, KQVEC(x), QVEC(r)) {
+    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
 int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
     int k;
     double val = *pval;
@@ -446,6 +523,27 @@ int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
 }
 
 
+int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
+    int k;
+    gsl_complex_float val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValQ");
+    switch (code) {
+        OPCA(0,gsl_complex_mul,val,xp[k])
+        OPCA(1,gsl_complex_div,val,xp[k])
+        OPCA(2,gsl_complex_add,val,xp[k])
+        OPCA(3,gsl_complex_sub,val,xp[k])
+        OPCA(4,gsl_complex_pow,val,xp[k])
+        OPCA(5,gsl_complex_pow,xp[k],val)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
+    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
 #define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
 #define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
 int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
@@ -519,6 +617,34 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
 }
 
 
+#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
+int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPCZE(0,"zipQ Add",gsl_complex_add)
+        OPCZE(1,"zipQ Sub",gsl_complex_sub)
+        OPCZE(2,"zipQ Mul",gsl_complex_mul)
+        OPCZE(3,"zipQ Div",gsl_complex_div)
+        OPCZE(4,"zipQ Pow",gsl_complex_pow)
+        //OPZE(5,"zipR ATan2",atan2)
+    }
+    //KCVVIEW(a);
+    //KCVVIEW(b);
+    //CVVIEW(r);
+    //gsl_vector_memcpy(V(r),V(a));
+    //int res;
+    switch(code) {
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
+    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
+}
+
+
 
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);