patch for finite field versions of matinverseimage

Peter Bruin on Fri, 28 Mar 2014 16:08:15 +0100

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

patch for finite field versions of matinverseimage

To: <pari-dev@pari.math.u-bordeaux.fr>
Subject: patch for finite field versions of matinverseimage
From: Peter Bruin <P.Bruin@warwick.ac.uk>
Date: Fri, 28 Mar 2014 15:08:00 +0000
Delivery-date: Fri, 28 Mar 2014 16:08:15 +0100
User-agent: Gnus/5.13 (Gnus v5.13) Emacs/24.3 (gnu/linux)

Hello,

Some time ago I mentioned having written some functions to speed up
matinverseimage for matrices/vectors with t_FFELT entries.  I polished
the code a bit and added some tests that should cover all new functions;
see the attached patch.

I realise that the functions gen_matinvimage and gen_matcolinvimage
(where the real work is done) is still written in a way that would
ideally be replaced by more standard primitives, as in Karim's message
quoted below.  Since I don't have that much time to spare at the moment,
I thought I would just send the patch as it is, hoping it is still
useful enough.

Thanks,

Peter

diff --git a/doc/develop.tex b/doc/develop.tex
index a711cec..0cafc15 100644
--- a/doc/develop.tex
+++ b/doc/develop.tex
@@ -432,10 +432,14 @@ fields:
 
 \fun{GEN}{gen_ker}{GEN x, long deplin, void *E, const struct bb_field *ff}
 
+\fun{GEN}{gen_matcolinvimage}{GEN a, GEN b, void *E, const struct bb_field *ff}
+
 \fun{GEN}{gen_matcolmul}{GEN a, GEN b, void *E, const struct bb_field *ff}
 
 \fun{GEN}{gen_matid}{long n, void *E, const struct bb_field *ff}
 
+\fun{GEN}{gen_matinvimage}{GEN a, GEN b, void *E, const struct bb_field *ff}
+
 \fun{GEN}{gen_matmul}{GEN a, GEN b, void *E, const struct bb_field *ff}
 
 \subsec{Functions returning black box fields}
diff --git a/doc/usersch5.tex b/doc/usersch5.tex
index d115321..81fca52 100644
--- a/doc/usersch5.tex
+++ b/doc/usersch5.tex
@@ -3454,6 +3454,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_invimage}{GEN a, GEN b, GEN T, GEN p}
+
 \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}
@@ -3463,6 +3465,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_invimage}{GEN a, GEN b, GEN T, GEN p}
+
 \fun{GEN}{FqM_mul}{GEN a, GEN b, GEN T, GEN p}
 
 \fun{long}{FqM_rank}{GEN x, GEN T, GEN p} as \kbd{rank}
@@ -3740,6 +3744,8 @@ and \kbd{y}
 
 \fun{GEN}{FlxqM_FlxqC_gauss}{GEN a, GEN b, GEN T, ulong p}
 
+\fun{GEN}{FlxqM_FlxqC_invimage}{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}
@@ -3750,6 +3756,8 @@ and \kbd{y}
 
 \fun{GEN}{FlxqM_inv}{GEN x, GEN T, ulong p}
 
+\fun{GEN}{FlxqM_invimage}{GEN a, GEN b, 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}
@@ -5120,6 +5128,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_invimage}{GEN a, GEN b, GEN T}
+
 \fun{GEN}{F2xqM_F2xqC_mul}{GEN a, GEN b, GEN T}
 
 \fun{GEN}{F2xqM_ker}{GEN x, GEN T}
@@ -5130,6 +5140,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_invimage}{GEN a, GEN b, GEN T}
+
 \fun{GEN}{F2xqM_mul}{GEN a, GEN b, GEN T}
 
 \fun{long}{F2xqM_rank}{GEN x, GEN T}
@@ -9296,6 +9308,11 @@ 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_FFC_invimage}{GEN M, GEN C, GEN ff} given a matrix
+\kbd{M} (\typ{MAT}) and a column vector \kbd{C} (\typ{COL}) over the
+finite field given by \kbd{ff} (\typ{FFELT}), return a column vector
+\kbd{X} such that $MX=C$, or \kbd{NULL} if no such vector exists.
+
 \fun{GEN}{FFM_FFC_mul}{GEN M, GEN C, GEN ff} returns the product of
 the matrix~\kbd{M} (\typ{MAT}) and the column vector~\kbd{C}
 (\typ{COL}) over the finite field given by \kbd{ff} (\typ{FFELT}).
