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