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