=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.c,v retrieving revision 1.1.1.1 retrieving revision 1.2 diff -u -p -r1.1.1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.c 2001/10/02 11:17:00 1.1.1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin2.c 2002/09/11 07:26:48 1.2 @@ -1,4 +1,4 @@ -/* $Id: alglin2.c,v 1.1.1.1 2001/10/02 11:17:00 noro Exp $ +/* $Id: alglin2.c,v 1.2 2002/09/11 07:26:48 noro Exp $ Copyright (C) 2000 The PARI group. @@ -43,7 +43,8 @@ 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)); @@ -55,7 +56,7 @@ 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); } @@ -75,7 +76,8 @@ 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)) @@ -96,11 +98,11 @@ easychar(GEN x, int v, GEN *py) case t_COMPLEX: case t_QUAD: if (py) err(typeer,"easychar"); - p1=cgetg(5,t_POL); - p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v); - p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma; - p1[3]=lpile(l,tetpil,gneg(p2)); - p1[4]=un; return p1; + p1 = cgetg(5,t_POL); + p1[1] = evalsigne(1) | evallgef(5) | evalvarn(v); + p1[2] = lnorm(x); av = avma; + p1[3] = lpileupto(av, gneg(gtrace(x))); + p1[4] = un; return p1; case t_POLMOD: if (py) err(typeer,"easychar"); @@ -123,24 +125,26 @@ easychar(GEN x, int v, GEN *py) GEN caract(GEN x, int v) { - long n,k,l=avma,tetpil; - GEN p1,p2,p3,p4,p5,p6; + long k, n; + gpmem_t av=avma; + GEN p1,p2,p3,p4,x_k; if ((p1 = easychar(x,v,NULL))) return p1; - p1=gzero; p2=gun; - n=lg(x)-1; if (n&1) p2=gneg_i(p2); - p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5; - p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3); + p1 = gzero; p2 = gun; + n = lg(x)-1; if (n&1) p2 = negi(p2); + x_k = dummycopy(polx[v]); + p4 = cgetg(3,t_RFRACN); p4[2] = (long)x_k; for (k=0; k<=n; k++) { - p3=det(gsub(gscalmat(stoi(k),n), x)); - p4[1]=lmul(p3,p2); p6[2]=k; - p1=gadd(p4,p1); p5[2]=(long)p6; - if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1); + p3 = det(gsub(gscalmat(stoi(k),n), x)); + p4[1] = lmul(p3,p2); x_k[2] = lstoi(-k); + p1 = gadd(p4,p1); + if (k == n) break; + + p2 = gdivgs(gmulsg(k-n,p2),k+1); } - p2=mpfact(n); tetpil=avma; - return gerepile(l,tetpil,gdiv((GEN) p1[1],p2)); + return gerepileupto(av, gdiv((GEN)p1[1], mpfact(n))); } GEN @@ -155,7 +159,8 @@ 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; @@ -217,7 +222,7 @@ GEN adj(GEN x) { GEN y; - caradj(x,MAXVARN,&y); return y; + (void)caradj(x,MAXVARN,&y); return y; } /*******************************************************************/ @@ -225,13 +230,13 @@ adj(GEN x) /* HESSENBERG FORM */ /* */ /*******************************************************************/ -#define swap(x,y) { long _t=x; x=y; y=_t; } -#define gswap(x,y) { GEN _t=x; x=y; y=_t; } +#define lswap(x,y) { long _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; @@ -246,8 +251,8 @@ hess(GEN x) p = gcoeff(x,i,m-1); if (gcmp0(p)) continue; - for (j=m-1; j=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 1 && lg(a[1]) != n) err(mattype1,"sqred2"); - av=avma; mun=negi(gun); r=new_chunk(n); - for (i=1; in) break; } } - if (no_signature) return gerepilecopy(av,a); - avma=av; - a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a; + if (no_signature) return gerepilecopy(av, a); + avma = av; a = cgetg(3,t_VEC); + a[1] = lstoi(sp); + a[2] = lstoi(sn); return a; } GEN @@ -824,7 +841,8 @@ 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"); @@ -840,9 +858,8 @@ 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); - for (i=1; i<=n; i++) - affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec))); + r[j] = lgetg(l,t_COL); + for (i=1; i<=n; i++) coeff(r,i,j) = (long)stor(i==j, prec); } av1=avma; @@ -940,63 +957,63 @@ 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 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]; - for (i=1; i<=m; i++) - { - 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[j] = (long)primpart((GEN)p1[j]); + if (!ZV_isin((GEN)p1[j])) err(talker, "not a rational matrix in matrixqz"); } - x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un; + /* x integral */ if (gcmp0(p)) { - p1=cgetg(n+1,t_MAT); p2=gtrans(x); - for (j=1; j<=n; j++) p1[j]=p2[j]; - p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1)); - if (!signe(p3)) + p1 = gtrans(x); setlg(p1,n+1); + p2 = det(p1); p1[n] = p1[n+1]; p2 = ggcd(p2,det(p1)); + if (!signe(p2)) 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 - { - p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1; - } - av1=avma; lim=stack_lim(av1,1); + else p1 = _vec(p); + nfact = lg(p1)-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; j1) err(warnmem,"matrixqz"); - x=gerepilecopy(av1,x); + x = gerepilecopy(av1,x); } } } @@ -1013,90 +1030,138 @@ 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); - long i,j,k,in[2]; - GEN p1; + if (gcmp1(u)) return A; + if (gcmp_1(u)) return gneg(A); + 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; in) 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])); + k = j+1; if (k == n) k = 1; + /* zero a = Aij using b = Aik */ + QV_elem(a, gcoeff(A,i,k), A, j,k); } - - j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++; - if (j<=n) + a = gcoeff(A,i,k); + if (!gcmp0(a)) { - p1=denom(gcoeff(x,i,j)); - if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]); + a = denom(a); + 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); - m=lg(x[1])-1; x=dummycopy(x); - return gerepileupto(av, matrixqz_aux(x,m,n)); + x = dummycopy(x); + return gerepileupto(av, matrixqz_aux(x)); } 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); - m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1); - for (j=1; j<=n; j++) c[j]=0; - av1=avma; lim=stack_lim(av1,1); - for (k=1; k<=m; k++) + n = lg(x); if (n==1) return gcopy(x); + m = lg(x[1]); x = dummycopy(x); + c = cgetg(n,t_VECSMALL); + for (j=1; j1) err(warnmem,"matrixqz3"); - x=gerepilecopy(av1,x); + GEN t = gcoeff(x,k,j1); + 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"); @@ -1128,7 +1193,7 @@ 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 */ @@ -1136,7 +1201,7 @@ 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 @@ -1182,13 +1247,14 @@ 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) @@ -1243,7 +1309,8 @@ 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"); @@ -1311,7 +1378,7 @@ hnf_special(GEN x, long remove) for (i=1,j=1; j1) err(warnmem,"hnffinal, i = %ld",i); - gerepilemany(av,gptr,2); + gerepileall(av, 2, &Cnew, &B); } } p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze; @@ -1452,16 +1521,16 @@ 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; - C[i1+col] = Cnew[j+zc]; + i1++; 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; - C[j1+zc] = Cnew[j+zc]; - if (nlze) depnew[j1] = dep[j]; + j1++; C[j1+zc] = Cnew[j+zc]; + p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1; + depnew[j1] = dep[j]; } for (i=k=1; k<=lnz; i++) if (!diagH1[i]) p1[k++] = z[i]; @@ -1484,30 +1553,21 @@ 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; - GEN p1, matj, matgen; - long co = lg(mat); - long li = lg(perm); - + gpmem_t av = avma; + perm = vecextract_i(perm, k+1, lg(perm)-1); fprintferr("Permutation: %Z\n",perm); - matgen = cgetg(co,t_MAT); - for (j=1; j 6) fprintferr("matgen = %Z\n",matgen); - avma=av; + if (DEBUGLEVEL > 6) + fprintferr("matgen = %Z\n", small_to_mat( rowextract_p((GEN)mat, perm) )); + avma = av; } static GEN @@ -1524,7 +1584,7 @@ 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] @@ -1536,11 +1596,11 @@ 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; - long n,s,t,lim,nlze,lnz,nr; + gpmem_t av=avma,av2,lim; + 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; @@ -1579,13 +1639,13 @@ 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]); - swap(T[n], T[col]); - gswap(mat[n], mat[col]); p = mat[col]; + lswap(perm[i], perm[lig]); + lswap(T[n], T[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 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]); - swap(T[n], T[col]); p1 = (GEN)T[col]; - gswap(mat[n], mat[col]); p = mat[col]; + lswap(perm[i], perm[lig]); + swap(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; j5) { fprintferr(" after phase2:\n"); - p_mat(mat,perm,k0); + p_mat(mat,perm,lk0); } for (i=li-2; i>lig; i--) { @@ -1748,10 +1804,10 @@ END2: /* clean up mat: remove everything to the right { if(DEBUGMEM>1) err(warnmem,"hnfspec[3], i = %ld", i); for (j=1; j5) { fprintferr(" matb cleaned up (using Id block)\n"); @@ -1829,10 +1885,7 @@ END2: /* clean up mat: remove everything to the right *ptdep = dep; *ptB = B; H = hnffinal(matbnew,perm,ptdep,ptB,ptC); - gptr[0]=ptC; - gptr[1]=ptdep; - gptr[2]=ptB; - gptr[3]=&H; gerepilemany(av,gptr,4); + gerepileall(av, 4, ptC, ptdep, ptB, &H); if (DEBUGLEVEL) msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1); return H; @@ -1891,15 +1944,13 @@ 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 B=*ptB, C=*ptC, dep=*ptdep, *gptr[4]; - long av = avma, i,j,lextra,colnew; - long li = lg(perm); - long co = lg(C); - long lB = lg(B); - long lig = li - lB; - long col = co - lB; + GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep; + gpmem_t av = avma; + long i; 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) @@ -1914,89 +1965,37 @@ hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC * [--------|-----] lig * [ 0 | Id ] * [ | ] li */ - lextra = lg(extramat)-1; - extratop = cgetg(lextra+1,t_MAT); /* [1..lig] part (top) */ - 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 ( ; i5) - { - fprintferr(" 1st phase done\n"); - if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat); - } + extramat = concatsp(extratop, vconcat(dep, H)); + Cnew = concatsp(extraC, vecextract_i(C, col-lH+1, co)); + if (DEBUGLEVEL>5) fprintferr(" 1st phase done\n"); permpro = imagecomplspec(extramat, &nlze); - p1 = new_chunk(lig+1); - for (i=1; i<=lig; i++) p1[i] = perm[permpro[i]]; - for (i=1; i<=lig; i++) perm[i] = p1[i]; + extramat = rowextract_p(extramat, permpro); + *ptB = rowextract_p(B, permpro); + /* perm o= permpro */ + permpro = vecextract_p(perm, permpro); + for (i=1; i<=lig; i++) perm[i] = permpro[i]; - matb = cgetg(colnew+1,t_MAT); - dep = cgetg(colnew+1,t_MAT); - 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; j5) fprintferr(" 2nd phase done\n"); - *ptdep = dep; - *ptB = B; H = hnffinal(matb,perm,ptdep,ptB,&Cnew); - p1 = cgetg(co+lextra,t_MAT); - for (j=1; j <= col-lH; j++) p1[j] = C[j]; - C = Cnew - (col-lH); - for ( ; j < co+lextra; j++) p1[j] = C[j]; + *ptC = concatsp(vecextract_i(C, 1, col-lH), Cnew); - 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); - msgtimer("hnfadd (%ld)",lextra); + if (DEBUGLEVEL>7) fprintferr("H = %Z\nC = %Z\n",H,*ptC); + msgtimer("hnfadd (%ld)", lg(extramat)-1); } + gerepileall(av, 4, ptC, ptdep, ptB, &H); return H; } @@ -2014,7 +2013,7 @@ 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)) @@ -2029,10 +2028,10 @@ 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; } @@ -2077,7 +2076,8 @@ 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); @@ -2100,7 +2100,7 @@ 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); @@ -2115,7 +2115,7 @@ 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) @@ -2124,9 +2124,8 @@ 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 @@ -2137,7 +2136,7 @@ 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; @@ -2184,7 +2183,7 @@ 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; @@ -2214,12 +2213,11 @@ 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)); - diag = gcoeff(w,i,i); + gerepileall(av, 2, &w, &dm); diag = gcoeff(w,i,i); } } } - tetpil=avma; return gerepile(av,tetpil,ZM_copy(w)); + return gerepilecopy(av, w); } GEN @@ -2238,7 +2236,7 @@ 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; @@ -2247,22 +2245,6 @@ 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) { @@ -2313,33 +2295,30 @@ 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]); - if (row0 && signe(gcoeff(A,row0,j)) < 0) - { - 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; + row[0] = row0 = findi_normalize((GEN)A[j], B,j,lambda); + row[1] = row1 = findi_normalize((GEN)A[k], B,k,lambda); 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; @@ -2371,17 +2350,15 @@ 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]); - if (B) swap(B[k],B[k-1]); - for (j=k-2; j; j--) SWAP(lambda[k-1][j], lambda[k][j]); + lswap(A[k],A[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 (i=k+1; i nz, A[j] = 0 for all j < nz. - * "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 nzcol) - {trivswap((GEN*)A,nzcol,k-1, lambda,D); nzcol++;} - 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]); + if (row[0]) + do_swap = (!row[1] || row[0] <= row[1]); + else if (!row[1]) + { /* row[0] == row[1] == 0 */ + gpmem_t 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; } else - { - 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; - } - } + do_swap = 0; 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 */ - gptr[0]=&A; gptr[1]=&a; gptr[2]=&b; gptr[3]=&B; - if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n); - gerepilemany(av,gptr,B? 4: 3); lambda = (GEN**)a; D = (GEN*)(b+1); + GEN b = (GEN)(D-1); + if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n-1); + gerepileall(av, B? 4: 3, &A, (GEN*)&lambda, &b, &B); + D = (GEN*)(b+1); } } - for (i=nzcol; i 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; i2) k--; + if (k > 2) k--; } else { @@ -2606,16 +2541,15 @@ 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; - p1 = cgetg(3,t_VEC); - p1[1]=lcopy((GEN)A[n-1]); - p1[2]=lcopy(B); - return gerepile(av,tetpil,p1); + z = cgetg(3,t_VEC); + z[1] = A[n-1]; + z[2] = (long)B; + return gerepilecopy(av, z); } /* HNF with permutation. TODO: obsolete, replace with hnflll. */ @@ -2623,7 +2557,8 @@ 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; @@ -2696,9 +2631,8 @@ 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) @@ -2735,8 +2669,8 @@ 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; - GEN *gptr[2]; + gpmem_t av = avma, av1, lim; + long m,n,r,i,j,k,li; if (typ(A)!=t_MAT) err(typeer,"hnfall"); n=lg(A)-1; @@ -2766,11 +2700,8 @@ 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; @@ -2780,8 +2711,8 @@ hnfall_i(GEN A, GEN *ptB, long remove) r--; if (j < r) /* A[j] != 0 */ { - swap(A[j], A[r]); - if (B) swap(B[j], B[r]); + lswap(A[j], A[r]); + if (B) lswap(B[j], B[r]); h[j]=h[r]; h[r]=li; c[li]=r; } p = gcoeff(A,li,r); @@ -2794,9 +2725,8 @@ 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: "); @@ -2812,18 +2742,14 @@ 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; - gerepilemany(av, gptr, B? 2: 1); + gerepileall(av, B? 2: 1, &A, &B); if (B) *ptB = B; return A; } @@ -2896,208 +2822,220 @@ 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], - * where d = u.x.v - */ -static GEN -smithall(GEN x, long all) +static void +snf_pile(gpmem_t av, GEN *x, GEN *U, GEN *V) { - long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim; - GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys; + GEN *gptr[3]; + 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); - n=lg(x)-1; - if (!n) return trivsmith(all); + if (DEBUGLEVEL>8) outerr(x); + n = lg(x)-1; + 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); - x = dummycopy(x); mdet = detint(x); - if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } } + U = ptU? gun: NULL; + V = ptV? gun: NULL; + 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); - if (all) { ml=idmat(n); mr=gauss(x,p1); } + p1 = hnfmod(x,mdet); + if (V) V = gauss(x,p1); } else { - if (all) - { - p1 = hnfall0(x,0); - ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1]; - } - else - p1 = hnf0(x,0); + p1 = hnfall_i(x, ptV, 0); + if (ptV) V = *ptV; } 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); - for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; } - ml=p3; mr=p4; + p1 = cgetg(n+1,t_MAT); + for (j=1; j<=n; j++) p1[j]=U[p2[j]]; + 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 = hnfall0(x,0); - mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1]; - } - else - p1 = hnf0(x,0); + p1 = hnfall_i(x,ptV,0); + if (ptV) V = gmul(V, *ptV); } - x=p1; mun = negi(gun); + x = p1; mun = negi(gun); - if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();} - for (i=n; i>=2; i--) + if (DEBUGLEVEL>7) fprintferr("starting SNF loop"); + 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) - { - p2=gcoeff(x,i,i); - if (!absi_cmp(p1,p2)) + 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 (s1 == s2) { d=p1; u=gun; p4=gun; } - else - { - if (s2>0) { u = gun; p4 = mun; } - else { u = mun; p4 = gun; } - d=(s1>0)? p1: absi(p1); - } - v = gzero; p3 = u; + if (s2>0) { u = gun; p4 = mun; } + else { u = mun; p4 = gun; } + d=(s1>0)? p1: absi(p1); } - else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); } - for (k=1; k<=i; k++) - { - b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j))); - coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)), - mulii(p4,gcoeff(x,k,i))); - coeff(x,k,i)=(long)b; - } - if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j)); - if (low_stack(lim, stack_lim(av,1))) - { - if (DEBUGMEM>1) err(warnmem,"[1]: smithall"); - if (all) - { - GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr; - gerepilemany(av,gptr,3); - } - else x=gerepileupto(av, ZM_copy(x)); - } - } + v = gzero; p3 = u; + } + else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); } + for (k=1; k<=i; k++) + { + b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j))); + coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)), + mulii(p4,gcoeff(x,k,i))); + coeff(x,k,i)=(long)b; + } + if (V) update(u,v,p3,p4,(GEN*)(V+i),(GEN*)(V+j)); + if (low_stack(lim, stack_lim(av,1))) + { + if (DEBUGMEM>1) err(warnmem,"[1]: smithall"); + snf_pile(av, &x,&U,&V); + } } - 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) - { - p2=gcoeff(x,i,i); - if (!absi_cmp(p1,p2)) + 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 (s1 == s2) { d=p1; u=gun; p4=gun; } - else - { - if (s2>0) { u = gun; p4 = mun; } - else { u = mun; p4 = gun; } - d=(s1>0)? p1: absi(p1); - } - v = gzero; p3 = u; + if (s2>0) { u = gun; p4 = mun; } + else { u = mun; p4 = gun; } + d=(s1>0)? p1: absi(p1); } - else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); } - for (k=1; k<=i; k++) - { - b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k))); - coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)), - mulii(p4,gcoeff(x,i,k))); - coeff(x,i,k)=(long)b; - } - if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j)); - c = 1; - } + v = gzero; p3 = u; + } + else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); } + for (k=1; k<=i; k++) + { + b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k))); + coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)), + mulii(p4,gcoeff(x,i,k))); + coeff(x,i,k)=(long)b; + } + if (U) update(u,v,p3,p4,(GEN*)(U+i),(GEN*)(U+j)); + c = 1; } if (!c) { - b=gcoeff(x,i,i); fl=1; - if (signe(b)) - { - for (k=1; k1) err(warnmem,"[2]: smithall"); - if (all) - { - GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr; - gerepilemany(av,gptr,3); - } - else x=gerepileupto(av,ZM_copy(x)); + if (DEBUGMEM>1) err(warnmem,"[2]: smithall"); + snf_pile(av, &x,&U,&V); } } } - if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();} - if (all) - { - for (k=1; k<=n; k++) - if (signe(gcoeff(x,k,k))<0) - { mr[k]=lneg((GEN)mr[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); - return gerepile(av,tetpil,z); - } - tetpil=avma; z=cgetg(n+1,t_VEC); j=n; - for (k=n; k; k--) - if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k)); - for ( ; j; j--) z[j]=zero; - return gerepile(av,tetpil,z); + if (DEBUGLEVEL>7) fprintferr("\n"); + for (k=1; k<=n; k++) + if (signe(gcoeff(x,k,k)) < 0) + { + if (V) V[k]=lneg((GEN)V[k]); + coeff(x,k,k) = lnegi(gcoeff(x,k,k)); + } + if (!U && !V) + return gerepileupto(av, mattodiagonal(x)); + + if (U) U = gtrans_i(U); + snf_pile(av, &x,&U,&V); + if (ptU) *ptU = U; + if (ptV) *ptV = V; + 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 @@ -3136,7 +3074,8 @@ 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"); @@ -3254,13 +3193,13 @@ 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); - else x = smithall(x, flag & 1); - if (flag & 4) x = smithclean(x); - return gerepileupto(av, x); + if (flag & 2) x = flag&1 ? gsmith2(x): gsmith(x); + else x = flag&1 ? smith2(x): smith(x); + if (flag & 4) x = gerepileupto(av, smithclean(x)); + return x; } GEN