3 files changed, 276 insertions, 131 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 063e601..7da6f6a 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -3,6 +3,8 @@
 #include <string.h>
 #include <math.h>
 #include <time.h>
+#include <inttypes.h>
+
 #include "lapack-aux.h"
 
 #define MACRO(B) do {B} while (0)
@@ -1300,6 +1302,18 @@ int multiplyI(KOIMAT(a), KOIMAT(b), OIMAT(r)) {
     OK
 }
 
+int multiplyL(KOLMAT(a), KOLMAT(b), OLMAT(r)) {
+    { TRAV(r,i,j) {
+        int k;
+        AT(r,i,j) = 0;
+        for (k=0;k<ac;k++) {
+            AT(r,i,j) += AT(a,i,k) * AT(b,k,j);
+        }
+      }
+    }
+    OK
+}
+
 
 ////////////////// sparse matrix-product ///////////////////////////////////////
 
@@ -1331,67 +1345,23 @@ int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
 
 //////////////////// transpose /////////////////////////
 
-int transF(KFMAT(x),FMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transF");
-    int i,j;
-    for (i=0; i<tr; i++) {
-        for (j=0; j<tc; j++) {
-        tp[i*tc+j] = xp[j*xc+i];
-        }
-    }
-    OK
-}
-
-int transR(KDMAT(x),DMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transR");
-    int i,j;
-    for (i=0; i<tr; i++) {
-        for (j=0; j<tc; j++) {
-        tp[i*tc+j] = xp[j*xc+i];
-        }
-    }
-    OK
-}
-
-int transQ(KQMAT(x),QMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transQ");
-    int i,j;
-    for (i=0; i<tr; i++) {
-        for (j=0; j<tc; j++) {
-        tp[i*tc+j] = xp[j*xc+i];
-        }
-    }
-    OK
-}
-
-int transC(KCMAT(x),CMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transC");
-    int i,j;
-    for (i=0; i<tr; i++) {
-        for (j=0; j<tc; j++) {
-        tp[i*tc+j] = xp[j*xc+i];
-        }
-    }
-    OK
-}
-
-
-int transI(KIMAT(x),IMAT(t)) {
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
-    DEBUGMSG("transI");
-    int i,j;
-    for (i=0; i<tr; i++) {
-        for (j=0; j<tc; j++) {
-            tp[i*tc+j] = xp[j*xc+i];
-        }
-    }
-    OK
-}
-
+#define TRANS_IMP {                       \
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);  \
+    DEBUGMSG("trans");                    \
+    int i,j;                              \
+    for (i=0; i<tr; i++) {                \
+        for (j=0; j<tc; j++) {            \
+            tp[i*tc+j] = xp[j*xc+i];      \
+        }                                 \
+    }                                     \
+    OK }
+
+int transF(KFMAT(x),FMAT(t)) TRANS_IMP
+int transR(KDMAT(x),DMAT(t)) TRANS_IMP
+int transQ(KQMAT(x),QMAT(t)) TRANS_IMP
+int transC(KCMAT(x),CMAT(t)) TRANS_IMP
+int transI(KIMAT(x),IMAT(t)) TRANS_IMP
+int transL(KLMAT(x),LMAT(t)) TRANS_IMP
 
 //////////////////////// extract /////////////////////////////////
 
@@ -1431,6 +1401,11 @@ int extractI(int modei, int modej, KIVEC(i), KIVEC(j), KOIMAT(m), OIMAT(r)) {
     EXTRACT_IMP
 }
 
+int extractL(int modei, int modej, KIVEC(i), KIVEC(j), KOLMAT(m), OLMAT(r)) {
+    EXTRACT_IMP
+}
+
+
 //////////////////////// remap /////////////////////////////////
 
 #define REMAP_IMP                                               \
@@ -1451,6 +1426,10 @@ int remapI(KOIMAT(i), KOIMAT(j), KOIMAT(m), OIMAT(r)) {
     REMAP_IMP
 }
 
+int remapL(KOIMAT(i), KOIMAT(j), KOLMAT(m), OLMAT(r)) {
+    REMAP_IMP
+}
+
 int remapC(KOIMAT(i), KOIMAT(j), KOCMAT(m), OCMAT(r)) {
     REMAP_IMP
 }
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
index 6ffbef1..1549bb5 100644
--- a/packages/base/src/Internal/C/lapack-aux.h
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -37,6 +37,7 @@ typedef short ftnlen;
 /********************************************************/
 
 #define IVEC(A) int A##n, int*A##p
