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