Annotation of OpenXM_contrib/pari/src/basemath/alglin1.c, Revision 1.1.1.1
1.1 maekawa 1: /********************************************************************/
2: /** **/
3: /** LINEAR ALGEBRA **/
4: /** (first part) **/
5: /** **/
6: /********************************************************************/
7: /* $Id: alglin1.c,v 1.1.1.1 1999/09/16 13:47:15 karim Exp $ */
8: #include "pari.h"
9:
10: /*******************************************************************/
11: /* */
12: /* TRANSPOSE */
13: /* */
14: /*******************************************************************/
15:
16: GEN
17: gtrans(GEN x)
18: {
19: long i,j,lx,dx, tx=typ(x);
20: GEN y,p1;
21:
22: if (! is_matvec_t(tx)) err(typeer,"gtrans");
23: switch(tx)
24: {
25: case t_VEC:
26: y=gcopy(x); settyp(y,t_COL); break;
27:
28: case t_COL:
29: y=gcopy(x); settyp(y,t_VEC); break;
30:
31: case t_MAT:
32: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
33: dx=lg(x[1]); y=cgetg(dx,tx);
34: for (i=1; i<dx; i++)
35: {
36: p1=cgetg(lx,t_COL); y[i]=(long)p1;
37: for (j=1; j<lx; j++) p1[j]=lcopy(gcoeff(x,i,j));
38: }
39: break;
40:
41: default: y=gcopy(x); break;
42: }
43: return y;
44: }
45:
46: /*******************************************************************/
47: /* */
48: /* CONCATENATION & EXTRACTION */
49: /* */
50: /*******************************************************************/
51:
52: static GEN
53: strconcat(GEN x, GEN y)
54: {
55: long flx=0,fly=0,l;
56: char *sx,*sy,*str;
57:
58: if (typ(x)==t_STR) sx = GSTR(x); else { flx=1; sx = GENtostr(x); }
59: if (typ(y)==t_STR) sy = GSTR(y); else { fly=1; sy = GENtostr(y); }
60: l = strlen(sx) + strlen(sy) + 1;
61: l = (l+BYTES_IN_LONG) >> TWOPOTBYTES_IN_LONG;
62: x = cgetg(l+1, t_STR); str = GSTR(x);
63: strcpy(str,sx);
64: strcat(str,sy);
65: if (flx) free(sx);
66: if (fly) free(sy);
67: return x;
68: }
69:
70: GEN
71: concatsp(GEN x, GEN y)
72: {
73: long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i;
74: GEN z,p1;
75:
76: if (tx==t_LIST || ty==t_LIST) return listconcat(x,y);
77: if (tx==t_STR || ty==t_STR) return strconcat(x,y);
78:
79: if (tx==t_MAT && lx==1)
80: {
81: if (ty!=t_VEC || ly==1) return gtomat(y);
82: err(concater);
83: }
84: if (ty==t_MAT && ly==1)
85: {
86: if (tx!=t_VEC || lx==1) return gtomat(x);
87: err(concater);
88: }
89:
90: if (! is_matvec_t(tx))
91: {
92: if (! is_matvec_t(ty))
93: {
94: z=cgetg(3,t_VEC); z[1]=(long)x; z[2]=(long)y;
95: return z;
96: }
97: z=cgetg(ly+1,ty);
98: if (ty != t_MAT) p1 = x;
99: else
100: {
101: if (lg(y[1])!=2) err(concater);
102: p1=cgetg(2,t_COL); p1[1]=(long)x;
103: }
104: for (i=2; i<=ly; i++) z[i]=y[i-1];
105: z[1]=(long)p1; return z;
106: }
107: if (! is_matvec_t(ty))
108: {
109: z=cgetg(lx+1,tx);
110: if (tx != t_MAT) p1 = y;
111: else
112: {
113: if (lg(x[1])!=2) err(concater);
114: p1=cgetg(2,t_COL); p1[1]=(long)y;
115: }
116: for (i=1; i<lx; i++) z[i]=x[i];
117: z[lx]=(long)p1; return z;
118: }
119:
120: if (tx == ty)
121: {
122: if (tx == t_MAT && lg(x[1]) != lg(y[1])) err(concater);
123: z=cgetg(lx+ly-1,tx);
124: for (i=1; i<lx; i++) z[i]=x[i];
125: for (i=1; i<ly; i++) z[lx+i-1]=y[i];
126: return z;
127: }
128:
129: switch(tx)
130: {
131: case t_VEC:
132: switch(ty)
133: {
134: case t_COL:
135: if (lx<=2) return (lx==1)? y: concatsp((GEN) x[1],y);
136: if (ly>=3) break;
137: return (ly==1)? x: concatsp(x,(GEN) y[1]);
138: case t_MAT:
139: z=cgetg(ly,ty); if (lx != ly) break;
140: for (i=1; i<ly; i++) z[i]=(long)concatsp((GEN) x[i],(GEN) y[i]);
141: return z;
142: }
143: break;
144:
145: case t_COL:
146: switch(ty)
147: {
148: case t_VEC:
149: if (lx<=2) return (lx==1)? y: concatsp((GEN) x[1],y);
150: if (ly>=3) break;
151: return (ly==1)? x: concatsp(x,(GEN) y[1]);
152: case t_MAT:
153: if (lx != lg(y[1])) break;
154: z=cgetg(ly+1,ty); z[1]=(long)x;
155: for (i=2; i<=ly; i++) z[i]=y[i-1];
156: return z;
157: }
158: break;
159:
160: case t_MAT:
161: switch(ty)
162: {
163: case t_VEC:
164: z=cgetg(lx,tx); if (ly != lx) break;
165: for (i=1; i<lx; i++) z[i]=(long)concatsp((GEN) x[i],(GEN) y[i]);
166: return z;
167: case t_COL:
168: if (ly != lg(x[1])) break;
169: z=cgetg(lx+1,tx); z[lx]=(long)y;
170: for (i=1; i<lx; i++) z[i]=x[i];
171: return z;
172: }
173: break;
174: }
175: err(concater);
176: return NULL; /* not reached */
177: }
178:
179: GEN
180: concat(GEN x, GEN y)
181: {
182: long tx = typ(x), lx,ty,ly,i;
183: GEN z,p1;
184:
185: if (!y)
186: {
187: long av = avma, tetpil;
188: if (tx == t_LIST)
189: { lx = lgef(x); i = 2; }
190: else if (tx == t_VEC)
191: { lx = lg(x); i = 1; }
192: else err(concater);
193: if (i>=lx) err(talker,"trying to concat elements of an empty vector");
194: y = (GEN)x[i++];
195: for (; i<lx; i++) y = concatsp(y, (GEN)x[i]);
196: tetpil = avma; return gerepile(av,tetpil,gcopy(y));
197: }
198: ty = typ(y);
199: if (tx==t_LIST || ty==t_LIST) return listconcat(x,y);
200: if (tx==t_STR || ty==t_STR) return strconcat(x,y);
201: lx=lg(x); ly=lg(y);
202:
203: if (tx==t_MAT && lx==1)
204: {
205: if (ty!=t_VEC || ly==1) return gtomat(y);
206: err(concater);
207: }
208: if (ty==t_MAT && ly==1)
209: {
210: if (tx!=t_VEC || lx==1) return gtomat(x);
211: err(concater);
212: }
213:
214: if (! is_matvec_t(tx))
215: {
216: if (! is_matvec_t(ty))
217: {
218: z=cgetg(3,t_VEC); z[1]=lcopy(x); z[2]=lcopy(y);
219: return z;
220: }
221: z=cgetg(ly+1,ty);
222: if (ty != t_MAT) p1 = gcopy(x);
223: else
224: {
225: if (lg(y[1])!=2) err(concater);
226: p1=cgetg(2,t_COL); p1[1]=lcopy(x);
227: }
228: for (i=2; i<=ly; i++) z[i]=lcopy((GEN) y[i-1]);
229: z[1]=(long)p1; return z;
230: }
231: if (! is_matvec_t(ty))
232: {
233: z=cgetg(lx+1,tx);
234: if (tx != t_MAT) p1 = gcopy(y);
235: else
236: {
237: if (lg(x[1])!=2) err(concater);
238: p1=cgetg(2,t_COL); p1[1]=lcopy(y);
239: }
240: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
241: z[lx]=(long)p1; return z;
242: }
243:
244: if (tx == ty)
245: {
246: if (tx == t_MAT && lg(x[1]) != lg(y[1])) err(concater);
247: z=cgetg(lx+ly-1,tx);
248: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
249: for (i=1; i<ly; i++) z[lx+i-1]=lcopy((GEN) y[i]);
250: return z;
251: }
252:
253: switch(tx)
254: {
255: case t_VEC:
256: switch(ty)
257: {
258: case t_COL:
259: if (lx<=2) return (lx==1)? gcopy(y): concat((GEN) x[1],y);
260: if (ly>=3) break;
261: return (ly==1)? gcopy(x): concat(x,(GEN) y[1]);
262: case t_MAT:
263: z=cgetg(ly,ty); if (lx != ly) break;
264: for (i=1; i<ly; i++) z[i]=lconcat((GEN) x[i],(GEN) y[i]);
265: return z;
266: }
267: break;
268:
269: case t_COL:
270: switch(ty)
271: {
272: case t_VEC:
273: if (lx<=2) return (lx==1)? gcopy(y): concat((GEN) x[1],y);
274: if (ly>=3) break;
275: return (ly==1)? gcopy(x): concat(x,(GEN) y[1]);
276: case t_MAT:
277: if (lx != lg(y[1])) break;
278: z=cgetg(ly+1,ty); z[1]=lcopy(x);
279: for (i=2; i<=ly; i++) z[i]=lcopy((GEN) y[i-1]);
280: return z;
281: }
282: break;
283:
284: case t_MAT:
285: switch(ty)
286: {
287: case t_VEC:
288: z=cgetg(lx,tx); if (ly != lx) break;
289: for (i=1; i<lx; i++) z[i]=lconcat((GEN) x[i],(GEN) y[i]);
290: return z;
291: case t_COL:
292: if (ly != lg(x[1])) break;
293: z=cgetg(lx+1,tx); z[lx]=lcopy(y);
294: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
295: return z;
296: }
297: break;
298: }
299: err(concater);
300: return NULL; /* not reached */
301: }
302:
303: static long
304: str_to_long(char *s, char **pt)
305: {
306: long a = atol(s);
307: while (isspace((int)*s)) s++;
308: if (*s == '-' || *s == '+') s++;
309: while (isdigit((int)*s) || isspace((int)*s)) s++;
310: *pt = s; return a;
311: }
312:
313: static int
314: get_range(char *s, long *a, long *b, long *compl, long lx)
315: {
316: long max = lx - 1;
317:
318: *a = 1; *b = max;
319: if (*s == '^') { *compl = 1; s++; } else *compl = 0;
320: if (*s == 0) return 0;
321: if (*s != '.')
322: {
323: *a = str_to_long(s, &s);
324: if (*a < 0) *a += lx;
325: if (*a<1 || *a>max) return 0;
326: }
327: if (*s == '.')
328: {
329: s++; if (*s != '.') return 0;
330: do s++; while (isspace((int)*s));
331: if (*s)
332: {
333: *b = str_to_long(s, &s);
334: if (*b < 0) *b += lx;
335: if (*b<1 || *b>max || *s) return 0;
336: }
337: return 1;
338: }
339: if (*s) return 0;
340: *b = *a; return 1;
341: }
342:
343: GEN
344: extract(GEN x, GEN l)
345: {
346: long av,i,j, tl = typ(l), tx = typ(x), lx = lg(x);
347: GEN y;
348:
349: if (! is_matvec_t(tx)) err(typeer,"extract");
350: if (tl==t_INT)
351: {
352: /* extract components of x as per the bits of mask l */
353: if (!signe(l)) return cgetg(1,tx);
354: av=avma; y = (GEN) gpmalloc(lx*sizeof(long));
355: i = j = 1; while (!mpodd(l)) { l=shifti(l,-1); i++; }
356: while (signe(l) && i<lx)
357: {
358: if (mod2(l)) y[j++] = x[i];
359: i++; l=shifti(l,-1);
360: }
361: if (signe(l)) err(talker,"mask too large in vecextract");
362: y[0] = evaltyp(tx) | evallg(j);
363: avma=av; x = gcopy(y); free(y); return x;
364: }
365: if (tl==t_STR)
366: {
367: char *s = GSTR(l);
368: long first, last, compl;
369: if (! get_range(s, &first, &last, &compl, lx))
370: err(talker, "incorrect range in extract");
371: if (lx == 1) return gcopy(x);
372: if (compl)
373: {
374: if (first <= last)
375: {
376: y = cgetg(lx - (last - first + 1),tx);
377: for (j=1; j<first; j++) y[j] = lcopy((GEN)x[j]);
378: for (i=last+1; i<lx; i++,j++) y[j] = lcopy((GEN)x[i]);
379: }
380: else
381: {
382: y = cgetg(lx - (first - last + 1),tx);
383: for (j=1,i=lx-1; i>first; i--,j++) y[j] = lcopy((GEN)x[i]);
384: for (i=last-1; i>0; i--,j++) y[j] = lcopy((GEN)x[i]);
385: }
386: }
387: else
388: {
389: if (first <= last)
390: {
391: y = cgetg(last-first+2,tx);
392: for (i=first,j=1; i<=last; i++,j++) y[j] = lcopy((GEN)x[i]);
393: }
394: else
395: {
396: y = cgetg(first-last+2,tx);
397: for (i=first,j=1; i>=last; i--,j++) y[j] = lcopy((GEN)x[i]);
398: }
399: }
400: return y;
401: }
402:
403: if (is_vec_t(tl))
404: {
405: long ll=lg(l); y=cgetg(ll,tx);
406: for (i=1; i<ll; i++)
407: {
408: j = itos((GEN) l[i]);
409: if (j>=lx || j<=0) err(talker,"no such component in vecextract");
410: y[i] = lcopy((GEN) x[j]);
411: }
412: return y;
413: }
414: if (tl == t_VECSMALL)
415: {
416: long ll=lg(l); y=cgetg(ll,tx);
417: for (i=1; i<ll; i++)
418: {
419: j = l[i];
420: if (j>=lx || j<=0) err(talker,"no such component in vecextract");
421: y[i] = lcopy((GEN) x[j]);
422: }
423: return y;
424: }
425: err(talker,"incorrect mask in vecextract");
426: return NULL; /* not reached */
427: }
428:
429: GEN
430: matextract(GEN x, GEN l1, GEN l2)
431: {
432: long av = avma, tetpil;
433:
434: if (typ(x)!=t_MAT) err(typeer,"matextract");
435: x = extract(gtrans(extract(x,l2)),l1); tetpil=avma;
436: return gerepile(av,tetpil, gtrans(x));
437: }
438:
439: GEN
440: extract0(GEN x, GEN l1, GEN l2)
441: {
442: if (! l2) return extract(x,l1);
443: return matextract(x,l1,l2);
444: }
445:
446: /*******************************************************************/
447: /* */
448: /* SCALAR-MATRIX OPERATIONS */
449: /* */
450: /*******************************************************************/
451:
452: /* create the square nxn matrix equal to z*Id */
453: static GEN
454: gscalmat_proto(GEN z, GEN myzero, long n, int flag)
455: {
456: long i,j;
457: GEN y = cgetg(n+1,t_MAT);
458: if (n < 0) err(talker,"negative size in scalmat");
459: if (flag) z = (flag==1)? stoi((long)z): gcopy(z);
460: for (i=1; i<=n; i++)
461: {
462: y[i]=lgetg(n+1,t_COL);
463: for (j=1; j<=n; j++)
464: coeff(y,j,i) = (i==j)? (long)z: (long)myzero;
465: }
466: return y;
467: }
468:
469: GEN
470: gscalmat(GEN x, long n) { return gscalmat_proto(x,gzero,n,2); }
471:
472: GEN
473: gscalsmat(long x, long n) { return gscalmat_proto((GEN)x,gzero,n,1); }
474:
475: GEN
476: idmat(long n) { return gscalmat_proto(gun,gzero,n,0); }
477:
478: GEN
479: idmat_intern(long n,GEN myun,GEN z) { return gscalmat_proto(myun,z,n,0); }
480:
481: GEN
482: gscalcol_proto(GEN z, GEN myzero, long n)
483: {
484: GEN y = cgetg(n+1,t_COL);
485: long i;
486:
487: if (n)
488: {
489: y[1]=(long)z;
490: for (i=2; i<=n; i++) y[i]=(long)myzero;
491: }
492: return y;
493: }
494:
495: GEN
496: zerocol(long n)
497: {
498: GEN y = cgetg(n+1,t_COL);
499: long i;
500: for (i=1; i<=n; i++) y[i]=zero;
501: return y;
502: }
503:
504: GEN
505: gscalcol(GEN x, long n) { return gscalcol_proto(gcopy(x),gzero,n); }
506:
507: GEN
508: gscalcol_i(GEN x, long n) { return gscalcol_proto(x,gzero,n); }
509:
510: GEN
511: gtomat(GEN x)
512: {
513: long tx,lx,i;
514: GEN y,p1;
515:
516: if (!x) return cgetg(1, t_MAT);
517: tx = typ(x);
518: if (! is_matvec_t(tx))
519: {
520: y=cgetg(2,t_MAT); p1=cgetg(2,t_COL); y[1]=(long)p1;
521: p1[1]=lcopy(x); return y;
522: }
523: switch(tx)
524: {
525: case t_VEC:
526: lx=lg(x); y=cgetg(lx,t_MAT);
527: for (i=1; i<lx; i++)
528: {
529: p1=cgetg(2,t_COL); y[i]=(long)p1;
530: p1[1]=lcopy((GEN) x[i]);
531: }
532: break;
533: case t_COL:
534: y=cgetg(2,t_MAT); y[1]=lcopy(x); break;
535: case t_MAT:
536: y=gcopy(x); break;
537: }
538: return y;
539: }
540:
541: long
542: isdiagonal(GEN x)
543: {
544: long nco,i,j;
545:
546: if (typ(x)!=t_MAT) err(typeer,"isdiagonal");
547: nco=lg(x)-1; if (!nco) return 1;
548: if (nco != lg(x[1])-1) return 0;
549:
550: for (j=1; j<=nco; j++)
551: {
552: GEN *col = (GEN*) x[j];
553: for (i=1; i<=nco; i++)
554: if (i!=j && !gcmp0(col[i])) return 0;
555: }
556: return 1;
557: }
558:
559: /* create the diagonal matrix, whose diagonal is given by x */
560: GEN
561: diagonal(GEN x)
562: {
563: long i,j,lx,tx=typ(x);
564: GEN y,p1;
565:
566: if (! is_matvec_t(tx)) return gscalmat(x,1);
567: if (tx==t_MAT)
568: {
569: if (isdiagonal(x)) return gcopy(x);
570: err(talker,"incorrect object in diagonal");
571: }
572: lx=lg(x); y=cgetg(lx,t_MAT);
573: for (j=1; j<lx; j++)
574: {
575: p1=cgetg(lx,t_COL); y[j]=(long)p1;
576: for (i=1; i<lx; i++)
577: p1[i] = (i==j)? lcopy((GEN) x[i]): zero;
578: }
579: return y;
580: }
581:
582: /* compute m*diagonal(d) */
583: GEN
584: matmuldiagonal(GEN m, GEN d)
585: {
586: long j=typ(d),lx=lg(m);
587: GEN y;
588:
589: if (typ(m)!=t_MAT) err(typeer,"matmuldiagonal");
590: if (! is_vec_t(j) || lg(d)!=lx)
591: err(talker,"incorrect vector in matmuldiagonal");
592: y=cgetg(lx,t_MAT);
593: for (j=1; j<lx; j++) y[j] = lmul((GEN) d[j],(GEN) m[j]);
594: return y;
595: }
596:
597: /* compute m*n assuming the result is a diagonal matrix */
598: GEN
599: matmultodiagonal(GEN m, GEN n)
600: {
601: long lx,i,j;
602: GEN s,y;
603:
604: if (typ(m)!=t_MAT || typ(n)!=t_MAT) err(typeer,"matmultodiagonal");
605: lx=lg(n); y=idmat(lx-1);
606: if (lx == 1)
607: { if (lg(m) != 1) err(consister,"matmultodiagonal"); }
608: else
609: { if (lg(m) != lg(n[1])) err(consister,"matmultodiagonal"); }
610: for (i=1; i<lx; i++)
611: {
612: s = gzero;
613: for (j=1; j<lx; j++)
614: s = gadd(s,gmul(gcoeff(m,i,j),gcoeff(n,j,i)));
615: coeff(y,i,i) = (long)s;
616: }
617: return y;
618: }
619:
620: /* [m[1,1], ..., m[l,l]] */
621: GEN
622: mattodiagonal(GEN m)
623: {
624: long i, lx = lg(m);
625: GEN y = cgetg(lx,t_VEC);
626:
627: if (typ(m)!=t_MAT) err(typeer,"mattodiagonal");
628: if (lx == 1) return y;
629: for (i=1; i<lx; i++) y[i] = lcopy(gcoeff(m,i,i));
630: return y;
631: }
632:
633: /*******************************************************************/
634: /* */
635: /* ADDITION SCALAR + MATRIX */
636: /* */
637: /*******************************************************************/
638:
639: /* create the square matrix x*Id + y */
640: GEN
641: gaddmat(GEN x, GEN y)
642: {
643: long ly,dy,i,j;
644: GEN z;
645:
646: ly=lg(y); if (ly==1) err(gadderf,"Scalar","t_MAT");
647: dy=lg(y[1]);
648: if (typ(y)!=t_MAT || ly!=dy) err(mattype1,"gaddmat");
649: z=cgetg(ly,t_MAT);
650: for (i=1; i<ly; i++)
651: {
652: z[i]=lgetg(dy,t_COL);
653: for (j=1; j<dy; j++)
654: coeff(z,j,i) = i==j? ladd(x,gcoeff(y,j,i)): lcopy(gcoeff(y,j,i));
655: }
656: return z;
657: }
658:
659: /*******************************************************************/
660: /* */
661: /* Solve A*X=B (Gauss pivot) */
662: /* */
663: /*******************************************************************/
664: #define swap(x,y) { long _t=x; x=y; y=_t; }
665:
666: /* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be
667: * used, and the matrix precision (min real precision of coeffs) otherwise.
668: */
669: static long
670: matprec(GEN x)
671: {
672: long tx,i,j,l, k = VERYBIGINT, lx = lg(x), ly = lg(x[1]);
673: GEN p1;
674: for (i=1; i<lx; i++)
675: for (j=1; j<ly; j++)
676: {
677: p1 = gmael(x,i,j); tx = typ(p1);
678: if (!is_scalar_t(tx)) return 0;
679: l = precision(p1); if (l && l<k) k = l;
680: }
681: return (k==VERYBIGINT)? 0: k;
682: }
683:
684: /* As above, returning 1 if the precision would be non-zero, 0 otherwise */
685: static long
686: use_maximal_pivot(GEN x)
687: {
688: long tx,i,j, lx = lg(x), ly = lg(x[1]);
689: GEN p1;
690: for (i=1; i<lx; i++)
691: for (j=1; j<ly; j++)
692: {
693: p1 = gmael(x,i,j); tx = typ(p1);
694: if (!is_scalar_t(tx)) return 0;
695: if (precision(p1)) return 1;
696: }
697: return 0;
698: }
699:
700: static GEN
701: check_b(GEN b, long nbli)
702: {
703: GEN col;
704: if (!b) return idmat(nbli);
705: b = dummycopy(b);
706: col = (typ(b) == t_MAT)? (GEN)b[1]: b;
707: if (nbli == lg(col)-1) return b;
708: err(talker,"incompatible matrix dimensions in gauss");
709: return NULL; /* not reached */
710: }
711:
712: GEN
713: gauss_get_col(GEN a, GEN b, GEN p, long nbli)
714: {
715: GEN m, u=cgetg(nbli+1,t_COL);
716: long i,j;
717:
718: u[nbli] = ldiv((GEN) b[nbli],p);
719: for (i=nbli-1; i>0; i--)
720: {
721: m = gneg_i((GEN)b[i]);
722: for (j=i+1; j<=nbli; j++)
723: m = gadd(m, gmul(gcoeff(a,i,j),(GEN) u[j]));
724: u[i] = ldiv(gneg_i(m), gcoeff(a,i,i));
725: }
726: return u;
727: }
728:
729: /* Gauss pivot.
730: * Compute a^(-1)*b, where nblig(a) = nbcol(a) = nblig(b).
731: * b is a matrix or column vector, NULL meaning: take the identity matrix
732: * Be careful, if a or b is empty, the result is the empty matrix...
733: */
734: GEN
735: gauss(GEN a, GEN b)
736: {
737: long inexact,ismat,nbli,nbco,i,j,k,av,tetpil,lim;
738: GEN p,m,u;
739: /* nbli: nb lines of b = nb columns of a */
740: /* nbco: nb columns of b (if matrix) */
741:
742: if (typ(a)!=t_MAT) err(mattype1,"gauss");
743: if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss");
744: if (lg(a) == 1)
745: {
746: if (b && lg(b)!=1) err(consister,"gauss");
747: if (DEBUGLEVEL)
748: err(warner,"in Gauss lg(a)=%ld lg(b)=%ld",lg(a),b?lg(b):-1);
749: return cgetg(1,t_MAT);
750: }
751: av=avma; lim=stack_lim(av,1);
752: nbli = lg(a)-1; if (nbli!=lg(a[1])-1) err(mattype1,"gauss");
753: a = dummycopy(a);
754: b = check_b(b,nbli);
755: nbco = lg(b)-1;
756: inexact = use_maximal_pivot(a);
757: ismat = (typ(b)==t_MAT);
758: if(DEBUGLEVEL>4)
759: fprintferr("Entering gauss with inexact=%ld ismat=%ld\n",inexact,ismat);
760:
761: for (i=1; i<nbli; i++)
762: {
763: /* k is the line where we find the pivot */
764: p=gcoeff(a,i,i); k=i;
765: if (inexact) /* maximal pivot */
766: {
767: long e, ex = gexpo(p);
768: for (j=i+1; j<=nbli; j++)
769: {
770: e = gexpo(gcoeff(a,j,i));
771: if (e > ex) { ex=e; k=j; }
772: }
773: if (gcmp0(gcoeff(a,k,i))) err(matinv1);
774: }
775: else if (gcmp0(p)) /* first non-zero pivot */
776: {
777: do k++; while (k<=nbli && gcmp0(gcoeff(a,k,i)));
778: if (k>nbli) err(matinv1);
779: }
780:
781: /* if (k!=i), exchange the lines s.t. k = i */
782: if (k != i)
783: {
784: for (j=i; j<=nbli; j++) swap(coeff(a,i,j), coeff(a,k,j));
785: if (ismat)
786: {
787: for (j=1; j<=nbco; j++) swap(coeff(b,i,j), coeff(b,k,j));
788: }
789: else
790: swap(b[i],b[k]);
791: p = gcoeff(a,i,i);
792: }
793:
794: for (k=i+1; k<=nbli; k++)
795: {
796: m=gcoeff(a,k,i);
797: if (!gcmp0(m))
798: {
799: m = gneg_i(gdiv(m,p));
800: for (j=i+1; j<=nbli; j++)
801: {
802: u = gmul(m,gcoeff(a,i,j));
803: coeff(a,k,j) = ladd(gcoeff(a,k,j),u);
804: }
805: if (ismat) for (j=1; j<=nbco; j++)
806: {
807: u = gmul(m,gcoeff(b,i,j));
808: coeff(b,k,j) = ladd(gcoeff(b,k,j),u);
809: }
810: else
811: {
812: u = gmul(m,(GEN) b[i]);
813: b[k] = ladd((GEN) b[k],u);
814: }
815: }
816: }
817: if (low_stack(lim, stack_lim(av,1)))
818: {
819: GEN *gptr[2];
820: if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i);
821: gptr[0]=&a; gptr[1]=&b;
822: gerepilemany(av,gptr,2);
823: }
824: }
825:
826: if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
827: p=gcoeff(a,nbli,nbli);
828: if (!inexact && gcmp0(p)) err(matinv1);
829: if (!ismat) u = gauss_get_col(a,b,p,nbli);
830: else
831: {
832: long av1 = avma;
833: lim = stack_lim(av1,1); u=cgetg(nbco+1,t_MAT);
834: for (j=2; j<=nbco; j++) u[j] = zero; /* dummy */
835: for (j=1; j<=nbco; j++)
836: {
837: u[j] = (long)gauss_get_col(a,(GEN)b[j],p,nbli);
838: if (low_stack(lim, stack_lim(av1,1)))
839: {
840: if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j);
841: tetpil=avma; u = gerepile(av1,tetpil,gcopy(u));
842: }
843: }
844: }
845: tetpil=avma; return gerepile(av,tetpil,gcopy(u));
846: }
847:
848: /* x a matrix with integer coefficients. Return a multiple of the determinant
849: * of the lattice generated by the columns of x (to be used with hnfmod)
850: */
851: GEN
852: detint(GEN x)
853: {
854: GEN pass,c,v,det1,piv,pivprec,vi,p1;
855: long i,j,k,rg,n,m,m1,av=avma,av1,lim;
856:
857: if (typ(x)!=t_MAT) err(typeer,"detint");
858: n=lg(x)-1; if (!n) return gun;
859: m1=lg(x[1]); m=m1-1; lim=stack_lim(av,1);
860: c=new_chunk(m1); for (k=1; k<=m; k++) c[k]=0;
861: av1=avma; pass=cgetg(m1,t_MAT);
862: for (j=1; j<=m; j++)
863: {
864: p1=cgetg(m1,t_COL); pass[j]=(long)p1;
865: for (i=1; i<=m; i++) p1[i]=zero;
866: }
867: v=cgetg(m1,t_COL);
868: det1=gzero; piv=pivprec=gun;
869: for (rg=0,k=1; k<=n; k++)
870: {
871: long t = 0;
872: for (i=1; i<=m; i++)
873: if (!c[i])
874: {
875: vi=mulii(piv,gcoeff(x,i,k));
876: for (j=1; j<=m; j++)
877: if (c[j]) vi=addii(vi,mulii(gcoeff(pass,i,j),gcoeff(x,j,k)));
878: v[i]=(long)vi; if (!t && signe(vi)) t=i;
879: }
880: if (t)
881: {
882: if (rg == m-1)
883: { det1=mppgcd((GEN)v[t],det1); c[t]=0; }
884: else
885: {
886: rg++; pivprec = piv; piv=(GEN)v[t]; c[t]=k;
887: for (i=1; i<=m; i++)
888: if (!c[i])
889: {
890: GEN p2 = negi((GEN)v[i]);
891: for (j=1; j<=m; j++)
892: if (c[j] && j!=t)
893: {
894: p1 = addii(mulii(gcoeff(pass,i,j), piv),
895: mulii(gcoeff(pass,t,j), p2));
896: if (rg>1) p1 = divii(p1,pivprec);
897: coeff(pass,i,j) = (long)p1;
898: }
899: coeff(pass,i,t) = (long)p2;
900: }
901: }
902: }
903: if (low_stack(lim, stack_lim(av,1)))
904: {
905: GEN *gptr[5];
906: if(DEBUGMEM>1) err(warnmem,"detint. k=%ld",k);
907: gptr[0]=&det1; gptr[1]=ϖ gptr[2]=&pivprec;
908: gptr[3]=&pass; gptr[4]=&v; gerepilemany(av1,gptr,5);
909: }
910: }
911: return gerepileupto(av, absi(det1));
912: }
913:
914: static void
915: gerepile_gauss_keep(GEN x, long m, long n, long k, long t, long av)
916: {
917: long tetpil = avma,dec,u,A,i;
918:
919: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_keep. k=%ld, n=%ld",k,n);
920: for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k));
921: for (i=k+1; i<=n; i++)
922: for (u=1; u<=m; u++) copyifstack(coeff(x,u,i), coeff(x,u,i));
923:
924: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
925: for (u=t+1; u<=m; u++)
926: {
927: A=coeff(x,u,k);
928: if (A<av && A>=bot) coeff(x,u,k)+=dec;
929: }
930: for (i=k+1; i<=n; i++)
931: for (u=1; u<=m; u++)
932: {
933: A=coeff(x,u,i);
934: if (A<av && A>=bot) coeff(x,u,i)+=dec;
935: }
936: }
937:
938: static void
939: gerepile_gauss_keep_mod_p(GEN x, GEN p, long m, long n, long k, long t, long av)
940: {
941: long tetpil = avma,dec,u,A,i;
942:
943: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_keep. k=%ld, n=%ld",k,n);
944: for (u=t+1; u<=m; u++)
945: if (isonstack(coeff(x,u,k))) coeff(x,u,k) = lmodii(gcoeff(x,u,k),p);
946: for (i=k+1; i<=n; i++)
947: for (u=1; u<=m; u++)
948: if (isonstack(coeff(x,u,i))) coeff(x,u,i) = lmodii(gcoeff(x,u,i),p);
949:
950: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
951: for (u=t+1; u<=m; u++)
952: {
953: A=coeff(x,u,k);
954: if (A<av && A>=bot) coeff(x,u,k)+=dec;
955: }
956: for (i=k+1; i<=n; i++)
957: for (u=1; u<=m; u++)
958: {
959: A=coeff(x,u,i);
960: if (A<av && A>=bot) coeff(x,u,i)+=dec;
961: }
962: }
963:
964: /* special gerepile for huge matrices */
965:
966: static void
967: gerepile_gauss(GEN x,long m,long n,long k,long t,long av, long j, GEN c)
968: {
969: long tetpil = avma,dec,u,A,i;
970:
971: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n);
972: for (u=t+1; u<=m; u++)
973: if (u==j || !c[u]) copyifstack(coeff(x,u,k), coeff(x,u,k));
974: for (u=1; u<=m; u++)
975: if (u==j || !c[u])
976: for (i=k+1; i<=n; i++) copyifstack(coeff(x,u,i), coeff(x,u,i));
977:
978: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
979: for (u=t+1; u<=m; u++)
980: if (u==j || !c[u])
981: {
982: A=coeff(x,u,k);
983: if (A<av && A>=bot) coeff(x,u,k)+=dec;
984: }
985: for (u=1; u<=m; u++)
986: if (u==j || !c[u])
987: for (i=k+1; i<=n; i++)
988: {
989: A=coeff(x,u,i);
990: if (A<av && A>=bot) coeff(x,u,i)+=dec;
991: }
992: }
993:
994: /*******************************************************************/
995: /* */
996: /* KERNEL of an m x n matrix */
997: /* return n - rk(x) linearly independant vectors */
998: /* */
999: /*******************************************************************/
1000:
1001: /* x has INTEGER coefficients */
1002: GEN
1003: keri(GEN x)
1004: {
1005: GEN c,d,y,p,pp;
1006: long i,j,k,r,t,n,m,av,av0,tetpil,lim;
1007:
1008: if (typ(x)!=t_MAT) err(typeer,"keri");
1009: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
1010:
1011: av0=avma; m=lg(x[1])-1; r=0;
1012: pp=cgetg(n+1,t_COL);
1013: x=dummycopy(x); p=gun;
1014: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1015: d=new_chunk(n+1); av=avma; lim=stack_lim(av,1);
1016: for (k=1; k<=n; k++)
1017: {
1018: j=1;
1019: while (j<=m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
1020: if (j>m)
1021: {
1022: r++; d[k]=0;
1023: for(j=1; j<k; j++)
1024: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
1025: pp[k]=lclone(p);
1026: }
1027: else
1028: {
1029: GEN p0 = p;
1030: long av1;
1031:
1032: c[j]=k; d[k]=j; p = gcoeff(x,j,k);
1033:
1034: for (t=1; t<=m; t++)
1035: if (t!=j)
1036: {
1037: GEN q=gcoeff(x,t,k), p1,p2;
1038: for (i=k+1; i<=n; i++)
1039: {
1040: av1=avma; (void)new_chunk(lgefint(p0));
1041: p1=mulii(q,gcoeff(x,j,i));
1042: p2=mulii(p,gcoeff(x,t,i));
1043: p1=subii(p2,p1); avma=av1;
1044: coeff(x,t,i) = ldivii(p1,p0);
1045: }
1046: if (low_stack(lim, stack_lim(av,1)))
1047: {
1048: p1 = gclone(p);
1049: gerepile_gauss_keep(x,m,n,k,t,av);
1050: p = gcopy(p1); gunclone(p1);
1051: }
1052: }
1053: }
1054: }
1055: if (!r) { avma=av0; y=cgetg(1,t_MAT); return y; }
1056:
1057: /* non trivial kernel */
1058: tetpil=avma; y=cgetg(r+1,t_MAT);
1059: for (j=k=1; j<=r; j++,k++)
1060: {
1061: p = cgetg(n+1, t_COL);
1062: y[j]=(long)p; while (d[k]) k++;
1063: for (i=1; i<k; i++)
1064: if (d[i])
1065: {
1066: c=gcoeff(x,d[i],k);
1067: p[i] = (long) forcecopy(c); gunclone(c);
1068: }
1069: else
1070: p[i] = zero;
1071: p[k]=lnegi((GEN)pp[k]); gunclone((GEN)pp[k]);
1072: for (i=k+1; i<=n; i++) p[i]=zero;
1073: }
1074: return gerepile(av0,tetpil,y);
1075: }
1076:
1077: GEN
1078: deplin(GEN x)
1079: {
1080: long i,j,k,t,nc,nl, av=avma;
1081: GEN y,q,c,l,d;
1082:
1083: if (typ(x) != t_MAT) err(typeer,"deplin");
1084: nc=lg(x)-1; if (!nc) err(talker,"empty matrix in deplin");
1085: nl=lg(x[1])-1;
1086: c=new_chunk(nl+1);
1087: l=new_chunk(nc+1);
1088: d=cgetg(nl+1,t_VEC);
1089: for (i=1; i<=nl; i++) { d[i]=un; c[i]=0; }
1090: k=1; t=1;
1091: while (t<=nl && k<=nc)
1092: {
1093: for (j=1; j<k; j++)
1094: for (i=1; i<=nl; i++)
1095: if (i!=l[j])
1096: coeff(x,i,k)=lsub(gmul((GEN) d[j],gcoeff(x,i,k)),
1097: gmul(gcoeff(x,i,j),gcoeff(x,l[j],k)));
1098: t=1;
1099: while ( t<=nl && (c[t] || gcmp0(gcoeff(x,t,k))) ) t++;
1100: if (t<=nl)
1101: {
1102: d[k]=coeff(x,t,k);
1103: c[t]=k; l[k++]=t;
1104: }
1105: }
1106: if (k>nc)
1107: {
1108: avma=av; y=cgetg(nc+1,t_COL);
1109: for (j=1; j<=nc; j++) y[j]=zero;
1110: return y;
1111: }
1112: y=cgetg(nc+1,t_COL);
1113: y[1]=(k>1)? coeff(x,l[1],k): un;
1114: for (q=gun,j=2; j<k; j++)
1115: {
1116: q=gmul(q,(GEN) d[j-1]);
1117: y[j]=lmul(gcoeff(x,l[j],k),q);
1118: }
1119: if (k>1) y[k]=lneg(gmul(q,(GEN) d[k-1]));
1120: for (j=k+1; j<=nc; j++) y[j]=zero;
1121: d=content(y); t=avma;
1122: return gerepile(av,t,gdiv(y,d));
1123: }
1124:
1125: /*******************************************************************/
1126: /* */
1127: /* GAUSS REDUCTION OF MATRICES (m lines x n cols) */
1128: /* (kernel, image, complementary image, rank) */
1129: /* */
1130: /*******************************************************************/
1131: static long gauss_ex;
1132: static int (*gauss_is_zero)(GEN);
1133:
1134: static int
1135: real0(GEN x)
1136: {
1137: return gcmp0(x) || (gexpo(x) < gauss_ex);
1138: }
1139:
1140: static void
1141: gauss_get_prec(GEN x, long prec)
1142: {
1143: long pr = matprec(x);
1144:
1145: if (!pr) { gauss_is_zero = &gcmp0; return; }
1146: if (pr > prec) prec = pr;
1147:
1148: gauss_ex = - (long)(0.85 * bit_accuracy(prec));
1149: gauss_is_zero = &real0;
1150: }
1151:
1152: /* return the transform of x under a standard Gauss pivot. r = dim ker(x).
1153: * d[k] contains the index of the first non-zero pivot in column k
1154: */
1155: static GEN
1156: gauss_pivot_keep(GEN x, long prec, GEN *dd, long *rr)
1157: {
1158: GEN c,d,p,mun;
1159: long i,j,k,r,t,n,m,av,lim;
1160:
1161: if (typ(x)!=t_MAT) err(typeer,"gauss_pivot");
1162: n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
1163:
1164: gauss_get_prec(x,prec); m=lg(x[1])-1; r=0;
1165: x=dummycopy(x); mun=negi(gun);
1166: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1167: d=(GEN)gpmalloc((n+1)*sizeof(long));
1168: av=avma; lim=stack_lim(av,1);
1169: for (k=1; k<=n; k++)
1170: {
1171: j=1; while (j<=m && (c[j] || gauss_is_zero(gcoeff(x,j,k)))) j++;
1172: if (j>m)
1173: {
1174: r++; d[k]=0;
1175: for(j=1; j<k; j++)
1176: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
1177: }
1178: else
1179: {
1180: c[j]=k; d[k]=j; p = gdiv(mun,gcoeff(x,j,k));
1181: coeff(x,j,k) = (long)mun;
1182: for (i=k+1; i<=n; i++)
1183: coeff(x,j,i) = lmul(p,gcoeff(x,j,i));
1184: for (t=1; t<=m; t++)
1185: if (t!=j)
1186: {
1187: p=gcoeff(x,t,k); coeff(x,t,k)=zero;
1188: for (i=k+1; i<=n; i++)
1189: coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
1190: if (low_stack(lim, stack_lim(av,1)))
1191: gerepile_gauss_keep(x,m,n,k,t,av);
1192: }
1193: }
1194: }
1195: *dd=d; *rr=r; return x;
1196: }
1197:
1198: /* r = dim ker(x).
1199: * d[k] contains the index of the first non-zero pivot in column k
1200: */
1201: static void
1202: gauss_pivot(GEN x, long prec, GEN *dd, long *rr)
1203: {
1204: GEN c,d,mun,p;
1205: long i,j,k,r,t,n,m,av,lim;
1206:
1207: if (typ(x)!=t_MAT) err(typeer,"gauss_pivot");
1208: n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; }
1209:
1210: gauss_get_prec(x,prec); m=lg(x[1])-1; r=0;
1211: x=dummycopy(x); mun=negi(gun);
1212: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1213: d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1);
1214: for (k=1; k<=n; k++)
1215: {
1216: j=1; while (j<=m && (c[j] || gauss_is_zero(gcoeff(x,j,k)))) j++;
1217: if (j>m) { r++; d[k]=0; }
1218: else
1219: {
1220: c[j]=k; d[k]=j; p = gdiv(mun,gcoeff(x,j,k));
1221: for (i=k+1; i<=n; i++)
1222: coeff(x,j,i) = lmul(p,gcoeff(x,j,i));
1223:
1224: for (t=1; t<=m; t++)
1225: if (!c[t]) /* no pivot on that line yet */
1226: {
1227: p=gcoeff(x,t,k); coeff(x,t,k)=zero;
1228: for (i=k+1; i<=n; i++)
1229: coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i)));
1230: if (low_stack(lim, stack_lim(av,1)))
1231: gerepile_gauss(x,m,n,k,t,av,j,c);
1232: }
1233: for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
1234: }
1235: }
1236: *dd=d; *rr=r;
1237: }
1238:
1239: static GEN
1240: ker0(GEN x, long prec)
1241: {
1242: GEN d,y;
1243: long i,j,k,r,n, av = avma, tetpil;
1244:
1245: x=gauss_pivot_keep(x,prec,&d,&r);
1246: if (!r)
1247: {
1248: avma=av; if (d) free(d);
1249: return cgetg(1,t_MAT);
1250: }
1251: n = lg(x)-1; tetpil=avma; y=cgetg(r+1,t_MAT);
1252: for (j=k=1; j<=r; j++,k++)
1253: {
1254: GEN p = cgetg(n+1,t_COL);
1255:
1256: y[j]=(long)p; while (d[k]) k++;
1257: for (i=1; i<k; i++)
1258: if (d[i])
1259: {
1260: GEN p1=gcoeff(x,d[i],k);
1261: p[i] = (long)forcecopy(p1); gunclone(p1);
1262: }
1263: else
1264: p[i] = zero;
1265: p[k]=un; for (i=k+1; i<=n; i++) p[i]=zero;
1266: }
1267: free(d); return gerepile(av,tetpil,y);
1268: }
1269:
1270: GEN
1271: ker(GEN x) /* Programme pour types exacts */
1272: {
1273: return ker0(x,0);
1274: }
1275:
1276: GEN
1277: matker0(GEN x,long flag)
1278: {
1279: return flag? keri(x): ker(x);
1280: }
1281:
1282: static GEN
1283: image0(GEN x, long prec)
1284: {
1285: GEN d,y;
1286: long j,k,r, av = avma;
1287:
1288: gauss_pivot(x,prec,&d,&r);
1289:
1290: /* r = dim ker(x) */
1291: if (!r) { avma=av; if (d) free(d); return gcopy(x); }
1292:
1293: /* r = dim Im(x) */
1294: r = lg(x)-1 - r; avma=av;
1295: y=cgetg(r+1,t_MAT);
1296: for (j=k=1; j<=r; k++)
1297: if (d[k]) y[j++] = lcopy((GEN)x[k]);
1298: free(d); return y;
1299: }
1300:
1301: GEN
1302: image(GEN x) /* Programme pour types exacts */
1303: {
1304: return image0(x,0);
1305: }
1306:
1307: GEN
1308: imagereel(GEN x, long prec) /* Programme pour types inexacts */
1309: {
1310: return image0(x,prec);
1311: }
1312:
1313: static GEN
1314: imagecompl0(GEN x, long prec)
1315: {
1316: GEN d,y;
1317: long j,i,r,av = avma;
1318:
1319: gauss_pivot(x,prec,&d,&r);
1320: avma=av; y=cgetg(r+1,t_VEC);
1321: for (i=j=1; j<=r; i++)
1322: if (!d[i]) y[j++]=lstoi(i);
1323: if (d) free(d); return y;
1324: }
1325:
1326: /* for hnfspec: imagecompl(trans(x)) + image(trans(x)) */
1327: static GEN
1328: imagecomplspec(GEN x, long *nlze)
1329: {
1330: GEN d,y;
1331: long i,j,k,l,r,av = avma;
1332:
1333: x = gtrans(x); l = lg(x);
1334: gauss_pivot(x,0,&d,&r);
1335: avma=av; y = cgetg(l,t_VECSMALL);
1336: for (i=j=1, k=r+1; i<l; i++)
1337: if (d[i]) y[k++]=i; else y[j++]=i;
1338: *nlze = r;
1339: if (d) free(d); return y;
1340: }
1341:
1342: GEN
1343: imagecompl(GEN x) /* Programme pour types exacts */
1344: {
1345: return imagecompl0(x,0);
1346: }
1347:
1348: static GEN
1349: sinverseimage(GEN mat, GEN y)
1350: {
1351: long av=avma,tetpil,i, nbcol = lg(mat);
1352: GEN col,p1 = cgetg(nbcol+1,t_MAT);
1353:
1354: if (nbcol==1) return NULL;
1355: if (lg(y) != lg(mat[1])) err(consister,"inverseimage");
1356:
1357: p1[nbcol] = (long)y;
1358: for (i=1; i<nbcol; i++) p1[i]=mat[i];
1359: p1 = ker(p1); i=lg(p1)-1;
1360: if (!i) return NULL;
1361:
1362: col = (GEN)p1[i]; p1 = (GEN) col[nbcol];
1363: if (gcmp0(p1)) return NULL;
1364:
1365: p1 = gneg_i(p1); setlg(col,nbcol); tetpil=avma;
1366: return gerepile(av,tetpil, gdiv(col, p1));
1367: }
1368:
1369: /* Calcule l'image reciproque de v par m */
1370: GEN
1371: inverseimage(GEN m,GEN v)
1372: {
1373: long av=avma,j,lv,tv=typ(v);
1374: GEN y,p1;
1375:
1376: if (typ(m)!=t_MAT) err(typeer,"inverseimage");
1377: if (tv==t_COL)
1378: {
1379: p1 = sinverseimage(m,v);
1380: if (p1) return p1;
1381: avma = av; return cgetg(1,t_MAT);
1382: }
1383: if (tv!=t_MAT) err(typeer,"inverseimage");
1384:
1385: lv=lg(v)-1; y=cgetg(lv+1,t_MAT);
1386: for (j=1; j<=lv; j++)
1387: {
1388: p1 = sinverseimage(m,(GEN)v[j]);
1389: if (!p1) { avma = av; return cgetg(1,t_MAT); }
1390: y[j] = (long)p1;
1391: }
1392: return y;
1393: }
1394:
1395: /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
1396: * whose first k columns are given by x. If rank(x)<k, the result may be wrong
1397: */
1398: GEN
1399: suppl_intern(GEN x, GEN myid)
1400: {
1401: long av = avma, lx = lg(x), n,i,j;
1402: GEN y,p1;
1403: stackzone *zone;
1404:
1405: if (typ(x) != t_MAT) err(typeer,"suppl");
1406: if (lx==1) err(talker,"empty matrix in suppl");
1407: n=lg(x[1]); if (lx>n) err(suppler2);
1408:
1409: zone = switch_stack(NULL, n*n);
1410: switch_stack(zone,1);
1411: y = myid? dummycopy(myid): idmat(n-1);
1412: switch_stack(zone,0);
1413: for (i=1; i<lx; i++)
1414: {
1415: p1=gauss(y,(GEN)x[i]); j=i;
1416: while (j<n && gcmp0((GEN)p1[j])) j++;
1417: if (j>=n) err(suppler2);
1418: p1=(GEN)y[i]; y[i]=x[i]; if (i!=j) y[j]=(long)p1;
1419: }
1420: avma = av; y = gcopy(y);
1421: free(zone); return y;
1422: }
1423:
1424: GEN
1425: suppl(GEN x)
1426: {
1427: return suppl_intern(x,NULL);
1428: }
1429:
1430: GEN
1431: image2(GEN x)
1432: {
1433: long av=avma,tetpil,k,n,i;
1434: GEN p1,p2;
1435:
1436: if (typ(x)!=t_MAT) err(typeer,"image2");
1437: k=lg(x)-1; if (!k) return gcopy(x);
1438: n=lg(x[1])-1; p1=ker(x); k=lg(p1)-1;
1439: if (k) { p1=suppl(p1); n=lg(p1)-1; }
1440: else p1=idmat(n);
1441:
1442: tetpil=avma; p2=cgetg(n-k+1,t_MAT);
1443: for (i=k+1; i<=n; i++) p2[i-k]=lmul(x,(GEN) p1[i]);
1444: return gerepile(av,tetpil,p2);
1445: }
1446:
1447: GEN
1448: matimage0(GEN x,long flag)
1449: {
1450: switch(flag)
1451: {
1452: case 0: return image(x);
1453: case 1: return image2(x);
1454: default: err(flagerr,"matimage");
1455: }
1456: return NULL; /* not reached */
1457: }
1458:
1459: long
1460: rank(GEN x)
1461: {
1462: long av = avma, r;
1463: GEN d;
1464:
1465: gauss_pivot(x,0,&d,&r);
1466: /* yield r = dim ker(x) */
1467:
1468: avma=av; if (d) free(d);
1469: return lg(x)-1 - r;
1470: }
1471:
1472: GEN
1473: indexrank(GEN x)
1474: {
1475: long av = avma, i,j,n,r;
1476: GEN res,d,p1,p2;
1477:
1478: /* yield r = dim ker(x) */
1479: gauss_pivot(x,0,&d,&r);
1480:
1481: /* now r = dim Im(x) */
1482: n = lg(x)-1; r = n - r;
1483:
1484: avma=av; res=cgetg(3,t_VEC);
1485: p1=cgetg(r+1,t_VEC); res[1]=(long)p1;
1486: p2=cgetg(r+1,t_VEC); res[2]=(long)p2;
1487: if (d)
1488: {
1489: for (i=0,j=1; j<=n; j++)
1490: if (d[j]) { i++; p1[i]=d[j]; p2[i]=j; }
1491: free(d);
1492: qsort(p1+1,r,sizeof(long),(QSCOMP)pari_compare_long);
1493: }
1494: for (i=1;i<=r;i++) { p1[i]=lstoi(p1[i]); p2[i]=lstoi(p2[i]); }
1495: return res;
1496: }
1497:
1498: /*******************************************************************/
1499: /* */
1500: /* LINEAR ALGEBRA MODULO P */
1501: /* */
1502: /*******************************************************************/
1503: #ifdef LONG_IS_64BIT
1504: # define MASK (0x7fffffff00000000UL)
1505: #else
1506: # define MASK (0x7fff0000UL)
1507: #endif
1508: static GEN
1509: ker_mod_p_small(GEN x, GEN pp, long nontriv)
1510: {
1511: GEN y,c,d;
1512: long a,i,j,k,r,t,n,m,av0,tetpil, p = pp[2], piv;
1513:
1514: n = lg(x)-1;
1515: m=lg(x[1])-1; r=0; av0 = avma;
1516: x = dummycopy(x);
1517: for (i=1; i<=n; i++)
1518: {
1519: GEN p1 = (GEN)x[i];
1520: for (j=1; j<=m; j++) p1[j] = itos((GEN)p1[j]);
1521: }
1522: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1523: d=new_chunk(n+1);
1524: for (k=1; k<=n; k++)
1525: {
1526: for (j=1; j<=m; j++)
1527: if (!c[j])
1528: {
1529: a = coeff(x,j,k) % p;
1530: if (a) break;
1531: }
1532: if (j>m)
1533: {
1534: if (nontriv) { avma=av0; return NULL; }
1535: r++; d[k]=0;
1536: }
1537: else
1538: {
1539: c[j]=k; d[k]=j;
1540: {
1541: long av1 = avma;
1542: GEN p1 = mpinvmod(stoi(a), pp);
1543: piv = -itos(p1); avma = av1;
1544: }
1545: coeff(x,j,k) = -1;
1546: for (i=k+1; i<=n; i++)
1547: coeff(x,j,i) = (piv * coeff(x,j,i)) % p;
1548: for (t=1; t<=m; t++)
1549: if (t!=j)
1550: {
1551: piv = coeff(x,t,k) % p;
1552: if (piv)
1553: {
1554: coeff(x,t,k) = 0;
1555: for (i=k+1; i<=n; i++)
1556: {
1557: a = coeff(x,t,i) + piv * coeff(x,j,i);
1558: if (a & MASK) a %= p;
1559: coeff(x,t,i) = a;
1560: }
1561: }
1562: }
1563: }
1564: }
1565:
1566: tetpil=avma; y=cgetg(r+1,t_MAT);
1567: for (j=k=1; j<=r; j++,k++)
1568: {
1569: GEN c = cgetg(n+1,t_COL);
1570:
1571: y[j]=(long)c; while (d[k]) k++;
1572: for (i=1; i<k; i++)
1573: if (d[i])
1574: {
1575: long a = coeff(x,d[i],k) % p;
1576: if (a < 0) a += p;
1577: c[i] = lstoi(a);
1578: }
1579: else
1580: c[i] = zero;
1581: c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
1582: }
1583: return gerepile(av0,tetpil,y);
1584: }
1585:
1586: static GEN
1587: ker_mod_p_i(GEN x, GEN p, long nontriv)
1588: {
1589: GEN y,c,d,piv,mun;
1590: long i,j,k,r,t,n,m,av0,av,lim,tetpil;
1591:
1592: if (typ(x)!=t_MAT) err(typeer,"ker_mod_p");
1593: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
1594: if (!is_bigint(p) && p[2] < (MAXHALFULONG>>1))
1595: return ker_mod_p_small(x, p, nontriv);
1596:
1597: m=lg(x[1])-1; r=0; av0 = avma;
1598: x=dummycopy(x); mun=negi(gun);
1599: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1600: d=new_chunk(n+1);
1601: av=avma; lim=stack_lim(av,1);
1602: for (k=1; k<=n; k++)
1603: {
1604: for (j=1; j<=m; j++)
1605: if (!c[j])
1606: {
1607: coeff(x,j,k) = lmodii(gcoeff(x,j,k), p);
1608: if (signe(coeff(x,j,k))) break;
1609: }
1610: if (j>m)
1611: {
1612: if (nontriv) { avma = av0; return NULL; }
1613: r++; d[k]=0;
1614: for(j=1; j<k; j++)
1615: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
1616: }
1617: else
1618: {
1619: c[j]=k; d[k]=j; piv = negi(mpinvmod(gcoeff(x,j,k), p));
1620: coeff(x,j,k) = (long)mun;
1621: for (i=k+1; i<=n; i++)
1622: coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);
1623: for (t=1; t<=m; t++)
1624: if (t!=j)
1625: {
1626: piv = modii(gcoeff(x,t,k), p);
1627: if (signe(piv))
1628: {
1629: coeff(x,t,k)=zero;
1630: for (i=k+1; i<=n; i++)
1631: coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i)));
1632: if (low_stack(lim, stack_lim(av,1)))
1633: gerepile_gauss_keep_mod_p(x,p,m,n,k,t,av);
1634: }
1635: }
1636: }
1637: }
1638:
1639: tetpil=avma; y=cgetg(r+1,t_MAT);
1640: for (j=k=1; j<=r; j++,k++)
1641: {
1642: GEN c = cgetg(n+1,t_COL);
1643:
1644: y[j]=(long)c; while (d[k]) k++;
1645: for (i=1; i<k; i++)
1646: if (d[i])
1647: {
1648: GEN p1=gcoeff(x,d[i],k);
1649: c[i] = lmodii(p1, p); gunclone(p1);
1650: }
1651: else
1652: c[i] = zero;
1653: c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
1654: }
1655: return gerepile(av0,tetpil,y);
1656: }
1657:
1658: static void
1659: gauss_pivot_mod_p(GEN x, GEN p, GEN *dd, long *rr)
1660: {
1661: GEN c,d,piv;
1662: long i,j,k,r,t,n,m,av,lim;
1663:
1664: if (typ(x)!=t_MAT) err(typeer,"gauss_pivot_mod_p");
1665: n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; }
1666:
1667: m=lg(x[1])-1; r=0;
1668: x=dummycopy(x);
1669: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
1670: d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1);
1671: for (k=1; k<=n; k++)
1672: {
1673: for (j=1; j<=m; j++)
1674: if (!c[j])
1675: {
1676: coeff(x,j,k) = lmodii(gcoeff(x,j,k), p);
1677: if (signe(coeff(x,j,k))) break;
1678: }
1679: if (j>m) { r++; d[k]=0; }
1680: else
1681: {
1682: c[j]=k; d[k]=j; piv = negi(mpinvmod(gcoeff(x,j,k), p));
1683: for (i=k+1; i<=n; i++)
1684: coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);
1685: for (t=1; t<=m; t++)
1686: if (!c[t]) /* no pivot on that line yet */
1687: {
1688: piv=gcoeff(x,t,k);
1689: if (signe(piv))
1690: {
1691: coeff(x,t,k)=zero;
1692: for (i=k+1; i<=n; i++)
1693: coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i)));
1694: if (low_stack(lim, stack_lim(av,1)))
1695: gerepile_gauss(x,m,n,k,t,av,j,c);
1696: }
1697: }
1698: for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
1699: }
1700: }
1701: *dd=d; *rr=r;
1702: }
1703:
1704: GEN
1705: ker_mod_p(GEN x, GEN p)
1706: {
1707: return ker_mod_p_i(x, p, 0);
1708: }
1709:
1710: int
1711: ker_trivial_mod_p(GEN x, GEN p)
1712: {
1713: return ker_mod_p_i(x, p, 1)==NULL;
1714: }
1715:
1716: GEN
1717: image_mod_p(GEN x, GEN p)
1718: {
1719: GEN d,y;
1720: long j,k,r, av = avma;
1721:
1722: gauss_pivot_mod_p(x,p,&d,&r);
1723:
1724: /* r = dim ker(x) */
1725: if (!r) { avma=av; if (d) free(d); return gcopy(x); }
1726:
1727: /* r = dim Im(x) */
1728: r = lg(x)-1 - r; avma=av;
1729: y=cgetg(r+1,t_MAT);
1730: for (j=k=1; j<=r; k++)
1731: if (d[k]) y[j++] = lcopy((GEN)x[k]);
1732: free(d); return y;
1733: }
1734:
1735: /*******************************************************************/
1736: /* */
1737: /* EIGENVECTORS */
1738: /* (independent eigenvectors, sorted by increasing eigenvalue) */
1739: /* */
1740: /*******************************************************************/
1741:
1742: GEN
1743: eigen(GEN x, long prec)
1744: {
1745: GEN y,z,rr,p,ssesp,r1,r2,r3;
1746: long i,k,l,ly,av,tetpil,nbrac,ex, n = lg(x);
1747:
1748: if (typ(x)!=t_MAT) err(typeer,"eigen");
1749: if (n != 1 && n != lg(x[1])) err(mattype1,"eigen");
1750: if (n<=2) return gcopy(x);
1751:
1752: av=avma; ex = 16 - bit_accuracy(prec);
1753: y=cgetg(n,t_MAT); z=dummycopy(x);
1754: p=caradj(x,0,NULL); rr=roots(p,prec); nbrac=lg(rr)-1;
1755: for (i=1; i<=nbrac; i++)
1756: {
1757: GEN p1 = (GEN)rr[i];
1758: if (!signe(p1[2])) rr[i]=p1[1];
1759: }
1760: ly=1; k=1; r2=(GEN)rr[1];
1761: for(;;)
1762: {
1763: r3 = ground(r2); if (gexpo(gsub(r2,r3)) < ex) r2 = r3;
1764: for (i=1; i<n; i++)
1765: coeff(z,i,i) = lsub(gcoeff(x,i,i),r2);
1766: ssesp=ker0(z,prec); l=lg(ssesp);
1767: if (l == 1)
1768: err(talker, "precision too low in eigen");
1769: for (i=1; i<l; ) y[ly++]=ssesp[i++]; /* done with this eigenspace */
1770:
1771: r1=r2; /* try to find a different eigenvalue */
1772: do
1773: {
1774: if (k==nbrac)
1775: {
1776: tetpil=avma; setlg(y,ly); /* x may not be diagonalizable */
1777: return gerepile(av,tetpil,gcopy(y));
1778: }
1779: k++; r2=(GEN)rr[k];
1780: }
1781: while (gexpo(gsub(r1,r2)) < ex);
1782: }
1783: }
1784:
1785: /*******************************************************************/
1786: /* */
1787: /* DETERMINANT */
1788: /* */
1789: /*******************************************************************/
1790:
1791: GEN
1792: det0(GEN a,long flag)
1793: {
1794: switch(flag)
1795: {
1796: case 0: return det(a);
1797: case 1: return det2(a);
1798: default: err(flagerr,"matdet");
1799: }
1800: return NULL; /* not reached */
1801: }
1802:
1803: /* Exact types: choose the first non-zero pivot. Otherwise: maximal pivot */
1804: static GEN
1805: det_simple_gauss(GEN a, long inexact)
1806: {
1807: long i,j,k,av,av1,s, nbco = lg(a)-1;
1808: GEN x,p;
1809:
1810: av=avma; s=1; x=gun; a=dummycopy(a);
1811: for (i=1; i<nbco; i++)
1812: {
1813: p=gcoeff(a,i,i); k=i;
1814: if (inexact)
1815: {
1816: long e, ex = gexpo(p);
1817: GEN p1;
1818:
1819: for (j=i+1; j<=nbco; j++)
1820: {
1821: e = gexpo(gcoeff(a,i,j));
1822: if (e > ex) { ex=e; k=j; }
1823: }
1824: p1 = gcoeff(a,i,k);
1825: if (gcmp0(p1)) return gerepileupto(av, gcopy(p1));
1826: }
1827: else if (gcmp0(p))
1828: {
1829: do k++; while(k<=nbco && gcmp0(gcoeff(a,i,k)));
1830: if (k>nbco) return gerepileupto(av, gcopy(p));
1831: }
1832: if (k != i)
1833: {
1834: swap(a[i],a[k]); s = -s;
1835: p = gcoeff(a,i,i);
1836: }
1837:
1838: x = gmul(x,p);
1839: for (k=i+1; k<=nbco; k++)
1840: {
1841: GEN m = gcoeff(a,i,k);
1842: if (!gcmp0(m))
1843: {
1844: m = gneg_i(gdiv(m,p));
1845: for (j=i+1; j<=nbco; j++)
1846: coeff(a,j,k) = ladd(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
1847: }
1848: }
1849: }
1850: if (s<0) x = gneg_i(x);
1851: av1=avma; return gerepile(av,av1,gmul(x,gcoeff(a,nbco,nbco)));
1852: }
1853:
1854: /* a has integer entries, N = P^n */
1855: GEN
1856: det_mod_P_n(GEN a, GEN N, GEN P)
1857: {
1858: long va,i,j,k,s, av = avma, nbco = lg(a)-1;
1859: GEN x,p;
1860:
1861: s=1; va=0; x=gun; a=dummycopy(a);
1862: for (i=1; i<nbco; i++)
1863: {
1864: long fl = 0;
1865: for(;;)
1866: {
1867: for (k=i; k<=nbco; k++)
1868: {
1869: p=gcoeff(a,i,k);
1870: if (signe(p))
1871: {
1872: fl = 1;
1873: if (resii(p,P) != gzero) break;
1874: }
1875: }
1876: if (k <= nbco) break;
1877: va++; N = divii(N, P);
1878: if (!fl || is_pm1(N)) { avma=av; return gzero; }
1879: for (k=i; k<=nbco; k++) coeff(a,i,k) = ldivii(gcoeff(a,i,k), P);
1880: }
1881: if (k != i) { swap(a[i],a[k]); s = -s; }
1882:
1883: x = gmul(x,p); p = mpinvmod(p,N);
1884: for (k=i+1; k<=nbco; k++)
1885: {
1886: GEN m = resii(gcoeff(a,i,k), N);
1887: coeff(a,i,k) = zero;
1888: if (signe(m))
1889: {
1890: m = negi(resii(mulii(m,p), N));
1891: for (j=i+1; j<=nbco; j++)
1892: coeff(a,j,k) = lresii(addii(gcoeff(a,j,k),
1893: mulii(m,gcoeff(a,j,i))), N);
1894: }
1895: }
1896: }
1897: if (s<0) x = negi(x);
1898: x = resii(mulii(x,gcoeff(a,nbco,nbco)), N);
1899: return gerepileuptoint(av, mulii(x, gpowgs(P,va)));
1900: }
1901:
1902: GEN
1903: det2(GEN a)
1904: {
1905: long nbco = lg(a)-1;
1906: if (typ(a)!=t_MAT) err(mattype1,"det2");
1907: if (!nbco) return gun;
1908: if (nbco != lg(a[1])-1) err(mattype1,"det2");
1909: return det_simple_gauss(a,use_maximal_pivot(a));
1910: }
1911:
1912: /* determinant in a ring A: all computations are done within A
1913: * (Gauss-Bareiss algorithm)
1914: */
1915: GEN
1916: det(GEN a)
1917: {
1918: long nbco = lg(a)-1,i,j,k,av,s;
1919: GEN p,pprec;
1920:
1921: if (typ(a)!=t_MAT) err(mattype1,"det");
1922: if (!nbco) return gun;
1923: if (nbco != lg(a[1])-1) err(mattype1,"det");
1924: if (use_maximal_pivot(a)) return det_simple_gauss(a,1);
1925:
1926: av=avma; a=dummycopy(a); s=1;
1927: if (DEBUGLEVEL > 7) timer2();
1928: for (pprec=gun,i=1; i<nbco; i++,pprec=p)
1929: {
1930: GEN *ci, *ck, m, p1;
1931: int diveuc = (gcmp1(pprec)==0);
1932:
1933: p = gcoeff(a,i,i);
1934: if (gcmp0(p))
1935: {
1936: k=i+1; while (k<=nbco && gcmp0(gcoeff(a,i,k))) k++;
1937: if (k>nbco) return gerepileupto(av, gcopy(p));
1938: swap(a[k], a[i]); s = -s;
1939: p=gcoeff(a,i,i);
1940: }
1941: ci = (GEN*)a[i];
1942: for (k=i+1; k<=nbco; k++)
1943: {
1944: ck = (GEN*)a[k]; m = (GEN)ck[i];
1945: if (gcmp0(m))
1946: {
1947: if (gcmp1(p))
1948: {
1949: if (!gcmp1(pprec))
1950: a[k] = (long)gdivexact((GEN)a[k], pprec);
1951: }
1952: else
1953: for (j=i+1; j<=nbco; j++)
1954: {
1955: p1 = gmul(p,ck[j]);
1956: if (diveuc) p1 = gdivexact(p1,pprec);
1957: ck[j] = p1;
1958: }
1959: }
1960: else
1961: {
1962: m = gneg_i(m);
1963: for (j=i+1; j<=nbco; j++)
1964: {
1965: p1 = gadd(gmul(p,ck[j]), gmul(m,ci[j]));
1966: if (diveuc) p1 = gdivexact(p1,pprec);
1967: ck[j] = p1;
1968: }
1969: }
1970: }
1971: if (DEBUGLEVEL > 7) msgtimer("det, col %ld / %ld",i,nbco-1);
1972: }
1973: p = gcoeff(a,nbco,nbco);
1974: if (s < 0) p = gneg(p); else if (nbco==1) p = gcopy(p);
1975: return gerepileupto(av, p);
1976: }
1977:
1978: /*******************************************************************/
1979: /* */
1980: /* SPECIAL HNF (FOR INTERNAL USE !!!) */
1981: /* */
1982: /*******************************************************************/
1983: GEN lincomb_integral(GEN u, GEN v, GEN X, GEN Y);
1984: GEN vconcat(GEN Q1, GEN Q2);
1985:
1986: static int
1987: count(long **mat, long row, long len, long *firstnonzero)
1988: {
1989: int j, n=0;
1990:
1991: for (j=1; j<=len; j++)
1992: {
1993: long p = mat[j][row];
1994: if (p)
1995: {
1996: if (labs(p)!=1) return -1;
1997: n++; *firstnonzero=j;
1998: }
1999: }
2000: return n;
2001: }
2002:
2003: static int
2004: count2(long **mat, long row, long len)
2005: {
2006: int j;
2007: for (j=len; j; j--)
2008: if (labs(mat[j][row]) == 1) return j;
2009: return 0;
2010: }
2011:
2012: static GEN
2013: hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
2014: {
2015: GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
2016: GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
2017: long av,i,j,k,s,i1,j1,lim,zc;
2018: long co = lg(C);
2019: long col = lg(matgen)-1;
2020: long lnz, nlze, lig;
2021:
2022: if (col == 0) return matgen;
2023: lnz = lg(matgen[1])-1; /* was called lnz-1 - nr in hnfspec */
2024: nlze = lg(dep[1])-1;
2025: lig = nlze + lnz;
2026: if (DEBUGLEVEL>5)
2027: {
2028: fprintferr("Entering hnffinal:\n");
2029: if (DEBUGLEVEL>6)
2030: {
2031: if (nlze) fprintferr("dep = %Z\n",dep);
2032: fprintferr("mit = %Z\n",matgen);
2033: fprintferr("B = %Z\n",B);
2034: }
2035: }
2036: /* [LLLKERIM]
2037: u1u2=lllkerim(matgen); u1=(GEN)u1u2[1]; u2=(GEN)u1u2[2];
2038: if (DEBUGLEVEL>6) fprintferr("lllkerim done\n");
2039: if (lg(u2)<=lnz)
2040: err(talker,"matrix not of maximal rank in hermite spec");
2041: p1=gmul(matgen,u2);
2042: detmat=absi(det(p1));
2043: if (DEBUGLEVEL>6) fprintferr("det done\n");
2044: H=hnfmod(p1,detmat);
2045: if (DEBUGLEVEL>6) fprintferr("hnfmod done\n");
2046: p2=gmul(u1,lllint(u1));
2047: if (DEBUGLEVEL>6) fprintferr("lllint done\n");
2048: p3=gmul(u2,gauss(p1,H));
2049: if (DEBUGLEVEL>6) fprintferr("gauss done\n");
2050: U=cgetg(col+1,t_MAT);
2051: for (j=1; j<lg(p2); j++) U[j]=p2[j];
2052: for (j=lg(p2); j<=col; j++) U[j]=p3[j+1-lg(p2)]; */
2053:
2054: /* [HNFHAVAS]
2055:
2056: p2=hnfhavas(matgen); p1=(GEN)p2[1]; U=(GEN)p2[2]; p5=(GEN)p2[3];
2057: if (DEBUGLEVEL>6) fprintferr("hnfhavas done\n");
2058: for (i=1; i < lg(p1) && gcmp0(p1[i]); i++);
2059: i1=i-1;
2060: u1=cgetg(i,t_MAT); for (j=1; j<i; j++) u1[j]=U[j];
2061: H=cgetg(j1=lg(p1)-i1,t_MAT); for (j=1; j<j1; j++) H[j]=p1[i1+j];
2062: p2=cgetg(lg(p5),t_VEC);
2063: for (i=1; i<lg(p5); i++) p2[i]=lstoi(perm[nlze+itos(p5[i])]);
2064: for (i=1; i<lg(p5); i++) perm[nlze+i]=itos(p2[i]);
2065: p2=u1;
2066: p1=cgetg(j1,t_MAT); for (j=1; j<j1; j++) p1[j]=U[i1+j];
2067: Bnew=cgetg(co-col,t_MAT);
2068: for (j=1; j<co-col; j++)
2069: {
2070: p3=cgetg(lig+1,t_COL); Bnew[j]=(long)p3;
2071: for (i=1; i<=nlze; i++) p3[i]=coeff(B,i,j);
2072: for (; i<=lig; i++) p3[i]=coeff(B,nlze+itos(p5[i-nlze]),j);
2073: }
2074: B=Bnew; */
2075:
2076: /* [HNFBATUT] */
2077: p1 = hnfall(matgen);
2078: H = (GEN)p1[1]; /* lnz x lnz */
2079: U = (GEN)p1[2]; /* col x col */
2080: /* Only keep the part above the H (above the 0s is 0 since the dep rows
2081: * are dependant from the ones in matgen) */
2082: zc = col - lnz; /* # of 0 columns, correspond to units */
2083: if (nlze) { dep = gmul(dep,U); dep += zc; }
2084:
2085: diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
2086:
2087: av = avma; lim = stack_lim(av,1);
2088: Cnew = cgetg(co,t_MAT);
2089: setlg(C, col+1);
2090: p1 = gmul(C,U); setlg(C, co);
2091: for (j=1; j<=col; j++) Cnew[j] = p1[j];
2092: for ( ; j<co ; j++) Cnew[j] = C[j];
2093: if (DEBUGLEVEL>5) fprintferr(" hnfall done\n");
2094:
2095: /* Clean up B using new H */
2096: for (s=0,i=lnz; i; i--)
2097: {
2098: GEN h = gcoeff(H,i,i);
2099: if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; }
2100: for (j=col+1; j<co; j++)
2101: {
2102: GEN z = (GEN)B[j-col];
2103: p1 = (GEN)z[i+nlze]; if (h) p1 = gdivent(p1,h);
2104: for (k=1; k<=nlze; k++)
2105: z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(dep,k,i)));
2106: for ( ; k<=lig; k++)
2107: z[k] = lsubii((GEN)z[k], mulii(p1, gcoeff(H,k-nlze,i)));
2108: Cnew[j] = lsub((GEN)Cnew[j], gmul(p1, (GEN)Cnew[i+zc]));
2109: }
2110: if (low_stack(lim, stack_lim(av,1)))
2111: {
2112: GEN *gptr[2]; gptr[0]=&Cnew; gptr[1]=&B;
2113: if(DEBUGMEM>1) err(warnmem,"hnffinal, i = %ld",i);
2114: gerepilemany(av,gptr,2);
2115: }
2116: }
2117: p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
2118: for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
2119: if (diagH1[i])
2120: p1[++j1] = p2[i];
2121: else
2122: p2[++i1] = p2[i];
2123: for (i=i1+1; i<=lnz; i++) p2[i] = p1[i];
2124: if (DEBUGLEVEL>5) fprintferr(" first pass in hnffinal done\n");
2125:
2126: /* s = # extra redundant generators taken from H
2127: * zc col-s co zc = col lnz
2128: * [ 0 |dep | ] i = lnze + lnz - s = lig - s
2129: * nlze [--------| B' ]
2130: * [ 0 | H' | ] H' = H minus the s rows with a 1 on diagonal
2131: * i [--------|-----] lig-s (= "1-rows")
2132: * [ 0 | Id ]
2133: * [ | ] li */
2134: lig -= s; col -= s; lnz -= s;
2135: Hnew = cgetg(lnz+1,t_MAT);
2136: if (nlze) depnew = cgetg(lnz+1,t_MAT);
2137: Bnew = cgetg(co-col,t_MAT);
2138: C = dummycopy(Cnew);
2139: for (j=1,i1=j1=0; j<=lnz+s; j++)
2140: {
2141: GEN z = (GEN)H[j];
2142: if (diagH1[j])
2143: { /* hit exactly s times */
2144: i1++; p1 = cgetg(lig+1,t_COL); Bnew[i1] = (long)p1;
2145: C[i1+col] = Cnew[j+zc];
2146: for (i=1; i<=nlze; i++) p1[i] = coeff(dep,i,j);
2147: p1 += nlze;
2148: }
2149: else
2150: {
2151: j1++; p1 = cgetg(lnz+1,t_COL); Hnew[j1] = (long)p1;
2152: C[j1+zc] = Cnew[j+zc];
2153: if (nlze) depnew[j1] = dep[j];
2154: }
2155: for (i=k=1; k<=lnz; i++)
2156: if (!diagH1[i]) p1[k++] = z[i];
2157: }
2158: for (j=s+1; j<co-col; j++)
2159: {
2160: GEN z = (GEN)B[j-s];
2161: p1 = cgetg(lig+1,t_COL); Bnew[j] = (long)p1;
2162: for (i=1; i<=nlze; i++) p1[i] = z[i];
2163: z += nlze; p1 += nlze;
2164: for (i=k=1; k<=lnz; i++)
2165: if (!diagH1[i]) p1[k++] = z[i];
2166: }
2167: if (DEBUGLEVEL>5)
2168: {
2169: fprintferr("Leaving hnffinal\n");
2170: if (DEBUGLEVEL>6)
2171: {
2172: if (nlze) fprintferr("dep = %Z\n",depnew);
2173: fprintferr("mit = %Z\n",Hnew); outerr(Hnew);
2174: fprintferr("B = %Z\n",Bnew);
2175: fprintferr("C = %Z\n",C);
2176: }
2177: }
2178: if (nlze) *ptdep = depnew;
2179: *ptC = C;
2180: *ptB = Bnew; return Hnew;
2181: }
2182:
2183: /* for debugging */
2184: static void
2185: p_mat(long **mat, long *perm, long k0)
2186: {
2187: long av=avma, i,j;
2188: GEN p1, matj, matgen;
2189: long co = lg(mat);
2190: long li = lg(perm);
2191:
2192: fprintferr("Permutation: %Z\n",perm);
2193: matgen = cgetg(co,t_MAT);
2194: for (j=1; j<co; j++)
2195: {
2196: p1 = cgetg(li-k0,t_COL); matgen[j]=(long)p1;
2197: p1 -= k0; matj = mat[j];
2198: for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);
2199: }
2200: if (DEBUGLEVEL > 6) fprintferr("matgen = %Z\n",matgen);
2201: avma=av;
2202: }
2203:
2204: #define gswap(x,y) { long *_t=x; x=y; y=_t; }
2205:
2206: /* HNF reduce a relation matrix (column operations + row permutation)
2207: ** Input:
2208: ** mat = (li-1) x (co-1) matrix of long
2209: ** C = r x (co-1) matrix of GEN
2210: ** perm= permutation vector (length li-1), indexing the rows of mat: easier
2211: ** to maintain perm than to copy rows. For columns we can do it directly
2212: ** using e.g. swap(mat[i], mat[j])
2213: ** k0 = integer. The k0 first lines of mat are dense, the others are sparse.
2214: ** Output: cf ASCII art in the function body
2215: **
2216: ** row permutations applied to perm
2217: ** column operations applied to C.
2218: **/
2219: GEN
2220: hnfspec(long** mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
2221: {
2222: long av=avma,av2,*p,i,j,k,lk0,col,lig,*matj;
2223: long n,s,t,lim,nlze,lnz,nr;
2224: GEN p1,p2,matb,matbnew,vmax,matt,T,extramat;
2225: GEN B,H,dep,permpro;
2226: GEN *gptr[4];
2227: long co = lg(mat);
2228: long li = lg(perm); /* = lg(mat[1]) */
2229: int updateT = 1;
2230:
2231: if (DEBUGLEVEL>5)
2232: {
2233: fprintferr("Entering hnfspec\n");
2234: p_mat(mat,perm,0);
2235: }
2236: matt = cgetg(co,t_MAT); /* dense part of mat (top) */
2237: for (j=1; j<co; j++)
2238: {
2239: p1=cgetg(k0+1,t_COL); matt[j]=(long)p1; matj = mat[j];
2240: for (i=1; i<=k0; i++) p1[i] = lstoi(matj[perm[i]]);
2241: }
2242: vmax = cgetg(co,t_VECSMALL);
2243: av2 = avma; lim = stack_lim(av2,1);
2244:
2245: i=lig=li-1; col=co-1; lk0=k0;
2246: if (k0 || (lg(*ptC) > 1 && lg((*ptC)[1]) > 1)) T = idmat(col);
2247: else
2248: { /* dummy ! */
2249: GEN z = cgetg(1,t_COL);
2250: T = cgetg(co, t_MAT); updateT = 0;
2251: for (j=1; j<co; j++) T[j] = (long)z;
2252: }
2253: /* Look for lines with a single non0 entry, equal to ±1 */
2254: while (i > lk0)
2255: switch( count(mat,perm[i],col,&n) )
2256: {
2257: case 0: /* move zero lines between k0+1 and lk0 */
2258: lk0++; swap(perm[i], perm[lk0]);
2259: i=lig; continue;
2260:
2261: case 1: /* move trivial generator between lig+1 and li */
2262: swap(perm[i], perm[lig]);
2263: swap(T[n], T[col]);
2264: gswap(mat[n], mat[col]); p = mat[col];
2265: if (p[perm[lig]] < 0) /* = -1 */
2266: { /* convert relation -g = 0 to g = 0 */
2267: for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
2268: if (updateT)
2269: {
2270: p1 = (GEN)T[col];
2271: for (i=1; ; i++)
2272: if (signe((GEN)p1[i])) { p1[i] = lnegi((GEN)p1[i]); break; }
2273: }
2274: }
2275: lig--; col--; i=lig; continue;
2276:
2277: default: i--;
2278: }
2279: if (DEBUGLEVEL>5)
2280: {
2281: fprintferr(" after phase1:\n");
2282: p_mat(mat,perm,0);
2283: }
2284:
2285: #define absmax(s,z) {long _z = labs(z); if (_z > s) s = _z;}
2286:
2287: #if 0 /* TODO: check, and put back in */
2288: /* Get rid of all lines containing only 0 and ± 1, keeping track of column
2289: * operations in T. Leave the rows 1..lk0 alone [up to k0, coeff
2290: * explosion, between k0+1 and lk0, row is 0]
2291: */
2292: s = 0;
2293: while (lig > lk0 && s < (HIGHBIT>>1))
2294: {
2295: for (i=lig; i>lk0; i--)
2296: if (count(mat,perm[i],col,&n) >= 0) break;
2297: if (i == lk0) break;
2298:
2299: /* only 0, ±1 entries, at least 2 of them non-zero */
2300: swap(perm[i], perm[lig]);
2301: swap(T[n], T[col]); p1 = (GEN)T[col];
2302: gswap(mat[n], mat[col]); p = mat[col];
2303: if (p[perm[lig]] < 0)
2304: {
2305: for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
2306: T[col] = lneg(p1);
2307: }
2308: for (j=1; j<n; j++)
2309: {
2310: matj = mat[j];
2311: if (! (t = matj[perm[lig]]) ) continue;
2312: if (t == 1)
2313: { /* t = 1 */
2314: for (i=lk0+1; i<=lig; i++)
2315: absmax(s, matj[perm[i]] -= p[perm[i]]);
2316: T[j] = lsub((GEN)T[j], p1);
2317: }
2318: else
2319: { /* t = -1 */
2320: for (i=lk0+1; i<=lig; i++)
2321: absmax(s, matj[perm[i]] += p[perm[i]]);
2322: T[j] = ladd((GEN)T[j], p1);
2323: }
2324: }
2325: lig--; col--;
2326: if (low_stack(lim, stack_lim(av2,1)))
2327: {
2328: if(DEBUGMEM>1) err(warnmem,"hnfspec[1]");
2329: T = gerepileupto(av2, gcopy(T));
2330: }
2331: }
2332: #endif
2333: /* As above with lines containing a ±1 (no other assumption).
2334: * Stop when single precision becomes dangerous */
2335: for (j=1; j<=col; j++)
2336: {
2337: matj = mat[j];
2338: for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
2339: vmax[j] = s;
2340: }
2341: while (lig > lk0)
2342: {
2343: for (i=lig; i>lk0; i--)
2344: if ( (n = count2(mat,perm[i],col)) ) break;
2345: if (i == lk0) break;
2346:
2347: swap(perm[i], perm[lig]);
2348: swap(vmax[n], vmax[col]);
2349: gswap(mat[n], mat[col]); p = mat[col];
2350: swap(T[n], T[col]); p1 = (GEN)T[col];
2351: if (p[perm[lig]] < 0)
2352: {
2353: for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
2354: p1 = gneg(p1); T[col] = (long)p1;
2355: }
2356: for (j=1; j<col; j++)
2357: {
2358: matj = mat[j];
2359: if (! (t = matj[perm[lig]]) ) continue;
2360: if (vmax[col] && labs(t) >= (HIGHBIT-vmax[j]) / vmax[col]) goto END2;
2361:
2362: for (s=0, i=lk0+1; i<=lig; i++)
2363: absmax(s, matj[perm[i]] -= t*p[perm[i]]);
2364: vmax[j] = s;
2365: T[j] = (long)lincomb_integral(gun,stoi(-t), (GEN)T[j],p1);
2366: }
2367: lig--; col--;
2368: if (low_stack(lim, stack_lim(av2,1)))
2369: {
2370: if(DEBUGMEM>1) err(warnmem,"hnfspec[2]");
2371: T = gerepileupto(av2,gcopy(T));
2372: }
2373: }
2374:
2375: END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
2376: /* go multiprecision first */
2377: matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
2378: for (j=1; j<co; j++)
2379: {
2380: p1 = cgetg(li-k0,t_COL); matb[j] = (long)p1;
2381: p1 -= k0; matj = mat[j];
2382: for (i=k0+1; i<li; i++) p1[i] = lstoi(matj[perm[i]]);
2383: }
2384: if (DEBUGLEVEL>5)
2385: {
2386: fprintferr(" after phase2:\n");
2387: p_mat(mat,perm,k0);
2388: }
2389: for (i=li-2; i>lig; i--)
2390: {
2391: long i1, i0 = i - k0, k = i + co-li;
2392: GEN Bk = (GEN)matb[k];
2393: GEN Tk = (GEN)T[k];
2394: for (j=k+1; j<co; j++)
2395: {
2396: p1=(GEN)matb[j]; p2=(GEN)p1[i0];
2397: if (! (s=signe(p2)) ) continue;
2398:
2399: p1[i0] = zero;
2400: if (is_pm1(p2))
2401: {
2402: if (s > 0)
2403: { /* p2 = 1 */
2404: for (i1=1; i1<i0; i1++)
2405: p1[i1] = lsubii((GEN)p1[i1], (GEN)Bk[i1]);
2406: T[j] = lsub((GEN)T[j], Tk);
2407: }
2408: else
2409: { /* p2 = -1 */
2410: for (i1=1; i1<i0; i1++)
2411: p1[i1] = laddii((GEN)p1[i1], (GEN)Bk[i1]);
2412: T[j] = ladd((GEN)T[j], Tk);
2413: }
2414: }
2415: else
2416: {
2417: for (i1=1; i1<i0; i1++)
2418: p1[i1] = lsubii((GEN)p1[i1], mulii(p2,(GEN) Bk[i1]));
2419: T[j] = (long)lincomb_integral(gun,negi(p2), (GEN)T[j],Tk);
2420: }
2421: }
2422: if (low_stack(lim, stack_lim(av2,1)))
2423: {
2424: if(DEBUGMEM>1) err(warnmem,"hnfspec[3], i = %ld", i);
2425: for (j=1; j<co; j++) setlg(matb[j], i0+1); /* bottom can be forgotten */
2426: gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
2427: }
2428: }
2429: gptr[0]=&T; gptr[1]=&matb; gerepilemany(av2,gptr,2);
2430: if (DEBUGLEVEL>5)
2431: {
2432: fprintferr(" matb cleaned up (using Id block)\n");
2433: if (DEBUGLEVEL>6) outerr(matb);
2434: }
2435:
2436: nlze = lk0 - k0; /* # of 0 rows */
2437: lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
2438: if (updateT) matt = gmul(matt,T); /* update top rows */
2439: extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
2440: for (j=1; j<=col; j++)
2441: {
2442: GEN z = (GEN)matt[j];
2443: GEN t = ((GEN)matb[j]) + nlze - k0;
2444: p2=cgetg(lnz,t_COL); extramat[j]=(long)p2;
2445: for (i=1; i<=k0; i++) p2[i] = z[i]; /* top k0 rows */
2446: for ( ; i<lnz; i++) p2[i] = t[i]; /* other non-0 rows */
2447: }
2448: permpro = imagecomplspec(extramat, &nr); /* lnz = lg(permpro) */
2449:
2450: if (nlze)
2451: { /* put the nlze 0 rows (trivial generators) at the top */
2452: p1 = new_chunk(lk0+1);
2453: for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
2454: for ( ; i<=lk0; i++) p1[i] = perm[i - nlze];
2455: for (i=1; i<=lk0; i++) perm[i] = p1[i];
2456: }
2457: /* sort other rows according to permpro (nr redundant generators first) */
2458: p1 = new_chunk(lnz); p2 = perm + nlze;
2459: for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
2460: for (i=1; i<lnz; i++) p2[i] = p1[i];
2461: /* perm indexes the rows of mat
2462: * |_0__|__redund__|__dense__|__too big__|_____done______|
2463: * 0 nlze lig li
2464: * \___nr___/ \___k0__/
2465: * \____________lnz ______________/
2466: *
2467: * col co
2468: * [dep | ]
2469: * i0 [--------| B ] (i0 = nlze + nr)
2470: * [matbnew | ] matbnew has maximal rank = lnz-1 - nr
2471: * mat = [--------|-----] lig
2472: * [ 0 | Id ]
2473: * [ | ] li */
2474:
2475: matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
2476: dep = cgetg(col+1,t_MAT); /* rows dependant from the ones in matbnew */
2477: for (j=1; j<=col; j++)
2478: {
2479: GEN z = (GEN)extramat[j];
2480: p1 = cgetg(nlze+nr+1,t_COL); dep[j]=(long)p1;
2481: p2 = cgetg(lnz-nr,t_COL); matbnew[j]=(long)p2;
2482: for (i=1; i<=nlze; i++) p1[i]=zero;
2483: p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
2484: p2 -= nr; for ( ; i<lnz; i++) p2[i] = z[permpro[i]];
2485: }
2486:
2487: /* redundant generators in terms of the genuine generators
2488: * (x_i) = - (g_i) B */
2489: B = cgetg(co-col,t_MAT);
2490: for (j=col+1; j<co; j++)
2491: {
2492: GEN y = (GEN)matt[j];
2493: GEN z = (GEN)matb[j];
2494: p1=cgetg(lig+1,t_COL); B[j-col]=(long)p1;
2495: for (i=1; i<=nlze; i++) p1[i] = z[i];
2496: p1 += nlze; z += nlze-k0;
2497: for (k=1; k<lnz; k++)
2498: {
2499: i = permpro[k];
2500: p1[k] = (i <= k0)? y[i]: z[i];
2501: }
2502: }
2503: if (updateT) *ptC = gmul(*ptC,T);
2504: *ptdep = dep;
2505: *ptB = B;
2506: H = hnffinal(matbnew,perm,ptdep,ptB,ptC);
2507: gptr[0]=ptC;
2508: gptr[1]=ptdep;
2509: gptr[2]=ptB;
2510: gptr[3]=&H; gerepilemany(av,gptr,4);
2511: if (DEBUGLEVEL)
2512: msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1);
2513: return H;
2514: }
2515:
2516: /* HNF reduce x, apply same transforms to C */
2517: GEN
2518: mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
2519: {
2520: long i,j,ly,lx = lg(x);
2521: GEN p1,p2,z,perm;
2522: if (lx == 1) return gcopy(x);
2523: ly = lg(x[1]);
2524: z = cgetg(lx,t_MAT);
2525: perm = cgetg(ly,t_VECSMALL); *ptperm = perm;
2526: for (i=1; i<ly; i++) perm[i] = i;
2527: for (i=1; i<lx; i++)
2528: {
2529: p1 = cgetg(ly,t_COL); z[i] = (long)p1;
2530: p2 = (GEN)x[i];
2531: for (j=1; j<ly; j++) p1[j] = itos((GEN)p2[j]);
2532: }
2533: /* [ dep | ]
2534: * [-----| B ]
2535: * [ H | ]
2536: * [-----|-----]
2537: * [ 0 | Id ] */
2538: return hnfspec((long**)z,perm, ptdep, ptB, ptC, 0);
2539: }
2540:
2541: /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
2542: GEN
2543: hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
2544: GEN extramat,GEN extraC)
2545: {
2546: GEN p1,p2,p3,matb,extratop,Cnew,permpro;
2547: GEN B=*ptB, C=*ptC, dep=*ptdep, *gptr[4];
2548: long av = avma, i,j,lextra,colnew;
2549: long li = lg(perm);
2550: long co = lg(C);
2551: long lB = lg(B);
2552: long lig = li - lB;
2553: long col = co - lB;
2554: long lH = lg(H)-1;
2555: long nlze = lH? lg(dep[1])-1: lg(B[1])-1;
2556:
2557: if (DEBUGLEVEL>5)
2558: {
2559: fprintferr("Entering hnfadd:\n");
2560: if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);
2561: }
2562: /* col co
2563: * [ 0 |dep | ]
2564: * nlze [--------| B ]
2565: * [ 0 | H | ]
2566: * [--------|-----] lig
2567: * [ 0 | Id ]
2568: * [ | ] li */
2569: lextra = lg(extramat)-1;
2570: extratop = cgetg(lextra+1,t_MAT); /* [1..lig] part (top) */
2571: p2 = cgetg(lextra+1,t_MAT); /* bottom */
2572: for (j=1; j<=lextra; j++)
2573: {
2574: GEN z = (GEN)extramat[j];
2575: p1=cgetg(lig+1,t_COL); extratop[j] = (long)p1;
2576: for (i=1; i<=lig; i++) p1[i] = z[i];
2577: p1=cgetg(lB,t_COL); p2[j] = (long)p1;
2578: p1 -= lig;
2579: for ( ; i<li; i++) p1[i] = z[i];
2580: }
2581: if (li-1 != lig)
2582: { /* zero out bottom part, using the Id block */
2583: GEN A = cgetg(lB,t_MAT);
2584: for (j=1; j<lB; j++) A[j] = C[j+col];
2585: extraC = gsub(extraC, gmul(A,p2));
2586: extratop = gsub(extratop,gmul(B,p2));
2587: }
2588:
2589: colnew = lH + lextra;
2590: extramat = cgetg(colnew+1,t_MAT);
2591: Cnew = cgetg(lB+colnew,t_MAT);
2592: for (j=1; j<=lextra; j++)
2593: {
2594: extramat[j] = extratop[j];
2595: Cnew[j] = extraC[j];
2596: }
2597: for ( ; j<=colnew; j++)
2598: {
2599: p1 = cgetg(lig+1,t_COL); extramat[j] = (long)p1;
2600: p2 = (GEN)dep[j-lextra]; for (i=1; i<=nlze; i++) p1[i] = p2[i];
2601: p2 = (GEN) H[j-lextra]; for ( ; i<=lig ; i++) p1[i] = p2[i-nlze];
2602: }
2603: for (j=lextra+1; j<lB+colnew; j++)
2604: Cnew[j] = C[j-lextra+col-lH];
2605: if (DEBUGLEVEL>5)
2606: {
2607: fprintferr(" 1st phase done\n");
2608: if (DEBUGLEVEL>6) fprintferr("extramat = %Z\n",extramat);
2609: }
2610: permpro = imagecomplspec(extramat, &nlze);
2611: p1 = new_chunk(lig+1);
2612: for (i=1; i<=lig; i++) p1[i] = perm[permpro[i]];
2613: for (i=1; i<=lig; i++) perm[i] = p1[i];
2614:
2615: matb = cgetg(colnew+1,t_MAT);
2616: dep = cgetg(colnew+1,t_MAT);
2617: for (j=1; j<=colnew; j++)
2618: {
2619: GEN z = (GEN)extramat[j];
2620: p1=cgetg(nlze+1,t_COL); dep[j]=(long)p1;
2621: p2=cgetg(lig+1-nlze,t_COL); matb[j]=(long)p2;
2622: p2 -= nlze;
2623: for (i=1; i<=nlze; i++) p1[i] = z[permpro[i]];
2624: for ( ; i<= lig; i++) p2[i] = z[permpro[i]];
2625: }
2626: p3 = cgetg(lB,t_MAT);
2627: for (j=1; j<lB; j++)
2628: {
2629: p2 = (GEN)B[j];
2630: p1 = cgetg(lig+1,t_COL); p3[j] = (long)p1;
2631: for (i=1; i<=lig; i++) p1[i] = p2[permpro[i]];
2632: }
2633: B = p3;
2634: if (DEBUGLEVEL>5) fprintferr(" 2nd phase done\n");
2635: *ptdep = dep;
2636: *ptB = B;
2637: H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
2638: p1 = cgetg(co+lextra,t_MAT);
2639: for (j=1; j <= col-lH; j++) p1[j] = C[j];
2640: C = Cnew - (col-lH);
2641: for ( ; j < co+lextra; j++) p1[j] = C[j];
2642:
2643: gptr[0]=ptC; *ptC=p1;
2644: gptr[1]=ptdep;
2645: gptr[2]=ptB;
2646: gptr[3]=&H; gerepilemany(av,gptr,4);
2647: if (DEBUGLEVEL)
2648: {
2649: if (DEBUGLEVEL>7)
2650: {
2651: fprintferr("mit = %Z\n",H);
2652: fprintferr("C = %Z\n",p1);
2653: }
2654: msgtimer("hnfadd (%d)",lextra);
2655: }
2656: return H;
2657: }
2658:
2659: /* return a solution of congruence system sum M_{ i,j } X_j = Y_i mod D_i
2660: * If ptu1 != NULL, put in *ptu1 a Z-basis of the homogeneous system
2661: */
2662: static GEN
2663: gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
2664: {
2665: long n,m,i,j,lM,av=avma,tetpil;
2666: GEN p1,delta,H,U,u1,u2,x;
2667:
2668: if (typ(M)!=t_MAT) err(typeer,"gaussmodulo");
2669: lM = lg(M); m = lM-1;
2670: if (!m)
2671: {
2672: if ((typ(Y)!=t_INT && lg(Y)!=1)
2673: || (typ(D)!=t_INT && lg(D)!=1)) err(consister,"gaussmodulo");
2674: return gzero;
2675: }
2676: n = lg(M[1])-1;
2677: switch(typ(D))
2678: {
2679: case t_VEC:
2680: case t_COL: delta=diagonal(D); break;
2681: case t_INT: delta=gscalmat(D,n); break;
2682: default: err(typeer,"gaussmodulo");
2683: }
2684: if (typ(Y) == t_INT)
2685: {
2686: p1 = cgetg(n+1,t_COL);
2687: for (i=1; i<=n; i++) p1[i]=(long)Y;
2688: Y = p1;
2689: }
2690: p1 = hnfall(concatsp(M,delta));
2691: H = (GEN)p1[1]; U = (GEN)p1[2];
2692: Y = gauss(H,Y);
2693: if (!gcmp1(denom(Y))) return gzero;
2694: u1 = cgetg(m+1,t_MAT);
2695: u2 = cgetg(n+1,t_MAT);
2696: for (j=1; j<=m; j++)
2697: {
2698: p1 = (GEN)U[j]; setlg(p1,lM);
2699: u1[j] = (long)p1;
2700: }
2701: U += m;
2702: for (j=1; j<=n; j++)
2703: {
2704: p1 = (GEN)U[j]; setlg(p1,lM);
2705: u2[j] = (long)p1;
2706: }
2707: x = gmul(u2,Y);
2708: tetpil=avma; x=lllreducemodmatrix(x,u1);
2709: if (!ptu1) x = gerepile(av,tetpil,x);
2710: else
2711: {
2712: GEN *gptr[2];
2713: *ptu1=gcopy(u1); gptr[0]=ptu1; gptr[1]=&x;
2714: gerepilemanysp(av,tetpil,gptr,2);
2715: }
2716: return x;
2717: }
2718:
2719: GEN
2720: matsolvemod0(GEN M, GEN D, GEN Y, long flag)
2721: {
2722: long av;
2723: GEN p1,y;
2724:
2725: if (!flag) return gaussmoduloall(M,D,Y,NULL);
2726:
2727: av=avma; y = cgetg(3,t_VEC);
2728: p1 = gaussmoduloall(M,D,Y, (GEN*)y+2);
2729: if (p1==gzero) { avma=av; return gzero; }
2730: y[1] = (long)p1; return y;
2731: }
2732:
2733: GEN
2734: gaussmodulo2(GEN M, GEN D, GEN Y)
2735: {
2736: return matsolvemod0(M,D,Y,1);
2737: }
2738:
2739: GEN
2740: gaussmodulo(GEN M, GEN D, GEN Y)
2741: {
2742: return matsolvemod0(M,D,Y,0);
2743: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to alglin1.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