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