new FFM_mul (and FlxqM_mul, FqM_mul, ...)

[Peter Bruin <P.Bruin@warwick.ac.uk>]

Hello,

Here are a couple of functions for multiplying matrices over finite
fields.  I would be very happy if these functions (with any appropriate
improvements) could be included in the PARI library.

I implemented F2xqM_mul, FlxqM_mul, FqM_mul, the variants where the
right-hand side is a column vector (F2xqM_F2xqC_mul, etc.), and FFM_mul.
All of the functions are pretty straightforward.  I personally didn't
need FFM_FFC_mul or anything involving FFV, so I didn't implement those.
The patch also modifies RgM_mul to call FFM_mul when appropriate.  I ran
a small number of examples; the new functions seem to be up to 2.5 times
as fast as the existing code.

One thing I'm not sure about is how much consistency checking to do in
the lower-level functions (FlxqM_FlxqC_mul, etc.); currently they check
if the dimensions are compatible but not if the input types (t_MAT etc.)
are correct.

I also ran gcov and discovered that the new functions are covered by the
existing test suite, so I didn't add any new tests.

Thanks,

Peter


diff --git a/doc/usersch5.tex b/doc/usersch5.tex
index 42680b3..55159fc 100644
--- a/doc/usersch5.tex
+++ b/doc/usersch5.tex
@@ -3757,6 +3757,8 @@ as \kbd{gauss}, where $b$ is a \kbd{FqM}.
 \fun{GEN}{FqM_FqC_gauss}{GEN a, GEN b, GEN T, GEN p}
 as \kbd{gauss}, where $b$ is a \kbd{FqC}.
 
+\fun{GEN}{FqM_FqC_mul}{GEN a, GEN b, GEN T, GEN p}
+
 \fun{GEN}{FqM_ker}{GEN x, GEN T, GEN p} as \kbd{ker}
 
 \fun{GEN}{FqM_image}{GEN x, GEN T, GEN p} as \kbd{image}
@@ -3764,6 +3766,8 @@ as \kbd{gauss}, where $b$ is a \kbd{FqC}.
 \fun{GEN}{FqM_inv}{GEN x, GEN T, GEN p} returns the inverse of \kbd{x}, or
 \kbd{NULL} if \kbd{x} is not invertible.
 
+\fun{GEN}{FqM_FqC_mul}{GEN a, GEN b, GEN T, GEN p}
+
 \fun{long}{FqM_rank}{GEN x, GEN T, GEN p} as \kbd{rank}
 
 \fun{GEN}{FqM_suppl}{GEN x, GEN T, GEN p} as \kbd{suppl}
@@ -4033,6 +4037,8 @@ and \kbd{y}
 
 \fun{GEN}{FlxqM_FlxqC_gauss}{GEN a, GEN b, GEN T, ulong p}
 
+\fun{GEN}{FlxqM_FlxqC_mul}{GEN a, GEN b, GEN T, ulong p}
+
 \fun{GEN}{FlxqM_ker}{GEN x, GEN T, ulong p}
 
 \fun{GEN}{FlxqM_image}{GEN x, GEN T, ulong p}
@@ -4041,6 +4047,8 @@ and \kbd{y}
 
 \fun{GEN}{FlxqM_inv}{GEN x, GEN T, ulong p}
 
+\fun{GEN}{FlxqM_mul}{GEN a, GEN b, GEN T, ulong p}
+
 \fun{long}{FlxqM_rank}{GEN x, GEN T, ulong p}
 
 \fun{GEN}{matid_FlxqM}{long n, GEN T, ulong p}
@@ -5395,6 +5403,8 @@ $a=\sigma(X)$ where $\sigma$ is an automorphism of the algebra $\F_2[X]/T(X)$.
 
 \subsec{\kbd{F2xqV}, \kbd{F2xqM}}. See \kbd{FqV}, \kbd{FqM} operations.
 
