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