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

Return to alglin2.c CVS log

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

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.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 43 charpoly0(GEN x, int v, long flag)

static GEN

caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,GEN*))

{

long av = avma, d;

gpmem_t av = avma;

long d;

GEN p1, p2 = leading_term(p);

if (!signe(x)) return gpowgs(polx[v], degpol(p));

Line 55 caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,

Line 56 caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,

else

p1 = gsubst(p1,MAXVARN,polx[v]);

if (!gcmp1(p2) && (d=degpol(x)) > 0) p1 = gdiv(p1, gpuigs(p2,d));

if (!gcmp1(p2) && (d=degpol(x)) > 0) p1 = gdiv(p1, gpowgs(p2,d));

return gerepileupto(av,p1);

}

Line 75 caractducos(GEN p, GEN x, int v)

Line 76 caractducos(GEN p, GEN x, int v)

static GEN

easychar(GEN x, int v, GEN *py)

{

long l,tetpil,lx;

gpmem_t av;

long lx;

GEN p1,p2;

switch(typ(x))

Line 96 easychar(GEN x, int v, GEN *py)

Line 98 easychar(GEN x, int v, GEN *py)

case t_COMPLEX: case t_QUAD:

if (py) err(typeer,"easychar");

p1=cgetg(5,t_POL);

p1 = cgetg(5,t_POL);

p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v);

p1[1] = evalsigne(1) | evallgef(5) | evalvarn(v);

p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma;

p1[2] = lnorm(x); av = avma;

p1[3]=lpile(l,tetpil,gneg(p2));

p1[3] = lpileupto(av, gneg(gtrace(x)));

p1[4]=un; return p1;

p1[4] = un; return p1;

case t_POLMOD:

if (py) err(typeer,"easychar");

Line 123 easychar(GEN x, int v, GEN *py)

Line 125 easychar(GEN x, int v, GEN *py)

GEN

caract(GEN x, int v)

{

long n,k,l=avma,tetpil;

long k, n;

GEN p1,p2,p3,p4,p5,p6;

gpmem_t av=avma;

GEN p1,p2,p3,p4,x_k;

if ((p1 = easychar(x,v,NULL))) return p1;

p1=gzero; p2=gun;

p1 = gzero; p2 = gun;

n=lg(x)-1; if (n&1) p2=gneg_i(p2);

n = lg(x)-1; if (n&1) p2 = negi(p2);

p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5;

x_k = dummycopy(polx[v]);

p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3);

p4 = cgetg(3,t_RFRACN); p4[2] = (long)x_k;

for (k=0; k<=n; k++)

{

p3=det(gsub(gscalmat(stoi(k),n), x));

p3 = det(gsub(gscalmat(stoi(k),n), x));

p4[1]=lmul(p3,p2); p6[2]=k;

p4[1] = lmul(p3,p2); x_k[2] = lstoi(-k);

p1=gadd(p4,p1); p5[2]=(long)p6;

p1 = gadd(p4,p1);

if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);

if (k == n) break;

p2 = gdivgs(gmulsg(k-n,p2),k+1);

}

p2=mpfact(n); tetpil=avma;

return gerepileupto(av, gdiv((GEN)p1[1], mpfact(n)));

return gerepile(l,tetpil,gdiv((GEN) p1[1],p2));

}

GEN

Line 155 caradj0(GEN x, long v)

Line 159 caradj0(GEN x, long v)

GEN

caradj(GEN x, long v, GEN *py)

{

long i,j,k,l,av,tetpil;

gpmem_t av,tetpil;

long i,j,k,l;

GEN p,y,t,*gptr[2];

if ((p = easychar(x,v,py))) return p;

Line 217 GEN

Line 222 GEN

adj(GEN x)

{

GEN y;

caradj(x,MAXVARN,&y); return y;

(void)caradj(x,MAXVARN,&y); return y;

}

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

Line 225 adj(GEN x)

Line 230 adj(GEN x)

/* HESSENBERG FORM */

/* */

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

#define swap(x,y) { long _t=x; x=y; y=_t; }

#define lswap(x,y) { long _t=x; x=y; y=_t; }

#define gswap(x,y) { GEN _t=x; x=y; y=_t; }

#define swap(x,y) { GEN _t=x; x=y; y=_t; }

GEN

hess(GEN x)

{

ulong av = avma;

gpmem_t av = avma;

long lx=lg(x),m,i,j;

GEN p,p1,p2;

Line 246 hess(GEN x)

Line 251 hess(GEN x)

p = gcoeff(x,i,m-1);

if (gcmp0(p)) continue;

for (j=m-1; j<lx; j++) swap(coeff(x,i,j), coeff(x,m,j));

for (j=m-1; j<lx; j++) lswap(coeff(x,i,j), coeff(x,m,j));

swap(x[i], x[m]); p = ginv(p);

lswap(x[i], x[m]); p = ginv(p);

for (i=m+1; i<lx; i++)

{

p1 = gcoeff(x,i,m-1);

Line 268 hess(GEN x)

Line 273 hess(GEN x)

GEN

carhess(GEN x, long v)

{

long av,tetpil,lx,r,i;

gpmem_t av,tetpil;

long lx,r,i;

GEN *y,p1,p2,p3,p4;

if ((p1 = easychar(x,v,NULL))) return p1;

Line 299 carhess(GEN x, long v)

Line 305 carhess(GEN x, long v)

GEN

gnorm(GEN x)

{

long l,lx,i,tetpil, tx=typ(x);

gpmem_t av;

long lx,i, tx=typ(x);

GEN p1,p2,y;

switch(tx)

Line 313 gnorm(GEN x)

Line 320 gnorm(GEN x)

case t_FRAC: case t_FRACN:

return gsqr(x);

case t_COMPLEX:

case t_COMPLEX: av = avma;

l=avma; p1=gsqr((GEN) x[1]); p2=gsqr((GEN) x[2]);

return gerepileupto(av, gadd(gsqr((GEN)x[1]), gsqr((GEN)x[2])));

tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));

case t_QUAD:

case t_QUAD: av = avma;

l=avma; p1=(GEN)x[1];

p1 = (GEN)x[1];

p2=gmul((GEN) p1[2], gsqr((GEN) x[3]));

p2 = gmul((GEN)p1[2], gsqr((GEN)x[3]));

p1 = gcmp0((GEN) p1[3])? gsqr((GEN) x[2])

p1 = gcmp0((GEN)p1[3])? gsqr((GEN)x[2])

: gmul((GEN) x[2], gadd((GEN) x[2],(GEN) x[3]));

: gmul((GEN)x[2], gadd((GEN)x[2],(GEN)x[3]));

tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));

return gerepileupto(av, gadd(p1,p2));

case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:

case t_POL: case t_SER: case t_RFRAC: case t_RFRACN: av = avma;

l=avma; p1=gmul(gconj(x),x); tetpil=avma;

return gerepileupto(av, greal(gmul(gconj(x),x)));

return gerepile(l,tetpil,greal(p1));

case t_POLMOD:

p1=(GEN)x[1]; p2=leading_term(p1);

{

if (gcmp1(p2) || gcmp0((GEN) x[2])) return subres(p1,(GEN) x[2]);

GEN T = (GEN)x[1], A = (GEN)x[2];

l=avma; y=subres(p1,(GEN)x[2]); p1=gpuigs(p2,degpol(x[2]));

if (typ(A) != t_POL) return gpowgs(A, degpol(T));

tetpil=avma; return gerepile(l,tetpil,gdiv(y,p1));

y = subres(T, A); p1 = leading_term(T);

if (gcmp1(p1) || gcmp0(A)) return y;

av = avma; T = gpowgs(p1, degpol(A));

return gerepileupto(av, gdiv(y,T));

}

case t_VEC: case t_COL: case t_MAT:

lx=lg(x); y=cgetg(lx,tx);

for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]);

Line 346 gnorm(GEN x)

Line 356 gnorm(GEN x)

GEN

gnorml2(GEN x)

