Re: Strassen multiplication over the integers

[Peter Bruin <P.J.Bruin@math.leidenuniv.nl>]

Hello Bill,

> I have applied your patch with very minor modifications.

Merci !

> However there are things that need to be done.
>
> - ZM_sqr should be implemented as well.
>
> - gmul or RgM_mul should be changed to call ZM_mul when appropriate.

The above two points are addressed by the attached patch.  (Note that
the assumption that A == B in A*B doesn't really seem to simplify either
the classical algorithm or Strassen's algorithm.  Even the following
could be reasonable: "GEN ZM_sqr(GEN x) { return ZM_mul(x, x); }" ...)

> - the tuning should depend on the size of the coefficients, I think.

Certainly, for larger coefficients it is advantageous to use Strassen
multiplication already for smaller matrices.  I haven't yet done any
good experiments to see how the bound should depend on the maximum
coefficient size, though.

> - ideally the program src/test/tune.c should deal with the tuning.
> This program deal with dependencies between tuning parameters.

OK, I hope to have some time soon to try to write an addition for
src/test/tune.c for this.

> - Most of the code is fairly generic, so maybe there could be a
> gen_matmul_sw function.

Indeed, and I actually already have an implementation of exactly this
function, but it is an older version with awkward conventions for the
indices.  I will try to update this soon.  (Unfortunately it will be
hard to tune that one...)

Thanks,

Peter


diff --git a/src/basemath/RgV.c b/src/basemath/RgV.c
index 8bfaeea..98ef4c8 100644
--- a/src/basemath/RgV.c
+++ b/src/basemath/RgV.c
@@ -526,6 +526,8 @@ RgM_mul(GEN x, GEN y)
   lx = lg(x);
   if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
   if (lx == 1) return zeromat(0,ly-1);
+  if (RgM_is_ZM(x) && RgM_is_ZM(y))
+    return ZM_mul(x, y);
   if (is_modular_mul(x,y,&z)) return gerepileupto(av, z);
   if (RgM_is_FFM(x, &ffx) && RgM_is_FFM(y, &ffy)) {
     if (!FF_samefield(ffx, ffy))
@@ -628,6 +630,7 @@ RgM_sqr(GEN x)
   GEN z, ffx = NULL;
   if (lx == 1) return cgetg(1, t_MAT);
   if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
+  if (RgM_is_ZM(x))         return ZM_sqr(x);
   if (is_modular_sqr(x,&z)) return gerepileupto(av, z);
   if (RgM_is_FFM(x, &ffx))  return FFM_mul(x, x, ffx);
   z = cgetg(lx, t_MAT);
diff --git a/src/basemath/ZV.c b/src/basemath/ZV.c
index 332b243..0b3cc3e 100644
--- a/src/basemath/ZV.c
+++ b/src/basemath/ZV.c
@@ -465,15 +465,22 @@ ZM_transmul(GEN x, GEN y)
   }
   return M;
 }
+
+static GEN
+ZM_sqr_i(GEN x, long l)
+{
+  if (l <= ZM_sw_bound)
+    return ZM_mul_classical(x, x, l, l, l);
+  else
+    return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
+}
+
 GEN
 ZM_sqr(GEN x)
 {
-  long j, l, lx=lg(x);
-  GEN z;
+  long lx=lg(x);
   if (lx==1) return cgetg(1,t_MAT);
-  l = lgcols(x); z = cgetg(lx,t_MAT);
-  for (j=1; j<lx; j++) gel(z,j) = ZM_ZC_mul_i(x, gel(x,j), lx, l);
-  return z;
+  return ZM_sqr_i(x, lx);
 }
 GEN
 ZM_ZC_mul(GEN x, GEN y)
@@ -564,7 +571,7 @@ _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
 { (void)data; return ZM_mul(x,y); }
 static GEN
 _ZM_sqr(void *data /*ignored*/, GEN x)
-{ (void)data; return ZM_mul(x,x); }
+{ (void)data; return ZM_sqr(x); }
 GEN
 ZM_pow(GEN x, GEN n)
 {