@@ -9310,6 +9327,11 @@ by \tet{RgM_is_FFM(M,\&ff)}).
 
 \fun{GEN}{FFM_inv}{GEN M, GEN ff}
 
+\fun{GEN}{FFM_invimage}{GEN M, GEN N, GEN ff} given two matrices
+\kbd{M} and~\kbd{N} (\typ{MAT}) over the finite field given by
+\kbd{ff} (\typ{FFELT}), return a matrix \kbd{X} such that $MX=N$, or
+\kbd{NULL} if no such matrix exists.
+
 \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}).
diff --git a/src/basemath/FF.c b/src/basemath/FF.c
index 1dfb9aa..1bf0a51 100644
--- a/src/basemath/FF.c
+++ b/src/basemath/FF.c
@@ -1681,6 +1681,7 @@ FqM_to_FpXQM(GEN x, GEN T)
   return y;
 }
 
+/* for functions t_MAT -> t_MAT */
 static GEN
 FFM_wrap(GEN M, GEN ff, GEN (*Fq)(GEN,GEN,GEN),
                        GEN (*Flxq)(GEN,GEN,ulong),
@@ -1698,6 +1699,53 @@ FFM_wrap(GEN M, GEN ff, GEN (*Fq)(GEN,GEN,GEN),
   }
   return gerepilecopy(av, raw_to_FFM(M, ff));
 }
+
+/* for functions (t_MAT, t_MAT) -> t_MAT */
+static GEN
+FFM_FFM_wrap(GEN M, GEN N, GEN ff,
+	     GEN (*Fq)(GEN, GEN, GEN, GEN),
+	     GEN (*Flxq)(GEN, GEN, GEN, ulong),
+	     GEN (*F2xq)(GEN, GEN, GEN))
+{
+  pari_sp av = avma;
+  ulong pp;
+  GEN T, p;
+  _getFF(ff, &T, &p, &pp);
+  M = FFM_to_raw(M);
+  N = FFM_to_raw(N);
+  switch(ff[1])
+  {
+  case t_FF_FpXQ: M = FqM_to_FpXQM(Fq(M, N, T, p), T); break;
+  case t_FF_F2xq: M = F2xq(M, N, T); break;
+  default: M = Flxq(M, N, T, pp); break;
+  }
+  if (M == NULL) return NULL;
+  return gerepilecopy(av, raw_to_FFM(M, ff));
+}
+
+/* for functions (t_MAT, t_COL) -> t_COL */
+static GEN
+FFM_FFC_wrap(GEN M, GEN C, GEN ff,
+	     GEN (*Fq)(GEN, GEN, GEN, GEN),
+	     GEN (*Flxq)(GEN, GEN, GEN, ulong),
+	     GEN (*F2xq)(GEN, GEN, GEN))
+{
+  pari_sp av = avma;
+  ulong pp;
+  GEN T, p;
+  _getFF(ff, &T, &p, &pp);
+  M = FFM_to_raw(M);
+  C = FFC_to_raw(C);
+  switch(ff[1])
+  {
+  case t_FF_FpXQ: C = FqC_to_FpXQC(Fq(M, C, T, p), T); break;
+  case t_FF_F2xq: C = F2xq(M, C, T); break;
+  default: C = Flxq(M, C, T, pp); break;
+  }
+  if (C == NULL) return NULL;
+  return gerepilecopy(av, raw_to_FFC(C, ff));
+}
+
 GEN
 FFM_ker(GEN M, GEN ff)
 { return FFM_wrap(M,ff, &FqM_ker,&FlxqM_ker,&F2xqM_ker); }
@@ -1740,37 +1788,26 @@ FFM_det(GEN M, GEN ff)
 }
 
 GEN
