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