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) |
|
Line 43 charpoly0(GEN x, int v, long flag) |
|
static GEN |
static GEN |
caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,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); |
GEN p1, p2 = leading_term(p); |
|
|
if (!signe(x)) return gpowgs(polx[v], degpol(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 |
else |
p1 = gsubst(p1,MAXVARN,polx[v]); |
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); |
return gerepileupto(av,p1); |
} |
} |
|
|
Line 75 caractducos(GEN p, GEN x, int v) |
|
Line 76 caractducos(GEN p, GEN x, int v) |
|
static GEN |
static GEN |
easychar(GEN x, int v, GEN *py) |
easychar(GEN x, int v, GEN *py) |
{ |
{ |
long l,tetpil,lx; |
gpmem_t av; |
|
long lx; |
GEN p1,p2; |
GEN p1,p2; |
|
|
switch(typ(x)) |
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: |
case t_COMPLEX: case t_QUAD: |
if (py) err(typeer,"easychar"); |
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: |
case t_POLMOD: |
if (py) err(typeer,"easychar"); |
if (py) err(typeer,"easychar"); |
Line 123 easychar(GEN x, int v, GEN *py) |
|
Line 125 easychar(GEN x, int v, GEN *py) |
|
GEN |
GEN |
caract(GEN x, int v) |
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; |
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++) |
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 |
GEN |
Line 155 caradj0(GEN x, long v) |
|
Line 159 caradj0(GEN x, long v) |
|
GEN |
GEN |
caradj(GEN x, long v, GEN *py) |
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]; |
GEN p,y,t,*gptr[2]; |
|
|
if ((p = easychar(x,v,py))) return p; |
if ((p = easychar(x,v,py))) return p; |
|
|
adj(GEN x) |
adj(GEN x) |
{ |
{ |
GEN y; |
GEN y; |
caradj(x,MAXVARN,&y); return y; |
(void)caradj(x,MAXVARN,&y); return y; |
} |
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
|
|
/* HESSENBERG FORM */ |
/* 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 |
GEN |
hess(GEN x) |
hess(GEN x) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long lx=lg(x),m,i,j; |
long lx=lg(x),m,i,j; |
GEN p,p1,p2; |
GEN p,p1,p2; |
|
|
|
|
p = gcoeff(x,i,m-1); |
p = gcoeff(x,i,m-1); |
if (gcmp0(p)) continue; |
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++) |
for (i=m+1; i<lx; i++) |
{ |
{ |
p1 = gcoeff(x,i,m-1); |
p1 = gcoeff(x,i,m-1); |
|
|
GEN |
GEN |
carhess(GEN x, long v) |
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; |
GEN *y,p1,p2,p3,p4; |
|
|
if ((p1 = easychar(x,v,NULL))) return p1; |
if ((p1 = easychar(x,v,NULL))) return p1; |
Line 299 carhess(GEN x, long v) |
|
Line 305 carhess(GEN x, long v) |
|
GEN |
GEN |
gnorm(GEN x) |
gnorm(GEN x) |
{ |
{ |
long l,lx,i,tetpil, tx=typ(x); |
gpmem_t av; |
|
long lx,i, tx=typ(x); |
GEN p1,p2,y; |
GEN p1,p2,y; |
|
|
switch(tx) |
switch(tx) |
|
|
case t_FRAC: case t_FRACN: |
case t_FRAC: case t_FRACN: |
return gsqr(x); |
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: |
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: |
case t_VEC: case t_COL: case t_MAT: |
lx=lg(x); y=cgetg(lx,tx); |
lx=lg(x); y=cgetg(lx,tx); |
for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]); |
for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]); |
|
|
GEN |
GEN |
gnorml2(GEN x) |
gnorml2(GEN x) |
{ |
{ |
|
gpmem_t av,lim; |
GEN y; |
GEN y; |
long i,tx=typ(x),lx,av,lim; |
long i,tx=typ(x),lx; |
|
|
if (! is_matvec_t(tx)) return gnorm(x); |
if (! is_matvec_t(tx)) return gnorm(x); |
lx=lg(x); if (lx==1) return gzero; |
lx=lg(x); if (lx==1) return gzero; |
|
|
GEN |
GEN |
QuickNormL2(GEN x, long prec) |
QuickNormL2(GEN x, long prec) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN y = gmul(x, realun(prec)); |
GEN y = gmul(x, realun(prec)); |
if (typ(x) == t_POL) |
if (typ(x) == t_POL) |
*++y = evaltyp(t_VEC) | evallg(lgef(x)-1); |
*++y = evaltyp(t_VEC) | evallg(lgef(x)-1); |
Line 378 QuickNormL2(GEN x, long prec) |
|
Line 389 QuickNormL2(GEN x, long prec) |
|
GEN |
GEN |
gnorml1(GEN x,long prec) |
gnorml1(GEN x,long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long lx,i; |
long lx,i; |
GEN s; |
GEN s; |
switch(typ(x)) |
switch(typ(x)) |
Line 406 gnorml1(GEN x,long prec) |
|
Line 417 gnorml1(GEN x,long prec) |
|
GEN |
GEN |
QuickNormL1(GEN x,long prec) |
QuickNormL1(GEN x,long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long lx,i; |
long lx,i; |
GEN p1,p2,s; |
GEN p1,p2,s; |
switch(typ(x)) |
switch(typ(x)) |
|
|
GEN |
GEN |
conjvec(GEN x,long prec) |
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; |
GEN z,y,p1,p2,p; |
|
|
switch(tx) |
switch(tx) |
Line 511 conjvec(GEN x,long prec) |
|
Line 523 conjvec(GEN x,long prec) |
|
case t_VEC: case t_COL: |
case t_VEC: case t_COL: |
lx=lg(x); z=cgetg(lx,t_MAT); |
lx=lg(x); z=cgetg(lx,t_MAT); |
for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec); |
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++) |
for (i=2; i<lx; i++) |
if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec"); |
if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec"); |
break; |
break; |
|
|
case t_POLMOD: |
case t_POLMOD: |
|
|
p1[j] = (j==(i+1))? un: zero; |
p1[j] = (j==(i+1))? un: zero; |
} |
} |
p1=cgetg(lx,t_COL); y[i]=(long)p1; |
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++) |
for (j=1; j<lx; j++) |
p1[j] = lneg((GEN)x[j+1]); |
p1[j] = lneg((GEN)x[j+1]); |
else |
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++) |
for (j=1; j<lx; j++) |
p1[j] = ldiv((GEN)x[j+1],p2); |
p1[j] = ldiv((GEN)x[j+1],p2); |
gnegz(p2,p2); |
gunclone(p2); |
} |
} |
return y; |
return y; |
} |
} |
|
|
GEN |
GEN |
gtrace(GEN x) |
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; |
GEN y,p1,p2; |
|
|
switch(tx) |
switch(tx) |
|
|
p1=(GEN)x[1]; |
p1=(GEN)x[1]; |
if (!gcmp0((GEN) p1[3])) |
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); |
return gmul2n((GEN)x[2],1); |
|
|
|
|
return y; |
return y; |
|
|
case t_POLMOD: |
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; |
p1=polsym((GEN)x[1],n); p2=gzero; |
for (i=0; i<=n; i++) |
for (i=0; i<=n; i++) |
p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1])); |
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: |
case t_RFRAC: case t_RFRACN: |
return gadd(x,gconj(x)); |
return gadd(x,gconj(x)); |
|
|
case t_MAT: |
case t_MAT: |
lx=lg(x); if (lx==1) return gzero;/*now lx>=2*/ |
lx=lg(x); if (lx==1) return gzero;/*now lx>=2*/ |
if (lx!=lg(x[1])) err(mattype1,"gtrace"); |
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++) |
for (i=2; i<lx-1; i++) |
p1=gadd(p1,gcoeff(x,i,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"); |
err(typeer,"gtrace"); |
return NULL; /* not reached */ |
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 a is not positive definite: |
* if flag is zero: raise an error |
* if flag is zero: raise an error |
* else: return NULL. |
* else: return NULL. */ |
*/ |
|
GEN |
GEN |
sqred1intern(GEN a,long flag) |
sqred1intern(GEN a,long flag) |
{ |
{ |
|
gpmem_t av = avma, lim=stack_lim(av,1); |
GEN b,p; |
GEN b,p; |
long i,j,k, n = lg(a); |
long i,j,k, n = lg(a); |
ulong av = avma, lim=stack_lim(av,1); |
|
|
|
if (typ(a)!=t_MAT) err(typeer,"sqred1"); |
if (typ(a)!=t_MAT) err(typeer,"sqred1"); |
if (n == 1) return cgetg(1, t_MAT); |
if (n == 1) return cgetg(1, t_MAT); |
|
|
GEN |
GEN |
sqred3(GEN a) |
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); |
long i,j,k,l, n = lg(a); |
GEN p1,b; |
GEN p1,b; |
|
|
|
|
sqred2(GEN a, long no_signature) |
sqred2(GEN a, long no_signature) |
{ |
{ |
GEN r,p,mun; |
GEN r,p,mun; |
ulong av,av1,lim; |
gpmem_t av,av1,lim; |
long n,i,j,k,l,sp,sn,t; |
long n,i,j,k,l,sp,sn,t; |
|
|
if (typ(a)!=t_MAT) err(typeer,"sqred2"); |
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) |
while (t) |
{ |
{ |
k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++; |
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 (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 |
GEN |
Line 824 signat(GEN a) { return sqred2(a,0); } |
|
Line 841 signat(GEN a) { return sqred2(a,0); } |
|
GEN |
GEN |
jacobi(GEN a, long prec) |
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; |
GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1; |
|
|
if (typ(a)!=t_MAT) err(mattype1,"jacobi"); |
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; |
r=cgetg(l,t_MAT); ja[2]=(long)r; |
for (j=1; j<=n; j++) |
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; |
av1=avma; |
|
|
Line 940 matrixqz0(GEN x,GEN p) |
|
Line 957 matrixqz0(GEN x,GEN p) |
|
err(flagerr,"matrixqz"); return NULL; /* not reached */ |
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 |
GEN |
matrixqz(GEN x, GEN p) |
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) err(talker,"more rows than columns in matrixqz"); |
if (n==m) |
if (n==m) |
{ |
{ |
p1=det(x); |
p1 = det(x); |
if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz"); |
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++) |
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)) |
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"); |
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++) |
for (i=1; i<=nfact; i++) |
{ |
{ |
p=(GEN)p1[i]; unmodp[1]=(long)p; |
p = (GEN)p1[i]; |
for(;;) |
for(;;) |
{ |
{ |
p2=ker(gmul(unmodp,x)); |
p2 = FpM_ker(x, p); |
if (lg(p2)==1) break; |
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++) |
for (j=1; j<lg(p2); j++) |
{ |
{ |
j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--; |
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 (low_stack(lim, stack_lim(av1,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"matrixqz"); |
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 |
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 (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"matrixqz_aux"); |
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 |
GEN |
matrixqz2(GEN x) |
matrixqz2(GEN x) |
{ |
{ |
long av = avma,m,n; |
gpmem_t av = avma; |
|
|
if (typ(x)!=t_MAT) err(typeer,"matrixqz2"); |
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 |
GEN |
matrixqz3(GEN x) |
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; |
GEN c; |
|
|
if (typ(x)!=t_MAT) err(typeer,"matrixqz3"); |
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)); |
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 |
GEN |
intersect(GEN x, GEN y) |
intersect(GEN x, GEN y) |
{ |
{ |
long av,tetpil,j, lx = lg(x); |
gpmem_t av,tetpil; |
|
long j, lx = lg(x); |
GEN z; |
GEN z; |
|
|
if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect"); |
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 |
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"); |
if (typ(x) != t_MAT) err(typeer,"mathnf"); |
*co=lg(x); if (*co==1) return NULL; /* empty matrix */ |
*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 */ |
if (gcmp1(*denx)) /* no denominator */ |
{ *denx = NULL; return dummycopy(x); } |
{ *denx = NULL; return dummycopy(x); } |
return gmul(x,*denx); |
return Q_muli_to_int(x, *denx); |
} |
} |
|
|
GEN |
GEN |
Line 1182 ZV_Z_mul(GEN c, GEN X) |
|
Line 1247 ZV_Z_mul(GEN c, GEN X) |
|
GEN |
GEN |
ZV_lincomb(GEN u, GEN v, GEN X, GEN Y) |
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; |
GEN p1,p2,A; |
|
|
if (!signe(u)) return ZV_Z_mul(v,Y); |
if (!signe(u)) return ZV_Z_mul(v,Y); |
if (!signe(v)) return ZV_Z_mul(u,X); |
if (!signe(v)) return ZV_Z_mul(u,X); |
lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4; |
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 (is_pm1(u)) |
{ |
{ |
if (signe(u) > 0) |
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 |
GEN |
hnf_special(GEN x, long remove) |
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; |
GEN p1,u,v,d,denx,a,b, x2,res; |
|
|
if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special"); |
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++) |
for (i=1,j=1; j<co; j++) |
if (!gcmp0((GEN)x[j])) |
if (!gcmp0((GEN)x[j])) |
{ |
{ |
x[i] = x[j]; |
x[i] = x[j]; |
x2[i] = x2[j]; i++; |
x2[i] = x2[j]; i++; |
} |
} |
setlg(x,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 p1,p2,U,H,Hnew,Bnew,Cnew,diagH1; |
GEN B = *ptB, C = *ptC, dep = *ptdep, depnew; |
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 co = lg(C); |
long col = lg(matgen)-1; |
long col = lg(matgen)-1; |
long lnz, nlze, lig; |
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); |
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 */ |
U = (GEN)p1[2]; /* col x col */ |
/* Only keep the part above the H (above the 0s is 0 since the dep rows |
/* Only keep the part above the H (above the 0s is 0 since the dep rows |
* are dependant from the ones in matgen) */ |
* 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) */ |
diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */ |
|
|
av = avma; lim = stack_lim(av,1); |
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); |
setlg(C, col+1); p1 = gmul(C,U); |
for (j=1; j<=col; j++) Cnew[j] = p1[j]; |
for (j=1; j<=col; j++) Cnew[j] = p1[j]; |
for ( ; j<co ; j++) Cnew[j] = C[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 */ |
/* Clean up B using new H */ |
for (s=0,i=lnz; i; i--) |
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++; } |
if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; } |
for (j=col+1; j<co; j++) |
for (j=col+1; j<co; j++) |
{ |
{ |
GEN z = (GEN)B[j-col]; |
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])); |
Cnew[j] = lsub((GEN)Cnew[j], gmul(p1, (GEN)Cnew[i+zc])); |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
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); |
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; |
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]; |
GEN z = (GEN)H[j]; |
if (diagH1[j]) |
if (diagH1[j]) |
{ /* hit exactly s times */ |
{ /* 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); |
for (i=1; i<=nlze; i++) p1[i] = coeff(dep,i,j); |
p1 += nlze; |
p1 += nlze; |
} |
} |
else |
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++) |
for (i=k=1; k<=lnz; i++) |
if (!diagH1[i]) p1[k++] = z[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); |
fprintferr("mit = %Z\nB = %Z\nC = %Z\n", Hnew, Bnew, C); |
} |
} |
} |
} |
if (nlze) *ptdep = depnew; |
*ptdep = depnew; |
*ptC = C; |
*ptC = C; |
*ptB = Bnew; return Hnew; |
*ptB = Bnew; return Hnew; |
} |
} |
|
|
/* for debugging */ |
/* for debugging */ |
static void |
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); |
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 |
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 |
** C = r x (co-1) matrix of GEN |
** perm= permutation vector (length li-1), indexing the rows of mat: easier |
** 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 |
** 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. |
** 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 |
** Also if k0=0, mat is modified in place [from mathnfspec], otherwise |
** it is left alone [from buchall] |
** it is left alone [from buchall] |
Line 1536 col_dup(long n, GEN col) |
|
Line 1596 col_dup(long n, GEN col) |
|
GEN |
GEN |
hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0) |
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 p1,p2,matb,matbnew,vmax,matt,T,extramat; |
GEN B,H,dep,permpro; |
GEN B,H,dep,permpro; |
GEN *gptr[4]; |
|
long co = lg(mat0); |
long co = lg(mat0); |
long li = lg(perm); /* = lg(mat0[1]) */ |
long li = lg(perm); /* = lg(mat0[1]) */ |
int updateT = 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) ) |
switch( count(mat,perm[i],col,&n) ) |
{ |
{ |
case 0: /* move zero lines between k0+1 and lk0 */ |
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; |
i=lig; continue; |
|
|
case 1: /* move trivial generator between lig+1 and li */ |
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 */ |
if (p[perm[lig]] < 0) /* = -1 */ |
{ /* convert relation -g = 0 to g = 0 */ |
{ /* convert relation -g = 0 to g = 0 */ |
for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]]; |
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;} |
#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 |
/* 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 |
* 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; |
s = 0; |
while (lig > lk0 && s < (HIGHBIT>>1)) |
while (lig > lk0 && s < (HIGHBIT>>1)) |
{ |
{ |
for (i=lig; i>lk0; i--) |
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; |
if (i == lk0) break; |
|
|
/* only 0, ±1 entries, at least 2 of them non-zero */ |
/* 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) |
if (p[perm[lig]] < 0) |
{ |
{ |
for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]]; |
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]; |
matj = mat[j]; |
if (! (t = matj[perm[lig]]) ) continue; |
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); |
T = gerepilecopy(av2, T); |
} |
} |
} |
} |
#endif |
|
/* As above with lines containing a ±1 (no other assumption). |
/* As above with lines containing a ±1 (no other assumption). |
* Stop when single precision becomes dangerous */ |
* Stop when single precision becomes dangerous */ |
for (j=1; j<=col; j++) |
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 ( (n = count2(mat,perm[i],col)) ) break; |
if (i == lk0) 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) |
if (p[perm[lig]] < 0) |
{ |
{ |
for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]]; |
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) |
if (DEBUGLEVEL>5) |
{ |
{ |
fprintferr(" after phase2:\n"); |
fprintferr(" after phase2:\n"); |
p_mat(mat,perm,k0); |
p_mat(mat,perm,lk0); |
} |
} |
for (i=li-2; i>lig; i--) |
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); |
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 */ |
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) |
if (DEBUGLEVEL>5) |
{ |
{ |
fprintferr(" matb cleaned up (using Id block)\n"); |
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; |
*ptdep = dep; |
*ptB = B; |
*ptB = B; |
H = hnffinal(matbnew,perm,ptdep,ptB,ptC); |
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) |
if (DEBUGLEVEL) |
msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1); |
msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1); |
return H; |
return H; |
|
|
hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */ |
hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */ |
GEN extramat,GEN extraC) |
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 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; |
long nlze = lH? lg(dep[1])-1: lg(B[1])-1; |
|
|
if (DEBUGLEVEL>5) |
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 |
* [--------|-----] lig |
* [ 0 | Id ] |
* [ 0 | Id ] |
* [ | ] li */ |
* [ | ] 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 */ |
{ /* 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); |
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"); |
if (DEBUGLEVEL>5) fprintferr(" 2nd phase done\n"); |
*ptdep = dep; |
|
*ptB = B; |
|
H = hnffinal(matb,perm,ptdep,ptB,&Cnew); |
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) |
{ |
{ |
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; |
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; |
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); |
d = bezout(aj,ak,&u,&v); |
/* frequent special case (u,v) = (1,0) or (0,1) */ |
/* frequent special case (u,v) = (1,0) or (0,1) */ |
if (!signe(u)) |
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 */ |
{ /* aj | ak */ |
p1 = negi(divii(ak,aj)); |
p1 = negi(divii(ak,aj)); |
A[k] = (long)ZV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]); |
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) { |
if (U) { |
U[k] = (long)ZV_lincomb(gun, p1, (GEN)U[k], (GEN)U[j]); |
U[k] = (long)ZV_lincomb(gun, p1, (GEN)U[k], (GEN)U[j]); |
swap(U[j], U[k]); |
lswap(U[j], U[k]); |
} |
} |
return; |
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 |
GEN |
hnf0(GEN A, int remove) |
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; |
GEN denx,a,p1; |
|
|
if (typ(A) == t_VEC) return hnf_special(A,remove); |
if (typ(A) == t_VEC) return hnf_special(A,remove); |
Line 2100 hnf0(GEN A, int remove) |
|
Line 2100 hnf0(GEN A, int remove) |
|
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i); |
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); |
p1 = gcoeff(A,i,def); s = signe(p1); |
Line 2115 hnf0(GEN A, int remove) |
|
Line 2115 hnf0(GEN A, int remove) |
|
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i); |
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) |
if (remove) |
Line 2124 hnf0(GEN A, int remove) |
|
Line 2124 hnf0(GEN A, int remove) |
|
if (!gcmp0((GEN)A[j])) A[i++] = A[j]; |
if (!gcmp0((GEN)A[j])) A[i++] = A[j]; |
setlg(A,i); |
setlg(A,i); |
} |
} |
tetpil=avma; |
|
A = denx? gdiv(A,denx): ZM_copy(A); |
A = denx? gdiv(A,denx): ZM_copy(A); |
return gerepile(av0,tetpil,A); |
return gerepileupto(av0, A); |
} |
} |
|
|
GEN |
GEN |
Line 2137 hnf(GEN x) { return hnf0(x,1); } |
|
Line 2136 hnf(GEN x) { return hnf0(x,1); } |
|
static GEN |
static GEN |
allhnfmod(GEN x,GEN dm,int flag) |
allhnfmod(GEN x,GEN dm,int flag) |
{ |
{ |
ulong av,tetpil,lim; |
gpmem_t av, lim; |
long li,co,i,j,k,def,ldef,ldm; |
long li,co,i,j,k,def,ldef,ldm; |
GEN a,b,w,p1,p2,d,u,v; |
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 (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i); |
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; |
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 (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i); |
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 |
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 |
static GEN |
mynegi(GEN x) |
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); |
long s = signe(x); |
|
|
if (!s) return gzero; |
if (!s) return gzero; |
|
|
setsigne(x,-s); return 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 |
static void |
Minus(long j, GEN **lambda) |
Minus(long j, GEN **lambda) |
{ |
{ |
|
|
return 0; |
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 |
static void |
reduce2(GEN A, GEN B, long k, long j, long *row, GEN **lambda, GEN *D) |
reduce2(GEN A, GEN B, long k, long j, long *row, GEN **lambda, GEN *D) |
{ |
{ |
GEN q; |
GEN q; |
long i, row0, row1; |
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) |
if (row0) |
q = truedvmdii(gcoeff(A,row0,k), gcoeff(A,row0,j), NULL); |
q = truedvmdii(gcoeff(A,row0,k), gcoeff(A,row0,j), NULL); |
else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0) |
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 |
else |
return; |
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 |
static void |
hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D) |
hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D) |
{ |
{ |
GEN t,p1,p2; |
GEN t,p1,p2; |
long i,j,n = lg(A); |
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++) |
for (i=k+1; i<n; i++) |
{ |
{ |
p1 = mulii(lambda[i][k-1], D[k]); |
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]); |
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 |
static GEN |
fix_rows(GEN A) |
fix_rows(GEN A) |
{ |
{ |
Line 2443 fix_rows(GEN A) |
|
Line 2394 fix_rows(GEN A) |
|
} |
} |
|
|
GEN |
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 m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */ |
long row[2], do_swap,i,n,k; |
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 z,B, **lambda, *D; |
GEN *gptr[4]; |
|
|
|
if (typ(A) != t_MAT) err(typeer,"hnflll"); |
if (typ(A) != t_MAT) err(typeer,"hnflll"); |
n = lg(A); |
n = lg(A); |
A = ZM_copy(fix_rows(A)); |
A = ZM_copy(fix_rows(A)); |
B = ptB? idmat(n-1): NULL; |
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) |
while (k < n) |
{ |
{ |
reduce2(A,B,k,k-1,row,lambda,D); |
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 |
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) |
if (do_swap) |
{ |
{ |
hnfswap(A,B,k,lambda,D); |
hnfswap(A,B,k,lambda,D); |
if (k > nzcol+1) k--; |
if (k > 2) k--; |
} |
} |
else |
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++; |
k++; |
} |
} |
if (low_stack(lim, stack_lim(av,3))) |
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); |
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; |
if (B) *ptB = B; |
return A; |
return A; |
} |
} |
|
|
hnflll(GEN A) |
hnflll(GEN A) |
{ |
{ |
GEN B, z = cgetg(3, t_VEC); |
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; |
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; |
long i; |
|
|
if (signe(A[j])) |
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) |
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 |
else |
return; |
return; |
|
|
if (! gcmp0(q)) |
if (signe(q)) |
{ |
{ |
|
GEN *Lk = lambda[k], *Lj = lambda[j]; |
q = mynegi(q); |
q = mynegi(q); |
A[k] = laddii((GEN)A[k], mulii(q,(GEN)A[j])); |
A[k] = laddii((GEN)A[k], mulii(q,(GEN)A[j])); |
elt_col((GEN)B[k],(GEN)B[j],q); |
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++) |
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])); |
|
} |
} |
} |
} |
|
|
|
|
extendedgcd(GEN A) |
extendedgcd(GEN A) |
{ |
{ |
long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */ |
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=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; |
k = 2; |
while (k < n) |
while (k < n) |
{ |
{ |
Line 2586 extendedgcd(GEN A) |
|
Line 2523 extendedgcd(GEN A) |
|
if (signe(A[k-1])) do_swap = 1; |
if (signe(A[k-1])) do_swap = 1; |
else if (!signe(A[k])) |
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; |
else do_swap = 0; |
|
|
if (do_swap) |
if (do_swap) |
{ |
{ |
hnfswap(A,B,k,lambda,D); |
hnfswap(A,B,k,lambda,D); |
if (k>2) k--; |
if (k > 2) k--; |
} |
} |
else |
else |
{ |
{ |
Line 2606 extendedgcd(GEN A) |
|
Line 2541 extendedgcd(GEN A) |
|
k++; |
k++; |
} |
} |
} |
} |
if (signe(A[n-1])<0) |
if (signe(A[n-1]) < 0) |
{ |
{ |
A[n-1] = (long)mynegi((GEN)A[n-1]); |
A[n-1] = (long)mynegi((GEN)A[n-1]); |
neg_col((GEN)B[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. */ |
/* HNF with permutation. TODO: obsolete, replace with hnflll. */ |
|
|
hnfperm(GEN A) |
hnfperm(GEN A) |
{ |
{ |
GEN U,c,l,perm,d,u,p,q,y,b; |
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"); |
if (typ(A) != t_MAT) err(typeer,"hnfperm"); |
n = lg(A)-1; |
n = lg(A)-1; |
|
|
} |
} |
if (low_stack(lim, stack_lim(av1,1))) |
if (low_stack(lim, stack_lim(av1,1))) |
{ |
{ |
GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U; |
|
if (DEBUGMEM>1) err(warnmem,"hnfperm"); |
if (DEBUGMEM>1) err(warnmem,"hnfperm"); |
gerepilemany(av1,gptr,2); |
gerepileall(av1, 2, &A, &U); |
} |
} |
} |
} |
if (r < m) |
if (r < m) |
|
|
hnfall_i(GEN A, GEN *ptB, long remove) |
hnfall_i(GEN A, GEN *ptB, long remove) |
{ |
{ |
GEN B,c,h,x,p,a; |
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"); |
if (typ(A)!=t_MAT) err(typeer,"hnfall"); |
n=lg(A)-1; |
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); |
ZM_reduce(A,B, i,k); |
if (low_stack(lim, stack_lim(av1,1))) |
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); |
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; |
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--; |
r--; |
if (j < r) /* A[j] != 0 */ |
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; |
h[j]=h[r]; h[r]=li; c[li]=r; |
} |
} |
p = gcoeff(A,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); |
ZM_reduce(A,B, li,r); |
if (low_stack(lim, stack_lim(av1,1))) |
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); |
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: "); |
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); |
ZM_reduce(A,B, i,k); |
if (low_stack(lim, stack_lim(av1,1))) |
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); |
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"); |
if (DEBUGLEVEL>5) fprintferr("\n"); |
/* remove the first r columns */ |
/* remove the first r columns */ |
if (remove) { A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1); } |
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; |
if (B) *ptB = B; |
return A; |
return A; |
} |
} |
Line 2896 trivsmith(long all) |
|
Line 2822 trivsmith(long all) |
|
z[3]=lgetg(1,t_MAT); return z; |
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 (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"); |
if (lg(x[1]) != n+1) err(mattype1,"smithall"); |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
for (j=1; j<=n; j++) |
for (j=1; j<=n; j++) |
if (typ(coeff(x,i,j)) != t_INT) |
if (typ(coeff(x,i,j)) != t_INT) |
err(talker,"non integral matrix in smithall"); |
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 |
else |
{ |
{ |
|
mdet = detint(x); |
if (signe(mdet)) |
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 |
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; |
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)); |
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++) |
for (j=1; j<=n; j++) |
{ |
{ |
p1=cgetg(n+1,t_COL); ys[j]=(long)p1; |
p1=cgetg(n+1,t_COL); ys[j]=(long)p1; |
for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]); |
for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]); |
} |
} |
x = ys; |
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)) |
if (signe(mdet)) |
{ |
{ |
p1 = hnfmod(x,mdet); |
p1 = hnfmod(x,mdet); |
if (all) mr=gmul(mr,gauss(x,p1)); |
if (V) V = gmul(V, gauss(x,p1)); |
} |
} |
else |
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(;;) |
for(;;) |
{ |
{ |
c = 0; |
c = 0; |
for (j=i-1; j>=1; j--) |
for (j=i-1; j>=1; j--) |
{ |
{ |
p1=gcoeff(x,i,j); s1 = signe(p1); |
p1=gcoeff(x,i,j); s1 = signe(p1); |
if (s1) |
if (!s1) continue; |
{ |
|
p2=gcoeff(x,i,i); |
p2=gcoeff(x,i,i); |
if (!absi_cmp(p1,p2)) |
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--) |
for (j=i-1; j>=1; j--) |
{ |
{ |
p1=gcoeff(x,j,i); s1 = signe(p1); |
p1=gcoeff(x,j,i); s1 = signe(p1); |
if (s1) |
if (!s1) continue; |
{ |
|
p2=gcoeff(x,i,i); |
p2=gcoeff(x,i,i); |
if (!absi_cmp(p1,p2)) |
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) |
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 (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if (DEBUGMEM>1) err(warnmem,"[2]: smithall"); |
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 |
GEN |
smith(GEN x) { return smithall(x,0); } |
smith(GEN x) { return smithall(x, NULL,NULL); } |
|
|
GEN |
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 */ |
/* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */ |
GEN |
GEN |
Line 3136 smithclean(GEN z) |
|
Line 3074 smithclean(GEN z) |
|
static GEN |
static GEN |
gsmithall(GEN x,long all) |
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; |
GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr; |
|
|
if (typ(x)!=t_MAT) err(typeer,"gsmithall"); |
if (typ(x)!=t_MAT) err(typeer,"gsmithall"); |
Line 3254 gsmithall(GEN x,long all) |
|
Line 3193 gsmithall(GEN x,long all) |
|
GEN |
GEN |
matsnf0(GEN x,long flag) |
matsnf0(GEN x,long flag) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
if (flag > 7) err(flagerr,"matsnf"); |
if (flag > 7) err(flagerr,"matsnf"); |
if (typ(x) == t_VEC && flag & 4) return smithclean(x); |
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 |
GEN |