{

gpmem_t av,lim;

GEN y;

long i,tx=typ(x),lx,av,lim;

long i,tx=typ(x),lx;

if (! is_matvec_t(tx)) return gnorm(x);

lx=lg(x); if (lx==1) return gzero;

Line 368 gnorml2(GEN x)

Line 379 gnorml2(GEN x)

GEN

QuickNormL2(GEN x, long prec)

{

long av = avma;

gpmem_t av = avma;

GEN y = gmul(x, realun(prec));

if (typ(x) == t_POL)

*++y = evaltyp(t_VEC) | evallg(lgef(x)-1);

Line 378 QuickNormL2(GEN x, long prec)

Line 389 QuickNormL2(GEN x, long prec)

GEN

gnorml1(GEN x,long prec)

{

ulong av = avma;

gpmem_t av = avma;

long lx,i;

GEN s;

switch(typ(x))

Line 406 gnorml1(GEN x,long prec)

Line 417 gnorml1(GEN x,long prec)

GEN

QuickNormL1(GEN x,long prec)

{

ulong av = avma;

gpmem_t av = avma;

long lx,i;

GEN p1,p2,s;

switch(typ(x))

Line 497 gconj(GEN x)

Line 508 gconj(GEN x)

GEN

conjvec(GEN x,long prec)

{

long lx,s,av,tetpil,i,tx=typ(x);

gpmem_t av,tetpil;

long lx,s,i,tx=typ(x);

GEN z,y,p1,p2,p;

switch(tx)

Line 511 conjvec(GEN x,long prec)

Line 523 conjvec(GEN x,long prec)

case t_VEC: case t_COL:

lx=lg(x); z=cgetg(lx,t_MAT);

for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec);

s=lg(z[1]);

if (lx == 1) break;

s = lg(z[1]);

for (i=2; i<lx; i++)

if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");

break;

case t_POLMOD:

Line 573 assmat(GEN x)

Line 586 assmat(GEN x)

p1[j] = (j==(i+1))? un: zero;

}

p1=cgetg(lx,t_COL); y[i]=(long)p1;

if (gcmp1((GEN) x[lx+1]))

if (gcmp1((GEN)x[lx+1]))

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

p1[j] = lneg((GEN)x[j+1]);

else

{

p2 = (GEN)x[lx+1]; gnegz(p2,p2);

gpmem_t av = avma;

p2 = gclone(gneg((GEN)x[lx+1]));

avma = av;

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

p1[j] = ldiv((GEN)x[j+1],p2);

gnegz(p2,p2);

gunclone(p2);

}

return y;

}

Line 589 assmat(GEN x)

Line 604 assmat(GEN x)

GEN

gtrace(GEN x)

{

long i,l,n,tx=typ(x),lx,tetpil;

gpmem_t av, tetpil;

long i,n,tx=typ(x),lx;

GEN y,p1,p2;

switch(tx)

Line 604 gtrace(GEN x)

Line 620 gtrace(GEN x)

p1=(GEN)x[1];

if (!gcmp0((GEN) p1[3]))

{

l=avma; p2=gmul2n((GEN)x[2],1);

av=avma; p2=gmul2n((GEN)x[2],1);

tetpil=avma; return gerepile(l,tetpil,gadd((GEN)x[3],p2));

tetpil=avma; return gerepile(av,tetpil,gadd((GEN)x[3],p2));

}

return gmul2n((GEN)x[2],1);

Line 620 gtrace(GEN x)

Line 636 gtrace(GEN x)

return y;

case t_POLMOD:

l=avma; n=(lgef(x[1])-4);

av=avma; n=(lgef(x[1])-4);

p1=polsym((GEN)x[1],n); p2=gzero;

for (i=0; i<=n; i++)

p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1]));

return gerepileupto(l,p2);

return gerepileupto(av,p2);

case t_RFRAC: case t_RFRACN:

return gadd(x,gconj(x));

Line 637 gtrace(GEN x)

Line 653 gtrace(GEN x)

case t_MAT:

lx=lg(x); if (lx==1) return gzero;/*now lx>=2*/

if (lx!=lg(x[1])) err(mattype1,"gtrace");

l=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);

av=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);

for (i=2; i<lx-1; i++)

p1=gadd(p1,gcoeff(x,i,i));

tetpil=avma; return gerepile(l,tetpil,gadd(p1,gcoeff(x,i,i)));

tetpil=avma; return gerepile(av,tetpil,gadd(p1,gcoeff(x,i,i)));

}

err(typeer,"gtrace");

return NULL; /* not reached */

}

/* Gauss reduction for positive definite matrix a

/* Cholesky Decomposition for positive definite matrix a

* If a is not positive definite:

* if flag is zero: raise an error

* else: return NULL.

* else: return NULL. */

GEN

sqred1intern(GEN a,long flag)

{

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

GEN b,p;

long i,j,k, n = lg(a);

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

if (typ(a)!=t_MAT) err(typeer,"sqred1");

if (n == 1) return cgetg(1, t_MAT);

Line 704 sqred1(GEN a)

Line 719 sqred1(GEN a)

GEN

sqred3(GEN a)

{

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

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

long i,j,k,l, n = lg(a);

GEN p1,b;

Line 746 static GEN

Line 761 static GEN

sqred2(GEN a, long no_signature)

{

GEN r,p,mun;

ulong av,av1,lim;

gpmem_t av,av1,lim;

long n,i,j,k,l,sp,sn,t;

if (typ(a)!=t_MAT) err(typeer,"sqred2");

n=lg(a); if (lg(a[1]) != n) err(mattype1,"sqred2");

n = lg(a); if (n > 1 && lg(a[1]) != n) err(mattype1,"sqred2");

av=avma; mun=negi(gun); r=new_chunk(n);

av = avma; mun = negi(gun);

for (i=1; i<n; i++) r[i]=1;

r = vecsmall_const(n-1, 1);

av1=avma; lim=stack_lim(av1,1); a=dummycopy(a);

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

n--; t=n; sp=sn=0;

a = dummycopy(a);

n--; t = n; sp = sn = 0;

while (t)

{

k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++;

Line 806 sqred2(GEN a, long no_signature)

Line 822 sqred2(GEN a, long no_signature)

if (k>n) break;

}

if (no_signature) return gerepilecopy(av,a);

if (no_signature) return gerepilecopy(av, a);

avma=av;

avma = av; a = cgetg(3,t_VEC);

a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a;

a[1] = lstoi(sp);

a[2] = lstoi(sn); return a;

}

GEN

Line 824 signat(GEN a) { return sqred2(a,0); }

Line 841 signat(GEN a) { return sqred2(a,0); }

GEN

jacobi(GEN a, long prec)

{

long de,e,e1,e2,l,n,i,j,p,q,av1,av2;

gpmem_t av1, av2;

long de,e,e1,e2,l,n,i,j,p,q;

GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1;

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

Line 840 jacobi(GEN a, long prec)

Line 858 jacobi(GEN a, long prec)

r=cgetg(l,t_MAT); ja[2]=(long)r;

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

{

r[j]=lgetg(l,t_COL);

r[j] = lgetg(l,t_COL);

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

for (i=1; i<=n; i++) coeff(r,i,j) = (long)stor(i==j, prec);

affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec)));

}

av1=avma;

Line 940 matrixqz0(GEN x,GEN p)

Line 957 matrixqz0(GEN x,GEN p)

err(flagerr,"matrixqz"); return NULL; /* not reached */

}

static int

ZV_isin(GEN x)

{

long i,l = lg(x);

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

if (typ(x[i]) != t_INT) return 0;

return 1;

}

GEN

matrixqz(GEN x, GEN p)

{

ulong av = avma, av1, lim;

gpmem_t av = avma, av1, lim;

long i,j,j1,m,n,t,nfact;

long i,j,j1,m,n,nfact;

GEN p1,p2,p3,unmodp;

GEN p1,p2,p3;

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

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

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

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

m=lg(x[1])-1;

m = lg(x[1])-1;

if (n > m) err(talker,"more rows than columns in matrixqz");

if (n==m)

{

p1=det(x);

p1 = det(x);

if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz");

avma=av; return idmat(n);

avma = av; return idmat(n);

}

p1=cgetg(n+1,t_MAT);

/* m > n */

p1 = x; x = cgetg(n+1,t_MAT);

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

{

p2=gun; p3=(GEN)x[j];

x[j] = (long)primpart((GEN)p1[j]);

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

if (!ZV_isin((GEN)p1[j])) err(talker, "not a rational matrix in matrixqz");

{

t=typ(p3[i]);

if (t != t_INT && !is_frac_t(t))

err(talker,"not a rational or integral matrix in matrixqz");

p2=ggcd(p2,(GEN)p3[i]);

}

p1[j]=ldiv(p3,p2);

}

x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un;

/* x integral */

if (gcmp0(p))

{

p1=cgetg(n+1,t_MAT); p2=gtrans(x);

p1 = gtrans(x); setlg(p1,n+1);

for (j=1; j<=n; j++) p1[j]=p2[j];

p2 = det(p1); p1[n] = p1[n+1]; p2 = ggcd(p2,det(p1));

p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1));

if (!signe(p2))

if (!signe(p3))

err(impl,"matrixqz when the first 2 dets are zero");

if (gcmp1(p3)) return gerepilecopy(av,x);

if (gcmp1(p2)) return gerepilecopy(av,x);

p3=factor(p3); p1=(GEN)p3[1]; nfact=lg(p1)-1;

p1 = (GEN)factor(p2)[1];

}

else

else p1 = _vec(p);

{

nfact = lg(p1)-1;

p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1;

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

}

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

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