+\fun{GEN}{F2xqM_F2xqC_mul}{GEN a, GEN b, GEN T}
+
 \fun{GEN}{F2xqM_ker}{GEN x, GEN T}
 
 \fun{GEN}{F2xqM_det}{GEN a, GEN T}
@@ -5403,6 +5413,8 @@ $a=\sigma(X)$ where $\sigma$ is an automorphism of the algebra $\F_2[X]/T(X)$.
 
 \fun{GEN}{F2xqM_inv}{GEN a, GEN T}
 
+\fun{GEN}{F2xqM_mul}{GEN a, GEN b, GEN T}
+
 \fun{long}{F2xqM_rank}{GEN x, GEN T}
 
 \fun{GEN}{matid_F2xqM}{long n, GEN T}
@@ -9090,9 +9102,9 @@ of the univariate polynomial \kbd{f} over the definition field of the
 \typ{FFELT} \kbd{a}. The coefficients of \kbd{f} must be of type \typ{INT},
 \typ{INTMOD} or \typ{FFELT} and compatible with \kbd{a}.
 
-\fun{GEN}{FFM_ker}{GEN M, GEN ff} return the kernel of the \typ{MAT} \kbd{M}
+\fun{GEN}{FFM_ker}{GEN M, GEN ff} returns the kernel of the \typ{MAT} \kbd{M}
 defined over the finite field given by the \typ{FFELT} \kbd{ff} (obtained
-by \tet{RgM_is_FFM(M,\&ff)}.
+by \tet{RgM_is_FFM(M,\&ff)}).
 
 \fun{GEN}{FFM_det}{GEN M, GEN ff}
 
@@ -9100,6 +9112,10 @@ by \tet{RgM_is_FFM(M,\&ff)}.
 
 \fun{GEN}{FFM_inv}{GEN M, GEN ff}
 
+\fun{GEN}{FFM_mul}{GEN M, GEN N, GEN ff} returns the product of the
+matrices \kbd{M} and~\kbd{N} (\typ{MAT}) over the finite field given
+by \kbd{ff} (\typ{FFELT}).
+
 \fun{long}{FFM_rank}{GEN M, GEN ff}
 
 \section{Transcendental functions}
diff --git a/src/basemath/FF.c b/src/basemath/FF.c
index 00aee57..2ec1395 100644
--- a/src/basemath/FF.c
+++ b/src/basemath/FF.c
@@ -1738,3 +1738,21 @@ FFM_det(GEN M, GEN ff)
   }
   return gerepilecopy(av, mkFF_i(ff, d));
 }
+
+GEN
+FFM_mul(GEN M, GEN N, GEN ff)
+{
+  pari_sp av = avma;
+  ulong pp;
+  GEN P, T, p;
+  _getFF(ff,&T,&p,&pp);
+  M = FFM_to_raw(M);
+  N = FFM_to_raw(N);
+  switch(ff[1])
+  {
+  case t_FF_FpXQ: P = FqM_mul(M,N,T,p); break;
+  case t_FF_F2xq: P = F2xqM_mul(M,N,T); break;
+  default: P = FlxqM_mul(M,N,T,pp); break;
+  }
+  return gerepilecopy(av, raw_to_FFM(P, ff));
+}
diff --git a/src/basemath/RgV.c b/src/basemath/RgV.c
index c5dac90..4fab6d3 100644
--- a/src/basemath/RgV.c
+++ b/src/basemath/RgV.c
@@ -464,12 +464,17 @@ RgM_mul(GEN x, GEN y)
 {
   pari_sp av = avma;
   long j, l, lx, ly = lg(y);
-  GEN z;
+  GEN z, ffx = NULL, ffy = NULL;
   if (ly == 1) return cgetg(1,t_MAT);
   lx = lg(x);
   if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
   if (lx == 1) return zeromat(0,lx-1);
   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))