+#define LVEC(A) int A##n, int64_t*A##p
 #define FVEC(A) int A##n, float*A##p
 #define DVEC(A) int A##n, double*A##p
 #define QVEC(A) int A##n, complex*A##p
@@ -44,6 +45,7 @@ typedef short ftnlen;
 #define PVEC(A) int A##n, void* A##p, int A##s
 
 #define IMAT(A) int A##r, int A##c, int* A##p
+#define LMAT(A) int A##r, int A##c, int64_t* A##p
 #define FMAT(A) int A##r, int A##c, float* A##p
 #define DMAT(A) int A##r, int A##c, double* A##p
 #define QMAT(A) int A##r, int A##c, complex* A##p
@@ -51,6 +53,7 @@ typedef short ftnlen;
 #define PMAT(A) int A##r, int A##c, void* A##p, int A##s
 
 #define OIMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, int* A##p
+#define OLMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, int64_t* A##p
 #define OFMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, float* A##p
 #define ODMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, double* A##p
 #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p
@@ -58,6 +61,7 @@ typedef short ftnlen;
 
 
 #define KIVEC(A) int A##n, const int*A##p
+#define KLVEC(A) int A##n, const int64_t*A##p
 #define KFVEC(A) int A##n, const float*A##p
 #define KDVEC(A) int A##n, const double*A##p
 #define KQVEC(A) int A##n, const complex*A##p
@@ -65,6 +69,7 @@ typedef short ftnlen;
 #define KPVEC(A) int A##n, const void* A##p, int A##s
 
 #define KIMAT(A) int A##r, int A##c, const int* A##p
+#define KLMAT(A) int A##r, int A##c, const int64_t* A##p
 #define KFMAT(A) int A##r, int A##c, const float* A##p
 #define KDMAT(A) int A##r, int A##c, const double* A##p
 #define KQMAT(A) int A##r, int A##c, const complex* A##p
@@ -72,6 +77,7 @@ typedef short ftnlen;
 #define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
 
 #define KOIMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const int* A##p
+#define KOLMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const int64_t* A##p
 #define KOFMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const float* A##p
 #define KODMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const double* A##p
 #define KOQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const complex* A##p
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
index fcc9dbd..0921c18 100644
--- a/packages/base/src/Internal/C/vector-aux.c
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -1,4 +1,5 @@
 #include <complex.h>
+#include <inttypes.h>
 
 typedef double complex TCD;
 typedef float  complex TCF;
@@ -66,6 +67,14 @@ int sumI(KIVEC(x),IVEC(r)) {
     OK
 }
 