{

p=(GEN)p1[i]; unmodp[1]=(long)p;

p = (GEN)p1[i];

for(;;)

{

p2=ker(gmul(unmodp,x));

p2 = FpM_ker(x, p);

if (lg(p2)==1) break;

p2=centerlift(p2); p3=gdiv(gmul(x,p2),p);

p2 = centermod(p2,p); p3 = gdiv(gmul(x,p2), p);

for (j=1; j<lg(p2); j++)

{

j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;

Line 1005 matrixqz(GEN x, GEN p)

Line 1022 matrixqz(GEN x, GEN p)

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

{

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

x=gerepilecopy(av1,x);

x = gerepilecopy(av1,x);

}

Line 1013 matrixqz(GEN x, GEN p)

Line 1030 matrixqz(GEN x, GEN p)

}

static GEN

matrixqz_aux(GEN x, long m, long n)

Z_V_mul(GEN u, GEN A)

{

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

if (gcmp1(u)) return A;

long i,j,k,in[2];

if (gcmp_1(u)) return gneg(A);

GEN p1;

if (gcmp0(u)) return zerocol(lg(A)-1);

return gmul(u,A);

}

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

static GEN

QV_lincomb(GEN u, GEN v, GEN A, GEN B)

{

if (!signe(u)) return Z_V_mul(v,B);

if (!signe(v)) return Z_V_mul(u,A);

return gadd(Z_V_mul(u,A), Z_V_mul(v,B));

}

/* cf ZV_elem */

/* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of

* A[j] and A[k] of determinant 1. */

static void

QV_elem(GEN aj, GEN ak, GEN A, long j, long k)

{

GEN p1,u,v,d, D;

if (gcmp0(ak)) { lswap(A[j],A[k]); return; }

D = mpppcm(denom(aj), denom(ak));

if (!is_pm1(D)) { aj = gmul(aj,D); ak = gmul(ak,D); }

d = bezout(aj,ak,&u,&v);

/* frequent special case (u,v) = (1,0) or (0,1) */

if (!signe(u))

{ /* ak | aj */

p1 = negi(divii(aj,ak));

A[j] = (long)QV_lincomb(gun, p1, (GEN)A[j], (GEN)A[k]);

return;

}

if (!signe(v))

{ /* aj | ak */

p1 = negi(divii(ak,aj));

A[k] = (long)QV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]);

lswap(A[j], A[k]);

return;

}

if (!is_pm1(d)) { aj = divii(aj,d); ak = divii(ak,d); }

p1 = (GEN)A[k]; aj = negi(aj);

A[k] = (long)QV_lincomb(u,v, (GEN)A[j],p1);

A[j] = (long)QV_lincomb(aj,ak, p1,(GEN)A[j]);

}

static GEN

matrixqz_aux(GEN A)

{

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

long i,j,k,n,m;

GEN a;

n = lg(A); if (n == 1) return cgetg(1,t_MAT);

m = lg(A[1]);

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

{

for(;;)

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

{

long fl=0;

GEN a = gcoeff(A,i,j);

if (gcmp0(a)) continue;

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

k = j+1; if (k == n) k = 1;

if (!gcmp0(gcoeff(x,i,j)))

/* zero a = Aij using b = Aik */

{ in[fl]=j; fl++; if (fl==2) break; }

QV_elem(a, gcoeff(A,i,k), A, j,k);

if (j>n) break;

j=(gcmp(gabs(gcoeff(x,i,in[0]),DEFAULTPREC),

gabs(gcoeff(x,i,in[1]),DEFAULTPREC)) > 0)? in[1]: in[0];

p1=gcoeff(x,i,j);

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

if (k!=j)

x[k]=lsub((GEN)x[k],gmul(ground(gdiv(gcoeff(x,i,k),p1)),(GEN)x[j]));

}

a = gcoeff(A,i,k);

j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++;

if (!gcmp0(a))

if (j<=n)

{

p1=denom(gcoeff(x,i,j));

a = denom(a);

if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]);

if (!is_pm1(a)) A[k] = lmul((GEN)A[k], a);

}

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

{

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

x=gerepilecopy(av,x);

A = gerepilecopy(av,A);

}

return hnf(x);

return m > 100? hnfall_i(A,NULL,1): hnf(A);

}

GEN

matrixqz2(GEN x)

{

long av = avma,m,n;

gpmem_t av = avma;

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

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

x = dummycopy(x);

m=lg(x[1])-1; x=dummycopy(x);

return gerepileupto(av, matrixqz_aux(x));

return gerepileupto(av, matrixqz_aux(x,m,n));

}

GEN

matrixqz3(GEN x)

{

long av=avma,av1,j,j1,k,m,n,lim;

gpmem_t av = avma, av1, lim;

long j,j1,k,m,n;

GEN c;

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

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

n = lg(x); if (n==1) return gcopy(x);

m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1);

m = lg(x[1]); x = dummycopy(x);

for (j=1; j<=n; j++) c[j]=0;

c = cgetg(n,t_VECSMALL);

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

for (j=1; j<n; j++) c[j] = 0;

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

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

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

{

j=1; while (j<=n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;

j=1; while (j<n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;

if (j<=n)

if (j==n) continue;

{

c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));

for (j1=1; j1<=n; j1++)

for (j1=1; j1<n; j1++)

if (j1!=j) x[j1]=lsub((GEN)x[j1],gmul(gcoeff(x,k,j1),(GEN)x[j]));

if (j1!=j)

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

{

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

GEN t = gcoeff(x,k,j1);

x=gerepilecopy(av1,x);

if (!gcmp0(t)) x[j1] = lsub((GEN)x[j1],gmul(t,(GEN)x[j]));

}

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

{

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

x = gerepilecopy(av1,x);

}

return gerepileupto(av, matrixqz_aux(x,m,n));

return gerepileupto(av, matrixqz_aux(x));

}

GEN

intersect(GEN x, GEN y)

{

long av,tetpil,j, lx = lg(x);

gpmem_t av,tetpil;

long j, lx = lg(x);

GEN z;

if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect");

Line 1128 mathnf0(GEN x, long flag)

Line 1193 mathnf0(GEN x, long flag)

}

static GEN

init_hnf(GEN x, GEN *denx, long *co, long *li, long *av)

init_hnf(GEN x, GEN *denx, long *co, long *li, gpmem_t *av)

{

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

*co=lg(x); if (*co==1) return NULL; /* empty matrix */

Line 1136 init_hnf(GEN x, GEN *denx, long *co, long *li, long *a

Line 1201 init_hnf(GEN x, GEN *denx, long *co, long *li, long *a

if (gcmp1(*denx)) /* no denominator */

{ *denx = NULL; return dummycopy(x); }

return gmul(x,*denx);

return Q_muli_to_int(x, *denx);

}

GEN

Line 1182 ZV_Z_mul(GEN c, GEN X)

Line 1247 ZV_Z_mul(GEN c, GEN X)

GEN

ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)

{

long av,i,lx,m;

gpmem_t av;

long i,lx,m;

GEN p1,p2,A;

if (!signe(u)) return ZV_Z_mul(v,Y);

if (!signe(v)) return ZV_Z_mul(u,X);

lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4;

if (is_pm1(v)) { gswap(u,v); gswap(X,Y); }

if (is_pm1(v)) { swap(u,v); swap(X,Y); }

if (is_pm1(u))

{

if (signe(u) > 0)

Line 1243 ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)

Line 1309 ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)

GEN

hnf_special(GEN x, long remove)

{

long av0, s,li,co,av,tetpil,i,j,k,def,ldef,lim;

gpmem_t av0,av,tetpil,lim;

long s,li,co,i,j,k,def,ldef;

GEN p1,u,v,d,denx,a,b, x2,res;

if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special");

Line 1311 hnf_special(GEN x, long remove)

Line 1378 hnf_special(GEN x, long remove)

for (i=1,j=1; j<co; j++)

if (!gcmp0((GEN)x[j]))

{

x[i] = x[j];

x2[i] = x2[j]; i++;

}

setlg(x,i);

Line 1367 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1434 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

{

GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;

GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;

long av,i,j,k,s,i1,j1,lim,zc;

gpmem_t av, lim;

long i,j,k,s,i1,j1,zc;

long co = lg(C);

long col = lg(matgen)-1;

long lnz, nlze, lig;

Line 1387 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1455 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

}

p1 = hnflll(matgen);

H = (GEN)p1[1]; /* lnz x lnz */

H = (GEN)p1[1]; /* lnz x lnz [disregarding initial 0 cols] */

H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);

U = (GEN)p1[2]; /* col x col */

/* Only keep the part above the H (above the 0s is 0 since the dep rows

* are dependant from the ones in matgen) */

Line 1397 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1466 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */

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

Cnew = cgetg(co,t_MAT);

Cnew = cgetg(co, typ(C));

setlg(C, col+1); p1 = gmul(C,U);

for (j=1; j<=col; j++) Cnew[j] = p1[j];

for ( ; j<co ; j++) Cnew[j] = C[j];

Line 1406 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1475 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

/* Clean up B using new H */

for (s=0,i=lnz; i; i--)

{

GEN h = gcoeff(H,i,i);

GEN Di = (GEN)dep[i], Hi = (GEN)H[i];

GEN h = (GEN)Hi[i]; /* H[i,i] */

if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; }

for (j=col+1; j<co; j++)

{

GEN z = (GEN)B[j-col];

p1 = (GEN)z[i+nlze]; if (h) p1 = gdivent(p1,h);

p1 = (GEN)z[i+nlze]; if (!signe(p1)) continue;

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

z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(dep,k,i)));

if (h) p1 = gdivent(p1,h);

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

for (k=1; k<=nlze; k++) z[k] = lsubii((GEN)z[k], mulii(p1, (GEN)Di[k]));

z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(H,k-nlze,i)));