+      pari_err_OP("*", ffx, ffy);
+    return FFM_mul(x, y, ffx);
+  }
   z = cgetg(ly, t_MAT);
   l = lgcols(x);
   for (j=1; j<ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
diff --git a/src/basemath/alglin1.c b/src/basemath/alglin1.c
index 43ccc91..ca4046f 100644
--- a/src/basemath/alglin1.c
+++ b/src/basemath/alglin1.c
@@ -917,6 +917,37 @@ FlxqM_det(GEN a, GEN T, ulong p)
   const struct bb_field *S = get_Flxq_field(&E, T, p);
   return gen_det(a, E, S);
 }
+
+GEN
+FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, unsigned long p) {
+  long i, l, lgA = lg(A);
+  GEN C;
+  if(lgA != lg(B))
+    pari_err_OP("* [FlxqM_FlxqC]", A, B);
+  if(lgA == 1)
+    return cgetg(1, t_COL);
+  l = lgcols(A);
+  C = cgetg(l, t_COL);
+  for(i = 1; i < l; i++) {
+    pari_sp av = avma;
+    GEN e = Flx_mul(gcoeff(A, i, 1), gel(B, 1), p);
+    long k;
+    for(k = 2; k < lgA; k++)
+      e = Flx_add(e, Flx_mul(gcoeff(A, i, k), gel(B, k), p), p);
+    gel(C, i) = gerepileuptoleaf(av, Flx_rem(e, T, p));
+  }
+  return C;
+}
+
+GEN
+FlxqM_mul(GEN A, GEN B, GEN T, unsigned long p) {
+  long j, lgB = lg(B);
+  GEN C = cgetg(lgB, t_MAT);
+  for(j = 1; j < lgB; j++)
+    gel(C, j) = FlxqM_FlxqC_mul(A, gel(B, j), T, p);
+  return C;
+}
+
 GEN
 F2xqM_det(GEN a, GEN T)
 {
@@ -944,6 +975,36 @@ F2xqM_inv(GEN a, GEN T)
   return gerepilecopy(av, u);
 }
 
+GEN
+F2xqM_F2xqC_mul(GEN A, GEN B, GEN T) {
+  long i, l, lgA = lg(A);
+  GEN C;
+  if(lgA != lg(B))
+    pari_err_OP("* [F2xqM_F2xqC]", A, B);
+  if(lgA == 1)
+    return cgetg(1, t_COL);
+  l = lgcols(A);
+  C = cgetg(l, t_COL);
+  for(i = 1; i < l; i++) {
+    pari_sp av = avma;
+    GEN e = F2x_mul(gcoeff(A, i, 1), gel(B, 1));
+    long k;
+    for(k = 2; k < lgA; k++)
+      e = F2x_add(e, F2x_mul(gcoeff(A, i, k), gel(B, k)));
+    gel(C, i) = gerepileuptoleaf(av, F2x_rem(e, T));
+  }
+  return C;
+}
+
+GEN
+F2xqM_mul(GEN A, GEN B, GEN T) {
+  long j, lgB = lg(B);
+  GEN C = cgetg(lgB, t_MAT);
+  for(j = 1; j < lgB; j++)
+    gel(C, j) = F2xqM_F2xqC_mul(A, gel(B, j), T);
+  return C;
+}
+
 static GEN
 FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
 {
@@ -992,6 +1053,36 @@ FqM_det(GEN x, GEN T, GEN p)
   return gen_det(x, E, S);
 }
 
