| Peter Bruin on Fri, 10 Feb 2017 17:02:30 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| new functions FFM_deplin, FFM_indexrank, FFM_suppl etc. |
Bonjour,
I am attaching a patch to add implementations of the following functions
using the bb_field framework:
FFM_deplin, FFM_indexrank, FFM_suppl
F2xqM_deplin, F2xqM_indexrank, F2xqM_suppl
FlxqM_deplin, FlxqM_indexrank, FlxqM_suppl
FqM_indexrank
This increases the range of basic linear algebra functions having
specialised implementations for all finite field types. (Note that
FqM_deplin and FqM_suppl already exist.)
The new functions are very straightforward; I did modify the logic of
indexrank very slightly. I also added a few tests so that the new code
is completely covered by the testsuite.
Here is an example showing the resulting speedup. Setup:
? t = ffgen(2017^3);
? M = matrix(120, 80, i, j, random(t)) * matrix(80, 160, i, j, random(t));
? default(timer, 1);
Before:
? matindexrank(M);
time = 272 ms.
? lindep(M);
*** lindep: Warning: increasing stack size to 64000000.
*** lindep: Warning: increasing stack size to 128000000.
*** lindep: Warning: increasing stack size to 256000000.
time = 208 ms.
? matsupplement(M);
time = 256 ms.
After:
? matindexrank(M);
time = 104 ms.
? lindep(M);
time = 152 ms.
? matsupplement(M);
time = 100 ms.
Thanks,
Peter
>From f0f65853f1422c2924fb7a6e2a577d1a4eada063 Mon Sep 17 00:00:00 2001
From: Peter Bruin <P.J.Bruin@math.leidenuniv.nl>
Date: Fri, 10 Feb 2017 11:16:23 +0100
Subject: [PATCH] add FFM_deplin, FFM_indexrank, FFM_suppl and supporting
functions
---
src/basemath/FF.c | 36 +++++++++++++++
src/basemath/alglin1.c | 119 ++++++++++++++++++++++++++++++++++++++++++++----
src/headers/paridecl.h | 10 ++++
src/test/32/ff | 11 ++++-
src/test/in/ff | 3 ++
5 files changed, 168 insertions(+), 11 deletions(-)
diff --git a/src/basemath/FF.c b/src/basemath/FF.c
index c0165ca..d4d7a4a 100644
--- a/src/basemath/FF.c
+++ b/src/basemath/FF.c
@@ -1896,6 +1896,42 @@ FFM_image(GEN M, GEN ff)
GEN
FFM_inv(GEN M, GEN ff)
{ return FFM_wrap(M,ff, &FqM_inv,&FlxqM_inv,&F2xqM_inv); }
+GEN
+FFM_suppl(GEN M, GEN ff)
+{ return FFM_wrap(M,ff, &FqM_suppl,&FlxqM_suppl,&F2xqM_suppl); }
+
+GEN
+FFM_deplin(GEN M, GEN ff)
+{
+ pari_sp av = avma;
+ ulong pp;
+ GEN C, T, p;
+ _getFF(ff, &T, &p, &pp); M = FFM_to_raw(M);
+ switch(ff[1]) {
+ case t_FF_FpXQ: C = FqM_deplin(M, T, p);
+ if (C) C = FqC_to_FpXQC(C, T, p); break;
+ case t_FF_F2xq: C = F2xqM_deplin(M, T); break;
+ default: C = FlxqM_deplin(M, T, pp); break;
+ }
+ if (!C) { avma = av; return NULL; }
+ return gerepilecopy(av, raw_to_FFC(C, ff));
+}
+
+GEN
+FFM_indexrank(GEN M, GEN ff)
+{
+ pari_sp av = avma;
+ ulong pp;
+ GEN R, T, p;
+ _getFF(ff,&T,&p,&pp); M = FFM_to_raw(M);
+ switch(ff[1]) {
+ case t_FF_FpXQ: R = FqM_indexrank(M,T,p); break;
+ case t_FF_F2xq: R = F2xqM_indexrank(M,T); break;
+ default: R = FlxqM_indexrank(M,T,pp); break;
+ }
+ return gerepileupto(av, R);
+}
+
long
FFM_rank(GEN M, GEN ff)
{
diff --git a/src/basemath/alglin1.c b/src/basemath/alglin1.c
index af888c9..e64be2b 100644
--- a/src/basemath/alglin1.c
+++ b/src/basemath/alglin1.c
@@ -1399,6 +1399,12 @@ FlxqM_ker(GEN x, GEN T, ulong p)
return FlxqM_ker_i(x, T, p, 0);
}
+GEN
+FlxqM_deplin(GEN x, GEN T, ulong p)
+{
+ return FlxqM_ker_i(x, T, p, 1);
+}
+
static GEN
F2xqM_ker_i(GEN x, GEN T, long deplin)
{
@@ -1415,6 +1421,13 @@ F2xqM_ker(GEN x, GEN T)
{
return F2xqM_ker_i(x, T, 0);
}
+
+GEN
+F2xqM_deplin(GEN x, GEN T)
+{
+ return F2xqM_ker_i(x, T, 1);
+}
+
static GEN
F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
{
@@ -2796,7 +2809,7 @@ RgV_deplin(GEN v)
GEN
deplin(GEN x)
{
- GEN p = NULL;
+ GEN p = NULL, ff = NULL;
switch(typ(x))
{
case t_MAT: break;
@@ -2828,6 +2841,7 @@ deplin(GEN x)
}
return gerepileupto(av, x);
}
+ if (RgM_is_FFM(x, &ff)) return FFM_deplin(x, ff);
return deplin_aux(x);
}
/*******************************************************************/
@@ -3462,7 +3476,7 @@ GEN
suppl(GEN x)
{
pari_sp av = avma;
- GEN d, X = x, p = NULL;
+ GEN d, X = x, p = NULL, ff = NULL;
long r;
if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
@@ -3478,6 +3492,7 @@ suppl(GEN x)
}
return gerepileupto(av, x);
}
+ if (RgM_is_FFM(x, &ff)) return FFM_suppl(x, ff);
avma = av; init_suppl(x);
d = gauss_pivot(X,&r);
avma = av; return get_suppl(X,d,nbrows(X),r,&col_ei);
@@ -3511,6 +3526,59 @@ F2m_suppl(GEN x)
avma = av; return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
}
+/* variable number to be filled in later */
+static GEN
+_FlxC_ei(long n, long i)
+{
+ GEN x = cgetg(n + 1, t_COL);
+ long j;
+ for (j = 1; j <= n; j++)
+ gel(x, j) = (j == i)? pol1_Flx(0): pol0_Flx(0);
+ return x;
+}
+
+GEN
+F2xqM_suppl(GEN x, GEN T)
+{
+ pari_sp av = avma;
+ GEN d, y;
+ long n = nbrows(x), r, sv = get_Flx_var(T);
+
+ init_suppl(x);
+ d = F2xqM_gauss_pivot(x, T, &r);
+ avma = av;
+ y = get_suppl(x, d, n, r, &_FlxC_ei);
+ if (sv) {
+ long i, j;
+ for (j = r + 1; j <= n; j++) {
+ for (i = 1; i <= n; i++)
+ gcoeff(y, i, j)[1] = sv;
+ }
+ }
+ return y;
+}
+
+GEN
+FlxqM_suppl(GEN x, GEN T, ulong p)
+{
+ pari_sp av = avma;
+ GEN d, y;
+ long n = nbrows(x), r, sv = get_Flx_var(T);
+
+ init_suppl(x);
+ d = FlxqM_gauss_pivot(x, T, p, &r);
+ avma = av;
+ y = get_suppl(x, d, n, r, &_FlxC_ei);
+ if (sv) {
+ long i, j;
+ for (j = r + 1; j <= n; j++) {
+ for (i = 1; i <= n; i++)
+ gcoeff(y, i, j)[1] = sv;
+ }
+ }
+ return y;
+}
+
GEN
FqM_suppl(GEN x, GEN T, GEN p)
{
@@ -3639,24 +3707,25 @@ init_indexrank(GEN x) {
GEN
indexrank(GEN x) {
- pari_sp av = avma;
+ pari_sp av;
long r;
- GEN d, p = NULL;
+ GEN d, p = NULL, ff = NULL;
if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
- init_indexrank(x);
if (RgM_is_FpM(x, &p) && p)
{
ulong pp;
x = RgM_Fp_init(x,p,&pp);
switch(pp)
{
- case 0: d = FpM_gauss_pivot(x,p,&r); break;
- case 2: d = F2m_gauss_pivot(x,&r); break;
- default:d = Flm_gauss_pivot(x,pp,&r); break;
+ case 0: return FpM_indexrank(x, p);
+ case 2: return F2m_indexrank(x);
+ default: return Flm_indexrank(x, pp);
}
}
- else
- d = gauss_pivot(x,&r);
+ if (RgM_is_FFM(x, &ff)) return FFM_indexrank(x, ff);
+ av = avma;
+ init_indexrank(x);
+ d = gauss_pivot(x, &r);
avma = av; return indexrank0(lg(x)-1, r, d);
}
@@ -3689,6 +3758,36 @@ F2m_indexrank(GEN x) {
}
GEN
+F2xqM_indexrank(GEN x, GEN T) {
+ pari_sp av = avma;
+ long r;
+ GEN d;
+ init_indexrank(x);
+ d = F2xqM_gauss_pivot(x, T, &r);
+ avma = av; return indexrank0(lg(x) - 1, r, d);
+}
+
+GEN
+FlxqM_indexrank(GEN x, GEN T, ulong p) {
+ pari_sp av = avma;
+ long r;
+ GEN d;
+ init_indexrank(x);
+ d = FlxqM_gauss_pivot(x, T, p, &r);
+ avma = av; return indexrank0(lg(x) - 1, r, d);
+}
+
+GEN
+FqM_indexrank(GEN x, GEN T, GEN p) {
+ pari_sp av = avma;
+ long r;
+ GEN d;
+ init_indexrank(x);
+ d = FqM_gauss_pivot(x, T, p, &r);
+ avma = av; return indexrank0(lg(x) - 1, r, d);
+}
+
+GEN
ZM_indeximage(GEN x) {
pari_sp av = avma;
long r;
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
index 0660119..2fc07ff 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -1360,14 +1360,17 @@ 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_deplin(GEN x, 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_indexrank(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 F2xqM_suppl(GEN x, GEN T);
GEN Flm_Flc_gauss(GEN a, GEN b, ulong p);
GEN Flm_Flc_invimage(GEN mat, GEN y, ulong p);
GEN Flm_deplin(GEN x, ulong p);
@@ -1386,14 +1389,17 @@ 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_deplin(GEN x, 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_indexrank(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 FlxqM_suppl(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);
GEN FpM_deplin(GEN x, GEN p);
@@ -1415,6 +1421,7 @@ 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_indexrank(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);
@@ -2969,14 +2976,17 @@ 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_deplin(GEN M, 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_indexrank(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);
+GEN FFM_suppl(GEN M, GEN ff);
GEN FFX_factor(GEN f, GEN x);
GEN FFX_roots(GEN f, GEN x);
GEN FqX_to_FFX(GEN x, GEN ff);
diff --git a/src/test/32/ff b/src/test/32/ff
index cf08fa6..101e9dc 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);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(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));print(matindexrank(N));print(matsupplement(t*[0;1]));print(lindep(t*[0,0;1,1]));
? test(2^5)
[t^4 + t^3; t^4 + t^3; 1]
[t, t^2; t + 1, t^2 + 1]
@@ -202,6 +202,9 @@ t^4 + t^2
[0, 1]~
[]~
[;]
+[Vecsmall([1]), Vecsmall([2])]
+[0, 1; t, 0]
+[1, 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]
@@ -216,6 +219,9 @@ t^4 + t^2
[0, 1]~
[]~
[;]
+[Vecsmall([1]), Vecsmall([2])]
+[0, 1; t, 0]
+[6, 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]
@@ -233,6 +239,9 @@ t^4 + t^2
[0, 1]~
[]~
[;]
+[Vecsmall([1]), Vecsmall([2])]
+[0, 1; t, 0]
+[18446744073709551628, 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 94bd6ab..ae4555e 100644
--- a/src/test/in/ff
+++ b/src/test/in/ff
@@ -82,6 +82,9 @@ test(q)=
print(matinverseimage(N, t*[1, 0]~));
print(matinverseimage(N, t*[0, 1]~));
print(matinverseimage(N, N^0));
+ print(matindexrank(N));
+ print(matsupplement(t*[0; 1]));
+ print(lindep(t*[0, 0; 1, 1]));
}
test(2^5)
test(7^5)
--
1.7.9.5