for ( ; k<=lig; k++) z[k] = lsubii((GEN)z[k], mulii(p1, (GEN)Hi[k-nlze]));

Cnew[j] = lsub((GEN)Cnew[j], gmul(p1, (GEN)Cnew[i+zc]));

}

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

{

GEN *gptr[2]; gptr[0]=&Cnew; gptr[1]=&B;

if(DEBUGMEM>1) err(warnmem,"hnffinal, i = %ld",i);

gerepilemany(av,gptr,2);

gerepileall(av, 2, &Cnew, &B);

}

p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;

Line 1452 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1521 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

GEN z = (GEN)H[j];

if (diagH1[j])

{ /* hit exactly s times */

i1++; p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;

i1++; C[i1+col] = Cnew[j+zc];

C[i1+col] = Cnew[j+zc];

p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;

for (i=1; i<=nlze; i++) p1[i] = coeff(dep,i,j);

p1 += nlze;

}

else

{

j1++; p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;

j1++; C[j1+zc] = Cnew[j+zc];

C[j1+zc] = Cnew[j+zc];

p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;

if (nlze) depnew[j1] = dep[j];

depnew[j1] = dep[j];

}

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

if (!diagH1[i]) p1[k++] = z[i];

Line 1484 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

Line 1553 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN*

fprintferr("mit = %Z\nB = %Z\nC = %Z\n", Hnew, Bnew, C);

}

if (nlze) *ptdep = depnew;

*ptdep = depnew;

*ptC = C;

*ptB = Bnew; return Hnew;

}

/* for debugging */

static void

p_mat(long **mat, long *perm, long k0)

p_mat(long **mat, GEN perm, long k)

{

long av=avma, i,j;

gpmem_t av = avma;

GEN p1, matj, matgen;

perm = vecextract_i(perm, k+1, lg(perm)-1);

long co = lg(mat);

long li = lg(perm);

fprintferr("Permutation: %Z\n",perm);

matgen = cgetg(co,t_MAT);

if (DEBUGLEVEL > 6)

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

fprintferr("matgen = %Z\n", small_to_mat( rowextract_p((GEN)mat, perm) ));

{

avma = av;

p1 = cgetg(li-k0,t_COL); matgen[j]=(long)p1;

p1 -= k0; matj = mat[j];

for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);

}

if (DEBUGLEVEL > 6) fprintferr("matgen = %Z\n",matgen);

avma=av;

}

static GEN

Line 1524 col_dup(long n, GEN col)

Line 1584 col_dup(long n, GEN col)

** C = r x (co-1) matrix of GEN

** perm= permutation vector (length li-1), indexing the rows of mat: easier

** to maintain perm than to copy rows. For columns we can do it directly

** using e.g. swap(mat[i], mat[j])

** using e.g. lswap(mat[i], mat[j])

** k0 = integer. The k0 first lines of mat are dense, the others are sparse.

** Also if k0=0, mat is modified in place [from mathnfspec], otherwise

** it is left alone [from buchall]

Line 1536 col_dup(long n, GEN col)

Line 1596 col_dup(long n, GEN col)

GEN

hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)

{

long av=avma,av2,*p,i,j,k,lk0,col,lig,*matj, **mat;

gpmem_t av=avma,av2,lim;

long n,s,t,lim,nlze,lnz,nr;

long *p,i,j,k,lk0,col,lig,*matj, **mat;

long n,s,t,nlze,lnz,nr;

GEN p1,p2,matb,matbnew,vmax,matt,T,extramat;

GEN B,H,dep,permpro;

GEN *gptr[4];

long co = lg(mat0);

long li = lg(perm); /* = lg(mat0[1]) */

int updateT = 1;

Line 1579 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

Line 1639 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

switch( count(mat,perm[i],col,&n) )

{

case 0: /* move zero lines between k0+1 and lk0 */

lk0++; swap(perm[i], perm[lk0]);

lk0++; lswap(perm[i], perm[lk0]);

i=lig; continue;

case 1: /* move trivial generator between lig+1 and li */

swap(perm[i], perm[lig]);

lswap(perm[i], perm[lig]);

swap(T[n], T[col]);

lswap(T[n], T[col]);

gswap(mat[n], mat[col]); p = mat[col];

swap(mat[n], mat[col]); p = mat[col];

if (p[perm[lig]] < 0) /* = -1 */

{ /* convert relation -g = 0 to g = 0 */

for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];

Line 1607 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

Line 1667 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

}

#define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}

#if 0 /* TODO: check, and put back in */

/* Get rid of all lines containing only 0 and � 1, keeping track of column

* operations in T. Leave the rows 1..lk0 alone [up to k0, coeff

* explosion, between k0+1 and lk0, row is 0]

* explosion, between k0+1 and lk0, row is 0] */

s = 0;

while (lig > lk0 && s < (HIGHBIT>>1))

{

for (i=lig; i>lk0; i--)

if (count(mat,perm[i],col,&n) >= 0) break;

if (count(mat,perm[i],col,&n) > 0) break;

if (i == lk0) break;

/* only 0, �1 entries, at least 2 of them non-zero */

swap(perm[i], perm[lig]);

lswap(perm[i], perm[lig]);

swap(T[n], T[col]); p1 = (GEN)T[col];

swap(mat[n], mat[col]); p = mat[col];

gswap(mat[n], mat[col]); p = mat[col];

lswap(T[n], T[col]); p1 = (GEN)T[col];

if (p[perm[lig]] < 0)

{

for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];

T[col] = lneg(p1);

p1 = gneg(p1); T[col] = (long)p1;

}

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

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

{

matj = mat[j];

if (! (t = matj[perm[lig]]) ) continue;

Line 1653 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

Line 1710 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

T = gerepilecopy(av2, T);

}

#endif

/* As above with lines containing a �1 (no other assumption).

* Stop when single precision becomes dangerous */

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

Line 1668 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

Line 1724 hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, G

if ( (n = count2(mat,perm[i],col)) ) break;

if (i == lk0) break;

swap(perm[i], perm[lig]);

lswap(vmax[n], vmax[col]);

swap(vmax[n], vmax[col]);

lswap(perm[i], perm[lig]);

gswap(mat[n], mat[col]); p = mat[col];

swap(mat[n], mat[col]); p = mat[col];

swap(T[n], T[col]); p1 = (GEN)T[col];

lswap(T[n], T[col]); p1 = (GEN)T[col];

if (p[perm[lig]] < 0)

{

for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];

Line 1709 END2: /* clean up mat: remove everything to the right

Line 1765 END2: /* clean up mat: remove everything to the right

if (DEBUGLEVEL>5)

{

fprintferr(" after phase2:\n");

p_mat(mat,perm,k0);

p_mat(mat,perm,lk0);

}

for (i=li-2; i>lig; i--)

{

Line 1748 END2: /* clean up mat: remove everything to the right

Line 1804 END2: /* clean up mat: remove everything to the right

{

if(DEBUGMEM>1) err(warnmem,"hnfspec[3], i = %ld", i);

for (j=1; j<co; j++) setlg(matb[j], i0+1); /* bottom can be forgotten */

gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);

gerepileall(av2, 2, &T, &matb);

}

gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);

gerepileall(av2, 2, &T, &matb);

if (DEBUGLEVEL>5)

{

fprintferr(" matb cleaned up (using Id block)\n");

Line 1829 END2: /* clean up mat: remove everything to the right

Line 1885 END2: /* clean up mat: remove everything to the right

*ptdep = dep;

*ptB = B;

H = hnffinal(matbnew,perm,ptdep,ptB,ptC);

gptr[0]=ptC;

gerepileall(av, 4, ptC, ptdep, ptB, &H);

gptr[1]=ptdep;

gptr[2]=ptB;

gptr[3]=&H; gerepilemany(av,gptr,4);

if (DEBUGLEVEL)

msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1);

return H;

Line 1891 GEN

Line 1944 GEN

hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */

GEN extramat,GEN extraC)