+GEN
+FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
+  long i, l, lgA = lg(A);
+  GEN C;
+  if(lgA != lg(B))
+    pari_err_OP("* [FqM_FqC]", A, B);
+  if(lgA == 1)
+    return cgetg(1, t_COL);
+  l = lgcols(A);
+  C = cgetg(l, t_COL);
+  for(i = 1; i < l; i++) {
+    pari_sp av = avma;
+    GEN e = FpX_mul(gcoeff(A, i, 1), gel(B, 1), p);
+    long k;
+    for(k = 2; k < lgA; k++)
+      e = FpX_add(e, FpX_mul(gcoeff(A, i, k), gel(B, k), p), p);
+    gel(C, i) = gerepileupto(av, FpX_rem(e, T, p));
+  }
+  return C;
+}
+
+GEN
+FqM_mul(GEN A, GEN B, GEN T, GEN p) {
+  long j, lgB = lg(B);
+  GEN C = cgetg(lgB, t_MAT);
+  for(j = 1; j < lgB; j++)
+    gel(C, j) = FqM_FqC_mul(A, gel(B, j), T, p);
+  return C;
+}
+
 static GEN
 FpM_ker_gen(GEN x, GEN p, long deplin)
 {
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
index 5ffb570..43ce0c0 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -954,10 +954,12 @@ GEN     F2m_ker(GEN x);
 GEN     F2m_ker_sp(GEN x, long deplin);
 long    F2m_rank(GEN x);
 GEN     F2m_suppl(GEN x);
+GEN     F2xqM_F2xqC_mul(GEN a, GEN b, GEN T);
 GEN     F2xqM_det(GEN a, GEN T);
 GEN     F2xqM_ker(GEN x, GEN T);
 GEN     F2xqM_image(GEN x, GEN T);
 GEN     F2xqM_inv(GEN a, GEN T);
+GEN     F2xqM_mul(GEN a, GEN b, GEN T);
 long    F2xqM_rank(GEN x, GEN T);
 GEN     Flm_Flc_gauss(GEN a, GEN b, ulong p);
 GEN     Flm_Flc_invimage(GEN mat, GEN y, ulong p);
@@ -974,11 +976,13 @@ GEN     Flm_ker_sp(GEN x, ulong p, long deplin);
 long    Flm_rank(GEN x, ulong p);
 GEN     Flm_suppl(GEN x, ulong p);
 GEN     FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p);
+GEN     FlxqM_FlxqC_mul(GEN a, GEN b, GEN T, ulong p);
 GEN     FlxqM_det(GEN a, GEN T, ulong p);
 GEN     FlxqM_gauss(GEN a, GEN b, GEN T, ulong p);
 GEN     FlxqM_ker(GEN x, GEN T, ulong p);
 GEN     FlxqM_image(GEN x, GEN T, ulong p);
 GEN     FlxqM_inv(GEN x, GEN T, ulong p);
+GEN     FlxqM_mul(GEN a, GEN b, GEN T, ulong p);
 long    FlxqM_rank(GEN x, GEN T, ulong p);
 GEN     FpM_FpC_gauss(GEN a, GEN b, GEN p);
 GEN     FpM_FpC_invimage(GEN m, GEN v, GEN p);
@@ -994,12 +998,14 @@ GEN     FpM_ker(GEN x, GEN p);
 long    FpM_rank(GEN x, GEN p);
 GEN     FpM_suppl(GEN x, GEN p);
 GEN     FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p);
+GEN     FqM_FqC_mul(GEN a, GEN b, GEN T, GEN p);
 GEN     FqM_deplin(GEN x, GEN T, GEN p);
 GEN     FqM_det(GEN x, GEN T, GEN p);
 GEN     FqM_gauss(GEN a, GEN b, GEN T, GEN p);
 GEN     FqM_ker(GEN x, GEN T, GEN p);
 GEN     FqM_image(GEN x, GEN T, GEN p);
 GEN     FqM_inv(GEN x, GEN T, GEN p);
+GEN     FqM_mul(GEN a, GEN b, GEN T, GEN p);
 long    FqM_rank(GEN a, GEN T, GEN p);
 GEN     FqM_suppl(GEN x, GEN T, GEN p);
 GEN     QM_inv(GEN M, GEN dM);
@@ -2194,6 +2200,7 @@ GEN     FFM_det(GEN M, GEN ff);
 GEN     FFM_image(GEN M, GEN ff);
 GEN     FFM_inv(GEN M, GEN ff);
 GEN     FFM_ker(GEN M, GEN ff);
+GEN     FFM_mul(GEN M, GEN N, GEN ff);
 long    FFM_rank(GEN M, GEN ff);
 GEN     FFX_factor(GEN f, GEN x);
 GEN     FFX_roots(GEN f, GEN x);