Annotation of OpenXM_contrib/pari/src/basemath/alglin2.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** LINEAR ALGEBRA **/
! 4: /** (second part) **/
! 5: /** **/
! 6: /********************************************************************/
! 7: /* $Id: alglin2.c,v 1.3 1999/09/23 17:50:55 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: /*******************************************************************/
! 11: /* */
! 12: /* CHARACTERISTIC POLYNOMIAL */
! 13: /* */
! 14: /*******************************************************************/
! 15:
! 16: GEN
! 17: charpoly0(GEN x, int v, long flag)
! 18: {
! 19: if (v<0) v = 0;
! 20: switch(flag)
! 21: {
! 22: case 0: return caradj(x,v,0);
! 23: case 1: return caract(x,v);
! 24: case 2: return carhess(x,v);
! 25: }
! 26: err(flagerr,"charpoly"); return NULL; /* not reached */
! 27: }
! 28:
! 29:
! 30: static GEN
! 31: caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,GEN*))
! 32: {
! 33: long av = avma, d;
! 34: GEN p1, p2 = leading_term(p);
! 35:
! 36: if (!signe(x)) return gpowgs(polx[v], lgef(p)-3);
! 37: x = gneg_i(x); x[2] = ladd((GEN)x[2], polx[MAXVARN]);
! 38: p1=subres_f(p, x, NULL);
! 39: if (varn(p1)==MAXVARN) setvarn(p1,v); else p1=gsubst(p1,MAXVARN,polx[v]);
! 40:
! 41: if (!gcmp1(p2) && (d=lgef(x)-3) > 0) p1 = gdiv(p1, gpuigs(p2,d));
! 42: return gerepileupto(av,p1);
! 43: }
! 44:
! 45: /* return caract(Mod(x,p)) in variable v */
! 46: GEN
! 47: caract2(GEN p, GEN x, int v)
! 48: {
! 49: return caract2_i(p,x,v, subresall);
! 50: }
! 51: GEN
! 52: caractducos(GEN p, GEN x, int v)
! 53: {
! 54: return caract2_i(p,x,v, (GEN (*)(GEN,GEN,GEN*))resultantducos);
! 55: }
! 56:
! 57: /* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
! 58: static GEN
! 59: easychar(GEN x, int v, GEN *py)
! 60: {
! 61: long l,tetpil,lx;
! 62: GEN p1,p2;
! 63:
! 64: switch(typ(x))
! 65: {
! 66: case t_INT: case t_REAL: case t_INTMOD:
! 67: case t_FRAC: case t_FRACN: case t_PADIC:
! 68: p1=cgetg(4,t_POL);
! 69: p1[1]=evalsigne(1) | evallgef(4) | evalvarn(v);
! 70: p1[2]=lneg(x); p1[3]=un;
! 71: if (py)
! 72: {
! 73: p2=cgetg(2,t_MAT);
! 74: p2[1]=lgetg(2,t_COL);
! 75: coeff(p2,1,1)=lcopy(x);
! 76: *py=p2;
! 77: }
! 78: return p1;
! 79:
! 80: case t_COMPLEX: case t_QUAD:
! 81: if (py) err(typeer,"easychar");
! 82: p1=cgetg(5,t_POL);
! 83: p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v);
! 84: p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma;
! 85: p1[3]=lpile(l,tetpil,gneg(p2));
! 86: p1[4]=un; return p1;
! 87:
! 88: case t_POLMOD:
! 89: if (py) err(typeer,"easychar");
! 90: return caract2((GEN)x[1], (GEN)x[2], v);
! 91:
! 92: case t_MAT:
! 93: lx=lg(x);
! 94: if (lx==1)
! 95: {
! 96: if (py) *py = cgetg(1,t_MAT);
! 97: return polun[v];
! 98: }
! 99: if (lg(x[1]) != lx) break;
! 100: return NULL;
! 101: }
! 102: err(mattype1,"easychar");
! 103: return NULL; /* not reached */
! 104: }
! 105:
! 106: GEN
! 107: caract(GEN x, int v)
! 108: {
! 109: long n,k,l=avma,tetpil;
! 110: GEN p1,p2,p3,p4,p5,p6;
! 111:
! 112: if ((p1 = easychar(x,v,NULL))) return p1;
! 113:
! 114: p1=gzero; p2=gun;
! 115: n=lg(x)-1; if (n&1) p2=gneg_i(p2);
! 116: p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5;
! 117: p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3);
! 118: for (k=0; k<=n; k++)
! 119: {
! 120: p3=det(gsub(gscalmat(stoi(k),n), x));
! 121: p4[1]=lmul(p3,p2); p6[2]=k;
! 122: p1=gadd(p4,p1); p5[2]=(long)p6;
! 123: if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
! 124: }
! 125: p2=mpfact(n); tetpil=avma;
! 126: return gerepile(l,tetpil,gdiv((GEN) p1[1],p2));
! 127: }
! 128:
! 129: GEN
! 130: caradj0(GEN x, long v)
! 131: {
! 132: return caradj(x,v,NULL);
! 133: }
! 134:
! 135: /* Using traces: return the characteristic polynomial of x (in variable v).
! 136: * If py != NULL, the adjoint matrix is put there.
! 137: */
! 138: GEN
! 139: caradj(GEN x, long v, GEN *py)
! 140: {
! 141: long i,j,k,l,av,tetpil;
! 142: GEN p,y,t,*gptr[2];
! 143:
! 144: if ((p = easychar(x,v,py))) return p;
! 145:
! 146: l=lg(x);
! 147: if (l==1) { p=polun[v]; if (py) *py=gcopy(x); return p; }
! 148: if (l==2) { p=gsub(polx[v],gtrace(x)); if (py) *py=idmat(1); return p; }
! 149:
! 150: p=cgetg(l+2,t_POL); p[1]=evalsigne(1) | evallgef(l+2) | evalvarn(v);
! 151: av=avma; t=gtrace(x); tetpil=avma;
! 152: t=gerepile(av,tetpil,gneg(t));
! 153: p[l]=(long)t; p[l+1]=un;
! 154:
! 155: av=avma; y=cgetg(l,t_MAT);
! 156: for (j=1; j<l; j++)
! 157: {
! 158: y[j]=lgetg(l,t_COL);
! 159: for (i=1; i<l; i++)
! 160: coeff(y,i,j) = (i==j) ? ladd(gcoeff(x,i,j),t) : coeff(x,i,j);
! 161: }
! 162:
! 163: for (k=2; k<l-1; k++)
! 164: {
! 165: GEN z=gmul(x,y);
! 166:
! 167: t=gtrace(z); tetpil=avma;
! 168: t=gdivgs(t,-k); y=cgetg(l,t_MAT);
! 169: for (j=1; j<l; j++)
! 170: {
! 171: y[j]=lgetg(l,t_COL);
! 172: for (i=1;i<l;i++)
! 173: coeff(y,i,j) = (i==j)? ladd(gcoeff(z,i,i),t): lcopy(gcoeff(z,i,j));
! 174: }
! 175: gptr[0]=&y; gptr[1]=&t;
! 176: gerepilemanysp(av,tetpil,gptr,2);
! 177: p[l-k+1]=(long)t; av=avma;
! 178: }
! 179:
! 180: t=gzero;
! 181: for (i=1; i<l; i++)
! 182: t=gadd(t,gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
! 183: tetpil=avma; t=gneg(t);
! 184:
! 185: if (py)
! 186: {
! 187: *py = (l&1) ? gneg(y) : gcopy(y);
! 188: gptr[0] = &t; gptr[1] = py;
! 189: gerepilemanysp(av,tetpil,gptr,2);
! 190: p[2]=(long)t;
! 191: }
! 192: else p[2]=lpile(av,tetpil,t);
! 193: i = gvar2(p);
! 194: if (i == v) err(talker,"incorrect variable in caradj");
! 195: if (i < v) p = poleval(p, polx[v]);
! 196: return p;
! 197: }
! 198:
! 199: GEN
! 200: adj(GEN x)
! 201: {
! 202: GEN y;
! 203: caradj(x,MAXVARN,&y); return y;
! 204: }
! 205:
! 206: /*******************************************************************/
! 207: /* */
! 208: /* HESSENBERG FORM */
! 209: /* */
! 210: /*******************************************************************/
! 211:
! 212: GEN
! 213: hess(GEN x)
! 214: {
! 215: long lx=lg(x),av=avma,tetpil,m,i,j;
! 216: GEN p,p1,p2;
! 217:
! 218: if (typ(x) != t_MAT) err(mattype1,"hess");
! 219: if (lx==1) return cgetg(1,t_MAT);
! 220: if (lg(x[1])!=lx) err(mattype1,"hess");
! 221: x = dummycopy(x);
! 222:
! 223: for (m=2; m<lx-1; m++)
! 224: for (i=m+1; i<lx; i++)
! 225: {
! 226: p = gcoeff(x,i,m-1);
! 227: if (!gcmp0(p))
! 228: {
! 229: for (j=m-1; j<lx; j++)
! 230: {
! 231: p1 = gcoeff(x,i,j);
! 232: coeff(x,i,j) = coeff(x,m,j);
! 233: coeff(x,m,j) = (long)p1;
! 234: }
! 235: p1=(GEN)x[i]; x[i]=x[m]; x[m]=(long)p1;
! 236: for (i=m+1; i<lx; i++)
! 237: {
! 238: p1 = gcoeff(x,i,m-1);
! 239: if (!gcmp0(p1))
! 240: {
! 241: p1 = gdiv(p1,p); p2 = gneg_i(p1);
! 242: coeff(x,i,m-1) = zero;
! 243: for (j=m; j<lx; j++)
! 244: coeff(x,i,j) = ladd(gcoeff(x,i,j), gmul(p2,gcoeff(x,m,j)));
! 245: for (j=1; j<lx; j++)
! 246: coeff(x,j,m) = ladd(gcoeff(x,j,m), gmul(p1,gcoeff(x,j,i)));
! 247: }
! 248: }
! 249: break;
! 250: }
! 251: }
! 252: tetpil=avma; return gerepile(av,tetpil,gcopy(x));
! 253: }
! 254:
! 255: GEN
! 256: carhess(GEN x, long v)
! 257: {
! 258: long av,tetpil,lx,r,i;
! 259: GEN *y,p1,p2,p3,p4;
! 260:
! 261: if ((p1 = easychar(x,v,NULL))) return p1;
! 262:
! 263: lx=lg(x); av=avma; y = (GEN*) new_chunk(lx);
! 264: y[0]=polun[v]; p1=hess(x); p2=polx[v];
! 265: tetpil=avma;
! 266: for (r=1; r<lx; r++)
! 267: {
! 268: y[r]=gmul(y[r-1], gsub(p2,gcoeff(p1,r,r)));
! 269: p3=gun; p4=gzero;
! 270: for (i=r-1; i; i--)
! 271: {
! 272: p3=gmul(p3,gcoeff(p1,i+1,i));
! 273: p4=gadd(p4,gmul(gmul(p3,gcoeff(p1,i,r)),y[i-1]));
! 274: }
! 275: tetpil=avma; y[r]=gsub(y[r],p4);
! 276: }
! 277: return gerepile(av,tetpil,y[lx-1]);
! 278: }
! 279:
! 280: /*******************************************************************/
! 281: /* */
! 282: /* NORM */
! 283: /* */
! 284: /*******************************************************************/
! 285:
! 286: GEN
! 287: gnorm(GEN x)
! 288: {
! 289: long l,lx,i,tetpil, tx=typ(x);
! 290: GEN p1,p2,y;
! 291:
! 292: switch(tx)
! 293: {
! 294: case t_INT:
! 295: return sqri(x);
! 296:
! 297: case t_REAL:
! 298: return mulrr(x,x);
! 299:
! 300: case t_FRAC: case t_FRACN:
! 301: return gsqr(x);
! 302:
! 303: case t_COMPLEX:
! 304: l=avma; p1=gsqr((GEN) x[1]); p2=gsqr((GEN) x[2]);
! 305: tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
! 306:
! 307: case t_QUAD:
! 308: l=avma; p1=(GEN)x[1];
! 309: p2=gmul((GEN) p1[2], gsqr((GEN) x[3]));
! 310: p1 = gcmp0((GEN) p1[3])? gsqr((GEN) x[2])
! 311: : gmul((GEN) x[2], gadd((GEN) x[2],(GEN) x[3]));
! 312: tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));
! 313:
! 314: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
! 315: l=avma; p1=gmul(gconj(x),x); tetpil=avma;
! 316: return gerepile(l,tetpil,greal(p1));
! 317:
! 318: case t_POLMOD:
! 319: p1=(GEN)x[1]; p2=leading_term(p1);
! 320: if (gcmp1(p2) || gcmp0((GEN) x[2])) return subres(p1,(GEN) x[2]);
! 321: l=avma; y=subres(p1,(GEN)x[2]); p1=gpuigs(p2,lgef(x[2])-3);
! 322: tetpil=avma; return gerepile(l,tetpil,gdiv(y,p1));
! 323:
! 324: case t_VEC: case t_COL: case t_MAT:
! 325: lx=lg(x); y=cgetg(lx,tx);
! 326: for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]);
! 327: return y;
! 328: }
! 329: err(typeer,"gnorm");
! 330: return NULL; /* not reached */
! 331: }
! 332:
! 333: GEN
! 334: gnorml2(GEN x)
! 335: {
! 336: GEN y;
! 337: long i,tx=typ(x),lx,av;
! 338:
! 339: if (! is_matvec_t(tx)) return gnorm(x);
! 340: lx=lg(x); if (lx==1) return gzero;
! 341:
! 342: av=avma; y = gnorml2((GEN) x[1]);
! 343: for (i=2; i<lx; i++)
! 344: y = gadd(gnorml2((GEN) x[i]), y);
! 345: return gerepileupto(av,y);
! 346: }
! 347:
! 348: /*******************************************************************/
! 349: /* */
! 350: /* CONJUGATION */
! 351: /* */
! 352: /*******************************************************************/
! 353:
! 354: GEN
! 355: gconj(GEN x)
! 356: {
! 357: long lx,i,tx=typ(x);
! 358: GEN z;
! 359:
! 360: switch(tx)
! 361: {
! 362: case t_INT: case t_REAL: case t_INTMOD:
! 363: case t_FRAC: case t_FRACN: case t_PADIC:
! 364: return gcopy(x);
! 365:
! 366: case t_COMPLEX:
! 367: z=cgetg(3,t_COMPLEX);
! 368: z[1]=lcopy((GEN) x[1]);
! 369: z[2]=lneg((GEN) x[2]);
! 370: break;
! 371:
! 372: case t_QUAD:
! 373: z=cgetg(4,t_QUAD);
! 374: copyifstack(x[1],z[1]);
! 375: z[2]=gcmp0(gmael(x,1,3))? lcopy((GEN) x[2])
! 376: : ladd((GEN) x[2],(GEN) x[3]);
! 377: z[3]=lneg((GEN) x[3]);
! 378: break;
! 379:
! 380: case t_POL:
! 381: lx=lgef(x); z=cgetg(lx,tx); z[1]=x[1];
! 382: for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 383: break;
! 384:
! 385: case t_SER:
! 386: lx=lg(x); z=cgetg(lx,tx); z[1]=x[1];
! 387: for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 388: break;
! 389:
! 390: case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 391: lx=lg(x); z=cgetg(lx,tx);
! 392: for (i=1; i<lx; i++) z[i]=lconj((GEN) x[i]);
! 393: break;
! 394:
! 395: default:
! 396: err(typeer,"gconj");
! 397: return NULL; /* not reached */
! 398: }
! 399: return z;
! 400: }
! 401:
! 402: GEN
! 403: conjvec(GEN x,long prec)
! 404: {
! 405: long lx,s,av,tetpil,i,tx=typ(x);
! 406: GEN z,y,p1,p2,p;
! 407:
! 408: switch(tx)
! 409: {
! 410: case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN:
! 411: z=cgetg(2,t_COL); z[1]=lcopy(x); break;
! 412:
! 413: case t_COMPLEX: case t_QUAD:
! 414: z=cgetg(3,t_COL); z[1]=lcopy(x); z[2]=lconj(x); break;
! 415:
! 416: case t_VEC: case t_COL:
! 417: lx=lg(x); z=cgetg(lx,t_MAT);
! 418: for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec);
! 419: s=lg(z[1]);
! 420: for (i=2; i<lx; i++)
! 421: if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");
! 422: break;
! 423:
! 424: case t_POLMOD:
! 425: y=(GEN)x[1]; lx=lgef(y);
! 426: if (lx<=3) return cgetg(1,t_COL);
! 427: av=avma; p=NULL;
! 428: for (i=2; i<lx; i++)
! 429: {
! 430: tx=typ(y[i]);
! 431: if (tx==t_INTMOD) p=gmael(y,i,1);
! 432: else
! 433: if (tx != t_INT && ! is_frac_t(tx))
! 434: err(polrationer,"conjvec");
! 435: }
! 436: if (!p)
! 437: {
! 438: p1=roots(y,prec); x = (GEN)x[2]; tetpil=avma;
! 439: z=cgetg(lx-2,t_COL);
! 440: for (i=1; i<=lx-3; i++)
! 441: {
! 442: p2=(GEN)p1[i]; if (gcmp0((GEN)p2[2])) p2 = (GEN)p2[1];
! 443: z[i] = (long)poleval(x, p2);
! 444: }
! 445: return gerepile(av,tetpil,z);
! 446: }
! 447: z=cgetg(lx-2,t_COL); z[1]=lcopy(x);
! 448: for (i=2; i<=lx-3; i++) z[i] = lpui((GEN) z[i-1],p,prec);
! 449: break;
! 450:
! 451: default:
! 452: err(typeer,"conjvec");
! 453: return NULL; /* not reached */
! 454: }
! 455: return z;
! 456: }
! 457:
! 458: /*******************************************************************/
! 459: /* */
! 460: /* TRACES */
! 461: /* */
! 462: /*******************************************************************/
! 463:
! 464: GEN
! 465: assmat(GEN x)
! 466: {
! 467: long lx,i,j;
! 468: GEN y,p1,p2;
! 469:
! 470: if (typ(x)!=t_POL) err(notpoler,"assmat");
! 471: if (gcmp0(x)) err(zeropoler,"assmat");
! 472:
! 473: lx=lgef(x)-2; y=cgetg(lx,t_MAT);
! 474: for (i=1; i<lx-1; i++)
! 475: {
! 476: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 477: for (j=1; j<lx; j++)
! 478: p1[j] = (j==(i+1))? un: zero;
! 479: }
! 480: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 481: if (gcmp1((GEN) x[lx+1]))
! 482: for (j=1; j<lx; j++)
! 483: p1[j] = lneg((GEN)x[j+1]);
! 484: else
! 485: {
! 486: p2 = (GEN)x[lx+1]; gnegz(p2,p2);
! 487: for (j=1; j<lx; j++)
! 488: p1[j] = ldiv((GEN)x[j+1],p2);
! 489: gnegz(p2,p2);
! 490: }
! 491: return y;
! 492: }
! 493:
! 494: GEN
! 495: gtrace(GEN x)
! 496: {
! 497: long i,l,n,tx=typ(x),lx,tetpil;
! 498: GEN y,p1,p2;
! 499:
! 500: switch(tx)
! 501: {
! 502: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 503: return gmul2n(x,1);
! 504:
! 505: case t_COMPLEX:
! 506: return gmul2n((GEN)x[1],1);
! 507:
! 508: case t_QUAD:
! 509: p1=(GEN)x[1];
! 510: if (!gcmp0((GEN) p1[3]))
! 511: {
! 512: l=avma; p2=gmul2n((GEN)x[2],1);
! 513: tetpil=avma; return gerepile(l,tetpil,gadd((GEN)x[3],p2));
! 514: }
! 515: return gmul2n((GEN)x[2],1);
! 516:
! 517: case t_POL:
! 518: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
! 519: for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 520: return y;
! 521:
! 522: case t_SER:
! 523: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 524: for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 525: return y;
! 526:
! 527: case t_POLMOD:
! 528: l=avma; n=(lgef(x[1])-4);
! 529: p1=polsym((GEN)x[1],n); p2=gzero;
! 530: for (i=0; i<=n; i++)
! 531: p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1]));
! 532: return gerepileupto(l,p2);
! 533:
! 534: case t_RFRAC: case t_RFRACN:
! 535: return gadd(x,gconj(x));
! 536:
! 537: case t_VEC: case t_COL:
! 538: lx=lg(x); y=cgetg(lx,tx);
! 539: for (i=1; i<lx; i++) y[i]=ltrace((GEN)x[i]);
! 540: return y;
! 541:
! 542: case t_MAT:
! 543: lx=lg(x); if (lx!=lg(x[1])) err(mattype1,"gtrace");
! 544: l=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);
! 545: for (i=2; i<lx-1; i++)
! 546: p1=gadd(p1,gcoeff(x,i,i));
! 547: tetpil=avma; return gerepile(l,tetpil,gadd(p1,gcoeff(x,i,i)));
! 548:
! 549: }
! 550: err(typeer,"gtrace");
! 551: return NULL; /* not reached */
! 552: }
! 553:
! 554: /* Gauss reduction for positive definite matrix a
! 555: * If a is not positive definite:
! 556: * if flag is zero: raise an error
! 557: * else: return NULL.
! 558: */
! 559: GEN
! 560: sqred1intern(GEN a,long flag)
! 561: {
! 562: GEN b,p;
! 563: long av = avma,tetpil,i,j,k, lim=stack_lim(av,1), n=lg(a);
! 564:
! 565: if (typ(a)!=t_MAT) err(typeer,"sqred1");
! 566: if (lg(a[1])!=n) err(mattype1,"sqred1");
! 567: b=cgetg(n,t_MAT);
! 568: for (j=1; j<n; j++)
! 569: {
! 570: GEN p1=cgetg(n,t_COL), p2=(GEN)a[j];
! 571:
! 572: b[j]=(long)p1;
! 573: for (i=1; i<=j; i++) p1[i] = p2[i];
! 574: for ( ; i<n ; i++) p1[i] = zero;
! 575: }
! 576: for (k=1; k<n; k++)
! 577: {
! 578: p = gcoeff(b,k,k);
! 579: if (gsigne(p)<=0) /* not positive definite */
! 580: {
! 581: if (flag) { avma=av; return NULL; }
! 582: err(talker,"not a positive definite matrix in sqred1");
! 583: }
! 584: p = ginv(p);
! 585: for (i=k+1; i<n; i++)
! 586: for (j=i; j<n; j++)
! 587: coeff(b,i,j)=lsub(gcoeff(b,i,j),
! 588: gmul(gmul(gcoeff(b,k,i),gcoeff(b,k,j)), p));
! 589: for (j=k+1; j<n; j++)
! 590: coeff(b,k,j)=lmul(gcoeff(b,k,j), p);
! 591: if (low_stack(lim, stack_lim(av,1)))
! 592: {
! 593: if (DEBUGMEM>1) err(warnmem,"sqred1");
! 594: tetpil=avma; b=gerepile(av,tetpil,gcopy(b));
! 595: }
! 596: }
! 597: tetpil=avma; return gerepile(av,tetpil,gcopy(b));
! 598: }
! 599:
! 600: GEN
! 601: sqred1(GEN a)
! 602: {
! 603: return sqred1intern(a,0);
! 604: }
! 605:
! 606: GEN
! 607: sqred3(GEN a)
! 608: {
! 609: long av=avma,tetpil,i,j,k,l, lim=stack_lim(av,1), n=lg(a);
! 610: GEN p1,b;
! 611:
! 612: if (typ(a)!=t_MAT) err(typeer,"sqred3");
! 613: if (lg(a[1])!=n) err(mattype1,"sqred3");
! 614: av=avma; b=cgetg(n,t_MAT);
! 615: for (j=1; j<n; j++)
! 616: {
! 617: p1=cgetg(n,t_COL); b[j]=(long)p1;
! 618: for (i=1; i<n; i++) p1[i]=zero;
! 619: }
! 620: for (i=1; i<n; i++)
! 621: {
! 622: for (k=1; k<i; k++)
! 623: {
! 624: p1=gzero;
! 625: for (l=1; l<k; l++)
! 626: p1=gadd(p1, gmul(gmul(gcoeff(b,l,l),gcoeff(b,k,l)), gcoeff(b,i,l)));
! 627: coeff(b,i,k)=ldiv(gsub(gcoeff(a,i,k),p1),gcoeff(b,k,k));
! 628: }
! 629: p1=gzero;
! 630: for (l=1; l<i; l++)
! 631: p1=gadd(p1, gmul(gcoeff(b,l,l), gsqr(gcoeff(b,i,l))));
! 632: coeff(b,i,k)=lsub(gcoeff(a,i,i),p1);
! 633: if (low_stack(lim, stack_lim(av,1)))
! 634: {
! 635: if (DEBUGMEM>1) err(warnmem,"sqred3");
! 636: tetpil=avma; b=gerepile(av,tetpil,gcopy(b));
! 637: }
! 638: }
! 639: tetpil=avma; return gerepile(av,tetpil,gcopy(b));
! 640: }
! 641:
! 642: /* Gauss reduction (arbitrary symetric matrix, only the part above the
! 643: * diagonal is considered). If no_signature is zero, return only the
! 644: * signature
! 645: */
! 646: static GEN
! 647: sqred2(GEN a, long no_signature)
! 648: {
! 649: GEN r,p,mun;
! 650: long av,tetpil,av1,lim,n,i,j,k,l,sp,sn,t;
! 651:
! 652: if (typ(a)!=t_MAT) err(typeer,"sqred2");
! 653: n=lg(a); if (lg(a[1]) != n) err(mattype1,"sqred2");
! 654:
! 655: av=avma; mun=negi(gun); r=new_chunk(n);
! 656: for (i=1; i<n; i++) r[i]=1;
! 657: av1=avma; lim=stack_lim(av1,1); a=dummycopy(a);
! 658: n--; t=n; sp=sn=0;
! 659: while (t)
! 660: {
! 661: k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++;
! 662: if (k<=n)
! 663: {
! 664: p=gcoeff(a,k,k); if (gsigne(p)>0) sp++; else sn++;
! 665: r[k]=0; t--;
! 666: for (j=1; j<=n; j++)
! 667: coeff(a,k,j)=r[j] ? ldiv(gcoeff(a,k,j),p) : zero;
! 668:
! 669: for (i=1; i<=n; i++) if (r[i])
! 670: for (j=1; j<=n; j++)
! 671: coeff(a,i,j) = r[j] ? lsub(gcoeff(a,i,j),
! 672: gmul(gmul(gcoeff(a,k,i),gcoeff(a,k,j)),p))
! 673: : zero;
! 674: coeff(a,k,k)=(long)p;
! 675: }
! 676: else
! 677: {
! 678: for (k=1; k<=n; k++) if (r[k])
! 679: {
! 680: l=k+1; while (l<=n && (gcmp0(gcoeff(a,k,l)) || !r[l])) l++;
! 681: if (l<=n)
! 682: {
! 683: p=gcoeff(a,k,l); r[k]=r[l]=0; sp++; sn++; t-=2;
! 684: for (i=1; i<=n; i++) if (r[i])
! 685: {
! 686: for (j=1; j<=n; j++)
! 687: coeff(a,i,j) =
! 688: r[j]? lsub(gcoeff(a,i,j),
! 689: gdiv(gadd(gmul(gcoeff(a,k,i),gcoeff(a,l,j)),
! 690: gmul(gcoeff(a,k,j),gcoeff(a,l,i))),
! 691: p))
! 692: : zero;
! 693: coeff(a,k,i)=ldiv(gadd(gcoeff(a,k,i),gcoeff(a,l,i)),p);
! 694: coeff(a,l,i)=ldiv(gsub(gcoeff(a,k,i),gcoeff(a,l,i)),p);
! 695: }
! 696: coeff(a,k,l)=un; coeff(a,l,k)=(long)mun;
! 697: coeff(a,k,k)=lmul2n(p,-1); coeff(a,l,l)=lneg(gcoeff(a,k,k));
! 698: if (low_stack(lim, stack_lim(av1,1)))
! 699: {
! 700: if(DEBUGMEM>1) err(warnmem,"sqred2");
! 701: tetpil=avma; a=gerepile(av1,tetpil,gcopy(a));
! 702: }
! 703: break;
! 704: }
! 705: }
! 706: if (k>n) break;
! 707: }
! 708: }
! 709: if (no_signature) { tetpil=avma; return gerepile(av,tetpil,gcopy(a)); }
! 710: avma=av;
! 711: a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a;
! 712: }
! 713:
! 714: GEN
! 715: sqred(GEN a) { return sqred2(a,1); }
! 716:
! 717: GEN
! 718: signat(GEN a) { return sqred2(a,0); }
! 719:
! 720: /* Diagonalization of a REAL symetric matrix. Return a 2-component vector:
! 721: * 1st comp = vector of eigenvalues
! 722: * 2nd comp = matrix of eigenvectors
! 723: */
! 724: GEN
! 725: jacobi(GEN a, long prec)
! 726: {
! 727: long de,e,e1,e2,l,n,i,j,p,q,av1,av2;
! 728: GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1;
! 729:
! 730: if (typ(a)!=t_MAT) err(mattype1,"jacobi");
! 731: ja=cgetg(3,t_VEC); l=lg(a); n=l-1;
! 732: e1=HIGHEXPOBIT-1;
! 733: lambda=cgetg(l,t_COL); ja[1]=(long)lambda;
! 734: for (j=1; j<=n; j++)
! 735: {
! 736: lambda[j]=lgetr(prec);
! 737: gaffect(gcoeff(a,j,j), (GEN)lambda[j]);
! 738: e=expo(lambda[j]); if (e<e1) e1=e;
! 739: }
! 740: r=cgetg(l,t_MAT); ja[2]=(long)r;
! 741: for (j=1; j<=n; j++)
! 742: {
! 743: r[j]=lgetg(l,t_COL);
! 744: for (i=1; i<=n; i++)
! 745: affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec)));
! 746: }
! 747: av1=avma;
! 748:
! 749: e2=-HIGHEXPOBIT;p=q=1;
! 750: c=cgetg(l,t_MAT);
! 751: for (j=1; j<=n; j++)
! 752: {
! 753: c[j]=lgetg(j,t_COL);
! 754: for (i=1; i<j; i++)
! 755: {
! 756: gaffect(gcoeff(a,i,j), (GEN)(coeff(c,i,j)=lgetr(prec)));
! 757: e=expo(gcoeff(c,i,j)); if (e>e2) { e2=e; p=i; q=j; }
! 758: }
! 759: }
! 760: a=c; unr = realun(prec);
! 761: de=bit_accuracy(prec);
! 762:
! 763: /* e1 = min des expo des coeff diagonaux
! 764: * e2 = max des expo des coeff extra-diagonaux
! 765: * Test d'arret: e2 < e1-precision
! 766: */
! 767: while (e1-e2<de)
! 768: {
! 769: /*calcul de la rotation associee dans le plan
! 770: des p et q-iemes vecteurs de base */
! 771: av2=avma;
! 772: x=divrr(subrr((GEN)lambda[q],(GEN)lambda[p]),shiftr(gcoeff(a,p,q),1));
! 773: y=mpsqrt(addrr(unr,mulrr(x,x)));
! 774: t=(gsigne(x)>0)? divrr(unr,addrr(x,y)) : divrr(unr,subrr(x,y));
! 775: c=divrr(unr,mpsqrt(addrr(unr,mulrr(t,t))));
! 776: s=mulrr(t,c); u=divrr(s,addrr(unr,c));
! 777:
! 778: /* Recalcul des transformees successives de a et de la matrice
! 779: cumulee (r) des rotations : */
! 780:
! 781: for (i=1; i<p; i++)
! 782: {
! 783: x=gcoeff(a,i,p); y=gcoeff(a,i,q);
! 784: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 785: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 786: affrr(x1,gcoeff(a,i,p)); affrr(y1,gcoeff(a,i,q));
! 787: }
! 788: for (i=p+1; i<q; i++)
! 789: {
! 790: x=gcoeff(a,p,i); y=gcoeff(a,i,q);
! 791: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 792: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 793: affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,i,q));
! 794: }
! 795: for (i=q+1; i<=n; i++)
! 796: {
! 797: x=gcoeff(a,p,i); y=gcoeff(a,q,i);
! 798: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 799: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 800: affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,q,i));
! 801: }
! 802: x=(GEN)lambda[p]; y=gcoeff(a,p,q); subrrz(x,mulrr(t,y),(GEN)lambda[p]);
! 803: x=y; y=(GEN)lambda[q]; addrrz(y,mulrr(t,x),y);
! 804: setexpo(x,expo(x)-de-1);
! 805:
! 806: for (i=1; i<=n; i++)
! 807: {
! 808: x=gcoeff(r,i,p); y=gcoeff(r,i,q);
! 809: x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
! 810: y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
! 811: affrr(x1,gcoeff(r,i,p)); affrr(y1,gcoeff(r,i,q));
! 812: }
! 813:
! 814: e2=expo(gcoeff(a,1,2)); p=1; q=2;
! 815: for (j=1; j<=n; j++)
! 816: {
! 817: for (i=1; i<j; i++)
! 818: if ((e=expo(gcoeff(a,i,j))) > e2) { e2=e; p=i; q=j; }
! 819: for (i=j+1; i<=n; i++)
! 820: if ((e=expo(gcoeff(a,j,i))) > e2) { e2=e; p=j; q=i; }
! 821: }
! 822: avma=av2;
! 823: }
! 824: avma=av1; return ja;
! 825: }
! 826:
! 827: /*************************************************************************/
! 828: /** **/
! 829: /** MATRICE RATIONNELLE --> ENTIERE **/
! 830: /** **/
! 831: /*************************************************************************/
! 832:
! 833: GEN
! 834: matrixqz0(GEN x,GEN p)
! 835: {
! 836: if (typ(p)!=t_INT) err(typeer,"matrixqz0");
! 837: if (signe(p)>=0) return matrixqz(x,p);
! 838: if (!cmpsi(-1,p)) return matrixqz2(x);
! 839: if (!cmpsi(-2,p)) return matrixqz3(x);
! 840: err(flagerr,"matrixqz"); return NULL; /* not reached */
! 841: }
! 842:
! 843: GEN
! 844: matrixqz(GEN x, GEN p)
! 845: {
! 846: long av=avma,av1,tetpil,i,j,j1,m,n,t,lim,nfact;
! 847: GEN p1,p2,p3,unmodp;
! 848:
! 849: if (typ(x)!=t_MAT) err(typeer,"matrixqz");
! 850: n=lg(x)-1; if (!n) return gcopy(x);
! 851: m=lg(x[1])-1;
! 852: if (n > m) err(talker,"more rows than columns in matrixqz");
! 853: if (n==m)
! 854: {
! 855: p1=det(x);
! 856: if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz");
! 857: avma=av; return idmat(n);
! 858: }
! 859: p1=cgetg(n+1,t_MAT);
! 860: for (j=1; j<=n; j++)
! 861: {
! 862: p2=gun; p3=(GEN)x[j];
! 863: for (i=1; i<=m; i++)
! 864: {
! 865: t=typ(p3[i]);
! 866: if (t != t_INT && !is_frac_t(t))
! 867: err(talker,"not a rational or integral matrix in matrixqz");
! 868: p2=ggcd(p2,(GEN)p3[i]);
! 869: }
! 870: p1[j]=ldiv(p3,p2);
! 871: }
! 872: x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un;
! 873:
! 874: if (gcmp0(p))
! 875: {
! 876: p1=cgetg(n+1,t_MAT); p2=gtrans(x);
! 877: for (j=1; j<=n; j++) p1[j]=p2[j];
! 878: p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1));
! 879: if (!signe(p3))
! 880: err(impl,"matrixqz when the first 2 dets are zero");
! 881: if (gcmp1(p3))
! 882: { tetpil=avma; return gerepile(av,tetpil,gcopy(x)); }
! 883:
! 884: p3=factor(p3); p1=(GEN)p3[1]; nfact=lg(p1)-1;
! 885: }
! 886: else
! 887: {
! 888: p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1;
! 889: }
! 890: av1=avma; lim=stack_lim(av1,1);
! 891: for (i=1; i<=nfact; i++)
! 892: {
! 893: p=(GEN)p1[i]; unmodp[1]=(long)p;
! 894: for(;;)
! 895: {
! 896: p2=ker(gmul(unmodp,x));
! 897: if (lg(p2)==1) break;
! 898:
! 899: p2=centerlift(p2); p3=gdiv(gmul(x,p2),p);
! 900: for (j=1; j<lg(p2); j++)
! 901: {
! 902: j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;
! 903: x[j1] = p3[j];
! 904: }
! 905: if (low_stack(lim, stack_lim(av1,1)))
! 906: {
! 907: if (DEBUGMEM>1) err(warnmem,"matrixqz");
! 908: tetpil=avma; x=gerepile(av1,tetpil,gcopy(x));
! 909: }
! 910: }
! 911: }
! 912: tetpil=avma; return gerepile(av,tetpil,gcopy(x));
! 913: }
! 914:
! 915: static GEN
! 916: matrixqz_aux(GEN x, long m, long n)
! 917: {
! 918: long av = avma, lim = stack_lim(av,1), tetpil,i,j,k,in[2];
! 919: GEN p1;
! 920:
! 921: for (i=1; i<=m; i++)
! 922: {
! 923: for(;;)
! 924: {
! 925: long fl=0;
! 926:
! 927: for (j=1; j<=n; j++)
! 928: if (!gcmp0(gcoeff(x,i,j)))
! 929: { in[fl]=j; fl++; if (fl==2) break; }
! 930: if (j>n) break;
! 931:
! 932: j=(gcmp(gabs(gcoeff(x,i,in[0]),DEFAULTPREC),
! 933: gabs(gcoeff(x,i,in[1]),DEFAULTPREC)) > 0)? in[1]: in[0];
! 934: p1=gcoeff(x,i,j);
! 935: for (k=1; k<=n; k++)
! 936: if (k!=j)
! 937: x[k]=lsub((GEN)x[k],gmul(ground(gdiv(gcoeff(x,i,k),p1)),(GEN)x[j]));
! 938: }
! 939:
! 940: j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++;
! 941: if (j<=n)
! 942: {
! 943: p1=denom(gcoeff(x,i,j));
! 944: if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]);
! 945: }
! 946: if (low_stack(lim, stack_lim(av,1)))
! 947: {
! 948: if(DEBUGMEM>1) err(warnmem,"matrixqz_aux");
! 949: tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
! 950: }
! 951: }
! 952: return hnf(x);
! 953: }
! 954:
! 955: GEN
! 956: matrixqz2(GEN x)
! 957: {
! 958: long av = avma,m,n;
! 959:
! 960: if (typ(x)!=t_MAT) err(typeer,"matrixqz2");
! 961: n=lg(x)-1; if (!n) return gcopy(x);
! 962: m=lg(x[1])-1; x=dummycopy(x);
! 963: return gerepileupto(av, matrixqz_aux(x,m,n));
! 964: }
! 965:
! 966: GEN
! 967: matrixqz3(GEN x)
! 968: {
! 969: long av=avma,av1,j,j1,k,m,n,lim;
! 970: GEN c;
! 971:
! 972: if (typ(x)!=t_MAT) err(typeer,"matrixqz3");
! 973: n=lg(x)-1; if (!n) return gcopy(x);
! 974: m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1);
! 975: for (j=1; j<=n; j++) c[j]=0;
! 976: av1=avma; lim=stack_lim(av1,1);
! 977: for (k=1; k<=m; k++)
! 978: {
! 979: j=1; while (j<=n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
! 980: if (j<=n)
! 981: {
! 982: c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));
! 983: for (j1=1; j1<=n; j1++)
! 984: if (j1!=j) x[j1]=lsub((GEN)x[j1],gmul(gcoeff(x,k,j1),(GEN)x[j]));
! 985: if (low_stack(lim, stack_lim(av1,1)))
! 986: {
! 987: long tetpil = avma;
! 988: if(DEBUGMEM>1) err(warnmem,"matrixqz3");
! 989: x=gerepile(av1,tetpil,gcopy(x));
! 990: }
! 991: }
! 992: }
! 993: return gerepileupto(av, matrixqz_aux(x,m,n));
! 994: }
! 995:
! 996: GEN
! 997: intersect(GEN x, GEN y)
! 998: {
! 999: long av,tetpil,j, lx = lg(x);
! 1000: GEN z;
! 1001:
! 1002: if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect");
! 1003: if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
! 1004:
! 1005: av=avma; z=ker(concatsp(x,y));
! 1006: for (j=lg(z)-1; j; j--) setlg(z[j],lx);
! 1007: tetpil=avma; return gerepile(av,tetpil,gmul(x,z));
! 1008: }
! 1009:
! 1010: /**************************************************************/
! 1011: /** **/
! 1012: /** HERMITE NORMAL FORM REDUCTION **/
! 1013: /** **/
! 1014: /**************************************************************/
! 1015:
! 1016: GEN
! 1017: mathnf0(GEN x, long flag)
! 1018: {
! 1019: switch(flag)
! 1020: {
! 1021: case 0: return hnf(x);
! 1022: case 1: return hnfall(x);
! 1023: case 2: return hnfhavas(x);
! 1024: case 3: return hnfperm(x);
! 1025: case 4: return hnflll(x);
! 1026: default: err(flagerr,"mathnf");
! 1027: }
! 1028: return NULL; /* not reached */
! 1029: }
! 1030:
! 1031: static GEN
! 1032: init_hnf(GEN x, GEN *denx, long *co, long *li, long *av)
! 1033: {
! 1034: if (typ(x) != t_MAT) err(typeer,"mathnf");
! 1035: *co=lg(x); if (*co==1) return NULL; /* empty matrix */
! 1036: *li=lg(x[1]); *denx=denom(x); *av=avma;
! 1037:
! 1038: if (gcmp1(*denx)) /* no denominator */
! 1039: { *denx = NULL; return dummycopy(x); }
! 1040: return gmul(x,*denx);
! 1041: }
! 1042:
! 1043: /* return c * X */
! 1044: #ifdef INLINE
! 1045: INLINE
! 1046: #endif
! 1047: GEN
! 1048: int_col_mul(GEN c, GEN X)
! 1049: {
! 1050: long i,m = lg(X);
! 1051: GEN A = new_chunk(m);
! 1052: for (i=1; i<m; i++) A[i] = lmulii(c,(GEN)X[i]);
! 1053: A[0] = X[0]; return A;
! 1054: }
! 1055:
! 1056: /* X,Y columns; u,v scalars; everybody is integral. return A = u*X + v*Y */
! 1057: GEN
! 1058: lincomb_integral(GEN u, GEN v, GEN X, GEN Y)
! 1059: {
! 1060: long av,i,lx,m;
! 1061: GEN p1,p2,A;
! 1062:
! 1063: if (!signe(u)) return int_col_mul(v,Y);
! 1064: if (!signe(v)) return int_col_mul(u,X);
! 1065: lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4;
! 1066: if (gcmp1(u))
! 1067: {
! 1068: for (i=1; i<lx; i++)
! 1069: {
! 1070: p1=(GEN)X[i]; p2=(GEN)Y[i];
! 1071: if (!signe(p1)) A[i] = lmulii(v,p2);
! 1072: else if (!signe(p2)) A[i] = licopy(p1);
! 1073: else
! 1074: {
! 1075: av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2));
! 1076: p2 = mulii(v,p2);
! 1077: avma = av; A[i] = laddii(p1,p2);
! 1078: }
! 1079: }
! 1080: }
! 1081: else
! 1082: for (i=1; i<lx; i++)
! 1083: {
! 1084: p1=(GEN)X[i]; p2=(GEN)Y[i];
! 1085: if (!signe(p1)) A[i] = lmulii(v,p2);
! 1086: else if (!signe(p2)) A[i] = lmulii(u,p1);
! 1087: else
! 1088: {
! 1089: av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2));
! 1090: p1 = mulii(u,p1);
! 1091: p2 = mulii(v,p2);
! 1092: avma = av; A[i] = laddii(p1,p2);
! 1093: }
! 1094: }
! 1095: return A;
! 1096: }
! 1097:
! 1098: GEN
! 1099: hnf_special(GEN x, long remove)
! 1100: {
! 1101: long av0, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
! 1102: GEN p1,u,v,d,denx,a,b, x2,res;
! 1103:
! 1104: if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special");
! 1105: res = cgetg(3,t_VEC);
! 1106: av0 = avma;
! 1107: x2 = (GEN)x[2];
! 1108: x = (GEN)x[1];
! 1109: x = init_hnf(x,&denx,&co,&li,&av);
! 1110: if (!x) return cgetg(1,t_MAT);
! 1111:
! 1112: lim = stack_lim(av,1);
! 1113: def=co-1; ldef=(li>co)? li-co: 0;
! 1114: if (lg(x2) != co) err(talker,"incompatible matrices in hnf_special");
! 1115: x2 = dummycopy(x2);
! 1116: for (i=li-1; i>ldef; i--)
! 1117: {
! 1118: for (j=def-1; j; j--)
! 1119: {
! 1120: while (j && !signe(gcoeff(x,i,j))) j--;
! 1121: if (!j) break;
! 1122: k = (j==1)? def: j-1;
! 1123: a = gcoeff(x,i,j);
! 1124: b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
! 1125: if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
! 1126: if (DEBUGLEVEL>5) { outerr(u); outerr(v); }
! 1127: p1 = (GEN)x[j]; b = negi(b);
! 1128: x[j] = (long)lincomb_integral(a,b, (GEN)x[k], p1);
! 1129: x[k] = (long)lincomb_integral(u,v, p1, (GEN)x[k]);
! 1130: p1 = (GEN)x2[j];
! 1131: x2[j]= ladd(gmul(a, (GEN)x2[k]), gmul(b,p1));
! 1132: x2[k] = ladd(gmul(u,p1), gmul(v, (GEN)x2[k]));
! 1133: if (low_stack(lim, stack_lim(av,1)))
! 1134: {
! 1135: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1136: if (DEBUGMEM>1) err(warnmem,"hnf_special[1]. i=%ld",i);
! 1137: gerepilemany(av,gptr,2);
! 1138: }
! 1139: }
! 1140: p1 = gcoeff(x,i,def); s = signe(p1);
! 1141: if (s)
! 1142: {
! 1143: if (s < 0)
! 1144: {
! 1145: x[def] = lneg((GEN)x[def]); p1 = gcoeff(x,i,def);
! 1146: x2[def]= lneg((GEN)x2[def]);
! 1147: }
! 1148: for (j=def+1; j<co; j++)
! 1149: {
! 1150: b = negi(gdivent(gcoeff(x,i,j),p1));
! 1151: x[j] = (long)lincomb_integral(gun,b, (GEN)x[j], (GEN)x[def]);
! 1152: x2[j]= ladd((GEN)x2[j], gmul(b, (GEN)x2[def]));
! 1153: }
! 1154: def--;
! 1155: }
! 1156: else
! 1157: if (ldef && i==ldef+1) ldef--;
! 1158: if (low_stack(lim, stack_lim(av,1)))
! 1159: {
! 1160: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1161: if (DEBUGMEM>1) err(warnmem,"hnf_special[2]. i=%ld",i);
! 1162: gerepilemany(av,gptr,2);
! 1163: }
! 1164: }
! 1165: if (remove)
! 1166: { /* remove null columns */
! 1167: for (i=1,j=1; j<co; j++)
! 1168: if (!gcmp0((GEN)x[j]))
! 1169: {
! 1170: x[i] = x[j];
! 1171: x2[i] = x2[j]; i++;
! 1172: }
! 1173: setlg(x,i);
! 1174: setlg(x2,i);
! 1175: }
! 1176: tetpil=avma;
! 1177: x = denx? gdiv(x,denx): gcopy(x);
! 1178: x2 = gcopy(x2);
! 1179: {
! 1180: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
! 1181: gerepilemanysp(av0,tetpil,gptr,2);
! 1182: }
! 1183: res[1] = (long)x;
! 1184: res[2] = (long)x2;
! 1185: return res;
! 1186: }
! 1187:
! 1188: GEN
! 1189: hnf0(GEN x, long remove) /* remove: throw away lin.dep.columns, GN */
! 1190: {
! 1191: long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
! 1192: GEN p1,u,v,d,denx,a,b;
! 1193:
! 1194: if (typ(x) == t_VEC) return hnf_special(x,remove);
! 1195: x = init_hnf(x,&denx,&co,&li,&av);
! 1196: if (!x) return cgetg(1,t_MAT);
! 1197:
! 1198: lim = stack_lim(av,1);
! 1199: def=co-1; ldef=(li>co)? li-co: 0;
! 1200: for (i=li-1; i>ldef; i--)
! 1201: {
! 1202: for (j=def-1; j; j--)
! 1203: {
! 1204: while (j && !signe(gcoeff(x,i,j))) j--;
! 1205: if (!j) break;
! 1206: k = (j==1)? def: j-1;
! 1207: a = gcoeff(x,i,j);
! 1208: b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
! 1209: if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
! 1210: if (DEBUGLEVEL>5) { outerr(u); outerr(v); }
! 1211: p1 = (GEN)x[j];
! 1212: x[j] = (long)lincomb_integral(a,negi(b), (GEN)x[k],p1);
! 1213: x[k] = (long)lincomb_integral(u,v, p1,(GEN)x[k]);
! 1214: if (low_stack(lim, stack_lim(av,1)))
! 1215: {
! 1216: if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i);
! 1217: tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
! 1218: }
! 1219: }
! 1220: p1 = gcoeff(x,i,def); s = signe(p1);
! 1221: if (s)
! 1222: {
! 1223: if (s < 0) { x[def] = lneg((GEN)x[def]); p1 = gcoeff(x,i,def); }
! 1224: for (j=def+1; j<co; j++)
! 1225: {
! 1226: b = negi(gdivent(gcoeff(x,i,j),p1));
! 1227: x[j] = (long)lincomb_integral(gun,b, (GEN)x[j], (GEN)x[def]);
! 1228: }
! 1229: def--;
! 1230: }
! 1231: else
! 1232: if (ldef && i==ldef+1) ldef--;
! 1233: if (low_stack(lim, stack_lim(av,1)))
! 1234: {
! 1235: if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i);
! 1236: tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
! 1237: }
! 1238: }
! 1239: if (remove)
! 1240: { /* remove null columns */
! 1241: for (i=1,j=1; j<co; j++)
! 1242: if (!gcmp0((GEN)x[j])) x[i++] = x[j];
! 1243: setlg(x,i);
! 1244: }
! 1245: tetpil=avma;
! 1246: x = denx? gdiv(x,denx): gcopy(x);
! 1247: return gerepile(av0,tetpil,x);
! 1248: }
! 1249:
! 1250: GEN
! 1251: hnf(GEN x) { return hnf0(x,1); }
! 1252:
! 1253: #define cmod(x,u,us2) \
! 1254: {GEN a=modii((GEN)x,u); if (cmpii(a,us2)>0) a=subii(a,u); x=(long)a;}
! 1255:
! 1256: /* dm = multiple of diag element (usually detint(x)) */
! 1257: /* flag: don't/do append dm*matid. */
! 1258: static GEN
! 1259: allhnfmod(GEN x,GEN dm,long flag)
! 1260: {
! 1261: long li,co,av0,av,tetpil,i,j,k,def,ldef,lim,ldm;
! 1262: GEN a,b,w,p1,p2,d,u,v,dms2;
! 1263:
! 1264: if (typ(dm)!=t_INT) err(typeer,"allhnfmod");
! 1265: if (!signe(dm)) return hnf(x);
! 1266: if (typ(x)!=t_MAT) err(typeer,"allhnfmod");
! 1267: if (DEBUGLEVEL>6) fprintferr("Enter hnfmod");
! 1268:
! 1269: co=lg(x); if (co==1) return cgetg(1,t_MAT);
! 1270: av0=avma; lim=stack_lim(av0,1);
! 1271: dms2=shifti(dm,-1); li=lg(x[1]);
! 1272: av=avma;
! 1273:
! 1274: if (flag)
! 1275: {
! 1276: p1 = cgetg(co,t_MAT); for (i=1; i<co; i++) p1[i]=x[i];
! 1277: if (li > co) err(talker,"nb lines > nb columns in hnfmod");
! 1278: x = p1;
! 1279: }
! 1280: else
! 1281: {
! 1282: /* concatenate dm * Id to x */
! 1283: x = concatsp(x, idmat_intern(li-1,dm,gzero));
! 1284: for (i=1; i<co; i++) x[i] = lmod((GEN)x[i], dm);
! 1285: co += li-1;
! 1286: }
! 1287: def=co-1; ldef=0;
! 1288: for (i=li-1; i>ldef; i--,def--)
! 1289: {
! 1290: if (DEBUGLEVEL>6) { fprintferr(" %ld",i); flusherr(); }
! 1291: for (j=def-1; j; j--)
! 1292: {
! 1293: while (j && !signe(gcoeff(x,i,j))) j--;
! 1294: if (!j) break;
! 1295: if (DEBUGLEVEL>8) { fprintferr(" %ld",j); flusherr(); }
! 1296: k = (j==1)? def: j-1;
! 1297: a = gcoeff(x,i,j);
! 1298: b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
! 1299: if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
! 1300: p1 = lincomb_integral(u,v, (GEN)x[j], (GEN)x[k]);
! 1301: p2 = lincomb_integral(a, negi(b), (GEN)x[k], (GEN)x[j]);
! 1302: x[k] = (long)p1;
! 1303: x[j] = (long)p2;
! 1304: for (k=1; k<=i; k++)
! 1305: {
! 1306: cmod(p1[k], dm, dms2);
! 1307: cmod(p2[k], dm, dms2);
! 1308: }
! 1309: if (low_stack(lim, stack_lim(av0,1)))
! 1310: {
! 1311: if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i);
! 1312: tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
! 1313: }
! 1314: }
! 1315: }
! 1316: w=cgetg(li,t_MAT); b=dm;
! 1317: for (i=li-1; i>=1; i--)
! 1318: {
! 1319: d = bezout(gcoeff(x,i,i+def),b,&u,&v);
! 1320: w[i] = lmod(gmul(u,(GEN)x[i+def]), b);
! 1321: if (!signe(gcoeff(w,i,i))) coeff(w,i,i)=(long)d;
! 1322: if (i > 1 && flag) b=divii(b,d);
! 1323: }
! 1324: ldm = lgefint(dm);
! 1325: for (i=li-2; i>=1; i--)
! 1326: {
! 1327: GEN diag = gcoeff(w,i,i);
! 1328: for (j=i+1; j<li; j++)
! 1329: {
! 1330: b = negi(gdivent(gcoeff(w,i,j), diag));
! 1331: p1 = lincomb_integral(gun,b, (GEN)w[j], (GEN)w[i]);
! 1332: w[j] = (long)p1;
! 1333: for (k=1; k<i; k++)
! 1334: {
! 1335: p2 = (GEN)p1[k];
! 1336: if (lgefint(p2) > ldm) p1[k] = lmodii(p2, dm);
! 1337: }
! 1338: if (low_stack(lim, stack_lim(av0,1)))
! 1339: {
! 1340: if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i);
! 1341: tetpil=avma; w=gerepile(av,tetpil,gcopy(w));
! 1342: diag = gcoeff(w,i,i);
! 1343: }
! 1344: }
! 1345: }
! 1346: if (DEBUGLEVEL>6) { fprintferr("\nEnd hnfmod\n"); flusherr(); }
! 1347: tetpil=avma; return gerepile(av0,tetpil,gcopy(w));
! 1348: }
! 1349: #undef cmod
! 1350:
! 1351: GEN
! 1352: hnfmod(GEN x, GEN detmat) { return allhnfmod(x,detmat,1); }
! 1353:
! 1354: GEN
! 1355: hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }
! 1356:
! 1357: /* Return [y,U,V] such that y=V.x.U, V permutation vector, U unimodular
! 1358: * matrix, and y in HNF form
! 1359: */
! 1360: GEN
! 1361: hnfhavas(GEN x)
! 1362: {
! 1363: long av0=avma, av,av1,tetpil,li,co,i,j,k,def,ldef,lim,imin,jmin,vpk;
! 1364: long jpro,com,vi;
! 1365: GEN p1,p2,z,u,denx,vperm,mat1,col2,lil2,s,pmin,apro,bpro,cpro;
! 1366:
! 1367: if (DEBUGLEVEL>6)
! 1368: { fprintferr("Entering hnfhavas: AVMA = %ld\n",avma); flusherr(); }
! 1369: if (typ(x)!=t_MAT) err(typeer,"hnfhavas");
! 1370: co=lg(x);
! 1371: if (co==1)
! 1372: {
! 1373: z=cgetg(4,t_VEC); z[1]=lgetg(1,t_MAT);
! 1374: z[2]=lgetg(1,t_MAT); z[3]=lgetg(1,t_VEC);
! 1375: return z;
! 1376: }
! 1377: li=lg(x[1]); denx=denom(x);
! 1378: vperm=new_chunk(li); for (i=1; i<li; i++) vperm[i]=i;
! 1379:
! 1380: av=avma; lim=stack_lim(av,1); u=idmat(co-1);
! 1381: x = gcmp1(denx)? dummycopy(x): gmul(denx,x);
! 1382: def=co; ldef=(li>co)?li-co+1:1;
! 1383: for (i=li-1; i>=ldef; i--)
! 1384: {
! 1385: def--; av1=avma; mat1=cgetg(def+1,t_MAT); col2=cgetg(def+1,t_COL);
! 1386: for (j=1; j<=def; j++)
! 1387: {
! 1388: p1=cgetg(i+1,t_COL); mat1[j]=(long)p1; s=gzero;
! 1389: for (k=1; k<=i; k++)
! 1390: {
! 1391: p2=gsqr(gcoeff(x,vperm[k],j));
! 1392: p1[k]=(long)p2; s=gadd(s,p2);
! 1393: }
! 1394: col2[j]=(long)s;
! 1395: }
! 1396: lil2=cgetg(i+1,t_COL);
! 1397: for (k=1; k<=i; k++)
! 1398: {
! 1399: s=gzero;
! 1400: for (j=1; j<=def; j++) s=gadd(s,gcoeff(mat1,k,j));
! 1401: lil2[k]=(long)s;
! 1402: }
! 1403:
! 1404: pmin = NULL;
! 1405: for (k=i; k>=1; k--)
! 1406: {
! 1407: while (k>=1 && !signe(lil2[k])) k--;
! 1408: if (!k) goto comterm;
! 1409: vpk=vperm[k];
! 1410: if (!pmin || cmpii((GEN)lil2[k],pmin) <= 0)
! 1411: {
! 1412: j=1; while (!signe(gcoeff(x,vpk,j))) j++;
! 1413: if (!pmin)
! 1414: {
! 1415: imin=k; jmin=j; pmin=mulii((GEN)lil2[k],(GEN)col2[j]);
! 1416: cpro=absi(gcoeff(x,vpk,j));
! 1417: }
! 1418: jpro=j; apro=absi(gcoeff(x,vpk,j)); j++;
! 1419: for (; j<=def; j++)
! 1420: {
! 1421: com=cmpii((GEN)col2[j],(GEN)col2[jpro]);
! 1422: if (signe(gcoeff(x,vpk,j)) && com <=0)
! 1423: {
! 1424: if (com<0) { jpro=j; apro=absi(gcoeff(x,vpk,j)); }
! 1425: else
! 1426: {
! 1427: bpro=absi(gcoeff(x,vpk,j));
! 1428: if (cmpii(bpro,apro)<0) { jpro=j; apro=bpro; }
! 1429: }
! 1430: }
! 1431: }
! 1432: p1=mulii((GEN)lil2[k],(GEN)col2[jpro]);
! 1433: com=cmpii(p1,pmin);
! 1434: if (com<0 || (com==0 && cmpii(apro,cpro)<0))
! 1435: {
! 1436: imin=k; jmin=jpro; pmin=p1; cpro=apro;
! 1437: }
! 1438: }
! 1439: }
! 1440: avma=av1;
! 1441: if (jmin!=def)
! 1442: {
! 1443: p1=(GEN)x[def]; x[def]=x[jmin]; x[jmin]=(long)p1;
! 1444: p1=(GEN)u[def]; u[def]=u[jmin]; u[jmin]=(long)p1;
! 1445: }
! 1446: if (imin!=i) { vpk=vperm[i]; vperm[i]=vperm[imin]; vperm[imin]=vpk; }
! 1447: vi=vperm[i];
! 1448: for(;;)
! 1449: {
! 1450: GEN p3,p12,p13;
! 1451:
! 1452: if (signe(gcoeff(x,vi,def))<0)
! 1453: {
! 1454: x[def]=lneg((GEN)x[def]); u[def]=lneg((GEN)u[def]);
! 1455: }
! 1456: p1=gcoeff(x,vi,def); p12=shifti(p1,-1); p13=negi(p12);
! 1457: for (j=1; j<def; j++)
! 1458: {
! 1459: p2=dvmdii(gcoeff(x,vi,j),p1,&p3);
! 1460: if (cmpii(p3,p13)<0) p2=addis(p2,-1);
! 1461: else { if (cmpii(p3,p12)>0) p2=addis(p2,1); }
! 1462: if (DEBUGLEVEL>5) outerr(p2);
! 1463: setsigne(p2,-signe(p2));
! 1464: x[j]=ladd((GEN)x[j],gmul(p2,(GEN)x[def]));
! 1465: u[j]=ladd((GEN)u[j],gmul(p2,(GEN)u[def]));
! 1466: }
! 1467: j=1; while (!signe(gcoeff(x,vi,j))) j++;
! 1468: if (j<def)
! 1469: {
! 1470: pmin=gnorml2((GEN)x[j]); jmin=j; apro=absi(gcoeff(x,vi,j));
! 1471: j++;
! 1472: for (; j<def; j++)
! 1473: {
! 1474: if (signe(gcoeff(x,vi,j)))
! 1475: {
! 1476: p1=gnorml2((GEN)x[j]);
! 1477: com=cmpii(p1,pmin);
! 1478: if (com<0)
! 1479: {
! 1480: pmin=p1; jmin=j;
! 1481: }
! 1482: else if (com==0)
! 1483: {
! 1484: bpro=absi(gcoeff(x,vi,j));
! 1485: if (cmpii(bpro,apro)<0)
! 1486: {
! 1487: pmin=p1; jmin=j; apro=bpro;
! 1488: }
! 1489: }
! 1490: }
! 1491: }
! 1492: p1=(GEN)x[def]; x[def]=x[jmin]; x[jmin]=(long)p1;
! 1493: p1=(GEN)u[def]; u[def]=u[jmin]; u[jmin]=(long)p1;
! 1494: }
! 1495: else break;
! 1496: }
! 1497: vi=vperm[i]; p1=gcoeff(x,vi,def);
! 1498: for (j=def+1; j<co; j++)
! 1499: {
! 1500: p2=gdivent(gcoeff(x,vi,j),p1); setsigne(p2,-signe(p2));
! 1501: if (DEBUGLEVEL>5) outerr(p2);
! 1502: x[j]=ladd((GEN)x[j],gmul(p2,(GEN)x[def]));
! 1503: u[j]=ladd((GEN)u[j],gmul(p2,(GEN)u[def]));
! 1504: }
! 1505:
! 1506: if (low_stack(lim, stack_lim(av,1)))
! 1507: {
! 1508: GEN *gptr[2];
! 1509: if (DEBUGMEM>1) err(warnmem,"hnfhavas");
! 1510: gptr[0]=&x; gptr[1]=&u;
! 1511: gerepilemany(av,gptr,2);
! 1512: }
! 1513: }
! 1514:
! 1515: comterm:
! 1516: tetpil=avma; z=cgetg(4,t_VEC); p1=cgetg(co,t_MAT);
! 1517: if (gcmp1(denx))
! 1518: {
! 1519: for (j=1; j<co; j++)
! 1520: {
! 1521: p2=cgetg(li,t_COL); p1[j]=(long)p2;
! 1522: for (i=1; i<li; i++)
! 1523: p2[i] = lcopy(gcoeff(x,vperm[i],j));
! 1524: }
! 1525: }
! 1526: else
! 1527: {
! 1528: for (j=1; j<co; j++)
! 1529: {
! 1530: p2=cgetg(li,t_COL); p1[j]=(long)p2;
! 1531: for (i=1; i<li; i++)
! 1532: p2[i] = ldiv(gcoeff(x,vperm[i],j),denx);
! 1533: }
! 1534: }
! 1535: z[1]=(long)p1; z[2]=lcopy(u);
! 1536: p1=cgetg(li,t_VEC); for (i=1; i<li; i++) p1[i]=lstoi(vperm[i]);
! 1537: z[3]=(long)p1; return gerepile(av0,tetpil,z);
! 1538: }
! 1539:
! 1540: /* HNF by Bo Majewski and Allan Steele */
! 1541:
! 1542: /* premier indice non nul de la j-eme ligne de la matrice b */
! 1543: static long
! 1544: depthvector(GEN b,long j)
! 1545: {
! 1546: long lv = lg(b), i = 1;
! 1547:
! 1548: while (i<lv && gcmp0(gcoeff(b,j,i))) i++;
! 1549: return (i==lv)?-1:i;
! 1550: }
! 1551:
! 1552: static GEN
! 1553: incompleteprod(GEN b,long k,long l,long deb,long fin)
! 1554: {
! 1555: GEN p1 = gzero;
! 1556: long j;
! 1557:
! 1558: for (j=deb; j<=fin; j++)
! 1559: p1 = addii(p1,mulii(gcoeff(b,k,j),gcoeff(b,l,j)));
! 1560: return p1;
! 1561: }
! 1562:
! 1563: static void
! 1564: redlll(GEN b,GEN mu,long l,long c,long k)
! 1565: {
! 1566: long i,j, lb;
! 1567: GEN q, p1=gcoeff(b,l,c);
! 1568:
! 1569: if (signe(p1)) q=gdivround(gcoeff(b,k,c),p1); else q=ground(gcoeff(mu,k,l));
! 1570: if (signe(q))
! 1571: {
! 1572: q=negi(q); lb=lg(b);
! 1573: for (j=1; j<lb; j++)
! 1574: coeff(b,k,j) = laddii(gcoeff(b,k,j),mulii(q,gcoeff(b,l,j)));
! 1575: coeff(mu,k,l)=ladd(gcoeff(mu,k,l),q);
! 1576: for (i=1; i<l; i++)
! 1577: {
! 1578: p1=gcoeff(mu,l,i);
! 1579: if (gsigne(p1))
! 1580: coeff(mu,k,i) = ladd(gcoeff(mu,k,i),gmul(q,p1));
! 1581: }
! 1582: }
! 1583: }
! 1584:
! 1585: GEN
! 1586: hnflll(GEN x)
! 1587: {
! 1588: long n,m,i,j,k,ii,jj,p,c,s,av=avma,tetpil,kmax,ok;
! 1589: GEN q,qneg,p1,bnew,B,mu,mmu,cst,E,U,y,t,BB;
! 1590:
! 1591: if (typ(x)!=t_MAT) err(typeer,"hnflll");
! 1592: n=lg(x)-1;
! 1593: if (!n)
! 1594: {
! 1595: y=cgetg(3,t_VEC);
! 1596: y[1]=lgetg(1,t_MAT);
! 1597: y[2]=lgetg(1,t_MAT); return y;
! 1598: }
! 1599: cst=gdivgs(stoi(9),10);
! 1600: x=gcopy(x); n=lg(x)-1; m=lg(x[1])-1; p=n+m;
! 1601: bnew=cgetg(p+1,t_MAT);
! 1602: for (j=1; j<=n; j++) bnew[j]=x[j];
! 1603: for (j=n+1; j<=p; j++)
! 1604: {
! 1605: p1=cgetg(m+1,t_COL); bnew[j]=(long)p1;
! 1606: for (i=1; i<=m; i++) p1[i]=(i==(j-n))?un:zero;
! 1607: }
! 1608: x=bnew; c=n+1;
! 1609: for (i=1; i<=m; i++) c=min(c,depthvector(x,i));
! 1610: s=0; mu=cgetg(m+2,t_MAT);
! 1611: for (j=1; j<=m; j++) mu[j]=lgetg(m+2,t_COL); B=cgetg(m+2,t_COL);
! 1612: while (c<=n)
! 1613: {
! 1614: k=2; kmax=1;
! 1615: B[1]=(long)incompleteprod(x,1,1,c+1,p);
! 1616: while (k<=m-c+1)
! 1617: {
! 1618: if (k>kmax)
! 1619: {
! 1620: kmax=k;
! 1621: for (j=1; j<k; j++)
! 1622: {
! 1623: mmu=incompleteprod(x,k,j,c+1,p);
! 1624: for (i=1; i<j; i++) mmu=gsub(mmu,gmul(gcoeff(mu,j,i),gcoeff(mu,k,i)));
! 1625: coeff(mu,k,j)=(long)mmu;
! 1626: }
! 1627: for (j=1; j<k; j++) coeff(mu,k,j)=ldiv(gcoeff(mu,k,j),(GEN)B[j]);
! 1628: B[k]=(long)incompleteprod(x,k,k,c+1,p);
! 1629: for (j=1; j<k; j++)
! 1630: B[k]=lsub((GEN)B[k],gmul((GEN)B[j],gsqr(gcoeff(mu,k,j))));
! 1631: }
! 1632: redlll(x,mu,k-1,c,k);
! 1633: ok = (absi(gcoeff(x,k-1,c))>absi(gcoeff(x,k,c))) ||
! 1634: (gegal(gcoeff(x,k-1,c),gcoeff(x,k,c)) &&
! 1635: (gcmp((GEN) B[k],
! 1636: gmul(gsub(cst,gsqr(gcoeff(mu,k,k-1))), (GEN) B[k-1])) < 0) );
! 1637: while (ok)
! 1638: {
! 1639: for (j=1; j<=p; j++)
! 1640: {
! 1641: t=gcoeff(x,k,j); coeff(x,k,j)=coeff(x,k-1,j);
! 1642: coeff(x,k-1,j)=(long)t;
! 1643: }
! 1644: for (j=1; j<=k-2; j++)
! 1645: {
! 1646: t=gcoeff(mu,k,j); coeff(mu,k,j)=coeff(mu,k-1,j);
! 1647: coeff(mu,k-1,j)=(long)t;
! 1648: }
! 1649: mmu=gcoeff(mu,k,k-1);
! 1650: BB=gadd((GEN)B[k],gmul(gmul(mmu,mmu),(GEN)B[k-1]));
! 1651: q=gdiv((GEN)B[k-1],BB);
! 1652: coeff(mu,k,k-1)=lmul(mmu,q);
! 1653: B[k]=lmul((GEN)B[k],q); B[k-1]=(long)BB;
! 1654: for (i=k+1; i<=kmax; i++)
! 1655: {
! 1656: t=gcoeff(mu,i,k);
! 1657: coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mmu,t));
! 1658: coeff(mu,i,k-1)=ladd(t,gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
! 1659: }
! 1660: if (k>2) k--;
! 1661: redlll(x,mu,k-1,c,k);
! 1662: ok=(absi(gcoeff(x,k-1,c))>absi(gcoeff(x,k,c))) ||
! 1663: (gegal(gcoeff(x,k-1,c),gcoeff(x,k,c)) &&
! 1664: (gcmp((GEN) B[k],
! 1665: gmul(gsub(cst,gsqr(gcoeff(mu,k,k-1))), (GEN) B[k-1])) < 0));
! 1666: }
! 1667: for (i=k-2; i; i--) redlll(x,mu,i,c,k);
! 1668: k++;
! 1669: }
! 1670: s++; c=n+1;
! 1671: for (i=1; i<=m-s; i++) c=min(c,depthvector(x,i));
! 1672: }
! 1673: s=m-s+1;
! 1674: if (signe(gcoeff(x,s,depthvector(x,s)))<0)
! 1675: for (j=1; j<=p; j++)
! 1676: coeff(x,s,j)=lnegi(gcoeff(x,s,j));
! 1677: for (i=s+1; i<=m; i++)
! 1678: {
! 1679: if (signe(gcoeff(x,i,depthvector(x,i)))<0)
! 1680: for (j=1; j<=p; j++)
! 1681: coeff(x,i,j)=lnegi(gcoeff(x,i,j));
! 1682: for (j=i-1; j>=s; j--)
! 1683: {
! 1684: k=depthvector(x,j);
! 1685: qneg=negi(gdivent(gcoeff(x,i,k),gcoeff(x,j,k)));
! 1686: for (jj=1; jj<=p; jj++)
! 1687: coeff(x,i,jj)=laddii(gcoeff(x,i,jj),mulii(qneg,gcoeff(x,j,jj)));
! 1688: }
! 1689: }
! 1690: for (k=s; k<=m; k++)
! 1691: {
! 1692: for (j=1; j<s; j++)
! 1693: {
! 1694: mmu=incompleteprod(x,k,j,n+1,p);
! 1695: for (i=1; i<j; i++) mmu=gsub(mmu,gmul(gcoeff(mu,j,i),gcoeff(mu,k,i)));
! 1696: coeff(mu,k,j)=(long)mmu;
! 1697: }
! 1698: for (j=1; j<s; j++) coeff(mu,k,j)=ldiv(gcoeff(mu,k,j),(GEN)B[j]);
! 1699: B[k]=(long)incompleteprod(x,k,k,n+1,p);
! 1700: for (j=1; j<s; j++)
! 1701: B[k]=lsub((GEN)B[k],gmul(gsqr(gcoeff(mu,k,j)),(GEN)B[j]));
! 1702: for (j=s-1; j; j--)
! 1703: {
! 1704: qneg=negi(ground(gcoeff(mu,k,j)));
! 1705: if (signe(qneg))
! 1706: {
! 1707: for (jj=1; jj<=p; jj++)
! 1708: coeff(x,k,jj)=laddii(gcoeff(x,k,jj),mulii(qneg,gcoeff(x,j,jj)));
! 1709: for (i=1; i<j; i++)
! 1710: if (gsigne(gcoeff(mu,j,i)))
! 1711: coeff(mu,k,i)=ladd(gcoeff(mu,k,i),gmul(qneg,gcoeff(mu,j,i)));
! 1712: }
! 1713: }
! 1714: }
! 1715: tetpil=avma; y=cgetg(3,t_VEC);
! 1716: E=cgetg(n+1,t_MAT);
! 1717: for (i=1; i<=n; i++)
! 1718: {
! 1719: p1=cgetg(m-s+2,t_COL); E[i]=(long)p1;
! 1720: for (ii=1; ii<=m-s+1; ii++)
! 1721: p1[ii]=lcopy(gcoeff(x,m-ii+1,i));
! 1722: }
! 1723: y[1]=(long)E; U=cgetg(m+1,t_MAT);
! 1724: for (i=1; i<=m; i++)
! 1725: {
! 1726: p1=cgetg(m+1,t_COL); U[i]=(long)p1;
! 1727: for (ii=m; ii>=1; ii--)
! 1728: p1[m-ii+1]=lcopy(gcoeff(x,ii,i+n));
! 1729: }
! 1730: y[2]=(long)U; return gerepile(av,tetpil,y);
! 1731: }
! 1732:
! 1733: /* HNF with permutation */
! 1734: GEN
! 1735: hnfperm(GEN A)
! 1736: {
! 1737: GEN U,c,l,perm,d,u,v,p,q,y,a,b,p1;
! 1738: long r,t,i,j,j1,k,m,n,av=avma,av1,tetpil,lim;
! 1739:
! 1740: if (typ(A)!=t_MAT) err(typeer,"hnfperm");
! 1741: n=lg(A)-1;
! 1742: if (!n)
! 1743: {
! 1744: y=cgetg(4,t_VEC);
! 1745: y[1]=lgetg(1,t_MAT);
! 1746: y[2]=lgetg(1,t_MAT);
! 1747: y[3]=lgetg(1,t_VEC); return y;
! 1748: }
! 1749: m=lg(A[1])-1;
! 1750: c=new_chunk(m+1); for (i=1; i<=m; i++) c[i]=0;
! 1751: l=new_chunk(n+1); for (j=1; j<=n; j++) l[j]=0;
! 1752: perm=new_chunk(m+1);
! 1753: av1=avma; lim=stack_lim(av1,1);
! 1754: U=idmat(n); A=dummycopy(A);
! 1755: /* U base change matrix : A0*U=A all along */
! 1756:
! 1757: for (r=0,k=1; k<=n; k++)
! 1758: {
! 1759: for (j=1; j<k; j++) if (l[j])
! 1760: {
! 1761: t=l[j]; b=gcoeff(A,t,k);
! 1762: if (signe(b) == 0) continue;
! 1763:
! 1764: a=gcoeff(A,t,j); d=bezout(a,b,&u,&v);
! 1765: if (!is_pm1(d)) { a=divii(a,d); b=divii(b,d); }
! 1766: p1 = (GEN)A[j]; b = negi(b);
! 1767: A[j] = (long)lincomb_integral(u,v, p1,(GEN)A[k]);
! 1768: A[k] = (long)lincomb_integral(a,b, (GEN)A[k],p1);
! 1769: p1 = (GEN)U[j];
! 1770: U[j] = (long)lincomb_integral(u,v, p1,(GEN)U[k]);
! 1771: U[k] = (long)lincomb_integral(a,b, (GEN)U[k],p1);
! 1772: for (j1=1; j1<j; j1++) if (l[j1])
! 1773: {
! 1774: q=truedvmdii(gcoeff(A,t,j1),d,NULL);
! 1775: if (signe(q))
! 1776: {
! 1777: q = negi(q);
! 1778: A[j1] = (long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[j]);
! 1779: U[j1] = (long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[j]);
! 1780: }
! 1781: }
! 1782: }
! 1783: t=m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
! 1784: if (t)
! 1785: {
! 1786: p = gcoeff(A,t,k);
! 1787: for (i=t-1; i; i--)
! 1788: {
! 1789: q = gcoeff(A,i,k);
! 1790: if (signe(q) && absi_cmp(p,q) > 0) { p=q; t=i; }
! 1791: }
! 1792: perm[++r]=l[k]=t; c[t]=k;
! 1793: if (signe(p) < 0)
! 1794: {
! 1795: for (i=1; i<=m; i++) coeff(A,i,k)= lnegi(gcoeff(A,i,k));
! 1796: for (i=1; i<=n; i++) coeff(U,i,k)= lnegi(gcoeff(U,i,k));
! 1797: p=gcoeff(A,t,k);
! 1798: }
! 1799: for (j=1; j<k; j++) if (l[j])
! 1800: {
! 1801: q=truedvmdii(gcoeff(A,t,j),p,NULL);
! 1802: if (signe(q))
! 1803: {
! 1804: q = negi(q);
! 1805: A[j]=(long)lincomb_integral(gun,q,(GEN)A[j],(GEN)A[k]);
! 1806: U[j]=(long)lincomb_integral(gun,q,(GEN)U[j],(GEN)U[k]);
! 1807: }
! 1808: }
! 1809: }
! 1810: if (low_stack(lim, stack_lim(av1,1)))
! 1811: {
! 1812: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
! 1813: if (DEBUGMEM>1) err(warnmem,"hnfperm");
! 1814: gerepilemany(av1,gptr,2);
! 1815: }
! 1816: }
! 1817: if (r < m)
! 1818: {
! 1819: for (i=1,k=r; i<=m; i++)
! 1820: if (c[i]==0) perm[++k] = i;
! 1821: }
! 1822:
! 1823: /* We have A0*U=A, U in Gl(n,Z)
! 1824: * basis for Im(A): columns of A s.t l[j]>0 (r cols)
! 1825: * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols)
! 1826: */
! 1827: tetpil=avma; y=cgetg(4,t_VEC);
! 1828: p=cgetg(r+1,t_MAT); u=cgetg(n+1,t_MAT);
! 1829: for (t=1,k=r,j=1; j<=n; j++)
! 1830: if (l[j])
! 1831: {
! 1832: q=cgetg(m+1,t_COL); p[k]=(long)q;
! 1833: for (i=1; i<=m; i++) q[i]=lcopy(gcoeff(A,perm[m-i+1],j));
! 1834: u[k+n-r]=lcopy((GEN)U[j]);
! 1835: k--;
! 1836: }
! 1837: else u[t++]=lcopy((GEN)U[j]);
! 1838: y[1]=(long)p; y[2]=(long)u;
! 1839: q = cgetg(m+1,t_VEC); y[3]=(long)q;
! 1840: for (i=1; i<=m; i++) q[m-i+1]=lstoi(perm[i]);
! 1841: return gerepile(av,tetpil,y);
! 1842: }
! 1843:
! 1844: GEN
! 1845: colint_copy(GEN x)
! 1846: {
! 1847: long i, lx = lg(x);
! 1848: GEN y = cgetg(lx, t_COL);
! 1849: for (i=1; i<lx; i++) y[i] = signe(x[i])? licopy((GEN)x[i]): zero;
! 1850: return y;
! 1851: }
! 1852:
! 1853: GEN
! 1854: matint_copy(GEN x)
! 1855: {
! 1856: long i, lx = lg(x);
! 1857: GEN y = cgetg(lx, t_MAT);
! 1858:
! 1859: for (i=1; i<lx; i++)
! 1860: y[i] = (long)colint_copy((GEN)x[i]);
! 1861: return y;
! 1862: }
! 1863:
! 1864: /*====================================================================
! 1865: * Forme Normale d'Hermite (Version par colonnes 31/01/94)
! 1866: *====================================================================*/
! 1867: GEN
! 1868: hnfall0(GEN A, long remove)
! 1869: {
! 1870: GEN U,c,h,x,y,u,v,p,q,d,a,b,p1;
! 1871: long m,n,r,i,j,j1,k,li,z,av=avma,av1,tetpil,lim;
! 1872:
! 1873: if (typ(A)!=t_MAT) err(typeer,"hnfall");
! 1874: n=lg(A)-1;
! 1875: if (!n)
! 1876: {
! 1877: y=cgetg(3,t_VEC);
! 1878: y[1]=lgetg(1,t_MAT);
! 1879: y[2]=lgetg(1,t_MAT); return y;
! 1880: }
! 1881: m=lg(A[1])-1;
! 1882: c=new_chunk(m+1); for (i=1; i<=m; i++) c[i]=0;
! 1883: h=new_chunk(n+1); for (j=1; j<=n; j++) h[j]=m;
! 1884: av1=avma; lim=stack_lim(av1,1);
! 1885: A=dummycopy(A); U=idmat(n); r=n+1;
! 1886: for (li=m; li; li--)
! 1887: {
! 1888: for (j=1; j<r; j++)
! 1889: {
! 1890: for (i=h[j]; i>li; i--)
! 1891: {
! 1892: b = gcoeff(A,i,j);
! 1893: if (signe(b))
! 1894: {
! 1895: k=c[i]; a=gcoeff(A,i,k); /* annuler bij A l'aide de p=bik */
! 1896: d=bezout(a,b,&u,&v);
! 1897: if (DEBUGLEVEL>5)
! 1898: { fprintferr("(u,v) = (%Z, %Z); ",u,v); flusherr(); }
! 1899: if (!is_pm1(d)) { a=divii(a,d); b =divii(b,d); }
! 1900: p1 = (GEN)A[k]; b = negi(b);
! 1901: A[k] = (long)lincomb_integral(u,v, p1,(GEN)A[j]);
! 1902: A[j] = (long)lincomb_integral(a,b, (GEN)A[j],p1);
! 1903: p1 = (GEN)U[k];
! 1904: U[k] = (long)lincomb_integral(u,v, p1,(GEN)U[j]);
! 1905: U[j] = (long)lincomb_integral(a,b, (GEN)U[j],p1);
! 1906: for (j1=k+1; j1<=n; j1++)
! 1907: {
! 1908: q=truedvmdii(gcoeff(A,i,j1),d,NULL);
! 1909: if (signe(q))
! 1910: {
! 1911: q = negi(q);
! 1912: A[j1]=(long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[k]);
! 1913: U[j1]=(long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[k]);
! 1914: }
! 1915: }
! 1916: if (low_stack(lim, stack_lim(av1,1)))
! 1917: {
! 1918: GEN *gptr[2]; tetpil = avma;
! 1919: A = matint_copy(A); gptr[0]=&A;
! 1920: U = matint_copy(U); gptr[1]=&U;
! 1921: if (DEBUGMEM>1) err(warnmem,"hnfall[1], li = %ld", li);
! 1922: gerepilemanysp(av1,tetpil,gptr,2);
! 1923: }
! 1924: }
! 1925: }
! 1926: x=gcoeff(A,li,j);
! 1927: if (signe(x))
! 1928: {
! 1929: r--;
! 1930: if (j<r)
! 1931: {
! 1932: z=A[j]; A[j]=A[r]; A[r]=z;
! 1933: z=U[j]; U[j]=U[r]; U[r]=z;
! 1934: h[j]=h[r]; h[r]=li; c[li]=r;
! 1935: }
! 1936: if (signe(gcoeff(A,li,r))<0)
! 1937: {
! 1938: p1=(GEN)A[r]; for (i=1; i<=li; i++) p1[i]=lnegi((GEN)p1[i]);
! 1939: p1=(GEN)U[r]; for (i=1; i<=n ; i++) p1[i]=lnegi((GEN)p1[i]);
! 1940: }
! 1941: p=gcoeff(A,li,r);
! 1942: for (j=r+1; j<=n; j++)
! 1943: {
! 1944: q = truedvmdii(gcoeff(A,li,j),p,NULL);
! 1945: if (signe(q))
! 1946: {
! 1947: q = negi(q);
! 1948: A[j]=(long)lincomb_integral(gun,q,(GEN)A[j],(GEN)A[r]);
! 1949: U[j]=(long)lincomb_integral(gun,q,(GEN)U[j],(GEN)U[r]);
! 1950: }
! 1951: }
! 1952: if (low_stack(lim, stack_lim(av1,1)))
! 1953: {
! 1954: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
! 1955: if (DEBUGMEM>1) err(warnmem,"hnfall[2], li = %ld", li);
! 1956: gerepilemany(av1,gptr,2);
! 1957: }
! 1958: break;
! 1959: }
! 1960: h[j]=li-1;
! 1961: }
! 1962: }
! 1963: if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");
! 1964: r--; /* first r cols are in the image the n-r (independent) last ones */
! 1965: for (j=1; j<=r; j++)
! 1966: for (i=h[j]; i; i--)
! 1967: if (signe(b=gcoeff(A,i,j)))
! 1968: {
! 1969: k=c[i]; a=gcoeff(A,i,k);
! 1970: d=bezout(a,b,&u,&v);
! 1971: if (!is_pm1(d)) { a=divii(a,d); b=divii(b,d); }
! 1972: if (DEBUGLEVEL>5)
! 1973: { fprintferr("(u,v) = (%Z, %Z); ",u,v); flusherr(); }
! 1974: p1 = (GEN)A[k]; b = negi(b);
! 1975: A[k] = (long)lincomb_integral(u,v, p1,(GEN)A[j]);
! 1976: A[j] = (long)lincomb_integral(a,b, (GEN)A[j],p1);
! 1977: p1 = (GEN)U[k];
! 1978: U[k] = (long)lincomb_integral(u,v, p1,(GEN)U[j]);
! 1979: U[j] = (long)lincomb_integral(a,b, (GEN)U[j],p1);
! 1980: for (j1=k+1; j1<=n; j1++)
! 1981: {
! 1982: q=truedvmdii(gcoeff(A,i,j1),d,NULL);
! 1983: if (signe(q))
! 1984: {
! 1985: q = negi(q);
! 1986: A[j1] = (long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[k]);
! 1987: U[j1] = (long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[k]);
! 1988: }
! 1989: }
! 1990: if (low_stack(lim, stack_lim(av1,1)))
! 1991: {
! 1992: GEN *gptr[2]; tetpil = avma;
! 1993: A = matint_copy(A); gptr[0]=&A;
! 1994: U = matint_copy(U); gptr[1]=&U;
! 1995: if (DEBUGMEM>1) err(warnmem,"hnfall[3], j = %ld", j);
! 1996: gerepilemanysp(av1,tetpil,gptr,2);
! 1997: }
! 1998: }
! 1999: if (DEBUGLEVEL>5) fprintferr("\n");
! 2000: tetpil=avma; y=cgetg(3,t_VEC);
! 2001: if (remove)
! 2002: { /* remove the first r columns */
! 2003: A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1);
! 2004: }
! 2005: y[1]=lcopy(A); y[2]=lcopy(U);
! 2006: return gerepile(av,tetpil,y);
! 2007: }
! 2008:
! 2009: GEN
! 2010: hnfall(GEN x) {return hnfall0(x,1);}
! 2011:
! 2012: /***************************************************************/
! 2013: /** **/
! 2014: /** SMITH NORMAL FORM REDUCTION **/
! 2015: /** **/
! 2016: /***************************************************************/
! 2017:
! 2018: static GEN
! 2019: col_mul(GEN x, GEN c)
! 2020: {
! 2021: long s = signe(x);
! 2022: GEN xc = NULL;
! 2023: if (s)
! 2024: {
! 2025: if (!is_pm1(x)) xc = gmul(x,c);
! 2026: else xc = (s>0)? c: gneg_i(c);
! 2027: }
! 2028: return xc;
! 2029: }
! 2030:
! 2031: static void
! 2032: do_zero(GEN x)
! 2033: {
! 2034: long i, lx = lg(x);
! 2035: for (i=1; i<lx; i++) x[i] = zero;
! 2036: }
! 2037:
! 2038: /* c1 <-- u.c1 + v.c2; c2 <-- a.c2 - b.c1 */
! 2039: static void
! 2040: update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
! 2041: {
! 2042: GEN p1,p2;
! 2043:
! 2044: u = col_mul(u,*c1);
! 2045: v = col_mul(v,*c2);
! 2046: if (u) p1 = v? gadd(u,v): u;
! 2047: else p1 = v? v: (GEN)NULL;
! 2048:
! 2049: a = col_mul(a,*c2);
! 2050: b = col_mul(gneg_i(b),*c1);
! 2051: if (a) p2 = b? gadd(a,b): a;
! 2052: else p2 = b? b: (GEN)NULL;
! 2053:
! 2054: if (!p1) do_zero(*c1); else *c1 = p1;
! 2055: if (!p2) do_zero(*c2); else *c2 = p2;
! 2056: }
! 2057:
! 2058: static GEN
! 2059: trivsmith(long all)
! 2060: {
! 2061: GEN z;
! 2062: if (!all) return cgetg(1,t_VEC);
! 2063: z=cgetg(4,t_VEC);
! 2064: z[1]=lgetg(1,t_MAT);
! 2065: z[2]=lgetg(1,t_MAT);
! 2066: z[3]=lgetg(1,t_VEC); return z;
! 2067: }
! 2068:
! 2069: /* Return the smith normal form d of matrix x. If all != 0 return [d,u,v],
! 2070: * where d = u.x.v
! 2071: */
! 2072: static GEN
! 2073: smithall(GEN x, long all)
! 2074: {
! 2075: long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim;
! 2076: GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys;
! 2077:
! 2078: if (typ(x)!=t_MAT) err(typeer,"smithall");
! 2079: if (DEBUGLEVEL>=9) outerr(x);
! 2080: n=lg(x)-1;
! 2081: if (!n) return trivsmith(all);
! 2082: if (lg(x[1]) != n+1) err(mattype1,"smithall");
! 2083: for (i=1; i<=n; i++)
! 2084: for (j=1; j<=n; j++)
! 2085: if (typ(coeff(x,i,j)) != t_INT)
! 2086: err(talker,"non integral matrix in smithall");
! 2087:
! 2088: av = avma; lim = stack_lim(av,1);
! 2089: x = dummycopy(x); mdet = detint(x);
! 2090: if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }
! 2091: else
! 2092: {
! 2093: if (signe(mdet))
! 2094: {
! 2095: p1=hnfmod(x,mdet);
! 2096: if (all) { ml=idmat(n); mr=gauss(x,p1); }
! 2097: }
! 2098: else
! 2099: {
! 2100: if (all)
! 2101: {
! 2102: p1 = hnfall0(x,0);
! 2103: ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1];
! 2104: }
! 2105: else
! 2106: p1 = hnf0(x,0);
! 2107: }
! 2108: x = p1;
! 2109: }
! 2110: p1=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) p1[i]=lnegi(gcoeff(x,i,i));
! 2111: p2=sindexsort(p1); ys=cgetg(n+1,t_MAT);
! 2112: for (j=1; j<=n; j++)
! 2113: {
! 2114: p1=cgetg(n+1,t_COL); ys[j]=(long)p1;
! 2115: for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]);
! 2116: }
! 2117: x = ys;
! 2118: if (all)
! 2119: {
! 2120: p3=cgetg(n+1,t_MAT); p4=cgetg(n+1,t_MAT);
! 2121: for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }
! 2122: ml=p3; mr=p4;
! 2123: }
! 2124: if (signe(mdet))
! 2125: {
! 2126: p1 = hnfmod(x,mdet);
! 2127: if (all) mr=gmul(mr,gauss(x,p1));
! 2128: }
! 2129: else
! 2130: {
! 2131: if (all)
! 2132: {
! 2133: p1 = hnfall0(x,0);
! 2134: mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1];
! 2135: }
! 2136: else
! 2137: p1 = hnf0(x,0);
! 2138: }
! 2139: x=p1; mun = negi(gun);
! 2140:
! 2141: if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();}
! 2142: for (i=n; i>=2; i--)
! 2143: {
! 2144: if (DEBUGLEVEL>7) {fprintferr("\ni = %ld: ",i);flusherr();}
! 2145: for(;;)
! 2146: {
! 2147: c = 0;
! 2148: for (j=i-1; j>=1; j--)
! 2149: {
! 2150: p1=gcoeff(x,i,j); s1 = signe(p1);
! 2151: if (s1)
! 2152: {
! 2153: p2=gcoeff(x,i,i);
! 2154: if (!absi_cmp(p1,p2))
! 2155: {
! 2156: s2=signe(p2);
! 2157: if (s1 == s2) { d=p1; u=gun; p4=gun; }
! 2158: else
! 2159: {
! 2160: if (s2>0) { u = gun; p4 = mun; }
! 2161: else { u = mun; p4 = gun; }
! 2162: d=(s1>0)? p1: absi(p1);
! 2163: }
! 2164: v = gzero; p3 = u;
! 2165: }
! 2166: else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
! 2167: for (k=1; k<=i; k++)
! 2168: {
! 2169: b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));
! 2170: coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),
! 2171: mulii(p4,gcoeff(x,k,i)));
! 2172: coeff(x,k,i)=(long)b;
! 2173: }
! 2174: if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
! 2175: if (low_stack(lim, stack_lim(av,1)))
! 2176: {
! 2177: if (DEBUGMEM>1) err(warnmem,"[1]: smithall");
! 2178: if (all)
! 2179: {
! 2180: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 2181: gerepilemany(av,gptr,3);
! 2182: }
! 2183: else x=gerepileupto(av,gcopy(x));
! 2184: }
! 2185: }
! 2186: }
! 2187: if (DEBUGLEVEL>=8) {fprintferr("; ");flusherr();}
! 2188: for (j=i-1; j>=1; j--)
! 2189: {
! 2190: p1=gcoeff(x,j,i); s1 = signe(p1);
! 2191: if (s1)
! 2192: {
! 2193: p2=gcoeff(x,i,i);
! 2194: if (!absi_cmp(p1,p2))
! 2195: {
! 2196: s2 = signe(p2);
! 2197: if (s1 == s2) { d=p1; u=gun; p4=gun; }
! 2198: else
! 2199: {
! 2200: if (s2>0) { u = gun; p4 = mun; }
! 2201: else { u = mun; p4 = gun; }
! 2202: d=(s1>0)? p1: absi(p1);
! 2203: }
! 2204: v = gzero; p3 = u;
! 2205: }
! 2206: else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
! 2207: for (k=1; k<=i; k++)
! 2208: {
! 2209: b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
! 2210: coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),
! 2211: mulii(p4,gcoeff(x,i,k)));
! 2212: coeff(x,i,k)=(long)b;
! 2213: }
! 2214: if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
! 2215: c = 1;
! 2216: }
! 2217: }
! 2218: if (!c)
! 2219: {
! 2220: b=gcoeff(x,i,i); fl=1;
! 2221: if (signe(b))
! 2222: {
! 2223: for (k=1; k<i && fl; k++)
! 2224: for (l=1; l<i && fl; l++)
! 2225: fl = (int)!signe(resii(gcoeff(x,k,l),b));
! 2226: /* cast to (int) necessary for gcc-2.95 on sparcv9-64 (IS) */
! 2227: if (!fl)
! 2228: {
! 2229: k--;
! 2230: for (l=1; l<=i; l++)
! 2231: coeff(x,i,l)=laddii(gcoeff(x,i,l),gcoeff(x,k,l));
! 2232: if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
! 2233: }
! 2234: }
! 2235: if (fl) break;
! 2236: }
! 2237: if (low_stack(lim, stack_lim(av,1)))
! 2238: {
! 2239: if (DEBUGMEM>1) err(warnmem,"[2]: smithall");
! 2240: if (all)
! 2241: {
! 2242: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 2243: gerepilemany(av,gptr,3);
! 2244: }
! 2245: else x=gerepileupto(av,gcopy(x));
! 2246: }
! 2247: }
! 2248: }
! 2249: if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();}
! 2250: if (all)
! 2251: {
! 2252: for (k=1; k<=n; k++)
! 2253: if (signe(gcoeff(x,k,k))<0)
! 2254: { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lnegi(gcoeff(x,k,k)); }
! 2255: tetpil=avma; z=cgetg(4,t_VEC);
! 2256: z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
! 2257: return gerepile(av,tetpil,z);
! 2258: }
! 2259: tetpil=avma; z=cgetg(n+1,t_VEC); j=n;
! 2260: for (k=n; k; k--)
! 2261: if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));
! 2262: for ( ; j; j--) z[j]=zero;
! 2263: return gerepile(av,tetpil,z);
! 2264: }
! 2265:
! 2266: GEN
! 2267: smith(GEN x) { return smithall(x,0); }
! 2268:
! 2269: GEN
! 2270: smith2(GEN x) { return smithall(x,1); }
! 2271:
! 2272: /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
! 2273: GEN
! 2274: smithclean(GEN z)
! 2275: {
! 2276: long i,j,l,c;
! 2277: GEN u,v,y,d,p1;
! 2278:
! 2279: if (typ(z) != t_VEC) err(typeer,"smithclean");
! 2280: l = lg(z); if (l == 1) return cgetg(1,t_VEC);
! 2281: u=(GEN)z[1];
! 2282: if (l != 4 || typ(u) != t_MAT)
! 2283: {
! 2284: if (typ(u) != t_INT) err(typeer,"smithclean");
! 2285: for (c=1; c<l; c++)
! 2286: if (gcmp1((GEN)z[c])) break;
! 2287: return gcopy_i(z, c);
! 2288: }
! 2289: v=(GEN)z[2]; d=(GEN)z[3]; l = lg(d);
! 2290: for (c=1; c<l; c++)
! 2291: if (gcmp1(gcoeff(d,c,c))) break;
! 2292:
! 2293: y=cgetg(4,t_VEC);
! 2294: y[1]=(long)(p1 = cgetg(l,t_MAT));
! 2295: for (i=1; i<l; i++) p1[i] = (long)gcopy_i((GEN)u[i], c);
! 2296: y[2]=(long)gcopy_i(v, c);
! 2297: y[3]=(long)(p1 = cgetg(c,t_MAT));
! 2298: for (i=1; i<c; i++)
! 2299: {
! 2300: GEN p2 = cgetg(c,t_COL); p1[i] = (long)p2;
! 2301: for (j=1; j<c; j++) p2[j] = i==j? lcopy(gcoeff(d,i,i)): zero;
! 2302: }
! 2303: return y;
! 2304: }
! 2305:
! 2306: static GEN
! 2307: gsmithall(GEN x,long all)
! 2308: {
! 2309: long av,tetpil,i,j,k,l,c,fl,n, lim;
! 2310: GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr;
! 2311:
! 2312: if (typ(x)!=t_MAT) err(typeer,"gsmithall");
! 2313: n=lg(x)-1;
! 2314: if (!n) return trivsmith(all);
! 2315: if (lg(x[1]) != n+1) err(mattype1,"gsmithall");
! 2316: av = avma; lim = stack_lim(av,1);
! 2317: x = dummycopy(x);
! 2318: if (all) { ml=idmat(n); mr=idmat(n); }
! 2319: for (i=n; i>=2; i--)
! 2320: {
! 2321: do
! 2322: {
! 2323: c=0;
! 2324: for (j=i-1; j>=1; j--)
! 2325: {
! 2326: p1=gcoeff(x,i,j);
! 2327: if (signe(p1))
! 2328: {
! 2329: p2=gcoeff(x,i,i); v=gdiventres(p1,p2);
! 2330: if (gcmp0((GEN)v[2])) { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
! 2331: else { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
! 2332: for (k=1; k<=i; k++)
! 2333: {
! 2334: b=gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
! 2335: coeff(x,k,j)=lsub(gmul(p3,gcoeff(x,k,j)),gmul(p4,gcoeff(x,k,i)));
! 2336: coeff(x,k,i)=(long)b;
! 2337: }
! 2338: if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
! 2339: }
! 2340: }
! 2341: for (j=i-1; j>=1; j--)
! 2342: {
! 2343: p1=gcoeff(x,j,i);
! 2344: if (signe(p1))
! 2345: {
! 2346: p2=gcoeff(x,i,i); v=gdiventres(p1,p2);
! 2347: if (gcmp0((GEN)v[2])) { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
! 2348: else { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
! 2349: for (k=1; k<=i; k++)
! 2350: {
! 2351: b=gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
! 2352: coeff(x,j,k)=lsub(gmul(p3,gcoeff(x,j,k)),gmul(p4,gcoeff(x,i,k)));
! 2353: coeff(x,i,k)=(long)b;
! 2354: }
! 2355: if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
! 2356: c = 1;
! 2357: }
! 2358: }
! 2359: if (!c)
! 2360: {
! 2361: b=gcoeff(x,i,i); fl=1;
! 2362: if (signe(b))
! 2363: {
! 2364: for (k=1; (k<i)&&fl; k++)
! 2365: for (l=1; (l<i)&&fl; l++)
! 2366: fl= !signe(gmod(gcoeff(x,k,l),b));
! 2367: if (!fl)
! 2368: {
! 2369: k--;
! 2370: for (l=1; l<=i; l++)
! 2371: coeff(x,i,l)=ladd(gcoeff(x,i,l),gcoeff(x,k,l));
! 2372: if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
! 2373: }
! 2374: }
! 2375: }
! 2376: if (low_stack(lim, stack_lim(av,1)))
! 2377: {
! 2378: if (DEBUGMEM>1) err(warnmem,"[5]: smithall");
! 2379: tetpil=avma;
! 2380: if (all)
! 2381: {
! 2382: GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
! 2383: gerepilemany(av,gptr,3);
! 2384: }
! 2385: else x=gerepile(av,tetpil,gcopy(x));
! 2386: }
! 2387: }
! 2388: while (c || !fl);
! 2389: }
! 2390: if (all)
! 2391: {
! 2392: for (k=1; k<=n; k++)
! 2393: if (signe(gcoeff(x,k,k))<0)
! 2394: { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lneg(gcoeff(x,k,k)); }
! 2395: tetpil=avma; z=cgetg(4,t_VEC);
! 2396: z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
! 2397: }
! 2398: else
! 2399: {
! 2400: tetpil=avma; z=cgetg(n+1,t_VEC);
! 2401: for (j=0,k=1; k<=n; k++) if (!signe(gcoeff(x,k,k))) z[++j]=zero;
! 2402: for (k=1; k<=n; k++)
! 2403: if (signe(p1=gcoeff(x,k,k))) z[++j]=(long)gabs(p1,0);
! 2404: }
! 2405: return gerepile(av,tetpil,z);
! 2406: }
! 2407:
! 2408: GEN
! 2409: matsnf0(GEN x,long flag)
! 2410: {
! 2411: long av = avma;
! 2412: if (flag > 7) err(flagerr,"matsnf");
! 2413: if (typ(x) == t_VEC && flag & 4) return smithclean(x);
! 2414: if (flag & 2) x = gsmithall(x,flag & 1);
! 2415: else x = smithall(x, flag & 1);
! 2416: if (flag & 4) x = smithclean(x);
! 2417: return gerepileupto(av, x);
! 2418: }
! 2419:
! 2420: GEN
! 2421: gsmith(GEN x) { return gsmithall(x,0); }
! 2422:
! 2423: GEN
! 2424: 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 / src / basemath