{

GEN p1,p2,p3,matb,extratop,Cnew,permpro;

GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;

GEN B=*ptB, C=*ptC, dep=*ptdep, *gptr[4];

gpmem_t av = avma;

long av = avma, i,j,lextra,colnew;

long i;

long li = lg(perm);

long co = lg(C);

long lB = lg(B);

long lig = li - lB;

long col = co - lB;

long lH = lg(H)-1;

long lB = lg(B)-1;

long li = lg(perm)-1, lig = li - lB;

long co = lg(C)-1, col = co - lB;

long nlze = lH? lg(dep[1])-1: lg(B[1])-1;

if (DEBUGLEVEL>5)

Line 1914 hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC

Line 1965 hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC

* [--------|-----] lig

* [ 0 | Id ]

* [ | ] li */

lextra = lg(extramat)-1;

extratop = rowextract_i(extramat, 1, lig);

extratop = cgetg(lextra+1,t_MAT); /* [1..lig] part (top) */

if (li != lig)

p2 = cgetg(lextra+1,t_MAT); /* bottom */

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

{

GEN z = (GEN)extramat[j];

p1=cgetg(lig+1,t_COL); extratop[j] = (long)p1;

for (i=1; i<=lig; i++) p1[i] = z[i];

p1=cgetg(lB,t_COL); p2[j] = (long)p1;

p1 -= lig;

for ( ; i<li; i++) p1[i] = z[i];

}

if (li-1 != lig)

{ /* zero out bottom part, using the Id block */

GEN A = cgetg(lB,t_MAT);

GEN A = vecextract_i(C, col+1, co);

for (j=1; j<lB; j++) A[j] = C[j+col];

GEN c = rowextract_i(extramat, lig+1, li);

extraC = gsub(extraC, gmul(A,p2));

extraC = gsub(extraC, gmul(A, c));

extratop = gsub(extratop,gmul(B,p2));

extratop = gsub(extratop, gmul(B, c));

}

colnew = lH + lextra;

extramat = concatsp(extratop, vconcat(dep, H));

extramat = cgetg(colnew+1,t_MAT);

Cnew = concatsp(extraC, vecextract_i(C, col-lH+1, co));

Cnew = cgetg(lB+colnew,t_MAT);

if (DEBUGLEVEL>5) fprintferr(" 1st phase done\n");

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

{

extramat[j] = extratop[j];

Cnew[j] = extraC[j];

}

for ( ; j<=colnew; j++)

{

p1 = cgetg(lig+1,t_COL); extramat[j] = (long)p1;

p2 = (GEN)dep[j-lextra]; for (i=1; i<=nlze; i++) p1[i] = p2[i];

p2 = (GEN) H[j-lextra]; for ( ; i<=lig ; i++) p1[i] = p2[i-nlze];

}

for (j=lextra+1; j<lB+colnew; j++)

Cnew[j] = C[j-lextra+col-lH];

if (DEBUGLEVEL>5)

{

fprintferr(" 1st phase done\n");

if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);

}

permpro = imagecomplspec(extramat, &nlze);

p1 = new_chunk(lig+1);

extramat = rowextract_p(extramat, permpro);

for (i=1; i<=lig; i++) p1[i] = perm[permpro[i]];

*ptB = rowextract_p(B, permpro);

for (i=1; i<=lig; i++) perm[i] = p1[i];

/* perm o= permpro */

permpro = vecextract_p(perm, permpro);

for (i=1; i<=lig; i++) perm[i] = permpro[i];

matb = cgetg(colnew+1,t_MAT);

*ptdep = rowextract_i(extramat, 1, nlze);

dep = cgetg(colnew+1,t_MAT);

matb = rowextract_i(extramat, nlze+1, lig);

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

{

GEN z = (GEN)extramat[j];

p1=cgetg(nlze+1,t_COL); dep[j]=(long)p1;

p2=cgetg(lig+1-nlze,t_COL); matb[j]=(long)p2;

p2 -= nlze;

for (i=1; i<=nlze; i++) p1[i] = z[permpro[i]];

for ( ; i<= lig; i++) p2[i] = z[permpro[i]];

}

p3 = cgetg(lB,t_MAT);

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

{

p2 = (GEN)B[j];

p1 = cgetg(lig+1,t_COL); p3[j] = (long)p1;

for (i=1; i<=lig; i++) p1[i] = p2[permpro[i]];

}

B = p3;

if (DEBUGLEVEL>5) fprintferr(" 2nd phase done\n");

*ptdep = dep;

*ptB = B;

H = hnffinal(matb,perm,ptdep,ptB,&Cnew);

p1 = cgetg(co+lextra,t_MAT);

*ptC = concatsp(vecextract_i(C, 1, col-lH), Cnew);

for (j=1; j <= col-lH; j++) p1[j] = C[j];

C = Cnew - (col-lH);

for ( ; j < co+lextra; j++) p1[j] = C[j];

gptr[0]=ptC; *ptC=p1;

gptr[1]=ptdep;

gptr[2]=ptB;

gptr[3]=&H; gerepilemany(av,gptr,4);

if (DEBUGLEVEL)

{

if (DEBUGLEVEL>7) fprintferr("mit = %Z\nC = %Z\n",H,*ptC);

if (DEBUGLEVEL>7) fprintferr("H = %Z\nC = %Z\n",H,*ptC);

msgtimer("hnfadd (%ld)",lextra);

msgtimer("hnfadd (%ld)", lg(extramat)-1);

}

gerepileall(av, 4, ptC, ptdep, ptB, &H);

return H;

}

Line 2014 ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)

Line 2013 ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)

{

GEN p1,u,v,d;

if (!signe(ak)) { swap(A[j],A[k]); if (U) swap(U[j],U[k]); return; }

if (!signe(ak)) { lswap(A[j],A[k]); if (U) lswap(U[j],U[k]); return; }

d = bezout(aj,ak,&u,&v);

/* frequent special case (u,v) = (1,0) or (0,1) */

if (!signe(u))

Line 2029 ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)

Line 2028 ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)

{ /* aj | ak */

p1 = negi(divii(ak,aj));

A[k] = (long)ZV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]);

swap(A[j], A[k]);

lswap(A[j], A[k]);

if (U) {

U[k] = (long)ZV_lincomb(gun, p1, (GEN)U[k], (GEN)U[j]);

swap(U[j], U[k]);

lswap(U[j], U[k]);

}

return;

}

Line 2077 ZM_reduce(GEN A, GEN U, long i, long j0)

Line 2076 ZM_reduce(GEN A, GEN U, long i, long j0)

GEN

hnf0(GEN A, int remove)

{

long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;

gpmem_t av0 = avma, av, lim;

long s,li,co,i,j,k,def,ldef;

GEN denx,a,p1;

if (typ(A) == t_VEC) return hnf_special(A,remove);

Line 2100 hnf0(GEN A, int remove)

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

{

if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i);

tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));

A = gerepilecopy(av, A);

}

p1 = gcoeff(A,i,def); s = signe(p1);

Line 2115 hnf0(GEN A, int remove)

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

{

if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i);

tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));

A = gerepilecopy(av, A);

}

if (remove)

Line 2124 hnf0(GEN A, int remove)

if (!gcmp0((GEN)A[j])) A[i++] = A[j];

setlg(A,i);

}

tetpil=avma;

A = denx? gdiv(A,denx): ZM_copy(A);

return gerepile(av0,tetpil,A);

return gerepileupto(av0, A);

}

GEN

Line 2137 hnf(GEN x) { return hnf0(x,1); }

Line 2136 hnf(GEN x) { return hnf0(x,1); }

static GEN

allhnfmod(GEN x,GEN dm,int flag)

{

ulong av,tetpil,lim;

gpmem_t av, lim;

long li,co,i,j,k,def,ldef,ldm;

GEN a,b,w,p1,p2,d,u,v;

Line 2184 allhnfmod(GEN x,GEN dm,int flag)

Line 2183 allhnfmod(GEN x,GEN dm,int flag)

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

{

if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i);

tetpil=avma; x=gerepile(av,tetpil,ZM_copy(x));

x = gerepilecopy(av, x);

}

w = cgetg(li,t_MAT); b = dm;

Line 2214 allhnfmod(GEN x,GEN dm,int flag)

Line 2213 allhnfmod(GEN x,GEN dm,int flag)

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

{

if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i);

tetpil=avma; w=gerepile(av,tetpil,ZM_copy(w));

gerepileall(av, 2, &w, &dm); diag = gcoeff(w,i,i);

diag = gcoeff(w,i,i);

}

tetpil=avma; return gerepile(av,tetpil,ZM_copy(w));

return gerepilecopy(av, w);

}

GEN

Line 2238 hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }

Line 2236 hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }

static GEN

mynegi(GEN x)