+int sumL(KLVEC(x),LVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    int res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
 
 int sumQ(KQVEC(x),QVEC(r)) {
     DEBUGMSG("sumQ");
@@ -127,7 +136,14 @@ int prodI(KIVEC(x),IVEC(r)) {
     OK
 }
 
-
+int prodL(KLVEC(x),LVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    int res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
 
 int prodQ(KQVEC(x),QVEC(r)) {
     DEBUGMSG("prodQ");
@@ -189,7 +205,7 @@ double vector_min(KDVEC(x)) {
     return r;
 }
 
-double vector_max_index(KDVEC(x)) {
+int vector_max_index(KDVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
         if(xp[k]>xp[r]) {
@@ -199,7 +215,7 @@ double vector_max_index(KDVEC(x)) {
     return r;
 }
 
-double vector_min_index(KDVEC(x)) {
+int vector_min_index(KDVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
         if(xp[k]<xp[r]) {
@@ -254,7 +270,7 @@ float vector_min_f(KFVEC(x)) {
     return r;
 }
 
-float vector_max_index_f(KFVEC(x)) {
+int vector_max_index_f(KFVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
         if(xp[k]>xp[r]) {
@@ -264,7 +280,7 @@ float vector_max_index_f(KFVEC(x)) {
     return r;
 }
 
-float vector_min_index_f(KFVEC(x)) {
+int vector_min_index_f(KFVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
         if(xp[k]<xp[r]) {
@@ -306,7 +322,7 @@ int vector_max_i(KIVEC(x)) {
 }
 
 int vector_min_i(KIVEC(x)) {
-    float r = xp[0];
+    int r = xp[0];
     int k;
     for (k = 1; k<xn; k++) {
         if(xp[k]<r) {
@@ -352,6 +368,64 @@ int toScalarI(int code, KIVEC(x), IVEC(r)) {
 }
 
 
+int64_t vector_max_l(KLVEC(x)) {
+    int64_t r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+int64_t vector_min_l(KLVEC(x)) {
+    int64_t r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+int vector_max_index_l(KLVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+int vector_min_index_l(KLVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+
+int toScalarL(int code, KLVEC(x), LVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int64_t res;
+    switch(code) {
+        case 2: { res = vector_max_index_l(V(x));  break; }
+        case 3: { res = vector_max_l(V(x));  break; }
+        case 4: { res = vector_min_index_l(V(x)); break; }
+        case 5: { res = vector_min_l(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
 double dznrm2_(integer*, const doublecomplex*, integer*);
 double dzasum_(integer*, const doublecomplex*, integer*);
 
@@ -477,6 +551,17 @@ int mapI(int code, KIVEC(x), IVEC(r)) {
 }
 
 
+int mapL(int code, KLVEC(x), LVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    switch (code) {
+        OP(3,abs)
+        OP(15,sign)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
 
 inline double abs_complex(doublecomplex z) {
     return sqrt(z.r*z.r + z.i*z.i);
@@ -641,6 +726,32 @@ int mapValI(int code, int* pval, KIVEC(x), IVEC(r)) {
 }
 
 
+int64_t mod_l (int64_t a, int64_t b) {
+    int64_t m = a % b;
+    if (b>0) {
+        return m >=0 ? m : m+b;
+    } else {
+        return m <=0 ? m : m+b;
+    }
+}
+
+int mapValL(int code, int64_t* pval, KLVEC(x), LVEC(r)) {
+    int k;
+    int64_t val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValL");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(6,mod(val,xp[k]))
+        OPV(7,mod(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
 
 inline doublecomplex complex_add(doublecomplex a, doublecomplex b) {
     doublecomplex r;
@@ -736,6 +847,19 @@ REQUIRES(an == bn && an == rn, BAD_SIZE);
 }
 
 
+int zipL(int code, KLVEC(a), KLVEC(b), LVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipI Add",+)
+        OPZO(1,"zipI Sub",-)
+        OPZO(2,"zipI Mul",*)
+        OPZO(3,"zipI Div",/)
+        OPZO(6,"zipI Mod",%)
+        default: ERROR(BAD_CODE);
+    }
+}
+
 
 #define OPZOb(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = a2p[k] O b2p[k]; OK }
 #define OPZEb(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = E(a2p[k],b2p[k]); OK }
@@ -990,6 +1114,18 @@ int sort_valuesI(KIVEC(v),IVEC(r)) {
     OK
 }
 
+int
+compare_longs(const void *a, const void *b) {
+  return *(int64_t*)a > *(int64_t*)b;
+}
+
+int sort_valuesL(KLVEC(v),LVEC(r)) {
+    memcpy(rp,vp,vn*sizeof(int64_t));
+    qsort(rp,rn,sizeof(int64_t),compare_ints);
+    OK
+}
+
+
 ////////////////////////////////////////
 
 
@@ -1010,7 +1146,7 @@ int sort_valuesI(KIVEC(v),IVEC(r)) {
     OK
 
 
-typedef struct SDI { int pos; double val;} DI;
+typedef struct DI { int pos; double val;} DI;
 
 int compare_doubles_i (const void *a, const void *b) {
   return ((DI*)a)->val > ((DI*)b)->val;
@@ -1043,6 +1179,17 @@ int sort_indexI(KIVEC(v),IVEC(r)) {
 }
 
 
+typedef struct LI { int pos; int64_t val;} LI;
+
+int compare_longs_i (const void *a, const void *b) {
+  return ((II*)a)->val > ((II*)b)->val;
+}
+
+int sort_indexL(KLVEC(v),LVEC(r)) {
+    SORTIDX_IMP(II,compare_longs_i)
+}
+
+
 ////////////////////////////////////////////////////////////////////////////////
 
 int round_vector(KDVEC(v),DVEC(r)) {
@@ -1088,6 +1235,44 @@ int range_vector(IVEC(r)) {
     OK
 }
 
+///////////////////////////
+
+
+int round_vector_l(KDVEC(v),LVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = round(vp[k]);
+    }
+    OK
+}
+
+
+int mod_vector_l(int64_t m, KLVEC(v), LVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] % m;
+    }
+    OK
+}
+
+int div_vector_l(int64_t m, KLVEC(v), LVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] / m;
+    }
+    OK
+}
+
+int range_vector_l(LVEC(r)) {
+    int k;
+    for(k=0; k<rn; k++) {
+        rp[k] = k;
+    }
+    OK
+}
+
+
+
 //////////////////// constant /////////////////////////
 
 int constantF(float * pval, FVEC(r)) {
@@ -1143,65 +1328,46 @@ int constantI(int * pval, IVEC(r)) {
 }
 
 
-//////////////////// float-double conversion /////////////////////////
 
-int float2double(FVEC(x),DVEC(y)) {
-    DEBUGMSG("float2double")
+int constantL(int64_t * pval, LVEC(r)) {
+    DEBUGMSG("constantL")
     int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
+    int64_t val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
     }
     OK
 }
 
-int float2int(KFVEC(x),IVEC(y)) {
-    DEBUGMSG("float2int")
-    int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
-    }
-    OK
-}
 
+//////////////////// type conversions /////////////////////////
 
-int double2float(DVEC(x),FVEC(y)) {
-    DEBUGMSG("double2float")
-    int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
-    }
-    OK
-}
+#define CONVERT_IMP {     \
+    int k;                \
+    for(k=0;k<xn;k++) {   \
+        yp[k]=xp[k];      \
+    }                     \
+    OK }
 
+int float2double(FVEC(x),DVEC(y)) CONVERT_IMP
 
-int double2int(KDVEC(x),IVEC(y)) {
-    DEBUGMSG("double2int")
-    int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
-    }
-    OK
-}
+int float2int(KFVEC(x),IVEC(y)) CONVERT_IMP
 
+int double2float(DVEC(x),FVEC(y)) CONVERT_IMP
 
-int int2float(KIVEC(x),FVEC(y)) {
-    DEBUGMSG("int2float")
-    int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
-    }
-    OK
-}
+int double2int(KDVEC(x),IVEC(y)) CONVERT_IMP
 
+int double2long(KDVEC(x),LVEC(y)) CONVERT_IMP
 
-int int2double(KIVEC(x),DVEC(y)) {
-    DEBUGMSG("int2double")
-    int k;
-    for(k=0;k<xn;k++) {
-        yp[k]=xp[k];
-    }
-    OK
-}
+int int2float(KIVEC(x),FVEC(y)) CONVERT_IMP
+
+int int2double(KIVEC(x),DVEC(y)) CONVERT_IMP
+
+int int2long(KIVEC(x),LVEC(y)) CONVERT_IMP
+
+int lont2int(KLVEC(x),IVEC(y)) CONVERT_IMP
+
+int long2double(KLVEC(x),DVEC(y)) CONVERT_IMP
 
 
 //////////////////// conjugate /////////////////////////
@@ -1249,6 +1415,11 @@ int stepI(KIVEC(x),IVEC(y)) {
     STEP_IMP
 }
 
+int stepL(KLVEC(x),LVEC(y)) {
+    STEP_IMP
+}
+
+
 //////////////////// cond /////////////////////////
 
 #define COMPARE_IMP                               \
@@ -1272,26 +1443,10 @@ int compareI(KIVEC(x),KIVEC(y),IVEC(r)) {
     COMPARE_IMP
 }
 
-
-#define COND_IMP                                                            \
-    REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);  \
-    int k;                                                                  \
-    for(k=0;k<xn;k++) {                                                     \
-        rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);             \
-    }                                                                       \
-    OK
-
-int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
-    COND_IMP
-}
-
-int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
-    COND_IMP
+int compareL(KLVEC(x),KLVEC(y),IVEC(r)) {
+    COMPARE_IMP
 }
 
-int condI(KIVEC(x),KIVEC(y),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
-    COND_IMP
-}
 
 
 #define CHOOSE_IMP                                                      \
@@ -1314,6 +1469,11 @@ int chooseI(KIVEC(cond),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
     CHOOSE_IMP
 }
 
+int chooseL(KIVEC(cond),KLVEC(lt),KLVEC(eq),KLVEC(gt),LVEC(r)) {
+    CHOOSE_IMP
+}
+
+
 int chooseC(KIVEC(cond),KCVEC(lt),KCVEC(eq),KCVEC(gt),CVEC(r)) {
     CHOOSE_IMP
 }