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