Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** LINEAR ALGEBRA **/
18 : /** (second part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mat
25 :
26 : /*******************************************************************/
27 : /* */
28 : /* CHARACTERISTIC POLYNOMIAL */
29 : /* */
30 : /*******************************************************************/
31 :
32 : static GEN
33 23351 : Flm_charpoly_i(GEN x, ulong p)
34 : {
35 23351 : long lx = lg(x), r, i;
36 23351 : GEN H, y = cgetg(lx+1, t_VEC);
37 23351 : gel(y,1) = pol1_Flx(0); H = Flm_hess(x, p);
38 150694 : for (r = 1; r < lx; r++)
39 : {
40 127342 : pari_sp av2 = avma;
41 127342 : ulong a = 1;
42 127342 : GEN z, b = zero_Flx(0);
43 378126 : for (i = r-1; i; i--)
44 : {
45 302073 : a = Fl_mul(a, ucoeff(H,i+1,i), p);
46 302069 : if (!a) break;
47 250812 : b = Flx_add(b, Flx_Fl_mul(gel(y,i), Fl_mul(a,ucoeff(H,i,r),p), p), p);
48 : }
49 127340 : z = Flx_sub(Flx_shift(gel(y,r), 1),
50 127310 : Flx_Fl_mul(gel(y,r), ucoeff(H,r,r), p), p);
51 : /* (X - H[r,r])y[r] - b */
52 127343 : gel(y,r+1) = gerepileuptoleaf(av2, Flx_sub(z, b, p));
53 : }
54 23352 : return gel(y,lx);
55 : }
56 :
57 : GEN
58 1733 : Flm_charpoly(GEN x, ulong p)
59 : {
60 1733 : pari_sp av = avma;
61 1733 : return gerepileuptoleaf(av, Flm_charpoly_i(x,p));
62 : }
63 :
64 : GEN
65 19131 : FpM_charpoly(GEN x, GEN p)
66 : {
67 19131 : pari_sp av = avma;
68 : long lx, r, i;
69 : GEN y, H;
70 :
71 19131 : if (lgefint(p) == 3)
72 : {
73 18395 : ulong pp = p[2];
74 18395 : y = Flx_to_ZX(Flm_charpoly_i(ZM_to_Flm(x,pp), pp));
75 18395 : return gerepileupto(av, y);
76 : }
77 736 : lx = lg(x); y = cgetg(lx+1, t_VEC);
78 736 : gel(y,1) = pol_1(0); H = FpM_hess(x, p);
79 4048 : for (r = 1; r < lx; r++)
80 : {
81 4048 : pari_sp av2 = avma;
82 4048 : GEN z, a = gen_1, b = pol_0(0);
83 9229 : for (i = r-1; i; i--)
84 : {
85 7008 : a = Fp_mul(a, gcoeff(H,i+1,i), p);
86 7008 : if (!signe(a)) break;
87 5181 : b = ZX_add(b, ZX_Z_mul(gel(y,i), Fp_mul(a,gcoeff(H,i,r),p)));
88 : }
89 4048 : b = FpX_red(b, p);
90 4048 : z = FpX_sub(RgX_shift_shallow(gel(y,r), 1),
91 4048 : FpX_Fp_mul(gel(y,r), gcoeff(H,r,r), p), p);
92 4048 : z = FpX_sub(z,b,p);
93 4048 : if (r+1 == lx) { gel(y,lx) = z; break; }
94 3312 : gel(y,r+1) = gerepileupto(av2, z); /* (X - H[r,r])y[r] - b */
95 : }
96 736 : return gerepileupto(av, gel(y,lx));
97 : }
98 :
99 : GEN
100 420 : charpoly0(GEN x, long v, long flag)
101 : {
102 420 : if (v<0) v = 0;
103 420 : switch(flag)
104 : {
105 14 : case 0: return caradj(x,v,NULL);
106 14 : case 1: return caract(x,v);
107 14 : case 2: return carhess(x,v);
108 14 : case 3: return carberkowitz(x,v);
109 7 : case 4:
110 7 : if (typ(x) != t_MAT) pari_err_TYPE("charpoly",x);
111 7 : RgM_check_ZM(x, "charpoly");
112 7 : x = ZM_charpoly(x); setvarn(x, v); return x;
113 357 : case 5:
114 357 : return charpoly(x, v);
115 : }
116 0 : pari_err_FLAG("charpoly");
117 : return NULL; /* LCOV_EXCL_LINE */
118 : }
119 :
120 : /* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
121 : static GEN
122 4794 : easychar(GEN x, long v)
123 : {
124 : pari_sp av;
125 : long lx;
126 : GEN p1;
127 :
128 4794 : switch(typ(x))
129 : {
130 35 : case t_INT: case t_REAL: case t_INTMOD:
131 : case t_FRAC: case t_PADIC:
132 35 : p1=cgetg(4,t_POL);
133 35 : p1[1]=evalsigne(1) | evalvarn(v);
134 35 : gel(p1,2) = gneg(x); gel(p1,3) = gen_1;
135 35 : return p1;
136 :
137 14 : case t_COMPLEX: case t_QUAD:
138 14 : p1 = cgetg(5,t_POL);
139 14 : p1[1] = evalsigne(1) | evalvarn(v);
140 14 : gel(p1,2) = gnorm(x); av = avma;
141 14 : gel(p1,3) = gerepileupto(av, gneg(gtrace(x)));
142 14 : gel(p1,4) = gen_1; return p1;
143 :
144 28 : case t_FFELT: {
145 28 : pari_sp ltop=avma;
146 28 : p1 = FpX_to_mod(FF_charpoly(x), FF_p_i(x));
147 28 : setvarn(p1,v); return gerepileupto(ltop,p1);
148 : }
149 :
150 70 : case t_POLMOD:
151 : {
152 70 : GEN A = gel(x,2), T = gel(x,1);
153 70 : if (typ(A)==t_POL && RgX_is_QX(A) && RgX_is_ZX(T))
154 28 : return QXQ_charpoly(A, T, v);
155 : else
156 42 : return RgXQ_charpoly(A, T, v);
157 : }
158 4647 : case t_MAT:
159 4647 : lx=lg(x);
160 4647 : if (lx==1) return pol_1(v);
161 4591 : if (lgcols(x) != lx) break;
162 4584 : return NULL;
163 : }
164 7 : pari_err_TYPE("easychar",x);
165 : return NULL; /* LCOV_EXCL_LINE */
166 : }
167 : /* compute charpoly by mapping to Fp first, return lift to Z */
168 : static GEN
169 35 : RgM_Fp_charpoly(GEN x, GEN p, long v)
170 : {
171 : GEN T;
172 35 : if (lgefint(p) == 3)
173 : {
174 21 : ulong pp = itou(p);
175 21 : T = Flm_charpoly_i(RgM_to_Flm(x, pp), pp);
176 21 : T = Flx_to_ZX(T);
177 : }
178 : else
179 14 : T = FpM_charpoly(RgM_to_FpM(x, p), p);
180 35 : setvarn(T, v); return T;
181 : }
182 : GEN
183 3009 : charpoly(GEN x, long v)
184 : {
185 3009 : GEN T, p = NULL;
186 : long prec;
187 3009 : if ((T = easychar(x,v))) return T;
188 2841 : switch(RgM_type(x, &p,&T,&prec))
189 : {
190 1910 : case t_INT:
191 1910 : T = ZM_charpoly(x); setvarn(T, v); break;
192 35 : case t_INTMOD:
193 35 : if (!BPSW_psp(p)) T = carberkowitz(x, v);
194 : else
195 : {
196 35 : pari_sp av = avma;
197 35 : T = RgM_Fp_charpoly(x,p,v);
198 35 : T = gerepileupto(av, FpX_to_mod(T,p));
199 : }
200 35 : break;
201 49 : case t_REAL:
202 : case t_COMPLEX:
203 : case t_PADIC:
204 49 : T = carhess(x, v);
205 49 : break;
206 847 : default:
207 847 : T = carberkowitz(x, v);
208 : }
209 2834 : return T;
210 : }
211 :
212 : /* We possibly worked with an "invalid" polynomial p, satisfying
213 : * varn(p) > gvar2(p). Fix this. */
214 : static GEN
215 1680 : fix_pol(pari_sp av, GEN p)
216 : {
217 1680 : long w = gvar2(p), v = varn(p);
218 1680 : if (w == v) pari_err_PRIORITY("charpoly", p, "=", w);
219 1673 : if (varncmp(w,v) < 0) p = gerepileupto(av, poleval(p, pol_x(v)));
220 1673 : return p;
221 : }
222 : GEN
223 14 : caract(GEN x, long v)
224 : {
225 14 : pari_sp av = avma;
226 : GEN T, C, x_k, Q;
227 : long k, n;
228 :
229 14 : if ((T = easychar(x,v))) return T;
230 :
231 14 : n = lg(x)-1;
232 14 : if (n == 1) return fix_pol(av, deg1pol(gen_1, gneg(gcoeff(x,1,1)), v));
233 :
234 14 : x_k = pol_x(v); /* to be modified in place */
235 14 : T = scalarpol(det(x), v); C = utoineg(n); Q = pol_x(v);
236 28 : for (k=1; k<=n; k++)
237 : {
238 28 : GEN mk = utoineg(k), d;
239 28 : gel(x_k,2) = mk;
240 28 : d = det(RgM_Rg_add_shallow(x, mk));
241 28 : T = RgX_add(RgX_mul(T, x_k), RgX_Rg_mul(Q, gmul(C, d)));
242 28 : if (k == n) break;
243 :
244 14 : Q = RgX_mul(Q, x_k);
245 14 : C = diviuexact(mulsi(k-n,C), k+1); /* (-1)^k binomial(n,k) */
246 : }
247 14 : return fix_pol(av, RgX_Rg_div(T, mpfact(n)));
248 : }
249 :
250 : /* C = charpoly(x, v) */
251 : static GEN
252 21 : RgM_adj_from_char(GEN x, long v, GEN C)
253 : {
254 21 : if (varn(C) != v) /* problem with variable priorities */
255 : {
256 7 : C = gdiv(gsub(C, gsubst(C, v, gen_0)), pol_x(v));
257 7 : if (odd(lg(x))) C = RgX_neg(C); /* even dimension */
258 7 : return gsubst(C, v, x);
259 : }
260 : else
261 : {
262 14 : C = RgX_shift_shallow(C, -1);
263 14 : if (odd(lg(x))) C = RgX_neg(C); /* even dimension */
264 14 : return RgX_RgM_eval(C, x);
265 : }
266 : }
267 : /* assume x square matrice */
268 : static GEN
269 1827 : mattrace(GEN x)
270 : {
271 1827 : long i, lx = lg(x);
272 : GEN t;
273 1827 : if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
274 1757 : t = gcoeff(x,1,1);
275 4151 : for (i = 2; i < lx; i++) t = gadd(t, gcoeff(x,i,i));
276 1757 : return t;
277 : }
278 : static int
279 56 : bad_char(GEN q, long n)
280 : {
281 : forprime_t S;
282 : ulong p;
283 56 : if (!signe(q)) return 0;
284 42 : (void)u_forprime_init(&S, 2, n);
285 98 : while ((p = u_forprime_next(&S)))
286 70 : if (!umodiu(q, p)) return 1;
287 28 : return 0;
288 : }
289 : /* Using traces: return the characteristic polynomial of x (in variable v).
290 : * If py != NULL, the adjoint matrix is put there. */
291 : GEN
292 119 : caradj(GEN x, long v, GEN *py)
293 : {
294 : pari_sp av, av0;
295 : long i, k, n;
296 : GEN T, y, t;
297 :
298 119 : if ((T = easychar(x, v)))
299 : {
300 42 : if (py)
301 : {
302 42 : if (typ(x) != t_MAT) pari_err_TYPE("matadjoint",x);
303 42 : *py = cgetg(1,t_MAT);
304 : }
305 42 : return T;
306 : }
307 :
308 77 : n = lg(x)-1; av0 = avma;
309 77 : T = cgetg(n+3,t_POL); T[1] = evalsigne(1) | evalvarn(v);
310 77 : gel(T,n+2) = gen_1;
311 77 : if (!n) { if (py) *py = cgetg(1,t_MAT); return T; }
312 77 : av = avma; t = gerepileupto(av, gneg(mattrace(x)));
313 77 : gel(T,n+1) = t;
314 77 : if (n == 1) {
315 7 : T = fix_pol(av0, T);
316 7 : if (py) *py = matid(1); return T;
317 : }
318 70 : if (n == 2) {
319 14 : GEN a = gcoeff(x,1,1), b = gcoeff(x,1,2);
320 14 : GEN c = gcoeff(x,2,1), d = gcoeff(x,2,2);
321 14 : av = avma;
322 14 : gel(T,2) = gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
323 14 : T = fix_pol(av0, T);
324 14 : if (py) {
325 7 : y = cgetg(3, t_MAT);
326 7 : gel(y,1) = mkcol2(gcopy(d), gneg(c));
327 7 : gel(y,2) = mkcol2(gneg(b), gcopy(a));
328 7 : *py = y;
329 : }
330 14 : return T;
331 : }
332 : /* l > 3 */
333 56 : if (bad_char(residual_characteristic(x), n))
334 : { /* n! not invertible in base ring */
335 14 : T = charpoly(x, v);
336 14 : if (!py) return gerepileupto(av, T);
337 14 : *py = RgM_adj_from_char(x, v, T); return gc_all(av, 2, &T,py);
338 : }
339 42 : av = avma; y = RgM_shallowcopy(x);
340 175 : for (i = 1; i <= n; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
341 91 : for (k = 2; k < n; k++)
342 : {
343 49 : GEN y0 = y;
344 49 : y = RgM_mul(y, x);
345 49 : t = gdivgs(mattrace(y), -k);
346 210 : for (i = 1; i <= n; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
347 49 : y = gclone(y);
348 49 : gel(T,n-k+2) = gerepilecopy(av, t); av = avma;
349 49 : if (k > 2) gunclone(y0);
350 : }
351 42 : t = gmul(gcoeff(x,1,1),gcoeff(y,1,1));
352 133 : for (i=2; i<=n; i++) t = gadd(t, gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
353 42 : gel(T,2) = gerepileupto(av, gneg(t));
354 42 : T = fix_pol(av0, T);
355 42 : if (py) *py = odd(n)? gcopy(y): RgM_neg(y);
356 42 : gunclone(y); return T;
357 : }
358 :
359 : GEN
360 105 : adj(GEN x)
361 : {
362 : GEN y;
363 105 : (void)caradj(x, fetch_var(), &y);
364 105 : (void)delete_var(); return y;
365 : }
366 :
367 : GEN
368 7 : adjsafe(GEN x)
369 : {
370 7 : const long v = fetch_var();
371 7 : pari_sp av = avma;
372 : GEN C, A;
373 7 : if (typ(x) != t_MAT) pari_err_TYPE("matadjoint",x);
374 7 : if (lg(x) < 3) return gcopy(x);
375 7 : C = charpoly(x,v);
376 7 : A = RgM_adj_from_char(x, v, C);
377 7 : (void)delete_var(); return gerepileupto(av, A);
378 : }
379 :
380 : GEN
381 112 : matadjoint0(GEN x, long flag)
382 : {
383 112 : switch(flag)
384 : {
385 105 : case 0: return adj(x);
386 7 : case 1: return adjsafe(x);
387 : }
388 0 : pari_err_FLAG("matadjoint");
389 : return NULL; /* LCOV_EXCL_LINE */
390 : }
391 :
392 : /*******************************************************************/
393 : /* */
394 : /* Frobenius form */
395 : /* */
396 : /*******************************************************************/
397 :
398 : /* The following section implement a mix of Ozello and Storjohann algorithms
399 :
400 : P. Ozello, doctoral thesis (in French):
401 : Calcul exact des formes de Jordan et de Frobenius d'une matrice, Chapitre 2
402 : http://tel.archives-ouvertes.fr/tel-00323705
403 :
404 : A. Storjohann, Diss. ETH No. 13922
405 : Algorithms for Matrix Canonical Forms, Chapter 9
406 : https://cs.uwaterloo.ca/~astorjoh/diss2up.pdf
407 :
408 : We use Storjohann Lemma 9.14 (step1, step2, step3) Ozello theorem 4,
409 : and Storjohann Lemma 9.18
410 : */
411 :
412 : /* Elementary transforms */
413 :
414 : /* M <- U^(-1) M U, U = E_{i,j}(k) => U^(-1) = E{i,j}(-k)
415 : * P = U * P */
416 : static void
417 14182 : transL(GEN M, GEN P, GEN k, long i, long j)
418 : {
419 14182 : long l, n = lg(M)-1;
420 173628 : for(l=1; l<=n; l++) /* M[,j]-=k*M[,i] */
421 159446 : gcoeff(M,l,j) = gsub(gcoeff(M,l,j), gmul(gcoeff(M,l,i), k));
422 173628 : for(l=1; l<=n; l++) /* M[i,]+=k*M[j,] */
423 159446 : gcoeff(M,i,l) = gadd(gcoeff(M,i,l), gmul(gcoeff(M,j,l), k));
424 14182 : if (P)
425 163653 : for(l=1; l<=n; l++)
426 150570 : gcoeff(P,i,l) = gadd(gcoeff(P,i,l), gmul(gcoeff(P,j,l), k));
427 14182 : }
428 :
429 : /* j = a or b */
430 : static void
431 2261 : transD(GEN M, GEN P, long a, long b, long j)
432 : {
433 : long l, n;
434 2261 : GEN k = gcoeff(M,a,b), ki;
435 :
436 2261 : if (gequal1(k)) return;
437 1127 : ki = ginv(k); n = lg(M)-1;
438 11802 : for(l=1; l<=n; l++)
439 10675 : if (l!=j)
440 : {
441 9548 : gcoeff(M,l,j) = gmul(gcoeff(M,l,j), k);
442 9548 : gcoeff(M,j,l) = (j==a && l==b)? gen_1: gmul(gcoeff(M,j,l), ki);
443 : }
444 1127 : if (P)
445 9457 : for(l=1; l<=n; l++)
446 8603 : gcoeff(P,j,l) = gmul(gcoeff(P,j,l), ki);
447 : }
448 :
449 : static void
450 378 : transS(GEN M, GEN P, long i, long j)
451 : {
452 378 : long l, n = lg(M)-1;
453 378 : swap(gel(M,i), gel(M,j));
454 5495 : for (l=1; l<=n; l++)
455 5117 : swap(gcoeff(M,i,l), gcoeff(M,j,l));
456 378 : if (P)
457 4060 : for (l=1; l<=n; l++)
458 3843 : swap(gcoeff(P,i,l), gcoeff(P,j,l));
459 378 : }
460 :
461 : /* Convert companion matrix to polynomial*/
462 : static GEN
463 280 : minpoly_polslice(GEN M, long i, long j, long v)
464 : {
465 280 : long k, d = j+1-i;
466 280 : GEN P = cgetg(d+3,t_POL);
467 280 : P[1] = evalsigne(1)|evalvarn(v);
468 1379 : for (k=0; k<d; k++)
469 1099 : gel(P,k+2) = gneg(gcoeff(M,i+k, j));
470 280 : gel(P,d+2) = gen_1;
471 280 : return P;
472 : }
473 :
474 : static GEN
475 49 : minpoly_listpolslice(GEN M, GEN V, long v)
476 : {
477 49 : long i, n = lg(M)-1, nb = lg(V)-1;
478 49 : GEN W = cgetg(nb+1, t_VEC);
479 147 : for (i=1; i<=nb; i++)
480 98 : gel(W,i) = minpoly_polslice(M, V[i], i < nb? V[i+1]-1: n, v);
481 49 : return W;
482 : }
483 :
484 : static int
485 91 : minpoly_dvdslice(GEN M, long i, long j, long k)
486 : {
487 91 : pari_sp av = avma;
488 91 : long r = signe(RgX_rem(minpoly_polslice(M, i, j-1, 0),
489 : minpoly_polslice(M, j, k, 0)));
490 91 : return gc_bool(av, r == 0);
491 : }
492 :
493 : static void
494 0 : RgM_replace(GEN M, GEN M2)
495 : {
496 0 : long n = lg(M)-1, m = nbrows(M), i, j;
497 0 : for(i=1; i<=n; i++)
498 0 : for(j=1; j<=m; j++)
499 0 : gcoeff(M, i, j) = gcoeff(M2, i, j);
500 0 : }
501 :
502 : static void
503 0 : gerepilemat2_inplace(pari_sp av, GEN M, GEN P)
504 : {
505 0 : GEN M2 = M, P2 = P;
506 0 : gerepileall(av, P ? 2: 1, &M2, &P2);
507 0 : RgM_replace(M, M2);
508 0 : if (P) RgM_replace(P, P2);
509 0 : }
510 :
511 : /* Lemma 9.14 */
512 : static long
513 581 : weakfrobenius_step1(GEN M, GEN P, long j0)
514 : {
515 581 : pari_sp av = avma;
516 581 : long n = lg(M)-1, k, j;
517 2821 : for (j = j0; j < n; ++j)
518 : {
519 2401 : if (gequal0(gcoeff(M, j+1, j)))
520 : {
521 1456 : for (k = j+2; k <= n; ++k)
522 1295 : if (!gequal0(gcoeff(M,k,j))) break;
523 518 : if (k > n) return j;
524 357 : transS(M, P, k, j+1);
525 : }
526 2240 : transD(M, P, j+1, j, j+1);
527 : /* Now M[j+1,j] = 1 */
528 25452 : for (k = 1; k <= n; ++k)
529 23212 : if (k != j+1 && !gequal0(gcoeff(M,k,j))) /* zero M[k,j] */
530 : {
531 13559 : transL(M, P, gneg(gcoeff(M,k,j)), k, j+1);
532 13559 : gcoeff(M,k,j) = gen_0; /* avoid approximate 0 */
533 : }
534 2240 : if (gc_needed(av,1))
535 : {
536 0 : if (DEBUGMEM > 1)
537 0 : pari_warn(warnmem,"RgM_minpoly stage 1: j0=%ld, j=%ld", j0, j);
538 0 : gerepilemat2_inplace(av, M, P);
539 : }
540 : }
541 420 : return n;
542 : }
543 :
544 : static void
545 581 : weakfrobenius_step2(GEN M, GEN P, long j)
546 : {
547 581 : pari_sp av = avma;
548 581 : long i, k, n = lg(M)-1;
549 3437 : for(i=j; i>=2; i--)
550 : {
551 6237 : for(k=j+1; k<=n; k++)
552 3381 : if (!gequal0(gcoeff(M,i,k)))
553 623 : transL(M, P, gcoeff(M,i,k), i-1, k);
554 2856 : if (gc_needed(av,1))
555 : {
556 0 : if (DEBUGMEM > 1)
557 0 : pari_warn(warnmem,"RgM_minpoly stage 2: j=%ld, i=%ld", j, i);
558 0 : gerepilemat2_inplace(av, M, P);
559 : }
560 : }
561 581 : }
562 :
563 : static long
564 581 : weakfrobenius_step3(GEN M, GEN P, long j0, long j)
565 : {
566 581 : long i, k, n = lg(M)-1;
567 581 : if (j == n) return 0;
568 161 : if (gequal0(gcoeff(M, j0, j+1)))
569 : {
570 588 : for (k=j+2; k<=n; k++)
571 448 : if (!gequal0(gcoeff(M, j0, k))) break;
572 140 : if (k > n) return 0;
573 0 : transS(M, P, k, j+1);
574 : }
575 21 : transD(M, P, j0, j+1, j+1);
576 21 : for (i=j+2; i<=n; i++)
577 0 : if (!gequal0(gcoeff(M, j0, i)))
578 0 : transL(M, P, gcoeff(M, j0, i),j+1, i);
579 21 : return 1;
580 : }
581 :
582 : /* flag: 0 -> full Frobenius from , 1 -> weak Frobenius form */
583 : static GEN
584 420 : RgM_Frobenius(GEN M, long flag, GEN *pt_P, GEN *pt_v)
585 : {
586 420 : pari_sp av = avma, av2, ltop;
587 420 : long n = lg(M)-1, eps, j0 = 1, nb = 0;
588 : GEN v, P;
589 420 : v = cgetg(n+1, t_VECSMALL);
590 420 : ltop = avma;
591 420 : P = pt_P ? matid(n): NULL;
592 420 : M = RgM_shallowcopy(M);
593 420 : av2 = avma;
594 1001 : while (j0 <= n)
595 : {
596 581 : long j = weakfrobenius_step1(M, P, j0);
597 581 : weakfrobenius_step2(M, P, j);
598 581 : eps = weakfrobenius_step3(M, P, j0, j);
599 581 : if (eps == 0)
600 : {
601 560 : v[++nb] = j0;
602 560 : if (flag == 0 && nb > 1 && !minpoly_dvdslice(M, v[nb-1], j0, j))
603 : {
604 0 : j = j0; j0 = v[nb-1]; nb -= 2;
605 0 : transL(M, P, gen_1, j, j0); /*lemma 9.18*/
606 : } else
607 560 : j0 = j+1;
608 : }
609 : else
610 21 : transS(M, P, j0, j+1); /*theorem 4*/
611 581 : if (gc_needed(av,1))
612 : {
613 0 : if (DEBUGMEM > 1)
614 0 : pari_warn(warnmem,"weakfrobenius j0=%ld",j0);
615 0 : gerepilemat2_inplace(av2, M, P);
616 : }
617 : }
618 420 : fixlg(v, nb+1);
619 420 : if (pt_v) *pt_v = v;
620 420 : gerepileall(pt_v ? ltop: av, P? 2: 1, &M, &P);
621 420 : if (pt_P) *pt_P = P;
622 420 : return M;
623 : }
624 :
625 : static GEN
626 49 : RgM_minpoly(GEN M, long v)
627 : {
628 49 : pari_sp av = avma;
629 : GEN V, W;
630 49 : if (lg(M) == 1) return pol_1(v);
631 49 : M = RgM_Frobenius(M, 1, NULL, &V);
632 49 : W = minpoly_listpolslice(M, V, v);
633 49 : if (varncmp(v,gvar2(W)) >= 0)
634 0 : pari_err_PRIORITY("matfrobenius", M, "<=", v);
635 49 : return gerepileupto(av, RgX_normalize(glcm0(W, NULL)));
636 : }
637 :
638 : GEN
639 0 : Frobeniusform(GEN V, long n)
640 : {
641 : long i, j, k;
642 0 : GEN M = zeromatcopy(n,n);
643 0 : for (k=1,i=1;i<lg(V);i++,k++)
644 : {
645 0 : GEN P = gel(V,i);
646 0 : long d = degpol(P);
647 0 : if (k+d-1 > n) pari_err_PREC("matfrobenius");
648 0 : for (j=0; j<d-1; j++, k++) gcoeff(M,k+1,k) = gen_1;
649 0 : for (j=0; j<d; j++) gcoeff(M,k-j,k) = gneg(gel(P, 1+d-j));
650 : }
651 0 : return M;
652 : }
653 :
654 : GEN
655 371 : matfrobenius(GEN M, long flag, long v)
656 : {
657 : long n;
658 371 : if (typ(M)!=t_MAT) pari_err_TYPE("matfrobenius",M);
659 371 : if (v < 0) v = 0;
660 371 : n = lg(M)-1;
661 371 : if (n && lgcols(M)!=n+1) pari_err_DIM("matfrobenius");
662 371 : if (flag > 2) pari_err_FLAG("matfrobenius");
663 371 : switch (flag)
664 : {
665 7 : case 0:
666 7 : return RgM_Frobenius(M, 0, NULL, NULL);
667 0 : case 1:
668 : {
669 0 : pari_sp av = avma;
670 : GEN V, W, F;
671 0 : F = RgM_Frobenius(M, 0, NULL, &V);
672 0 : W = minpoly_listpolslice(F, V, v);
673 0 : if (varncmp(v, gvar2(W)) >= 0)
674 0 : pari_err_PRIORITY("matfrobenius", M, "<=", v);
675 0 : return gerepileupto(av, W);
676 : }
677 364 : case 2:
678 : {
679 364 : GEN P, F, R = cgetg(3, t_VEC);
680 364 : F = RgM_Frobenius(M, 0, &P, NULL);
681 364 : gel(R,1) = F; gel(R,2) = P;
682 364 : return R;
683 : }
684 0 : default:
685 0 : pari_err_FLAG("matfrobenius");
686 : }
687 : return NULL; /*LCOV_EXCL_LINE*/
688 : }
689 :
690 : /*******************************************************************/
691 : /* */
692 : /* MINIMAL POLYNOMIAL */
693 : /* */
694 : /*******************************************************************/
695 :
696 : static GEN
697 63 : RgXQ_minpoly_naive(GEN y, GEN P)
698 : {
699 63 : long n = lgpol(P);
700 63 : GEN M = ker(RgXQ_matrix_pow(y,n,n,P));
701 63 : return content(RgM_to_RgXV(M,varn(P)));
702 : }
703 :
704 : static GEN
705 0 : RgXQ_minpoly_FpXQ(GEN x, GEN y, GEN p)
706 : {
707 0 : pari_sp av = avma;
708 : GEN r;
709 0 : if (lgefint(p) == 3)
710 : {
711 0 : ulong pp = uel(p, 2);
712 0 : r = Flx_to_ZX_inplace(Flxq_minpoly(RgX_to_Flx(x, pp),
713 : RgX_to_Flx(y, pp), pp));
714 : }
715 : else
716 0 : r = FpXQ_minpoly(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
717 0 : return gerepileupto(av, FpX_to_mod(r, p));
718 : }
719 :
720 : static GEN
721 21 : RgXQ_minpoly_FpXQXQ(GEN x, GEN S, GEN pol, GEN p)
722 : {
723 21 : pari_sp av = avma;
724 : GEN r;
725 21 : GEN T = RgX_to_FpX(pol, p);
726 21 : if (signe(T)==0) pari_err_OP("minpoly",x,S);
727 21 : if (lgefint(p) == 3)
728 : {
729 13 : ulong pp = uel(p, 2);
730 13 : GEN Tp = ZX_to_Flx(T, pp);
731 13 : r = FlxX_to_ZXX(FlxqXQ_minpoly(RgX_to_FlxqX(x, Tp, pp),
732 : RgX_to_FlxqX(S, Tp, pp), Tp, pp));
733 : }
734 : else
735 8 : r = FpXQXQ_minpoly(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(S, T, p), T, p);
736 21 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
737 : }
738 :
739 : #define code(t1,t2) ((t1 << 6) | t2)
740 :
741 : static GEN
742 1015 : RgXQ_minpoly_fast(GEN x, GEN y)
743 : {
744 : GEN p, pol;
745 : long pa;
746 1015 : long t = RgX_type2(x,y, &p,&pol,&pa);
747 1015 : switch(t)
748 : {
749 0 : case t_INTMOD: return RgXQ_minpoly_FpXQ(x, y, p);
750 77 : case t_FFELT: return FFXQ_minpoly(x, y, pol);
751 21 : case code(t_POLMOD, t_INTMOD):
752 21 : return RgXQ_minpoly_FpXQXQ(x, y, pol, p);
753 917 : default: return NULL;
754 : }
755 : }
756 :
757 : #undef code
758 :
759 : /* return caract(Mod(x,T)) in variable v */
760 : GEN
761 1015 : RgXQ_minpoly(GEN x, GEN T, long v)
762 : {
763 1015 : pari_sp av = avma;
764 1015 : GEN R = RgXQ_minpoly_fast(x,T);
765 1015 : if (R) { setvarn(R, v); return R; }
766 917 : if (!issquarefree(T))
767 : {
768 63 : R = RgXQ_minpoly_naive(x, T);
769 63 : setvarn(R,v); return R;
770 : }
771 854 : R = RgXQ_charpoly(x, T, v);
772 854 : return gerepileupto(av, RgX_div(R,RgX_gcd(R, RgX_deriv(R))));
773 : }
774 :
775 : static GEN
776 1323 : easymin(GEN x, long v)
777 : {
778 : GEN G, R, dR;
779 1323 : long tx = typ(x);
780 1323 : if (tx == t_FFELT)
781 : {
782 154 : R = FpX_to_mod(FF_minpoly(x), FF_p_i(x));
783 154 : setvarn(R,v); return R;
784 : }
785 1169 : if (tx == t_POLMOD)
786 : {
787 1120 : GEN a = gel(x,2), b = gel(x,1);
788 1120 : if (degpol(b)==0) return pol_1(v);
789 1092 : if (typ(a) != t_POL || varn(a) != varn(b))
790 : {
791 77 : if (varncmp(gvar(a), v) <= 0) pari_err_PRIORITY("minpoly", x, "<", v);
792 70 : return deg1pol(gen_1, gneg_i(a), v);
793 : }
794 1015 : return RgXQ_minpoly(a, b, v);
795 : }
796 49 : R = easychar(x, v); if (!R) return NULL;
797 0 : dR = RgX_deriv(R); if (!lgpol(dR)) return NULL;
798 0 : G = RgX_normalize(RgX_gcd(R,dR));
799 0 : return RgX_div(R,G);
800 : }
801 : GEN
802 1323 : minpoly(GEN x, long v)
803 : {
804 1323 : pari_sp av = avma;
805 : GEN P;
806 1323 : if (v < 0) v = 0;
807 1323 : P = easymin(x,v);
808 1316 : if (P) return gerepileupto(av,P);
809 : /* typ(x) = t_MAT */
810 49 : set_avma(av); return RgM_minpoly(x,v);
811 : }
812 :
813 : /*******************************************************************/
814 : /* */
815 : /* HESSENBERG FORM */
816 : /* */
817 : /*******************************************************************/
818 : static int
819 462 : relative0(GEN a, GEN a0, long bit)
820 462 : { return gequal0(a) || (bit && gexpo(a)-gexpo(a0) < bit); }
821 : /* assume x a nonempty square t_MAT */
822 : static GEN
823 70 : RgM_hess(GEN x0, long prec)
824 : {
825 : pari_sp av;
826 70 : long lx = lg(x0), bit = prec? 8 - bit_accuracy(prec): 0, m, i, j;
827 : GEN x;
828 :
829 70 : if (bit) x0 = RgM_shallowcopy(x0);
830 70 : av = avma; x = RgM_shallowcopy(x0);
831 287 : for (m=2; m<lx-1; m++)
832 : {
833 217 : GEN t = NULL;
834 560 : for (i=m+1; i<lx; i++)
835 : {
836 462 : t = gcoeff(x,i,m-1);
837 462 : if (!relative0(t, gcoeff(x0,i,m-1), bit)) break;
838 : }
839 217 : if (i == lx) continue;
840 665 : for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
841 119 : swap(gel(x,i), gel(x,m));
842 119 : if (bit)
843 : {
844 329 : for (j=m-1; j<lx; j++) swap(gcoeff(x0,i,j), gcoeff(x0,m,j));
845 63 : swap(gel(x0,i), gel(x0,m));
846 : }
847 119 : t = ginv(t);
848 :
849 427 : for (i=m+1; i<lx; i++)
850 : {
851 308 : GEN c = gcoeff(x,i,m-1);
852 308 : if (gequal0(c)) continue;
853 :
854 266 : c = gmul(c,t); gcoeff(x,i,m-1) = gen_0;
855 1330 : for (j=m; j<lx; j++)
856 1064 : gcoeff(x,i,j) = gsub(gcoeff(x,i,j), gmul(c,gcoeff(x,m,j)));
857 1981 : for (j=1; j<lx; j++)
858 1715 : gcoeff(x,j,m) = gadd(gcoeff(x,j,m), gmul(c,gcoeff(x,j,i)));
859 266 : if (gc_needed(av,2))
860 : {
861 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
862 0 : gerepileall(av,2, &x, &t);
863 : }
864 : }
865 : }
866 70 : return x;
867 : }
868 :
869 : GEN
870 70 : hess(GEN x)
871 : {
872 70 : pari_sp av = avma;
873 70 : GEN p = NULL, T = NULL;
874 70 : long prec, lx = lg(x);
875 70 : if (typ(x) != t_MAT) pari_err_TYPE("hess",x);
876 70 : if (lx == 1) return cgetg(1,t_MAT);
877 70 : if (lgcols(x) != lx) pari_err_DIM("hess");
878 70 : switch(RgM_type(x, &p, &T, &prec))
879 : {
880 42 : case t_REAL:
881 42 : case t_COMPLEX: break;
882 28 : default: prec = 0;
883 : }
884 70 : return gerepilecopy(av, RgM_hess(x,prec));
885 : }
886 :
887 : GEN
888 23351 : Flm_hess(GEN x, ulong p)
889 : {
890 23351 : long lx = lg(x), m, i, j;
891 23351 : if (lx == 1) return cgetg(1,t_MAT);
892 23351 : if (lgcols(x) != lx) pari_err_DIM("hess");
893 :
894 23351 : x = Flm_copy(x);
895 104364 : for (m=2; m<lx-1; m++)
896 : {
897 81262 : ulong t = 0;
898 286143 : for (i=m+1; i<lx; i++) { t = ucoeff(x,i,m-1); if (t) break; }
899 81262 : if (i == lx) continue;
900 454375 : for (j=m-1; j<lx; j++) lswap(ucoeff(x,i,j), ucoeff(x,m,j));
901 51164 : swap(gel(x,i), gel(x,m)); t = Fl_inv(t, p);
902 :
903 351270 : for (i=m+1; i<lx; i++)
904 : {
905 300356 : ulong c = ucoeff(x,i,m-1);
906 300356 : if (!c) continue;
907 :
908 174660 : c = Fl_mul(c,t,p); ucoeff(x,i,m-1) = 0;
909 2014050 : for (j=m; j<lx; j++)
910 1839442 : ucoeff(x,i,j) = Fl_sub(ucoeff(x,i,j), Fl_mul(c,ucoeff(x,m,j), p), p);
911 3154981 : for (j=1; j<lx; j++)
912 2980584 : ucoeff(x,j,m) = Fl_add(ucoeff(x,j,m), Fl_mul(c,ucoeff(x,j,i), p), p);
913 : }
914 : }
915 23102 : return x;
916 : }
917 : GEN
918 736 : FpM_hess(GEN x, GEN p)
919 : {
920 736 : pari_sp av = avma;
921 736 : long lx = lg(x), m, i, j;
922 736 : if (lx == 1) return cgetg(1,t_MAT);
923 736 : if (lgcols(x) != lx) pari_err_DIM("hess");
924 736 : if (lgefint(p) == 3)
925 : {
926 0 : ulong pp = p[2];
927 0 : x = Flm_hess(ZM_to_Flm(x, pp), pp);
928 0 : return gerepileupto(av, Flm_to_ZM(x));
929 : }
930 736 : x = RgM_shallowcopy(x);
931 3312 : for (m=2; m<lx-1; m++)
932 : {
933 2576 : GEN t = NULL;
934 5481 : for (i=m+1; i<lx; i++) { t = gcoeff(x,i,m-1); if (signe(t)) break; }
935 2576 : if (i == lx) continue;
936 12019 : for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
937 1890 : swap(gel(x,i), gel(x,m)); t = Fp_inv(t, p);
938 :
939 8239 : for (i=m+1; i<lx; i++)
940 : {
941 6349 : GEN c = gcoeff(x,i,m-1);
942 6349 : if (!signe(c)) continue;
943 :
944 5061 : c = Fp_mul(c,t, p); gcoeff(x,i,m-1) = gen_0;
945 29785 : for (j=m; j<lx; j++)
946 24724 : gcoeff(x,i,j) = Fp_sub(gcoeff(x,i,j), Fp_mul(c,gcoeff(x,m,j),p), p);
947 45465 : for (j=1; j<lx; j++)
948 40404 : gcoeff(x,j,m) = Fp_add(gcoeff(x,j,m), Fp_mul(c,gcoeff(x,j,i),p), p);
949 5061 : if (gc_needed(av,2))
950 : {
951 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
952 0 : gerepileall(av,2, &x, &t);
953 : }
954 : }
955 : }
956 736 : return gerepilecopy(av,x);
957 : }
958 : /* H in Hessenberg form */
959 : static GEN
960 63 : RgM_hess_charpoly(GEN H, long v)
961 : {
962 63 : long n = lg(H), r, i;
963 63 : GEN y = cgetg(n+1, t_VEC);
964 63 : gel(y,1) = pol_1(v);
965 371 : for (r = 1; r < n; r++)
966 : {
967 308 : pari_sp av2 = avma;
968 308 : GEN z, a = gen_1, b = pol_0(v);
969 812 : for (i = r-1; i; i--)
970 : {
971 581 : a = gmul(a, gcoeff(H,i+1,i));
972 581 : if (gequal0(a)) break;
973 504 : b = RgX_add(b, RgX_Rg_mul(gel(y,i), gmul(a,gcoeff(H,i,r))));
974 : }
975 308 : z = RgX_sub(RgX_shift_shallow(gel(y,r), 1),
976 308 : RgX_Rg_mul(gel(y,r), gcoeff(H,r,r)));
977 308 : gel(y,r+1) = gerepileupto(av2, RgX_sub(z, b)); /* (X - H[r,r])y[r] - b */
978 : }
979 63 : return gel(y,n);
980 : }
981 : GEN
982 63 : carhess(GEN x, long v)
983 : {
984 : pari_sp av;
985 : GEN y;
986 63 : if ((y = easychar(x,v))) return y;
987 63 : av = avma; y = RgM_hess_charpoly(hess(x), v);
988 63 : return fix_pol(av, y);
989 : }
990 :
991 : /* Bound for sup norm of charpoly(M/dM), M integral: let B = |M|oo / |dM|,
992 : * s = max_k binomial(n,k) (kB^2)^(k/2),
993 : * return ceil(log2(s)) */
994 : static long
995 3436 : charpoly_bound(GEN M, GEN dM, GEN N)
996 : {
997 3436 : pari_sp av = avma;
998 3436 : GEN B = itor(N, LOWDEFAULTPREC);
999 3436 : GEN s = real_0(LOWDEFAULTPREC), bin, B2;
1000 3436 : long n = lg(M)-1, k;
1001 3436 : bin = gen_1;
1002 3436 : if (dM) B = divri(B, dM);
1003 3436 : B2 = sqrr(B);
1004 15123 : for (k = n; k >= (n+1)>>1; k--)
1005 : {
1006 11687 : GEN t = mulri(powruhalf(mulur(k, B2), k), bin);
1007 11687 : if (abscmprr(t, s) > 0) s = t;
1008 11687 : bin = diviuexact(muliu(bin, k), n-k+1);
1009 : }
1010 3436 : return gc_long(av, ceil(dbllog2(s)));
1011 : }
1012 :
1013 : static GEN
1014 3817 : QM_charpoly_ZX_slice(GEN A, GEN dM, GEN P, GEN *mod)
1015 : {
1016 3817 : pari_sp av = avma;
1017 3817 : long i, n = lg(P)-1;
1018 : GEN H, T;
1019 3817 : if (n == 1)
1020 : {
1021 3202 : ulong p = uel(P,1), dp = dM ? umodiu(dM, p): 1;
1022 3202 : GEN Hp, a = ZM_to_Flm(A, p);
1023 3202 : Hp = Flm_charpoly_i(a, p);
1024 3202 : if (dp != 1) Hp = Flx_rescale(Hp, Fl_inv(dp, p), p);
1025 3202 : Hp = gerepileupto(av, Flx_to_ZX(Hp));
1026 3202 : *mod = utoipos(p); return Hp;
1027 : }
1028 615 : T = ZV_producttree(P);
1029 615 : A = ZM_nv_mod_tree(A, P, T);
1030 615 : H = cgetg(n+1, t_VEC);
1031 2348 : for(i=1; i <= n; i++)
1032 : {
1033 1733 : ulong p = uel(P,i), dp = dM ? umodiu(dM, p): 1;
1034 1733 : gel(H,i) = Flm_charpoly(gel(A, i), p);
1035 1734 : if (dp != 1) gel(H,i) = Flx_rescale(gel(H,i), Fl_inv(dp, p), p);
1036 : }
1037 615 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P,T));
1038 615 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
1039 : }
1040 :
1041 : GEN
1042 3817 : QM_charpoly_ZX_worker(GEN P, GEN M, GEN dM)
1043 : {
1044 3817 : GEN V = cgetg(3, t_VEC);
1045 3817 : gel(V,1) = QM_charpoly_ZX_slice(M, equali1(dM) ? NULL:dM, P, &gel(V,2));
1046 3817 : return V;
1047 : }
1048 :
1049 : /* Assume M a square ZM, dM integer. Return charpoly(M / dM) in Z[X] */
1050 : static GEN
1051 4003 : QM_charpoly_ZX_i(GEN M, GEN dM, long bound)
1052 : {
1053 4003 : long n = lg(M)-1;
1054 : forprime_t S;
1055 4003 : GEN worker = snm_closure(is_entry("_QM_charpoly_ZX_worker"),
1056 : mkvec2(M, dM? dM: gen_1));
1057 4003 : if (!n) return pol_1(0);
1058 4003 : if (bound < 0)
1059 : {
1060 3737 : GEN N = ZM_supnorm(M);
1061 3737 : if (signe(N) == 0) return monomial(gen_1, n, 0);
1062 3436 : bound = charpoly_bound(M, dM, N) + 1;
1063 : }
1064 3702 : if (DEBUGLEVEL>5) err_printf("ZM_charpoly: bound 2^%ld\n", bound);
1065 3702 : init_modular_big(&S);
1066 3702 : return gen_crt("QM_charpoly_ZX", worker, &S, dM, bound, 0, NULL,
1067 : nxV_chinese_center, FpX_center);
1068 : }
1069 :
1070 : GEN
1071 266 : QM_charpoly_ZX_bound(GEN M, long bit)
1072 : {
1073 266 : pari_sp av = avma;
1074 266 : GEN dM; M = Q_remove_denom(M, &dM);
1075 266 : return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, bit));
1076 : }
1077 : GEN
1078 1820 : QM_charpoly_ZX(GEN M)
1079 : {
1080 1820 : pari_sp av = avma;
1081 1820 : GEN dM; M = Q_remove_denom(M, &dM);
1082 1820 : return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, -1));
1083 : }
1084 : GEN
1085 1917 : ZM_charpoly(GEN M)
1086 : {
1087 1917 : pari_sp av = avma;
1088 1917 : return gerepilecopy(av, QM_charpoly_ZX_i(M, NULL, -1));
1089 : }
1090 :
1091 : /*******************************************************************/
1092 : /* */
1093 : /* CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM) */
1094 : /* */
1095 : /*******************************************************************/
1096 : GEN
1097 1540 : carberkowitz(GEN x, long v)
1098 : {
1099 : long lx, i, j, k, r;
1100 : GEN V, S, C, Q;
1101 : pari_sp av0, av;
1102 1540 : if ((V = easychar(x,v))) return V;
1103 1540 : lx = lg(x); av0 = avma;
1104 1540 : V = cgetg(lx+1, t_VEC);
1105 1540 : S = cgetg(lx+1, t_VEC);
1106 1540 : C = cgetg(lx+1, t_VEC);
1107 1540 : Q = cgetg(lx+1, t_VEC);
1108 1540 : av = avma;
1109 1540 : gel(C,1) = gen_m1;
1110 1540 : gel(V,1) = gen_m1;
1111 11711 : for (i=2;i<=lx; i++) gel(C,i) = gel(Q,i) = gel(S,i) = gel(V,i) = gen_0;
1112 1540 : gel(V,2) = gcoeff(x,1,1);
1113 10171 : for (r = 2; r < lx; r++)
1114 : {
1115 : pari_sp av2;
1116 : GEN t;
1117 :
1118 69104 : for (i = 1; i < r; i++) gel(S,i) = gcoeff(x,i,r);
1119 8631 : gel(C,2) = gcoeff(x,r,r);
1120 60473 : for (i = 1; i < r-1; i++)
1121 : {
1122 51842 : av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
1123 760088 : for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
1124 51842 : gel(C,i+2) = gerepileupto(av2, t);
1125 811930 : for (j = 1; j < r; j++)
1126 : {
1127 760088 : av2 = avma; t = gmul(gcoeff(x,j,1), gel(S,1));
1128 15031912 : for (k = 2; k < r; k++) t = gadd(t, gmul(gcoeff(x,j,k), gel(S,k)));
1129 760088 : gel(Q,j) = gerepileupto(av2, t);
1130 : }
1131 811930 : for (j = 1; j < r; j++) gel(S,j) = gel(Q,j);
1132 : }
1133 8631 : av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
1134 60473 : for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
1135 8631 : gel(C,r+1) = gerepileupto(av2, t);
1136 8631 : if (gc_needed(av0,1))
1137 : {
1138 0 : if (DEBUGMEM>1) pari_warn(warnmem,"carberkowitz");
1139 0 : gerepileall(av, 2, &C, &V);
1140 : }
1141 86366 : for (i = 1; i <= r+1; i++)
1142 : {
1143 77735 : av2 = avma; t = gmul(gel(C,i), gel(V,1));
1144 578725 : for (j = 2; j <= minss(r,i); j++)
1145 500990 : t = gadd(t, gmul(gel(C,i+1-j), gel(V,j)));
1146 77735 : gel(Q,i) = gerepileupto(av2, t);
1147 : }
1148 86366 : for (i = 1; i <= r+1; i++) gel(V,i) = gel(Q,i);
1149 : }
1150 1540 : V = RgV_to_RgX(vecreverse(V), v); /* not gtopoly: fail if v > gvar(V) */
1151 1540 : V = odd(lx)? gcopy(V): RgX_neg(V);
1152 1540 : return fix_pol(av0, V);
1153 : }
1154 :
1155 : /*******************************************************************/
1156 : /* */
1157 : /* NORMS */
1158 : /* */
1159 : /*******************************************************************/
1160 : GEN
1161 3286003 : gnorm(GEN x)
1162 : {
1163 : pari_sp av;
1164 : long lx, i;
1165 : GEN y;
1166 :
1167 3286003 : switch(typ(x))
1168 : {
1169 42533 : case t_INT: return sqri(x);
1170 334031 : case t_REAL: return sqrr(x);
1171 1337 : case t_FRAC: return sqrfrac(x);
1172 2834864 : case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
1173 69055 : case t_QUAD: av = avma; return gerepileupto(av, quadnorm(x));
1174 :
1175 14 : case t_POL: case t_SER: case t_RFRAC: av = avma;
1176 14 : return gerepileupto(av, greal(gmul(conj_i(x),x)));
1177 :
1178 28 : case t_FFELT:
1179 28 : y = cgetg(3, t_INTMOD);
1180 28 : gel(y,1) = FF_p(x);
1181 28 : gel(y,2) = FF_norm(x); return y;
1182 :
1183 4144 : case t_POLMOD:
1184 : {
1185 4144 : GEN T = gel(x,1), a = gel(x,2);
1186 4144 : if (typ(a) != t_POL || varn(a) != varn(T)) return gpowgs(a, degpol(T));
1187 3962 : return RgXQ_norm(a, T);
1188 : }
1189 0 : case t_VEC: case t_COL: case t_MAT:
1190 0 : y = cgetg_copy(x, &lx);
1191 0 : for (i=1; i<lx; i++) gel(y,i) = gnorm(gel(x,i));
1192 0 : return y;
1193 : }
1194 0 : pari_err_TYPE("gnorm",x);
1195 : return NULL; /* LCOV_EXCL_LINE */
1196 : }
1197 :
1198 : /* return |q|^2, complex modulus */
1199 : static GEN
1200 28 : cxquadnorm(GEN q, long prec)
1201 : {
1202 28 : GEN X = gel(q,1), c = gel(X,2); /* (1-D)/4, -D/4 */
1203 28 : if (signe(c) > 0) return quadnorm(q); /* imaginary */
1204 21 : if (!prec) pari_err_TYPE("gnorml2", q);
1205 7 : return sqrr(quadtofp(q, prec));
1206 : }
1207 :
1208 : static GEN
1209 55312312 : gnorml2_i(GEN x, long prec)
1210 : {
1211 : pari_sp av;
1212 : long i, lx;
1213 : GEN s;
1214 :
1215 55312312 : switch(typ(x))
1216 : {
1217 20265174 : case t_INT: return sqri(x);
1218 22827266 : case t_REAL: return sqrr(x);
1219 7 : case t_FRAC: return sqrfrac(x);
1220 4114770 : case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
1221 21 : case t_QUAD: av = avma; return gerepileupto(av, cxquadnorm(x,prec));
1222 :
1223 57611 : case t_POL: lx = lg(x)-1; x++; break;
1224 :
1225 8054992 : case t_VEC:
1226 : case t_COL:
1227 8054992 : case t_MAT: lx = lg(x); break;
1228 :
1229 0 : default: pari_err_TYPE("gnorml2",x);
1230 : return NULL; /* LCOV_EXCL_LINE */
1231 : }
1232 8112603 : if (lx == 1) return gen_0;
1233 8112603 : av = avma;
1234 8112603 : s = gnorml2(gel(x,1));
1235 47268218 : for (i=2; i<lx; i++)
1236 : {
1237 39157479 : s = gadd(s, gnorml2(gel(x,i)));
1238 39156157 : if (gc_needed(av,1))
1239 : {
1240 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gnorml2");
1241 0 : s = gerepileupto(av, s);
1242 : }
1243 : }
1244 8110739 : return gerepileupto(av,s);
1245 : }
1246 : GEN
1247 55309794 : gnorml2(GEN x) { return gnorml2_i(x, 0); }
1248 :
1249 : static GEN pnormlp(GEN,GEN,long);
1250 : static GEN
1251 63 : pnormlpvec(long i0, GEN x, GEN p, long prec)
1252 : {
1253 63 : pari_sp av = avma;
1254 63 : long i, lx = lg(x);
1255 63 : GEN s = gen_0;
1256 224 : for (i=i0; i<lx; i++)
1257 : {
1258 161 : s = gadd(s, pnormlp(gel(x,i),p,prec));
1259 161 : if (gc_needed(av,1))
1260 : {
1261 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gnormlp, i = %ld", i);
1262 0 : s = gerepileupto(av, s);
1263 : }
1264 : }
1265 63 : return s;
1266 : }
1267 : /* (||x||_p)^p */
1268 : static GEN
1269 196 : pnormlp(GEN x, GEN p, long prec)
1270 : {
1271 196 : switch(typ(x))
1272 : {
1273 119 : case t_INT: case t_REAL: x = mpabs(x); break;
1274 0 : case t_FRAC: x = absfrac(x); break;
1275 14 : case t_COMPLEX: case t_QUAD: x = gabs(x,prec); break;
1276 7 : case t_POL: return pnormlpvec(2, x, p, prec);
1277 56 : case t_VEC: case t_COL: case t_MAT: return pnormlpvec(1, x, p, prec);
1278 0 : default: pari_err_TYPE("gnormlp",x);
1279 : }
1280 133 : return gpow(x, p, prec);
1281 : }
1282 :
1283 : GEN
1284 343 : gnormlp(GEN x, GEN p, long prec)
1285 : {
1286 343 : pari_sp av = avma;
1287 343 : if (!p || (typ(p) == t_INFINITY && inf_get_sign(p) > 0))
1288 182 : return gsupnorm(x, prec);
1289 161 : if (gsigne(p) <= 0) pari_err_DOMAIN("normlp", "p", "<=", gen_0, p);
1290 154 : if (is_scalar_t(typ(x))) return gabs(x, prec);
1291 91 : if (typ(p) == t_INT)
1292 : {
1293 63 : ulong pp = itou_or_0(p);
1294 63 : switch(pp)
1295 : {
1296 28 : case 1: return gnorml1(x, prec);
1297 28 : case 2: x = gnorml2_i(x, prec); break;
1298 7 : default: x = pnormlp(x, p, prec); break;
1299 : }
1300 35 : if (pp && typ(x) == t_INT && Z_ispowerall(x, pp, &x))
1301 7 : return gerepileuptoleaf(av, x);
1302 28 : if (pp == 2) return gerepileupto(av, gsqrt(x, prec));
1303 : }
1304 : else
1305 28 : x = pnormlp(x, p, prec);
1306 28 : x = gpow(x, ginv(p), prec);
1307 28 : return gerepileupto(av, x);
1308 : }
1309 :
1310 : GEN
1311 168 : gnorml1(GEN x,long prec)
1312 : {
1313 168 : pari_sp av = avma;
1314 : long lx,i;
1315 : GEN s;
1316 168 : switch(typ(x))
1317 : {
1318 98 : case t_INT: case t_REAL: return mpabs(x);
1319 0 : case t_FRAC: return absfrac(x);
1320 :
1321 14 : case t_COMPLEX: case t_QUAD:
1322 14 : return gabs(x,prec);
1323 :
1324 7 : case t_POL:
1325 7 : lx = lg(x); s = gen_0;
1326 28 : for (i=2; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
1327 7 : break;
1328 :
1329 49 : case t_VEC: case t_COL: case t_MAT:
1330 49 : lx = lg(x); s = gen_0;
1331 168 : for (i=1; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
1332 49 : break;
1333 :
1334 0 : default: pari_err_TYPE("gnorml1",x);
1335 : return NULL; /* LCOV_EXCL_LINE */
1336 : }
1337 56 : return gerepileupto(av, s);
1338 : }
1339 : /* As gnorml1, except for t_QUAD and t_COMPLEX: |x + wy| := |x| + |y|
1340 : * Still a norm of R-vector spaces, and can be cheaply computed without
1341 : * square roots */
1342 : GEN
1343 142513 : gnorml1_fake(GEN x)
1344 : {
1345 142513 : pari_sp av = avma;
1346 : long lx, i;
1347 : GEN s;
1348 142513 : switch(typ(x))
1349 : {
1350 126742 : case t_INT: case t_REAL: return mpabs(x);
1351 0 : case t_FRAC: return absfrac(x);
1352 :
1353 0 : case t_COMPLEX:
1354 0 : s = gadd(gnorml1_fake(gel(x,1)), gnorml1_fake(gel(x,2)));
1355 0 : break;
1356 0 : case t_QUAD:
1357 0 : s = gadd(gnorml1_fake(gel(x,2)), gnorml1_fake(gel(x,3)));
1358 0 : break;
1359 :
1360 15771 : case t_POL:
1361 15771 : lx = lg(x); s = gen_0;
1362 102893 : for (i=2; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
1363 15771 : break;
1364 :
1365 0 : case t_VEC: case t_COL: case t_MAT:
1366 0 : lx = lg(x); s = gen_0;
1367 0 : for (i=1; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
1368 0 : break;
1369 :
1370 0 : default: pari_err_TYPE("gnorml1_fake",x);
1371 : return NULL; /* LCOV_EXCL_LINE */
1372 : }
1373 15771 : return gerepileupto(av, s);
1374 : }
1375 :
1376 : static void
1377 207585 : store(GEN z, GEN *m) { if (!*m || gcmp(z, *m) > 0) *m = z; }
1378 : /* compare |x| to *m or |x|^2 to *msq, whichever is easiest, and update
1379 : * the pointed value if x is larger */
1380 : void
1381 245553 : gsupnorm_aux(GEN x, GEN *m, GEN *msq, long prec)
1382 : {
1383 : long i, lx;
1384 : GEN z;
1385 245553 : switch(typ(x))
1386 : {
1387 92513 : case t_COMPLEX: z = cxnorm(x); store(z, msq); return;
1388 7 : case t_QUAD: z = cxquadnorm(x,prec); store(z, msq); return;
1389 115067 : case t_INT: case t_REAL: z = mpabs(x); store(z,m); return;
1390 0 : case t_FRAC: z = absfrac(x); store(z,m); return;
1391 :
1392 6034 : case t_POL: lx = lg(x)-1; x++; break;
1393 :
1394 31935 : case t_VEC:
1395 : case t_COL:
1396 31935 : case t_MAT: lx = lg(x); break;
1397 :
1398 0 : default: pari_err_TYPE("gsupnorm",x);
1399 : return; /* LCOV_EXCL_LINE */
1400 : }
1401 245932 : for (i=1; i<lx; i++) gsupnorm_aux(gel(x,i), m, msq, prec);
1402 : }
1403 : GEN
1404 37591 : gsupnorm(GEN x, long prec)
1405 : {
1406 37591 : GEN m = NULL, msq = NULL;
1407 37591 : pari_sp av = avma;
1408 37591 : gsupnorm_aux(x, &m, &msq, prec);
1409 : /* now set m = max (m, sqrt(msq)) */
1410 37593 : if (msq) {
1411 15181 : msq = gsqrt(msq, prec);
1412 15183 : if (!m || gcmp(m, msq) < 0) m = msq;
1413 22412 : } else if (!m) m = gen_0;
1414 37595 : return gerepilecopy(av, m);
1415 : }
1416 :
1417 : /*******************************************************************/
1418 : /* */
1419 : /* TRACES */
1420 : /* */
1421 : /*******************************************************************/
1422 : GEN
1423 35 : matcompanion(GEN x)
1424 : {
1425 35 : long n = degpol(x), j;
1426 : GEN y, c;
1427 :
1428 35 : if (typ(x)!=t_POL) pari_err_TYPE("matcompanion",x);
1429 35 : if (!signe(x)) pari_err_DOMAIN("matcompanion","polynomial","=",gen_0,x);
1430 28 : if (n == 0) return cgetg(1, t_MAT);
1431 :
1432 28 : y = cgetg(n+1,t_MAT);
1433 105 : for (j=1; j < n; j++) gel(y,j) = col_ei(n, j+1);
1434 28 : c = cgetg(n+1,t_COL); gel(y,n) = c;
1435 28 : if (gequal1(gel(x, n+2)))
1436 112 : for (j=1; j<=n; j++) gel(c,j) = gneg(gel(x,j+1));
1437 : else
1438 : { /* not monic. Hardly ever used */
1439 7 : pari_sp av = avma;
1440 7 : GEN d = gclone(gneg(gel(x,n+2)));
1441 7 : set_avma(av);
1442 21 : for (j=1; j<=n; j++) gel(c,j) = gdiv(gel(x,j+1), d);
1443 7 : gunclone(d);
1444 : }
1445 28 : return y;
1446 : }
1447 :
1448 : GEN
1449 712759 : gtrace(GEN x)
1450 : {
1451 : pari_sp av;
1452 712759 : long i, lx, tx = typ(x);
1453 : GEN y, z;
1454 :
1455 712759 : switch(tx)
1456 : {
1457 23113 : case t_INT: case t_REAL: case t_FRAC:
1458 23113 : return gmul2n(x,1);
1459 :
1460 687080 : case t_COMPLEX:
1461 687080 : return gmul2n(gel(x,1),1);
1462 :
1463 154 : case t_QUAD:
1464 154 : y = gel(x,1);
1465 154 : if (!gequal0(gel(y,3)))
1466 : { /* assume quad. polynomial is either x^2 + d or x^2 - x + d */
1467 154 : av = avma;
1468 154 : return gerepileupto(av, gadd(gel(x,3), gmul2n(gel(x,2),1)));
1469 : }
1470 0 : return gmul2n(gel(x,2),1);
1471 :
1472 7 : case t_POL:
1473 7 : y = cgetg_copy(x, &lx); y[1] = x[1];
1474 21 : for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
1475 7 : return normalizepol_lg(y, lx);
1476 :
1477 14 : case t_SER:
1478 14 : if (ser_isexactzero(x)) return gcopy(x);
1479 7 : y = cgetg_copy(x, &lx); y[1] = x[1];
1480 21 : for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
1481 7 : return normalizeser(y);
1482 :
1483 644 : case t_POLMOD:
1484 644 : y = gel(x,1); z = gel(x,2);
1485 644 : if (typ(z) != t_POL || varn(y) != varn(z)) return gmulsg(degpol(y), z);
1486 413 : av = avma;
1487 413 : return gerepileupto(av, RgXQ_trace(z, y));
1488 :
1489 28 : case t_FFELT:
1490 28 : y=cgetg(3, t_INTMOD);
1491 28 : gel(y,1) = FF_p(x);
1492 28 : gel(y,2) = FF_trace(x);
1493 28 : return y;
1494 :
1495 7 : case t_RFRAC:
1496 7 : av = avma; return gerepileupto(av, gadd(x, conj_i(x)));
1497 :
1498 0 : case t_VEC: case t_COL:
1499 0 : y = cgetg_copy(x, &lx);
1500 0 : for (i=1; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
1501 0 : return y;
1502 :
1503 1715 : case t_MAT:
1504 1715 : lx = lg(x); if (lx == 1) return gen_0;
1505 : /*now lx >= 2*/
1506 1708 : if (lx != lgcols(x)) pari_err_DIM("gtrace");
1507 1701 : av = avma; return gerepileupto(av, mattrace(x));
1508 : }
1509 0 : pari_err_TYPE("gtrace",x);
1510 : return NULL; /* LCOV_EXCL_LINE */
1511 : }
1512 :
1513 : /* Cholesky decomposition for positive definite matrix a
1514 : * [GTM138, Algo 2.7.6, matrix Q]
1515 : * If a is not positive definite return NULL. */
1516 : GEN
1517 1328 : qfgaussred_positive(GEN a)
1518 : {
1519 1328 : pari_sp av = avma;
1520 : GEN b;
1521 1328 : long i,j,k, n = lg(a);
1522 :
1523 1328 : if (typ(a)!=t_MAT) pari_err_TYPE("qfgaussred_positive",a);
1524 1328 : if (n == 1) return cgetg(1, t_MAT);
1525 1321 : if (lgcols(a)!=n) pari_err_DIM("qfgaussred_positive");
1526 1321 : b = cgetg(n,t_MAT);
1527 8232 : for (j=1; j<n; j++)
1528 : {
1529 6911 : GEN p1=cgetg(n,t_COL), p2=gel(a,j);
1530 :
1531 6911 : gel(b,j) = p1;
1532 41321 : for (i=1; i<=j; i++) gel(p1,i) = gel(p2,i);
1533 34410 : for ( ; i<n ; i++) gel(p1,i) = gen_0;
1534 : }
1535 8232 : for (k=1; k<n; k++)
1536 : {
1537 6911 : GEN bk, p = gcoeff(b,k,k), invp;
1538 6911 : if (gsigne(p)<=0) return gc_NULL(av); /* not positive definite */
1539 6911 : invp = ginv(p);
1540 6911 : bk = row(b, k);
1541 34410 : for (i=k+1; i<n; i++) gcoeff(b,k,i) = gmul(gel(bk,i), invp);
1542 34410 : for (i=k+1; i<n; i++)
1543 : {
1544 27499 : GEN c = gel(bk, i);
1545 155461 : for (j=i; j<n; j++)
1546 127962 : gcoeff(b,i,j) = gsub(gcoeff(b,i,j), gmul(c,gcoeff(b,k,j)));
1547 : }
1548 6911 : if (gc_needed(av,1))
1549 : {
1550 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfgaussred_positive");
1551 0 : b=gerepilecopy(av,b);
1552 : }
1553 : }
1554 1321 : return gerepilecopy(av,b);
1555 : }
1556 :
1557 : /* Maximal pivot strategy: x is a suitable pivot if it is non zero and either
1558 : * - an exact type, or
1559 : * - it is maximal among remaining nonzero (t_REAL) pivots */
1560 : static int
1561 46830 : suitable(GEN x, long k, GEN *pp, long *pi)
1562 : {
1563 46830 : long t = typ(x);
1564 46830 : switch(t)
1565 : {
1566 24374 : case t_INT: return signe(x) != 0;
1567 22302 : case t_FRAC: return 1;
1568 154 : case t_REAL: {
1569 154 : GEN p = *pp;
1570 154 : if (signe(x) && (!p || abscmprr(p, x) < 0)) { *pp = x; *pi = k; }
1571 154 : return 0;
1572 : }
1573 0 : default: return !gequal0(x);
1574 : }
1575 : }
1576 :
1577 : /* Gauss reduction (arbitrary symetric matrix, only the part above the
1578 : * diagonal is considered). If signature is nonzero, return only the
1579 : * signature, in which case gsigne() should be defined for elements of a. */
1580 : static GEN
1581 12012 : gaussred(GEN a, long signature)
1582 : {
1583 : GEN r, ak, al;
1584 12012 : pari_sp av = avma, av1;
1585 12012 : long n = lg(a), i, j, k, l, sp, sn, t;
1586 :
1587 12012 : if (typ(a) != t_MAT) pari_err_TYPE("gaussred",a);
1588 12012 : if (n == 1) return signature? mkvec2(gen_0, gen_0): cgetg(1, t_MAT);
1589 12012 : if (lgcols(a) != n) pari_err_DIM("gaussred");
1590 12012 : n--;
1591 12012 : r = const_vecsmall(n, 1); av1= avma;
1592 12012 : a = RgM_shallowcopy(a);
1593 12012 : t = n; sp = sn = 0;
1594 58702 : while (t)
1595 : {
1596 46690 : long pind = 0;
1597 46690 : GEN invp, p = NULL;
1598 129381 : k=1; while (k<=n && (!r[k] || !suitable(gcoeff(a,k,k), k, &p, &pind))) k++;
1599 46690 : if (k > n && p) k = pind;
1600 46690 : if (k <= n)
1601 : {
1602 46683 : p = gcoeff(a,k,k); invp = ginv(p); /* != 0 */
1603 46683 : if (signature) { /* skip if (!signature): gsigne may fail ! */
1604 46634 : if (gsigne(p) > 0) sp++; else sn++;
1605 : }
1606 46683 : r[k] = 0; t--;
1607 46683 : ak = row(a, k);
1608 258356 : for (i=1; i<=n; i++)
1609 211673 : gcoeff(a,k,i) = r[i]? gmul(gcoeff(a,k,i), invp): gen_0;
1610 :
1611 258356 : for (i=1; i<=n; i++) if (r[i])
1612 : {
1613 82481 : GEN c = gel(ak,i); /* - p * a[k,i] */
1614 82481 : if (gequal0(c)) continue;
1615 461405 : for (j=1; j<=n; j++) if (r[j])
1616 244048 : gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
1617 : }
1618 46683 : gcoeff(a,k,k) = p;
1619 46683 : if (gc_needed(av1,1))
1620 : {
1621 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gaussred (t = %ld)", t);
1622 0 : a = gerepilecopy(av1, a);
1623 : }
1624 : }
1625 : else
1626 : { /* all remaining diagonal coeffs are currently 0 */
1627 7 : for (k=1; k<=n; k++) if (r[k])
1628 : {
1629 7 : l=k+1; while (l<=n && (!r[l] || !suitable(gcoeff(a,k,l), l, &p, &pind))) l++;
1630 7 : if (l > n && p) l = pind;
1631 7 : if (l > n) continue;
1632 :
1633 7 : p = gcoeff(a,k,l); invp = ginv(p);
1634 7 : sp++; sn++;
1635 7 : r[k] = r[l] = 0; t -= 2;
1636 7 : ak = row(a, k);
1637 7 : al = row(a, l);
1638 35 : for (i=1; i<=n; i++) if (r[i])
1639 : {
1640 14 : gcoeff(a,k,i) = gmul(gcoeff(a,k,i), invp);
1641 14 : gcoeff(a,l,i) = gmul(gcoeff(a,l,i), invp);
1642 : } else {
1643 14 : gcoeff(a,k,i) = gen_0;
1644 14 : gcoeff(a,l,i) = gen_0;
1645 : }
1646 :
1647 35 : for (i=1; i<=n; i++) if (r[i])
1648 : { /* c = a[k,i] * p, d = a[l,i] * p; */
1649 14 : GEN c = gel(ak,i), d = gel(al,i);
1650 70 : for (j=1; j<=n; j++) if (r[j])
1651 28 : gcoeff(a,i,j) = gsub(gcoeff(a,i,j),
1652 28 : gadd(gmul(gcoeff(a,l,j), c),
1653 28 : gmul(gcoeff(a,k,j), d)));
1654 : }
1655 35 : for (i=1; i<=n; i++) if (r[i])
1656 : {
1657 14 : GEN c = gcoeff(a,k,i), d = gcoeff(a,l,i);
1658 14 : gcoeff(a,k,i) = gadd(c, d);
1659 14 : gcoeff(a,l,i) = gsub(c, d);
1660 : }
1661 7 : gcoeff(a,k,l) = gen_1;
1662 7 : gcoeff(a,l,k) = gen_m1;
1663 7 : gcoeff(a,k,k) = gmul2n(p,-1);
1664 7 : gcoeff(a,l,l) = gneg(gcoeff(a,k,k));
1665 7 : if (gc_needed(av1,1))
1666 : {
1667 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gaussred");
1668 0 : a = gerepilecopy(av1, a);
1669 : }
1670 7 : break;
1671 : }
1672 7 : if (k > n) break;
1673 : }
1674 : }
1675 12012 : if (!signature) return gerepilecopy(av, a);
1676 11998 : set_avma(av); return mkvec2s(sp, sn);
1677 : }
1678 :
1679 : GEN
1680 14 : qfgaussred(GEN a) { return gaussred(a,0); }
1681 :
1682 : GEN
1683 11998 : qfsign(GEN a) { return gaussred(a,1); }
1684 :
1685 : /* x -= s(y+u*x) */
1686 : /* y += s(x-u*y), simultaneously */
1687 : static void
1688 19180 : rot(GEN x, GEN y, GEN s, GEN u) {
1689 19180 : GEN x1 = subrr(x, mulrr(s,addrr(y,mulrr(u,x))));
1690 19180 : GEN y1 = addrr(y, mulrr(s,subrr(x,mulrr(u,y))));
1691 19180 : affrr(x1,x);
1692 19180 : affrr(y1,y);
1693 19180 : }
1694 :
1695 : /* Diagonalization of a REAL symetric matrix. Return a vector [L, r]:
1696 : * L = vector of eigenvalues
1697 : * r = matrix of eigenvectors */
1698 : GEN
1699 28 : jacobi(GEN a, long prec)
1700 : {
1701 : pari_sp av;
1702 28 : long de, e, e1, e2, i, j, p, q, l = lg(a);
1703 : GEN c, ja, L, r, L2, r2, unr;
1704 :
1705 28 : if (typ(a) != t_MAT) pari_err_TYPE("jacobi",a);
1706 28 : ja = cgetg(3,t_VEC);
1707 28 : L = cgetg(l,t_COL); gel(ja,1) = L;
1708 28 : r = cgetg(l,t_MAT); gel(ja,2) = r;
1709 28 : if (l == 1) return ja;
1710 28 : if (lgcols(a) != l) pari_err_DIM("jacobi");
1711 :
1712 28 : e1 = HIGHEXPOBIT-1;
1713 224 : for (j=1; j<l; j++)
1714 : {
1715 196 : GEN z = gtofp(gcoeff(a,j,j), prec);
1716 196 : gel(L,j) = z;
1717 196 : e = expo(z); if (e < e1) e1 = e;
1718 : }
1719 224 : for (j=1; j<l; j++)
1720 : {
1721 196 : gel(r,j) = cgetg(l,t_COL);
1722 1582 : for (i=1; i<l; i++) gcoeff(r,i,j) = utor(i==j? 1: 0, prec);
1723 : }
1724 28 : av = avma;
1725 :
1726 28 : e2 = -(long)HIGHEXPOBIT; p = q = 1;
1727 28 : c = cgetg(l,t_MAT);
1728 224 : for (j=1; j<l; j++)
1729 : {
1730 196 : gel(c,j) = cgetg(j,t_COL);
1731 791 : for (i=1; i<j; i++)
1732 : {
1733 595 : GEN z = gtofp(gcoeff(a,i,j), prec);
1734 595 : gcoeff(c,i,j) = z;
1735 595 : if (!signe(z)) continue;
1736 308 : e = expo(z); if (e > e2) { e2 = e; p = i; q = j; }
1737 : }
1738 : }
1739 28 : a = c; unr = real_1(prec);
1740 28 : de = prec2nbits(prec);
1741 :
1742 : /* e1 = min expo(a[i,i])
1743 : * e2 = max expo(a[i,j]), i != j */
1744 1568 : while (e1-e2 < de)
1745 : {
1746 1540 : pari_sp av2 = avma;
1747 : GEN x, y, t, c, s, u;
1748 : /* compute attached rotation in the plane formed by basis vectors number
1749 : * p and q */
1750 1540 : x = subrr(gel(L,q),gel(L,p));
1751 1540 : if (signe(x))
1752 : {
1753 1512 : x = divrr(x, shiftr(gcoeff(a,p,q),1));
1754 1512 : y = sqrtr(addrr(unr, sqrr(x)));
1755 1512 : t = invr((signe(x)>0)? addrr(x,y): subrr(x,y));
1756 : }
1757 : else
1758 28 : y = t = unr;
1759 1540 : c = sqrtr(addrr(unr,sqrr(t)));
1760 1540 : s = divrr(t,c);
1761 1540 : u = divrr(t,addrr(unr,c));
1762 :
1763 : /* compute successive transforms of a and the matrix of accumulated
1764 : * rotations (r) */
1765 4144 : for (i=1; i<p; i++) rot(gcoeff(a,i,p), gcoeff(a,i,q), s,u);
1766 4039 : for (i=p+1; i<q; i++) rot(gcoeff(a,p,i), gcoeff(a,i,q), s,u);
1767 4487 : for (i=q+1; i<l; i++) rot(gcoeff(a,p,i), gcoeff(a,q,i), s,u);
1768 1540 : y = gcoeff(a,p,q);
1769 1540 : t = mulrr(t, y); shiftr_inplace(y, -de - 1);
1770 1540 : x = gel(L,p); subrrz(x,t, x);
1771 1540 : y = gel(L,q); addrrz(y,t, y);
1772 12670 : for (i=1; i<l; i++) rot(gcoeff(r,i,p), gcoeff(r,i,q), s,u);
1773 :
1774 1540 : e2 = -(long)HIGHEXPOBIT; p = q = 1;
1775 12670 : for (j=1; j<l; j++)
1776 : {
1777 46396 : for (i=1; i<j; i++)
1778 : {
1779 35266 : GEN z = gcoeff(a,i,j);
1780 35266 : if (!signe(z)) continue;
1781 31080 : e = expo(z); if (e > e2) { e2=e; p=i; q=j; }
1782 : }
1783 46396 : for (i=j+1; i<l; i++)
1784 : {
1785 35266 : GEN z = gcoeff(a,j,i);
1786 35266 : if (!signe(z)) continue;
1787 31080 : e = expo(z); if (e > e2) { e2=e; p=j; q=i; }
1788 : }
1789 : }
1790 1540 : set_avma(av2);
1791 : }
1792 : /* sort eigenvalues from smallest to largest */
1793 28 : c = indexsort(L);
1794 224 : r2 = vecpermute(r, c); for (i=1; i<l; i++) gel(r,i) = gel(r2,i);
1795 224 : L2 = vecpermute(L, c); for (i=1; i<l; i++) gel(L,i) = gel(L2,i);
1796 28 : set_avma(av); return ja;
1797 : }
1798 :
1799 : /*************************************************************************/
1800 : /** **/
1801 : /** Q-vector space -> Z-modules **/
1802 : /** **/
1803 : /*************************************************************************/
1804 :
1805 : GEN
1806 133 : matrixqz0(GEN x,GEN p)
1807 : {
1808 133 : if (typ(x) != t_MAT) pari_err_TYPE("matrixqz",x);
1809 133 : if (!p) return QM_minors_coprime(x,NULL);
1810 98 : if (typ(p) != t_INT) pari_err_TYPE("matrixqz",p);
1811 98 : if (signe(p)>=0) return QM_minors_coprime(x,p);
1812 91 : if (!RgM_is_QM(x)) pari_err_TYPE("matrixqz", x);
1813 91 : if (absequaliu(p,1)) return QM_ImZ(x); /* p = -1 */
1814 63 : if (absequaliu(p,2)) return QM_ImQ(x); /* p = -2 */
1815 7 : pari_err_FLAG("QM_minors_coprime");
1816 : return NULL; /* LCOV_EXCL_LINE */
1817 : }
1818 :
1819 : GEN
1820 42 : QM_minors_coprime(GEN x, GEN D)
1821 : {
1822 42 : pari_sp av = avma, av1;
1823 : long i, j, m, n, lP;
1824 : GEN P, y;
1825 :
1826 42 : n = lg(x)-1; if (!n) return gcopy(x);
1827 42 : m = nbrows(x);
1828 42 : if (n > m) pari_err_DOMAIN("QM_minors_coprime","n",">",strtoGENstr("m"),x);
1829 35 : y = x; x = cgetg(n+1,t_MAT);
1830 112 : for (j=1; j<=n; j++)
1831 : {
1832 77 : gel(x,j) = Q_primpart(gel(y,j));
1833 77 : RgV_check_ZV(gel(x,j), "QM_minors_coprime");
1834 : }
1835 : /* x now a ZM */
1836 35 : if (n==m)
1837 : {
1838 21 : if (gequal0(ZM_det(x)))
1839 14 : pari_err_DOMAIN("QM_minors_coprime", "rank(A)", "<",stoi(n),x);
1840 7 : set_avma(av); return matid(n);
1841 : }
1842 : /* m > n */
1843 14 : if (!D || gequal0(D))
1844 : {
1845 14 : pari_sp av2 = avma;
1846 14 : D = ZM_detmult(shallowtrans(x));
1847 14 : if (is_pm1(D)) { set_avma(av2); return ZM_copy(x); }
1848 : }
1849 14 : P = gel(Z_factor(D), 1); lP = lg(P);
1850 14 : av1 = avma;
1851 56 : for (i=1; i < lP; i++)
1852 : {
1853 42 : GEN p = gel(P,i), pov2 = shifti(p, -1);
1854 : for(;;)
1855 42 : {
1856 84 : GEN N, M = FpM_ker(x, p);
1857 84 : long lM = lg(M);
1858 84 : if (lM==1) break;
1859 :
1860 42 : FpM_center_inplace(M, p, pov2);
1861 42 : N = ZM_Z_divexact(ZM_mul(x,M), p);
1862 126 : for (j=1; j<lM; j++)
1863 : {
1864 147 : long k = n; while (!signe(gcoeff(M,k,j))) k--;
1865 84 : gel(x,k) = gel(N,j);
1866 : }
1867 42 : if (gc_needed(av1,1))
1868 : {
1869 0 : if (DEBUGMEM>1) pari_warn(warnmem,"QM_minors_coprime, p = %Ps", p);
1870 0 : x = gerepilecopy(av1, x); pov2 = shifti(p, -1);
1871 : }
1872 : }
1873 : }
1874 14 : return gerepilecopy(av, x);
1875 : }
1876 :
1877 : static GEN
1878 1162 : QM_ImZ_all_i(GEN A, GEN *U, long remove, long hnf, long linindep)
1879 : {
1880 1162 : GEN V = NULL, D;
1881 1162 : A = Q_remove_denom(A,&D);
1882 1162 : if (D)
1883 : {
1884 : long l, lA;
1885 126 : V = matkermod(A,D,NULL);
1886 126 : l = lg(V); lA = lg(A);
1887 126 : if (l == 1) V = scalarmat_shallow(D, lA-1);
1888 : else
1889 : {
1890 84 : if (l < lA) V = hnfmodid(V,D);
1891 84 : A = ZM_Z_divexact(ZM_mul(A, V), D);
1892 : }
1893 : }
1894 1162 : if (!linindep && ZM_rank(A)==lg(A)-1) linindep = 1;
1895 1162 : if (hnf || !linindep) A = ZM_hnflll(A, U, remove);
1896 1162 : if (U && V)
1897 : {
1898 35 : if (hnf) *U = ZM_mul(V,*U);
1899 0 : else *U = V;
1900 : }
1901 1162 : return A;
1902 : }
1903 : GEN
1904 28 : QM_ImZ_all(GEN x, GEN *U, long remove, long hnf)
1905 : {
1906 28 : pari_sp av = avma;
1907 28 : x = QM_ImZ_all_i(x, U, remove, hnf, 0);
1908 28 : return gc_all(av, U?2:1, &x, &U);
1909 : }
1910 : GEN
1911 0 : QM_ImZ_hnfall(GEN x, GEN *U, long remove) { return QM_ImZ_all(x, U, remove, 1); }
1912 : GEN
1913 0 : QM_ImZ_hnf(GEN x) { return QM_ImZ_hnfall(x, NULL, 1); }
1914 : GEN
1915 28 : QM_ImZ(GEN x) { return QM_ImZ_all(x, NULL, 1, 0); }
1916 :
1917 : GEN
1918 1141 : QM_ImQ_all(GEN x, GEN *U, long remove, long hnf)
1919 : {
1920 1141 : pari_sp av = avma;
1921 1141 : long i, n = lg(x), m;
1922 1141 : GEN ir, V, D, c, K = NULL;
1923 :
1924 1141 : if (U) *U = matid(n-1);
1925 1141 : if (n==1) return gcopy(x);
1926 1134 : m = lg(gel(x,1));
1927 :
1928 1134 : x = RgM_shallowcopy(x);
1929 7287 : for (i=1; i<n; i++)
1930 : {
1931 6153 : gel(x,i) = Q_primitive_part(gel(x,i), &c);
1932 6153 : if (U && c && signe(c)) gcoeff(*U,i,i) = ginv(c);
1933 : }
1934 :
1935 1134 : ir = ZM_indexrank(x);
1936 1134 : if (U)
1937 : {
1938 77 : *U = vecpermute(*U, gel(ir,2));
1939 77 : if (remove < 2) K = ZM_ker(x);
1940 : }
1941 1134 : x = vecpermute(x, gel(ir,2));
1942 :
1943 1134 : D = absi(ZM_det(rowpermute(x,gel(ir,1))));
1944 1134 : x = RgM_Rg_div(x, D);
1945 1134 : x = QM_ImZ_all_i(x, U? &V: NULL, remove, hnf, 1);
1946 :
1947 1134 : if (U)
1948 : {
1949 77 : *U = RgM_Rg_div(RgM_mul(*U,V),D);
1950 77 : if (remove < 2) *U = shallowconcat(K,*U);
1951 77 : if (!remove) x = shallowconcat(zeromatcopy(m-1,lg(K)-1), x);
1952 77 : gerepileall(av, 2, &x, U);
1953 : }
1954 1057 : else x = gerepilecopy(av,x);
1955 1134 : return x;
1956 : }
1957 : GEN
1958 1085 : QM_ImQ_hnfall(GEN x, GEN *U, long remove) { return QM_ImQ_all(x, U, remove, 1); }
1959 : GEN
1960 1008 : QM_ImQ_hnf(GEN x) { return QM_ImQ_hnfall(x, NULL, 1); }
1961 : GEN
1962 56 : QM_ImQ(GEN x) { return QM_ImQ_all(x, NULL, 1, 0); }
1963 :
1964 : GEN
1965 3752 : intersect(GEN x, GEN y)
1966 : {
1967 3752 : long j, lx = lg(x);
1968 : pari_sp av;
1969 : GEN z;
1970 :
1971 3752 : if (typ(x)!=t_MAT) pari_err_TYPE("intersect",x);
1972 3752 : if (typ(y)!=t_MAT) pari_err_TYPE("intersect",y);
1973 3752 : if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
1974 :
1975 3752 : av = avma; z = ker(shallowconcat(x,y));
1976 16324 : for (j=lg(z)-1; j; j--) setlg(z[j], lx);
1977 3752 : return gerepileupto(av, image(RgM_mul(x,z)));
1978 : }
|
