OpenXM_contrib/pari-2.2/src/basemath/alglin1.c - diff

Return to alglin1.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.c between version 1.1 and 1.2

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,

/* for linear algebra mod p */

#ifdef LONG_IS_64BIT

# define MASK (0x7fffffff00000000UL)

# define MASK (0xffffffff00000000UL)

#else

# define MASK (0x7fff0000UL)

# define MASK (0xffff0000UL)

#endif

/* 2p^2 < 2^BITS_IN_LONG */

Line 36 Foundation, Inc., 59 Temple Place - Suite 330, Boston,

#endif

#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))

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);

/*******************************************************************/

/* */

Line 75 gtrans_i(GEN x)

return y;

}

GEN

gtrans(GEN x)

{

Line 135 GEN

Line 134 GEN

concatsp3(GEN x, GEN y, GEN 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<ly; i++) *++t = *++y;

for (i=1; i<lz; i++) *++t = *++z;

Line 149 vconcat(GEN A, GEN B)

long la,ha,hb,hc,i,j;

GEN M,a,b,c;

if (!A) return B;

if (!B) return A;

la = lg(A); if (la==1) return A;

ha = lg(A[1]); M = cgetg(la,t_MAT);

hb = lg(B[1]); hc = ha+hb-1;

Line 162 vconcat(GEN A, GEN B)

Line 164 vconcat(GEN A, GEN B)

}

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; }

GEN

_col(GEN x) { GEN v = cgetg(2, t_COL); v[1] = (long)x; return v; }

/* routines for naive growarrays */

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)

{

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)

{

ulong av = avma;

gpmem_t av = avma;

if (tx == t_LIST)

{ lx = lgef(x); i = 2; }

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;

if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;

if (*s == 0) return 0;

if (!*s) return 0;

if (*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;

}

/* A[x0,] */

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)

{

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

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;

if (! is_matvec_t(tx)) err(typeer,"extract");

Line 585 extract(GEN x, GEN l)

Line 623 extract(GEN x, GEN l)

GEN

matextract(GEN x, GEN l1, GEN l2)

{

long av = avma, tetpil;

gpmem_t av = avma, tetpil;

if (typ(x)!=t_MAT) err(typeer,"matextract");

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);

}

/* 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 */

Line 668 GEN

Line 719 GEN

zeromat(long m, long n)

{

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;

}

GEN

gscalcol(GEN x, long n)

{

GEN y=gscalcol_proto(gzero,gzero,n);

if (n) y[1]=lcopy(x);

return y;

}

Line 847 mattodiagonal(GEN m)

Line 899 mattodiagonal(GEN m)

GEN

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;

}

/* 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) */

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; }

/* 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;

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;

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

check_b(GEN b, long nbli)

check_b(GEN b, long nbli, int *iscol)

{

GEN col;

*iscol = (b && typ(b) == t_COL);

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");

return dummycopy(b);

return b;

}

/* 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++)

{

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);

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--)

{

av = avma; m = mulii(negi((GEN)b[i]),t);

for (j=i+1; j<=n; 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;

}

/* 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

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)

{

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);

for (i=li-1; i>0; i--)

{

gpmem_t av = avma;

m = gneg_i((GEN)b[i]);

for (j=i+1; j<=li; 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;

}

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);

for (i=li-1; i>0; i--)

{

gpmem_t av = avma;

m = (GEN)b[i];

for (j=i+1; j<=li; j++)

m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j]));

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;

}

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 */

static long

u_Fp_inv(long a, long p)

{

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

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;

}

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]));

}

/* same, reduce mod (T,p) */

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

* 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.

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

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;

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);

return cgetg(1,t_MAT);

}

av=avma; lim=stack_lim(av,1);

av = avma; lim = stack_lim(av,1);

li = lg(a[1])-1;

if (li != aco && (li < aco || b)) err(mattype1,"gauss");

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);

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 */

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)

{

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);

}

if (i == aco) break;

for (k=i+1; k<=li; k++)

{

m=gcoeff(a,k,i);

m = gcoeff(a,k,i);

if (!gcmp0(m))

{

m = gneg_i(gdiv(m,p));

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)))

