version 1.1, 2001/10/02 11:17:00 |
version 1.2, 2002/09/11 07:26:48 |
Line 23 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 23 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
|
|
/* for linear algebra mod p */ |
/* for linear algebra mod p */ |
#ifdef LONG_IS_64BIT |
#ifdef LONG_IS_64BIT |
# define MASK (0x7fffffff00000000UL) |
# define MASK (0xffffffff00000000UL) |
#else |
#else |
# define MASK (0x7fff0000UL) |
# define MASK (0xffff0000UL) |
#endif |
#endif |
|
|
/* 2p^2 < 2^BITS_IN_LONG */ |
/* 2p^2 < 2^BITS_IN_LONG */ |
Line 36 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 36 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
#endif |
#endif |
#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2])) |
#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2])) |
extern GEN ZM_init_CRT(GEN Hp, ulong p); |
extern GEN ZM_init_CRT(GEN Hp, ulong p); |
extern ulong u_invmod(ulong x, ulong p); |
extern int ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p); |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
|
|
return y; |
return y; |
} |
} |
|
|
|
|
GEN |
GEN |
gtrans(GEN x) |
gtrans(GEN x) |
{ |
{ |
|
|
concatsp3(GEN x, GEN y, GEN z) |
concatsp3(GEN x, GEN y, GEN z) |
{ |
{ |
long i, lx = lg(x), ly = lg(y), lz = lg(z); |
long i, lx = lg(x), ly = lg(y), lz = lg(z); |
GEN r = cgetg(lx+ly+lz-2, t_MAT), t = r; |
GEN t, r = cgetg(lx+ly+lz-2, t_MAT); |
|
t = r; |
for (i=1; i<lx; i++) *++t = *++x; |
for (i=1; i<lx; i++) *++t = *++x; |
for (i=1; i<ly; i++) *++t = *++y; |
for (i=1; i<ly; i++) *++t = *++y; |
for (i=1; i<lz; i++) *++t = *++z; |
for (i=1; i<lz; i++) *++t = *++z; |
Line 149 vconcat(GEN A, GEN B) |
|
Line 149 vconcat(GEN A, GEN B) |
|
long la,ha,hb,hc,i,j; |
long la,ha,hb,hc,i,j; |
GEN M,a,b,c; |
GEN M,a,b,c; |
|
|
|
if (!A) return B; |
|
if (!B) return A; |
la = lg(A); if (la==1) return A; |
la = lg(A); if (la==1) return A; |
ha = lg(A[1]); M = cgetg(la,t_MAT); |
ha = lg(A[1]); M = cgetg(la,t_MAT); |
hb = lg(B[1]); hc = ha+hb-1; |
hb = lg(B[1]); hc = ha+hb-1; |
Line 162 vconcat(GEN A, GEN B) |
|
Line 164 vconcat(GEN A, GEN B) |
|
} |
} |
|
|
GEN |
GEN |
|
_veccopy(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = lcopy(x); return v; } |
|
GEN |
_vec(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = (long)x; return v; } |
_vec(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = (long)x; return v; } |
GEN |
GEN |
_col(GEN x) { GEN v = cgetg(2, t_COL); v[1] = (long)x; return v; } |
_col(GEN x) { GEN v = cgetg(2, t_COL); v[1] = (long)x; return v; } |
|
|
|
/* routines for naive growarrays */ |
GEN |
GEN |
|
cget1(long l, long t) |
|
{ |
|
GEN z = new_chunk(l); |
|
z[0] = evaltyp(t) | evallg(1); return z; |
|
} |
|
void |
|
appendL(GEN x, GEN t) |
|
{ |
|
long l = lg(x); x[l] = (long)t; setlg(x, l+1); |
|
} |
|
|
|
GEN |
concatsp(GEN x, GEN y) |
concatsp(GEN x, GEN y) |
{ |
{ |
long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i; |
long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i; |
Line 283 concat(GEN x, GEN y) |
|
Line 300 concat(GEN x, GEN y) |
|
|
|
if (!y) |
if (!y) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
if (tx == t_LIST) |
if (tx == t_LIST) |
{ lx = lgef(x); i = 2; } |
{ lx = lgef(x); i = 2; } |
else if (tx == t_VEC) |
else if (tx == t_VEC) |
Line 420 get_range(char *s, long *a, long *b, long *cmpl, long |
|
Line 437 get_range(char *s, long *a, long *b, long *cmpl, long |
|
|
|
*a = 1; *b = max; |
*a = 1; *b = max; |
if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0; |
if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0; |
if (*s == 0) return 0; |
if (!*s) return 0; |
if (*s != '.') |
if (*s != '.') |
{ |
{ |
*a = str_to_long(s, &s); |
*a = str_to_long(s, &s); |
Line 461 rowextract_i(GEN A, long x1, long x2) |
|
Line 478 rowextract_i(GEN A, long x1, long x2) |
|
return B; |
return B; |
} |
} |
|
|
|
/* A[x0,] */ |
GEN |
GEN |
|
row(GEN A, long x0) |
|
{ |
|
long i, lB = lg(A); |
|
GEN B = cgetg(lB, t_VEC); |
|
for (i=1; i<lB; i++) B[i] = coeff(A, x0, i); |
|
return B; |
|
} |
|
|
|
/* A[x0, x1..x2] */ |
|
GEN |
|
row_i(GEN A, long x0, long x1, long x2) |
|
{ |
|
long i, lB = x2 - x1 + 2; |
|
GEN B = cgetg(lB, t_VEC); |
|
for (i=x1; i<=x2; i++) B[i] = coeff(A, x0, i); |
|
return B; |
|
} |
|
|
|
GEN |
vecextract_p(GEN A, GEN p) |
vecextract_p(GEN A, GEN p) |
{ |
{ |
long i,lB = lg(p); |
long i,lB = lg(p); |
Line 499 rowselect_p(GEN A, GEN B, GEN p, long init) |
|
Line 536 rowselect_p(GEN A, GEN B, GEN p, long init) |
|
GEN |
GEN |
extract(GEN x, GEN l) |
extract(GEN x, GEN l) |
{ |
{ |
long av,i,j, tl = typ(l), tx = typ(x), lx = lg(x); |
gpmem_t av; |
|
long i,j, tl = typ(l), tx = typ(x), lx = lg(x); |
GEN y; |
GEN y; |
|
|
if (! is_matvec_t(tx)) err(typeer,"extract"); |
if (! is_matvec_t(tx)) err(typeer,"extract"); |
Line 585 extract(GEN x, GEN l) |
|
Line 623 extract(GEN x, GEN l) |
|
GEN |
GEN |
matextract(GEN x, GEN l1, GEN l2) |
matextract(GEN x, GEN l1, GEN l2) |
{ |
{ |
long av = avma, tetpil; |
gpmem_t av = avma, tetpil; |
|
|
if (typ(x)!=t_MAT) err(typeer,"matextract"); |
if (typ(x)!=t_MAT) err(typeer,"matextract"); |
x = extract(gtrans(extract(x,l2)),l1); tetpil=avma; |
x = extract(gtrans(extract(x,l2)),l1); tetpil=avma; |
Line 599 extract0(GEN x, GEN l1, GEN l2) |
|
Line 637 extract0(GEN x, GEN l1, GEN l2) |
|
return matextract(x,l1,l2); |
return matextract(x,l1,l2); |
} |
} |
|
|
|
/* v[a] + ... + v[b] */ |
|
GEN |
|
sum(GEN v, long a, long b) |
|
{ |
|
GEN p; |
|
long i; |
|
if (a > b) return gzero; |
|
if (b > lg(v)-1) err(talker,"incorrect length in sum"); |
|
p = (GEN)v[a]; |
|
for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]); |
|
return p; |
|
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* SCALAR-MATRIX OPERATIONS */ |
/* SCALAR-MATRIX OPERATIONS */ |
|
|
zeromat(long m, long n) |
zeromat(long m, long n) |
{ |
{ |
GEN y = cgetg(n+1,t_MAT); |
GEN y = cgetg(n+1,t_MAT); |
long i; for (i=1; i<=n; i++) y[i]=(long)zerocol(m); |
GEN v = zerocol(m); |
|
long i; for (i=1; i<=n; i++) y[i]=(long)v; |
return y; |
return y; |
} |
} |
|
|
GEN |
GEN |
gscalcol(GEN x, long n) |
gscalcol(GEN x, long n) |
{ |
{ |
GEN y=gscalcol_proto(gzero,gzero,n); |
GEN y=gscalcol_proto(gzero,gzero,n); |
if (n) y[1]=lcopy(x); |
if (n) y[1]=lcopy(x); |
return y; |
return y; |
} |
} |
|
|
Line 847 mattodiagonal(GEN m) |
|
Line 899 mattodiagonal(GEN m) |
|
GEN |
GEN |
gaddmat(GEN x, GEN y) |
gaddmat(GEN x, GEN y) |
{ |
{ |
long ly,dy,i,j; |
long l,d,i,j; |
GEN z; |
GEN z, cz,cy; |
|
|
ly=lg(y); if (ly==1) return cgetg(1,t_MAT); |
l = lg(y); if (l==1) return cgetg(1,t_MAT); |
dy=lg(y[1]); |
d = lg(y[1]); |
if (typ(y)!=t_MAT || ly!=dy) err(mattype1,"gaddmat"); |
if (typ(y)!=t_MAT || l!=d) err(mattype1,"gaddmat"); |
z=cgetg(ly,t_MAT); |
z=cgetg(l,t_MAT); |
for (i=1; i<ly; i++) |
for (i=1; i<l; i++) |
{ |
{ |
z[i]=lgetg(dy,t_COL); |
cz = cgetg(d,t_COL); z[i] = (long)cz; cy = (GEN)y[i]; |
for (j=1; j<dy; j++) |
for (j=1; j<d; j++) |
coeff(z,j,i) = i==j? ladd(x,gcoeff(y,j,i)): lcopy(gcoeff(y,j,i)); |
cz[j] = i==j? ladd(x,(GEN)cy[j]): lcopy((GEN)cy[j]); |
} |
} |
return z; |
return z; |
} |
} |
|
|
|
/* same, no copy */ |
|
GEN |
|
gaddmat_i(GEN x, GEN y) |
|
{ |
|
long l,d,i,j; |
|
GEN z, cz,cy; |
|
|
|
l = lg(y); if (l==1) return cgetg(1,t_MAT); |
|
d = lg(y[1]); |
|
if (typ(y)!=t_MAT || l!=d) err(mattype1,"gaddmat"); |
|
z = cgetg(l,t_MAT); |
|
for (i=1; i<l; i++) |
|
{ |
|
cz = cgetg(d,t_COL); z[i] = (long)cz; cy = (GEN)y[i]; |
|
for (j=1; j<d; j++) |
|
cz[j] = i==j? ladd(x,(GEN)cy[j]): cy[j]; |
|
} |
|
return z; |
|
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* Solve A*X=B (Gauss pivot) */ |
/* Solve A*X=B (Gauss pivot) */ |
Line 871 gaddmat(GEN x, GEN y) |
|
Line 943 gaddmat(GEN x, GEN y) |
|
#define swap(x,y) { long _t=x; x=y; y=_t; } |
#define swap(x,y) { long _t=x; x=y; y=_t; } |
|
|
/* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be |
/* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be |
* used, and the matrix precision (min real precision of coeffs) otherwise. |
* used, 1 otherwise */ |
*/ |
static int |
static long |
use_maximal_pivot(GEN x) |
matprec(GEN x) |
|
{ |
{ |
long tx,i,j,l, k = VERYBIGINT, lx = lg(x), ly = lg(x[1]); |
int res = 0; |
|
long tx,i,j, lx = lg(x), ly = lg(x[1]); |
GEN p1; |
GEN p1; |
for (i=1; i<lx; i++) |
for (i=1; i<lx; i++) |
for (j=1; j<ly; j++) |
for (j=1; j<ly; j++) |
{ |
{ |
p1 = gmael(x,i,j); tx = typ(p1); |
p1 = gmael(x,i,j); tx = typ(p1); |
if (!is_scalar_t(tx)) return 0; |
if (!is_scalar_t(tx)) return 0; |
l = precision(p1); if (l && l<k) k = l; |
if (precision(p1)) res = 1; |
} |
} |
return (k==VERYBIGINT)? 0: k; |
return res; |
} |
} |
|
|
/* As above, returning 1 if the precision would be non-zero, 0 otherwise */ |
static GEN |
static long |
col_to_mat(GEN b) |
use_maximal_pivot(GEN x) |
|
{ |
{ |
long tx,i,j, lx = lg(x), ly = lg(x[1]); |
GEN B = cgetg(2, t_MAT); |
GEN p1; |
B[1] = (long)b; return B; |
for (i=1; i<lx; i++) |
|
for (j=1; j<ly; j++) |
|
{ |
|
p1 = gmael(x,i,j); tx = typ(p1); |
|
if (!is_scalar_t(tx)) return 0; |
|
if (precision(p1)) return 1; |
|
} |
|
return 0; |
|
} |
} |
|
|
static GEN |
static GEN |
check_b(GEN b, long nbli) |
check_b(GEN b, long nbli, int *iscol) |
{ |
{ |
GEN col; |
GEN col; |
|
*iscol = (b && typ(b) == t_COL); |
if (!b) return idmat(nbli); |
if (!b) return idmat(nbli); |
col = (typ(b) == t_MAT)? (GEN)b[1]: b; |
b = dummycopy(b); |
|
if (*iscol) { col = b; b = col_to_mat(b); } |
|
else |
|
{ |
|
if (lg(b) == 1) err(consister,"gauss"); |
|
col = (GEN)b[1]; |
|
} |
if (nbli != lg(col)-1) err(talker,"incompatible matrix dimensions in gauss"); |
if (nbli != lg(col)-1) err(talker,"incompatible matrix dimensions in gauss"); |
return dummycopy(b); |
return b; |
} |
} |
|
|
/* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */ |
/* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */ |
Line 925 gauss_triangle_i(GEN A, GEN B, GEN t) |
|
Line 995 gauss_triangle_i(GEN A, GEN B, GEN t) |
|
for (k=1; k<=n; k++) |
for (k=1; k<=n; k++) |
{ |
{ |
GEN u = cgetg(n+1, t_COL), b = (GEN)B[k]; |
GEN u = cgetg(n+1, t_COL), b = (GEN)B[k]; |
ulong av = avma; |
gpmem_t av = avma; |
c[k] = (long)u; m = mulii((GEN)b[n],t); |
c[k] = (long)u; m = mulii((GEN)b[n],t); |
u[n] = (long)gerepileuptoint(av, divii(m, gcoeff(A,n,n))); |
u[n] = lpileuptoint(av, divii(m, gcoeff(A,n,n))); |
for (i=n-1; i>0; i--) |
for (i=n-1; i>0; i--) |
{ |
{ |
av = avma; m = mulii(negi((GEN)b[i]),t); |
av = avma; m = mulii(negi((GEN)b[i]),t); |
for (j=i+1; j<=n; j++) |
for (j=i+1; j<=n; j++) |
m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); |
m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); |
u[i] = (long)gerepileuptoint(av, divii(negi(m), gcoeff(A,i,i))); |
u[i] = lpileuptoint(av, divii(negi(m), gcoeff(A,i,i))); |
} |
} |
} |
} |
return c; |
return c; |
} |
} |
|
|
|
/* A, B integral upper HNF. A^(-1) B integral ? */ |
|
int |
|
hnfdivide(GEN A, GEN B) |
|
{ |
|
gpmem_t av = avma; |
|
long n = lg(A)-1, i,j,k; |
|
GEN u, b, m, r; |
|
|
|
if (!n) return 1; |
|
u = cgetg(n+1, t_COL); |
|
for (k=1; k<=n; k++) |
|
{ |
|
b = (GEN)B[k]; |
|
m = (GEN)b[k]; |
|
u[k] = ldvmdii(m, gcoeff(A,k,k), &r); |
|
if (r != gzero) { avma = av; return 0; } |
|
for (i=k-1; i>0; i--) |
|
{ |
|
m = negi((GEN)b[i]); |
|
for (j=i+1; j<=k; j++) |
|
m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); |
|
m = dvmdii(m, gcoeff(A,i,i), &r); |
|
if (r != gzero) { avma = av; return 0; } |
|
u[i] = lnegi(m); |
|
} |
|
} |
|
avma = av; return 1; |
|
} |
|
|
|
/* A upper HNF, b integral vector. Return A^(-1) b if integral, |
|
* NULL otherwise. */ |
GEN |
GEN |
|
hnf_invimage(GEN A, GEN b) |
|
{ |
|
gpmem_t av = avma, av2; |
|
long n = lg(A)-1, i,j; |
|
GEN u, m, r; |
|
|
|
if (!n) return NULL; |
|
u = cgetg(n+1, t_COL); |
|
m = (GEN)b[n]; |
|
u[n] = ldvmdii(m, gcoeff(A,n,n), &r); |
|
if (r != gzero) { avma = av; return NULL; } |
|
for (i=n-1; i>0; i--) |
|
{ |
|
av2 = avma; |
|
m = negi((GEN)b[i]); |
|
for (j=i+1; j<=n; j++) |
|
m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); |
|
m = dvmdii(m, gcoeff(A,i,i), &r); |
|
if (r != gzero) { avma = av; return NULL; } |
|
u[i] = lpileuptoint(av2, negi(m)); |
|
} |
|
return u; |
|
} |
|
|
|
GEN |
gauss_get_col(GEN a, GEN b, GEN p, long li) |
gauss_get_col(GEN a, GEN b, GEN p, long li) |
{ |
{ |
GEN m, u=cgetg(li+1,t_COL); |
GEN m, u=cgetg(li+1,t_COL); |
Line 948 gauss_get_col(GEN a, GEN b, GEN p, long li) |
|
Line 1074 gauss_get_col(GEN a, GEN b, GEN p, long li) |
|
u[li] = ldiv((GEN)b[li], p); |
u[li] = ldiv((GEN)b[li], p); |
for (i=li-1; i>0; i--) |
for (i=li-1; i>0; i--) |
{ |
{ |
|
gpmem_t av = avma; |
m = gneg_i((GEN)b[i]); |
m = gneg_i((GEN)b[i]); |
for (j=i+1; j<=li; j++) |
for (j=i+1; j<=li; j++) |
m = gadd(m, gmul(gcoeff(a,i,j), (GEN)u[j])); |
m = gadd(m, gmul(gcoeff(a,i,j), (GEN)u[j])); |
u[i] = ldiv(gneg_i(m), gcoeff(a,i,i)); |
u[i] = lpileupto(av, gdiv(gneg_i(m), gcoeff(a,i,i))); |
} |
} |
return u; |
return u; |
} |
} |
Line 965 Fp_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN p |
|
Line 1092 Fp_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN p |
|
u[li] = lresii(mulii((GEN)b[li], mpinvmod(piv,p)), p); |
u[li] = lresii(mulii((GEN)b[li], mpinvmod(piv,p)), p); |
for (i=li-1; i>0; i--) |
for (i=li-1; i>0; i--) |
{ |
{ |
|
gpmem_t av = avma; |
m = (GEN)b[i]; |
m = (GEN)b[i]; |
for (j=i+1; j<=li; j++) |
for (j=i+1; j<=li; j++) |
m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j])); |
m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j])); |
m = resii(m, p); |
m = resii(m, p); |
u[i] = lresii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p); |
u[i] = lpileuptoint(av, resii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p)); |
} |
} |
return u; |
return u; |
} |
} |
|
GEN |
|
Fq_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN T, GEN p) |
|
{ |
|
GEN m, u=cgetg(li+1,t_COL); |
|
long i,j; |
|
|
|
u[li] = (long)FpXQ_mul((GEN)b[li], FpXQ_inv(piv,T,p), T,p); |
|
for (i=li-1; i>0; i--) |
|
{ |
|
gpmem_t av = avma; |
|
m = (GEN)b[i]; |
|
for (j=i+1; j<=li; j++) |
|
m = gsub(m, gmul(gcoeff(a,i,j), (GEN)u[j])); |
|
m = FpX_res(m, T,p); |
|
u[i] = lpileupto(av, FpXQ_mul(m, FpXQ_inv(gcoeff(a,i,i), T,p), T,p)); |
|
} |
|
return u; |
|
} |
|
|
/* assume -p < a < p, return 1/a mod p */ |
/* assume -p < a < p, return 1/a mod p */ |
static long |
static long |
u_Fp_inv(long a, long p) |
u_Fp_inv(long a, long p) |
{ |
{ |
if (a < 0) a = p + a; /* pb with ulongs < 0 */ |
if (a < 0) a = p + a; /* pb with ulongs < 0 */ |
return u_invmod(a,p); |
return (long)u_invmod((ulong)a,(ulong)p); |
} |
} |
|
|
|
/* assume 0 <= a[i,i] < p */ |
GEN |
GEN |
u_Fp_gauss_get_col(GEN a, GEN b, long piv, long li, long p) |
u_Fp_gauss_get_col(GEN a, GEN b, ulong piv, long li, ulong p) |
{ |
{ |
GEN u=cgetg(li+1,t_VECSMALL); |
GEN u = cgetg(li+1,t_VECSMALL); |
long i,j, m; |
ulong m; |
|
long i,j; |
|
|
u[li] = (b[li] * u_Fp_inv(piv,p)) % p; |
m = b[li] % p; |
if (u[li] < 0) u[li] += p; |
if (u_OK_ULONG(p)) |
for (i=li-1; i>0; i--) |
|
{ |
{ |
m = b[i]; |
u[li] = (m * u_Fp_inv(piv,p)) % p; |
for (j=i+1; j<=li; j++) { m -= coeff(a,i,j) * u[j]; if (m & MASK) m %= p; } |
for (i=li-1; i>0; i--) |
u[i] = ((m%p) * u_Fp_inv(coeff(a,i,i), p)) % p; |
{ |
if (u[i] < 0) u[i] += p; |
m = p - b[i]%p; |
|
for (j=i+1; j<=li; j++) { |
|
if (m & HIGHBIT) m %= p; |
|
m += coeff(a,i,j) * u[j]; |
|
} |
|
m %= p; |
|
if (m) m = ((p-m) * u_invmod((ulong)coeff(a,i,i), p)) % p; |
|
u[i] = m; |
|
} |
} |
} |
|
else |
|
{ |
|
u[li] = mulssmod(m, u_Fp_inv(piv,p), p); |
|
for (i=li-1; i>0; i--) |
|
{ |
|
m = p - b[i]%p; |
|
for (j=i+1; j<=li; j++) |
|
m += mulssmod(coeff(a,i,j), u[j], p); |
|
m %= p; |
|
if (m) m = mulssmod(p-m, u_invmod((ulong)coeff(a,i,i), p), p); |
|
u[i] = m; |
|
} |
|
} |
return u; |
return u; |
} |
} |
|
|
Line 1015 _Fp_addmul(GEN b, long k, long i, GEN m, GEN p) |
|
Line 1183 _Fp_addmul(GEN b, long k, long i, GEN m, GEN p) |
|
b[k] = laddii((GEN)b[k], mulii(m, (GEN)b[i])); |
b[k] = laddii((GEN)b[k], mulii(m, (GEN)b[i])); |
} |
} |
|
|
|
/* same, reduce mod (T,p) */ |
static void |
static void |
_u_Fp_addmul(GEN b, long k, long i, long m, long p) |
_Fq_addmul(GEN b, long k, long i, GEN m, GEN T, GEN p) |
{ |
{ |
long a; |
b[i] = (long)FpX_res((GEN)b[i], T,p); |
if (b[i] & MASK) b[i] %= p; |
b[k] = (long)ladd((GEN)b[k], gmul(m, (GEN)b[i])); |
a = b[k] + m * b[i]; |
|
if (a & MASK) a %= p; |
|
b[k] = a; |
|
} |
} |
|
|
|
/* assume m < p && u_OK_ULONG(p) && (! (b[i] & MASK)) */ |
|
static void |
|
_u_Fp_addmul_OK(GEN b, long k, long i, ulong m, ulong p) |
|
{ |
|
b[k] += m * b[i]; |
|
if (b[k] & MASK) b[k] %= p; |
|
} |
|
/* assume m < p */ |
|
static void |
|
_u_Fp_addmul(GEN b, long k, long i, ulong m, ulong p) |
|
{ |
|
b[i] %= p; |
|
b[k] += mulssmod(m, b[i], p); |
|
if (b[k] & MASK) b[k] %= p; |
|
} |
|
/* same m = 1 */ |
|
static void |
|
_u_Fp_add(GEN b, long k, long i, ulong p) |
|
{ |
|
b[k] += b[i]; |
|
if (b[k] & MASK) b[k] %= p; |
|
} |
|
|
/* Gaussan Elimination. Compute a^(-1)*b |
/* Gaussan Elimination. Compute a^(-1)*b |
* b is a matrix or column vector, NULL meaning: take the identity matrix |
* b is a matrix or column vector, NULL meaning: take the identity matrix |
* If a and b are empty, the result is the empty matrix. |
* If a and b are empty, the result is the empty matrix. |
Line 1037 _u_Fp_addmul(GEN b, long k, long i, long m, long p) |
|
Line 1226 _u_Fp_addmul(GEN b, long k, long i, long m, long p) |
|
GEN |
GEN |
gauss_intern(GEN a, GEN b) |
gauss_intern(GEN a, GEN b) |
{ |
{ |
long inexact,iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1; |
gpmem_t av, lim; |
|
long i,j,k,li,bco, aco = lg(a)-1; |
|
int inexact, iscol; |
GEN p,m,u; |
GEN p,m,u; |
|
|
if (typ(a)!=t_MAT) err(mattype1,"gauss"); |
if (typ(a)!=t_MAT) err(mattype1,"gauss"); |
Line 1048 gauss_intern(GEN a, GEN b) |
|
Line 1239 gauss_intern(GEN a, GEN b) |
|
if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1); |
if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1); |
return cgetg(1,t_MAT); |
return cgetg(1,t_MAT); |
} |
} |
av=avma; lim=stack_lim(av,1); |
av = avma; lim = stack_lim(av,1); |
li = lg(a[1])-1; |
li = lg(a[1])-1; |
if (li != aco && (li < aco || b)) err(mattype1,"gauss"); |
if (li != aco && (li < aco || b)) err(mattype1,"gauss"); |
a = dummycopy(a); |
a = dummycopy(a); |
b = check_b(b,li); bco = lg(b)-1; |
b = check_b(b,li, &iscol); bco = lg(b)-1; |
inexact = use_maximal_pivot(a); |
inexact = use_maximal_pivot(a); |
iscol = (typ(b)==t_COL); |
if(DEBUGLEVEL>4) fprintferr("Entering gauss with inexact=%ld\n",inexact); |
if(DEBUGLEVEL>4) |
|
fprintferr("Entering gauss with inexact=%ld iscol=%ld\n",inexact,iscol); |
|
|
|
p = NULL; /* gcc -Wall */ |
p = NULL; /* gcc -Wall */ |
for (i=1; i<=aco; i++) |
for (i=1; i<=aco; i++) |
Line 1083 gauss_intern(GEN a, GEN b) |
|
Line 1272 gauss_intern(GEN a, GEN b) |
|
if (k != i) |
if (k != i) |
{ |
{ |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
if (iscol) { swap(b[i],b[k]); } |
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
else |
|
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
|
p = gcoeff(a,i,i); |
p = gcoeff(a,i,i); |
} |
} |
if (i == aco) break; |
if (i == aco) break; |
|
|
for (k=i+1; k<=li; k++) |
for (k=i+1; k<=li; k++) |
{ |
{ |
m=gcoeff(a,k,i); |
m = gcoeff(a,k,i); |
if (!gcmp0(m)) |
if (!gcmp0(m)) |
{ |
{ |
m = gneg_i(gdiv(m,p)); |
m = gneg_i(gdiv(m,p)); |
for (j=i+1; j<=aco; j++) _addmul((GEN)a[j],k,i,m); |
for (j=i+1; j<=aco; j++) _addmul((GEN)a[j],k,i,m); |
if (iscol) _addmul(b,k,i,m); |
for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m); |
else |
|
for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m); |
|
} |
} |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b; |
|
if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i); |
if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i); |
gerepilemany(av,gptr,2); |
gerepileall(av,2, &a,&b); |
} |
} |
} |
} |
|
|
if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); |
if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); |
if (iscol) u = gauss_get_col(a,b,p,aco); |
u = cgetg(bco+1,t_MAT); |
else |
for (j=1; j<=bco; j++) u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco); |
{ |
return gerepilecopy(av, iscol? (GEN)u[1]: u); |
long av1 = avma; |
|
lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT); |
|
for (j=2; j<=bco; j++) u[j] = zero; /* dummy */ |
|
for (j=1; j<=bco; j++) |
|
{ |
|
u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco); |
|
if (low_stack(lim, stack_lim(av1,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j); |
|
u = gerepilecopy(av1, u); |
|
} |
|
} |
|
} |
|
return gerepilecopy(av, u); |
|
} |
} |
|
|
GEN |
GEN |
Line 1138 gauss(GEN a, GEN b) |
|
Line 1308 gauss(GEN a, GEN b) |
|
return z; |
return z; |
} |
} |
|
|
GEN |
/* destroy a, b */ |
u_FpM_gauss(GEN a, GEN b, long p) |
static GEN |
|
u_FpM_gauss_sp(GEN a, GEN b, ulong p) |
{ |
{ |
long piv,m,iscol,i,j,k,li,bco, aco = lg(a)-1; |
long iscol,i,j,k,li,bco, aco = lg(a)-1; |
|
ulong piv, m; |
GEN u; |
GEN u; |
|
|
if (!aco) return cgetg(1,t_MAT); |
if (!aco) return cgetg(1,t_MAT); |
li = lg(a[1])-1; |
li = lg(a[1])-1; |
bco = lg(b)-1; |
bco = lg(b)-1; |
iscol = (typ(b)==t_COL); |
iscol = (typ(b)!=t_MAT); |
|
if (iscol) b = col_to_mat(b); |
piv = 0; /* gcc -Wall */ |
piv = 0; /* gcc -Wall */ |
for (i=1; i<=aco; i++) |
for (i=1; i<=aco; i++) |
{ |
{ |
/* k is the line where we find the pivot */ |
/* k is the line where we find the pivot */ |
coeff(a,i,i) = coeff(a,i,i) % p; |
piv = coeff(a,i,i) % p; |
piv = coeff(a,i,i); k = i; |
coeff(a,i,i) = piv; k = i; |
if (!piv) |
if (!piv) |
{ |
{ |
for (k++; k <= li; k++) |
for (k++; k <= li; k++) |
{ |
{ |
coeff(a,k,i) %= p; |
coeff(a,k,i) %= p; |
if (coeff(a,k,i)) break; |
if (coeff(a,k,i)) break; |
} |
} |
if (k>li) return NULL; |
if (k>li) return NULL; |
} |
} |
|
|
Line 1168 u_FpM_gauss(GEN a, GEN b, long p) |
|
Line 1341 u_FpM_gauss(GEN a, GEN b, long p) |
|
if (k != i) |
if (k != i) |
{ |
{ |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
if (iscol) { swap(b[i],b[k]); } |
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
else |
|
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
|
piv = coeff(a,i,i); |
piv = coeff(a,i,i); |
} |
} |
if (i == aco) break; |
if (i == aco) break; |
|
|
for (k=i+1; k<=li; k++) |
for (k=i+1; k<=li; k++) |
{ |
{ |
coeff(a,k,i) %= p; |
m = coeff(a,k,i) % p; coeff(a,k,i) = 0; |
m = coeff(a,k,i); coeff(a,k,i) = 0; |
|
if (m) |
if (m) |
{ |
{ |
m = - (m * u_Fp_inv(piv,p)) % p; |
m = mulssmod(m, u_invmod(piv,p), p); |
for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p); |
m = p - m; |
if (iscol) _u_Fp_addmul(b,k,i,m, p); |
if (u_OK_ULONG(p)) |
else |
{ |
for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p); |
for (j=i+1; j<=aco; j++) _u_Fp_addmul_OK((GEN)a[j],k,i,m, p); |
|
for (j=1; j<=bco; j++) _u_Fp_addmul_OK((GEN)b[j],k,i,m, p); |
|
} |
|
else |
|
{ |
|
for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p); |
|
for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p); |
|
} |
} |
} |
} |
} |
} |
} |
if (iscol) u = u_Fp_gauss_get_col(a,b,piv,aco,p); |
u = cgetg(bco+1,t_MAT); |
else |
for (j=1; j<=bco; j++) u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); |
{ |
return iscol? (GEN)u[1]: u; |
u=cgetg(bco+1,t_MAT); |
|
for (j=1; j<=bco; j++) |
|
u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); |
|
} |
|
return u; |
|
} |
} |
|
|
|
static GEN |
|
u_idmat(long n) |
|
{ |
|
GEN y = cgetg(n+1,t_MAT); |
|
long i; |
|
if (n < 0) err(talker,"negative size in u_idmat"); |
|
for (i=1; i<=n; i++) { y[i] = (long)vecsmall_const(n, 0); coeff(y, i,i) = 1; } |
|
return y; |
|
} |
|
|
GEN |
GEN |
|
u_FpM_gauss(GEN a, GEN b, ulong p) { |
|
return u_FpM_gauss_sp(dummycopy(a), dummycopy(b), p); |
|
} |
|
static GEN |
|
u_FpM_inv_sp(GEN a, ulong p) { |
|
return u_FpM_gauss_sp(a, u_idmat(lg(a)-1), p); |
|
} |
|
GEN |
|
u_FpM_inv(GEN a, ulong p) { |
|
return u_FpM_inv_sp(dummycopy(a), p); |
|
} |
|
|
|
GEN |
FpM_gauss(GEN a, GEN b, GEN p) |
FpM_gauss(GEN a, GEN b, GEN p) |
{ |
{ |
long iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1; |
gpmem_t av,lim; |
|
long i,j,k,li,bco, aco = lg(a)-1; |
|
int iscol; |
GEN piv,m,u; |
GEN piv,m,u; |
|
|
if (typ(a)!=t_MAT) err(mattype1,"gauss"); |
if (typ(a)!=t_MAT) err(mattype1,"gauss"); |
Line 1215 FpM_gauss(GEN a, GEN b, GEN p) |
|
Line 1412 FpM_gauss(GEN a, GEN b, GEN p) |
|
} |
} |
li = lg(a[1])-1; |
li = lg(a[1])-1; |
if (li != aco && (li < aco || b)) err(mattype1,"gauss"); |
if (li != aco && (li < aco || b)) err(mattype1,"gauss"); |
b = check_b(b,li); av = avma; |
b = check_b(b,li, &iscol); av = avma; |
if (OK_ULONG(p)) |
if (OK_ULONG(p)) |
{ |
{ |
ulong pp=p[2]; |
ulong pp=p[2]; |
a = u_Fp_FpM(a, pp); |
a = u_Fp_FpM(a, pp); |
b = u_Fp_FpM(b, pp); |
b = u_Fp_FpM(b, pp); |
u = u_FpM_gauss(a,b, pp); |
u = u_FpM_gauss_sp(a,b, pp); |
return gerepileupto(av, small_to_mat(u)); |
u = iscol? small_to_col((GEN)u[1]): small_to_mat(u); |
|
return gerepileupto(av, u); |
} |
} |
lim = stack_lim(av,1); |
lim = stack_lim(av,1); |
a = dummycopy(a); |
a = dummycopy(a); |
bco = lg(b)-1; |
bco = lg(b)-1; |
iscol = (typ(b)==t_COL); |
|
piv = NULL; /* gcc -Wall */ |
piv = NULL; /* gcc -Wall */ |
for (i=1; i<=aco; i++) |
for (i=1; i<=aco; i++) |
{ |
{ |
Line 1240 FpM_gauss(GEN a, GEN b, GEN p) |
|
Line 1437 FpM_gauss(GEN a, GEN b, GEN p) |
|
{ |
{ |
coeff(a,k,i) = lresii(gcoeff(a,k,i), p); |
coeff(a,k,i) = lresii(gcoeff(a,k,i), p); |
if (signe(coeff(a,k,i))) break; |
if (signe(coeff(a,k,i))) break; |
} |
} |
if (k>li) return NULL; |
if (k>li) return NULL; |
} |
} |
|
|
Line 1248 FpM_gauss(GEN a, GEN b, GEN p) |
|
Line 1445 FpM_gauss(GEN a, GEN b, GEN p) |
|
if (k != i) |
if (k != i) |
{ |
{ |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
if (iscol) { swap(b[i],b[k]); } |
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
else |
|
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
|
piv = gcoeff(a,i,i); |
piv = gcoeff(a,i,i); |
} |
} |
if (i == aco) break; |
if (i == aco) break; |
Line 1264 FpM_gauss(GEN a, GEN b, GEN p) |
|
Line 1459 FpM_gauss(GEN a, GEN b, GEN p) |
|
m = mulii(m, mpinvmod(piv,p)); |
m = mulii(m, mpinvmod(piv,p)); |
m = resii(negi(m), p); |
m = resii(negi(m), p); |
for (j=i+1; j<=aco; j++) _Fp_addmul((GEN)a[j],k,i,m, p); |
for (j=i+1; j<=aco; j++) _Fp_addmul((GEN)a[j],k,i,m, p); |
if (iscol) _Fp_addmul(b,k,i,m, p); |
for (j=1 ; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p); |
else |
|
for (j=1; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p); |
|
} |
} |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
Line 1278 FpM_gauss(GEN a, GEN b, GEN p) |
|
Line 1471 FpM_gauss(GEN a, GEN b, GEN p) |
|
} |
} |
|
|
if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); |
if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); |
if (iscol) u = Fp_gauss_get_col(a,b,piv,aco,p); |
u = cgetg(bco+1,t_MAT); |
else |
for (j=1; j<=bco; j++) u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); |
|
return gerepilecopy(av, iscol? (GEN)u[1]: u); |
|
} |
|
GEN |
|
FqM_gauss(GEN a, GEN b, GEN T, GEN p) |
|
{ |
|
gpmem_t av,lim; |
|
long i,j,k,li,bco, aco = lg(a)-1; |
|
int iscol; |
|
GEN piv,m,u; |
|
|
|
if (!T) return FpM_gauss(a,b,p); |
|
|
|
if (typ(a)!=t_MAT) err(mattype1,"gauss"); |
|
if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss"); |
|
if (!aco) |
{ |
{ |
long av1 = avma; |
if (b && lg(b)!=1) err(consister,"gauss"); |
lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT); |
if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1); |
for (j=2; j<=bco; j++) u[j] = zero; /* dummy */ |
return cgetg(1,t_MAT); |
for (j=1; j<=bco; j++) |
} |
|
li = lg(a[1])-1; |
|
if (li != aco && (li < aco || b)) err(mattype1,"gauss"); |
|
b = check_b(b,li,&iscol); av = avma; |
|
lim = stack_lim(av,1); |
|
a = dummycopy(a); |
|
bco = lg(b)-1; |
|
piv = NULL; /* gcc -Wall */ |
|
for (i=1; i<=aco; i++) |
|
{ |
|
/* k is the line where we find the pivot */ |
|
coeff(a,i,i) = (long)FpX_res(gcoeff(a,i,i), T,p); |
|
piv = gcoeff(a,i,i); k = i; |
|
if (!signe(piv)) |
{ |
{ |
u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); |
for (k++; k <= li; k++) |
if (low_stack(lim, stack_lim(av1,1))) |
|
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j); |
coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p); |
u = gerepilecopy(av1, u); |
if (signe(coeff(a,k,i))) break; |
} |
} |
|
if (k>li) return NULL; |
} |
} |
|
|
|
/* if (k!=i), exchange the lines s.t. k = i */ |
|
if (k != i) |
|
{ |
|
for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); |
|
for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); |
|
piv = gcoeff(a,i,i); |
|
} |
|
if (i == aco) break; |
|
|
|
for (k=i+1; k<=li; k++) |
|
{ |
|
coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p); |
|
m = gcoeff(a,k,i); coeff(a,k,i) = zero; |
|
if (signe(m)) |
|
{ |
|
m = FpXQ_mul(m, FpXQ_inv(piv,T,p),T,p); |
|
m = resii(negi(m), p); |
|
for (j=i+1; j<=aco; j++) _Fq_addmul((GEN)a[j],k,i,m, T,p); |
|
for (j=1; j<=bco; j++) _Fq_addmul((GEN)b[j],k,i,m, T,p); |
|
} |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b; |
|
if(DEBUGMEM>1) err(warnmem,"FpM_gauss. i=%ld",i); |
|
gerepilemany(av,gptr,2); |
|
} |
} |
} |
return gerepilecopy(av, u); |
|
|
if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); |
|
u = cgetg(bco+1,t_MAT); |
|
for (j=1; j<=bco; j++) u[j] = (long)Fq_gauss_get_col(a,(GEN)b[j],piv,aco,T,p); |
|
return gerepilecopy(av, iscol? (GEN)u[1]: u); |
} |
} |
|
|
|
|
GEN |
GEN |
FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); } |
FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); } |
|
|
static GEN |
|
u_idmat(long n) |
|
{ |
|
GEN y = cgetg(n+1,t_MAT); |
|
long i,j; |
|
if (n < 0) err(talker,"negative size in u_idmat"); |
|
for (i=1; i<=n; i++) |
|
{ |
|
y[i]=lgetg(n+1,t_VECSMALL); |
|
for (j=1; j<=n; j++) coeff(y,j,i) = (i==j)? 1: 0; |
|
} |
|
return y; |
|
} |
|
|
|
/* set y = x * y */ |
/* set y = x * y */ |
#define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG) |
#define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG) |
static GEN |
static GEN |
Line 1332 u_FpM_Fp_mul_ip(GEN y, long x, long p) |
|
Line 1572 u_FpM_Fp_mul_ip(GEN y, long x, long p) |
|
} |
} |
|
|
/* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */ |
/* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */ |
GEN |
GEN |
ZM_inv(GEN M, GEN dM) |
ZM_inv(GEN M, GEN dM) |
{ |
{ |
GEN ID,Hp,q,H; |
gpmem_t av2, av = avma, lim = stack_lim(av,1); |
ulong p,dMp, av2, av = avma, lim = stack_lim(av,1); |
GEN Hp,q,H; |
|
ulong p,dMp; |
byteptr d = diffptr; |
byteptr d = diffptr; |
long lM = lg(M), stable = 0; |
long lM = lg(M), stable = 0; |
|
|
Line 1346 ZM_inv(GEN M, GEN dM) |
|
Line 1587 ZM_inv(GEN M, GEN dM) |
|
av2 = avma; |
av2 = avma; |
H = NULL; |
H = NULL; |
d += 3000; /* 27449 = prime(3000) */ |
d += 3000; /* 27449 = prime(3000) */ |
for(p = 27449; ; p+= *d++) |
for(p = 27449; ; ) |
{ |
{ |
if (!*d) err(primer1); |
|
dMp = umodiu(dM,p); |
dMp = umodiu(dM,p); |
if (!dMp) continue; |
if (!dMp) goto repeat; |
ID = u_idmat(lM-1); |
Hp = u_FpM_inv_sp(u_Fp_FpM(M,p), p); |
Hp = u_FpM_gauss(u_Fp_FpM(M,p), ID, p); |
|
if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p); |
if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p); |
|
|
if (!H) |
if (!H) |
{ |
{ |
H = ZM_init_CRT(Hp, p); |
H = ZM_init_CRT(Hp, p); |
q = utoi(p); |
q = utoi(p); |
Line 1368 ZM_inv(GEN M, GEN dM) |
|
Line 1607 ZM_inv(GEN M, GEN dM) |
|
} |
} |
if (DEBUGLEVEL>5) msgtimer("inverse mod %ld (stable=%ld)", p,stable); |
if (DEBUGLEVEL>5) msgtimer("inverse mod %ld (stable=%ld)", p,stable); |
if (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */ |
if (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */ |
|
|
if (low_stack(lim, stack_lim(av,2))) |
if (low_stack(lim, stack_lim(av,2))) |
{ |
{ |
GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q; |
GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q; |
if (DEBUGMEM>1) err(warnmem,"ZM_inv"); |
if (DEBUGMEM>1) err(warnmem,"ZM_inv"); |
gerepilemany(av2,gptr, 2); |
gerepilemany(av2,gptr, 2); |
} |
} |
|
repeat: |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
if (DEBUGLEVEL>5) msgtimer("ZM_inv done"); |
if (DEBUGLEVEL>5) msgtimer("ZM_inv done"); |
return gerepilecopy(av, H); |
return gerepilecopy(av, H); |
} |
} |
|
|
/* same as above, M rational */ |
/* same as above, M rational */ |
GEN |
GEN |
QM_inv(GEN M, GEN dM) |
QM_inv(GEN M, GEN dM) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN cM = content(M); |
GEN cM, pM = Q_primitive_part(M, &cM); |
if (is_pm1(cM)) { avma = av; return ZM_inv(M,dM); } |
if (!cM) return ZM_inv(pM,dM); |
return gerepileupto(av, ZM_inv(gdiv(M,cM), gdiv(dM,cM))); |
return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM))); |
} |
} |
|
|
/* x a matrix with integer coefficients. Return a multiple of the determinant |
/* x a matrix with integer coefficients. Return a multiple of the determinant |
Line 1397 QM_inv(GEN M, GEN dM) |
|
Line 1637 QM_inv(GEN M, GEN dM) |
|
GEN |
GEN |
detint(GEN x) |
detint(GEN x) |
{ |
{ |
|
gpmem_t av = avma, av1, lim; |
GEN pass,c,v,det1,piv,pivprec,vi,p1; |
GEN pass,c,v,det1,piv,pivprec,vi,p1; |
long i,j,k,rg,n,m,m1,av=avma,av1,lim; |
long i,j,k,rg,n,m,m1; |
|
|
if (typ(x)!=t_MAT) err(typeer,"detint"); |
if (typ(x)!=t_MAT) err(typeer,"detint"); |
n=lg(x)-1; if (!n) return gun; |
n=lg(x)-1; if (!n) return gun; |
|
|
} |
} |
|
|
static void |
static void |
gerepile_gauss_ker(GEN x, long m, long n, long k, long t, ulong av) |
gerepile_gauss_ker(GEN x, long k, long t, gpmem_t av) |
{ |
{ |
ulong tetpil = avma,A; |
gpmem_t tetpil = avma, A; |
long dec,u,i; |
long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; |
|
size_t dec; |
|
|
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k)); |
for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k)); |
Line 1475 gerepile_gauss_ker(GEN x, long m, long n, long k, long |
|
Line 1717 gerepile_gauss_ker(GEN x, long m, long n, long k, long |
|
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,k); |
A=(gpmem_t)coeff(x,u,k); |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
} |
} |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
for (u=1; u<=m; u++) |
for (u=1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,i); |
A=(gpmem_t)coeff(x,u,i); |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
} |
} |
} |
} |
|
|
static void |
static void |
gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, long k, long t, ulong av) |
gerepile_gauss_FpM_ker(GEN x, GEN p, long k, long t, gpmem_t av) |
{ |
{ |
ulong tetpil = avma,A; |
gpmem_t tetpil = avma, A; |
long dec,u,i; |
long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; |
|
size_t dec; |
|
|
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
Line 1502 gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l |
|
Line 1745 gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l |
|
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,k); |
A=(gpmem_t)coeff(x,u,k); |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
} |
} |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
for (u=1; u<=m; u++) |
for (u=1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,i); |
A=(gpmem_t)coeff(x,u,i); |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
} |
} |
} |
} |
Line 1516 gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l |
|
Line 1759 gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l |
|
/* special gerepile for huge matrices */ |
/* special gerepile for huge matrices */ |
|
|
static void |
static void |
gerepile_gauss(GEN x,long m,long n,long k,long t,ulong av, long j, GEN c) |
gerepile_gauss(GEN x,long k,long t,gpmem_t av, long j, GEN c) |
{ |
{ |
ulong tetpil = avma,A; |
gpmem_t tetpil = avma, A; |
long dec,u,i; |
long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; |
|
size_t dec; |
|
|
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n); |
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n); |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
Line 1532 gerepile_gauss(GEN x,long m,long n,long k,long t,ulong |
|
Line 1776 gerepile_gauss(GEN x,long m,long n,long k,long t,ulong |
|
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
if (u==j || !c[u]) |
if (u==j || !c[u]) |
{ |
{ |
A=coeff(x,u,k); |
A=(gpmem_t)coeff(x,u,k); |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
} |
} |
for (u=1; u<=m; u++) |
for (u=1; u<=m; u++) |
if (u==j || !c[u]) |
if (u==j || !c[u]) |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
{ |
{ |
A=coeff(x,u,i); |
A=(gpmem_t)coeff(x,u,i); |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
} |
} |
} |
} |
Line 1550 gerepile_gauss(GEN x,long m,long n,long k,long t,ulong |
|
Line 1794 gerepile_gauss(GEN x,long m,long n,long k,long t,ulong |
|
/* return n - rk(x) linearly independant vectors */ |
/* return n - rk(x) linearly independant vectors */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
|
static long |
|
gauss_get_pivot_NZ(GEN x, GEN x0/*unused*/, GEN c, long i0) |
|
{ |
|
long i,lx = lg(x); |
|
(void)x0; |
|
if (c) |
|
for (i=i0; i<lx; i++) |
|
{ |
|
if (c[i]) continue; /* already a pivot in line i */ |
|
if (!gcmp0((GEN)x[i])) break; |
|
} |
|
else |
|
for (i=i0; i<lx; i++) |
|
if (!gcmp0((GEN)x[i])) break; |
|
return i; |
|
|
|
} |
|
|
|
/* x ~ 0 compared to reference y */ |
|
int |
|
approx_0(GEN x, GEN y) |
|
{ |
|
long tx = typ(x); |
|
if (tx == t_COMPLEX) |
|
return approx_0((GEN)x[1], y) && approx_0((GEN)x[2], y); |
|
return gcmp0(x) || |
|
(tx == t_REAL && gexpo(y) - gexpo(x) > bit_accuracy(lg(x))); |
|
} |
|
|
|
static long |
|
gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0) |
|
{ |
|
long i,e, k = i0, ex = - (long)HIGHEXPOBIT, lx = lg(x); |
|
GEN p,r; |
|
if (c) |
|
for (i=i0; i<lx; i++) |
|
{ |
|
if (c[i]) continue; |
|
e = gexpo((GEN)x[i]); |
|
if (e > ex) { ex=e; k=i; } |
|
} |
|
else |
|
for (i=i0; i<lx; i++) |
|
{ |
|
e = gexpo((GEN)x[i]); |
|
if (e > ex) { ex=e; k=i; } |
|
} |
|
p = (GEN)x[k]; |
|
r = (GEN)x0[k]; if (isexactzero(r)) r = x0; |
|
return approx_0(p, r)? i: k; |
|
} |
|
|
/* x has INTEGER coefficients */ |
/* x has INTEGER coefficients */ |
GEN |
GEN |
keri(GEN x) |
keri(GEN x) |
{ |
{ |
GEN c,d,y,p,pp; |
gpmem_t av, av0, tetpil, lim; |
long i,j,k,r,t,n,m,av,av0,tetpil,lim; |
GEN c,l,y,p,pp; |
|
long i,j,k,r,t,n,m; |
|
|
if (typ(x)!=t_MAT) err(typeer,"keri"); |
if (typ(x)!=t_MAT) err(typeer,"keri"); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
|
|
av0=avma; m=lg(x[1])-1; r=0; |
av0=avma; m=lg(x[1])-1; r=0; |
pp=cgetg(n+1,t_COL); |
pp=cgetg(n+1,t_COL); |
x=dummycopy(x); p=gun; |
x=dummycopy(x); p=gun; |
c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; |
c=cgetg(m+1, t_VECSMALL); for (k=1; k<=m; k++) c[k]=0; |
d=new_chunk(n+1); av=avma; lim=stack_lim(av,1); |
l=cgetg(n+1, t_VECSMALL); |
|
av = avma; lim = stack_lim(av,1); |
for (k=1; k<=n; k++) |
for (k=1; k<=n; k++) |
{ |
{ |
j=1; |
j = 1; |
while (j<=m && (c[j] || !signe(gcoeff(x,j,k))) ) j++; |
while ( j <= m && (c[j] || !signe(gcoeff(x,j,k))) ) j++; |
if (j>m) |
if (j > m) |
{ |
{ |
r++; d[k]=0; |
r++; l[k] = 0; |
for(j=1; j<k; j++) |
for(j=1; j<k; j++) |
if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k)); |
if (l[j]) coeff(x,l[j],k) = lclone(gcoeff(x,l[j],k)); |
pp[k]=lclone(p); |
pp[k]=lclone(p); |
} |
} |
else |
else |
{ |
{ |
GEN p0 = p; |
GEN p0 = p; |
long av1; |
c[j]=k; l[k]=j; p = gcoeff(x,j,k); |
|
|
c[j]=k; d[k]=j; p = gcoeff(x,j,k); |
|
|
|
for (t=1; t<=m; t++) |
for (t=1; t<=m; t++) |
if (t!=j) |
if (t!=j) |
{ |
{ |
GEN q=gcoeff(x,t,k), p1,p2; |
GEN q=gcoeff(x,t,k), p1; |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
{ |
{ |
av1=avma; (void)new_chunk(lgefint(p0)); |
gpmem_t av1 = avma; |
p1=mulii(q,gcoeff(x,j,i)); |
p1 = subii(mulii(p,gcoeff(x,t,i)), mulii(q,gcoeff(x,j,i))); |
p2=mulii(p,gcoeff(x,t,i)); |
p1 = gerepileuptoint(av1, diviiexact(p1,p0)); |
p1=subii(p2,p1); avma=av1; |
coeff(x,t,i) = (long)p1; |
coeff(x,t,i) = ldivii(p1,p0); |
|
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
p1 = gclone(p); |
p1 = gclone(p); |
gerepile_gauss_ker(x,m,n,k,t,av); |
gerepile_gauss_ker(x,k,t,av); |
p = gcopy(p1); gunclone(p1); |
p = gcopy(p1); gunclone(p1); |
} |
} |
} |
} |
|
|
for (j=k=1; j<=r; j++,k++) |
for (j=k=1; j<=r; j++,k++) |
{ |
{ |
p = cgetg(n+1, t_COL); |
p = cgetg(n+1, t_COL); |
y[j]=(long)p; while (d[k]) k++; |
y[j]=(long)p; while (l[k]) k++; |
for (i=1; i<k; i++) |
for (i=1; i<k; i++) |
if (d[i]) |
if (l[i]) |
{ |
{ |
c=gcoeff(x,d[i],k); |
c=gcoeff(x,l[i],k); |
p[i] = (long) forcecopy(c); gunclone(c); |
p[i] = (long) forcecopy(c); gunclone(c); |
} |
} |
else |
else |
|
|
} |
} |
|
|
GEN |
GEN |
deplin(GEN x) |
deplin(GEN x0) |
{ |
{ |
long i,j,k,t,nc,nl, av=avma; |
gpmem_t av = avma; |
GEN y,q,c,l,d; |
long i,j,k,t,nc,nl; |
|
GEN x,y,piv,q,c,l,*d,*ck,*cj; |
|
|
if (typ(x) != t_MAT) err(typeer,"deplin"); |
t = typ(x0); |
nc=lg(x)-1; if (!nc) err(talker,"empty matrix in deplin"); |
if (t == t_MAT) x = dummycopy(x0); |
nl=lg(x[1])-1; |
else |
c=new_chunk(nl+1); |
|
l=new_chunk(nc+1); |
|
d=cgetg(nl+1,t_VEC); |
|
for (i=1; i<=nl; i++) { d[i]=un; c[i]=0; } |
|
k=1; t=1; |
|
while (t<=nl && k<=nc) |
|
{ |
{ |
|
if (t != t_VEC) err(typeer,"deplin"); |
|
x = gtomat(x0); |
|
} |
|
nc = lg(x)-1; if (!nc) err(talker,"empty matrix in deplin"); |
|
nl = lg(x[1])-1; |
|
d = (GEN*)cgetg(nl+1,t_VEC); /* pivot list */ |
|
c = cgetg(nl+1, t_VECSMALL); |
|
l = cgetg(nc+1, t_VECSMALL); /* not initialized */ |
|
for (i=1; i<=nl; i++) { d[i] = gun; c[i] = 0; } |
|
ck = NULL; /* gcc -Wall */ |
|
for (k=1; k<=nc; k++) |
|
{ |
|
ck = (GEN*)x[k]; |
for (j=1; j<k; j++) |
for (j=1; j<k; j++) |
for (i=1; i<=nl; i++) |
|
if (i!=l[j]) |
|
coeff(x,i,k)=lsub(gmul((GEN) d[j],gcoeff(x,i,k)), |
|
gmul(gcoeff(x,i,j),gcoeff(x,l[j],k))); |
|
t=1; |
|
while ( t<=nl && (c[t] || gcmp0(gcoeff(x,t,k))) ) t++; |
|
if (t<=nl) |
|
{ |
{ |
d[k]=coeff(x,t,k); |
cj = (GEN*)x[j]; piv = d[j]; q = gneg(ck[l[j]]); |
c[t]=k; l[k++]=t; |
for (i=1; i<=nl; i++) |
|
if (i != l[j]) ck[i] = gadd(gmul(piv, ck[i]), gmul(q, cj[i])); |
} |
} |
|
|
|
i = gauss_get_pivot_NZ((GEN)ck, NULL, c, 1); |
|
if (i > nl) break; |
|
|
|
d[k] = ck[i]; |
|
c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */ |
} |
} |
if (k>nc) |
if (k > nc) { avma = av; return zerocol(nc); } |
|
if (k == 1) { avma = av; return gscalcol_i(gun,nc); } |
|
y = cgetg(nc+1,t_COL); |
|
y[1] = (long)ck[ l[1] ]; |
|
for (q=d[1],j=2; j<k; j++) |
{ |
{ |
avma=av; y=cgetg(nc+1,t_COL); |
y[j] = lmul(ck[ l[j] ], q); |
for (j=1; j<=nc; j++) y[j]=zero; |
q = gmul(q, d[j]); |
return y; |
|
} |
} |
y=cgetg(nc+1,t_COL); |
y[j] = lneg(q); |
y[1]=(k>1)? coeff(x,l[1],k): un; |
for (j++; j<=nc; j++) y[j] = zero; |
for (q=gun,j=2; j<k; j++) |
return gerepileupto(av, gdiv(y,content(y))); |
{ |
|
q=gmul(q,(GEN) d[j-1]); |
|
y[j]=lmul(gcoeff(x,l[j],k),q); |
|
} |
|
if (k>1) y[k]=lneg(gmul(q,(GEN) d[k-1])); |
|
for (j=k+1; j<=nc; j++) y[j]=zero; |
|
d=content(y); t=avma; |
|
return gerepile(av,t,gdiv(y,d)); |
|
} |
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
|
|
/* (kernel, image, complementary image, rank) */ |
/* (kernel, image, complementary image, rank) */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
static long gauss_ex; |
|
static int (*gauss_is_zero)(GEN); |
|
|
|
static int |
|
real0(GEN x) |
|
{ |
|
return gcmp0(x) || (gexpo(x) < gauss_ex); |
|
} |
|
|
|
static void |
|
gauss_get_prec(GEN x, long prec) |
|
{ |
|
long pr = matprec(x); |
|
|
|
if (!pr) { gauss_is_zero = &gcmp0; return; } |
|
if (pr > prec) prec = pr; |
|
|
|
gauss_ex = - (long)(0.85 * bit_accuracy(prec)); |
|
gauss_is_zero = &real0; |
|
} |
|
|
|
static long |
|
gauss_get_pivot_NZ(GEN x, GEN x0/* unused */, GEN c, long i0) |
|
{ |
|
long i,lx = lg(x); |
|
if (c) |
|
for (i=i0; i<lx; i++) |
|
{ |
|
if (c[i]) continue; |
|
if (!gcmp0((GEN)x[i])) break; |
|
} |
|
else |
|
for (i=i0; i<lx; i++) |
|
if (!gcmp0((GEN)x[i])) break; |
|
return i; |
|
|
|
} |
|
|
|
/* ~ 0 compared to reference y */ |
|
static int |
|
approx_0(GEN x, GEN y) |
|
{ |
|
long tx = typ(x); |
|
if (tx == t_COMPLEX) |
|
return approx_0((GEN)x[1], y) && approx_0((GEN)x[2], y); |
|
return gcmp0(x) || |
|
(tx == t_REAL && gexpo(y) - gexpo(x) > bit_accuracy(lg(x))); |
|
} |
|
|
|
static long |
|
gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0) |
|
{ |
|
long i,e, k = i0, ex = -HIGHEXPOBIT, lx = lg(x); |
|
GEN p,r; |
|
if (c) |
|
for (i=i0; i<lx; i++) |
|
{ |
|
if (c[i]) continue; |
|
e = gexpo((GEN)x[i]); |
|
if (e > ex) { ex=e; k=i; } |
|
} |
|
else |
|
for (i=i0; i<lx; i++) |
|
{ |
|
e = gexpo((GEN)x[i]); |
|
if (e > ex) { ex=e; k=i; } |
|
} |
|
p = (GEN)x[k]; |
|
r = (GEN)x0[k]; if (isexactzero(r)) r = x0; |
|
return approx_0(p, r)? i: k; |
|
} |
|
|
|
/* return the transform of x under a standard Gauss pivot. r = dim ker(x). |
/* return the transform of x under a standard Gauss pivot. r = dim ker(x). |
* d[k] contains the index of the first non-zero pivot in column k |
* d[k] contains the index of the first non-zero pivot in column k |
* If a != NULL, use x - a Id instead (for eigen) |
* If a != NULL, use x - a Id instead (for eigen) |
|
|
gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) |
gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) |
{ |
{ |
GEN x,c,d,p,mun; |
GEN x,c,d,p,mun; |
long i,j,k,r,t,n,m,av,lim; |
gpmem_t av, lim; |
|
long i,j,k,r,t,n,m; |
long (*get_pivot)(GEN,GEN,GEN,long); |
long (*get_pivot)(GEN,GEN,GEN,long); |
|
|
if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); |
if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); |
Line 1801 gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) |
|
Line 2027 gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) |
|
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i))); |
coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i))); |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
gerepile_gauss_ker(x,m,n,k,t,av); |
gerepile_gauss_ker(x,k,t,av); |
} |
} |
} |
} |
} |
} |
|
|
gauss_pivot(GEN x0, GEN *dd, long *rr) |
gauss_pivot(GEN x0, GEN *dd, long *rr) |
{ |
{ |
GEN x,c,d,d0,mun,p; |
GEN x,c,d,d0,mun,p; |
long i,j,k,r,t,n,m,av,lim; |
long i, j, k, r, t, n, m; |
|
gpmem_t av, lim; |
long (*get_pivot)(GEN,GEN,GEN,long); |
long (*get_pivot)(GEN,GEN,GEN,long); |
|
|
if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); |
if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); |
Line 1856 gauss_pivot(GEN x0, GEN *dd, long *rr) |
|
Line 2083 gauss_pivot(GEN x0, GEN *dd, long *rr) |
|
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i))); |
coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i))); |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
gerepile_gauss(x,m,n,k,t,av,j,c); |
gerepile_gauss(x,k,t,av,j,c); |
} |
} |
for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ |
for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ |
} |
} |
Line 1868 gauss_pivot(GEN x0, GEN *dd, long *rr) |
|
Line 2095 gauss_pivot(GEN x0, GEN *dd, long *rr) |
|
static GEN |
static GEN |
ker0(GEN x, GEN a) |
ker0(GEN x, GEN a) |
{ |
{ |
|
gpmem_t av = avma, tetpil; |
GEN d,y; |
GEN d,y; |
long i,j,k,r,n, av = avma, tetpil; |
long i,j,k,r,n; |
|
|
x = gauss_pivot_ker(x,a,&d,&r); |
x = gauss_pivot_ker(x,a,&d,&r); |
if (!r) { avma=av; return cgetg(1,t_MAT); } |
if (!r) { avma=av; return cgetg(1,t_MAT); } |
Line 1907 matker0(GEN x,long flag) |
|
Line 2135 matker0(GEN x,long flag) |
|
GEN |
GEN |
image(GEN x) |
image(GEN x) |
{ |
{ |
|
gpmem_t av = avma; |
GEN d,y; |
GEN d,y; |
long j,k,r, av = avma; |
long j,k,r; |
|
|
gauss_pivot(x,&d,&r); |
gauss_pivot(x,&d,&r); |
|
|
|
|
GEN |
GEN |
imagecompl(GEN x) |
imagecompl(GEN x) |
{ |
{ |
|
gpmem_t av = avma; |
GEN d,y; |
GEN d,y; |
long j,i,r,av = avma; |
long j,i,r; |
|
|
gauss_pivot(x,&d,&r); |
gauss_pivot(x,&d,&r); |
avma=av; y=cgetg(r+1,t_VEC); |
avma=av; y=cgetg(r+1,t_VEC); |
Line 1940 imagecompl(GEN x) |
|
Line 2170 imagecompl(GEN x) |
|
GEN |
GEN |
imagecomplspec(GEN x, long *nlze) |
imagecomplspec(GEN x, long *nlze) |
{ |
{ |
|
gpmem_t av = avma; |
GEN d,y; |
GEN d,y; |
long i,j,k,l,r,av = avma; |
long i,j,k,l,r; |
|
|
x = gtrans(x); l = lg(x); |
x = gtrans_i(x); l = lg(x); |
gauss_pivot(x,&d,&r); |
gauss_pivot(x,&d,&r); |
avma=av; y = cgetg(l,t_VECSMALL); |
avma=av; y = cgetg(l,t_VECSMALL); |
for (i=j=1, k=r+1; i<l; i++) |
for (i=j=1, k=r+1; i<l; i++) |
Line 1955 imagecomplspec(GEN x, long *nlze) |
|
Line 2186 imagecomplspec(GEN x, long *nlze) |
|
static GEN |
static GEN |
sinverseimage(GEN mat, GEN y) |
sinverseimage(GEN mat, GEN y) |
{ |
{ |
long av=avma,tetpil,i, nbcol = lg(mat); |
gpmem_t av=avma,tetpil; |
|
long i, nbcol = lg(mat); |
GEN col,p1 = cgetg(nbcol+1,t_MAT); |
GEN col,p1 = cgetg(nbcol+1,t_MAT); |
|
|
if (nbcol==1) return NULL; |
if (nbcol==1) return NULL; |
Line 1977 sinverseimage(GEN mat, GEN y) |
|
Line 2209 sinverseimage(GEN mat, GEN y) |
|
GEN |
GEN |
inverseimage(GEN m,GEN v) |
inverseimage(GEN m,GEN v) |
{ |
{ |
long av=avma,j,lv,tv=typ(v); |
gpmem_t av=avma; |
|
long j,lv,tv=typ(v); |
GEN y,p1; |
GEN y,p1; |
|
|
if (typ(m)!=t_MAT) err(typeer,"inverseimage"); |
if (typ(m)!=t_MAT) err(typeer,"inverseimage"); |
Line 1999 inverseimage(GEN m,GEN v) |
|
Line 2232 inverseimage(GEN m,GEN v) |
|
return y; |
return y; |
} |
} |
|
|
/* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix |
/* i-th vector in the standard basis */ |
* whose first k columns are given by x. If rank(x)<k, the result may be wrong |
|
*/ |
|
GEN |
GEN |
suppl_intern(GEN x, GEN myid) |
_ei(long n, long i) |
{ |
{ |
long av = avma, lx = lg(x), n,i,j; |
GEN e = zerocol(n); e[i] = un; return e; |
GEN y,p1; |
} |
stackzone *zone; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"suppl"); |
/* NB: d freed */ |
if (lx==1) err(talker,"empty matrix in suppl"); |
static GEN |
n=lg(x[1]); if (lx>n) err(suppler2); |
get_suppl(GEN x, GEN d, long r) |
if (lx == n) return gcopy(x); |
{ |
|
gpmem_t av; |
|
GEN y,c; |
|
long j,k,rx,n; |
|
|
zone = switch_stack(NULL, n*n); |
rx = lg(x)-1; |
switch_stack(zone,1); |
if (!rx) err(talker,"empty matrix in suppl"); |
y = myid? dummycopy(myid): idmat(n-1); |
n = lg(x[1])-1; |
switch_stack(zone,0); |
if (rx == n && r == 0) { free(d); return gcopy(x); } |
gauss_get_prec(x,0); |
y = cgetg(n+1, t_MAT); |
for (i=1; i<lx; i++) |
av = avma; |
{ |
c = vecsmall_const(n,0); |
p1 = gauss(y,(GEN)x[i]); j=i; |
k = 1; |
while (j<n && gauss_is_zero((GEN)p1[j])) j++; |
/* c = lines containing pivots (could get it from gauss_pivot, but cheap) |
if (j>=n) err(suppler2); |
* In theory r = 0 and d[j] > 0 for all j, but why take chances? */ |
p1=(GEN)y[i]; y[i]=x[i]; if (i!=j) y[j]=(long)p1; |
for (j=1; j<=rx; j++) |
} |
if (d[j]) |
avma = av; y = gcopy(y); |
{ |
free(zone); return y; |
c[ d[j] ] = 1; |
|
y[k++] = x[j]; |
|
} |
|
for (j=1; j<=n; j++) |
|
if (!c[j]) y[k++] = j; |
|
avma = av; |
|
|
|
rx -= r; |
|
for (j=1; j<=rx; j++) |
|
y[j] = lcopy((GEN)y[j]); |
|
for ( ; j<=n; j++) |
|
y[j] = (long)_ei(n, y[j]); |
|
free(d); return y; |
} |
} |
|
|
|
/* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix |
|
* whose first k columns are given by x. If rank(x) < k, undefined result. */ |
GEN |
GEN |
suppl(GEN x) |
suppl(GEN x) |
{ |
{ |
return suppl_intern(x,NULL); |
gpmem_t av = avma; |
|
GEN d; |
|
long r; |
|
|
|
gauss_pivot(x,&d,&r); |
|
avma = av; return get_suppl(x,d,r); |
} |
} |
|
|
|
static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr); |
|
static void FqM_gauss_pivot(GEN x, GEN T, GEN p, GEN *dd, long *rr); |
|
|
GEN |
GEN |
|
FpM_suppl(GEN x, GEN p) |
|
{ |
|
gpmem_t av = avma; |
|
GEN d; |
|
long r; |
|
|
|
FpM_gauss_pivot(x,p, &d,&r); |
|
avma = av; return get_suppl(x,d,r); |
|
} |
|
GEN |
|
FqM_suppl(GEN x, GEN T, GEN p) |
|
{ |
|
gpmem_t av = avma; |
|
GEN d; |
|
long r; |
|
|
|
if (!T) return FpM_suppl(x,p); |
|
FqM_gauss_pivot(x,T,p, &d,&r); |
|
avma = av; return get_suppl(x,d,r); |
|
} |
|
|
|
GEN |
image2(GEN x) |
image2(GEN x) |
{ |
{ |
long av=avma,tetpil,k,n,i; |
gpmem_t av=avma,tetpil; |
|
long k,n,i; |
GEN p1,p2; |
GEN p1,p2; |
|
|
if (typ(x)!=t_MAT) err(typeer,"image2"); |
if (typ(x)!=t_MAT) err(typeer,"image2"); |
Line 2068 matimage0(GEN x,long flag) |
|
Line 2346 matimage0(GEN x,long flag) |
|
long |
long |
rank(GEN x) |
rank(GEN x) |
{ |
{ |
long av = avma, r; |
gpmem_t av = avma; |
|
long r; |
GEN d; |
GEN d; |
|
|
gauss_pivot(x,&d,&r); |
gauss_pivot(x,&d,&r); |
|
|
return lg(x)-1 - r; |
return lg(x)-1 - r; |
} |
} |
|
|
static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr); |
|
|
|
/* if p != NULL, assume x integral and compute rank over Fp */ |
/* if p != NULL, assume x integral and compute rank over Fp */ |
static GEN |
static GEN |
indexrank0(GEN x, GEN p, int small) |
indexrank0(GEN x, GEN p, int small) |
{ |
{ |
long av = avma, i,j,n,r; |
gpmem_t av = avma; |
|
long i,j,n,r; |
GEN res,d,p1,p2; |
GEN res,d,p1,p2; |
|
|
/* yield r = dim ker(x) */ |
/* yield r = dim ker(x) */ |
Line 2108 indexrank0(GEN x, GEN p, int small) |
|
Line 2386 indexrank0(GEN x, GEN p, int small) |
|
return res; |
return res; |
} |
} |
|
|
GEN |
GEN |
indexrank(GEN x) { return indexrank0(x,NULL,0); } |
indexrank(GEN x) { return indexrank0(x,NULL,0); } |
|
|
GEN |
GEN |
sindexrank(GEN x) { return indexrank0(x,NULL,1); } |
sindexrank(GEN x) { return indexrank0(x,NULL,1); } |
|
|
GEN |
GEN |
FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1); } |
FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1); } |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
Line 2122 FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1 |
|
Line 2400 FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1 |
|
/* LINEAR ALGEBRA MODULO P */ |
/* LINEAR ALGEBRA MODULO P */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
|
|
/*If p is NULL no reduction is performed.*/ |
/*If p is NULL no reduction is performed.*/ |
GEN |
GEN |
FpM_mul(GEN x, GEN y, GEN p) |
FpM_mul(GEN x, GEN y, GEN p) |
Line 2130 FpM_mul(GEN x, GEN y, GEN p) |
|
Line 2408 FpM_mul(GEN x, GEN y, GEN p) |
|
long i,j,l,lx=lg(x), ly=lg(y); |
long i,j,l,lx=lg(x), ly=lg(y); |
GEN z; |
GEN z; |
if (ly==1) return cgetg(ly,t_MAT); |
if (ly==1) return cgetg(ly,t_MAT); |
if (lx != lg(y[1])) err(operi,"* [mod p]",t_MAT,t_MAT); |
if (lx != lg(y[1])) err(operi,"* [mod p]",x,y); |
z=cgetg(ly,t_MAT); |
z=cgetg(ly,t_MAT); |
if (lx==1) |
if (lx==1) |
{ |
{ |
Line 2142 FpM_mul(GEN x, GEN y, GEN p) |
|
Line 2420 FpM_mul(GEN x, GEN y, GEN p) |
|
{ |
{ |
z[j] = lgetg(l,t_COL); |
z[j] = lgetg(l,t_COL); |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
{ |
ulong av; |
gpmem_t av; |
GEN p1,p2; |
GEN p1,p2; |
int k; |
int k; |
p1=gzero; av=avma; |
p1=gzero; av=avma; |
Line 2162 FpM_mul(GEN x, GEN y, GEN p) |
|
Line 2440 FpM_mul(GEN x, GEN y, GEN p) |
|
GEN |
GEN |
FpM_FpV_mul(GEN x, GEN y, GEN p) |
FpM_FpV_mul(GEN x, GEN y, GEN p) |
{ |
{ |
long i,l,lx=lg(x), ly=lg(y); |
long i,k,l,lx=lg(x), ly=lg(y); |
GEN z; |
GEN z; |
if (lx != ly) err(operi,"* [mod p]",t_MAT,t_COL); |
if (lx != ly) err(operi,"* [mod p]",x,y); |
z=cgetg(ly,t_COL); |
|
if (lx==1) return cgetg(1,t_COL); |
if (lx==1) return cgetg(1,t_COL); |
l=lg(x[1]); |
l = lg(x[1]); |
|
z = cgetg(l,t_COL); |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
{ |
ulong av; |
gpmem_t av = avma; |
GEN p1,p2; |
GEN p1 = gzero; |
int k; |
|
p1=gzero; av=avma; |
|
for (k=1; k<lx; k++) |
for (k=1; k<lx; k++) |
{ |
p1 = addii(p1, mulii(gcoeff(x,i,k),(GEN)y[k])); |
p2=mulii(gcoeff(x,i,k),(GEN)y[k]); |
z[i] = lpileupto(av,p?modii(p1,p):p1); |
p1=addii(p1,p2); |
|
} |
|
z[i]=lpileupto(av,p?modii(p1,p):p1); |
|
} |
} |
return z; |
return z; |
} |
} |
|
|
|
/* in place, destroy x */ |
static GEN |
static GEN |
u_FpM_ker(GEN x, GEN pp, long nontriv) |
u_FpM_ker_sp(GEN x, ulong p, long deplin) |
{ |
{ |
GEN y,c,d; |
GEN y,c,d; |
long a,i,j,k,r,t,n,m,av0,tetpil, p = pp[2], piv; |
long i,j,k,r,t,n; |
|
ulong a, piv, m; |
|
|
n = lg(x)-1; |
n = lg(x)-1; |
m=lg(x[1])-1; r=0; av0 = avma; |
m=lg(x[1])-1; r=0; |
x = dummycopy(x); |
|
for (i=1; i<=n; i++) |
c = new_chunk(m+1); for (k=1; k<=m; k++) c[k] = 0; |
{ |
d = new_chunk(n+1); |
GEN p1 = (GEN)x[i]; |
a = 0; /* for gcc -Wall */ |
for (j=1; j<=m; j++) p1[j] = itos((GEN)p1[j]); |
|
} |
|
c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; |
|
d=new_chunk(n+1); |
|
a=0; /* for gcc -Wall */ |
|
for (k=1; k<=n; k++) |
for (k=1; k<=n; k++) |
{ |
{ |
for (j=1; j<=m; j++) |
for (j=1; j<=m; j++) |
Line 2209 u_FpM_ker(GEN x, GEN pp, long nontriv) |
|
Line 2479 u_FpM_ker(GEN x, GEN pp, long nontriv) |
|
a = coeff(x,j,k) % p; |
a = coeff(x,j,k) % p; |
if (a) break; |
if (a) break; |
} |
} |
if (j>m) |
if (j > m) |
{ |
{ |
if (nontriv) { avma=av0; return NULL; } |
if (deplin) { |
r++; d[k]=0; |
c = cgetg(n+1, t_VECSMALL); |
|
for (i=1; i<k; i++) c[i] = coeff(x,d[i],k) % p; |
|
c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0; |
|
return c; |
|
} |
|
r++; d[k] = 0; |
} |
} |
else |
else |
{ |
{ |
c[j]=k; d[k]=j; |
c[j] = k; d[k] = j; |
piv = - u_Fp_inv(a, p); |
piv = p - u_invmod(a, p); /* -1/a */ |
coeff(x,j,k) = -1; |
coeff(x,j,k) = p-1; |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,j,i) = (piv * coeff(x,j,i)) % p; |
coeff(x,j,i) = (piv * coeff(x,j,i)) % p; |
for (t=1; t<=m; t++) |
for (t=1; t<=m; t++) |
if (t!=j) |
{ |
{ |
if (t==j) continue; |
piv = coeff(x,t,k) % p; |
|
if (piv) |
piv = coeff(x,t,k) % p; |
{ |
if (!piv) continue; |
coeff(x,t,k) = 0; |
|
|
coeff(x,t,k) = 0; |
|
if (piv == 1) |
|
for (i=k+1; i<=n; i++) _u_Fp_add((GEN)x[i],t,j,p); |
|
else |
|
{ |
|
if (u_OK_ULONG(p)) |
|
for (i=k+1; i<=n; i++) _u_Fp_addmul_OK((GEN)x[i],t,j,piv,p); |
|
else |
for (i=k+1; i<=n; i++) _u_Fp_addmul((GEN)x[i],t,j,piv,p); |
for (i=k+1; i<=n; i++) _u_Fp_addmul((GEN)x[i],t,j,piv,p); |
} |
} |
} |
} |
} |
} |
} |
} |
|
if (deplin) return NULL; |
|
|
tetpil=avma; y=cgetg(r+1,t_MAT); |
y = cgetg(r+1, t_MAT); |
for (j=k=1; j<=r; j++,k++) |
for (j=k=1; j<=r; j++,k++) |
{ |
{ |
GEN c = cgetg(n+1,t_COL); |
GEN c = cgetg(n+1, t_VECSMALL); |
|
|
y[j]=(long)c; while (d[k]) k++; |
y[j] = (long)c; while (d[k]) k++; |
for (i=1; i<k; i++) |
for (i=1; i<k; i++) |
if (d[i]) |
if (d[i]) |
{ |
c[i] = coeff(x,d[i],k) % p; |
long a = coeff(x,d[i],k) % p; |
|
if (a < 0) a += p; |
|
c[i] = lstoi(a); |
|
} |
|
else |
else |
c[i] = zero; |
c[i] = 0; |
c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero; |
c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0; |
} |
} |
return gerepile(av0,tetpil,y); |
return y; |
} |
} |
|
|
|
/* assume x has integer entries */ |
static GEN |
static GEN |
FpM_ker_i(GEN x, GEN p, long nontriv) |
FpM_ker_i(GEN x, GEN p, long deplin) |
{ |
{ |
|
gpmem_t av0 = avma, av,lim,tetpil; |
GEN y,c,d,piv,mun; |
GEN y,c,d,piv,mun; |
long i,j,k,r,t,n,m,av0,av,lim,tetpil; |
long i,j,k,r,t,n,m; |
|
|
if (typ(x)!=t_MAT) err(typeer,"FpM_ker"); |
if (typ(x)!=t_MAT) err(typeer,"FpM_ker"); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
if (OK_ULONG(p)) return u_FpM_ker(x, p, nontriv); |
if (OK_ULONG(p)) |
|
{ |
|
ulong pp = p[2]; |
|
y = u_Fp_FpM(x, pp); |
|
y = u_FpM_ker_sp(y, pp, deplin); |
|
if (!y) return y; |
|
y = deplin? small_to_col(y): small_to_mat(y); |
|
return gerepileupto(av0, y); |
|
} |
|
|
m=lg(x[1])-1; r=0; av0 = avma; |
m=lg(x[1])-1; r=0; |
x=dummycopy(x); mun=negi(gun); |
x=dummycopy(x); mun=negi(gun); |
c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; |
c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; |
d=new_chunk(n+1); |
d=new_chunk(n+1); |
Line 2279 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
Line 2569 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
} |
} |
if (j>m) |
if (j>m) |
{ |
{ |
if (nontriv) { avma = av0; return NULL; } |
if (deplin) { |
|
c = cgetg(n+1, t_COL); |
|
for (i=1; i<k; i++) c[i] = lmodii(gcoeff(x,d[i],k), p); |
|
c[k] = un; for (i=k+1; i<=n; i++) c[i] = zero; |
|
return gerepileupto(av0, c); |
|
} |
r++; d[k]=0; |
r++; d[k]=0; |
for(j=1; j<k; j++) |
for(j=1; j<k; j++) |
if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k)); |
if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k)); |
Line 2291 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
Line 2586 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p); |
coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p); |
for (t=1; t<=m; t++) |
for (t=1; t<=m; t++) |
if (t!=j) |
{ |
{ |
if (t==j) continue; |
piv = modii(gcoeff(x,t,k), p); |
|
if (signe(piv)) |
piv = modii(gcoeff(x,t,k), p); |
{ |
if (!signe(piv)) continue; |
coeff(x,t,k)=zero; |
|
for (i=k+1; i<=n; i++) |
coeff(x,t,k)=zero; |
coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i))); |
for (i=k+1; i<=n; i++) |
if (low_stack(lim, stack_lim(av,1))) |
coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i))); |
gerepile_gauss_FpM_ker(x,p,m,n,k,t,av); |
if (low_stack(lim, stack_lim(av,1))) |
} |
gerepile_gauss_FpM_ker(x,p,k,t,av); |
} |
} |
} |
} |
} |
} |
|
if (deplin) { avma = av0; return NULL; } |
|
|
tetpil=avma; y=cgetg(r+1,t_MAT); |
tetpil=avma; y=cgetg(r+1,t_MAT); |
for (j=k=1; j<=r; j++,k++) |
for (j=k=1; j<=r; j++,k++) |
Line 2325 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
Line 2621 FpM_ker_i(GEN x, GEN p, long nontriv) |
|
return gerepile(av0,tetpil,y); |
return gerepile(av0,tetpil,y); |
} |
} |
|
|
|
GEN |
|
FpM_intersect(GEN x, GEN y, GEN p) |
|
{ |
|
gpmem_t av = avma; |
|
long j, lx = lg(x); |
|
GEN z; |
|
|
|
if (lx==1 || lg(y)==1) return cgetg(1,t_MAT); |
|
z = FpM_ker(concatsp(x,y), p); |
|
for (j=lg(z)-1; j; j--) setlg(z[j],lx); |
|
return gerepileupto(av, FpM_mul(x,z,p)); |
|
} |
|
|
|
GEN |
|
FpM_ker(GEN x, GEN p) { return FpM_ker_i(x, p, 0); } |
|
GEN |
|
FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x, p, 1); } |
|
/* not memory clean */ |
|
GEN |
|
u_FpM_ker(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 0); } |
|
GEN |
|
u_FpM_deplin(GEN x, ulong p) { return u_FpM_ker_sp(dummycopy(x), p, 1); } |
|
|
static void |
static void |
FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
{ |
{ |
|
gpmem_t av,lim; |
GEN c,d,piv; |
GEN c,d,piv; |
long i,j,k,r,t,n,m,av,lim; |
long i,j,k,r,t,n,m; |
|
|
if (!p) return gauss_pivot(x,dd,rr); |
if (!p) { gauss_pivot(x,dd,rr); return; } |
if (typ(x)!=t_MAT) err(typeer,"FpM_gauss_pivot"); |
if (typ(x)!=t_MAT) err(typeer,"FpM_gauss_pivot"); |
n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; } |
n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; } |
|
|
Line 2363 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
|
Line 2683 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
|
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i))); |
coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i))); |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
gerepile_gauss(x,m,n,k,t,av,j,c); |
gerepile_gauss(x,k,t,av,j,c); |
} |
} |
} |
} |
for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ |
for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ |
Line 2371 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
|
Line 2691 FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr) |
|
} |
} |
*dd=d; *rr=r; |
*dd=d; *rr=r; |
} |
} |
|
static void |
GEN |
FqM_gauss_pivot(GEN x, GEN T, GEN p, GEN *dd, long *rr) |
FpM_ker(GEN x, GEN p) |
|
{ |
{ |
return FpM_ker_i(x, p, 0); |
gpmem_t av,lim; |
} |
GEN c,d,piv; |
|
long i,j,k,r,t,n,m; |
|
|
int |
if (typ(x)!=t_MAT) err(typeer,"FqM_gauss_pivot"); |
ker_trivial_mod_p(GEN x, GEN p) |
n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; } |
{ |
|
return FpM_ker_i(x, p, 1)==NULL; |
m=lg(x[1])-1; r=0; |
|
x=dummycopy(x); |
|
c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; |
|
d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1); |
|
for (k=1; k<=n; k++) |
|
{ |
|
for (j=1; j<=m; j++) |
|
if (!c[j]) |
|
{ |
|
coeff(x,j,k) = (long)FpX_res(gcoeff(x,j,k), T,p); |
|
if (signe(coeff(x,j,k))) break; |
|
} |
|
if (j>m) { r++; d[k]=0; } |
|
else |
|
{ |
|
c[j]=k; d[k]=j; piv = gneg(FpXQ_inv(gcoeff(x,j,k), T,p)); |
|
for (i=k+1; i<=n; i++) |
|
coeff(x,j,i) = (long)Fq_mul(piv,gcoeff(x,j,i), T, p); |
|
for (t=1; t<=m; t++) |
|
if (!c[t]) /* no pivot on that line yet */ |
|
{ |
|
piv=gcoeff(x,t,k); |
|
if (signe(piv)) |
|
{ |
|
coeff(x,t,k)=zero; |
|
for (i=k+1; i<=n; i++) |
|
coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(piv,gcoeff(x,j,i))); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
gerepile_gauss(x,k,t,av,j,c); |
|
} |
|
} |
|
for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ |
|
} |
|
} |
|
*dd=d; *rr=r; |
} |
} |
|
|
GEN |
GEN |
FpM_image(GEN x, GEN p) |
FpM_image(GEN x, GEN p) |
{ |
{ |
|
gpmem_t av = avma; |
GEN d,y; |
GEN d,y; |
long j,k,r, av = avma; |
long j,k,r; |
|
|
FpM_gauss_pivot(x,p,&d,&r); |
FpM_gauss_pivot(x,p,&d,&r); |
|
|
Line 2406 FpM_image(GEN x, GEN p) |
|
Line 2761 FpM_image(GEN x, GEN p) |
|
static GEN |
static GEN |
sFpM_invimage(GEN mat, GEN y, GEN p) |
sFpM_invimage(GEN mat, GEN y, GEN p) |
{ |
{ |
long av=avma,i, nbcol = lg(mat); |
gpmem_t av=avma; |
|
long i, nbcol = lg(mat); |
GEN col,p1 = cgetg(nbcol+1,t_MAT),res; |
GEN col,p1 = cgetg(nbcol+1,t_MAT),res; |
|
|
if (nbcol==1) return NULL; |
if (nbcol==1) return NULL; |
Line 2434 sFpM_invimage(GEN mat, GEN y, GEN p) |
|
Line 2790 sFpM_invimage(GEN mat, GEN y, GEN p) |
|
GEN |
GEN |
FpM_invimage(GEN m, GEN v, GEN p) |
FpM_invimage(GEN m, GEN v, GEN p) |
{ |
{ |
long av=avma,j,lv,tv=typ(v); |
gpmem_t av=avma; |
|
long j,lv,tv=typ(v); |
GEN y,p1; |
GEN y,p1; |
|
|
if (typ(m)!=t_MAT) err(typeer,"inverseimage"); |
if (typ(m)!=t_MAT) err(typeer,"inverseimage"); |
Line 2455 FpM_invimage(GEN m, GEN v, GEN p) |
|
Line 2812 FpM_invimage(GEN m, GEN v, GEN p) |
|
} |
} |
return y; |
return y; |
} |
} |
/************************************************************** |
/************************************************************** |
Rather stupid implementation of linear algebra in finite fields |
Rather stupid implementation of linear algebra in finite fields |
**************************************************************/ |
**************************************************************/ |
static GEN |
static GEN |
Line 2467 Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) |
|
Line 2824 Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) |
|
case 1: return FpX_Fp_add(x,y,p); |
case 1: return FpX_Fp_add(x,y,p); |
case 2: return FpX_Fp_add(y,x,p); |
case 2: return FpX_Fp_add(y,x,p); |
case 3: return FpX_add(x,y,p); |
case 3: return FpX_add(x,y,p); |
} |
} |
return NULL; |
return NULL; |
} |
} |
#if 0 |
#if 0 |
Line 2475 Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) |
|
Line 2832 Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) |
|
* For consistency we write the code there. |
* For consistency we write the code there. |
* To avoid having to remove static status, we rewrite it in polarit1.c |
* To avoid having to remove static status, we rewrite it in polarit1.c |
*/ |
*/ |
static GEN |
static GEN |
Fq_neg(GEN x, GEN T, GEN p) |
Fq_neg(GEN x, GEN T, GEN p) |
{ |
{ |
switch(typ(x)==t_POL) |
switch(typ(x)==t_POL) |
Line 2487 Fq_neg(GEN x, GEN T, GEN p) |
|
Line 2844 Fq_neg(GEN x, GEN T, GEN p) |
|
} |
} |
#endif |
#endif |
|
|
static GEN |
GEN |
Fq_mul(GEN x, GEN y, GEN T, GEN p) |
Fq_mul(GEN x, GEN y, GEN T, GEN p) |
{ |
{ |
switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1)) |
switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1)) |
Line 2514 Fq_neg_inv(GEN x, GEN T, GEN p) |
|
Line 2871 Fq_neg_inv(GEN x, GEN T, GEN p) |
|
static GEN |
static GEN |
Fq_res(GEN x, GEN T, GEN p) |
Fq_res(GEN x, GEN T, GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
switch(typ(x)==t_POL) |
switch(typ(x)==t_POL) |
{ |
{ |
case 0: return modii(x,p); |
case 0: return modii(x,p); |
Line 2524 Fq_res(GEN x, GEN T, GEN p) |
|
Line 2881 Fq_res(GEN x, GEN T, GEN p) |
|
} |
} |
|
|
static void |
static void |
Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, long n, long k, long t, |
Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long k, long t, gpmem_t av) |
ulong av) |
|
{ |
{ |
ulong tetpil = avma,A; |
gpmem_t tetpil = avma,A; |
long dec,u,i; |
long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; |
|
size_t dec; |
|
|
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
Line 2540 Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, lon |
|
Line 2897 Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, lon |
|
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
(void)gerepile(av,tetpil,NULL); dec = av-tetpil; |
for (u=t+1; u<=m; u++) |
for (u=t+1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,k); |
A=(gpmem_t)coeff(x,u,k); |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
if (A<av && A>=bot) coeff(x,u,k)+=dec; |
} |
} |
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
for (u=1; u<=m; u++) |
for (u=1; u<=m; u++) |
{ |
{ |
A=coeff(x,u,i); |
A=(gpmem_t)coeff(x,u,i); |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
if (A<av && A>=bot) coeff(x,u,i)+=dec; |
} |
} |
} |
} |
|
|
static GEN |
static GEN |
FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) |
FqM_ker_i(GEN x, GEN T, GEN p, long deplin) |
{ |
{ |
|
gpmem_t av0,av,lim,tetpil; |
GEN y,c,d,piv,mun; |
GEN y,c,d,piv,mun; |
long i,j,k,r,t,n,m,av0,av,lim,tetpil; |
long i,j,k,r,t,n,m; |
|
|
|
if (!T) return FpM_ker_i(x,p,deplin); |
|
|
if (typ(x)!=t_MAT) err(typeer,"FpM_ker"); |
if (typ(x)!=t_MAT) err(typeer,"FpM_ker"); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
n=lg(x)-1; if (!n) return cgetg(1,t_MAT); |
|
|
Line 2575 FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) |
|
Line 2935 FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) |
|
} |
} |
if (j>m) |
if (j>m) |
{ |
{ |
if (nontriv) { avma = av0; return NULL; } |
if (deplin) { |
|
c = cgetg(n+1, t_COL); |
|
for (i=1; i<k; i++) c[i] = (long)Fq_res(gcoeff(x,d[i],k), T, p); |
|
c[k] = un; for (i=k+1; i<=n; i++) c[i] = zero; |
|
return gerepileupto(av0, c); |
|
} |
r++; d[k]=0; |
r++; d[k]=0; |
for(j=1; j<k; j++) |
for(j=1; j<k; j++) |
if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k)); |
if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k)); |
Line 2587 FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) |
|
Line 2952 FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) |
|
for (i=k+1; i<=n; i++) |
for (i=k+1; i<=n; i++) |
coeff(x,j,i) = (long) Fq_mul(piv,gcoeff(x,j,i), T, p); |
coeff(x,j,i) = (long) Fq_mul(piv,gcoeff(x,j,i), T, p); |
for (t=1; t<=m; t++) |
for (t=1; t<=m; t++) |
if (t!=j) |
{ |
{ |
if (t==j) continue; |
piv = Fq_res(gcoeff(x,t,k), T, p); |
|
if (signe(piv)) |
piv = Fq_res(gcoeff(x,t,k), T, p); |
{ |
if (!signe(piv)) continue; |
coeff(x,t,k)=zero; |
|
for (i=k+1; i<=n; i++) |
coeff(x,t,k)=zero; |
coeff(x,t,i) = (long)Fq_add(gcoeff(x,t,i),Fq_mul(piv,gcoeff(x,j,i), T, p), T, p); |
for (i=k+1; i<=n; i++) |
if (low_stack(lim, stack_lim(av,1))) |
coeff(x,t,i) = (long)Fq_add(gcoeff(x,t,i),Fq_mul(piv,gcoeff(x,j,i), T, p), T, p); |
Fq_gerepile_gauss_ker(x,T,p,m,n,k,t,av); |
if (low_stack(lim, stack_lim(av,1))) |
} |
Fq_gerepile_gauss_ker(x,T,p,k,t,av); |
} |
} |
} |
} |
} |
} |
|
if (deplin) { avma = av0; return NULL; } |
|
|
tetpil=avma; y=cgetg(r+1,t_MAT); |
tetpil=avma; y=cgetg(r+1,t_MAT); |
for (j=k=1; j<=r; j++,k++) |
for (j=k=1; j<=r; j++,k++) |
Line 2638 eigen(GEN x, long prec) |
|
Line 3004 eigen(GEN x, long prec) |
|
{ |
{ |
GEN y,rr,p,ssesp,r1,r2,r3; |
GEN y,rr,p,ssesp,r1,r2,r3; |
long e,i,k,l,ly,ex, n = lg(x); |
long e,i,k,l,ly,ex, n = lg(x); |
ulong av = avma; |
gpmem_t av = avma; |
|
|
if (typ(x)!=t_MAT) err(typeer,"eigen"); |
if (typ(x)!=t_MAT) err(typeer,"eigen"); |
if (n != 1 && n != lg(x[1])) err(mattype1,"eigen"); |
if (n != 1 && n != lg(x[1])) err(mattype1,"eigen"); |
Line 2696 det0(GEN a,long flag) |
|
Line 3062 det0(GEN a,long flag) |
|
|
|
/* Exact types: choose the first non-zero pivot. Otherwise: maximal pivot */ |
/* Exact types: choose the first non-zero pivot. Otherwise: maximal pivot */ |
static GEN |
static GEN |
det_simple_gauss(GEN a, long inexact) |
det_simple_gauss(GEN a, int inexact) |
{ |
{ |
long i,j,k,av,av1,s, nbco = lg(a)-1; |
gpmem_t av, av1; |
|
long i,j,k,s, nbco = lg(a)-1; |
GEN x,p; |
GEN x,p; |
|
|
av=avma; s=1; x=gun; a=dummycopy(a); |
av=avma; s=1; x=gun; a=dummycopy(a); |
Line 2745 det_simple_gauss(GEN a, long inexact) |
|
Line 3112 det_simple_gauss(GEN a, long inexact) |
|
av1=avma; return gerepile(av,av1,gmul(x,gcoeff(a,nbco,nbco))); |
av1=avma; return gerepile(av,av1,gmul(x,gcoeff(a,nbco,nbco))); |
} |
} |
|
|
/* a has integer entries, N = P^n */ |
|
GEN |
GEN |
det_mod_P_n(GEN a, GEN N, GEN P) |
|
{ |
|
long va,i,j,k,s, av = avma, nbco = lg(a)-1; |
|
GEN x,p; |
|
|
|
s=1; va=0; x=gun; a=dummycopy(a); |
|
p = NULL; /* for gcc -Wall */ |
|
for (i=1; i<nbco; i++) |
|
{ |
|
long fl = 0; |
|
for(;;) |
|
{ |
|
for (k=i; k<=nbco; k++) |
|
{ |
|
p=gcoeff(a,i,k); |
|
if (signe(p)) |
|
{ |
|
fl = 1; |
|
if (resii(p,P) != gzero) break; |
|
} |
|
} |
|
if (k <= nbco) break; |
|
va++; N = divii(N, P); |
|
if (!fl || is_pm1(N)) { avma=av; return gzero; } |
|
for (k=i; k<=nbco; k++) coeff(a,i,k) = ldivii(gcoeff(a,i,k), P); |
|
} |
|
if (k != i) { swap(a[i],a[k]); s = -s; } |
|
|
|
x = gmul(x,p); p = mpinvmod(p,N); |
|
for (k=i+1; k<=nbco; k++) |
|
{ |
|
GEN m = resii(gcoeff(a,i,k), N); |
|
coeff(a,i,k) = zero; |
|
if (signe(m)) |
|
{ |
|
m = negi(resii(mulii(m,p), N)); |
|
for (j=i+1; j<=nbco; j++) |
|
coeff(a,j,k) = lresii(addii(gcoeff(a,j,k), |
|
mulii(m,gcoeff(a,j,i))), N); |
|
} |
|
} |
|
} |
|
if (s<0) x = negi(x); |
|
x = resii(mulii(x,gcoeff(a,nbco,nbco)), N); |
|
return gerepileuptoint(av, mulii(x, gpowgs(P,va))); |
|
} |
|
|
|
GEN |
|
det2(GEN a) |
det2(GEN a) |
{ |
{ |
long nbco = lg(a)-1; |
long nbco = lg(a)-1; |
if (typ(a)!=t_MAT) err(mattype1,"det2"); |
if (typ(a)!=t_MAT) err(mattype1,"det2"); |
if (!nbco) return gun; |
if (!nbco) return gun; |
if (nbco != lg(a[1])-1) err(mattype1,"det2"); |
if (nbco != lg(a[1])-1) err(mattype1,"det2"); |
return det_simple_gauss(a,use_maximal_pivot(a)); |
return det_simple_gauss(a, use_maximal_pivot(a)); |
} |
} |
|
|
static GEN |
static GEN |
Line 2818 mydiv(GEN x, GEN y) |
|
Line 3136 mydiv(GEN x, GEN y) |
|
GEN |
GEN |
det(GEN a) |
det(GEN a) |
{ |
{ |
ulong av, lim; |
gpmem_t av, lim; |
long nbco = lg(a)-1,i,j,k,s; |
long nbco = lg(a)-1,i,j,k,s; |
GEN p,pprec; |
GEN p,pprec; |
|
|
|
|
if (!nbco) return gun; |
if (!nbco) return gun; |
if (nbco != lg(a[1])-1) err(mattype1,"det"); |
if (nbco != lg(a[1])-1) err(mattype1,"det"); |
if (use_maximal_pivot(a)) return det_simple_gauss(a,1); |
if (use_maximal_pivot(a)) return det_simple_gauss(a,1); |
if (DEBUGLEVEL > 7) timer2(); |
if (DEBUGLEVEL > 7) (void)timer2(); |
|
|
av = avma; lim = stack_lim(av,2); |
av = avma; lim = stack_lim(av,2); |
a = dummycopy(a); s = 1; |
a = dummycopy(a); s = 1; |
for (pprec=gun,i=1; i<nbco; i++,pprec=p) |
for (pprec=gun,i=1; i<nbco; i++,pprec=p) |
{ |
{ |
|
|
static GEN |
static GEN |
gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1) |
gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1) |
{ |
{ |
long n,m,i,j,lM,av=avma,tetpil; |
gpmem_t av = avma, tetpil; |
|
long n,m,i,j,lM; |
GEN p1,delta,H,U,u1,u2,x; |
GEN p1,delta,H,U,u1,u2,x; |
|
|
if (typ(M)!=t_MAT) err(typeer,"gaussmodulo"); |
if (typ(M)!=t_MAT) err(typeer,"gaussmodulo"); |
Line 2950 gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1) |
|
Line 3269 gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1) |
|
GEN |
GEN |
matsolvemod0(GEN M, GEN D, GEN Y, long flag) |
matsolvemod0(GEN M, GEN D, GEN Y, long flag) |
{ |
{ |
long av; |
gpmem_t av; |
GEN p1,y; |
GEN p1,y; |
|
|
if (!flag) return gaussmoduloall(M,D,Y,NULL); |
if (!flag) return gaussmoduloall(M,D,Y,NULL); |