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