modular C matrix product

2 files changed, 41 insertions, 29 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 7da6f6a..1601bef 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1290,29 +1290,25 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
 }
 
 
-int multiplyI(KOIMAT(a), KOIMAT(b), OIMAT(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);
-        }
-      }
+#define MULT_IMP_VER(OP)        \
+    { TRAV(r,i,j) {             \
+        int k;                  \
+        AT(r,i,j) = 0;          \
+        for (k=0;k<ac;k++) {    \
+            OP                  \
+        }                       \
+      }                         \
     }
-    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
-}
+#define MULT_IMP {                                                                   \
+    if (m==1) {                                                                      \
+        MULT_IMP_VER( AT(r,i,j) += AT(a,i,k) * AT(b,k,j); )                          \
+    } else {                                                                         \
+        MULT_IMP_VER( AT(r,i,j) = (AT(r,i,j) + (AT(a,i,k) * AT(b,k,j)) % m) % m ; )  \
+    } OK }
+
+int multiplyI(int     m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP
+int multiplyL(int32_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP
 
 
 ////////////////// sparse matrix-product ///////////////////////////////////////
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
index 70e46bc..580aa1c 100644
--- a/packages/base/src/Internal/C/vector-aux.c
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -58,20 +58,28 @@ int sumR(KDVEC(x),DVEC(r)) {
     OK
 }
 
-int sumI(KIVEC(x),IVEC(r)) {
+int sumI(int m, KIVEC(x),IVEC(r)) {
     REQUIRES(rn==1,BAD_SIZE);
     int i;
     int res = 0;
-    for (i = 0; i < xn; i++) res += xp[i];
+    if (m==1) {
+        for (i = 0; i < xn; i++) res += xp[i];
+    } else {
+        for (i = 0; i < xn; i++) res = (res + xp[i]) % m;
+    }
     rp[0] = res;
     OK
 }
 
-int sumL(KLVEC(x),LVEC(r)) {
+int sumL(int32_t m, KLVEC(x),LVEC(r)) {
     REQUIRES(rn==1,BAD_SIZE);
     int i;
     int res = 0;
-    for (i = 0; i < xn; i++) res += xp[i];
+    if (m==1) {
+        for (i = 0; i < xn; i++) res += xp[i];
+    } else {
+        for (i = 0; i < xn; i++) res = (res + xp[i]) % m;
+    }
     rp[0] = res;
     OK
 }
@@ -127,20 +135,28 @@ int prodR(KDVEC(x),DVEC(r)) {
     OK
 }
 
-int prodI(KIVEC(x),IVEC(r)) {
+int prodI(int m, KIVEC(x),IVEC(r)) {
     REQUIRES(rn==1,BAD_SIZE);
     int i;
     int res = 1;
-    for (i = 0; i < xn; i++) res *= xp[i];
+    if (m==1) {
+        for (i = 0; i < xn; i++) res *= xp[i];
+    } else {
+        for (i = 0; i < xn; i++) res = (res * xp[i]) % m;
+    }
     rp[0] = res;
     OK
 }
 
-int prodL(KLVEC(x),LVEC(r)) {
+int prodL(int32_t m, KLVEC(x),LVEC(r)) {
     REQUIRES(rn==1,BAD_SIZE);
     int i;
     int res = 1;
-    for (i = 0; i < xn; i++) res *= xp[i];
+    if (m==1) {
+        for (i = 0; i < xn; i++) res *= xp[i];
+    } else {
+        for (i = 0; i < xn; i++) res = (res * xp[i]) % m;
+    }
     rp[0] = res;
     OK
 }