{

GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b;

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 (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

Line 1138 gauss(GEN a, GEN b)

Line 1308 gauss(GEN a, GEN b)

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;

if (!aco) return cgetg(1,t_MAT);

li = lg(a[1])-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 */

for (i=1; i<=aco; i++)

{

/* 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)

{

for (k++; k <= li; k++)

{

coeff(a,k,i) %= p;

if (coeff(a,k,i)) break;

}

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)

{

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);

}

if (i == aco) break;

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)

{

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

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)

{

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;

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;

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))

{

ulong pp=p[2];

a = u_Fp_FpM(a, 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);

a = dummycopy(a);

bco = lg(b)-1;

iscol = (typ(b)==t_COL);

piv = NULL; /* gcc -Wall */

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);

if (signe(coeff(a,k,i))) break;

}

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)

{

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);

}

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 = resii(negi(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)))

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 (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

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 */

#define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG)

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' */

GEN

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;

long lM = lg(M), stable = 0;

Line 1346 ZM_inv(GEN M, GEN dM)

Line 1587 ZM_inv(GEN M, GEN dM)

av2 = avma;

H = NULL;

d += 3000; /* 27449 = prime(3000) */

for(p = 27449; ; p+= *d++)

for(p = 27449; ; )

{

if (!*d) err(primer1);

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 (!H)

{

H = ZM_init_CRT(Hp, 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 (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */

if (low_stack(lim, stack_lim(av,2)))

{

GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q;

if (DEBUGMEM>1) err(warnmem,"ZM_inv");

gerepilemany(av2,gptr, 2);

}

repeat:

NEXT_PRIME_VIADIFF_CHECK(p,d);

}

if (DEBUGLEVEL>5) msgtimer("ZM_inv done");

return gerepilecopy(av, H);

}

/* same as above, M rational */

GEN

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

Line 1397 QM_inv(GEN M, GEN dM)

Line 1637 QM_inv(GEN M, GEN dM)

GEN

detint(GEN x)

{

gpmem_t av = avma, av1, lim;

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");

n=lg(x)-1; if (!n) return gun;

Line 1462 detint(GEN x)

Line 1703 detint(GEN x)

}

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);

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;

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;

}

for (i=k+1; i<=n; i++)

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;

}

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);

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;

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;

}

for (i=k+1; i<=n; i++)

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;

}

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 */

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);

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++)

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;

}

for (u=1; u<=m; u++)

if (u==j || !c[u])

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;

}

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 */

/* */

/*******************************************************************/

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 */

GEN

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");

n=lg(x)-1; if (!n) return cgetg(1,t_MAT);

Line 1564 keri(GEN x)

Line 1860 keri(GEN x)

av0=avma; m=lg(x[1])-1; r=0;

pp=cgetg(n+1,t_COL);

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++)

{

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++)

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);

}

else

{

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++)

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++)

{

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)))

{

p1 = gclone(p);

gerepile_gauss_ker(x,m,n,k,t,av);

gerepile_gauss_ker(x,k,t,av);

p = gcopy(p1); gunclone(p1);

}

Line 1612 keri(GEN x)

Line 1906 keri(GEN x)

for (j=k=1; j<=r; j++,k++)

{

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++)

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);

}

else

Line 1628 keri(GEN x)

Line 1922 keri(GEN x)

}

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 (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));

}

/*******************************************************************/

Line 1681 deplin(GEN x)

Line 1978 deplin(GEN x)

/* (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).

* d[k] contains the index of the first non-zero pivot in column k

* If a != NULL, use x - a Id instead (for eigen)

Line 1761 static GEN

Line 1986 static GEN

gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr)

{

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);

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++)

coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));

if (low_stack(lim, stack_lim(av,1)))

gerepile_gauss_ker(x,m,n,k,t,av);

gerepile_gauss_ker(x,k,t,av);

}

Line 1815 static void

Line 2041 static void

gauss_pivot(GEN x0, GEN *dd, long *rr)

{

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);

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++)

coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i)));

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 */

}

Line 1868 gauss_pivot(GEN x0, GEN *dd, long *rr)

Line 2095 gauss_pivot(GEN x0, GEN *dd, long *rr)

static GEN

ker0(GEN x, GEN a)

{

gpmem_t av = avma, tetpil;

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);

if (!r) { avma=av; return cgetg(1,t_MAT); }

Line 1907 matker0(GEN x,long flag)

Line 2135 matker0(GEN x,long flag)

GEN

image(GEN x)