{

static long mun[]={evaltyp(t_INT)|m_evallg(3),evalsigne(-1)|evallgefint(3),1};

static long mun[]={evaltyp(t_INT)|_evallg(3),evalsigne(-1)|evallgefint(3),1};

long s = signe(x);

if (!s) return gzero;

Line 2247 mynegi(GEN x)

Line 2245 mynegi(GEN x)

setsigne(x,-s); return x;

}

/* We assume y > 0 */

static GEN

divnearest(GEN x, GEN y)

{

GEN r, q = dvmdii(x,y,&r);

long av = avma, s=signe(r), t;

if (!s) { cgiv(r); return q; }

if (s<0) r = mynegi(r);

y = shifti(y,-1); t = cmpii(r,y);

avma=av; cgiv(r);

if (t < 0 || (t == 0 && s > 0)) return q;

return addsi(s,q);

}

static void

Minus(long j, GEN **lambda)

{

Line 2313 findi(GEN M)

Line 2295 findi(GEN M)

return 0;

}

static long

findi_normalize(GEN Aj, GEN B, long j, GEN **lambda)

{

long r = findi(Aj);

if (r && signe(Aj[r]) < 0)

{

neg_col(Aj); if (B) neg_col((GEN)B[j]);

Minus(j,lambda);

}

return r;

}

static void

reduce2(GEN A, GEN B, long k, long j, long *row, GEN **lambda, GEN *D)

{

GEN q;

long i, row0, row1;

row0 = findi((GEN)A[j]);

row[0] = row0 = findi_normalize((GEN)A[j], B,j,lambda);

if (row0 && signe(gcoeff(A,row0,j)) < 0)

row[1] = row1 = findi_normalize((GEN)A[k], B,k,lambda);

{

neg_col((GEN)A[j]);

if (B) neg_col((GEN)B[j]);

Minus(j,lambda);

}

row1 = findi((GEN)A[k]);

if (row1 && signe(gcoeff(A,row1,k)) < 0)

{

neg_col((GEN)A[k]);

if (B) neg_col((GEN)B[k]);

Minus(k,lambda);

}

row[0]=row0; row[1]=row1;

if (row0)

q = truedvmdii(gcoeff(A,row0,k), gcoeff(A,row0,j), NULL);

else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)

q = divnearest(lambda[k][j], D[j]);

q = diviiround(lambda[k][j], D[j]);

else

return;

Line 2371 reduce2(GEN A, GEN B, long k, long j, long *row, GEN *

Line 2350 reduce2(GEN A, GEN B, long k, long j, long *row, GEN *

}

#define SWAP(x,y) {GEN _t=y; y=x; x=_t;}

static void

hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)

{

GEN t,p1,p2;

long i,j,n = lg(A);

swap(A[k],A[k-1]);

lswap(A[k],A[k-1]);

if (B) swap(B[k],B[k-1]);

if (B) lswap(B[k],B[k-1]);

for (j=k-2; j; j--) SWAP(lambda[k-1][j], lambda[k][j]);

for (j=k-2; j; j--) swap(lambda[k-1][j], lambda[k][j]);

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

{

p1 = mulii(lambda[i][k-1], D[k]);

Line 2399 hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)

Line 2376 hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)

D[k-1] = divii(addii(p1,p2), D[k-1]);

}

/* A[k] = 0, A[nz] != 0, k > nz, A[j] = 0 for all j < nz.

/* reverse row order in matrix A */

* "Deep insert" A[k] at nz */

static void

trivswap(GEN *A, long nz, long k, GEN **lambda, GEN *D)

{

GEN p1;

long i,j,n = lg(A);

p1 = A[nz]; /* cycle A */

for (j = nz; j < k; j++) SWAP(A[j+1], p1);

A[nz] = p1; /* = [0...0] */

p1 = D[nz]; /* cycle D */

for (j = nz; j < k; j++) SWAP(D[j+1], p1);

D[nz] = gun;

for (j=k-1; j>=nz; j--) /* cycle lambda */

for (i=k-1; i>=nz; i--) lambda[i+1][j+1] = lambda[i][j];

for (j=n-1; j>k; j--)

for (i=k-1; i>=nz; i--)

{

lambda[i+1][j] = lambda[i][j];

lambda[j][i+1] = lambda[j][i];

}

for (i=1; i<n; i++) lambda[nz][i] = lambda[i][nz] = gzero;

}

static GEN

fix_rows(GEN A)

{

Line 2443 fix_rows(GEN A)

Line 2394 fix_rows(GEN A)

}

GEN

hnflll_i(GEN A, GEN *ptB)

hnflll_i(GEN A, GEN *ptB, int remove)

{

ulong av = avma, lim = stack_lim(av,3);

gpmem_t av = avma, lim = stack_lim(av,3);

long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */

long row[2], do_swap,i,n,k;

long nzcol = 1; /* index of 1st (maybe) non-0 col [only updated when !B] */

GEN z,B, **lambda, *D;

GEN *gptr[4];

if (typ(A) != t_MAT) err(typeer,"hnflll");

n = lg(A);

A = ZM_copy(fix_rows(A));

B = ptB? idmat(n-1): NULL;

D = (GEN*) cgetg(n+1, t_VEC); D++; /* hack: need a "sentinel" D[0] */

D = (GEN*)cgetg(n+1,t_VEC); lambda = (GEN**) cgetg(n,t_MAT);

if (n == 2) /* handle trivial case: return negative diag coeff otherwise */

D++;

{

for (i=0; i<n; i++) D[i] = gun;

i = findi((GEN)A[1]);

for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);

if (i && signe(gcoeff(A,i,1)) < 0)

k = 2;

{

neg_col((GEN)A[1]);

if (B) neg_col((GEN)B[1]);

}

lambda = (GEN**) cgetg(n,t_MAT);

for (i=1; i<n; i++) { D[i] = gun; lambda[i] = (GEN*)zerocol(n-1); }

D[0] = gun; k = 2;

while (k < n)

{

reduce2(A,B,k,k-1,row,lambda,D);

if (!B)

if (row[0])

{ /* not interested in B. Eliminate 0 cols */

do_swap = (!row[1] || row[0] <= row[1]);

if (!row[0] || !row[1])

else if (!row[1])

{

{ /* row[0] == row[1] == 0 */

while (!findi((GEN)A[nzcol]) && nzcol < n) nzcol++;

gpmem_t av1 = avma;

/* A[nzcol] != 0, A[i] = 0 for i < nzcol */

z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));

if (!row[0] && k-1 > nzcol)

do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);

{trivswap((GEN*)A,nzcol,k-1, lambda,D); nzcol++;}

avma = av1;

if (!row[1] && k > nzcol)

{trivswap((GEN*)A,nzcol,k , lambda,D); nzcol++;}

if (k <= nzcol) k = nzcol+1;

continue;

}

do_swap = (row[0] && row[0] <= row[1]);

}

else

{

do_swap = 0;

if (row[0])

do_swap = (!row[1] || row[0] <= row[1]);

else if (row[1])

do_swap = 0;

else

{ /* row[0] == row[1] == 0 */

ulong av1 = avma;

z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));

do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);

avma = av1;

}

if (do_swap)

{

hnfswap(A,B,k,lambda,D);

if (k > nzcol+1) k--;

if (k > 2) k--;

}

else

{

for (i=k-2; i>=nzcol; i--) reduce2(A,B,k,i,row,lambda,D);

for (i=k-2; i; i--)

{

reduce2(A,B,k,i,row,lambda,D);

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

{

GEN b = (GEN)(D-1);

if (DEBUGMEM) err(warnmem,"hnflll (reducing), i = %ld",i);

gerepileall(av, B? 4: 3, &A, (GEN*)&lambda, &b, &B);

D = (GEN*)(b+1);

}

k++;

}

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

{

GEN a = (GEN)lambda, b = (GEN)(D-1); /* gcc -Wall */

GEN b = (GEN)(D-1);

gptr[0]=&A; gptr[1]=&a; gptr[2]=&b; gptr[3]=&B;

if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n-1);

if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n);

gerepileall(av, B? 4: 3, &A, (GEN*)&lambda, &b, &B);

gerepilemany(av,gptr,B? 4: 3); lambda = (GEN**)a; D = (GEN*)(b+1);

D = (GEN*)(b+1);

}

for (i=nzcol; i<n; i++)

/* handle trivial case: return negative diag coeff otherwise */

if (findi((GEN)A[i])) break;

if (n == 2) (void)findi_normalize((GEN)A[1], B,1,lambda);

i--; A += i; A[0] = evaltyp(t_MAT)|evallg(n-i);

A = fix_rows(A);

gptr[0] = &A; gptr[1] = &B;

if (remove)

gerepilemany(av, gptr, B? 2: 1);

{

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

if (findi((GEN)A[i])) break;

i--;

A += i; A[0] = evaltyp(t_MAT) | evallg(n-i);

}

gerepileall(av, B? 2: 1, &A, &B);

if (B) *ptB = B;

return A;

}