+FFM_FFC_invimage(GEN M, GEN C, GEN ff)
+{
+  return FFM_FFC_wrap(M, C, ff, &FqM_FqC_invimage,
+		      &FlxqM_FlxqC_invimage, &F2xqM_F2xqC_invimage);
+}
+
+GEN
+FFM_invimage(GEN M, GEN N, GEN ff)
+{
+  return FFM_FFM_wrap(M, N, ff, &FqM_invimage, &FlxqM_invimage, &F2xqM_invimage);
+}
+
+GEN
 FFM_FFC_mul(GEN M, GEN C, GEN ff)
 {
-  pari_sp av = avma;
-  ulong pp;
-  GEN P, T, p;
-  _getFF(ff, &T, &p, &pp);
-  M = FFM_to_raw(M);
-  C = FFC_to_raw(C);
-  switch (ff[1])
-  {
-  case t_FF_FpXQ: P = FqM_FqC_mul(M, C, T, p); break;
-  case t_FF_F2xq: P = F2xqM_F2xqC_mul(M, C, T); break;
-  default: P = FlxqM_FlxqC_mul(M, C, T, pp); break;
-  }
-  return gerepilecopy(av, raw_to_FFC(P, ff));
+  return FFM_FFC_wrap(M, C, ff, &FqM_FqC_mul, &FlxqM_FlxqC_mul, &F2xqM_F2xqC_mul);
 }
 
 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));
+  return FFM_FFM_wrap(M, N, ff, &FqM_mul, &FlxqM_mul, &F2xqM_mul);
 }
diff --git a/src/basemath/alglin1.c b/src/basemath/alglin1.c
index d3868dc..dad10a8 100644
--- a/src/basemath/alglin1.c
+++ b/src/basemath/alglin1.c
@@ -408,6 +408,107 @@ gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
 }
 
 static GEN
+gen_colneg(GEN A, void *E, const struct bb_field *ff)
+{
+  long i, l;
+  GEN B = cgetg_copy(A, &l);
+  for (i = 1; i < l; i++)
+    gel(B, i) = ff->neg(E, gel(A, i));
+  return B;
+}
+
+static GEN
+gen_matneg(GEN A, void *E, const struct bb_field *ff)
+{
+  long i, l;
+  GEN B = cgetg_copy(A, &l);
+  for (i = 1; i < l; i++)
+    gel(B, i) = gen_colneg(gel(A, i), E, ff);
+  return B;
+}
+
+/* assume dim A >= 1, A invertible + upper triangular  */
+static GEN
+gen_matinv_upper_ind(GEN A, long index, void *E, const struct bb_field *ff)
+{
+  long n = lg(A) - 1, i, j;
+  GEN u = cgetg(n + 1, t_COL);
+  for (i = n; i > index; i--)
+    gel(u, i) = ff->s(E, 0);
+  gel(u, i) = ff->inv(E, gcoeff(A, i, i));
+  for (i--; i > 0; i--) {
+    pari_sp av = avma;
+    GEN m = ff->neg(E, ff->mul(E, gcoeff(A, i, i + 1), gel(u, i + 1)));
+    for (j = i + 2; j <= n; j++)
+      m = ff->add(E, m, ff->neg(E, ff->mul(E, gcoeff(A, i, j), gel(u, j))));
+    gel(u, i) = gerepileupto(av, ff->red(E, ff->mul(E, m, ff->inv(E, gcoeff(A, i, i)))));
+  }
+  return u;
+}
+
+static GEN
+gen_matinv_upper(GEN A, void *E, const struct bb_field *ff)
+{
+  long i, l;
+  GEN B = cgetg_copy(A, &l);
+  for (i = 1; i < l; i++)
+    gel(B,i) = gen_matinv_upper_ind(A, i, E, ff);
+  return B;
+}
+
+/* find z such that A z = y. Return NULL if no solution */
+GEN
+gen_matcolinvimage(GEN A, GEN y, void *E, const struct bb_field *ff)
+{
+  pari_sp av = avma;
+  long i, l = lg(A);
+  GEN M, x, t;
+
+  M = gen_ker(shallowconcat(A, y), 0, E, ff);
+  i = lg(M) - 1;
+  if (!i) { avma = av; return NULL; }
+
+  x = gel(M, i);
+  t = gel(x, l);
+  if (ff->equal0(t)) { avma = av; return NULL; }
+
+  t = ff->neg(E, ff->inv(E, t));
+  setlg(x, l);
+  for (i = 1; i < l; i++)
+    gel(x, i) = ff->red(E, ff->mul(E, t, gel(x, i)));
+  return gerepilecopy(av, x);
+}
+
+/* find Z such that A Z = B. Return NULL if no solution */
+GEN
+gen_matinvimage(GEN A, GEN B, void *E, const struct bb_field *ff)
+{
+  pari_sp av = avma;
+  GEN d, x, X, Y;
+  long i, j, nY, nA, nB;
+  x = gen_ker(shallowconcat(gen_matneg(A, E, ff), B), 0, E, ff);
+  /* AX = BY, Y in strict upper echelon form with pivots = 1.
+   * We must find T such that Y T = Id_nB then X T = Z. This exists
+   * iff Y has at least nB columns and full rank. */
+  nY = lg(x) - 1;
+  nB = lg(B) - 1;
+  if (nY < nB) { avma = av; return NULL; }
+  nA = lg(A) - 1;
+  Y = rowslice(x, nA + 1, nA + nB); /* nB rows */
+  d = cgetg(nB + 1, t_VECSMALL);
+  for (i = nB, j = nY; i >= 1; i--) {
+    for (; j >= 1; j--)
+      if (!ff->equal0(gcoeff(Y, i, j))) { d[i] = j; break; }
+    if (!j) { avma = av; return NULL; }
+  }
+  /* reduce to the case Y square, upper triangular with 1s on diagonal */
+  Y = vecpermute(Y, d);
+  x = vecpermute(x, d);
+  X = rowslice(x, 1, nA);
+  return gerepileupto(av, gen_matmul(X, gen_matinv_upper(Y, E, ff), E, ff));
+}
+
+static GEN
 image_from_pivot(GEN x, GEN d, long r)
 {
   GEN y;
@@ -966,20 +1067,34 @@ FlxqM_det(GEN a, GEN T, ulong p)
 }
 
 GEN
-FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, unsigned long p) {
+FlxqM_FlxqC_invimage(GEN A, GEN B, GEN T, ulong p) {
+  void *E;
+  const struct bb_field *ff = get_Flxq_field(&E, T, p);
+  return gen_matcolinvimage(A, B, E, ff);
+}
+
+GEN
+FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, ulong p) {
   void *E;
   const struct bb_field *ff = get_Flxq_field(&E, T, p);
   return gen_matcolmul(A, B, E, ff);
 }
 
 GEN
-FlxqM_mul(GEN A, GEN B, GEN T, unsigned long p) {
+FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
   void *E;
   const struct bb_field *ff = get_Flxq_field(&E, T, p);
   return gen_matmul(A, B, E, ff);
 }
 
 GEN
+FlxqM_invimage(GEN A, GEN B, GEN T, ulong p) {
+  void *E;
+  const struct bb_field *ff = get_Flxq_field(&E, T, p);
+  return gen_matinvimage(A, B, E, ff);
+}
+
+GEN
 F2xqM_det(GEN a, GEN T)
 {
   void *E;
@@ -1007,6 +1122,13 @@ F2xqM_inv(GEN a, GEN T)
 }
 
 GEN
+F2xqM_F2xqC_invimage(GEN A, GEN B, GEN T) {
+  void *E;
+  const struct bb_field *ff = get_F2xq_field(&E, T);
+  return gen_matcolinvimage(A, B, E, ff);
+}
+
+GEN
 F2xqM_F2xqC_mul(GEN A, GEN B, GEN T) {
   void *E;
   const struct bb_field *ff = get_F2xq_field(&E, T);
@@ -1020,6 +1142,13 @@ F2xqM_mul(GEN A, GEN B, GEN T) {
   return gen_matmul(A, B, E, ff);
 }
 
+GEN
+F2xqM_invimage(GEN A, GEN B, GEN T) {
+  void *E;
+  const struct bb_field *ff = get_F2xq_field(&E, T);
+  return gen_matinvimage(A, B, E, ff);
+}
+
 static GEN
 FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
 {
@@ -1069,6 +1198,13 @@ FqM_det(GEN x, GEN T, GEN p)
 }
 
 GEN
+FqM_FqC_invimage(GEN A, GEN B, GEN T, GEN p) {
+  void *E;
+  const struct bb_field *ff = get_Fq_field(&E, T, p);
+  return gen_matcolinvimage(A, B, E, ff);
+}
+
+GEN
 FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
   void *E;
   const struct bb_field *ff = get_Fq_field(&E, T, p);
@@ -1082,6 +1218,13 @@ FqM_mul(GEN A, GEN B, GEN T, GEN p) {
   return gen_matmul(A, B, E, ff);
 }
 
+GEN
+FqM_invimage(GEN A, GEN B, GEN T, GEN p) {
+  void *E;
+  const struct bb_field *ff = get_Fq_field(&E, T, p);
+  return gen_matinvimage(A, B, E, ff);
+}
+
 static GEN
 FpM_ker_gen(GEN x, GEN p, long deplin)
 {
@@ -2867,6 +3010,8 @@ RgM_RgC_invimage(GEN A, GEN y)
   long i, l = lg(A);
   GEN M, x, t, p = NULL;
 
+  if (l==1) return NULL;
+  if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
   if (RgM_is_FpM(A, &p) && RgV_is_FpV(y, &p) && p)
   {
     ulong pp;
@@ -2891,9 +3036,11 @@ RgM_RgC_invimage(GEN A, GEN y)
     if (!x) { avma = av; return NULL; }
     return gerepileupto(av, x);
   }
-
-  if (l==1) return NULL;
-  if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
+  {
+    GEN ff = NULL;
+    if (RgM_is_FFM(A, &ff) && RgC_is_FFC(y, &ff))
+      return FFM_FFC_invimage(A, y, ff);
+  }
   M = ker(shallowconcat(A, y));
   i = lg(M)-1;
   if (!i) { avma = av; return NULL; }
@@ -3137,6 +3284,11 @@ RgM_invimage(GEN A, GEN B)
     if (!x) { avma = av; return NULL; }
     return gerepileupto(av, x);
   }
+  {
+    GEN ff = NULL;
+    if (RgM_is_FFM(A, &ff) && RgM_is_FFM(B, &ff))
+      return FFM_invimage(A, B, ff);
+  }
   x = ker(shallowconcat(RgM_neg(A), B));
   /* AX = BY, Y in strict upper echelon form with pivots = 1.
    * We must find T such that Y T = Id_nB then X T = Z. This exists iff
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
index f360963..2cbf0b7 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -969,11 +969,13 @@ 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_invimage(GEN a, GEN b, GEN T);
 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_invimage(GEN a, GEN b, 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);
@@ -991,12 +993,14 @@ 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_invimage(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_invimage(GEN a, GEN b, 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);
@@ -1013,6 +1017,7 @@ 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_invimage(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);
@@ -1020,6 +1025,7 @@ 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_invimage(GEN a, GEN b, 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);
@@ -1056,7 +1062,9 @@ GEN     gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff);
 GEN     gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff);
 GEN     gen_det(GEN a, void *E, const struct bb_field *ff);
 GEN     gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff);
+GEN     gen_matcolinvimage(GEN a, GEN b, void *E, const struct bb_field *ff);
 GEN     gen_matcolmul(GEN a, GEN b, void *E, const struct bb_field *ff);
+GEN     gen_matinvimage(GEN a, GEN b, void *E, const struct bb_field *ff);
 GEN     gen_matmul(GEN a, GEN b, void *E, const struct bb_field *ff);
 GEN     image(GEN x);
 GEN     image2(GEN x);
@@ -2227,10 +2235,12 @@ GEN     FF_to_FpXQ(GEN x);
 GEN     FF_to_FpXQ_i(GEN x);
 GEN     FF_trace(GEN x);
 GEN     FF_zero(GEN a);
+GEN     FFM_FFC_invimage(GEN M, GEN C, GEN ff);
 GEN     FFM_FFC_mul(GEN M, GEN C, GEN ff);
 GEN     FFM_det(GEN M, GEN ff);
 GEN     FFM_image(GEN M, GEN ff);
 GEN     FFM_inv(GEN M, GEN ff);
+GEN     FFM_invimage(GEN M, GEN N, 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);
diff --git a/src/test/32/ff b/src/test/32/ff
index d6333d6..6755757 100644
--- a/src/test/32/ff
+++ b/src/test/32/ff
@@ -240,21 +240,26 @@ t^2 + Mod(1, 5)*t + Mod(3, 5))*x^3 + Mod(Mod(0, 5), Mod(1, 5)*t^4 + Mod(1, 5
 od(4, 5)*t, Mod(1, 5)*t^4 + Mod(1, 5)*t^3 + Mod(2, 5)*t^2 + Mod(1, 5)*t + Mo
 d(3, 5)) 1]
 
-? test(q)=my(t=ffgen(q,'t),m=[t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matker(m));print(matimage(m));print(matrank(m));my(M=[t,2*t^0,3*t^0;t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matdet(M));print(M^(-1)*M);my(v=[t^0,t^1,t^2]~);print(M*v);
+? test(q)=my(t=ffgen(q,'t),m=[t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matker(m));print(matimage(m));print(matrank(m));my(M=[t,2*t^0,3*t^0;t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matdet(M));print(M^(-1)*M);print(matinverseimage(M,matid(3)*t^0)==M^(-1));my(v=[t^0,t^1,t^2]~);print(M*v);print(matinverseimage(M,v));
 ? test(2^5)
 [t^4 + t^3; t^4 + t^3; 1]
 [t, t^2; t + 1, t^2 + 1]
 2
 t^4 + t^2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [t^2 + t, t^4 + t^3 + 1, t^4 + t^3 + t]~
+[t^4, t^4 + t^2 + t, t^4 + t^2 + t]~
 ? test(7^5)
 [3*t^4 + 5*t^3 + 6*t^2 + 2*t; 4*t^4 + 2*t^3 + t^2 + 5*t; 1]
 [t, t^2; t + 1, t^2 + 1]
 2
 6*t^4 + 2*t^3 + 4*t^2 + 2*t + 2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [3*t^2 + 3*t, 6*t^4 + 5*t^3 + 4*t^2 + 5*t + 6, 6*t^4 + 5*t^3 + 4*t^2 + 6*t]~
+[6*t^4 + 3*t^2 + 4*t + 1, t^4 + 5*t^2 + 2*t + 6, 6*t^4 + 5*t^3 + t^2 + 2*t +
+ 3]~
 ? test((2^64+13)^5)
 [3*t^4 + 5*t^3 + 18446744073709551621*t^2 + 18446744073709551617*t; 18446744
 073709551626*t^4 + 18446744073709551624*t^3 + 8*t^2 + 12*t; 1]
@@ -262,9 +267,15 @@ t^4 + t^2
 2
 18446744073709551628*t^4 + 2*t^3 + 18446744073709551626*t^2 + 2*t + 2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [3*t^2 + 3*t, 18446744073709551628*t^4 + 5*t^3 + 4*t^2 + 1844674407370955162
 7*t + 18446744073709551628, 18446744073709551628*t^4 + 5*t^3 + 4*t^2 + 18446
 744073709551628*t]~
+[3488353224204417687*t^4 + 4892000116800957185*t^3 + 9053107177101941144*t^2
+ + 13820523250181311982*t + 390363336994303883, 14958390849505133942*t^4 + 1
+3554743956908594444*t^3 + 9393636896607610486*t^2 + 4626220823528239646*t + 
+18056380736715247746, 8006601209840615835*t^4 + 7890322769033801912*t^3 + 16
+386954550845990710*t^2 + 274084896187489962*t + 13720856015204042904]~
 ? p=2^64+13;g=ffprimroot(ffgen(p^2),&o);a=2*g^0;
 ? v=[I,-1,Mat(1),matid(2)/2];
 ? for(i=1,#v,print(iferr(fflog(a,g,v[i]),E,E)));
diff --git a/src/test/64/ff b/src/test/64/ff
index 2b163ca..69ea16b 100644
--- a/src/test/64/ff
+++ b/src/test/64/ff
@@ -243,21 +243,26 @@ t^2 + Mod(1, 5)*t + Mod(3, 5))*x^3 + Mod(Mod(0, 5), Mod(1, 5)*t^4 + Mod(1, 5
 od(4, 5)*t, Mod(1, 5)*t^4 + Mod(1, 5)*t^3 + Mod(2, 5)*t^2 + Mod(1, 5)*t + Mo
 d(3, 5)) 1]
 
-? test(q)=my(t=ffgen(q,'t),m=[t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matker(m));print(matimage(m));print(matrank(m));my(M=[t,2*t^0,3*t^0;t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matdet(M));print(M^(-1)*M);my(v=[t^0,t^1,t^2]~);print(M*v);
+? test(q)=my(t=ffgen(q,'t),m=[t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matker(m));print(matimage(m));print(matrank(m));my(M=[t,2*t^0,3*t^0;t,t^2,1+t^3;1+t,1+t^2,1+t^3]);print(matdet(M));print(M^(-1)*M);print(matinverseimage(M,matid(3)*t^0)==M^(-1));my(v=[t^0,t^1,t^2]~);print(M*v);print(matinverseimage(M,v));
 ? test(2^5)
 [t^4 + t^3; t^4 + t^3; 1]
 [t, t^2; t + 1, t^2 + 1]
 2
 t^4 + t^2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [t^2 + t, t^4 + t^3 + 1, t^4 + t^3 + t]~
+[t^4, t^4 + t^2 + t, t^4 + t^2 + t]~
 ? test(7^5)
 [3*t^4 + 5*t^3 + 6*t^2 + 2*t; 4*t^4 + 2*t^3 + t^2 + 5*t; 1]
 [t, t^2; t + 1, t^2 + 1]
 2
 6*t^4 + 2*t^3 + 4*t^2 + 2*t + 2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [3*t^2 + 3*t, 6*t^4 + 5*t^3 + 4*t^2 + 5*t + 6, 6*t^4 + 5*t^3 + 4*t^2 + 6*t]~
+[6*t^4 + 3*t^2 + 4*t + 1, t^4 + 5*t^2 + 2*t + 6, 6*t^4 + 5*t^3 + t^2 + 2*t +
+ 3]~
 ? test((2^64+13)^5)
 [3*t^4 + 5*t^3 + 18446744073709551621*t^2 + 18446744073709551617*t; 18446744
 073709551626*t^4 + 18446744073709551624*t^3 + 8*t^2 + 12*t; 1]
@@ -265,9 +270,15 @@ t^4 + t^2
 2
 18446744073709551628*t^4 + 2*t^3 + 18446744073709551626*t^2 + 2*t + 2
 [1, 0, 0; 0, 1, 0; 0, 0, 1]
+1
 [3*t^2 + 3*t, 18446744073709551628*t^4 + 5*t^3 + 4*t^2 + 1844674407370955162
 7*t + 18446744073709551628, 18446744073709551628*t^4 + 5*t^3 + 4*t^2 + 18446
 744073709551628*t]~
+[3488353224204417687*t^4 + 4892000116800957185*t^3 + 9053107177101941144*t^2
+ + 13820523250181311982*t + 390363336994303883, 14958390849505133942*t^4 + 1
+3554743956908594444*t^3 + 9393636896607610486*t^2 + 4626220823528239646*t + 
+18056380736715247746, 8006601209840615835*t^4 + 7890322769033801912*t^3 + 16
+386954550845990710*t^2 + 274084896187489962*t + 13720856015204042904]~
 ? p=2^64+13;g=ffprimroot(ffgen(p^2),&o);a=2*g^0;
 ? v=[I,-1,Mat(1),matid(2)/2];
 ? for(i=1,#v,print(iferr(fflog(a,g,v[i]),E,E)));
diff --git a/src/test/in/ff b/src/test/in/ff
index 6404029..2f6b98e 100644
--- a/src/test/in/ff
+++ b/src/test/in/ff
@@ -95,8 +95,10 @@ test(q)=
   my(M = [t,2*t^0,3*t^0; t,t^2,1+t^3; 1+t,1+t^2,1+t^3]);
   print(matdet(M));
   print(M^(-1)*M);
+  print(matinverseimage(M, matid(3)*t^0) == M^(-1));
   my(v = [t^0, t^1, t^2]~);
   print(M*v);
+  print(matinverseimage(M, v));
 }
 test(2^5)
 test(7^5)


Karim Belabas wrote:

> As a matter of fact, it's a little sad that we have so many "Gauss
> pivot" variants as primitives, instead of standard algorithms,
> e.g. LUP factorization. Currently
>
> - RgM_pivots (allow image / supplement / rank profile)
> - gauss_pivot_ker (analogous but more expensive, allow kernel)
> - gauss (compute A^(-1) B)
>
> Other operations are then implemented in terms of these non-standard
> primitives. (e.g. inverseimage reduces to kernel). These 3 are not too
> complicated, and share most non-trivial code (pivot search), but
> various minor optimizations related to memory use made them diverge
> (e.g. avoid a few operations, specialized gerepile).
>
> General linear algebra over fields has never been a strong focus of
> PARI, but cleaning up this set of ad hoc historical functions would be
> useful.  Again it's a bit sad that domain specialization (allowing
> large efficiency gains!) occured before that cleanup. As I see it, we
> could
>
> - merge RgM_pivots / gauss_pivot_ker
>
> - replace gauss(A,B) by a proper PA = LU + (independent)
>   back-substitution to solve triangular linear systems associated to
>   B. This would allow to implement matinverseimage directly (more
>   efficiently for fixed A and varying B), without reducing to ker().

Prev by Date: Re: third release candidate for 2.7.0
Previous by thread: RgXn class function for power series
Index(es):
- Date
- Thread