| Peter Bruin on Thu, 12 Jan 2017 16:56:04 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| FFM_gauss, FFM_invimage and analogous functions |
Bonjour, A long time ago I implemented *_invimage functions for finite fields using the bb_field framework, but so far these have not been included in the PARI library. During the Atelier PARI/GP I discovered that also the *_gauss functions for finite fields of characteristic 2 and the FFM_* variants were still missing. I made a patch containing implementations of the following tests: F2xqM_F2xqC_invimage F2xqM_F2xqC_gauss F2xqM_gauss F2xqM_invimage FlxqM_FlxqC_invimage FlxqM_invimage FqM_FqC_invimage FqM_invimage FFM_FFC_gauss FFM_FFC_invimage FFM_gauss FFM_invimage gen_matcolvinvimage gen_matinvimage There are some new tests that cover all the new code as far as I can see, and an entry in usersch5.tex for each function. Here is an example of the resulting speed improvement. The set-up is as follows: ? t = ffgen(2017^3); ? M = matrix(100, 100, i, j, random(t)); ? N = matrix(100, 100, i, j, random(t)); Before: ? matsolve(M, N); ? ## *** last result computed in 475 ms. ? matinverseimage(M, N); ? ## *** last result computed in 343 ms. After: ? matsolve(M, N); ? ## *** last result computed in 190 ms. ? matinverseimage(M, N); ? ## *** last result computed in 225 ms. These functions simply use classical algorithms. The next step will hopefully be to use LQUP factorisation for all standard linear algebra operations on matrices above a certain size. Thanks, Peter
>From c80c85859df6a86ef618384349af961d582a7120 Mon Sep 17 00:00:00 2001
From: Peter Bruin <P.J.Bruin@math.leidenuniv.nl>
Date: Thu, 12 Jan 2017 16:10:46 +0100
Subject: [PATCH] add FFM_gauss, FFM_invimage, FFM_FFC_gauss, FFM_FFC_invimage
and supporting functions
---
doc/usersch5.tex | 40 +++++++++++
src/basemath/FF.c | 109 ++++++++++++++++++++++-------
src/basemath/alglin1.c | 186 ++++++++++++++++++++++++++++++++++++++++++++++++-
src/headers/paridecl.h | 14 ++++
src/test/32/ff | 23 +++++-
src/test/in/ff | 8 +++
6 files changed, 351 insertions(+), 29 deletions(-)
diff --git a/doc/usersch5.tex b/doc/usersch5.tex
index 5de91e2..b74cd7f 100644
--- a/doc/usersch5.tex
+++ b/doc/usersch5.tex
@@ -3902,6 +3902,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}
@@ -3911,6 +3913,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}
@@ -4229,6 +4233,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}
@@ -4239,6 +4245,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}
@@ -5891,16 +5899,24 @@ $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_gauss}{GEN a, GEN b, GEN T}
+
+\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}
\fun{GEN}{F2xqM_det}{GEN a, GEN T}
+\fun{GEN}{F2xqM_gauss}{GEN a, GEN b, GEN T}
+
\fun{GEN}{F2xqM_image}{GEN x, GEN T}
\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}
@@ -9149,10 +9165,14 @@ operate on black box 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}
@@ -11196,6 +11216,16 @@ 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_gauss}{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}) such that $M$ is invertible,
+return the unique column vector $X$ such that $MX=C$.
+
+\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}).
@@ -11206,10 +11236,20 @@ by \tet{RgM_is_FFM(M,\&ff)}).
\fun{GEN}{FFM_det}{GEN M, GEN ff}
+\fun{GEN}{FFM_gauss}{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}) such that $M$ is invertible, return the unique matrix
+$X$ such that $MX=N$.
+
\fun{GEN}{FFM_image}{GEN M, GEN 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 887bfe7..c0165ca 100644
--- a/src/basemath/FF.c
+++ b/src/basemath/FF.c
@@ -1817,6 +1817,7 @@ FqM_to_FpXQM(GEN x, GEN T, GEN p)
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),
@@ -1836,6 +1837,56 @@ FFM_wrap(GEN M, GEN ff, GEN (*Fq)(GEN,GEN,GEN),
if (!M) { avma = av; return NULL; }
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;
+ int is_sqr = M==N;
+ _getFF(ff, &T, &p, &pp);
+ M = FFM_to_raw(M);
+ N = is_sqr? M: FFM_to_raw(N);
+ switch(ff[1])
+ {
+ case t_FF_FpXQ: M = Fq(M, N, T, p); if (M) M = FqM_to_FpXQM(M, T, p);
+ break;
+ case t_FF_F2xq: M = F2xq(M, N, T); break;
+ default: M = Flxq(M, N, T, pp); break;
+ }
+ if (!M) { avma = av; 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 = Fq(M, C, T, p); if (C) C = FqC_to_FpXQC(C, T, p);
+ break;
+ case t_FF_F2xq: C = F2xq(M, C, T); break;
+ default: C = Flxq(M, C, T, pp); break;
+ }
+ if (!C) { avma = av; 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); }
@@ -1878,38 +1929,42 @@ FFM_det(GEN M, GEN ff)
}
GEN
+FFM_FFC_gauss(GEN M, GEN C, GEN ff)
+{
+ return FFM_FFC_wrap(M, C, ff, FqM_FqC_gauss,
+ FlxqM_FlxqC_gauss, F2xqM_F2xqC_gauss);
+}
+
+GEN
+FFM_gauss(GEN M, GEN N, GEN ff)
+{
+ return FFM_FFM_wrap(M, N, ff, FqM_gauss,
+ FlxqM_gauss, F2xqM_gauss);
+}
+
+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;
- int is_sqr = M==N;
- _getFF(ff, &T, &p, &pp);
- M = FFM_to_raw(M);
- N = is_sqr? M: 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 b31b2fe..4f451ec 100644
--- a/src/basemath/alglin1.c
+++ b/src/basemath/alglin1.c
@@ -410,6 +410,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--, j--) {
+ 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;
@@ -1024,6 +1125,13 @@ FlxqM_det(GEN a, GEN T, ulong p)
}
GEN
+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);
@@ -1045,6 +1153,13 @@ FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
}
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;
@@ -1059,6 +1174,19 @@ F2xqM_gauss_gen(GEN a, GEN b, GEN T)
const struct bb_field *S = get_F2xq_field(&E, T);
return gen_Gauss(a, b, E, S);
}
+
+GEN
+F2xqM_gauss(GEN a, GEN b, GEN T)
+{
+ pari_sp av = avma;
+ long n = lg(a)-1;
+ GEN u;
+ if (!n || lg(b)==1) { avma = av; return cgetg(1, t_MAT); }
+ u = F2xqM_gauss_gen(a, b, T);
+ if (!u) { avma = av; return NULL; }
+ return gerepilecopy(av, u);
+}
+
GEN
F2xqM_inv(GEN a, GEN T)
{
@@ -1071,6 +1199,24 @@ F2xqM_inv(GEN a, GEN T)
}
GEN
+F2xqM_F2xqC_gauss(GEN a, GEN b, GEN T)
+{
+ pari_sp av = avma;
+ GEN u;
+ if (lg(a) == 1) return cgetg(1, t_COL);
+ u = F2xqM_gauss_gen(a, mkmat(b), T);
+ if (!u) { avma = av; return NULL; }
+ return gerepilecopy(av, gel(u,1));
+}
+
+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);
@@ -1084,6 +1230,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)
{
@@ -1133,6 +1286,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);
@@ -1146,6 +1306,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)
{
@@ -1597,7 +1764,14 @@ RgM_solve(GEN a, GEN b)
GEN p, u, data, ff = NULL;
if (is_modular_solve(a,b,&u)) return gerepileupto(av, u);
- if (!b && RgM_is_FFM(a, &ff)) return FFM_inv(a, ff);
+ if (RgM_is_FFM(a, &ff)) {
+ if (!b)
+ return FFM_inv(a, ff);
+ if (typ(b) == t_COL && RgC_is_FFC(b, &ff))
+ return FFM_FFC_gauss(a, b, ff);
+ if (typ(b) == t_MAT && RgM_is_FFM(b, &ff))
+ return FFM_gauss(a, b, ff);
+ }
avma = av;
if (lg(a)-1 == 2 && nbrows(a) == 2) {
@@ -3048,6 +3222,11 @@ RgM_RgC_invimage(GEN A, GEN y)
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; }
@@ -3291,6 +3470,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 9cd3290..ae978bc 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -1347,11 +1347,15 @@ 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_gauss(GEN a, GEN b, GEN T);
+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_gauss(GEN a, GEN b, 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);
@@ -1370,12 +1374,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);
@@ -1392,6 +1398,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);
@@ -1399,6 +1406,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);
@@ -1438,7 +1446,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);
@@ -2934,10 +2944,14 @@ 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_gauss(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_gauss(GEN M, GEN N, 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 65c979c..cf08fa6 100644
--- a/src/test/32/ff
+++ b/src/test/32/ff
@@ -187,7 +187,7 @@ Mod(6687681666819568403, 18446744073944432641)
? conjvec(Mod(x,x^2+Mod(1,3)))
[Mod(Mod(1, 3)*x, Mod(1, 3)*x^2 + Mod(1, 3)), Mod(Mod(2, 3)*x, Mod(1, 3)*x^2
+ Mod(1, 3))]~
-? 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);my(N=t*[0,1;0,0]);iferr(N^-1,e,,errname(e)=="e_INV");
+? 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);print(M*matsolve(M,v)==v);print(M*matinverseimage(M,v)==v);print(matsolve(M,matid(3)*t^0)==M^(-1));print(matinverseimage(M,matid(3)*t^0)==M^(-1));my(N=t*[0,1;0,0]);iferr(N^-1,e,,errname(e)=="e_INV");iferr(matsolve(N,t*[0,1]~),e,,errname(e)=="e_INV");print(matinverseimage(N,t*[1,0]~));print(matinverseimage(N,t*[0,1]~));print(matinverseimage(N,N^0));
? test(2^5)
[t^4 + t^3; t^4 + t^3; 1]
[t, t^2; t + 1, t^2 + 1]
@@ -195,6 +195,13 @@ Mod(6687681666819568403, 18446744073944432641)
t^4 + t^2
[1, 0, 0; 0, 1, 0; 0, 0, 1]
[t^2 + t, t^4 + t^3 + 1, t^4 + t^3 + t]~
+1
+1
+1
+1
+[0, 1]~
+[]~
+[;]
? 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]
@@ -202,6 +209,13 @@ t^4 + t^2
6*t^4 + 2*t^3 + 4*t^2 + 2*t + 2
[1, 0, 0; 0, 1, 0; 0, 0, 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]~
+1
+1
+1
+1
+[0, 1]~
+[]~
+[;]
? 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]
@@ -212,6 +226,13 @@ t^4 + t^2
[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]~
+1
+1
+1
+1
+[0, 1]~
+[]~
+[;]
? test(q)=my(t=ffgen(q,'t),M=matrix(10,10,i,j,random(t)));subst(charpoly(M),'x,M)==0;
? test(nextprime(2^7)^5)
1
diff --git a/src/test/in/ff b/src/test/in/ff
index 7e18974..94bd6ab 100644
--- a/src/test/in/ff
+++ b/src/test/in/ff
@@ -72,8 +72,16 @@ test(q)=
print(M^(-1)*M);
my(v = [t^0, t^1, t^2]~);
print(M*v);
+ print(M*matsolve(M, v) == v);
+ print(M*matinverseimage(M, v) == v);
+ print(matsolve(M, matid(3)*t^0) == M^(-1));
+ print(matinverseimage(M, matid(3)*t^0) == M^(-1));
my(N = t*[0, 1; 0, 0]);
iferr(N^-1, e,, errname(e) == "e_INV");
+ iferr(matsolve(N, t*[0, 1]~), e,, errname(e) == "e_INV");
+ print(matinverseimage(N, t*[1, 0]~));
+ print(matinverseimage(N, t*[0, 1]~));
+ print(matinverseimage(N, N^0));
}
test(2^5)
test(7^5)
--
2.7.3