Annotation of OpenXM_contrib/pari-2.2/src/basemath/alglin2.c, Revision 1.1
1.1 ! noro 1: /* $Id: alglin2.c,v 1.32 2001/10/01 12:11:28 karim Exp $
! 2:
! 3: Copyright (C) 2000 The PARI group.
! 4:
! 5: This file is part of the PARI/GP package.
! 6:
! 7: PARI/GP is free software; you can redistribute it and/or modify it under the
! 8: terms of the GNU General Public License as published by the Free Software
! 9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
! 10: ANY WARRANTY WHATSOEVER.
! 11:
! 12: Check the License for details. You should have received a copy of it, along
! 13: with the package; see the file 'COPYING'. If not, write to the Free Software
! 14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
! 15:
! 16: /********************************************************************/
! 17: /** **/
! 18: /** LINEAR ALGEBRA **/
! 19: /** (second part) **/
! 20: /** **/
! 21: /********************************************************************/
! 22: #include "pari.h"
! 23:
! 24: /*******************************************************************/
! 25: /* */
! 26: /* CHARACTERISTIC POLYNOMIAL */
! 27: /* */
! 28: /*******************************************************************/
! 29:
! 30: GEN
! 31: charpoly0(GEN x, int v, long flag)
! 32: {
! 33: if (v<0) v = 0;
! 34: switch(flag)
! 35: {
! 36: case 0: return caradj(x,v,0);
! 37: case 1: return caract(x,v);
! 38: case 2: return carhess(x,v);
! 39: }
! 40: err(flagerr,"charpoly"); return NULL; /* not reached */
! 41: }
! 42:
! 43: static GEN
! 44: caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,GEN*))
! 45: {
! 46: long av = avma, d;
! 47: GEN p1, p2 = leading_term(p);
! 48:
! 49: if (!signe(x)) return gpowgs(polx[v], degpol(p));
! 50: if (typ(x) != t_POL) x = scalarpol(x,v);
! 51: x = gneg_i(x); x[2] = ladd((GEN)x[2], polx[MAXVARN]);
! 52: p1=subres_f(p, x, NULL);
! 53: if (typ(p1) == t_POL && varn(p1)==MAXVARN)
! 54: setvarn(p1,v);
! 55: else
! 56: p1 = gsubst(p1,MAXVARN,polx[v]);
! 57:
! 58: if (!gcmp1(p2) && (d=degpol(x)) > 0) p1 = gdiv(p1, gpuigs(p2,d));
! 59: return gerepileupto(av,p1);
! 60: }
! 61:
! 62: /* return caract(Mod(x,p)) in variable v */
! 63: GEN
! 64: caract2(GEN p, GEN x, int v)
! 65: {
! 66: return caract2_i(p,x,v, subresall);
! 67: }
! 68: GEN
! 69: caractducos(GEN p, GEN x, int v)
! 70: {
! 71: return caract2_i(p,x,v, (GEN (*)(GEN,GEN,GEN*))resultantducos);
! 72: }
! 73:
! 74: /* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
! 75: static GEN
! 76: easychar(GEN x, int v, GEN *py)
! 77: {
! 78: long l,tetpil,lx;
! 79: GEN p1,p2;
! 80:
! 81: switch(typ(x))
! 82: {
! 83: case t_INT: case t_REAL: case t_INTMOD:
! 84: case t_FRAC: case t_FRACN: case t_PADIC:
! 85: p1=cgetg(4,t_POL);
! 86: p1[1]=evalsigne(1) | evallgef(4) | evalvarn(v);
! 87: p1[2]=lneg(x); p1[3]=un;
! 88: if (py)
! 89: {
! 90: p2=cgetg(2,t_MAT);
! 91: p2[1]=lgetg(2,t_COL);
! 92: coeff(p2,1,1)=lcopy(x);
! 93: *py=p2;
! 94: }
! 95: return p1;
! 96:
! 97: case t_COMPLEX: case t_QUAD:
! 98: if (py) err(typeer,"easychar");
! 99: p1=cgetg(5,t_POL);
! 100: p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v);
! 101: p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma;
! 102: p1[3]=lpile(l,tetpil,gneg(p2));
! 103: p1[4]=un; return p1;
! 104:
! 105: case t_POLMOD:
! 106: if (py) err(typeer,"easychar");
! 107: return caract2((GEN)x[1], (GEN)x[2], v);
! 108:
! 109: case t_MAT:
! 110: lx=lg(x);
! 111: if (lx==1)
! 112: {
! 113: if (py) *py = cgetg(1,t_MAT);
! 114: return polun[v];
! 115: }
! 116: if (lg(x[1]) != lx) break;
! 117: return NULL;
! 118: }
! 119: err(mattype1,"easychar");
! 120: return NULL; /* not reached */
! 121: }
! 122:
! 123: GEN
! 124: caract(GEN x, int v)
! 125: {
! 126: long n,k,l=avma,tetpil;
! 127: GEN p1,p2,p3,p4,p5,p6;
! 128:
! 129: if ((p1 = easychar(x,v,NULL))) return p1;
! 130:
! 131: p1=gzero; p2=gun;
! 132: n=lg(x)-1; if (n&1) p2=gneg_i(p2);
! 133: p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5;
! 134: p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3);
! 135: for (k=0; k<=n; k++)
! 136: {
! 137: p3=det(gsub(gscalmat(stoi(k),n), x));
! 138: p4[1]=lmul(p3,p2); p6[2]=k;
! 139: p1=gadd(p4,p1); p5[2]=(long)p6;
! 140: if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
! 141: }
! 142: p2=mpfact(n); tetpil=avma;
! 143: return gerepile(l,tetpil,gdiv((GEN) p1[1],p2));
! 144: }
! 145:
! 146: GEN
! 147: caradj0(GEN x, long v)
! 148: {
! 149: return caradj(x,v,NULL);
! 150: }
! 151:
! 152: /* Using traces: return the characteristic polynomial of x (in variable v).
! 153: * If py != NULL, the adjoint matrix is put there.
! 154: */
! 155: GEN
! 156: caradj(GEN x, long v, GEN *py)
! 157: {
! 158: long i,j,k,l,av,tetpil;
! 159: GEN p,y,t,*gptr[2];
! 160:
! 161: if ((p = easychar(x,v,py))) return p;
! 162:
! 163: l=lg(x);
! 164: if (l==1) { p=polun[v]; if (py) *py=gcopy(x); return p; }
! 165: if (l==2) { p=gsub(polx[v],gtrace(x)); if (py) *py=idmat(1); return p; }
! 166:
! 167: p=cgetg(l+2,t_POL); p[1]=evalsigne(1) | evallgef(l+2) | evalvarn(v);
! 168: av=avma; t=gtrace(x); tetpil=avma;
! 169: t=gerepile(av,tetpil,gneg(t));
! 170: p[l]=(long)t; p[l+1]=un;
! 171:
! 172: av=avma; y=cgetg(l,t_MAT);
! 173: for (j=1; j<l; j++)
! 174: {
! 175: y[j]=lgetg(l,t_COL);
! 176: for (i=1; i<l; i++)
! 177: coeff(y,i,j) = (i==j) ? ladd(gcoeff(x,i,j),t) : coeff(x,i,j);
! 178: }
! 179:
! 180: for (k=2; k<l-1; k++)
! 181: {
! 182: GEN z=gmul(x,y);
! 183:
! 184: t=gtrace(z); tetpil=avma;
! 185: t=gdivgs(t,-k); y=cgetg(l,t_MAT);
! 186: for (j=1; j<l; j++)
! 187: {
! 188: y[j]=lgetg(l,t_COL);
! 189: for (i=1;i<l;i++)
! 190: coeff(y,i,j) = (i==j)? ladd(gcoeff(z,i,i),t): lcopy(gcoeff(z,i,j));
! 191: }
! 192: gptr[0]=&y; gptr[1]=&t;
! 193: gerepilemanysp(av,tetpil,gptr,2);
! 194: p[l-k+1]=(long)t; av=avma;
! 195: }
! 196:
! 197: t=gzero;
! 198: for (i=1; i<l; i++)
! 199: t=gadd(t,gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
! 200: tetpil=avma; t=gneg(t);
! 201:
! 202: if (py)
! 203: {
! 204: *py = (l&1) ? gneg(y) : gcopy(y);
! 205: gptr[0] = &t; gptr[1] = py;
! 206: gerepilemanysp(av,tetpil,gptr,2);
! 207: p[2]=(long)t;
! 208: }
! 209: else p[2]=lpile(av,tetpil,t);
! 210: i = gvar2(p);
! 211: if (i == v) err(talker,"incorrect variable in caradj");
! 212: if (i < v) p = poleval(p, polx[v]);
! 213: return p;
! 214: }
! 215:
! 216: GEN
! 217: adj(GEN x)
! 218: {
! 219: GEN y;
! 220: caradj(x,MAXVARN,&y); return y;
! 221: }
! 222:
! 223: /*******************************************************************/
! 224: /* */
! 225: /* HESSENBERG FORM */
! 226: /* */
! 227: /*******************************************************************/
! 228: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 229: #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
! 230:
! 231: GEN
! 232: hess(GEN x)
! 233: {
! 234: ulong av = avma;
! 235: long lx=lg(x),m,i,j;
! 236: GEN p,p1,p2;
! 237:
! 238: if (typ(x) != t_MAT) err(mattype1,"hess");
! 239: if (lx==1) return cgetg(1,t_MAT);
! 240: if (lg(x[1])!=lx) err(mattype1,"hess");
! 241: x = dummycopy(x);
! 242:
! 243: for (m=2; m<lx-1; m++)
! 244: for (i=m+1; i<lx; i++)
! 245: {
! 246: p = gcoeff(x,i,m-1);
! 247: if (gcmp0(p)) continue;
! 248:
! 249: for (j=m-1; j<lx; j++) swap(coeff(x,i,j), coeff(x,m,j));
! 250: swap(x[i], x[m]); p = ginv(p);
! 251: for (i=m+1; i<lx; i++)
! 252: {
! 253: p1 = gcoeff(x,i,m-1);
! 254: if (gcmp0(p1)) continue;
! 255:
! 256: p1 = gmul(p1,p); p2 = gneg_i(p1);
! 257: coeff(x,i,m-1) = zero;
! 258: for (j=m; j<lx; j++)
! 259: coeff(x,i,j) = ladd(gcoeff(x,i,j), gmul(p2,gcoeff(x,m,j)));
! 260: for (j=1; j<lx; j++)
! 261: coeff(x,j,m) = ladd(gcoeff(x,j,m), gmul(p1,gcoeff(x,j,i)));
! 262: }
! 263: break;
! 264: }
! 265: return gerepilecopy(av,x);
! 266: }
! 267:
! 268: GEN
! 269: carhess(GEN x, long v)
! 270: {
! 271: long av,tetpil,lx,r,i;
! 272: GEN *y,p1,p2,p3,p4;
! 273:
! 274: if ((p1 = easychar(x,v,NULL))) return p1;
! 275:
! 276: lx=lg(x); av=avma; y = (GEN*) new_chunk(lx);
! 277: y[0]=polun[v]; p1=hess(x); p2=polx[v];
! 278: tetpil=avma;
! 279: for (r=1; r<lx; r++)
! 280: {
! 281: y[r]=gmul(y[r-1], gsub(p2,gcoeff(p1,r,r)));
! 282: p3=gun; p4=gzero;
! 283: for (i=r-1; i; i--)
! 284: {
! 285: p3=gmul(p3,gcoeff(p1,i+1,i));
! 286: p4=gadd(p4,gmul(gmul(p3,gcoeff(p1,i,r)),y[i-1]));
! 287: }
! 288: tetpil=avma; y[r]=gsub(y[r],p4);
! 289: }
! 290: return gerepile(av,tetpil,y[lx-1]);
! 291: }
! 292:
! 293: /*******************************************************************/
! 294: /* */
! 295: /* NORM */
! 296: /* */
! 297: /*******************************************************************/
! 298:
! 299: GEN
! 300: gnorm(GEN x)
! 301: {
! 302: long l,lx,i,tetpil, tx=typ(x);
! 303: GEN p1,p2,y;
! 304:
! 305: switch(tx)
! 306: {
! 307: case t_INT:
! 308: return sqri(x);
! 309:
! 310: case t_REAL:
! 311: return mulrr(x,x);
! 312:
! 313: case t_FRAC: case t_FRACN:
! 314: return gsqr(x);
! 315:
! 316: case t_COMPLEX:
! 317: l=avma; p1=gsqr((GEN) x[1]); p2=gsqr((GEN) x[2]);
! 318: tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
! 319:
! 320: case t_QUAD:
! 321: l=avma; p1=(GEN)x[1];
! 322: p2=gmul((GEN) p1[2], gsqr((GEN) x[3]));
! 323: p1 = gcmp0((GEN) p1[3])? gsqr((GEN) x[2])
! 324: : gmul((GEN) x[2], gadd((GEN) x[2],(GEN) x[3]));
! 325: tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
! 326:
! 327: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
! 328: l=avma; p1=gmul(gconj(x),x); tetpil=avma;
! 329: return gerepile(l,tetpil,greal(p1));
! 330:
! 331: case t_POLMOD:
! 332: p1=(GEN)x[1]; p2=leading_term(p1);
! 333: if (gcmp1(p2) || gcmp0((GEN) x[2])) return subres(p1,(GEN) x[2]);
! 334: l=avma; y=subres(p1,(GEN)x[2]); p1=gpuigs(p2,degpol(x[2]));
! 335: tetpil=avma; return gerepile(l,tetpil,gdiv(y,p1));
! 336:
! 337: case t_VEC: case t_COL: case t_MAT:
! 338: lx=lg(x); y=cgetg(lx,tx);
! 339: for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]);
! 340: return y;
! 341: }
! 342: err(typeer,"gnorm");
! 343: return NULL; /* not reached */
! 344: }
! 345:
! 346: GEN
! 347: gnorml2(GEN x)
! 348: {
! 349: GEN y;
! 350: long i,tx=typ(x),lx,av,lim;
! 351:
! 352: if (! is_matvec_t(tx)) return gnorm(x);
! 353: lx=lg(x); if (lx==1) return gzero;
! 354:
! 355: av=avma; lim = stack_lim(av,1); y = gnorml2((GEN) x[1]);
! 356: for (i=2; i<lx; i++)
! 357: {
! 358: y = gadd(y, gnorml2((GEN) x[i]));
! 359: if (low_stack(lim, stack_lim(av,1)))
! 360: {
! 361: if (DEBUGMEM>1) err(warnmem,"gnorml2");
! 362: y = gerepileupto(av, y);
! 363: }
! 364: }
! 365: return gerepileupto(av,y);
! 366: }
! 367:
! 368: GEN
! 369: QuickNormL2(GEN x, long prec)
! 370: {
! 371: long av = avma;
! 372: GEN y = gmul(x, realun(prec));
! 373: if (typ(x) == t_POL)
! 374: *++y = evaltyp(t_VEC) | evallg(lgef(x)-1);
! 375: return gerepileupto(av, gnorml2(y));
! 376: }
! 377:
! 378: GEN
! 379: gnorml1(GEN x,long prec)
! 380: {
! 381: ulong av = avma;
! 382: long lx,i;
! 383: GEN s;
! 384: switch(typ(x))
! 385: {
! 386: case t_INT: case t_REAL: case t_COMPLEX: case t_FRAC:
! 387: case t_FRACN: case t_QUAD:
! 388: return gabs(x,prec);
! 389:
! 390: case t_POL:
! 391: lx = lgef(x); s = gzero;
! 392: for (i=2; i<lx; i++) s = gadd(s, gabs((GEN)x[i],prec));
! 393: break;
! 394:
! 395: case t_VEC: case t_COL: case t_MAT:
! 396: lx = lg(x); s = gzero;
! 397: for (i=1; i<lx; i++) s = gadd(s, gabs((GEN)x[i],prec));
! 398: break;
! 399:
! 400: default: err(typeer,"gnorml1");
! 401: return NULL; /* not reached */
! 402: }
! 403: return gerepileupto(av, s);
! 404: }
! 405:
! 406: GEN
! 407: QuickNormL1(GEN x,long prec)
! 408: {
! 409: ulong av = avma;
! 410: long lx,i;
! 411: GEN p1,p2,s;
! 412: switch(typ(x))
! 413: {
! 414: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 415: return gabs(x,prec);
! 416:
! 417: case t_INTMOD: case t_PADIC: case t_POLMOD:
! 418: case t_SER: case t_RFRAC: case t_RFRACN:
! 419: return gcopy(x);
! 420:
! 421: case t_COMPLEX:
! 422: p1=gabs((GEN)x[1],prec); p2=gabs((GEN)x[2],prec);
! 423: return gerepileupto(av, gadd(p1,p2));
! 424:
! 425: case t_QUAD:
! 426: p1=gabs((GEN)x[2],prec); p2=gabs((GEN)x[3],prec);
! 427: return gerepileupto(av, gadd(p1,p2));
! 428:
! 429: case t_POL:
! 430: lx=lg(x); s=gzero;
! 431: for (i=2; i<lx; i++) s=gadd(s,QuickNormL1((GEN)x[i],prec));
! 432: return gerepileupto(av, s);
! 433:
! 434: case t_VEC: case t_COL: case t_MAT:
! 435: lx=lg(x); s=gzero;
! 436: for (i=1; i<lx; i++) s=gadd(s,QuickNormL1((GEN)x[i],prec));
! 437: return gerepileupto(av, s);
! 438: }
! 439: err(typeer,"QuickNormL1");
! 440: return NULL; /* not reached */
! 441: }
! 442:
! 443: /*******************************************************************/
! 444: /* */
! 445: /* CONJUGATION */
! 446: /* */
! 447: /*******************************************************************/
! 448:
! 449: GEN
! 450: gconj(GEN x)
! 451: {
! 452: long lx,i,tx=typ(x);
! 453: GEN z;
! 454:
! 455: switch(tx)
! 456: {
! 457: case t_INT: case t_REAL: case t_INTMOD:
! 458: case t_FRAC: case t_FRACN: case t_PADIC:
! 459: return gcopy(x);
! 460:
! 461: case t_COMPLEX:
! 462: z=cgetg(3,t_COMPLEX);
! 463: z[1]=lcopy((GEN) x[1]);
! 464: z[2]=lneg((GEN) x[2]);
! 465: break;
! 466:
! 467: case t_QUAD:
! 468: z=cgetg(4,t_QUAD);
! 469: copyifstack(x[1],z[1]);
! 470: z[2]=gcmp0(gmael(x,1,3))? lcopy((GEN) x[2])
! 471: : ladd((GEN) x[2],(GEN) x[3]);
! 472: z[3]=lneg((GEN) x[3]);
! 473: break;
! 474:
! 475: case t_POL:
! 476: lx=lgef(x); z=cgetg(lx,tx); z[1]=x[1];
! 477: for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 478: break;
! 479:
! 480: case t_SER:
! 481: lx=lg(x); z=cgetg(lx,tx); z[1]=x[1];
! 482: for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 483: break;
! 484:
! 485: case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 486: lx=lg(x); z=cgetg(lx,tx);
! 487: for (i=1; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 488: break;
! 489:
! 490: default:
! 491: err(typeer,"gconj");
! 492: return NULL; /* not reached */
! 493: }
! 494: return z;
! 495: }
! 496:
! 497: GEN
! 498: conjvec(GEN x,long prec)
! 499: {
! 500: long lx,s,av,tetpil,i,tx=typ(x);
! 501: GEN z,y,p1,p2,p;
! 502:
! 503: switch(tx)
! 504: {
! 505: case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN:
! 506: z=cgetg(2,t_COL); z[1]=lcopy(x); break;
! 507:
! 508: case t_COMPLEX: case t_QUAD:
! 509: z=cgetg(3,t_COL); z[1]=lcopy(x); z[2]=lconj(x); break;
! 510:
! 511: case t_VEC: case t_COL:
! 512: lx=lg(x); z=cgetg(lx,t_MAT);
! 513: for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec);
! 514: s=lg(z[1]);
! 515: for (i=2; i<lx; i++)
! 516: if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");
! 517: break;
! 518:
! 519: case t_POLMOD:
! 520: y=(GEN)x[1]; lx=lgef(y);
! 521: if (lx<=3) return cgetg(1,t_COL);
! 522: av=avma; p=NULL;
! 523: for (i=2; i<lx; i++)
! 524: {
! 525: tx=typ(y[i]);
! 526: if (tx==t_INTMOD) p=gmael(y,i,1);
! 527: else
! 528: if (tx != t_INT && ! is_frac_t(tx))
! 529: err(polrationer,"conjvec");
! 530: }
! 531: if (!p)
! 532: {
! 533: p1=roots(y,prec); x = (GEN)x[2]; tetpil=avma;
! 534: z=cgetg(lx-2,t_COL);
! 535: for (i=1; i<=lx-3; i++)
! 536: {
! 537: p2=(GEN)p1[i]; if (gcmp0((GEN)p2[2])) p2 = (GEN)p2[1];
! 538: z[i] = (long)poleval(x, p2);
! 539: }
! 540: return gerepile(av,tetpil,z);
! 541: }
! 542: z=cgetg(lx-2,t_COL); z[1]=lcopy(x);
! 543: for (i=2; i<=lx-3; i++) z[i] = lpui((GEN) z[i-1],p,prec);
! 544: break;
! 545:
! 546: default:
! 547: err(typeer,"conjvec");
! 548: return NULL; /* not reached */
! 549: }
! 550: return z;
! 551: }
! 552:
! 553: /*******************************************************************/
! 554: /* */
! 555: /* TRACES */
! 556: /* */
! 557: /*******************************************************************/
! 558:
! 559: GEN
! 560: assmat(GEN x)
! 561: {
! 562: long lx,i,j;
! 563: GEN y,p1,p2;
! 564:
! 565: if (typ(x)!=t_POL) err(notpoler,"assmat");
! 566: if (gcmp0(x)) err(zeropoler,"assmat");
! 567:
! 568: lx=lgef(x)-2; y=cgetg(lx,t_MAT);
! 569: for (i=1; i<lx-1; i++)
! 570: {
! 571: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 572: for (j=1; j<lx; j++)
! 573: p1[j] = (j==(i+1))? un: zero;
! 574: }
! 575: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 576: if (gcmp1((GEN) x[lx+1]))
! 577: for (j=1; j<lx; j++)
! 578: p1[j] = lneg((GEN)x[j+1]);
! 579: else
! 580: {
! 581: p2 = (GEN)x[lx+1]; gnegz(p2,p2);
! 582: for (j=1; j<lx; j++)
! 583: p1[j] = ldiv((GEN)x[j+1],p2);
! 584: gnegz(p2,p2);
! 585: }
! 586: return y;
! 587: }
! 588:
! 589: GEN
! 590: gtrace(GEN x)
! 591: {
! 592: long i,l,n,tx=typ(x),lx,tetpil;
! 593: GEN y,p1,p2;
! 594:
! 595: switch(tx)
! 596: {
! 597: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 598: return gmul2n(x,1);
! 599:
! 600: case t_COMPLEX:
! 601: return gmul2n((GEN)x[1],1);
! 602:
! 603: case t_QUAD:
! 604: p1=(GEN)x[1];
! 605: if (!gcmp0((GEN) p1[3]))
! 606: {
! 607: l=avma; p2=gmul2n((GEN)x[2],1);
! 608: tetpil=avma; return gerepile(l,tetpil,gadd((GEN)x[3],p2));
! 609: }
! 610: return gmul2n((GEN)x[2],1);
! 611:
! 612: case t_POL:
! 613: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
! 614: for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 615: return y;
! 616:
! 617: case t_SER:
! 618: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 619: for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 620: return y;
! 621:
! 622: case t_POLMOD:
! 623: l=avma; n=(lgef(x[1])-4);
! 624: p1=polsym((GEN)x[1],n); p2=gzero;
! 625: for (i=0; i<=n; i++)
! 626: p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1]));
! 627: return gerepileupto(l,p2);
! 628:
! 629: case t_RFRAC: case t_RFRACN:
! 630: return gadd(x,gconj(x));
! 631:
! 632: case t_VEC: case t_COL:
! 633: lx=lg(x); y=cgetg(lx,tx);
! 634: for (i=1; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 635: return y;
! 636:
! 637: case t_MAT:
! 638: lx=lg(x); if (lx==1) return gzero;/*now lx>=2*/
! 639: if (lx!=lg(x[1])) err(mattype1,"gtrace");
! 640: l=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);
! 641: for (i=2; i<lx-1; i++)
! 642: p1=gadd(p1,gcoeff(x,i,i));
! 643: tetpil=avma; return gerepile(l,tetpil,gadd(p1,gcoeff(x,i,i)));
! 644:
! 645: }
! 646: err(typeer,"gtrace");
! 647: return NULL; /* not reached */
! 648: }
! 649:
! 650: /* Gauss reduction for positive definite matrix a
! 651: * If a is not positive definite:
! 652: * if flag is zero: raise an error
! 653: * else: return NULL.
! 654: */
! 655: GEN
! 656: sqred1intern(GEN a,long flag)
! 657: {
! 658: GEN b,p;
! 659: long i,j,k, n = lg(a);
! 660: ulong av = avma, lim=stack_lim(av,1);
! 661:
! 662: if (typ(a)!=t_MAT) err(typeer,"sqred1");
! 663: if (n == 1) return cgetg(1, t_MAT);
! 664: if (lg(a[1])!=n) err(mattype1,"sqred1");
! 665: b=cgetg(n,t_MAT);
! 666: for (j=1; j<n; j++)
! 667: {
! 668: GEN p1=cgetg(n,t_COL), p2=(GEN)a[j];
! 669:
! 670: b[j]=(long)p1;
! 671: for (i=1; i<=j; i++) p1[i] = p2[i];
! 672: for ( ; i<n ; i++) p1[i] = zero;
! 673: }
! 674: for (k=1; k<n; k++)
! 675: {
! 676: p = gcoeff(b,k,k);
! 677: if (gsigne(p)<=0) /* not positive definite */
! 678: {
! 679: if (flag) { avma=av; return NULL; }
! 680: err(talker,"not a positive definite matrix in sqred1");
! 681: }
! 682: p = ginv(p);
! 683: for (i=k+1; i<n; i++)
! 684: for (j=i; j<n; j++)
! 685: coeff(b,i,j)=lsub(gcoeff(b,i,j),
! 686: gmul(gmul(gcoeff(b,k,i),gcoeff(b,k,j)), p));
! 687: for (j=k+1; j<n; j++)
! 688: coeff(b,k,j)=lmul(gcoeff(b,k,j), p);
! 689: if (low_stack(lim, stack_lim(av,1)))
! 690: {
! 691: if (DEBUGMEM>1) err(warnmem,"sqred1");
! 692: b=gerepilecopy(av,b);
! 693: }
! 694: }
! 695: return gerepilecopy(av,b);
! 696: }
! 697:
! 698: GEN
! 699: sqred1(GEN a)
! 700: {
! 701: return sqred1intern(a,0);
! 702: }
! 703:
! 704: GEN
! 705: sqred3(GEN a)
! 706: {
! 707: ulong av = avma, lim = stack_lim(av,1);
! 708: long i,j,k,l, n = lg(a);
! 709: GEN p1,b;
! 710:
! 711: if (typ(a)!=t_MAT) err(typeer,"sqred3");
! 712: if (lg(a[1])!=n) err(mattype1,"sqred3");
! 713: av=avma; b=cgetg(n,t_MAT);
! 714: for (j=1; j<n; j++)
! 715: {
! 716: p1=cgetg(n,t_COL); b[j]=(long)p1;
! 717: for (i=1; i<n; i++) p1[i]=zero;
! 718: }
! 719: for (i=1; i<n; i++)
! 720: {
! 721: for (k=1; k<i; k++)
! 722: {
! 723: p1=gzero;
! 724: for (l=1; l<k; l++)
! 725: p1=gadd(p1, gmul(gmul(gcoeff(b,l,l),gcoeff(b,k,l)), gcoeff(b,i,l)));
! 726: coeff(b,i,k)=ldiv(gsub(gcoeff(a,i,k),p1),gcoeff(b,k,k));
! 727: }
! 728: p1=gzero;
! 729: for (l=1; l<i; l++)
! 730: p1=gadd(p1, gmul(gcoeff(b,l,l), gsqr(gcoeff(b,i,l))));
! 731: coeff(b,i,k)=lsub(gcoeff(a,i,i),p1);
! 732: if (low_stack(lim, stack_lim(av,1)))
! 733: {
! 734: if (DEBUGMEM>1) err(warnmem,"sqred3");
! 735: b=gerepilecopy(av,b);
! 736: }
! 737: }
! 738: return gerepilecopy(av,b);
! 739: }
! 740:
! 741: /* Gauss reduction (arbitrary symetric matrix, only the part above the
! 742: * diagonal is considered). If no_signature is zero, return only the
! 743: * signature
! 744: */
! 745: static GEN
! 746: sqred2(GEN a, long no_signature)
! 747: {
! 748: GEN r,p,mun;
! 749: ulong av,av1,lim;
! 750: long n,i,j,k,l,sp,sn,t;
! 751:
! 752: if (typ(a)!=t_MAT) err(typeer,"sqred2");
! 753: n=lg(a); if (lg(a[1]) != n) err(mattype1,"sqred2");
! 754:
! 755: av=avma; mun=negi(gun); r=new_chunk(n);
! 756: for (i=1; i<n; i++) r[i]=1;
! 757: av1=avma; lim=stack_lim(av1,1); a=dummycopy(a);
! 758: n--; t=n; sp=sn=0;
! 759: while (t)
! 760: {
! 761: k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++;
! 762: if (k<=n)
! 763: {
! 764: p=gcoeff(a,k,k); if (gsigne(p)>0) sp++; else sn++;
! 765: r[k]=0; t--;
! 766: for (j=1; j<=n; j++)
! 767: coeff(a,k,j)=r[j] ? ldiv(gcoeff(a,k,j),p) : zero;
! 768:
! 769: for (i=1; i<=n; i++) if (r[i])
! 770: for (j=1; j<=n; j++)
! 771: coeff(a,i,j) = r[j] ? lsub(gcoeff(a,i,j),
! 772: gmul(gmul(gcoeff(a,k,i),gcoeff(a,k,j)),p))
! 773: : zero;
! 774: coeff(a,k,k)=(long)p;
! 775: }
! 776: else
! 777: {
! 778: for (k=1; k<=n; k++) if (r[k])
! 779: {
! 780: l=k+1; while (l<=n && (gcmp0(gcoeff(a,k,l)) || !r[l])) l++;
! 781: if (l<=n)
! 782: {
! 783: p=gcoeff(a,k,l); r[k]=r[l]=0; sp++; sn++; t-=2;
! 784: for (i=1; i<=n; i++) if (r[i])
! 785: {
! 786: for (j=1; j<=n; j++)
! 787: coeff(a,i,j) =
! 788: r[j]? lsub(gcoeff(a,i,j),
! 789: gdiv(gadd(gmul(gcoeff(a,k,i),gcoeff(a,l,j)),
! 790: gmul(gcoeff(a,k,j),gcoeff(a,l,i))),
! 791: p))
! 792: : zero;
! 793: coeff(a,k,i)=ldiv(gadd(gcoeff(a,k,i),gcoeff(a,l,i)),p);
! 794: coeff(a,l,i)=ldiv(gsub(gcoeff(a,k,i),gcoeff(a,l,i)),p);
! 795: }
! 796: coeff(a,k,l)=un; coeff(a,l,k)=(long)mun;
! 797: coeff(a,k,k)=lmul2n(p,-1); coeff(a,l,l)=lneg(gcoeff(a,k,k));
! 798: if (low_stack(lim, stack_lim(av1,1)))
! 799: {
! 800: if(DEBUGMEM>1) err(warnmem,"sqred2");
! 801: a=gerepilecopy(av1,a);
! 802: }
! 803: break;
! 804: }
! 805: }
! 806: if (k>n) break;
! 807: }
! 808: }
! 809: if (no_signature) return gerepilecopy(av,a);
! 810: avma=av;
! 811: a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a;
! 812: }
! 813:
! 814: GEN
! 815: sqred(GEN a) { return sqred2(a,1); }
! 816:
! 817: GEN
! 818: signat(GEN a) { return sqred2(a,0); }
! 819:
! 820: /* Diagonalization of a REAL symetric matrix. Return a 2-component vector:
! 821: * 1st comp = vector of eigenvalues
! 822: * 2nd comp = matrix of eigenvectors
! 823: */
! 824: GEN
! 825: jacobi(GEN a, long prec)
! 826: {
! 827: long de,e,e1,e2,l,n,i,j,p,q,av1,av2;
! 828: GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1;
! 829:
! 830: if (typ(a)!=t_MAT) err(mattype1,"jacobi");
! 831: ja=cgetg(3,t_VEC); l=lg(a); n=l-1;
! 832: e1=HIGHEXPOBIT-1;
! 833: lambda=cgetg(l,t_COL); ja[1]=(long)lambda;
! 834: for (j=1; j<=n; j++)
! 835: {
! 836: lambda[j]=lgetr(prec);
! 837: gaffect(gcoeff(a,j,j), (GEN)lambda[j]);
! 838: e=expo(lambda[j]); if (e<e1) e1=e;
! 839: }
! 840: r=cgetg(l,t_MAT); ja[2]=(long)r;
! 841: for (j=1; j<=n; j++)
! 842: {
! 843: r[j]=lgetg(l,t_COL);
! 844: for (i=1; i<=n; i++)
! 845: affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec)));
! 846: }
! 847: av1=avma;
! 848:
! 849: e2=-HIGHEXPOBIT;p=q=1;
! 850: c=cgetg(l,t_MAT);
! 851: for (j=1; j<=n; j++)
! 852: {
! 853: c[j]=lgetg(j,t_COL);
! 854: for (i=1; i<j; i++)
! 855: {
! 856: gaffect(gcoeff(a,i,j), (GEN)(coeff(c,i,j)=lgetr(prec)));
! 857: e=expo(gcoeff(c,i,j)); if (e>e2) { e2=e; p=i; q=j; }
! 858: }
! 859: }
! 860: a=c; unr = realun(prec);
! 861: de=bit_accuracy(prec);
! 862:
! 863: /* e1 = min des expo des coeff diagonaux
! 864: * e2 = max des expo des coeff extra-diagonaux
! 865: * Test d'arret: e2 < e1-precision
! 866: */
! 867: while (e1-e2<de)
! 868: {
! 869: /*calcul de la rotation associee dans le plan
! 870: des p et q-iemes vecteurs de base */
! 871: av2=avma;
! 872: x=divrr(subrr((GEN)lambda[q],(GEN)lambda[p]),shiftr(gcoeff(a,p,q),1));
! 873: y=mpsqrt(addrr(unr,mulrr(x,x)));
! 874: t=(gsigne(x)>0)? divrr(unr,addrr(x,y)) : divrr(unr,subrr(x,y));
! 875: c=divrr(unr,mpsqrt(addrr(unr,mulrr(t,t))));
! 876: s=mulrr(t,c); u=divrr(s,addrr(unr,c));
! 877:
! 878: /* Recalcul des transformees successives de a et de la matrice
! 879: cumulee (r) des rotations : */
! 880:
! 881: for (i=1; i<p; i++)
! 882: {
! 883: x=gcoeff(a,i,p); y=gcoeff(a,i,q);
! 884: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 885: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 886: affrr(x1,gcoeff(a,i,p)); affrr(y1,gcoeff(a,i,q));
! 887: }
! 888: for (i=p+1; i<q; i++)
! 889: {
! 890: x=gcoeff(a,p,i); y=gcoeff(a,i,q);
! 891: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 892: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 893: affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,i,q));
! 894: }
! 895: for (i=q+1; i<=n; i++)
! 896: {
! 897: x=gcoeff(a,p,i); y=gcoeff(a,q,i);
! 898: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 899: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 900: affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,q,i));
! 901: }
! 902: x=(GEN)lambda[p]; y=gcoeff(a,p,q); subrrz(x,mulrr(t,y),(GEN)lambda[p]);
! 903: x=y; y=(GEN)lambda[q]; addrrz(y,mulrr(t,x),y);
! 904: setexpo(x,expo(x)-de-1);
! 905:
! 906: for (i=1; i<=n; i++)
! 907: {
! 908: x=gcoeff(r,i,p); y=gcoeff(r,i,q);
! 909: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 910: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 911: affrr(x1,gcoeff(r,i,p)); affrr(y1,gcoeff(r,i,q));
! 912: }
! 913:
! 914: e2=expo(gcoeff(a,1,2)); p=1; q=2;
! 915: for (j=1; j<=n; j++)
! 916: {
! 917: for (i=1; i<j; i++)
! 918: if ((e=expo(gcoeff(a,i,j))) > e2) { e2=e; p=i; q=j; }
! 919: for (i=j+1; i<=n; i++)
! 920: if ((e=expo(gcoeff(a,j,i))) > e2) { e2=e; p=j; q=i; }
! 921: }
! 922: avma=av2;
! 923: }
! 924: avma=av1; return ja;
! 925: }
! 926:
! 927: /*************************************************************************/
! 928: /** **/
! 929: /** MATRICE RATIONNELLE --> ENTIERE **/
! 930: /** **/
! 931: /*************************************************************************/
! 932:
! 933: GEN
! 934: matrixqz0(GEN x,GEN p)
! 935: {
! 936: if (typ(p)!=t_INT) err(typeer,"matrixqz0");
! 937: if (signe(p)>=0) return matrixqz(x,p);
! 938: if (!cmpsi(-1,p)) return matrixqz2(x);
! 939: if (!cmpsi(-2,p)) return matrixqz3(x);
! 940: err(flagerr,"matrixqz"); return NULL; /* not reached */
! 941: }
! 942:
! 943: GEN
! 944: matrixqz(GEN x, GEN p)
! 945: {
! 946: ulong av = avma, av1, lim;
! 947: long i,j,j1,m,n,t,nfact;
! 948: GEN p1,p2,p3,unmodp;
! 949:
! 950: if (typ(x)!=t_MAT) err(typeer,"matrixqz");
! 951: n=lg(x)-1; if (!n) return gcopy(x);
! 952: m=lg(x[1])-1;
! 953: if (n > m) err(talker,"more rows than columns in matrixqz");
! 954: if (n==m)
! 955: {
! 956: p1=det(x);
! 957: if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz");
! 958: avma=av; return idmat(n);
! 959: }
! 960: p1=cgetg(n+1,t_MAT);
! 961: for (j=1; j<=n; j++)
! 962: {
! 963: p2=gun; p3=(GEN)x[j];
! 964: for (i=1; i<=m; i++)
! 965: {
! 966: t=typ(p3[i]);
! 967: if (t != t_INT && !is_frac_t(t))
! 968: err(talker,"not a rational or integral matrix in matrixqz");
! 969: p2=ggcd(p2,(GEN)p3[i]);
! 970: }
! 971: p1[j]=ldiv(p3,p2);
! 972: }
! 973: x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un;
! 974:
! 975: if (gcmp0(p))
! 976: {
! 977: p1=cgetg(n+1,t_MAT); p2=gtrans(x);
! 978: for (j=1; j<=n; j++) p1[j]=p2[j];
! 979: p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1));
! 980: if (!signe(p3))
! 981: err(impl,"matrixqz when the first 2 dets are zero");
! 982: if (gcmp1(p3)) return gerepilecopy(av,x);
! 983:
! 984: p3=factor(p3); p1=(GEN)p3[1]; nfact=lg(p1)-1;
! 985: }
! 986: else
! 987: {
! 988: p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1;
! 989: }
! 990: av1=avma; lim=stack_lim(av1,1);
! 991: for (i=1; i<=nfact; i++)
! 992: {
! 993: p=(GEN)p1[i]; unmodp[1]=(long)p;
! 994: for(;;)
! 995: {
! 996: p2=ker(gmul(unmodp,x));
! 997: if (lg(p2)==1) break;
! 998:
! 999: p2=centerlift(p2); p3=gdiv(gmul(x,p2),p);
! 1000: for (j=1; j<lg(p2); j++)
! 1001: {
! 1002: j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;
! 1003: x[j1] = p3[j];
! 1004: }
! 1005: if (low_stack(lim, stack_lim(av1,1)))
! 1006: {
! 1007: if (DEBUGMEM>1) err(warnmem,"matrixqz");
! 1008: x=gerepilecopy(av1,x);
! 1009: }
! 1010: }
! 1011: }
! 1012: return gerepilecopy(av,x);
! 1013: }
! 1014:
! 1015: static GEN
! 1016: matrixqz_aux(GEN x, long m, long n)
! 1017: {
! 1018: ulong av = avma, lim = stack_lim(av,1);
! 1019: long i,j,k,in[2];
! 1020: GEN p1;
! 1021:
! 1022: for (i=1; i<=m; i++)
! 1023: {
! 1024: for(;;)
! 1025: {
! 1026: long fl=0;
! 1027:
! 1028: for (j=1; j<=n; j++)
! 1029: if (!gcmp0(gcoeff(x,i,j)))
! 1030: { in[fl]=j; fl++; if (fl==2) break; }
! 1031: if (j>n) break;
! 1032:
! 1033: j=(gcmp(gabs(gcoeff(x,i,in[0]),DEFAULTPREC),
! 1034: gabs(gcoeff(x,i,in[1]),DEFAULTPREC)) > 0)? in[1]: in[0];
! 1035: p1=gcoeff(x,i,j);
! 1036: for (k=1; k<=n; k++)
! 1037: if (k!=j)
! 1038: x[k]=lsub((GEN)x[k],gmul(ground(gdiv(gcoeff(x,i,k),p1)),(GEN)x[j]));
! 1039: }
! 1040:
! 1041: j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++;
! 1042: if (j<=n)
! 1043: {
! 1044: p1=denom(gcoeff(x,i,j));
! 1045: if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]);
! 1046: }
! 1047: if (low_stack(lim, stack_lim(av,1)))
! 1048: {
! 1049: if(DEBUGMEM>1) err(warnmem,"matrixqz_aux");
! 1050: x=gerepilecopy(av,x);
! 1051: }
! 1052: }
! 1053: return hnf(x);
! 1054: }
! 1055:
! 1056: GEN
! 1057: matrixqz2(GEN x)
! 1058: {
! 1059: long av = avma,m,n;
! 1060:
! 1061: if (typ(x)!=t_MAT) err(typeer,"matrixqz2");
! 1062: n=lg(x)-1; if (!n) return gcopy(x);
! 1063: m=lg(x[1])-1; x=dummycopy(x);
! 1064: return gerepileupto(av, matrixqz_aux(x,m,n));
! 1065: }
! 1066:
! 1067: GEN
! 1068: matrixqz3(GEN x)
! 1069: {
! 1070: long av=avma,av1,j,j1,k,m,n,lim;
! 1071: GEN c;
! 1072:
! 1073: if (typ(x)!=t_MAT) err(typeer,"matrixqz3");
! 1074: n=lg(x)-1; if (!n) return gcopy(x);
! 1075: m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1);
! 1076: for (j=1; j<=n; j++) c[j]=0;
! 1077: av1=avma; lim=stack_lim(av1,1);
! 1078: for (k=1; k<=m; k++)
! 1079: {
! 1080: j=1; while (j<=n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
! 1081: if (j<=n)
! 1082: {
! 1083: c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));
! 1084: for (j1=1; j1<=n; j1++)
! 1085: if (j1!=j) x[j1]=lsub((GEN)x[j1],gmul(gcoeff(x,k,j1),(GEN)x[j]));
! 1086: if (low_stack(lim, stack_lim(av1,1)))
! 1087: {
! 1088: if(DEBUGMEM>1) err(warnmem,"matrixqz3");
! 1089: x=gerepilecopy(av1,x);
! 1090: }
! 1091: }
! 1092: }
! 1093: return gerepileupto(av, matrixqz_aux(x,m,n));
! 1094: }
! 1095:
! 1096: GEN
! 1097: intersect(GEN x, GEN y)
! 1098: {
! 1099: long av,tetpil,j, lx = lg(x);
! 1100: GEN z;
! 1101:
! 1102: if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect");
! 1103: if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
! 1104:
! 1105: av=avma; z=ker(concatsp(x,y));
! 1106: for (j=lg(z)-1; j; j--) setlg(z[j],lx);
! 1107: tetpil=avma; return gerepile(av,tetpil,gmul(x,z));
! 1108: }
! 1109:
! 1110: /**************************************************************/
! 1111: /** **/
! 1112: /** HERMITE NORMAL FORM REDUCTION **/
! 1113: /** **/
! 1114: /**************************************************************/
! 1115:
! 1116: GEN
! 1117: mathnf0(GEN x, long flag)
! 1118: {
! 1119: switch(flag)
! 1120: {
! 1121: case 0: return hnf(x);
! 1122: case 1: return hnfall(x);
! 1123: case 3: return hnfperm(x);
! 1124: case 4: return hnflll(x);
! 1125: default: err(flagerr,"mathnf");
! 1126: }
! 1127: return NULL; /* not reached */
! 1128: }
! 1129:
! 1130: static GEN
! 1131: init_hnf(GEN x, GEN *denx, long *co, long *li, long *av)
! 1132: {
! 1133: if (typ(x) != t_MAT) err(typeer,"mathnf");
! 1134: *co=lg(x); if (*co==1) return NULL; /* empty matrix */
! 1135: *li=lg(x[1]); *denx=denom(x); *av=avma;
! 1136:
! 1137: if (gcmp1(*denx)) /* no denominator */
! 1138: { *denx = NULL; return dummycopy(x); }
! 1139: return gmul(x,*denx);
! 1140: }
! 1141:
! 1142: GEN
! 1143: ZV_copy(GEN x)
! 1144: {
! 1145: long i, lx = lg(x);
! 1146: GEN y = cgetg(lx, t_COL);
! 1147: for (i=1; i<lx; i++) y[i] = signe(x[i])? licopy((GEN)x[i]): zero;
! 1148: return y;
! 1149: }
! 1150:
! 1151: GEN
! 1152: ZM_copy(GEN x)
! 1153: {
! 1154: long i, lx = lg(x);
! 1155: GEN y = cgetg(lx, t_MAT);
! 1156:
! 1157: for (i=1; i<lx; i++)
! 1158: y[i] = (long)ZV_copy((GEN)x[i]);
! 1159: return y;
! 1160: }
! 1161:
! 1162: /* return c * X. Not memory clean if c = 1 */
! 1163: GEN
! 1164: ZV_Z_mul(GEN c, GEN X)
! 1165: {
! 1166: long i,m = lg(X);
! 1167: GEN A = new_chunk(m);
! 1168: if (is_pm1(c))
! 1169: {
! 1170: if (signe(c) > 0)
! 1171: { for (i=1; i<m; i++) A[i] = X[i]; }
! 1172: else
! 1173: { for (i=1; i<m; i++) A[i] = lnegi((GEN)X[i]); }
! 1174: }
! 1175: else /* c = 0 should be exceedingly rare */
! 1176: { for (i=1; i<m; i++) A[i] = lmulii(c,(GEN)X[i]); }
! 1177: A[0] = X[0]; return A;
! 1178: }
! 1179:
! 1180: /* X,Y columns; u,v scalars; everybody is integral. return A = u*X + v*Y
! 1181: * Not memory clean if (u,v) = (1,0) or (0,1) */
! 1182: GEN
! 1183: ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)
! 1184: {
! 1185: long av,i,lx,m;
! 1186: GEN p1,p2,A;
! 1187:
! 1188: if (!signe(u)) return ZV_Z_mul(v,Y);
! 1189: if (!signe(v)) return ZV_Z_mul(u,X);
! 1190: lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4;
! 1191: if (is_pm1(v)) { gswap(u,v); gswap(X,Y); }
! 1192: if (is_pm1(u))
! 1193: {
! 1194: if (signe(u) > 0)
! 1195: {
! 1196: for (i=1; i<lx; i++)
! 1197: {
! 1198: p1=(GEN)X[i]; p2=(GEN)Y[i];
! 1199: if (!signe(p1)) A[i] = lmulii(v,p2);
! 1200: else if (!signe(p2)) A[i] = licopy(p1);
! 1201: else
! 1202: {
! 1203: av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /* HACK */
! 1204: p2 = mulii(v,p2);
! 1205: avma = av; A[i] = laddii(p1,p2);
! 1206: }
! 1207: }
! 1208: }
! 1209: else
! 1210: {
! 1211: for (i=1; i<lx; i++)
! 1212: {
! 1213: p1=(GEN)X[i]; p2=(GEN)Y[i];
! 1214: if (!signe(p1)) A[i] = lmulii(v,p2);
! 1215: else if (!signe(p2)) A[i] = lnegi(p1);
! 1216: else
! 1217: {
! 1218: av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /* HACK */
! 1219: p2 = mulii(v,p2);
! 1220: avma = av; A[i] = lsubii(p2,p1);
! 1221: }
! 1222: }
! 1223: }
! 1224: }
! 1225: else
! 1226: for (i=1; i<lx; i++)
! 1227: {
! 1228: p1=(GEN)X[i]; p2=(GEN)Y[i];
! 1229: if (!signe(p1)) A[i] = lmulii(v,p2);
! 1230: else if (!signe(p2)) A[i] = lmulii(u,p1);
! 1231: else
! 1232: {
! 1233: av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /* HACK */
! 1234: p1 = mulii(u,p1);
! 1235: p2 = mulii(v,p2);
! 1236: avma = av; A[i] = laddii(p1,p2);
! 1237: }
! 1238: }
! 1239: return A;
! 1240: }
! 1241:
! 1242: /* x = [A,U], nbcol(A) = nbcol(U), A integral. Return [AV, UV], with AV HNF */
! 1243: GEN
! 1244: hnf_special(GEN x, long remove)
! 1245: {
! 1246: long av0, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
! 1247: GEN p1,u,v,d,denx,a,b, x2,res;
! 1248:
! 1249: if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special");
! 1250: res = cgetg(3,t_VEC);
! 1251: av0 = avma;
! 1252: x2 = (GEN)x[2];
! 1253: x = (GEN)x[1];
! 1254: x = init_hnf(x,&denx,&co,&li,&av);
! 1255: if (!x) return cgetg(1,t_MAT);
! 1256:
! 1257: lim = stack_lim(av,1);
! 1258: def=co-1; ldef=(li>co)? li-co: 0;
! 1259: if (lg(x2) != co) err(talker,"incompatible matrices in hnf_special");
! 1260: x2 = dummycopy(x2);
! 1261: for (i=li-1; i>ldef; i--)
! 1262: {
! 1263: for (j=def-1; j; j--)
! 1264: {
! 1265: a = gcoeff(x,i,j);
! 1266: if (!signe(a)) continue;
! 1267:
! 1268: k = (j==1)? def: j-1;
! 1269: b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
! 1270: if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
! 1271: p1 = (GEN)x[j]; b = negi(b);
! 1272: x[j] = (long)ZV_lincomb(a,b, (GEN)x[k], p1);
! 1273: x[k] = (long)ZV_lincomb(u,v, p1, (GEN)x[k]);
! 1274: p1 = (GEN)x2[j];
! 1275: x2[j]= ladd(gmul(a, (GEN)x2[k]), gmul(b,p1));
! 1276: x2[k] = ladd(gmul(u,p1), gmul(v, (GEN)x2[k]));
! 1277: if (low_stack(lim, stack_lim(av,1)))
! 1278: {
! 1279: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1280: if (DEBUGMEM>1) err(warnmem,"hnf_special[1]. i=%ld",i);
! 1281: gerepilemany(av,gptr,2);
! 1282: }
! 1283: }
! 1284: p1 = gcoeff(x,i,def); s = signe(p1);
! 1285: if (s)
! 1286: {
! 1287: if (s < 0)
! 1288: {
! 1289: x[def] = lneg((GEN)x[def]); p1 = gcoeff(x,i,def);
! 1290: x2[def]= lneg((GEN)x2[def]);
! 1291: }
! 1292: for (j=def+1; j<co; j++)
! 1293: {
! 1294: b = negi(gdivent(gcoeff(x,i,j),p1));
! 1295: x[j] = (long)ZV_lincomb(gun,b, (GEN)x[j], (GEN)x[def]);
! 1296: x2[j]= ladd((GEN)x2[j], gmul(b, (GEN)x2[def]));
! 1297: }
! 1298: def--;
! 1299: }
! 1300: else
! 1301: if (ldef && i==ldef+1) ldef--;
! 1302: if (low_stack(lim, stack_lim(av,1)))
! 1303: {
! 1304: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1305: if (DEBUGMEM>1) err(warnmem,"hnf_special[2]. i=%ld",i);
! 1306: gerepilemany(av,gptr,2);
! 1307: }
! 1308: }
! 1309: if (remove)
! 1310: { /* remove null columns */
! 1311: for (i=1,j=1; j<co; j++)
! 1312: if (!gcmp0((GEN)x[j]))
! 1313: {
! 1314: x[i] = x[j];
! 1315: x2[i] = x2[j]; i++;
! 1316: }
! 1317: setlg(x,i);
! 1318: setlg(x2,i);
! 1319: }
! 1320: tetpil=avma;
! 1321: x = denx? gdiv(x,denx): ZM_copy(x);
! 1322: x2 = gcopy(x2);
! 1323: {
! 1324: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1325: gerepilemanysp(av0,tetpil,gptr,2);
! 1326: }
! 1327: res[1] = (long)x;
! 1328: res[2] = (long)x2;
! 1329: return res;
! 1330: }
! 1331:
! 1332: /*******************************************************************/
! 1333: /* */
! 1334: /* SPECIAL HNF (FOR INTERNAL USE !!!) */
! 1335: /* */
! 1336: /*******************************************************************/
! 1337: extern GEN imagecomplspec(GEN x, long *nlze);
! 1338:
! 1339: static int
! 1340: count(long **mat, long row, long len, long *firstnonzero)
! 1341: {
! 1342: int j, n=0;
! 1343:
! 1344: for (j=1; j<=len; j++)
! 1345: {
! 1346: long p = mat[j][row];
! 1347: if (p)
! 1348: {
! 1349: if (labs(p)!=1) return -1;
! 1350: n++; *firstnonzero=j;
! 1351: }
! 1352: }
! 1353: return n;
! 1354: }
! 1355:
! 1356: static int
! 1357: count2(long **mat, long row, long len)
! 1358: {
! 1359: int j;
! 1360: for (j=len; j; j--)
! 1361: if (labs(mat[j][row]) == 1) return j;
! 1362: return 0;
! 1363: }
! 1364:
! 1365: static GEN
! 1366: hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
! 1367: {
! 1368: GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
! 1369: GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
! 1370: long av,i,j,k,s,i1,j1,lim,zc;
! 1371: long co = lg(C);
! 1372: long col = lg(matgen)-1;
! 1373: long lnz, nlze, lig;
! 1374:
! 1375: if (col == 0) return matgen;
! 1376: lnz = lg(matgen[1])-1; /* was called lnz-1 - nr in hnfspec */
! 1377: nlze = lg(dep[1])-1;
! 1378: lig = nlze + lnz;
! 1379: if (DEBUGLEVEL>5)
! 1380: {
! 1381: fprintferr("Entering hnffinal:\n");
! 1382: if (DEBUGLEVEL>6)
! 1383: {
! 1384: if (nlze) fprintferr("dep = %Z\n",dep);
! 1385: fprintferr("mit = %Z\n",matgen);
! 1386: fprintferr("B = %Z\n",B);
! 1387: }
! 1388: }
! 1389: p1 = hnflll(matgen);
! 1390: H = (GEN)p1[1]; /* lnz x lnz */
! 1391: U = (GEN)p1[2]; /* col x col */
! 1392: /* Only keep the part above the H (above the 0s is 0 since the dep rows
! 1393: * are dependant from the ones in matgen) */
! 1394: zc = col - lnz; /* # of 0 columns, correspond to units */
! 1395: if (nlze) { dep = gmul(dep,U); dep += zc; }
! 1396:
! 1397: diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
! 1398:
! 1399: av = avma; lim = stack_lim(av,1);
! 1400: Cnew = cgetg(co,t_MAT);
! 1401: setlg(C, col+1); p1 = gmul(C,U);
! 1402: for (j=1; j<=col; j++) Cnew[j] = p1[j];
! 1403: for ( ; j<co ; j++) Cnew[j] = C[j];
! 1404: if (DEBUGLEVEL>5) fprintferr(" hnflll done\n");
! 1405:
! 1406: /* Clean up B using new H */
! 1407: for (s=0,i=lnz; i; i--)
! 1408: {
! 1409: GEN h = gcoeff(H,i,i);
! 1410: if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; }
! 1411: for (j=col+1; j<co; j++)
! 1412: {
! 1413: GEN z = (GEN)B[j-col];
! 1414: p1 = (GEN)z[i+nlze]; if (h) p1 = gdivent(p1,h);
! 1415: for (k=1; k<=nlze; k++)
! 1416: z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(dep,k,i)));
! 1417: for ( ; k<=lig; k++)
! 1418: z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(H,k-nlze,i)));
! 1419: Cnew[j] = lsub((GEN)Cnew[j], gmul(p1, (GEN)Cnew[i+zc]));
! 1420: }
! 1421: if (low_stack(lim, stack_lim(av,1)))
! 1422: {
! 1423: GEN *gptr[2]; gptr[0]=&Cnew; gptr[1]=&B;
! 1424: if(DEBUGMEM>1) err(warnmem,"hnffinal, i = %ld",i);
! 1425: gerepilemany(av,gptr,2);
! 1426: }
! 1427: }
! 1428: p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
! 1429: for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
! 1430: if (diagH1[i])
! 1431: p1[++j1] = p2[i];
! 1432: else
! 1433: p2[++i1] = p2[i];
! 1434: for (i=i1+1; i<=lnz; i++) p2[i] = p1[i];
! 1435: if (DEBUGLEVEL>5) fprintferr(" first pass in hnffinal done\n");
! 1436:
! 1437: /* s = # extra redundant generators taken from H
! 1438: * zc col-s co zc = col lnz
! 1439: * [ 0 |dep | ] i = lnze + lnz - s = lig - s
! 1440: * nlze [--------| B' ]
! 1441: * [ 0 | H' | ] H' = H minus the s rows with a 1 on diagonal
! 1442: * i [--------|-----] lig-s (= "1-rows")
! 1443: * [ 0 | Id ]
! 1444: * [ | ] li */
! 1445: lig -= s; col -= s; lnz -= s;
! 1446: Hnew = cgetg(lnz+1,t_MAT);
! 1447: depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
! 1448: Bnew = cgetg(co-col,t_MAT);
! 1449: C = dummycopy(Cnew);
! 1450: for (j=1,i1=j1=0; j<=lnz+s; j++)
! 1451: {
! 1452: GEN z = (GEN)H[j];
! 1453: if (diagH1[j])
! 1454: { /* hit exactly s times */
! 1455: i1++; p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;
! 1456: C[i1+col] = Cnew[j+zc];
! 1457: for (i=1; i<=nlze; i++) p1[i] = coeff(dep,i,j);
! 1458: p1 += nlze;
! 1459: }
! 1460: else
! 1461: {
! 1462: j1++; p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;
! 1463: C[j1+zc] = Cnew[j+zc];
! 1464: if (nlze) depnew[j1] = dep[j];
! 1465: }
! 1466: for (i=k=1; k<=lnz; i++)
! 1467: if (!diagH1[i]) p1[k++] = z[i];
! 1468: }
! 1469: for (j=s+1; j<co-col; j++)
! 1470: {
! 1471: GEN z = (GEN)B[j-s];
! 1472: p1 = cgetg(lig+1,t_COL); Bnew[j] = (long)p1;
! 1473: for (i=1; i<=nlze; i++) p1[i] = z[i];
! 1474: z += nlze; p1 += nlze;
! 1475: for (i=k=1; k<=lnz; i++)
! 1476: if (!diagH1[i]) p1[k++] = z[i];
! 1477: }
! 1478: if (DEBUGLEVEL>5)
! 1479: {
! 1480: fprintferr("Leaving hnffinal\n");
! 1481: if (DEBUGLEVEL>6)
! 1482: {
! 1483: if (nlze) fprintferr("dep = %Z\n",depnew);
! 1484: fprintferr("mit = %Z\nB = %Z\nC = %Z\n", Hnew, Bnew, C);
! 1485: }
! 1486: }
! 1487: if (nlze) *ptdep = depnew;
! 1488: *ptC = C;
! 1489: *ptB = Bnew; return Hnew;
! 1490: }
! 1491:
! 1492: /* for debugging */
! 1493: static void
! 1494: p_mat(long **mat, long *perm, long k0)
! 1495: {
! 1496: long av=avma, i,j;
! 1497: GEN p1, matj, matgen;
! 1498: long co = lg(mat);
! 1499: long li = lg(perm);
! 1500:
! 1501: fprintferr("Permutation: %Z\n",perm);
! 1502: matgen = cgetg(co,t_MAT);
! 1503: for (j=1; j<co; j++)
! 1504: {
! 1505: p1 = cgetg(li-k0,t_COL); matgen[j]=(long)p1;
! 1506: p1 -= k0; matj = mat[j];
! 1507: for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);
! 1508: }
! 1509: if (DEBUGLEVEL > 6) fprintferr("matgen = %Z\n",matgen);
! 1510: avma=av;
! 1511: }
! 1512:
! 1513: static GEN
! 1514: col_dup(long n, GEN col)
! 1515: {
! 1516: GEN c = new_chunk(n+1);
! 1517: memcpy(c,col,(n+1)*sizeof(long));
! 1518: return c;
! 1519: }
! 1520:
! 1521: /* HNF reduce a relation matrix (column operations + row permutation)
! 1522: ** Input:
! 1523: ** mat = (li-1) x (co-1) matrix of long
! 1524: ** C = r x (co-1) matrix of GEN
! 1525: ** perm= permutation vector (length li-1), indexing the rows of mat: easier
! 1526: ** to maintain perm than to copy rows. For columns we can do it directly
! 1527: ** using e.g. swap(mat[i], mat[j])
! 1528: ** k0 = integer. The k0 first lines of mat are dense, the others are sparse.
! 1529: ** Also if k0=0, mat is modified in place [from mathnfspec], otherwise
! 1530: ** it is left alone [from buchall]
! 1531: ** Output: cf ASCII art in the function body
! 1532: **
! 1533: ** row permutations applied to perm
! 1534: ** column operations applied to C.
! 1535: **/
! 1536: GEN
! 1537: hnfspec(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
! 1538: {
! 1539: long av=avma,av2,*p,i,j,k,lk0,col,lig,*matj, **mat;
! 1540: long n,s,t,lim,nlze,lnz,nr;
! 1541: GEN p1,p2,matb,matbnew,vmax,matt,T,extramat;
! 1542: GEN B,H,dep,permpro;
! 1543: GEN *gptr[4];
! 1544: long co = lg(mat0);
! 1545: long li = lg(perm); /* = lg(mat0[1]) */
! 1546: int updateT = 1;
! 1547:
! 1548: if (!k0) mat = mat0; /* in place */
! 1549: else
! 1550: { /* keep original mat0 safe, modify a copy */
! 1551: mat = (long**)new_chunk(co); mat[0] = mat0[0];
! 1552: for (j=1; j<co; j++) mat[j] = col_dup(li,mat0[j]);
! 1553: }
! 1554:
! 1555: if (DEBUGLEVEL>5)
! 1556: {
! 1557: fprintferr("Entering hnfspec\n");
! 1558: p_mat(mat,perm,0);
! 1559: }
! 1560: matt = cgetg(co,t_MAT); /* dense part of mat (top) */
! 1561: for (j=1; j<co; j++)
! 1562: {
! 1563: p1=cgetg(k0+1,t_COL); matt[j]=(long)p1; matj = mat[j];
! 1564: for (i=1; i<=k0; i++) p1[i] = lstoi(matj[perm[i]]);
! 1565: }
! 1566: vmax = cgetg(co,t_VECSMALL);
! 1567: av2 = avma; lim = stack_lim(av2,1);
! 1568:
! 1569: i=lig=li-1; col=co-1; lk0=k0;
! 1570: if (k0 || (lg(*ptC) > 1 && lg((*ptC)[1]) > 1)) T = idmat(col);
! 1571: else
! 1572: { /* dummy ! */
! 1573: GEN z = cgetg(1,t_COL);
! 1574: T = cgetg(co, t_MAT); updateT = 0;
! 1575: for (j=1; j<co; j++) T[j] = (long)z;
! 1576: }
! 1577: /* Look for lines with a single non0 entry, equal to ±1 */
! 1578: while (i > lk0)
! 1579: switch( count(mat,perm[i],col,&n) )
! 1580: {
! 1581: case 0: /* move zero lines between k0+1 and lk0 */
! 1582: lk0++; swap(perm[i], perm[lk0]);
! 1583: i=lig; continue;
! 1584:
! 1585: case 1: /* move trivial generator between lig+1 and li */
! 1586: swap(perm[i], perm[lig]);
! 1587: swap(T[n], T[col]);
! 1588: gswap(mat[n], mat[col]); p = mat[col];
! 1589: if (p[perm[lig]] < 0) /* = -1 */
! 1590: { /* convert relation -g = 0 to g = 0 */
! 1591: for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
! 1592: if (updateT)
! 1593: {
! 1594: p1 = (GEN)T[col];
! 1595: for (i=1; ; i++)
! 1596: if (signe((GEN)p1[i])) { p1[i] = lnegi((GEN)p1[i]); break; }
! 1597: }
! 1598: }
! 1599: lig--; col--; i=lig; continue;
! 1600:
! 1601: default: i--;
! 1602: }
! 1603: if (DEBUGLEVEL>5)
! 1604: {
! 1605: fprintferr(" after phase1:\n");
! 1606: p_mat(mat,perm,0);
! 1607: }
! 1608:
! 1609: #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
! 1610:
! 1611: #if 0 /* TODO: check, and put back in */
! 1612: /* Get rid of all lines containing only 0 and ± 1, keeping track of column
! 1613: * operations in T. Leave the rows 1..lk0 alone [up to k0, coeff
! 1614: * explosion, between k0+1 and lk0, row is 0]
! 1615: */
! 1616: s = 0;
! 1617: while (lig > lk0 && s < (HIGHBIT>>1))
! 1618: {
! 1619: for (i=lig; i>lk0; i--)
! 1620: if (count(mat,perm[i],col,&n) >= 0) break;
! 1621: if (i == lk0) break;
! 1622:
! 1623: /* only 0, ±1 entries, at least 2 of them non-zero */
! 1624: swap(perm[i], perm[lig]);
! 1625: swap(T[n], T[col]); p1 = (GEN)T[col];
! 1626: gswap(mat[n], mat[col]); p = mat[col];
! 1627: if (p[perm[lig]] < 0)
! 1628: {
! 1629: for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
! 1630: T[col] = lneg(p1);
! 1631: }
! 1632: for (j=1; j<n; j++)
! 1633: {
! 1634: matj = mat[j];
! 1635: if (! (t = matj[perm[lig]]) ) continue;
! 1636: if (t == 1)
! 1637: { /* t = 1 */
! 1638: for (i=lk0+1; i<=lig; i++)
! 1639: absmax(s, matj[perm[i]] -= p[perm[i]]);
! 1640: T[j] = lsub((GEN)T[j], p1);
! 1641: }
! 1642: else
! 1643: { /* t = -1 */
! 1644: for (i=lk0+1; i<=lig; i++)
! 1645: absmax(s, matj[perm[i]] += p[perm[i]]);
! 1646: T[j] = ladd((GEN)T[j], p1);
! 1647: }
! 1648: }
! 1649: lig--; col--;
! 1650: if (low_stack(lim, stack_lim(av2,1)))
! 1651: {
! 1652: if(DEBUGMEM>1) err(warnmem,"hnfspec[1]");
! 1653: T = gerepilecopy(av2, T);
! 1654: }
! 1655: }
! 1656: #endif
! 1657: /* As above with lines containing a ±1 (no other assumption).
! 1658: * Stop when single precision becomes dangerous */
! 1659: for (j=1; j<=col; j++)
! 1660: {
! 1661: matj = mat[j];
! 1662: for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
! 1663: vmax[j] = s;
! 1664: }
! 1665: while (lig > lk0)
! 1666: {
! 1667: for (i=lig; i>lk0; i--)
! 1668: if ( (n = count2(mat,perm[i],col)) ) break;
! 1669: if (i == lk0) break;
! 1670:
! 1671: swap(perm[i], perm[lig]);
! 1672: swap(vmax[n], vmax[col]);
! 1673: gswap(mat[n], mat[col]); p = mat[col];
! 1674: swap(T[n], T[col]); p1 = (GEN)T[col];
! 1675: if (p[perm[lig]] < 0)
! 1676: {
! 1677: for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
! 1678: p1 = gneg(p1); T[col] = (long)p1;
! 1679: }
! 1680: for (j=1; j<col; j++)
! 1681: {
! 1682: matj = mat[j];
! 1683: if (! (t = matj[perm[lig]]) ) continue;
! 1684: if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
! 1685: goto END2;
! 1686:
! 1687: for (s=0, i=lk0+1; i<=lig; i++)
! 1688: absmax(s, matj[perm[i]] -= t*p[perm[i]]);
! 1689: vmax[j] = s;
! 1690: T[j] = (long)ZV_lincomb(gun,stoi(-t), (GEN)T[j],p1);
! 1691: }
! 1692: lig--; col--;
! 1693: if (low_stack(lim, stack_lim(av2,1)))
! 1694: {
! 1695: if(DEBUGMEM>1) err(warnmem,"hnfspec[2]");
! 1696: T = gerepilecopy(av2,T);
! 1697: }
! 1698: }
! 1699:
! 1700: END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
! 1701: /* go multiprecision first */
! 1702: matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
! 1703: for (j=1; j<co; j++)
! 1704: {
! 1705: p1 = cgetg(li-k0,t_COL); matb[j] = (long)p1;
! 1706: p1 -= k0; matj = mat[j];
! 1707: for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);
! 1708: }
! 1709: if (DEBUGLEVEL>5)
! 1710: {
! 1711: fprintferr(" after phase2:\n");
! 1712: p_mat(mat,perm,k0);
! 1713: }
! 1714: for (i=li-2; i>lig; i--)
! 1715: {
! 1716: long i1, i0 = i - k0, k = i + co-li;
! 1717: GEN Bk = (GEN)matb[k];
! 1718: GEN Tk = (GEN)T[k];
! 1719: for (j=k+1; j<co; j++)
! 1720: {
! 1721: p1=(GEN)matb[j]; p2=(GEN)p1[i0];
! 1722: if (! (s=signe(p2)) ) continue;
! 1723:
! 1724: p1[i0] = zero;
! 1725: if (is_pm1(p2))
! 1726: {
! 1727: if (s > 0)
! 1728: { /* p2 = 1 */
! 1729: for (i1=1; i1<i0; i1++)
! 1730: p1[i1] = lsubii((GEN)p1[i1], (GEN)Bk[i1]);
! 1731: T[j] = lsub((GEN)T[j], Tk);
! 1732: }
! 1733: else
! 1734: { /* p2 = -1 */
! 1735: for (i1=1; i1<i0; i1++)
! 1736: p1[i1] = laddii((GEN)p1[i1], (GEN)Bk[i1]);
! 1737: T[j] = ladd((GEN)T[j], Tk);
! 1738: }
! 1739: }
! 1740: else
! 1741: {
! 1742: for (i1=1; i1<i0; i1++)
! 1743: p1[i1] = lsubii((GEN)p1[i1], mulii(p2,(GEN) Bk[i1]));
! 1744: T[j] = (long)ZV_lincomb(gun,negi(p2), (GEN)T[j],Tk);
! 1745: }
! 1746: }
! 1747: if (low_stack(lim, stack_lim(av2,1)))
! 1748: {
! 1749: if(DEBUGMEM>1) err(warnmem,"hnfspec[3], i = %ld", i);
! 1750: for (j=1; j<co; j++) setlg(matb[j], i0+1); /* bottom can be forgotten */
! 1751: gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
! 1752: }
! 1753: }
! 1754: gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
! 1755: if (DEBUGLEVEL>5)
! 1756: {
! 1757: fprintferr(" matb cleaned up (using Id block)\n");
! 1758: if (DEBUGLEVEL>6) outerr(matb);
! 1759: }
! 1760:
! 1761: nlze = lk0 - k0; /* # of 0 rows */
! 1762: lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
! 1763: if (updateT) matt = gmul(matt,T); /* update top rows */
! 1764: extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
! 1765: for (j=1; j<=col; j++)
! 1766: {
! 1767: GEN z = (GEN)matt[j];
! 1768: GEN t = ((GEN)matb[j]) + nlze - k0;
! 1769: p2=cgetg(lnz,t_COL); extramat[j]=(long)p2;
! 1770: for (i=1; i<=k0; i++) p2[i] = z[i]; /* top k0 rows */
! 1771: for ( ; i<lnz; i++) p2[i] = t[i]; /* other non-0 rows */
! 1772: }
! 1773: permpro = imagecomplspec(extramat, &nr); /* lnz = lg(permpro) */
! 1774:
! 1775: if (nlze)
! 1776: { /* put the nlze 0 rows (trivial generators) at the top */
! 1777: p1 = new_chunk(lk0+1);
! 1778: for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
! 1779: for ( ; i<=lk0; i++) p1[i] = perm[i - nlze];
! 1780: for (i=1; i<=lk0; i++) perm[i] = p1[i];
! 1781: }
! 1782: /* sort other rows according to permpro (nr redundant generators first) */
! 1783: p1 = new_chunk(lnz); p2 = perm + nlze;
! 1784: for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
! 1785: for (i=1; i<lnz; i++) p2[i] = p1[i];
! 1786: /* perm indexes the rows of mat
! 1787: * |_0__|__redund__|__dense__|__too big__|_____done______|
! 1788: * 0 nlze lig li
! 1789: * \___nr___/ \___k0__/
! 1790: * \____________lnz ______________/
! 1791: *
! 1792: * col co
! 1793: * [dep | ]
! 1794: * i0 [--------| B ] (i0 = nlze + nr)
! 1795: * [matbnew | ] matbnew has maximal rank = lnz-1 - nr
! 1796: * mat = [--------|-----] lig
! 1797: * [ 0 | Id ]
! 1798: * [ | ] li */
! 1799:
! 1800: matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
! 1801: dep = cgetg(col+1,t_MAT); /* rows dependant from the ones in matbnew */
! 1802: for (j=1; j<=col; j++)
! 1803: {
! 1804: GEN z = (GEN)extramat[j];
! 1805: p1 = cgetg(nlze+nr+1,t_COL); dep[j]=(long)p1;
! 1806: p2 = cgetg(lnz-nr,t_COL); matbnew[j]=(long)p2;
! 1807: for (i=1; i<=nlze; i++) p1[i]=zero;
! 1808: p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
! 1809: p2 -= nr; for ( ; i<lnz; i++) p2[i] = z[permpro[i]];
! 1810: }
! 1811:
! 1812: /* redundant generators in terms of the genuine generators
! 1813: * (x_i) = - (g_i) B */
! 1814: B = cgetg(co-col,t_MAT);
! 1815: for (j=col+1; j<co; j++)
! 1816: {
! 1817: GEN y = (GEN)matt[j];
! 1818: GEN z = (GEN)matb[j];
! 1819: p1=cgetg(lig+1,t_COL); B[j-col]=(long)p1;
! 1820: for (i=1; i<=nlze; i++) p1[i] = z[i];
! 1821: p1 += nlze; z += nlze-k0;
! 1822: for (k=1; k<lnz; k++)
! 1823: {
! 1824: i = permpro[k];
! 1825: p1[k] = (i <= k0)? y[i]: z[i];
! 1826: }
! 1827: }
! 1828: if (updateT) *ptC = gmul(*ptC,T);
! 1829: *ptdep = dep;
! 1830: *ptB = B;
! 1831: H = hnffinal(matbnew,perm,ptdep,ptB,ptC);
! 1832: gptr[0]=ptC;
! 1833: gptr[1]=ptdep;
! 1834: gptr[2]=ptB;
! 1835: gptr[3]=&H; gerepilemany(av,gptr,4);
! 1836: if (DEBUGLEVEL)
! 1837: msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1);
! 1838: return H;
! 1839: }
! 1840:
! 1841: /* HNF reduce x, apply same transforms to C */
! 1842: GEN
! 1843: mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
! 1844: {
! 1845: long i,j,k,ly,lx = lg(x);
! 1846: GEN p1,p2,z,perm;
! 1847: if (lx == 1) return gcopy(x);
! 1848: ly = lg(x[1]);
! 1849: z = cgetg(lx,t_MAT);
! 1850: perm = cgetg(ly,t_VECSMALL); *ptperm = perm;
! 1851: for (i=1; i<ly; i++) perm[i] = i;
! 1852: for (i=1; i<lx; i++)
! 1853: {
! 1854: p1 = cgetg(ly,t_COL); z[i] = (long)p1;
! 1855: p2 = (GEN)x[i];
! 1856: for (j=1; j<ly; j++)
! 1857: {
! 1858: if (is_bigint(p2[j])) goto TOOLARGE;
! 1859: p1[j] = itos((GEN)p2[j]);
! 1860: }
! 1861: }
! 1862: /* [ dep | ]
! 1863: * [-----| B ]
! 1864: * [ H | ]
! 1865: * [-----|-----]
! 1866: * [ 0 | Id ] */
! 1867: return hnfspec((long**)z,perm, ptdep, ptB, ptC, 0);
! 1868:
! 1869: TOOLARGE:
! 1870: if (lg(*ptC) > 1 && lg((*ptC)[1]) > 1)
! 1871: err(impl,"mathnfspec with large entries");
! 1872: x = hnf(x); lx = lg(x); j = ly; k = 0;
! 1873: for (i=1; i<ly; i++)
! 1874: {
! 1875: if (gcmp1(gcoeff(x,i,i + lx-ly)))
! 1876: perm[--j] = i;
! 1877: else
! 1878: perm[++k] = i;
! 1879: }
! 1880: setlg(perm,k+1);
! 1881: x = rowextract_p(x, perm); /* upper part */
! 1882: setlg(perm,ly);
! 1883: *ptB = vecextract_i(x, j+lx-ly, lx-1);
! 1884: setlg(x, j);
! 1885: *ptdep = rowextract_i(x, 1, lx-ly);
! 1886: return rowextract_i(x, lx-ly+1, k); /* H */
! 1887: }
! 1888:
! 1889: /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
! 1890: GEN
! 1891: hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
! 1892: GEN extramat,GEN extraC)
! 1893: {
! 1894: GEN p1,p2,p3,matb,extratop,Cnew,permpro;
! 1895: GEN B=*ptB, C=*ptC, dep=*ptdep, *gptr[4];
! 1896: long av = avma, i,j,lextra,colnew;
! 1897: long li = lg(perm);
! 1898: long co = lg(C);
! 1899: long lB = lg(B);
! 1900: long lig = li - lB;
! 1901: long col = co - lB;
! 1902: long lH = lg(H)-1;
! 1903: long nlze = lH? lg(dep[1])-1: lg(B[1])-1;
! 1904:
! 1905: if (DEBUGLEVEL>5)
! 1906: {
! 1907: fprintferr("Entering hnfadd:\n");
! 1908: if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);
! 1909: }
! 1910: /* col co
! 1911: * [ 0 |dep | ]
! 1912: * nlze [--------| B ]
! 1913: * [ 0 | H | ]
! 1914: * [--------|-----] lig
! 1915: * [ 0 | Id ]
! 1916: * [ | ] li */
! 1917: lextra = lg(extramat)-1;
! 1918: extratop = cgetg(lextra+1,t_MAT); /* [1..lig] part (top) */
! 1919: p2 = cgetg(lextra+1,t_MAT); /* bottom */
! 1920: for (j=1; j<=lextra; j++)
! 1921: {
! 1922: GEN z = (GEN)extramat[j];
! 1923: p1=cgetg(lig+1,t_COL); extratop[j] = (long)p1;
! 1924: for (i=1; i<=lig; i++) p1[i] = z[i];
! 1925: p1=cgetg(lB,t_COL); p2[j] = (long)p1;
! 1926: p1 -= lig;
! 1927: for ( ; i<li; i++) p1[i] = z[i];
! 1928: }
! 1929: if (li-1 != lig)
! 1930: { /* zero out bottom part, using the Id block */
! 1931: GEN A = cgetg(lB,t_MAT);
! 1932: for (j=1; j<lB; j++) A[j] = C[j+col];
! 1933: extraC = gsub(extraC, gmul(A,p2));
! 1934: extratop = gsub(extratop,gmul(B,p2));
! 1935: }
! 1936:
! 1937: colnew = lH + lextra;
! 1938: extramat = cgetg(colnew+1,t_MAT);
! 1939: Cnew = cgetg(lB+colnew,t_MAT);
! 1940: for (j=1; j<=lextra; j++)
! 1941: {
! 1942: extramat[j] = extratop[j];
! 1943: Cnew[j] = extraC[j];
! 1944: }
! 1945: for ( ; j<=colnew; j++)
! 1946: {
! 1947: p1 = cgetg(lig+1,t_COL); extramat[j] = (long)p1;
! 1948: p2 = (GEN)dep[j-lextra]; for (i=1; i<=nlze; i++) p1[i] = p2[i];
! 1949: p2 = (GEN) H[j-lextra]; for ( ; i<=lig ; i++) p1[i] = p2[i-nlze];
! 1950: }
! 1951: for (j=lextra+1; j<lB+colnew; j++)
! 1952: Cnew[j] = C[j-lextra+col-lH];
! 1953: if (DEBUGLEVEL>5)
! 1954: {
! 1955: fprintferr(" 1st phase done\n");
! 1956: if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);
! 1957: }
! 1958: permpro = imagecomplspec(extramat, &nlze);
! 1959: p1 = new_chunk(lig+1);
! 1960: for (i=1; i<=lig; i++) p1[i] = perm[permpro[i]];
! 1961: for (i=1; i<=lig; i++) perm[i] = p1[i];
! 1962:
! 1963: matb = cgetg(colnew+1,t_MAT);
! 1964: dep = cgetg(colnew+1,t_MAT);
! 1965: for (j=1; j<=colnew; j++)
! 1966: {
! 1967: GEN z = (GEN)extramat[j];
! 1968: p1=cgetg(nlze+1,t_COL); dep[j]=(long)p1;
! 1969: p2=cgetg(lig+1-nlze,t_COL); matb[j]=(long)p2;
! 1970: p2 -= nlze;
! 1971: for (i=1; i<=nlze; i++) p1[i] = z[permpro[i]];
! 1972: for ( ; i<= lig; i++) p2[i] = z[permpro[i]];
! 1973: }
! 1974: p3 = cgetg(lB,t_MAT);
! 1975: for (j=1; j<lB; j++)
! 1976: {
! 1977: p2 = (GEN)B[j];
! 1978: p1 = cgetg(lig+1,t_COL); p3[j] = (long)p1;
! 1979: for (i=1; i<=lig; i++) p1[i] = p2[permpro[i]];
! 1980: }
! 1981: B = p3;
! 1982: if (DEBUGLEVEL>5) fprintferr(" 2nd phase done\n");
! 1983: *ptdep = dep;
! 1984: *ptB = B;
! 1985: H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
! 1986: p1 = cgetg(co+lextra,t_MAT);
! 1987: for (j=1; j <= col-lH; j++) p1[j] = C[j];
! 1988: C = Cnew - (col-lH);
! 1989: for ( ; j < co+lextra; j++) p1[j] = C[j];
! 1990:
! 1991: gptr[0]=ptC; *ptC=p1;
! 1992: gptr[1]=ptdep;
! 1993: gptr[2]=ptB;
! 1994: gptr[3]=&H; gerepilemany(av,gptr,4);
! 1995: if (DEBUGLEVEL)
! 1996: {
! 1997: if (DEBUGLEVEL>7) fprintferr("mit = %Z\nC = %Z\n",H,*ptC);
! 1998: msgtimer("hnfadd (%ld)",lextra);
! 1999: }
! 2000: return H;
! 2001: }
! 2002:
! 2003: static void
! 2004: ZV_neg(GEN x)
! 2005: {
! 2006: long i, lx = lg(x);
! 2007: for (i=1; i<lx ; i++) x[i]=lnegi((GEN)x[i]);
! 2008: }
! 2009:
! 2010: /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
! 2011: * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
! 2012: static void
! 2013: ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
! 2014: {
! 2015: GEN p1,u,v,d;
! 2016:
! 2017: if (!signe(ak)) { swap(A[j],A[k]); if (U) swap(U[j],U[k]); return; }
! 2018: d = bezout(aj,ak,&u,&v);
! 2019: /* frequent special case (u,v) = (1,0) or (0,1) */
! 2020: if (!signe(u))
! 2021: { /* ak | aj */
! 2022: p1 = negi(divii(aj,ak));
! 2023: A[j] = (long)ZV_lincomb(gun, p1, (GEN)A[j], (GEN)A[k]);
! 2024: if (U)
! 2025: U[j] = (long)ZV_lincomb(gun, p1, (GEN)U[j], (GEN)U[k]);
! 2026: return;
! 2027: }
! 2028: if (!signe(v))
! 2029: { /* aj | ak */
! 2030: p1 = negi(divii(ak,aj));
! 2031: A[k] = (long)ZV_lincomb(gun, p1, (GEN)A[k], (GEN)A[j]);
! 2032: swap(A[j], A[k]);
! 2033: if (U) {
! 2034: U[k] = (long)ZV_lincomb(gun, p1, (GEN)U[k], (GEN)U[j]);
! 2035: swap(U[j], U[k]);
! 2036: }
! 2037: return;
! 2038: }
! 2039:
! 2040: if (!is_pm1(d)) { aj = divii(aj,d); ak = divii(ak,d); }
! 2041: p1 = (GEN)A[k]; aj = negi(aj);
! 2042: A[k] = (long)ZV_lincomb(u,v, (GEN)A[j],p1);
! 2043: A[j] = (long)ZV_lincomb(aj,ak, p1,(GEN)A[j]);
! 2044: if (U)
! 2045: {
! 2046: p1 = (GEN)U[k];
! 2047: U[k] = (long)ZV_lincomb(u,v, (GEN)U[j],p1);
! 2048: U[j] = (long)ZV_lincomb(aj,ak, p1,(GEN)U[j]);
! 2049: }
! 2050: }
! 2051:
! 2052: /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
! 2053: static void
! 2054: ZM_reduce(GEN A, GEN U, long i, long j0)
! 2055: {
! 2056: long j, lA = lg(A);
! 2057: GEN d = gcoeff(A,i,j0);
! 2058: if (signe(d) < 0)
! 2059: {
! 2060: ZV_neg((GEN)A[j0]);
! 2061: if (U) ZV_neg((GEN)U[j0]);
! 2062: d = gcoeff(A,i,j0);
! 2063: }
! 2064: for (j=j0+1; j<lA; j++)
! 2065: {
! 2066: GEN q = truedvmdii(gcoeff(A,i,j), d, NULL);
! 2067: if (!signe(q)) continue;
! 2068:
! 2069: q = negi(q);
! 2070: A[j] = (long)ZV_lincomb(gun,q, (GEN)A[j], (GEN)A[j0]);
! 2071: if (U)
! 2072: U[j] = (long)ZV_lincomb(gun,q,(GEN)U[j],(GEN)U[j0]);
! 2073: }
! 2074: }
! 2075:
! 2076: /* remove: throw away lin.dep.columns */
! 2077: GEN
! 2078: hnf0(GEN A, int remove)
! 2079: {
! 2080: long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
! 2081: GEN denx,a,p1;
! 2082:
! 2083: if (typ(A) == t_VEC) return hnf_special(A,remove);
! 2084: A = init_hnf(A,&denx,&co,&li,&av);
! 2085: if (!A) return cgetg(1,t_MAT);
! 2086:
! 2087: lim = stack_lim(av,1);
! 2088: def=co-1; ldef=(li>co)? li-co: 0;
! 2089: for (i=li-1; i>ldef; i--)
! 2090: {
! 2091: for (j=def-1; j; j--)
! 2092: {
! 2093: a = gcoeff(A,i,j);
! 2094: if (!signe(a)) continue;
! 2095:
! 2096: /* zero a = Aij using b = Aik */
! 2097: k = (j==1)? def: j-1;
! 2098: ZV_elem(a,gcoeff(A,i,k), A,NULL, j,k);
! 2099:
! 2100: if (low_stack(lim, stack_lim(av,1)))
! 2101: {
! 2102: if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i);
! 2103: tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));
! 2104: }
! 2105: }
! 2106: p1 = gcoeff(A,i,def); s = signe(p1);
! 2107: if (s)
! 2108: {
! 2109: if (s < 0) ZV_neg((GEN)A[def]);
! 2110: ZM_reduce(A, NULL, i,def);
! 2111: def--;
! 2112: }
! 2113: else
! 2114: if (ldef) ldef--;
! 2115: if (low_stack(lim, stack_lim(av,1)))
! 2116: {
! 2117: if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i);
! 2118: tetpil=avma; A=gerepile(av,tetpil,ZM_copy(A));
! 2119: }
! 2120: }
! 2121: if (remove)
! 2122: { /* remove null columns */
! 2123: for (i=1,j=1; j<co; j++)
! 2124: if (!gcmp0((GEN)A[j])) A[i++] = A[j];
! 2125: setlg(A,i);
! 2126: }
! 2127: tetpil=avma;
! 2128: A = denx? gdiv(A,denx): ZM_copy(A);
! 2129: return gerepile(av0,tetpil,A);
! 2130: }
! 2131:
! 2132: GEN
! 2133: hnf(GEN x) { return hnf0(x,1); }
! 2134:
! 2135: /* dm = multiple of diag element (usually detint(x)) */
! 2136: /* flag: don't/do append dm*matid. */
! 2137: static GEN
! 2138: allhnfmod(GEN x,GEN dm,int flag)
! 2139: {
! 2140: ulong av,tetpil,lim;
! 2141: long li,co,i,j,k,def,ldef,ldm;
! 2142: GEN a,b,w,p1,p2,d,u,v;
! 2143:
! 2144: if (typ(dm)!=t_INT) err(typeer,"allhnfmod");
! 2145: if (!signe(dm)) return hnf(x);
! 2146: if (typ(x)!=t_MAT) err(typeer,"allhnfmod");
! 2147:
! 2148: co=lg(x); if (co==1) return cgetg(1,t_MAT);
! 2149: li=lg(x[1]);
! 2150: av = avma; lim = stack_lim(av,1);
! 2151: x = dummycopy(x);
! 2152:
! 2153: if (flag)
! 2154: {
! 2155: if (li > co) err(talker,"nb lines > nb columns in hnfmod");
! 2156: }
! 2157: else
! 2158: { /* concatenate dm * Id to x */
! 2159: x = concatsp(x, idmat_intern(li-1,dm,gzero));
! 2160: co += li-1;
! 2161: }
! 2162: /* Avoid wasteful divisions. we only want to prevent coeff explosion, so
! 2163: * only reduce mod dm when lg(coeff) > ldm */
! 2164: ldm = lgefint(dm);
! 2165: def = co-1; ldef = 0;
! 2166: for (i=li-1; i>ldef; i--,def--)
! 2167: for (j=def-1; j; j--)
! 2168: {
! 2169: coeff(x,i,j) = lresii(gcoeff(x,i,j), dm);
! 2170: a = gcoeff(x,i,j);
! 2171: if (!signe(a)) continue;
! 2172:
! 2173: k = (j==1)? def: j-1;
! 2174: /* do not reduce the appended dm [hnfmodid] */
! 2175: if (flag || j != 1) coeff(x,i,k) = lresii(gcoeff(x,i,k), dm);
! 2176: ZV_elem(a,gcoeff(x,i,k), x,NULL, j,k);
! 2177: p1 = (GEN)x[j];
! 2178: p2 = (GEN)x[k];
! 2179: for (k=1; k<i; k++)
! 2180: {
! 2181: if (lgefint(p1[k]) > ldm) p1[k] = lresii((GEN)p1[k], dm);
! 2182: if (lgefint(p2[k]) > ldm) p2[k] = lresii((GEN)p2[k], dm);
! 2183: }
! 2184: if (low_stack(lim, stack_lim(av,1)))
! 2185: {
! 2186: if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i);
! 2187: tetpil=avma; x=gerepile(av,tetpil,ZM_copy(x));
! 2188: }
! 2189: }
! 2190: w = cgetg(li,t_MAT); b = dm;
! 2191: for (i=li-1; i>=1; i--)
! 2192: {
! 2193: d = bezout(gcoeff(x,i,i+def),b,&u,&v);
! 2194: w[i] = lmod(gmul(u,(GEN)x[i+def]), b);
! 2195: if (!signe(gcoeff(w,i,i))) coeff(w,i,i) = (long)d;
! 2196: if (flag && i > 1) b = diviiexact(b,d);
! 2197: }
! 2198: if (flag)
! 2199: { /* compute optimal value for dm */
! 2200: dm = gcoeff(w,1,1);
! 2201: for (i=2; i<li; i++) dm = mpppcm(dm, gcoeff(w,i,i));
! 2202: ldm = lgefint(dm);
! 2203: }
! 2204: for (i=li-2; i>=1; i--)
! 2205: {
! 2206: GEN diag = gcoeff(w,i,i);
! 2207: for (j=i+1; j<li; j++)
! 2208: {
! 2209: b = negi(truedvmdii(gcoeff(w,i,j), diag, NULL));
! 2210: p1 = ZV_lincomb(gun,b, (GEN)w[j], (GEN)w[i]);
! 2211: w[j] = (long)p1;
! 2212: for (k=1; k<i; k++)
! 2213: if (lgefint(p1[k]) > ldm) p1[k] = lresii((GEN)p1[k], dm);
! 2214: if (low_stack(lim, stack_lim(av,1)))
! 2215: {
! 2216: if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i);
! 2217: tetpil=avma; w=gerepile(av,tetpil,ZM_copy(w));
! 2218: diag = gcoeff(w,i,i);
! 2219: }
! 2220: }
! 2221: }
! 2222: tetpil=avma; return gerepile(av,tetpil,ZM_copy(w));
! 2223: }
! 2224:
! 2225: GEN
! 2226: hnfmod(GEN x, GEN detmat) { return allhnfmod(x,detmat,1); }
! 2227:
! 2228: GEN
! 2229: hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }
! 2230:
! 2231: /***********************************************************************/
! 2232: /* */
! 2233: /* HNFLLL (Havas, Majewski, Mathews) */
! 2234: /* */
! 2235: /***********************************************************************/
! 2236:
! 2237: /* negate in place, except universal constants */
! 2238: static GEN
! 2239: mynegi(GEN x)
! 2240: {
! 2241: static long mun[]={evaltyp(t_INT)|m_evallg(3),evalsigne(-1)|evallgefint(3),1};
! 2242: long s = signe(x);
! 2243:
! 2244: if (!s) return gzero;
! 2245: if (is_pm1(x))
! 2246: return (s>0)? mun: gun;
! 2247: setsigne(x,-s); return x;
! 2248: }
! 2249:
! 2250: /* We assume y > 0 */
! 2251: static GEN
! 2252: divnearest(GEN x, GEN y)
! 2253: {
! 2254: GEN r, q = dvmdii(x,y,&r);
! 2255: long av = avma, s=signe(r), t;
! 2256:
! 2257: if (!s) { cgiv(r); return q; }
! 2258: if (s<0) r = mynegi(r);
! 2259:
! 2260: y = shifti(y,-1); t = cmpii(r,y);
! 2261: avma=av; cgiv(r);
! 2262: if (t < 0 || (t == 0 && s > 0)) return q;
! 2263: return addsi(s,q);
! 2264: }
! 2265:
! 2266: static void
! 2267: Minus(long j, GEN **lambda)
! 2268: {
! 2269: long k, n = lg(lambda);
! 2270:
! 2271: for (k=1 ; k<j; k++) lambda[j][k] = mynegi(lambda[j][k]);
! 2272: for (k=j+1; k<n; k++) lambda[k][j] = mynegi(lambda[k][j]);
! 2273: }
! 2274:
! 2275: static void
! 2276: neg_col(GEN M)
! 2277: {
! 2278: long i;
! 2279: for (i = lg(M)-1; i; i--)
! 2280: M[i] = (long)mynegi((GEN)M[i]);
! 2281: }
! 2282:
! 2283: /* M_k = M_k + q M_i (col operations) */
! 2284: static void
! 2285: elt_col(GEN Mk, GEN Mi, GEN q)
! 2286: {
! 2287: long j;
! 2288: if (is_pm1(q))
! 2289: {
! 2290: if (signe(q) > 0)
! 2291: {
! 2292: for (j = lg(Mk)-1; j; j--)
! 2293: if (signe(Mi[j])) Mk[j] = laddii((GEN)Mk[j], (GEN)Mi[j]);
! 2294: }
! 2295: else
! 2296: {
! 2297: for (j = lg(Mk)-1; j; j--)
! 2298: if (signe(Mi[j])) Mk[j] = lsubii((GEN)Mk[j], (GEN)Mi[j]);
! 2299: }
! 2300: }
! 2301: else
! 2302: for (j = lg(Mk)-1; j; j--)
! 2303: if (signe(Mi[j])) Mk[j] = laddii((GEN)Mk[j], mulii(q, (GEN)Mi[j]));
! 2304: }
! 2305:
! 2306: /* index of first non-zero entry */
! 2307: static long
! 2308: findi(GEN M)
! 2309: {
! 2310: long i, n = lg(M);
! 2311: for (i=1; i<n; i++)
! 2312: if (signe(M[i])) return i;
! 2313: return 0;
! 2314: }
! 2315:
! 2316: static void
! 2317: reduce2(GEN A, GEN B, long k, long j, long *row, GEN **lambda, GEN *D)
! 2318: {
! 2319: GEN q;
! 2320: long i, row0, row1;
! 2321:
! 2322: row0 = findi((GEN)A[j]);
! 2323: if (row0 && signe(gcoeff(A,row0,j)) < 0)
! 2324: {
! 2325: neg_col((GEN)A[j]);
! 2326: if (B) neg_col((GEN)B[j]);
! 2327: Minus(j,lambda);
! 2328: }
! 2329:
! 2330: row1 = findi((GEN)A[k]);
! 2331: if (row1 && signe(gcoeff(A,row1,k)) < 0)
! 2332: {
! 2333: neg_col((GEN)A[k]);
! 2334: if (B) neg_col((GEN)B[k]);
! 2335: Minus(k,lambda);
! 2336: }
! 2337:
! 2338: row[0]=row0; row[1]=row1;
! 2339: if (row0)
! 2340: q = truedvmdii(gcoeff(A,row0,k), gcoeff(A,row0,j), NULL);
! 2341: else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
! 2342: q = divnearest(lambda[k][j], D[j]);
! 2343: else
! 2344: return;
! 2345:
! 2346: if (signe(q))
! 2347: {
! 2348: GEN *Lk = lambda[k], *Lj = lambda[j];
! 2349: q = mynegi(q);
! 2350: if (row0) elt_col((GEN)A[k],(GEN)A[j],q);
! 2351: if (B) elt_col((GEN)B[k],(GEN)B[j],q);
! 2352: Lk[j] = addii(Lk[j], mulii(q,D[j]));
! 2353: if (is_pm1(q))
! 2354: {
! 2355: if (signe(q) > 0)
! 2356: {
! 2357: for (i=1; i<j; i++)
! 2358: if (signe(Lj[i])) Lk[i] = addii(Lk[i], Lj[i]);
! 2359: }
! 2360: else
! 2361: {
! 2362: for (i=1; i<j; i++)
! 2363: if (signe(Lj[i])) Lk[i] = subii(Lk[i], Lj[i]);
! 2364: }
! 2365: }
! 2366: else
! 2367: {
! 2368: for (i=1; i<j; i++)
! 2369: if (signe(Lj[i])) Lk[i] = addii(Lk[i], mulii(q,Lj[i]));
! 2370: }
! 2371: }
! 2372: }
! 2373:
! 2374: #define SWAP(x,y) {GEN _t=y; y=x; x=_t;}
! 2375:
! 2376: static void
! 2377: hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)
! 2378: {
! 2379: GEN t,p1,p2;
! 2380: long i,j,n = lg(A);
! 2381:
! 2382: swap(A[k],A[k-1]);
! 2383: if (B) swap(B[k],B[k-1]);
! 2384: for (j=k-2; j; j--) SWAP(lambda[k-1][j], lambda[k][j]);
! 2385: for (i=k+1; i<n; i++)
! 2386: {
! 2387: p1 = mulii(lambda[i][k-1], D[k]);
! 2388: p2 = mulii(lambda[i][k], lambda[k][k-1]);
! 2389: t = subii(p1,p2);
! 2390:
! 2391: p1 = mulii(lambda[i][k], D[k-2]);
! 2392: p2 = mulii(lambda[i][k-1], lambda[k][k-1]);
! 2393: lambda[i][k-1] = divii(addii(p1,p2), D[k-1]);
! 2394:
! 2395: lambda[i][k] = divii(t, D[k-1]);
! 2396: }
! 2397: p1 = mulii(D[k-2],D[k]);
! 2398: p2 = sqri(lambda[k][k-1]);
! 2399: D[k-1] = divii(addii(p1,p2), D[k-1]);
! 2400: }
! 2401:
! 2402: /* A[k] = 0, A[nz] != 0, k > nz, A[j] = 0 for all j < nz.
! 2403: * "Deep insert" A[k] at nz */
! 2404: static void
! 2405: trivswap(GEN *A, long nz, long k, GEN **lambda, GEN *D)
! 2406: {
! 2407: GEN p1;
! 2408: long i,j,n = lg(A);
! 2409:
! 2410: p1 = A[nz]; /* cycle A */
! 2411: for (j = nz; j < k; j++) SWAP(A[j+1], p1);
! 2412: A[nz] = p1; /* = [0...0] */
! 2413:
! 2414: p1 = D[nz]; /* cycle D */
! 2415: for (j = nz; j < k; j++) SWAP(D[j+1], p1);
! 2416: D[nz] = gun;
! 2417:
! 2418: for (j=k-1; j>=nz; j--) /* cycle lambda */
! 2419: for (i=k-1; i>=nz; i--) lambda[i+1][j+1] = lambda[i][j];
! 2420: for (j=n-1; j>k; j--)
! 2421: for (i=k-1; i>=nz; i--)
! 2422: {
! 2423: lambda[i+1][j] = lambda[i][j];
! 2424: lambda[j][i+1] = lambda[j][i];
! 2425: }
! 2426: for (i=1; i<n; i++) lambda[nz][i] = lambda[i][nz] = gzero;
! 2427: }
! 2428:
! 2429: static GEN
! 2430: fix_rows(GEN A)
! 2431: {
! 2432: long i,j, h,n = lg(A);
! 2433: GEN cB,cA,B = cgetg(n,t_MAT);
! 2434: if (n == 1) return B;
! 2435: h = lg(A[1]);
! 2436: for (j=1; j<n; j++)
! 2437: {
! 2438: cB = cgetg(h, t_COL);
! 2439: cA = (GEN)A[j]; B[j] = (long)cB;
! 2440: for (i=h>>1; i; i--) { cB[h-i] = cA[i]; cB[i] = cA[h-i]; }
! 2441: }
! 2442: return B;
! 2443: }
! 2444:
! 2445: GEN
! 2446: hnflll_i(GEN A, GEN *ptB)
! 2447: {
! 2448: ulong av = avma, lim = stack_lim(av,3);
! 2449: long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
! 2450: long row[2], do_swap,i,n,k;
! 2451: long nzcol = 1; /* index of 1st (maybe) non-0 col [only updated when !B] */
! 2452: GEN z,B, **lambda, *D;
! 2453: GEN *gptr[4];
! 2454:
! 2455: if (typ(A) != t_MAT) err(typeer,"hnflll");
! 2456: n = lg(A);
! 2457: A = ZM_copy(fix_rows(A));
! 2458: B = ptB? idmat(n-1): NULL;
! 2459: D = (GEN*) cgetg(n+1, t_VEC); D++; /* hack: need a "sentinel" D[0] */
! 2460: if (n == 2) /* handle trivial case: return negative diag coeff otherwise */
! 2461: {
! 2462: i = findi((GEN)A[1]);
! 2463: if (i && signe(gcoeff(A,i,1)) < 0)
! 2464: {
! 2465: neg_col((GEN)A[1]);
! 2466: if (B) neg_col((GEN)B[1]);
! 2467: }
! 2468: }
! 2469: lambda = (GEN**) cgetg(n,t_MAT);
! 2470: for (i=1; i<n; i++) { D[i] = gun; lambda[i] = (GEN*)zerocol(n-1); }
! 2471: D[0] = gun; k = 2;
! 2472: while (k < n)
! 2473: {
! 2474: reduce2(A,B,k,k-1,row,lambda,D);
! 2475: if (!B)
! 2476: { /* not interested in B. Eliminate 0 cols */
! 2477: if (!row[0] || !row[1])
! 2478: {
! 2479: while (!findi((GEN)A[nzcol]) && nzcol < n) nzcol++;
! 2480: /* A[nzcol] != 0, A[i] = 0 for i < nzcol */
! 2481: if (!row[0] && k-1 > nzcol)
! 2482: {trivswap((GEN*)A,nzcol,k-1, lambda,D); nzcol++;}
! 2483: if (!row[1] && k > nzcol)
! 2484: {trivswap((GEN*)A,nzcol,k , lambda,D); nzcol++;}
! 2485: if (k <= nzcol) k = nzcol+1;
! 2486: continue;
! 2487: }
! 2488: do_swap = (row[0] && row[0] <= row[1]);
! 2489: }
! 2490: else
! 2491: {
! 2492: if (row[0])
! 2493: do_swap = (!row[1] || row[0] <= row[1]);
! 2494: else if (row[1])
! 2495: do_swap = 0;
! 2496: else
! 2497: { /* row[0] == row[1] == 0 */
! 2498: ulong av1 = avma;
! 2499: z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
! 2500: do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
! 2501: avma = av1;
! 2502: }
! 2503: }
! 2504: if (do_swap)
! 2505: {
! 2506: hnfswap(A,B,k,lambda,D);
! 2507: if (k > nzcol+1) k--;
! 2508: }
! 2509: else
! 2510: {
! 2511: for (i=k-2; i>=nzcol; i--) reduce2(A,B,k,i,row,lambda,D);
! 2512: k++;
! 2513: }
! 2514: if (low_stack(lim, stack_lim(av,3)))
! 2515: {
! 2516: GEN a = (GEN)lambda, b = (GEN)(D-1); /* gcc -Wall */
! 2517: gptr[0]=&A; gptr[1]=&a; gptr[2]=&b; gptr[3]=&B;
! 2518: if (DEBUGMEM) err(warnmem,"hnflll, k = %ld / %ld",k,n);
! 2519: gerepilemany(av,gptr,B? 4: 3); lambda = (GEN**)a; D = (GEN*)(b+1);
! 2520: }
! 2521: }
! 2522: for (i=nzcol; i<n; i++)
! 2523: if (findi((GEN)A[i])) break;
! 2524: i--; A += i; A[0] = evaltyp(t_MAT)|evallg(n-i);
! 2525: A = fix_rows(A);
! 2526: gptr[0] = &A; gptr[1] = &B;
! 2527: gerepilemany(av, gptr, B? 2: 1);
! 2528: if (B) *ptB = B;
! 2529: return A;
! 2530: }
! 2531:
! 2532: GEN
! 2533: hnflll(GEN A)
! 2534: {
! 2535: GEN B, z = cgetg(3, t_VEC);
! 2536: z[1] = (long)hnflll_i(A, &B);
! 2537: z[2] = (long)B; return z;
! 2538: }
! 2539:
! 2540: /* Variation on HNFLLL: Extended GCD */
! 2541:
! 2542: static void
! 2543: reduce1(GEN A, GEN B, long k, long j, GEN **lambda, GEN *D)
! 2544: {
! 2545: GEN q;
! 2546: long i;
! 2547:
! 2548: if (signe(A[j]))
! 2549: q = divnearest((GEN)A[k],(GEN)A[j]);
! 2550: else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
! 2551: q = divnearest(lambda[k][j], D[j]);
! 2552: else
! 2553: return;
! 2554:
! 2555: if (! gcmp0(q))
! 2556: {
! 2557: q = mynegi(q);
! 2558: A[k] = laddii((GEN)A[k], mulii(q,(GEN)A[j]));
! 2559: elt_col((GEN)B[k],(GEN)B[j],q);
! 2560: lambda[k][j] = addii(lambda[k][j], mulii(q,D[j]));
! 2561: for (i=1; i<j; i++)
! 2562: if (signe(lambda[j][i]))
! 2563: lambda[k][i] = addii(lambda[k][i], mulii(q,lambda[j][i]));
! 2564: }
! 2565: }
! 2566:
! 2567: GEN
! 2568: extendedgcd(GEN A)
! 2569: {
! 2570: long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
! 2571: long av = avma, tetpil, do_swap,i,j,n,k;
! 2572: GEN p1,p2,B, **lambda, *D;
! 2573:
! 2574: n = lg(A); B = idmat(n-1); A = ZM_copy(A);
! 2575: D = (GEN*) cgeti(n); lambda = (GEN**) cgetg(n,t_MAT);
! 2576: for (i=0; i<n; i++) D[i] = gun;
! 2577: for (i=1; i<n; i++)
! 2578: {
! 2579: lambda[i] = (GEN*) cgetg(n,t_COL);
! 2580: for(j=1;j<n;j++) lambda[i][j] = gzero;
! 2581: }
! 2582: k = 2;
! 2583: while (k < n)
! 2584: {
! 2585: reduce1(A,B,k,k-1,lambda,D);
! 2586: if (signe(A[k-1])) do_swap = 1;
! 2587: else if (!signe(A[k]))
! 2588: {
! 2589: long av1=avma;
! 2590: p1 = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
! 2591: p1 = mulsi(n1,p1);
! 2592: p2 = mulsi(m1,sqri(D[k-1]));
! 2593: do_swap = (cmpii(p1,p2) < 0);
! 2594: avma=av1;
! 2595: }
! 2596: else do_swap = 0;
! 2597:
! 2598: if (do_swap)
! 2599: {
! 2600: hnfswap(A,B,k,lambda,D);
! 2601: if (k>2) k--;
! 2602: }
! 2603: else
! 2604: {
! 2605: for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
! 2606: k++;
! 2607: }
! 2608: }
! 2609: if (signe(A[n-1])<0)
! 2610: {
! 2611: A[n-1] = (long)mynegi((GEN)A[n-1]);
! 2612: neg_col((GEN)B[n-1]);
! 2613: }
! 2614: tetpil = avma;
! 2615: p1 = cgetg(3,t_VEC);
! 2616: p1[1]=lcopy((GEN)A[n-1]);
! 2617: p1[2]=lcopy(B);
! 2618: return gerepile(av,tetpil,p1);
! 2619: }
! 2620:
! 2621: /* HNF with permutation. TODO: obsolete, replace with hnflll. */
! 2622: GEN
! 2623: hnfperm(GEN A)
! 2624: {
! 2625: GEN U,c,l,perm,d,u,p,q,y,b;
! 2626: long r,t,i,j,j1,k,m,n,av=avma,av1,tetpil,lim;
! 2627:
! 2628: if (typ(A) != t_MAT) err(typeer,"hnfperm");
! 2629: n = lg(A)-1;
! 2630: if (!n)
! 2631: {
! 2632: y = cgetg(4,t_VEC);
! 2633: y[1] = lgetg(1,t_MAT);
! 2634: y[2] = lgetg(1,t_MAT);
! 2635: y[3] = lgetg(1,t_VEC); return y;
! 2636: }
! 2637: m = lg(A[1])-1;
! 2638: c = cgetg(m+1,t_VECSMALL); for (i=1; i<=m; i++) c[i]=0;
! 2639: l = cgetg(n+1,t_VECSMALL); for (j=1; j<=n; j++) l[j]=0;
! 2640: perm = cgetg(m+1, t_VECSMALL);
! 2641: av1 = avma; lim = stack_lim(av1,1);
! 2642: U = idmat(n); A = dummycopy(A); /* U base change matrix : A0*U=A all along */
! 2643: for (r=0,k=1; k<=n; k++)
! 2644: {
! 2645: for (j=1; j<k; j++)
! 2646: {
! 2647: if (!l[j]) continue;
! 2648: t=l[j]; b=gcoeff(A,t,k);
! 2649: if (!signe(b)) continue;
! 2650:
! 2651: ZV_elem(b,gcoeff(A,t,j), A,U,k,j);
! 2652: d = gcoeff(A,t,j);
! 2653: if (signe(d) < 0)
! 2654: {
! 2655: ZV_neg((GEN)A[j]);
! 2656: ZV_neg((GEN)U[j]);
! 2657: d = gcoeff(A,t,j);
! 2658: }
! 2659: for (j1=1; j1<j; j1++)
! 2660: {
! 2661: if (!l[j1]) continue;
! 2662: q = truedvmdii(gcoeff(A,t,j1),d,NULL);
! 2663: if (!signe(q)) continue;
! 2664:
! 2665: q = negi(q);
! 2666: A[j1] = (long)ZV_lincomb(gun,q,(GEN)A[j1],(GEN)A[j]);
! 2667: U[j1] = (long)ZV_lincomb(gun,q,(GEN)U[j1],(GEN)U[j]);
! 2668: }
! 2669: }
! 2670: t=m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
! 2671: if (t)
! 2672: {
! 2673: p = gcoeff(A,t,k);
! 2674: for (i=t-1; i; i--)
! 2675: {
! 2676: q = gcoeff(A,i,k);
! 2677: if (signe(q) && absi_cmp(p,q) > 0) { p=q; t=i; }
! 2678: }
! 2679: perm[++r]=l[k]=t; c[t]=k;
! 2680: if (signe(p) < 0)
! 2681: {
! 2682: ZV_neg((GEN)A[k]);
! 2683: ZV_neg((GEN)U[k]);
! 2684: p = gcoeff(A,t,k);
! 2685: }
! 2686: for (j=1; j<k; j++)
! 2687: {
! 2688: if (!l[j]) continue;
! 2689: q = truedvmdii(gcoeff(A,t,j),p,NULL);
! 2690: if (!signe(q)) continue;
! 2691:
! 2692: q = negi(q);
! 2693: A[j] = (long)ZV_lincomb(gun,q,(GEN)A[j],(GEN)A[k]);
! 2694: U[j] = (long)ZV_lincomb(gun,q,(GEN)U[j],(GEN)U[k]);
! 2695: }
! 2696: }
! 2697: if (low_stack(lim, stack_lim(av1,1)))
! 2698: {
! 2699: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
! 2700: if (DEBUGMEM>1) err(warnmem,"hnfperm");
! 2701: gerepilemany(av1,gptr,2);
! 2702: }
! 2703: }
! 2704: if (r < m)
! 2705: {
! 2706: for (i=1,k=r; i<=m; i++)
! 2707: if (c[i]==0) perm[++k] = i;
! 2708: }
! 2709:
! 2710: /* We have A0*U=A, U in Gl(n,Z)
! 2711: * basis for Im(A): columns of A s.t l[j]>0 (r cols)
! 2712: * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols)
! 2713: */
! 2714: tetpil=avma; y=cgetg(4,t_VEC);
! 2715: p=cgetg(r+1,t_MAT); u=cgetg(n+1,t_MAT);
! 2716: for (t=1,k=r,j=1; j<=n; j++)
! 2717: if (l[j])
! 2718: {
! 2719: q=cgetg(m+1,t_COL); p[k]=(long)q;
! 2720: for (i=1; i<=m; i++) q[i]=lcopy(gcoeff(A,perm[m-i+1],j));
! 2721: u[k+n-r]=lcopy((GEN)U[j]);
! 2722: k--;
! 2723: }
! 2724: else u[t++]=lcopy((GEN)U[j]);
! 2725: y[1]=(long)p; y[2]=(long)u;
! 2726: q = cgetg(m+1,t_VEC); y[3]=(long)q;
! 2727: for (i=1; i<=m; i++) q[m-i+1]=lstoi(perm[i]);
! 2728: return gerepile(av,tetpil,y);
! 2729: }
! 2730:
! 2731: /*====================================================================
! 2732: * Forme Normale d'Hermite (Version par colonnes 31/01/94)
! 2733: *====================================================================*/
! 2734: GEN
! 2735: hnfall_i(GEN A, GEN *ptB, long remove)
! 2736: {
! 2737: GEN B,c,h,x,p,a;
! 2738: long m,n,r,i,j,k,li,av=avma,av1,tetpil,lim;
! 2739: GEN *gptr[2];
! 2740:
! 2741: if (typ(A)!=t_MAT) err(typeer,"hnfall");
! 2742: n=lg(A)-1;
! 2743: if (!n)
! 2744: {
! 2745: if (ptB) *ptB = cgetg(1,t_MAT);
! 2746: return cgetg(1,t_MAT);
! 2747: }
! 2748: m = lg(A[1])-1;
! 2749: c = cgetg(m+1,t_VECSMALL); for (i=1; i<=m; i++) c[i]=0;
! 2750: h = cgetg(n+1,t_VECSMALL); for (j=1; j<=n; j++) h[j]=m;
! 2751: av1 = avma; lim = stack_lim(av1,1);
! 2752: A = dummycopy(A);
! 2753: B = ptB? idmat(n): NULL;
! 2754: r = n+1;
! 2755: for (li=m; li; li--)
! 2756: {
! 2757: for (j=1; j<r; j++)
! 2758: {
! 2759: for (i=h[j]; i>li; i--)
! 2760: {
! 2761: a = gcoeff(A,i,j);
! 2762: if (!signe(a)) continue;
! 2763:
! 2764: k = c[i]; /* zero a = Aij using Aik */
! 2765: ZV_elem(a,gcoeff(A,i,k), A,B,j,k);
! 2766: ZM_reduce(A,B, i,k);
! 2767: if (low_stack(lim, stack_lim(av1,1)))
! 2768: {
! 2769: tetpil = avma;
! 2770: A = ZM_copy(A); gptr[0]=&A;
! 2771: if (B) { B = ZM_copy(B); gptr[1]=&B; }
! 2772: if (DEBUGMEM>1) err(warnmem,"hnfall[1], li = %ld", li);
! 2773: gerepilemanysp(av1,tetpil,gptr,B? 2: 1);
! 2774: }
! 2775: }
! 2776: x = gcoeff(A,li,j); if (signe(x)) break;
! 2777: h[j] = li-1;
! 2778: }
! 2779: if (j == r) continue;
! 2780: r--;
! 2781: if (j < r) /* A[j] != 0 */
! 2782: {
! 2783: swap(A[j], A[r]);
! 2784: if (B) swap(B[j], B[r]);
! 2785: h[j]=h[r]; h[r]=li; c[li]=r;
! 2786: }
! 2787: p = gcoeff(A,li,r);
! 2788: if (signe(p) < 0)
! 2789: {
! 2790: ZV_neg((GEN)A[r]);
! 2791: if (B) ZV_neg((GEN)B[r]);
! 2792: p = gcoeff(A,li,r);
! 2793: }
! 2794: ZM_reduce(A,B, li,r);
! 2795: if (low_stack(lim, stack_lim(av1,1)))
! 2796: {
! 2797: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&B;
! 2798: if (DEBUGMEM>1) err(warnmem,"hnfall[2], li = %ld", li);
! 2799: gerepilemany(av1,gptr,B? 2: 1);
! 2800: }
! 2801: }
! 2802: if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");
! 2803: r--; /* first r cols are in the image the n-r (independent) last ones */
! 2804: for (j=1; j<=r; j++)
! 2805: for (i=h[j]; i; i--)
! 2806: {
! 2807: a = gcoeff(A,i,j);
! 2808: if (!signe(a)) continue;
! 2809:
! 2810: k = c[i];
! 2811: ZV_elem(a,gcoeff(A,i,k), A,B, j,k);
! 2812: ZM_reduce(A,B, i,k);
! 2813: if (low_stack(lim, stack_lim(av1,1)))
! 2814: {
! 2815: tetpil = avma;
! 2816: A = ZM_copy(A); gptr[0]=&A;
! 2817: if (B) { B = ZM_copy(B); gptr[1]=&B; }
! 2818: if (DEBUGMEM>1) err(warnmem,"hnfall[3], j = %ld", j);
! 2819: gerepilemanysp(av1,tetpil,gptr, B? 2: 1);
! 2820: }
! 2821: }
! 2822: if (DEBUGLEVEL>5) fprintferr("\n");
! 2823: /* remove the first r columns */
! 2824: if (remove) { A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1); }
! 2825: gptr[0] = &A; gptr[1] = &B;
! 2826: gerepilemany(av, gptr, B? 2: 1);
! 2827: if (B) *ptB = B;
! 2828: return A;
! 2829: }
! 2830:
! 2831: GEN
! 2832: hnfall0(GEN A, long remove)
! 2833: {
! 2834: GEN B, z = cgetg(3, t_VEC);
! 2835: z[1] = (long)hnfall_i(A, &B, remove);
! 2836: z[2] = (long)B; return z;
! 2837: }
! 2838:
! 2839: GEN
! 2840: hnfall(GEN x) {return hnfall0(x,1);}
! 2841:
! 2842: /***************************************************************/
! 2843: /** **/
! 2844: /** SMITH NORMAL FORM REDUCTION **/
! 2845: /** **/
! 2846: /***************************************************************/
! 2847:
! 2848: static GEN
! 2849: col_mul(GEN x, GEN c)
! 2850: {
! 2851: long s = signe(x);
! 2852: GEN xc = NULL;
! 2853: if (s)
! 2854: {
! 2855: if (!is_pm1(x)) xc = gmul(x,c);
! 2856: else xc = (s>0)? c: gneg_i(c);
! 2857: }
! 2858: return xc;
! 2859: }
! 2860:
! 2861: static void
! 2862: do_zero(GEN x)
! 2863: {
! 2864: long i, lx = lg(x);
! 2865: for (i=1; i<lx; i++) x[i] = zero;
! 2866: }
! 2867:
! 2868: /* c1 <-- u.c1 + v.c2; c2 <-- a.c2 - b.c1 */
! 2869: static void
! 2870: update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
! 2871: {
! 2872: GEN p1,p2;
! 2873:
! 2874: u = col_mul(u,*c1);
! 2875: v = col_mul(v,*c2);
! 2876: if (u) p1 = v? gadd(u,v): u;
! 2877: else p1 = v? v: (GEN)NULL;
! 2878:
! 2879: a = col_mul(a,*c2);
! 2880: b = col_mul(gneg_i(b),*c1);
! 2881: if (a) p2 = b? gadd(a,b): a;
! 2882: else p2 = b? b: (GEN)NULL;
! 2883:
! 2884: if (!p1) do_zero(*c1); else *c1 = p1;
! 2885: if (!p2) do_zero(*c2); else *c2 = p2;
! 2886: }
! 2887:
! 2888: static GEN
! 2889: trivsmith(long all)
! 2890: {
! 2891: GEN z;
! 2892: if (!all) return cgetg(1,t_VEC);
! 2893: z=cgetg(4,t_VEC);
! 2894: z[1]=lgetg(1,t_MAT);
! 2895: z[2]=lgetg(1,t_MAT);
! 2896: z[3]=lgetg(1,t_MAT); return z;
! 2897: }
! 2898:
! 2899: /* Return the smith normal form d of matrix x. If all != 0 return [d,u,v],
! 2900: * where d = u.x.v
! 2901: */
! 2902: static GEN
! 2903: smithall(GEN x, long all)
! 2904: {
! 2905: long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim;
! 2906: GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys;
! 2907:
! 2908: if (typ(x)!=t_MAT) err(typeer,"smithall");
! 2909: if (DEBUGLEVEL>=9) outerr(x);
! 2910: n=lg(x)-1;
! 2911: if (!n) return trivsmith(all);
! 2912: if (lg(x[1]) != n+1) err(mattype1,"smithall");
! 2913: for (i=1; i<=n; i++)
! 2914: for (j=1; j<=n; j++)
! 2915: if (typ(coeff(x,i,j)) != t_INT)
! 2916: err(talker,"non integral matrix in smithall");
! 2917:
! 2918: av = avma; lim = stack_lim(av,1);
! 2919: x = dummycopy(x); mdet = detint(x);
! 2920: if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }
! 2921: else
! 2922: {
! 2923: if (signe(mdet))
! 2924: {
! 2925: p1=hnfmod(x,mdet);
! 2926: if (all) { ml=idmat(n); mr=gauss(x,p1); }
! 2927: }
! 2928: else
! 2929: {
! 2930: if (all)
! 2931: {
! 2932: p1 = hnfall0(x,0);
! 2933: ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1];
! 2934: }
! 2935: else
! 2936: p1 = hnf0(x,0);
! 2937: }
! 2938: x = p1;
! 2939: }
! 2940: p1=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) p1[i]=lnegi(gcoeff(x,i,i));
! 2941: p2=sindexsort(p1); ys=cgetg(n+1,t_MAT);
! 2942: for (j=1; j<=n; j++)
! 2943: {
! 2944: p1=cgetg(n+1,t_COL); ys[j]=(long)p1;
! 2945: for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]);
! 2946: }
! 2947: x = ys;
! 2948: if (all)
! 2949: {
! 2950: p3=cgetg(n+1,t_MAT); p4=cgetg(n+1,t_MAT);
! 2951: for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }
! 2952: ml=p3; mr=p4;
! 2953: }
! 2954: if (signe(mdet))
! 2955: {
! 2956: p1 = hnfmod(x,mdet);
! 2957: if (all) mr=gmul(mr,gauss(x,p1));
! 2958: }
! 2959: else
! 2960: {
! 2961: if (all)
! 2962: {
! 2963: p1 = hnfall0(x,0);
! 2964: mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1];
! 2965: }
! 2966: else
! 2967: p1 = hnf0(x,0);
! 2968: }
! 2969: x=p1; mun = negi(gun);
! 2970:
! 2971: if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();}
! 2972: for (i=n; i>=2; i--)
! 2973: {
! 2974: if (DEBUGLEVEL>7) {fprintferr("\ni = %ld: ",i);flusherr();}
! 2975: for(;;)
! 2976: {
! 2977: c = 0;
! 2978: for (j=i-1; j>=1; j--)
! 2979: {
! 2980: p1=gcoeff(x,i,j); s1 = signe(p1);
! 2981: if (s1)
! 2982: {
! 2983: p2=gcoeff(x,i,i);
! 2984: if (!absi_cmp(p1,p2))
! 2985: {
! 2986: s2=signe(p2);
! 2987: if (s1 == s2) { d=p1; u=gun; p4=gun; }
! 2988: else
! 2989: {
! 2990: if (s2>0) { u = gun; p4 = mun; }
! 2991: else { u = mun; p4 = gun; }
! 2992: d=(s1>0)? p1: absi(p1);
! 2993: }
! 2994: v = gzero; p3 = u;
! 2995: }
! 2996: else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
! 2997: for (k=1; k<=i; k++)
! 2998: {
! 2999: b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));
! 3000: coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),
! 3001: mulii(p4,gcoeff(x,k,i)));
! 3002: coeff(x,k,i)=(long)b;
! 3003: }
! 3004: if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
! 3005: if (low_stack(lim, stack_lim(av,1)))
! 3006: {
! 3007: if (DEBUGMEM>1) err(warnmem,"[1]: smithall");
! 3008: if (all)
! 3009: {
! 3010: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 3011: gerepilemany(av,gptr,3);
! 3012: }
! 3013: else x=gerepileupto(av, ZM_copy(x));
! 3014: }
! 3015: }
! 3016: }
! 3017: if (DEBUGLEVEL>=8) {fprintferr("; ");flusherr();}
! 3018: for (j=i-1; j>=1; j--)
! 3019: {
! 3020: p1=gcoeff(x,j,i); s1 = signe(p1);
! 3021: if (s1)
! 3022: {
! 3023: p2=gcoeff(x,i,i);
! 3024: if (!absi_cmp(p1,p2))
! 3025: {
! 3026: s2 = signe(p2);
! 3027: if (s1 == s2) { d=p1; u=gun; p4=gun; }
! 3028: else
! 3029: {
! 3030: if (s2>0) { u = gun; p4 = mun; }
! 3031: else { u = mun; p4 = gun; }
! 3032: d=(s1>0)? p1: absi(p1);
! 3033: }
! 3034: v = gzero; p3 = u;
! 3035: }
! 3036: else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
! 3037: for (k=1; k<=i; k++)
! 3038: {
! 3039: b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
! 3040: coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),
! 3041: mulii(p4,gcoeff(x,i,k)));
! 3042: coeff(x,i,k)=(long)b;
! 3043: }
! 3044: if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
! 3045: c = 1;
! 3046: }
! 3047: }
! 3048: if (!c)
! 3049: {
! 3050: b=gcoeff(x,i,i); fl=1;
! 3051: if (signe(b))
! 3052: {
! 3053: for (k=1; k<i && fl; k++)
! 3054: for (l=1; l<i && fl; l++)
! 3055: fl = (int)!signe(resii(gcoeff(x,k,l),b));
! 3056: /* cast to (int) necessary for gcc-2.95 on sparcv9-64 (IS) */
! 3057: if (!fl)
! 3058: {
! 3059: k--;
! 3060: for (l=1; l<=i; l++)
! 3061: coeff(x,i,l)=laddii(gcoeff(x,i,l),gcoeff(x,k,l));
! 3062: if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
! 3063: }
! 3064: }
! 3065: if (fl) break;
! 3066: }
! 3067: if (low_stack(lim, stack_lim(av,1)))
! 3068: {
! 3069: if (DEBUGMEM>1) err(warnmem,"[2]: smithall");
! 3070: if (all)
! 3071: {
! 3072: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 3073: gerepilemany(av,gptr,3);
! 3074: }
! 3075: else x=gerepileupto(av,ZM_copy(x));
! 3076: }
! 3077: }
! 3078: }
! 3079: if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();}
! 3080: if (all)
! 3081: {
! 3082: for (k=1; k<=n; k++)
! 3083: if (signe(gcoeff(x,k,k))<0)
! 3084: { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lnegi(gcoeff(x,k,k)); }
! 3085: tetpil=avma; z=cgetg(4,t_VEC);
! 3086: z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
! 3087: return gerepile(av,tetpil,z);
! 3088: }
! 3089: tetpil=avma; z=cgetg(n+1,t_VEC); j=n;
! 3090: for (k=n; k; k--)
! 3091: if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));
! 3092: for ( ; j; j--) z[j]=zero;
! 3093: return gerepile(av,tetpil,z);
! 3094: }
! 3095:
! 3096: GEN
! 3097: smith(GEN x) { return smithall(x,0); }
! 3098:
! 3099: GEN
! 3100: smith2(GEN x) { return smithall(x,1); }
! 3101:
! 3102: /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
! 3103: GEN
! 3104: smithclean(GEN z)
! 3105: {
! 3106: long i,j,l,c;
! 3107: GEN u,v,y,d,p1;
! 3108:
! 3109: if (typ(z) != t_VEC) err(typeer,"smithclean");
! 3110: l = lg(z); if (l == 1) return cgetg(1,t_VEC);
! 3111: u=(GEN)z[1];
! 3112: if (l != 4 || typ(u) != t_MAT)
! 3113: {
! 3114: if (typ(u) != t_INT) err(typeer,"smithclean");
! 3115: for (c=1; c<l; c++)
! 3116: if (gcmp1((GEN)z[c])) break;
! 3117: return gcopy_i(z, c);
! 3118: }
! 3119: v=(GEN)z[2]; d=(GEN)z[3]; l = lg(d);
! 3120: for (c=1; c<l; c++)
! 3121: if (gcmp1(gcoeff(d,c,c))) break;
! 3122:
! 3123: y=cgetg(4,t_VEC);
! 3124: y[1]=(long)(p1 = cgetg(l,t_MAT));
! 3125: for (i=1; i<l; i++) p1[i] = (long)gcopy_i((GEN)u[i], c);
! 3126: y[2]=(long)gcopy_i(v, c);
! 3127: y[3]=(long)(p1 = cgetg(c,t_MAT));
! 3128: for (i=1; i<c; i++)
! 3129: {
! 3130: GEN p2 = cgetg(c,t_COL); p1[i] = (long)p2;
! 3131: for (j=1; j<c; j++) p2[j] = i==j? lcopy(gcoeff(d,i,i)): zero;
! 3132: }
! 3133: return y;
! 3134: }
! 3135:
! 3136: static GEN
! 3137: gsmithall(GEN x,long all)
! 3138: {
! 3139: long av,tetpil,i,j,k,l,c,n, lim;
! 3140: GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr;
! 3141:
! 3142: if (typ(x)!=t_MAT) err(typeer,"gsmithall");
! 3143: n=lg(x)-1;
! 3144: if (!n) return trivsmith(all);
! 3145: if (lg(x[1]) != n+1) err(mattype1,"gsmithall");
! 3146: av = avma; lim = stack_lim(av,1);
! 3147: x = dummycopy(x);
! 3148: if (all) { ml=idmat(n); mr=idmat(n); }
! 3149: for (i=n; i>=2; i--)
! 3150: {
! 3151: do
! 3152: {
! 3153: c=0;
! 3154: for (j=i-1; j>=1; j--)
! 3155: {
! 3156: p1=gcoeff(x,i,j);
! 3157: if (signe(p1))
! 3158: {
! 3159: p2=gcoeff(x,i,i);
! 3160: if (!signe(p2))
! 3161: {
! 3162: p3 = gzero; p4 = gun; u = gzero; v = gun;
! 3163: }
! 3164: else
! 3165: {
! 3166: v = gdiventres(p1,p2);
! 3167: if (gcmp0((GEN)v[2]))
! 3168: { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
! 3169: else
! 3170: { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
! 3171: }
! 3172: for (k=1; k<=i; k++)
! 3173: {
! 3174: b=gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
! 3175: coeff(x,k,j)=lsub(gmul(p3,gcoeff(x,k,j)),gmul(p4,gcoeff(x,k,i)));
! 3176: coeff(x,k,i)=(long)b;
! 3177: }
! 3178: if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
! 3179: }
! 3180: }
! 3181: for (j=i-1; j>=1; j--)
! 3182: {
! 3183: p1=gcoeff(x,j,i);
! 3184: if (signe(p1))
! 3185: {
! 3186: p2 = gcoeff(x,i,i);
! 3187: if (!signe(p2))
! 3188: {
! 3189: p3 = gzero; p4 = gun; u = gzero; v = gun;
! 3190: }
! 3191: else
! 3192: {
! 3193: v = gdiventres(p1,p2);
! 3194: if (gcmp0((GEN)v[2]))
! 3195: { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
! 3196: else
! 3197: { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
! 3198: }
! 3199: for (k=1; k<=i; k++)
! 3200: {
! 3201: b=gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
! 3202: coeff(x,j,k)=lsub(gmul(p3,gcoeff(x,j,k)),gmul(p4,gcoeff(x,i,k)));
! 3203: coeff(x,i,k)=(long)b;
! 3204: }
! 3205: if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
! 3206: c = 1;
! 3207: }
! 3208: }
! 3209: if (!c)
! 3210: {
! 3211: b = gcoeff(x,i,i);
! 3212: if (signe(b))
! 3213: for (k=1; k<i; k++)
! 3214: for (l=1; l<i; l++)
! 3215: if (signe(gmod(gcoeff(x,k,l),b)))
! 3216: {
! 3217: for (l=1; l<=i; l++)
! 3218: coeff(x,i,l) = ladd(gcoeff(x,i,l),gcoeff(x,k,l));
! 3219: if (all) ml[i] = ladd((GEN)ml[i],(GEN)ml[k]);
! 3220: k = l = i; c = 1;
! 3221: }
! 3222: }
! 3223: if (low_stack(lim, stack_lim(av,1)))
! 3224: {
! 3225: if (DEBUGMEM>1) err(warnmem,"[5]: smithall");
! 3226: if (all)
! 3227: {
! 3228: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 3229: gerepilemany(av,gptr,3);
! 3230: }
! 3231: else x=gerepilecopy(av,x);
! 3232: }
! 3233: }
! 3234: while (c);
! 3235: }
! 3236: if (all)
! 3237: {
! 3238: for (k=1; k<=n; k++)
! 3239: if (signe(gcoeff(x,k,k)) < 0)
! 3240: { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lneg(gcoeff(x,k,k)); }
! 3241: tetpil=avma; z=cgetg(4,t_VEC);
! 3242: z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
! 3243: }
! 3244: else
! 3245: {
! 3246: tetpil=avma; z=cgetg(n+1,t_VEC);
! 3247: for (j=0,k=1; k<=n; k++) if (!signe(gcoeff(x,k,k))) z[++j]=zero;
! 3248: for (k=1; k<=n; k++)
! 3249: if (signe(p1=gcoeff(x,k,k))) z[++j]=(long)gabs(p1,0);
! 3250: }
! 3251: return gerepile(av,tetpil,z);
! 3252: }
! 3253:
! 3254: GEN
! 3255: matsnf0(GEN x,long flag)
! 3256: {
! 3257: long av = avma;
! 3258: if (flag > 7) err(flagerr,"matsnf");
! 3259: if (typ(x) == t_VEC && flag & 4) return smithclean(x);
! 3260: if (flag & 2) x = gsmithall(x,flag & 1);
! 3261: else x = smithall(x, flag & 1);
! 3262: if (flag & 4) x = smithclean(x);
! 3263: return gerepileupto(av, x);
! 3264: }
! 3265:
! 3266: GEN
! 3267: gsmith(GEN x) { return gsmithall(x,0); }
! 3268:
! 3269: GEN
! 3270: gsmith2(GEN x) { return gsmithall(x,1); }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to alglin2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath