various matrix multiplication procedures

3 files changed, 38 insertions, 4 deletions

diff --git a/test/monniaux/int_mat.c b/test/monniaux/int_mat.c
index 936f45d6..5bea64fc 100644
--- a/test/monniaux/int_mat.c
+++ b/test/monniaux/int_mat.c
@@ -38,6 +38,26 @@ void modint_mat_mul2(unsigned m, unsigned n, unsigned p,
   }
 }
 
+void modint_mat_mul3(unsigned m, unsigned n, unsigned p,
+		     modint * restrict c, unsigned stride_c,
+		     const modint *a, unsigned stride_a,
+		     const modint *b, unsigned stride_b) {
+  for(unsigned i=0; i<m; i++) {
+    for(unsigned k=0; k<p; k++) {
+      modint total0 = 0, total1 = 0;
+      unsigned j;
+      for(j=0; j<n; j+=2) {
+	total0 += a[i*stride_a + j] * b[j*stride_b + k];
+	total1 += a[i*stride_a + (j+1)] * b[(j+1)*stride_b + k];
+      }
+      if (j < n) {
+	total0 += a[i*stride_a + j] * b[j*stride_b + k];
+      }
+      c[i*stride_c+k] = (total0+total1) % MODULUS;
+    }
+  }
+}
+
 modint modint_random(void) {
   static uint64_t next = 1325997111;
   next = next * 1103515245 + 12345;
diff --git a/test/monniaux/int_mat_run.c b/test/monniaux/int_mat_run.c
index 3091f929..b8d836af 100644
--- a/test/monniaux/int_mat_run.c
+++ b/test/monniaux/int_mat_run.c
@@ -32,7 +32,7 @@ static inline cycle_t get_cycle(void) { return 0; }
 #endif
 
 int main() {
-  const unsigned m = 40, n = 20, p = 30;
+  const unsigned m = 40, n = 21, p = 30;
   cycle_count_config();
   modint *a = malloc(sizeof(modint) * m * n);
   modint_mat_random(m, n, a, n);
@@ -49,12 +49,21 @@ int main() {
   modint_mat_mul2(m, n, p, c2, p, a, n, b, p);
   c2_time = get_cycle()-c2_time;
   
-  printf("equal = %s\n"
+  modint *c3 = malloc(sizeof(modint) * m * p);
+  cycle_t c3_time = get_cycle();
+  modint_mat_mul3(m, n, p, c3, p, a, n, b, p);
+  c3_time = get_cycle()-c3_time;
+  
+  printf("c1==c2: %s\n"
+	 "c1==c3: %s\n"
 	 "c1_time = %" PRIu64 "\n"
-	 "c2_time = %" PRIu64 "\n",
+	 "c2_time = %" PRIu64 "\n"
+	 "c3_time = %" PRIu64 "\n",
 	 modint_mat_equal(m, n, c1, p, c2, p)?"true":"false",
+	 modint_mat_equal(m, n, c1, p, c3, p)?"true":"false",
 	 c1_time,
-	 c2_time);
+	 c2_time,
+	 c3_time);
   
   free(a);
   free(b);
diff --git a/test/monniaux/modint.h b/test/monniaux/modint.h
index f1c591f8..3636657d 100644
--- a/test/monniaux/modint.h
+++ b/test/monniaux/modint.h
@@ -14,6 +14,11 @@ void modint_mat_mul2(unsigned m, unsigned n, unsigned p,
 		     const modint *a, unsigned stride_a,
 		     const modint *b, unsigned stride_b);
 
+void modint_mat_mul3(unsigned m, unsigned n, unsigned p,
+		     modint * restrict c, unsigned stride_c,
+		     const modint *a, unsigned stride_a,
+		     const modint *b, unsigned stride_b);
+
 modint modint_random(void);
 
 void modint_mat_random(unsigned m,