Line 2533 GEN

Line 2471 GEN

hnflll(GEN A)

{

GEN B, z = cgetg(3, t_VEC);

z[1] = (long)hnflll_i(A, &B);

z[1] = (long)hnflll_i(A, &B, 0);

z[2] = (long)B; return z;

}

Line 2546 reduce1(GEN A, GEN B, long k, long j, GEN **lambda, GE

Line 2484 reduce1(GEN A, GEN B, long k, long j, GEN **lambda, GE

long i;

if (signe(A[j]))

q = divnearest((GEN)A[k],(GEN)A[j]);

q = diviiround((GEN)A[k],(GEN)A[j]);

else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)

q = divnearest(lambda[k][j], D[j]);

q = diviiround(lambda[k][j], D[j]);

else

return;

if (! gcmp0(q))

if (signe(q))

{

GEN *Lk = lambda[k], *Lj = lambda[j];

q = mynegi(q);

A[k] = laddii((GEN)A[k], mulii(q,(GEN)A[j]));

elt_col((GEN)B[k],(GEN)B[j],q);

lambda[k][j] = addii(lambda[k][j], mulii(q,D[j]));

Lk[j] = addii(Lk[j], mulii(q,D[j]));

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

if (signe(lambda[j][i]))

if (signe(Lj[i])) Lk[i] = addii(Lk[i], mulii(q,Lj[i]));

lambda[k][i] = addii(lambda[k][i], mulii(q,lambda[j][i]));

}

Line 2568 GEN

Line 2506 GEN

extendedgcd(GEN A)

{

long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */

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

gpmem_t av = avma;

GEN p1,p2,B, **lambda, *D;

long do_swap,i,n,k;

GEN z,B, **lambda, *D;

n = lg(A); B = idmat(n-1); A = ZM_copy(A);

n = lg(A);

D = (GEN*) cgeti(n); lambda = (GEN**) cgetg(n,t_MAT);

A = dummycopy(A);

B = idmat(n-1);

D = (GEN*)new_chunk(n); lambda = (GEN**) cgetg(n,t_MAT);

for (i=0; i<n; i++) D[i] = gun;

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

for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);

{

lambda[i] = (GEN*) cgetg(n,t_COL);

for(j=1;j<n;j++) lambda[i][j] = gzero;

}

k = 2;

while (k < n)

{

Line 2586 extendedgcd(GEN A)

Line 2523 extendedgcd(GEN A)

if (signe(A[k-1])) do_swap = 1;

else if (!signe(A[k]))

{

long av1=avma;

gpmem_t av1 = avma;

p1 = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));

z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));

p1 = mulsi(n1,p1);

do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);

p2 = mulsi(m1,sqri(D[k-1]));

avma = av1;

do_swap = (cmpii(p1,p2) < 0);

avma=av1;

}

else do_swap = 0;

if (do_swap)

{

hnfswap(A,B,k,lambda,D);

if (k>2) k--;

if (k > 2) k--;

}

else

{

Line 2606 extendedgcd(GEN A)

Line 2541 extendedgcd(GEN A)

k++;

}

if (signe(A[n-1])<0)

if (signe(A[n-1]) < 0)

{

A[n-1] = (long)mynegi((GEN)A[n-1]);

neg_col((GEN)B[n-1]);

}

tetpil = avma;

z = cgetg(3,t_VEC);

p1 = cgetg(3,t_VEC);

z[1] = A[n-1];

p1[1]=lcopy((GEN)A[n-1]);

z[2] = (long)B;

p1[2]=lcopy(B);

return gerepilecopy(av, z);

return gerepile(av,tetpil,p1);

}

/* HNF with permutation. TODO: obsolete, replace with hnflll. */

Line 2623 GEN

Line 2557 GEN

hnfperm(GEN A)

{

GEN U,c,l,perm,d,u,p,q,y,b;

long r,t,i,j,j1,k,m,n,av=avma,av1,tetpil,lim;

gpmem_t av=avma,av1,tetpil,lim;

long r,t,i,j,j1,k,m,n;

if (typ(A) != t_MAT) err(typeer,"hnfperm");

n = lg(A)-1;

Line 2696 hnfperm(GEN A)

Line 2631 hnfperm(GEN A)

}

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

{

GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;

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

gerepilemany(av1,gptr,2);

gerepileall(av1, 2, &A, &U);

}

if (r < m)

Line 2735 GEN

Line 2669 GEN

hnfall_i(GEN A, GEN *ptB, long remove)

{

GEN B,c,h,x,p,a;

long m,n,r,i,j,k,li,av=avma,av1,tetpil,lim;

gpmem_t av = avma, av1, lim;

GEN *gptr[2];

long m,n,r,i,j,k,li;

if (typ(A)!=t_MAT) err(typeer,"hnfall");

n=lg(A)-1;

Line 2766 hnfall_i(GEN A, GEN *ptB, long remove)

Line 2700 hnfall_i(GEN A, GEN *ptB, long remove)

ZM_reduce(A,B, i,k);

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

{

tetpil = avma;

A = ZM_copy(A); gptr[0]=&A;

if (B) { B = ZM_copy(B); gptr[1]=&B; }

if (DEBUGMEM>1) err(warnmem,"hnfall[1], li = %ld", li);

gerepilemanysp(av1,tetpil,gptr,B? 2: 1);

gerepileall(av1, B? 2: 1, &A, &B);

}

x = gcoeff(A,li,j); if (signe(x)) break;

Line 2780 hnfall_i(GEN A, GEN *ptB, long remove)

Line 2711 hnfall_i(GEN A, GEN *ptB, long remove)

r--;

if (j < r) /* A[j] != 0 */

{

swap(A[j], A[r]);

lswap(A[j], A[r]);

if (B) swap(B[j], B[r]);

if (B) lswap(B[j], B[r]);

h[j]=h[r]; h[r]=li; c[li]=r;

}

p = gcoeff(A,li,r);

Line 2794 hnfall_i(GEN A, GEN *ptB, long remove)

Line 2725 hnfall_i(GEN A, GEN *ptB, long remove)

ZM_reduce(A,B, li,r);

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

{

GEN *gptr[2]; gptr[0]=&A; gptr[1]=&B;

if (DEBUGMEM>1) err(warnmem,"hnfall[2], li = %ld", li);

gerepilemany(av1,gptr,B? 2: 1);

gerepileall(av1, B? 2: 1, &A, &B);

}

if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");

Line 2812 hnfall_i(GEN A, GEN *ptB, long remove)

Line 2742 hnfall_i(GEN A, GEN *ptB, long remove)

ZM_reduce(A,B, i,k);

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

{

tetpil = avma;

A = ZM_copy(A); gptr[0]=&A;

if (B) { B = ZM_copy(B); gptr[1]=&B; }

if (DEBUGMEM>1) err(warnmem,"hnfall[3], j = %ld", j);

gerepilemanysp(av1,tetpil,gptr, B? 2: 1);

gerepileall(av1, B? 2: 1, &A, &B);

}

if (DEBUGLEVEL>5) fprintferr("\n");

/* remove the first r columns */

if (remove) { A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1); }

gptr[0] = &A; gptr[1] = &B;

gerepileall(av, B? 2: 1, &A, &B);

gerepilemany(av, gptr, B? 2: 1);

if (B) *ptB = B;

return A;

}

Line 2896 trivsmith(long all)

Line 2822 trivsmith(long all)

z[3]=lgetg(1,t_MAT); return z;

}

/* Return the smith normal form d of matrix x. If all != 0 return [d,u,v],

static void

* where d = u.x.v

snf_pile(gpmem_t av, GEN *x, GEN *U, GEN *V)

static GEN

smithall(GEN x, long all)

{

long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim;

GEN *gptr[3];

GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys;

int c = 1; gptr[0]=x;

if (*U) gptr[c++] = U;

if (*V) gptr[c++] = V;

gerepilemany(av,gptr,c);

}

/* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V

* to that D = UXV */

GEN

smithall(GEN x, GEN *ptU, GEN *ptV)

{

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

long i,j,k,l,c,n,s1,s2;

GEN p1,p2,p3,p4,b,u,v,d,U,V,mun,mdet,ys;

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

if (DEBUGLEVEL>=9) outerr(x);

if (DEBUGLEVEL>8) outerr(x);

n=lg(x)-1;

n = lg(x)-1;

if (!n) return trivsmith(all);

if (!n) {

if (ptU) *ptU = cgetg(1,t_MAT);

if (ptV) *ptV = cgetg(1,t_MAT);

return cgetg(1,t_MAT);

}

if (lg(x[1]) != n+1) err(mattype1,"smithall");

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

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

if (typ(coeff(x,i,j)) != t_INT)

err(talker,"non integral matrix in smithall");

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

U = ptU? gun: NULL;

x = dummycopy(x); mdet = detint(x);

V = ptV? gun: NULL;

if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }

