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