{

gpmem_t av = avma;

GEN d,y;

long j,k,r, av = avma;

long j,k,r;

gauss_pivot(x,&d,&r);

Line 1926 image(GEN x)

Line 2155 image(GEN x)

GEN

imagecompl(GEN x)

{

gpmem_t av = avma;

GEN d,y;

long j,i,r,av = avma;

long j,i,r;

gauss_pivot(x,&d,&r);

avma=av; y=cgetg(r+1,t_VEC);

Line 1940 imagecompl(GEN x)

Line 2170 imagecompl(GEN x)

GEN

imagecomplspec(GEN x, long *nlze)

{

gpmem_t av = avma;

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);

avma=av; y = cgetg(l,t_VECSMALL);

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

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);

if (nbcol==1) return NULL;

Line 1977 sinverseimage(GEN mat, GEN y)

Line 2209 sinverseimage(GEN mat, GEN y)

GEN

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;

if (typ(m)!=t_MAT) err(typeer,"inverseimage");

Line 1999 inverseimage(GEN m,GEN v)

Line 2232 inverseimage(GEN m,GEN v)

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

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

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

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)

{

long av=avma,tetpil,k,n,i;

gpmem_t av=avma,tetpil;

long k,n,i;

GEN p1,p2;

if (typ(x)!=t_MAT) err(typeer,"image2");

Line 2068 matimage0(GEN x,long flag)

Line 2346 matimage0(GEN x,long flag)

long

rank(GEN x)

{

long av = avma, r;

gpmem_t av = avma;

long r;

GEN d;

gauss_pivot(x,&d,&r);

Line 2078 rank(GEN x)

Line 2357 rank(GEN x)

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 */

static GEN

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;

/* 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;

}

GEN

indexrank(GEN x) { return indexrank0(x,NULL,0); }

GEN

sindexrank(GEN x) { return indexrank0(x,NULL,1); }

GEN

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 */

/* */

/*******************************************************************/

/*If p is NULL no reduction is performed.*/

GEN

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);

GEN z;

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);

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);

for (i=1; i<l; i++)

{

ulong av;

gpmem_t av;

GEN p1,p2;

int k;

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

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;

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);

l=lg(x[1]);

l = lg(x[1]);

z = cgetg(l,t_COL);

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++)

{

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;

}

/* in place, destroy x */

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;

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;

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 (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;

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

{

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++)

coeff(x,j,i) = (piv * coeff(x,j,i)) % p;

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;

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);

}

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++)

{

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++)

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

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

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;

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");

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);

c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;

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 (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;

for(j=1; j<k; j++)

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++)

coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);

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);

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);

}

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

FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)

{

gpmem_t av,lim;

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");

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++)

coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i)));

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 */

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;

}

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

FpM_image(GEN x, GEN p)

{

gpmem_t av = avma;

GEN d,y;

long j,k,r, av = avma;

long j,k,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

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;

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

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;

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;

}

/**************************************************************

Rather stupid implementation of linear algebra in finite fields

**************************************************************/

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 2: return FpX_Fp_add(y,x,p);

case 3: return FpX_add(x,y,p);

}

return NULL;

}

#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.

* To avoid having to remove static status, we rewrite it in polarit1.c

static GEN

Fq_neg(GEN x, GEN T, GEN p)

{

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

static GEN

GEN

Fq_mul(GEN x, GEN y, GEN T, GEN p)

{

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

Fq_res(GEN x, GEN T, GEN p)

{

ulong ltop=avma;

gpmem_t ltop=avma;

switch(typ(x)==t_POL)

{

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

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);

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;

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;

}

for (i=k+1; i<=n; i++)

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;

}

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;

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");

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 (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;

for(j=1; j<k; j++)

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++)

coeff(x,j,i) = (long) Fq_mul(piv,gcoeff(x,j,i), T, p);

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);

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;

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 (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 */

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;

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)));

}

/* a has integer entries, N = P^n */

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)

{

long nbco = lg(a)-1;

if (typ(a)!=t_MAT) err(mattype1,"det2");

if (!nbco) return gun;

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

Line 2818 mydiv(GEN x, GEN y)

Line 3136 mydiv(GEN x, GEN y)

GEN

det(GEN a)

{

ulong av, lim;

gpmem_t av, lim;

long nbco = lg(a)-1,i,j,k,s;

GEN p,pprec;

Line 2826 det(GEN a)

Line 3144 det(GEN a)

if (!nbco) return gun;

if (nbco != lg(a[1])-1) err(mattype1,"det");

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);

a = dummycopy(a); s = 1;

for (pprec=gun,i=1; i<nbco; i++,pprec=p)

{

Line 2892 det(GEN a)

Line 3210 det(GEN a)

static GEN

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;

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

matsolvemod0(GEN M, GEN D, GEN Y, long flag)

{

long av;

gpmem_t av;

GEN p1,y;

if (!flag) return gaussmoduloall(M,D,Y,NULL);

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>