x = dummycopy(x);

if (ishnfall(x))

{

mdet = dethnf_i(x);

if (V) V = idmat(n);

}

else

{

mdet = detint(x);

if (signe(mdet))

{

p1=hnfmod(x,mdet);

p1 = hnfmod(x,mdet);

if (all) { ml=idmat(n); mr=gauss(x,p1); }

if (V) V = gauss(x,p1);

}

else

{

if (all)

p1 = hnfall_i(x, ptV, 0);

{

if (ptV) V = *ptV;

p1 = hnfall0(x,0);

ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1];

}

else

p1 = hnf0(x,0);

}

x = p1;

}

if (U) U = idmat(n);

p1=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) p1[i]=lnegi(gcoeff(x,i,i));

p2=sindexsort(p1); ys=cgetg(n+1,t_MAT);

p2=sindexsort(p1);

ys = cgetg(n+1,t_MAT);

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

{

p1=cgetg(n+1,t_COL); ys[j]=(long)p1;

for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]);

}

x = ys;

if (all)

if (U)

{

p3=cgetg(n+1,t_MAT); p4=cgetg(n+1,t_MAT);

p1 = cgetg(n+1,t_MAT);

for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }

for (j=1; j<=n; j++) p1[j]=U[p2[j]];

ml=p3; mr=p4;

U = p1;

}

if (V)

{

p1 = cgetg(n+1,t_MAT);

for (j=1; j<=n; j++) p1[j]=V[p2[j]];

V = p1;

}

if (signe(mdet))

{

p1 = hnfmod(x,mdet);

if (all) mr=gmul(mr,gauss(x,p1));

if (V) V = gmul(V, gauss(x,p1));

}

else

{

if (all)

p1 = hnfall_i(x,ptV,0);

{

if (ptV) V = gmul(V, *ptV);

p1 = hnfall0(x,0);

mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1];

}

else

p1 = hnf0(x,0);

}

x=p1; mun = negi(gun);

x = p1; mun = negi(gun);

if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();}

if (DEBUGLEVEL>7) fprintferr("starting SNF loop");

for (i=n; i>=2; i--)

for (i=n; i>1; i--)

{

if (DEBUGLEVEL>7) {fprintferr("\ni = %ld: ",i);flusherr();}

if (DEBUGLEVEL>7) fprintferr("\ni = %ld: ",i);

for(;;)

{

c = 0;

for (j=i-1; j>=1; j--)

{

p1=gcoeff(x,i,j); s1 = signe(p1);

if (s1)

if (!s1) continue;

{

p2=gcoeff(x,i,i);

if (!absi_cmp(p1,p2))

{

s2=signe(p2);

if (s1 == s2) { d=p1; u=gun; p4=gun; }

else

{

s2=signe(p2);

if (s2>0) { u = gun; p4 = mun; }

if (s1 == s2) { d=p1; u=gun; p4=gun; }

else { u = mun; p4 = gun; }

else

d=(s1>0)? p1: absi(p1);

{

if (s2>0) { u = gun; p4 = mun; }

else { u = mun; p4 = gun; }

d=(s1>0)? p1: absi(p1);

}

v = gzero; p3 = u;

}

else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }

v = gzero; p3 = u;

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

}

{

else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }

b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));

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

coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),

{

mulii(p4,gcoeff(x,k,i)));

b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));

coeff(x,k,i)=(long)b;

coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),

}

mulii(p4,gcoeff(x,k,i)));

if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));

coeff(x,k,i)=(long)b;

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

}

{

if (V) update(u,v,p3,p4,(GEN*)(V+i),(GEN*)(V+j));

if (DEBUGMEM>1) err(warnmem,"[1]: smithall");

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

if (all)

{

if (DEBUGMEM>1) err(warnmem,"[1]: smithall");

GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;

snf_pile(av, &x,&U,&V);

gerepilemany(av,gptr,3);

}

else x=gerepileupto(av, ZM_copy(x));

}

if (DEBUGLEVEL>=8) {fprintferr("; ");flusherr();}

if (DEBUGLEVEL>7) fprintferr("; ");

for (j=i-1; j>=1; j--)

{

p1=gcoeff(x,j,i); s1 = signe(p1);

if (s1)

if (!s1) continue;

{

p2=gcoeff(x,i,i);

if (!absi_cmp(p1,p2))

{

s2 = signe(p2);

if (s1 == s2) { d=p1; u=gun; p4=gun; }

else

{

s2 = signe(p2);

if (s2>0) { u = gun; p4 = mun; }

if (s1 == s2) { d=p1; u=gun; p4=gun; }

else { u = mun; p4 = gun; }

else

d=(s1>0)? p1: absi(p1);

{

if (s2>0) { u = gun; p4 = mun; }

else { u = mun; p4 = gun; }

d=(s1>0)? p1: absi(p1);

}

v = gzero; p3 = u;

}

else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }

v = gzero; p3 = u;

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

}

{

else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }

b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));

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

coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),

{

mulii(p4,gcoeff(x,i,k)));

b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));

coeff(x,i,k)=(long)b;

coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),

}

mulii(p4,gcoeff(x,i,k)));

if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));

coeff(x,i,k)=(long)b;

c = 1;

}

if (U) update(u,v,p3,p4,(GEN*)(U+i),(GEN*)(U+j));

c = 1;

}

if (!c)

{

b=gcoeff(x,i,i); fl=1;

b = gcoeff(x,i,i);

if (signe(b))

if (!signe(b)) break;

{

for (k=1; k<i && fl; k++)

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

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

{

fl = (int)!signe(resii(gcoeff(x,k,l),b));

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

/* cast to (int) necessary for gcc-2.95 on sparcv9-64 (IS) */

if (signe(resii(gcoeff(x,k,l),b))) break;

if (!fl)

if (l != i) break;

{

}

k--;

if (k == i) break;

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

coeff(x,i,l)=laddii(gcoeff(x,i,l),gcoeff(x,k,l));

/* x[k,l] != 0 mod b */

if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);

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

}

coeff(x,i,l) = laddii(gcoeff(x,i,l),gcoeff(x,k,l));

}

if (U) U[i] = ladd((GEN)U[i],(GEN)U[k]);

if (fl) break;

}

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

{

if (DEBUGMEM>1) err(warnmem,"[2]: smithall");

if (all)

snf_pile(av, &x,&U,&V);

{

GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;

gerepilemany(av,gptr,3);

}

else x=gerepileupto(av,ZM_copy(x));

}

if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();}

if (DEBUGLEVEL>7) fprintferr("\n");

if (all)

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

{

if (signe(gcoeff(x,k,k)) < 0)

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

{

if (signe(gcoeff(x,k,k))<0)

if (V) V[k]=lneg((GEN)V[k]);

{ mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lnegi(gcoeff(x,k,k)); }

coeff(x,k,k) = lnegi(gcoeff(x,k,k));

tetpil=avma; z=cgetg(4,t_VEC);

}

z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);

if (!U && !V)

return gerepile(av,tetpil,z);

return gerepileupto(av, mattodiagonal(x));

}

tetpil=avma; z=cgetg(n+1,t_VEC); j=n;

if (U) U = gtrans_i(U);

for (k=n; k; k--)

snf_pile(av, &x,&U,&V);

if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));

if (ptU) *ptU = U;

for ( ; j; j--) z[j]=zero;

if (ptV) *ptV = V;

return gerepile(av,tetpil,z);

return x;

}

GEN

smith(GEN x) { return smithall(x,0); }

smith(GEN x) { return smithall(x, NULL,NULL); }

GEN

smith2(GEN x) { return smithall(x,1); }

smith2(GEN x)

{

GEN z = cgetg(4, t_VEC);

z[3] = (long)smithall(x, (GEN*)(z+1),(GEN*)(z+2));

return z;

}

/* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */

GEN

Line 3136 smithclean(GEN z)

Line 3074 smithclean(GEN z)

static GEN

gsmithall(GEN x,long all)

{

long av,tetpil,i,j,k,l,c,n, lim;

gpmem_t av,tetpil,lim;

long i,j,k,l,c,n;

GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr;

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

Line 3254 gsmithall(GEN x,long all)

Line 3193 gsmithall(GEN x,long all)

GEN

matsnf0(GEN x,long flag)

{

long av = avma;

gpmem_t av = avma;

if (flag > 7) err(flagerr,"matsnf");

if (typ(x) == t_VEC && flag & 4) return smithclean(x);

if (flag & 2) x = gsmithall(x,flag & 1);

if (flag & 2) x = flag&1 ? gsmith2(x): gsmith(x);

else x = smithall(x, flag & 1);

else x = flag&1 ? smith2(x): smith(x);

if (flag & 4) x = smithclean(x);

if (flag & 4) x = gerepileupto(av, smithclean(x));

return gerepileupto(av, x);

return x;

}

GEN

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