Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - alglin2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30347-cb65b7994e) Lines: 1075 1189 90.4 %
Date: 2025-06-27 09:22:08 Functions: 89 96 92.7 %
Legend: Lines: hit not hit

          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       34573 : Flm_charpoly_i(GEN x, ulong p)
      34             : {
      35       34573 :   long lx = lg(x), r, i;
      36       34573 :   GEN H, y = cgetg(lx+1, t_VEC);
      37       34573 :   gel(y,1) = pol1_Flx(0); H = Flm_hess(x, p);
      38      255413 :   for (r = 1; r < lx; r++)
      39             :   {
      40      220840 :     pari_sp av2 = avma;
      41      220840 :     ulong a = 1;
      42      220840 :     GEN z, b = zero_Flx(0);
      43      659278 :     for (i = r-1; i; i--)
      44             :     {
      45      550986 :       a = Fl_mul(a, ucoeff(H,i+1,i), p);
      46      550968 :       if (!a) break;
      47      438391 :       b = Flx_add(b, Flx_Fl_mul(gel(y,i), Fl_mul(a,ucoeff(H,i,r),p), p), p);
      48             :     }
      49      220840 :     z = Flx_sub(Flx_shift(gel(y,r), 1),
      50      220869 :                 Flx_Fl_mul(gel(y,r), ucoeff(H,r,r), p), p);
      51             :     /* (X - H[r,r])y[r] - b */
      52      220837 :     gel(y,r+1) = gc_leaf(av2, Flx_sub(z, b, p));
      53             :   }
      54       34573 :   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       29584 : FpM_charpoly(GEN x, GEN p)
      66             : {
      67       29584 :   pari_sp av = avma;
      68             :   long lx, r, i;
      69             :   GEN y, H;
      70             : 
      71       29584 :   if (lgefint(p) == 3)
      72             :   {
      73       28782 :     ulong pp = p[2];
      74       28782 :     y = Flx_to_ZX(Flm_charpoly_i(ZM_to_Flm(x,pp), pp));
      75       28782 :     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(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        5954 : easychar(GEN x, long v)
     128             : {
     129             :   pari_sp av;
     130             :   long lx;
     131             :   GEN p1;
     132             : 
     133        5954 :   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) = 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        5744 :     case t_MAT:
     167        5744 :       lx=lg(x);
     168        5744 :       if (lx==1) return pol_1(v);
     169        5688 :       if (lgcols(x) != lx) break;
     170        5681 :       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        3868 : charpoly(GEN x, long v)
     192             : {
     193        3868 :   GEN T, p = NULL;
     194             :   long prec;
     195        3868 :   if ((T = easychar(x,v))) return T;
     196        3637 :   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         147 :     case t_REAL:
     210             :     case t_COMPLEX:
     211             :     case t_PADIC:
     212         147 :       T = carhess(x, v);
     213         147 :       break;
     214        1050 :     default:
     215        1050 :       T = carberkowitz(x, v);
     216             :   }
     217        3637 :   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        1974 : fix_pol(GEN p, long v)
     223             : {
     224        1974 :   long w = gvar2(p);
     225        1974 :   if (w != NO_VARIABLE && varncmp(w, v) <= 0)
     226          56 :     p = poleval(p, pol_x(v));
     227             :   else
     228        1918 :     setvarn(p, v);
     229        1974 :   (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      219685 : FpM_trace(GEN x, GEN p)
     287             : {
     288      219685 :   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 from , 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        3249 : RgXQ_minpoly_fast(GEN x, GEN y)
     774             : {
     775             :   GEN p, pol;
     776             :   long pa;
     777        3249 :   long t = RgX_type2(x,y, &p,&pol,&pa);
     778        3249 :   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        3151 :     default:       return NULL;
     785             :   }
     786             : }
     787             : 
     788             : /* return caract(Mod(x,T)) in variable v */
     789             : GEN
     790        3249 : RgXQ_minpoly(GEN x, GEN T, long v)
     791             : {
     792        3249 :   pari_sp av = avma;
     793        3249 :   GEN R = RgXQ_minpoly_fast(x,T);
     794        3249 :   if (R) { setvarn(R, v); return R; }
     795        3151 :   if (!issquarefree(T))
     796             :   {
     797          63 :     R = RgXQ_minpoly_naive(x, T);
     798          63 :     setvarn(R,v); return R;
     799             :   }
     800        3088 :   R = RgXQ_charpoly(x, T, v);
     801        3088 :   return gc_upto(av, RgX_div(R,RgX_gcd(R, RgX_deriv(R))));
     802             : }
     803             : 
     804             : static GEN
     805        3668 : easymin(GEN x, long v)
     806             : {
     807             :   GEN G, R, dR;
     808        3668 :   long tx = typ(x);
     809        3668 :   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        3514 :   if (tx == t_POLMOD)
     815             :   {
     816        3465 :     GEN a = gel(x,2), b = gel(x,1);
     817        3465 :     if (degpol(b)==0) return pol_1(v);
     818        3437 :     if (typ(a) != t_POL || varn(a) != varn(b))
     819             :     {
     820         223 :       if (varncmp(gvar(a), v) <= 0) pari_err_PRIORITY("minpoly", x, "<", v);
     821         216 :       return deg1pol(gen_1, gneg_i(a), v);
     822             :     }
     823        3214 :     return RgXQ_minpoly(a, b, v);
     824             :   }
     825          49 :   R = easychar(x, v); if (!R) return NULL;
     826           0 :   dR = RgX_deriv(R);  if (!lgpol(dR)) return NULL;
     827           0 :   G = RgX_normalize(RgX_gcd(R,dR));
     828           0 :   return RgX_div(R,G);
     829             : }
     830             : GEN
     831        3668 : minpoly(GEN x, long v)
     832             : {
     833        3668 :   pari_sp av = avma;
     834             :   GEN P;
     835        3668 :   if (v < 0) v = 0;
     836        3668 :   P = easymin(x,v);
     837        3661 :   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             : static int
     848         364 : relative0(GEN a, GEN a0, long bit)
     849             : {
     850         364 :   if (gequal0(a)) return 1;
     851         343 :   if (gequal0(a0)) return gexpo(a) < bit;
     852         203 :   return (gexpo(a)-gexpo(a0) < bit);
     853             : }
     854             : /* x0 a nonempty square t_MAT that can be written to */
     855             : static GEN
     856         168 : RgM_hess(GEN x0, long prec)
     857             : {
     858         168 :   pari_sp av = avma;
     859         168 :   long lx = lg(x0), bit = prec? 8 - prec2nbits(prec): 0, m, i, j;
     860         168 :   GEN x = bit? RgM_shallowcopy(x0): x0;
     861             : 
     862         665 :   for (m=2; m<lx-1; m++)
     863             :   {
     864         497 :     GEN t = NULL;
     865         497 :     if (!bit)
     866             :     { /* first non-zero pivot */
     867          84 :       for (i=m+1; i<lx; i++)
     868          77 :         if (!gequal0(t = gcoeff(x,i,m-1))) break;
     869             :     }
     870             :     else
     871             :     { /* maximal pivot */
     872         434 :       long E = -(long)HIGHEXPOBIT, k = lx;
     873        3906 :       for (i=m+1; i<lx; i++)
     874             :       {
     875        3472 :         long e = gexpo(gcoeff(x,i,m-1));
     876        3472 :         if (e > E) { E = e; k = i; t = gcoeff(x,i,m-1); }
     877             :       }
     878         434 :       if (k != lx && relative0(t, gcoeff(x0,k,m-1), bit)) k = lx;
     879         434 :       i = k;
     880             :     }
     881         497 :     if (i == lx) continue;
     882        4438 :     for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
     883         385 :     swap(gel(x,i), gel(x,m));
     884         385 :     if (bit)
     885             :     {
     886        4102 :       for (j=m-1; j<lx; j++) swap(gcoeff(x0,i,j), gcoeff(x0,m,j));
     887         329 :       swap(gel(x0,i), gel(x0,m));
     888             :     }
     889         385 :     t = ginv(t);
     890             : 
     891        3668 :     for (i=m+1; i<lx; i++)
     892             :     {
     893        3283 :       GEN c = gcoeff(x,i,m-1);
     894        3283 :       if (gequal0(c)) continue;
     895             : 
     896        2520 :       c = gmul(c,t); gcoeff(x,i,m-1) = gen_0;
     897       41062 :       for (j=m; j<lx; j++)
     898       38542 :         gcoeff(x,i,j) = gsub(gcoeff(x,i,j), gmul(c,gcoeff(x,m,j)));
     899       67417 :       for (j=1; j<lx; j++)
     900       64897 :         gcoeff(x,j,m) = gadd(gcoeff(x,j,m), gmul(c,gcoeff(x,j,i)));
     901        2520 :       if (gc_needed(av,2))
     902             :       {
     903           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
     904           0 :         (void)gc_all(av,2, &x, &t);
     905             :       }
     906             :     }
     907             :   }
     908         168 :   return x;
     909             : }
     910             : 
     911             : GEN
     912         168 : hess(GEN x)
     913             : {
     914         168 :   pari_sp av = avma;
     915         168 :   GEN p = NULL, T = NULL;
     916         168 :   long prec, lx = lg(x);
     917         168 :   if (typ(x) != t_MAT) pari_err_TYPE("hess",x);
     918         168 :   if (lx == 1) return cgetg(1,t_MAT);
     919         168 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     920         168 :   switch(RgM_type(x, &p, &T, &prec))
     921             :   {
     922         140 :     case t_REAL:
     923         140 :     case t_COMPLEX: break;
     924          28 :     default: prec = 0;
     925             :   }
     926         168 :   return gc_GEN(av, RgM_hess(RgM_shallowcopy(x),prec));
     927             : }
     928             : 
     929             : GEN
     930       34573 : Flm_hess_pre(GEN x, ulong p, ulong pi)
     931             : {
     932       34573 :   long lx = lg(x), m, i, j;
     933       34573 :   if (lx == 1) return cgetg(1,t_MAT);
     934       34573 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     935             : 
     936       34573 :   x = Flm_copy(x);
     937      186729 :   for (m=2; m<lx-1; m++)
     938             :   {
     939      152316 :     ulong t = 0;
     940      718370 :     for (i=m+1; i<lx; i++) { t = ucoeff(x,i,m-1); if (t) break; }
     941      152316 :     if (i == lx) continue;
     942     1249223 :     for (j=m-1; j<lx; j++) lswap(ucoeff(x,i,j), ucoeff(x,m,j));
     943       91216 :     swap(gel(x,i), gel(x,m)); t = Fl_inv(t, p);
     944             : 
     945     1066356 :     for (i=m+1; i<lx; i++)
     946             :     {
     947      975300 :       ulong c = ucoeff(x,i,m-1);
     948      975300 :       if (!c) continue;
     949      344814 :       if (pi)
     950             :       {
     951      124791 :         c = Fl_mul_pre(c,t,p,pi); ucoeff(x,i,m-1) = 0;
     952     2015485 :         for (j=m; j<lx; j++)
     953     1890739 :           ucoeff(x,i,j) = Fl_sub(ucoeff(x,i,j), Fl_mul_pre(c,ucoeff(x,m,j), p, pi), p);
     954     3133794 :         for (j=1; j<lx; j++)
     955     3009476 :           ucoeff(x,j,m) = Fl_add(ucoeff(x,j,m), Fl_mul_pre(c,ucoeff(x,j,i), p, pi), p);
     956             :       } else
     957             :       {
     958      220023 :         c = Fl_mul(c,t,p); ucoeff(x,i,m-1) = 0;
     959     4327014 :         for (j=m; j<lx; j++)
     960     4106990 :           ucoeff(x,i,j) = Fl_sub(ucoeff(x,i,j), Fl_mul(c,ucoeff(x,m,j), p), p);
     961     7378743 :         for (j=1; j<lx; j++)
     962     7158414 :           ucoeff(x,j,m) = Fl_add(ucoeff(x,j,m), Fl_mul(c,ucoeff(x,j,i), p), p);
     963             :       }
     964             :     }
     965             :   }
     966       34413 :   return x;
     967             : }
     968             : 
     969             : GEN
     970       34573 : Flm_hess(GEN x, ulong p)
     971       34573 : { return Flm_hess_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
     972             : 
     973             : GEN
     974         802 : FpM_hess(GEN x, GEN p)
     975             : {
     976         802 :   pari_sp av = avma;
     977         802 :   long lx = lg(x), m, i, j;
     978         802 :   if (lx == 1) return cgetg(1,t_MAT);
     979         802 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     980         802 :   if (lgefint(p) == 3)
     981             :   {
     982           0 :     ulong pp = p[2];
     983           0 :     x = Flm_hess(ZM_to_Flm(x, pp), pp);
     984           0 :     return gc_upto(av, Flm_to_ZM(x));
     985             :   }
     986         802 :   x = RgM_shallowcopy(x);
     987        3392 :   for (m=2; m<lx-1; m++)
     988             :   {
     989        2590 :     GEN t = NULL;
     990        5509 :     for (i=m+1; i<lx; i++) { t = gcoeff(x,i,m-1); if (signe(t)) break; }
     991        2590 :     if (i == lx) continue;
     992       12019 :     for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
     993        1890 :     swap(gel(x,i), gel(x,m)); t = Fp_inv(t, p);
     994             : 
     995        8239 :     for (i=m+1; i<lx; i++)
     996             :     {
     997        6349 :       GEN c = gcoeff(x,i,m-1);
     998        6349 :       if (!signe(c)) continue;
     999             : 
    1000        5061 :       c = Fp_mul(c,t, p); gcoeff(x,i,m-1) = gen_0;
    1001       29785 :       for (j=m; j<lx; j++)
    1002       24724 :         gcoeff(x,i,j) = Fp_sub(gcoeff(x,i,j), Fp_mul(c,gcoeff(x,m,j),p), p);
    1003       45465 :       for (j=1; j<lx; j++)
    1004       40404 :         gcoeff(x,j,m) = Fp_add(gcoeff(x,j,m), Fp_mul(c,gcoeff(x,j,i),p), p);
    1005        5061 :       if (gc_needed(av,2))
    1006             :       {
    1007           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
    1008           0 :         (void)gc_all(av,2, &x, &t);
    1009             :       }
    1010             :     }
    1011             :   }
    1012         802 :   return gc_GEN(av,x);
    1013             : }
    1014             : /* H in Hessenberg form */
    1015             : static GEN
    1016         161 : RgM_hess_charpoly(GEN H, long v)
    1017             : {
    1018         161 :   long n = lg(H), r, i;
    1019         161 :   GEN y = cgetg(n+1, t_VEC);
    1020         161 :   gel(y,1) = pol_1(v);
    1021         945 :   for (r = 1; r < n; r++)
    1022             :   {
    1023         784 :     pari_sp av2 = avma;
    1024         784 :     GEN z, a = gen_1, b = pol_0(v);
    1025        4690 :     for (i = r-1; i; i--)
    1026             :     {
    1027        3983 :       a = gmul(a, gcoeff(H,i+1,i));
    1028        3983 :       if (gequal0(a)) break;
    1029        3906 :       b = RgX_add(b, RgX_Rg_mul(gel(y,i), gmul(a,gcoeff(H,i,r))));
    1030             :     }
    1031         784 :     z = RgX_sub(RgX_shift_shallow(gel(y,r), 1),
    1032         784 :                 RgX_Rg_mul(gel(y,r), gcoeff(H,r,r)));
    1033         784 :     gel(y,r+1) = gc_upto(av2, RgX_sub(z, b)); /* (X - H[r,r])y[r] - b */
    1034             :   }
    1035         161 :   return gel(y,n);
    1036             : }
    1037             : GEN
    1038         161 : carhess(GEN x, long v)
    1039             : {
    1040             :   pari_sp av;
    1041             :   GEN y;
    1042         161 :   if ((y = easychar(x,v))) return y;
    1043         161 :   av = avma; y = RgM_hess_charpoly(hess(x), fetch_var_higher());
    1044         161 :   return gc_upto(av, fix_pol(y, v));
    1045             : }
    1046             : 
    1047             : /* Bound for sup norm of charpoly(M/dM), M integral: let B = |M|oo / |dM|,
    1048             :  *   s = max_k binomial(n,k) (kB^2)^(k/2),
    1049             :  * return ceil(log2(s)) */
    1050             : static long
    1051        3938 : charpoly_bound(GEN M, GEN dM, GEN N)
    1052             : {
    1053        3938 :   pari_sp av = avma;
    1054        3938 :   GEN B = itor(N, LOWDEFAULTPREC);
    1055        3938 :   GEN s = real_0(LOWDEFAULTPREC), bin, B2;
    1056        3938 :   long n = lg(M)-1, k;
    1057        3938 :   bin = gen_1;
    1058        3938 :   if (dM) B = divri(B, dM);
    1059        3938 :   B2 = sqrr(B);
    1060       17845 :   for (k = n; k >= (n+1)>>1; k--)
    1061             :   {
    1062       13907 :     GEN t = mulri(powruhalf(mulur(k, B2), k), bin);
    1063       13907 :     if (abscmprr(t, s) > 0) s = t;
    1064       13907 :     bin = diviuexact(muliu(bin, k), n-k+1);
    1065             :   }
    1066        3938 :   return gc_long(av, ceil(dbllog2(s)));
    1067             : }
    1068             : 
    1069             : static GEN
    1070        4355 : QM_charpoly_ZX_slice(GEN A, GEN dM, GEN P, GEN *mod)
    1071             : {
    1072        4355 :   pari_sp av = avma;
    1073        4355 :   long i, n = lg(P)-1;
    1074             :   GEN H, T;
    1075        4355 :   if (n == 1)
    1076             :   {
    1077        3613 :     ulong p = uel(P,1), dp = dM ? umodiu(dM, p): 1;
    1078        3613 :     GEN Hp, a = ZM_to_Flm(A, p);
    1079        3613 :     Hp = Flm_charpoly_i(a, p);
    1080        3613 :     if (dp != 1) Hp = Flx_rescale(Hp, Fl_inv(dp, p), p);
    1081        3613 :     Hp = gc_upto(av, Flx_to_ZX(Hp));
    1082        3613 :     *mod = utoipos(p); return Hp;
    1083             :   }
    1084         742 :   T = ZV_producttree(P);
    1085         742 :   A = ZM_nv_mod_tree(A, P, T);
    1086         742 :   H = cgetg(n+1, t_VEC);
    1087        2899 :   for(i=1; i <= n; i++)
    1088             :   {
    1089        2157 :     ulong p = uel(P,i), dp = dM ? umodiu(dM, p): 1;
    1090        2157 :     gel(H,i) = Flm_charpoly(gel(A, i), p);
    1091        2157 :     if (dp != 1) gel(H,i) = Flx_rescale(gel(H,i), Fl_inv(dp, p), p);
    1092             :   }
    1093         742 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P,T));
    1094         742 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    1095             : }
    1096             : 
    1097             : GEN
    1098        4355 : QM_charpoly_ZX_worker(GEN P, GEN M, GEN dM)
    1099             : {
    1100        4355 :   GEN V = cgetg(3, t_VEC);
    1101        4355 :   gel(V,1) = QM_charpoly_ZX_slice(M, equali1(dM) ? NULL:dM, P, &gel(V,2));
    1102        4355 :   return V;
    1103             : }
    1104             : 
    1105             : /* Assume M a square ZM, dM integer. Return charpoly(M / dM) in Z[X] */
    1106             : static GEN
    1107        4526 : QM_charpoly_ZX_i(GEN M, GEN dM, long bound)
    1108             : {
    1109        4526 :   long n = lg(M)-1;
    1110             :   forprime_t S;
    1111        4526 :   GEN worker = snm_closure(is_entry("_QM_charpoly_ZX_worker"),
    1112             :                            mkvec2(M, dM? dM: gen_1));
    1113        4526 :   if (!n) return pol_1(0);
    1114        4526 :   if (bound < 0)
    1115             :   {
    1116        4260 :     GEN N = ZM_supnorm(M);
    1117        4260 :     if (signe(N) == 0) return monomial(gen_1, n, 0);
    1118        3938 :     bound = charpoly_bound(M, dM, N) + 1;
    1119             :   }
    1120        4204 :   if (DEBUGLEVEL>5) err_printf("ZM_charpoly: bound 2^%ld\n", bound);
    1121        4204 :   init_modular_big(&S);
    1122        4204 :   return gen_crt("QM_charpoly_ZX", worker, &S, dM, bound, 0, NULL,
    1123             :                  nxV_chinese_center, FpX_center);
    1124             : }
    1125             : 
    1126             : GEN
    1127         266 : QM_charpoly_ZX_bound(GEN M, long bit)
    1128             : {
    1129         266 :   pari_sp av = avma;
    1130         266 :   GEN dM; M = Q_remove_denom(M, &dM);
    1131         266 :   return gc_GEN(av, QM_charpoly_ZX_i(M, dM, bit));
    1132             : }
    1133             : GEN
    1134        1848 : QM_charpoly_ZX(GEN M)
    1135             : {
    1136        1848 :   pari_sp av = avma;
    1137        1848 :   GEN dM; M = Q_remove_denom(M, &dM);
    1138        1848 :   return gc_GEN(av, QM_charpoly_ZX_i(M, dM, -1));
    1139             : }
    1140             : GEN
    1141        2412 : ZM_charpoly(GEN M)
    1142             : {
    1143        2412 :   pari_sp av = avma;
    1144        2412 :   return gc_GEN(av, QM_charpoly_ZX_i(M, NULL, -1));
    1145             : }
    1146             : 
    1147             : GEN
    1148      289419 : ZM_trace(GEN x)
    1149             : {
    1150      289419 :   pari_sp av = avma;
    1151      289419 :   long i, lx = lg(x);
    1152             :   GEN t;
    1153      289419 :   if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
    1154      276341 :   t = gcoeff(x,1,1);
    1155     2857971 :   for (i = 2; i < lx; i++) t = addii(t, gcoeff(x,i,i));
    1156      276341 :   return gc_INT(av, t);
    1157             : }
    1158             : 
    1159             : /*******************************************************************/
    1160             : /*                                                                 */
    1161             : /*        CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)        */
    1162             : /*                                                                 */
    1163             : /*******************************************************************/
    1164             : GEN
    1165        1743 : carberkowitz(GEN x, long v)
    1166             : {
    1167             :   long lx, i, j, k, r;
    1168             :   GEN V, S, C, Q;
    1169             :   pari_sp av0, av;
    1170        1743 :   if ((V = easychar(x,v))) return V;
    1171        1743 :   lx = lg(x); av0 = avma;
    1172        1743 :   V = cgetg(lx+1, t_VEC);
    1173        1743 :   S = cgetg(lx+1, t_VEC);
    1174        1743 :   C = cgetg(lx+1, t_VEC);
    1175        1743 :   Q = cgetg(lx+1, t_VEC);
    1176        1743 :   av = avma;
    1177        1743 :   gel(C,1) = gen_m1;
    1178        1743 :   gel(V,1) = gen_m1;
    1179       12306 :   for (i=2;i<=lx; i++) gel(C,i) = gel(Q,i) = gel(S,i) = gel(V,i) = gen_0;
    1180        1743 :   gel(V,2) = gcoeff(x,1,1);
    1181       10563 :   for (r = 2; r < lx; r++)
    1182             :   {
    1183             :     pari_sp av2;
    1184             :     GEN t;
    1185             : 
    1186       68068 :     for (i = 1; i < r; i++) gel(S,i) = gcoeff(x,i,r);
    1187        8820 :     gel(C,2) = gcoeff(x,r,r);
    1188       59248 :     for (i = 1; i < r-1; i++)
    1189             :     {
    1190       50428 :       av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1191      726376 :       for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1192       50428 :       gel(C,i+2) = gc_upto(av2, t);
    1193      776804 :       for (j = 1; j < r; j++)
    1194             :       {
    1195      726376 :         av2 = avma; t = gmul(gcoeff(x,j,1), gel(S,1));
    1196    14286608 :         for (k = 2; k < r; k++) t = gadd(t, gmul(gcoeff(x,j,k), gel(S,k)));
    1197      726376 :         gel(Q,j) = gc_upto(av2, t);
    1198             :       }
    1199      776804 :       for (j = 1; j < r; j++) gel(S,j) = gel(Q,j);
    1200             :     }
    1201        8820 :     av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1202       59248 :     for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1203        8820 :     gel(C,r+1) = gc_upto(av2, t);
    1204        8820 :     if (gc_needed(av0,1))
    1205             :     {
    1206           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"carberkowitz");
    1207           0 :       (void)gc_all(av, 2, &C, &V);
    1208             :     }
    1209       85708 :     for (i = 1; i <= r+1; i++)
    1210             :     {
    1211       76888 :       av2 = avma; t = gmul(gel(C,i), gel(V,1));
    1212      558572 :       for (j = 2; j <= minss(r,i); j++)
    1213      481684 :         t = gadd(t, gmul(gel(C,i+1-j), gel(V,j)));
    1214       76888 :       gel(Q,i) = gc_upto(av2, t);
    1215             :     }
    1216       85708 :     for (i = 1; i <= r+1; i++) gel(V,i) = gel(Q,i);
    1217             :   }
    1218        1743 :   V = gtopoly(V, fetch_var_higher());
    1219        1743 :   if (!odd(lx)) V = RgX_neg(V);
    1220        1743 :   return gc_upto(av0, fix_pol(V, v));
    1221             : }
    1222             : 
    1223             : /*******************************************************************/
    1224             : /*                                                                 */
    1225             : /*                            NORMS                                */
    1226             : /*                                                                 */
    1227             : /*******************************************************************/
    1228             : GEN
    1229     3540275 : gnorm(GEN x)
    1230             : {
    1231             :   pari_sp av;
    1232     3540275 :   switch(typ(x))
    1233             :   {
    1234       50856 :     case t_INT:  return sqri(x);
    1235      537528 :     case t_REAL: return sqrr(x);
    1236        1393 :     case t_FRAC: return sqrfrac(x);
    1237     2879700 :     case t_COMPLEX: av = avma; return gc_upto(av, cxnorm(x));
    1238       69041 :     case t_QUAD:    av = avma; return gc_upto(av, quadnorm(x));
    1239             : 
    1240          14 :     case t_POL: case t_SER: case t_RFRAC: av = avma;
    1241          14 :       return gc_upto(av, greal(gmul(conj_i(x),x)));
    1242             : 
    1243          28 :     case t_FFELT:
    1244          28 :       retmkintmod(FF_norm(x), FF_p(x));
    1245             : 
    1246        1715 :     case t_POLMOD:
    1247             :     {
    1248        1715 :       GEN T = gel(x,1), a = gel(x,2);
    1249        1715 :       if (typ(a) != t_POL || varn(a) != varn(T)) return gpowgs(a, degpol(T));
    1250        1533 :       return RgXQ_norm(a, T);
    1251             :     }
    1252           0 :     case t_VEC: case t_COL: case t_MAT:
    1253           0 :       pari_APPLY_same(gnorm(gel(x,i)));
    1254             :   }
    1255           0 :   pari_err_TYPE("gnorm",x);
    1256             :   return NULL; /* LCOV_EXCL_LINE */
    1257             : }
    1258             : 
    1259             : /* return |q|^2, complex modulus */
    1260             : static GEN
    1261          28 : cxquadnorm(GEN q, long prec)
    1262             : {
    1263          28 :   GEN X = gel(q,1), c = gel(X,2); /* (1-D)/4, -D/4 */
    1264          28 :   if (signe(c) > 0) return quadnorm(q); /* imaginary */
    1265          21 :   if (!prec) pari_err_TYPE("gnorml2", q);
    1266           7 :   return sqrr(quadtofp(q, prec));
    1267             : }
    1268             : 
    1269             : static GEN
    1270    37213970 : gnorml2_i(GEN x, long prec)
    1271             : {
    1272             :   pari_sp av;
    1273             :   long i, lx;
    1274             :   GEN s;
    1275             : 
    1276    37213970 :   switch(typ(x))
    1277             :   {
    1278     1353634 :     case t_INT:  return sqri(x);
    1279    27049187 :     case t_REAL: return sqrr(x);
    1280           7 :     case t_FRAC: return sqrfrac(x);
    1281     4150822 :     case t_COMPLEX: av = avma; return gc_upto(av, cxnorm(x));
    1282          21 :     case t_QUAD:    av = avma; return gc_upto(av, cxquadnorm(x,prec));
    1283             : 
    1284       61614 :     case t_POL: lx = lg(x)-1; x++; break;
    1285             : 
    1286     4601385 :     case t_VEC:
    1287             :     case t_COL:
    1288     4601385 :     case t_MAT: lx = lg(x); break;
    1289             : 
    1290           0 :     default: pari_err_TYPE("gnorml2",x);
    1291             :       return NULL; /* LCOV_EXCL_LINE */
    1292             :   }
    1293     4662999 :   if (lx == 1) return gen_0;
    1294     4662999 :   av = avma;
    1295     4662999 :   s = gnorml2(gel(x,1));
    1296    32665647 :   for (i=2; i<lx; i++)
    1297             :   {
    1298    28003389 :     s = gadd(s, gnorml2(gel(x,i)));
    1299    28002817 :     if (gc_needed(av,1))
    1300             :     {
    1301           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnorml2");
    1302           0 :       s = gc_upto(av, s);
    1303             :     }
    1304             :   }
    1305     4662258 :   return gc_upto(av,s);
    1306             : }
    1307             : GEN
    1308    37213303 : gnorml2(GEN x) { return gnorml2_i(x, 0); }
    1309             : 
    1310             : static GEN pnormlp(GEN,GEN,long);
    1311             : static GEN
    1312          63 : pnormlpvec(long i0, GEN x, GEN p, long prec)
    1313             : {
    1314          63 :   pari_sp av = avma;
    1315          63 :   long i, lx = lg(x);
    1316          63 :   GEN s = gen_0;
    1317         224 :   for (i=i0; i<lx; i++)
    1318             :   {
    1319         161 :     s = gadd(s, pnormlp(gel(x,i),p,prec));
    1320         161 :     if (gc_needed(av,1))
    1321             :     {
    1322           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnormlp, i = %ld", i);
    1323           0 :       s = gc_upto(av, s);
    1324             :     }
    1325             :   }
    1326          63 :   return s;
    1327             : }
    1328             : /* (||x||_p)^p */
    1329             : static GEN
    1330         196 : pnormlp(GEN x, GEN p, long prec)
    1331             : {
    1332         196 :   switch(typ(x))
    1333             :   {
    1334         119 :     case t_INT: case t_REAL: x = mpabs(x); break;
    1335           0 :     case t_FRAC: x = absfrac(x); break;
    1336          14 :     case t_COMPLEX: case t_QUAD: x = gabs(x,prec); break;
    1337           7 :     case t_POL: return pnormlpvec(2, x, p, prec);
    1338          56 :     case t_VEC: case t_COL: case t_MAT: return pnormlpvec(1, x, p, prec);
    1339           0 :     default: pari_err_TYPE("gnormlp",x);
    1340             :   }
    1341         133 :   return gpow(x, p, prec);
    1342             : }
    1343             : 
    1344             : GEN
    1345         371 : gnormlp(GEN x, GEN p, long prec)
    1346             : {
    1347         371 :   pari_sp av = avma;
    1348         371 :   if (!p || (typ(p) == t_INFINITY && inf_get_sign(p) > 0))
    1349         210 :     return gsupnorm(x, prec);
    1350         161 :   if (gsigne(p) <= 0) pari_err_DOMAIN("normlp", "p", "<=", gen_0, p);
    1351         154 :   if (is_scalar_t(typ(x))) return gabs(x, prec);
    1352          91 :   if (typ(p) == t_INT)
    1353             :   {
    1354          63 :     ulong pp = itou_or_0(p);
    1355          63 :     switch(pp)
    1356             :     {
    1357          28 :       case 1: return gnorml1(x, prec);
    1358          28 :       case 2: x = gnorml2_i(x, prec); break;
    1359           7 :       default: x = pnormlp(x, p, prec); break;
    1360             :     }
    1361          35 :     if (pp && typ(x) == t_INT && Z_ispowerall(x, pp, &x))
    1362           7 :       return gc_leaf(av, x);
    1363          28 :     if (pp == 2) return gc_upto(av, gsqrt(x, prec));
    1364             :   }
    1365             :   else
    1366          28 :     x = pnormlp(x, p, prec);
    1367          28 :   x = gpow(x, ginv(p), prec);
    1368          28 :   return gc_upto(av, x);
    1369             : }
    1370             : 
    1371             : GEN
    1372         168 : gnorml1(GEN x,long prec)
    1373             : {
    1374         168 :   pari_sp av = avma;
    1375             :   long lx,i;
    1376             :   GEN s;
    1377         168 :   switch(typ(x))
    1378             :   {
    1379          98 :     case t_INT: case t_REAL: return mpabs(x);
    1380           0 :     case t_FRAC: return absfrac(x);
    1381             : 
    1382          14 :     case t_COMPLEX: case t_QUAD:
    1383          14 :       return gabs(x,prec);
    1384             : 
    1385           7 :     case t_POL:
    1386           7 :       lx = lg(x); s = gen_0;
    1387          28 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1388           7 :       break;
    1389             : 
    1390          49 :     case t_VEC: case t_COL: case t_MAT:
    1391          49 :       lx = lg(x); s = gen_0;
    1392         168 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1393          49 :       break;
    1394             : 
    1395           0 :     default: pari_err_TYPE("gnorml1",x);
    1396             :       return NULL; /* LCOV_EXCL_LINE */
    1397             :   }
    1398          56 :   return gc_upto(av, s);
    1399             : }
    1400             : /* As gnorml1, except for t_QUAD and t_COMPLEX: |x + wy| := |x| + |y|
    1401             :  * Still a norm of R-vector spaces, and can be cheaply computed without
    1402             :  * square roots */
    1403             : GEN
    1404      152173 : gnorml1_fake(GEN x)
    1405             : {
    1406      152173 :   pari_sp av = avma;
    1407             :   long lx, i;
    1408             :   GEN s;
    1409      152173 :   switch(typ(x))
    1410             :   {
    1411      133378 :     case t_INT: case t_REAL: return mpabs(x);
    1412           0 :     case t_FRAC: return absfrac(x);
    1413             : 
    1414           0 :     case t_COMPLEX:
    1415           0 :       s = gadd(gnorml1_fake(gel(x,1)), gnorml1_fake(gel(x,2)));
    1416           0 :       break;
    1417           0 :     case t_QUAD:
    1418           0 :       s = gadd(gnorml1_fake(gel(x,2)), gnorml1_fake(gel(x,3)));
    1419           0 :       break;
    1420             : 
    1421       18795 :     case t_POL:
    1422       18795 :       lx = lg(x); s = gen_0;
    1423      110789 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1424       18795 :       break;
    1425             : 
    1426           0 :     case t_VEC: case t_COL: case t_MAT:
    1427           0 :       lx = lg(x); s = gen_0;
    1428           0 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1429           0 :       break;
    1430             : 
    1431           0 :     default: pari_err_TYPE("gnorml1_fake",x);
    1432             :       return NULL; /* LCOV_EXCL_LINE */
    1433             :   }
    1434       18795 :   return gc_upto(av, s);
    1435             : }
    1436             : 
    1437             : static void
    1438    29518930 : store(GEN z, GEN *m) { if (!*m || gcmp(z, *m) > 0) *m = z; }
    1439             : /* compare |x| to *m or |x|^2 to *msq, whichever is easiest, and update
    1440             :  * the pointed value if x is larger */
    1441             : void
    1442    35518196 : gsupnorm_aux(GEN x, GEN *m, GEN *msq, long prec)
    1443             : {
    1444             :   long i, lx;
    1445             :   GEN z;
    1446    35518196 :   switch(typ(x))
    1447             :   {
    1448       92370 :     case t_COMPLEX: z = cxnorm(x); store(z, msq); return;
    1449           7 :     case t_QUAD:  z = cxquadnorm(x,prec); store(z, msq); return;
    1450    29426605 :     case t_INT: case t_REAL: z = mpabs(x); store(z,m); return;
    1451          35 :     case t_FRAC: z = absfrac(x); store(z,m); return;
    1452          14 :     case t_INFINITY: store(mkoo(), m);
    1453             : 
    1454       81322 :     case t_POL: lx = lg(x)-1; x++; break;
    1455             : 
    1456     5917934 :     case t_VEC:
    1457             :     case t_COL:
    1458     5917934 :     case t_MAT: lx = lg(x); break;
    1459             : 
    1460           0 :     default: pari_err_TYPE("gsupnorm",x);
    1461             :       return; /* LCOV_EXCL_LINE */
    1462             :   }
    1463    40279346 :   for (i=1; i<lx; i++) gsupnorm_aux(gel(x,i), m, msq, prec);
    1464             : }
    1465             : GEN
    1466     1238043 : gsupnorm(GEN x, long prec)
    1467             : {
    1468     1238043 :   GEN m = NULL, msq = NULL;
    1469     1238043 :   pari_sp av = avma;
    1470     1238043 :   gsupnorm_aux(x, &m, &msq, prec);
    1471             :   /* now set m = max (m, sqrt(msq)) */
    1472     1238053 :   if (msq) {
    1473       15094 :     msq = gsqrt(msq, prec);
    1474       15093 :     if (!m || gcmp(m, msq) < 0) m = msq;
    1475     1222959 :   } else if (!m) m = gen_0;
    1476     1238052 :   return gc_GEN(av, m);
    1477             : }
    1478             : 
    1479             : /*******************************************************************/
    1480             : /*                                                                 */
    1481             : /*                            TRACES                               */
    1482             : /*                                                                 */
    1483             : /*******************************************************************/
    1484             : GEN
    1485          35 : matcompanion(GEN x)
    1486             : {
    1487          35 :   long n = degpol(x), j;
    1488             :   GEN y, c;
    1489             : 
    1490          35 :   if (typ(x)!=t_POL) pari_err_TYPE("matcompanion",x);
    1491          35 :   if (!signe(x)) pari_err_DOMAIN("matcompanion","polynomial","=",gen_0,x);
    1492          28 :   if (n == 0) return cgetg(1, t_MAT);
    1493             : 
    1494          28 :   y = cgetg(n+1,t_MAT);
    1495         105 :   for (j=1; j < n; j++) gel(y,j) = col_ei(n, j+1);
    1496          28 :   c = cgetg(n+1,t_COL); gel(y,n) = c;
    1497          28 :   if (gequal1(gel(x, n+2)))
    1498         112 :     for (j=1; j<=n; j++) gel(c,j) = gneg(gel(x,j+1));
    1499             :   else
    1500             :   { /* not monic. Hardly ever used */
    1501           7 :     pari_sp av = avma;
    1502           7 :     GEN d = gclone(gneg(gel(x,n+2)));
    1503           7 :     set_avma(av);
    1504          21 :     for (j=1; j<=n; j++) gel(c,j) = gdiv(gel(x,j+1), d);
    1505           7 :     gunclone(d);
    1506             :   }
    1507          28 :   return y;
    1508             : }
    1509             : 
    1510             : GEN
    1511      762291 : gtrace(GEN x)
    1512             : {
    1513             :   pari_sp av;
    1514      762291 :   long lx, tx = typ(x);
    1515             :   GEN y, z;
    1516             : 
    1517      762291 :   switch(tx)
    1518             :   {
    1519       23630 :     case t_INT: case t_REAL: case t_FRAC:
    1520       23630 :       return gmul2n(x,1);
    1521             : 
    1522      735808 :     case t_COMPLEX:
    1523      735808 :       return gmul2n(gel(x,1),1);
    1524             : 
    1525         154 :     case t_QUAD:
    1526         154 :       y = gel(x,1);
    1527         154 :       if (!gequal0(gel(y,3)))
    1528             :       { /* assume quad. polynomial is either x^2 + d or x^2 - x + d */
    1529         154 :         av = avma;
    1530         154 :         return gc_upto(av, gadd(gel(x,3), gmul2n(gel(x,2),1)));
    1531             :       }
    1532           0 :       return gmul2n(gel(x,2),1);
    1533             : 
    1534           7 :     case t_POL:
    1535          21 :       pari_APPLY_pol(gtrace(gel(x,i)));
    1536          14 :     case t_SER:
    1537          14 :       if (ser_isexactzero(x)) return gcopy(x);
    1538          21 :       pari_APPLY_ser(gtrace(gel(x,i)));
    1539             : 
    1540         784 :     case t_POLMOD:
    1541         784 :       y = gel(x,1); z = gel(x,2);
    1542         784 :       if (typ(z) != t_POL || varn(y) != varn(z)) return gmulsg(degpol(y), z);
    1543         476 :       av = avma;
    1544         476 :       return gc_upto(av, RgXQ_trace(z, y));
    1545             : 
    1546          28 :     case t_FFELT:
    1547          28 :       retmkintmod(FF_trace(x), FF_p(x));
    1548             : 
    1549           7 :     case t_RFRAC:
    1550           7 :       av = avma; return gc_upto(av, gadd(x, conj_i(x)));
    1551             : 
    1552           0 :     case t_VEC: case t_COL:
    1553           0 :       pari_APPLY_same(gtrace(gel(x,i)));
    1554             : 
    1555        1862 :     case t_MAT:
    1556        1862 :       lx = lg(x); if (lx == 1) return gen_0;
    1557             :       /*now lx >= 2*/
    1558        1855 :       if (lx != lgcols(x)) pari_err_DIM("gtrace");
    1559        1848 :       av = avma; return gc_upto(av, mattrace(x));
    1560             :   }
    1561           0 :   pari_err_TYPE("gtrace",x);
    1562             :   return NULL; /* LCOV_EXCL_LINE */
    1563             : }
    1564             : 
    1565             : /* Cholesky decomposition for positive definite matrix a
    1566             :  * [GTM138, Algo 2.7.6, matrix Q]
    1567             :  * If a is not positive definite return NULL. */
    1568             : GEN
    1569       70142 : qfgaussred_positive(GEN a)
    1570             : {
    1571       70142 :   pari_sp av = avma;
    1572             :   GEN b;
    1573       70142 :   long i,j,k, n = lg(a);
    1574             : 
    1575       70142 :   if (typ(a)!=t_MAT) pari_err_TYPE("qfgaussred_positive",a);
    1576       70142 :   if (n == 1) return cgetg(1, t_MAT);
    1577       70135 :   if (lgcols(a)!=n) pari_err_DIM("qfgaussred_positive");
    1578       70135 :   b = cgetg(n,t_MAT);
    1579      371149 :   for (j=1; j<n; j++)
    1580             :   {
    1581      301014 :     GEN p1=cgetg(n,t_COL), p2=gel(a,j);
    1582             : 
    1583      301014 :     gel(b,j) = p1;
    1584     1518194 :     for (i=1; i<=j; i++) gel(p1,i) = gel(p2,i);
    1585     1217180 :     for (   ; i<n ; i++) gel(p1,i) = gen_0;
    1586             :   }
    1587      365364 :   for (k=1; k<n; k++)
    1588             :   {
    1589      297664 :     GEN bk, p = gcoeff(b,k,k), invp;
    1590      297664 :     if (gsigne(p)<=0) return gc_NULL(av); /* not positive definite */
    1591      295230 :     invp = ginv(p);
    1592      295225 :     bk = row(b, k);
    1593     1196945 :     for (i=k+1; i<n; i++) gcoeff(b,k,i) = gmul(gel(bk,i), invp);
    1594     1196943 :     for (i=k+1; i<n; i++)
    1595             :     {
    1596      901715 :       GEN c = gel(bk, i);
    1597     5563873 :       for (j=i; j<n; j++)
    1598     4662160 :         gcoeff(b,i,j) = gsub(gcoeff(b,i,j), gmul(c,gcoeff(b,k,j)));
    1599             :     }
    1600      295228 :     if (gc_needed(av,1))
    1601             :     {
    1602           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfgaussred_positive");
    1603           0 :       b=gc_GEN(av,b);
    1604             :     }
    1605             :   }
    1606       67700 :   return gc_GEN(av,b);
    1607             : }
    1608             : 
    1609             : GEN
    1610       68786 : RgM_Cholesky(GEN M, long prec)
    1611             : {
    1612       68786 :   pari_sp av = avma;
    1613       68786 :   long i, j, lM = lg(M);
    1614       68786 :   GEN R, L = qfgaussred_positive(M);
    1615       68786 :   if (!L) return gc_NULL(av);
    1616      340201 :   R = cgetg(lM, t_MAT); for (j = 1; j < lM; j++) gel(R,j) = cgetg(lM, t_COL);
    1617      340198 :   for (i = 1; i < lM; i++)
    1618             :   {
    1619      273848 :     GEN r = gsqrt(gcoeff(L,i,i), prec);
    1620     2034357 :     for (j = 1; j < lM; j++)
    1621     1760510 :       gcoeff(R, i, j) = (i == j) ? r: gmul(r, gcoeff(L, i, j));
    1622             :   }
    1623       66350 :   return gc_upto(av, R);
    1624             : }
    1625             : 
    1626             : /* Maximal pivot strategy: x is a suitable pivot if it is non zero and either
    1627             :  * - an exact type, or
    1628             :  * - it is maximal among remaining nonzero (t_REAL) pivots */
    1629             : static int
    1630       47429 : suitable(GEN x, long k, GEN *pp, long *pi)
    1631             : {
    1632       47429 :   long t = typ(x);
    1633       47429 :   switch(t)
    1634             :   {
    1635       24785 :     case t_INT: return signe(x) != 0;
    1636       22490 :     case t_FRAC: return 1;
    1637         154 :     case t_REAL: {
    1638         154 :       GEN p = *pp;
    1639         154 :       if (signe(x) && (!p || abscmprr(p, x) < 0)) { *pp = x; *pi = k; }
    1640         154 :       return 0;
    1641             :     }
    1642           0 :     default: return !gequal0(x);
    1643             :   }
    1644             : }
    1645             : 
    1646             : /* Gauss reduction (arbitrary symmetric matrix, only the part above the
    1647             :  * diagonal is considered). If signature is nonzero, return only the
    1648             :  * signature, in which case gsigne() should be defined for elements of a. */
    1649             : static GEN
    1650       12200 : gaussred(GEN a, long signature)
    1651             : {
    1652             :   GEN r, ak, al;
    1653       12200 :   pari_sp av = avma, av1;
    1654       12200 :   long n = lg(a), i, j, k, l, sp, sn, t;
    1655             : 
    1656       12200 :   if (typ(a) != t_MAT) pari_err_TYPE("gaussred",a);
    1657       12200 :   if (n == 1) return signature? mkvec2(gen_0, gen_0): cgetg(1, t_MAT);
    1658       12200 :   if (lgcols(a) != n) pari_err_DIM("gaussred");
    1659       12200 :   n--;
    1660       12200 :   r = const_vecsmall(n, 1); av1= avma;
    1661       12200 :   a = RgM_shallowcopy(a);
    1662       12200 :   t = n; sp = sn = 0;
    1663       59461 :   while (t)
    1664             :   {
    1665       47261 :     long pind = 0;
    1666       47261 :     GEN invp, p = NULL;
    1667      130579 :     k=1; while (k<=n && (!r[k] || !suitable(gcoeff(a,k,k), k, &p, &pind))) k++;
    1668       47261 :     if (k > n && p) k = pind;
    1669       47261 :     if (k <= n)
    1670             :     {
    1671       47247 :       p = gcoeff(a,k,k); invp = ginv(p); /* != 0 */
    1672       47247 :       if (signature) { /* skip if (!signature): gsigne may fail ! */
    1673       47184 :         if (gsigne(p) > 0) sp++; else sn++;
    1674             :       }
    1675       47247 :       r[k] = 0; t--;
    1676       47247 :       ak = row(a, k);
    1677      260654 :       for (i=1; i<=n; i++)
    1678      213407 :         gcoeff(a,k,i) = r[i]? gmul(gcoeff(a,k,i), invp): gen_0;
    1679             : 
    1680      260654 :       for (i=1; i<=n; i++) if (r[i])
    1681             :       {
    1682       83052 :         GEN c = gel(ak,i); /* - p * a[k,i] */
    1683       83052 :         if (gequal0(c)) continue;
    1684      462136 :         for (j=1; j<=n; j++) if (r[j])
    1685      244348 :           gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
    1686             :       }
    1687       47247 :       gcoeff(a,k,k) = p;
    1688       47247 :       if (gc_needed(av1,1))
    1689             :       {
    1690           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"gaussred (t = %ld)", t);
    1691           0 :         a = gc_GEN(av1, a);
    1692             :       }
    1693             :     }
    1694             :     else
    1695             :     { /* all remaining diagonal coeffs are currently 0 */
    1696          14 :       for (k=1; k<=n; k++) if (r[k])
    1697             :       {
    1698          14 :         l=k+1; while (l<=n && (!r[l] || !suitable(gcoeff(a,k,l), l, &p, &pind))) l++;
    1699          14 :         if (l > n && p) l = pind;
    1700          14 :         if (l > n) continue;
    1701             : 
    1702          14 :         p = gcoeff(a,k,l); invp = ginv(p);
    1703          14 :         sp++; sn++;
    1704          14 :         r[k] = r[l] = 0; t -= 2;
    1705          14 :         ak = row(a, k);
    1706          14 :         al = row(a, l);
    1707          70 :         for (i=1; i<=n; i++) if (r[i])
    1708             :         {
    1709          28 :           gcoeff(a,k,i) = gmul(gcoeff(a,k,i), invp);
    1710          28 :           gcoeff(a,l,i) = gmul(gcoeff(a,l,i), invp);
    1711             :         } else {
    1712          28 :           gcoeff(a,k,i) = gen_0;
    1713          28 :           gcoeff(a,l,i) = gen_0;
    1714             :         }
    1715             : 
    1716          70 :         for (i=1; i<=n; i++) if (r[i])
    1717             :         { /* c = a[k,i] * p, d = a[l,i] * p; */
    1718          28 :           GEN c = gel(ak,i), d = gel(al,i);
    1719         140 :           for (j=1; j<=n; j++) if (r[j])
    1720          56 :             gcoeff(a,i,j) = gsub(gcoeff(a,i,j),
    1721          56 :                                  gadd(gmul(gcoeff(a,l,j), c),
    1722          56 :                                       gmul(gcoeff(a,k,j), d)));
    1723             :         }
    1724          70 :         for (i=1; i<=n; i++) if (r[i])
    1725             :         {
    1726          28 :           GEN c = gcoeff(a,k,i), d = gcoeff(a,l,i);
    1727          28 :           gcoeff(a,k,i) = gadd(c, d);
    1728          28 :           gcoeff(a,l,i) = gsub(c, d);
    1729             :         }
    1730          14 :         gcoeff(a,k,l) = gen_1;
    1731          14 :         gcoeff(a,l,k) = gen_m1;
    1732          14 :         gcoeff(a,k,k) = gmul2n(p,-1);
    1733          14 :         gcoeff(a,l,l) = gneg(gcoeff(a,k,k));
    1734          14 :         if (gc_needed(av1,1))
    1735             :         {
    1736           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gaussred");
    1737           0 :           a = gc_GEN(av1, a);
    1738             :         }
    1739          14 :         break;
    1740             :       }
    1741          14 :       if (k > n) break;
    1742             :     }
    1743             :   }
    1744       12200 :   if (!signature) return gc_GEN(av, a);
    1745       12179 :   set_avma(av); return mkvec2s(sp, sn);
    1746             : }
    1747             : 
    1748             : GEN
    1749          21 : qfgaussred(GEN a) { return gaussred(a,0); }
    1750             : 
    1751             : GEN
    1752           7 : qfgaussred2(GEN a)
    1753             : {
    1754           7 :   pari_sp av = avma;
    1755           7 :   GEN M = qfgaussred(a);
    1756           7 :   long i, l = lg(M);
    1757           7 :   GEN D = cgetg(l, t_VEC);
    1758          35 :   for (i = 1; i < l; i++)
    1759             :   {
    1760          28 :     gel(D,i) = gcoeff(M,i,i);
    1761          28 :     gcoeff(M,i,i) = gen_1;
    1762             :   }
    1763           7 :   return gc_GEN(av, mkvec2(M,D));
    1764             : }
    1765             : 
    1766             : GEN
    1767          21 : qfgaussred0(GEN a, long flag)
    1768             : {
    1769          21 :   switch (flag)
    1770             :   {
    1771          14 :     case 0: return qfgaussred(a);
    1772           7 :     case 1: return qfgaussred2(a);
    1773           0 :     default: pari_err_FLAG("qfgaussred");
    1774             :   }
    1775             :   return NULL; /* LCOV_EXCL_LINE */
    1776             : }
    1777             : 
    1778             : GEN
    1779          14 : qfcholesky(GEN a, long prec)
    1780             : {
    1781             :   GEN M;
    1782          14 :   long n = lg(a);
    1783          14 :   if (typ(a) != t_MAT) pari_err_TYPE("qfcholesky",a);
    1784          14 :   if (n == 1) return cgetg(1, t_MAT);
    1785          14 :   if (lgcols(a) != lg(a)) pari_err_DIM("qfcholesky");
    1786          14 :   M =  RgM_Cholesky(a, prec);
    1787          14 :   if (!M) return cgetg(1, t_VEC);
    1788           7 :   return M;
    1789             : }
    1790             : 
    1791             : GEN
    1792       12179 : qfsign(GEN a) { return gaussred(a,1); }
    1793             : 
    1794             : /* x -= s(y+u*x) */
    1795             : /* y += s(x-u*y), simultaneously */
    1796             : static void
    1797       19180 : rot(GEN x, GEN y, GEN s, GEN u) {
    1798       19180 :   GEN x1 = subrr(x, mulrr(s,addrr(y,mulrr(u,x))));
    1799       19180 :   GEN y1 = addrr(y, mulrr(s,subrr(x,mulrr(u,y))));
    1800       19180 :   affrr(x1,x);
    1801       19180 :   affrr(y1,y);
    1802       19180 : }
    1803             : 
    1804             : /* Diagonalization of a REAL symmetric matrix. Return a vector [L, r]:
    1805             :  * L = vector of eigenvalues
    1806             :  * r = matrix of eigenvectors */
    1807             : GEN
    1808          28 : jacobi(GEN a, long prec)
    1809             : {
    1810             :   pari_sp av;
    1811          28 :   long de, e, e1, e2, i, j, p, q, l = lg(a);
    1812             :   GEN c, ja, L, r, L2, r2, unr, sqrt2;
    1813             : 
    1814          28 :   if (typ(a) != t_MAT) pari_err_TYPE("jacobi",a);
    1815          28 :   ja = cgetg(3,t_VEC);
    1816          28 :   L = cgetg(l,t_COL); gel(ja,1) = L;
    1817          28 :   r = cgetg(l,t_MAT); gel(ja,2) = r;
    1818          28 :   if (l == 1) return ja;
    1819          28 :   if (lgcols(a) != l) pari_err_DIM("jacobi");
    1820             : 
    1821          28 :   e1 = HIGHEXPOBIT-1;
    1822         224 :   for (j=1; j<l; j++)
    1823             :   {
    1824         196 :     GEN z = gtofp(gcoeff(a,j,j), prec);
    1825         196 :     gel(L,j) = z;
    1826         196 :     e = expo(z); if (e < e1) e1 = e;
    1827             :   }
    1828         224 :   for (j=1; j<l; j++)
    1829             :   {
    1830         196 :     gel(r,j) = cgetg(l,t_COL);
    1831        1582 :     for (i=1; i<l; i++) gcoeff(r,i,j) = utor(i==j? 1: 0, prec);
    1832             :   }
    1833          28 :   av = avma;
    1834             : 
    1835          28 :   e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1836          28 :   c = cgetg(l,t_MAT);
    1837         224 :   for (j=1; j<l; j++)
    1838             :   {
    1839         196 :     gel(c,j) = cgetg(j,t_COL);
    1840         791 :     for (i=1; i<j; i++)
    1841             :     {
    1842         595 :       GEN z = gtofp(gcoeff(a,i,j), prec);
    1843         595 :       gcoeff(c,i,j) = z;
    1844         595 :       if (!signe(z)) continue;
    1845         308 :       e = expo(z); if (e > e2) { e2 = e; p = i; q = j; }
    1846             :     }
    1847             :   }
    1848          28 :   a = c; unr = real_1(prec);
    1849          28 :   sqrt2 = sqrtr_abs(shiftr(unr, 1));
    1850          28 :   de = prec2nbits(prec);
    1851             : 
    1852             :  /* e1 = min expo(a[i,i])
    1853             :   * e2 = max expo(a[i,j]), i < j, occurs at a[p,q] (p < q)*/
    1854        1568 :   while (e1-e2 < de)
    1855             :   {
    1856        1540 :     pari_sp av2 = avma;
    1857             :     GEN x, y, t, c, s, u;
    1858             :     /* compute attached rotation in the plane formed by basis vectors p and q */
    1859        1540 :     x = subrr(gel(L,q),gel(L,p));
    1860        1540 :     if (signe(x))
    1861             :     {
    1862        1512 :       x = divrr(x, shiftr(gcoeff(a,p,q),1));
    1863        1512 :       y = sqrtr(addrr(unr, sqrr(x)));
    1864        1512 :       t = invr((signe(x)>0)? addrr(x,y): subrr(x,y));
    1865        1512 :       c = sqrtr(addrr(unr,sqrr(t)));
    1866        1512 :       s = divrr(t,c);
    1867        1512 :       u = divrr(t,addrr(unr,c));
    1868             :     }
    1869             :     else /* same formulas for t = 1.0 */
    1870             :     {
    1871          28 :       t = NULL; /* 1.0 */
    1872          28 :       c = sqrt2;
    1873          28 :       s = shiftr(c, -1);
    1874          28 :       u = subrr(c, unr);
    1875             :     }
    1876             : 
    1877             :     /* compute successive transforms of a and the matrix of accumulated
    1878             :      * rotations (r) */
    1879        4158 :     for (i=1;   i<p; i++) rot(gcoeff(a,i,p), gcoeff(a,i,q), s,u);
    1880        4025 :     for (i=p+1; i<q; i++) rot(gcoeff(a,p,i), gcoeff(a,i,q), s,u);
    1881        4487 :     for (i=q+1; i<l; i++) rot(gcoeff(a,p,i), gcoeff(a,q,i), s,u);
    1882        1540 :     y = gcoeff(a,p,q); t = t? mulrr(t, y): rcopy(y);
    1883        1540 :     shiftr_inplace(y, -de - 1);
    1884        1540 :     affrr(subrr(gel(L,p),t), gel(L,p));
    1885        1540 :     affrr(addrr(gel(L,q),t), gel(L,q));
    1886       12670 :     for (i=1; i<l; i++) rot(gcoeff(r,i,p), gcoeff(r,i,q), s,u);
    1887             : 
    1888        1540 :     e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1889       12670 :     for (j=1; j<l; j++)
    1890       46396 :       for (i=1; i<j; i++)
    1891             :       {
    1892       35266 :         GEN z = gcoeff(a,i,j);
    1893       35266 :         if (!signe(z)) continue;
    1894       31066 :         e = expo(z); if (e > e2) { e2 = e; p = i; q = j; }
    1895             :       }
    1896        1540 :     set_avma(av2);
    1897             :   }
    1898             :   /* sort eigenvalues from smallest to largest */
    1899          28 :   c = indexsort(L);
    1900         224 :   r2 = vecpermute(r, c); for (i=1; i<l; i++) gel(r,i) = gel(r2,i);
    1901         224 :   L2 = vecpermute(L, c); for (i=1; i<l; i++) gel(L,i) = gel(L2,i);
    1902          28 :   set_avma(av); return ja;
    1903             : }
    1904             : 
    1905             : /*************************************************************************/
    1906             : /**                                                                     **/
    1907             : /**                   Q-vector space -> Z-modules                       **/
    1908             : /**                                                                     **/
    1909             : /*************************************************************************/
    1910             : 
    1911             : GEN
    1912         133 : matrixqz0(GEN x,GEN p)
    1913             : {
    1914         133 :   if (typ(x) != t_MAT) pari_err_TYPE("matrixqz",x);
    1915         133 :   if (!p) return QM_minors_coprime(x,NULL);
    1916          98 :   if (typ(p) != t_INT) pari_err_TYPE("matrixqz",p);
    1917          98 :   if (signe(p)>=0) return QM_minors_coprime(x,p);
    1918          91 :   if (!RgM_is_QM(x)) pari_err_TYPE("matrixqz", x);
    1919          91 :   if (absequaliu(p,1)) return QM_ImZ(x); /* p = -1 */
    1920          63 :   if (absequaliu(p,2)) return QM_ImQ(x); /* p = -2 */
    1921           7 :   pari_err_FLAG("QM_minors_coprime");
    1922             :   return NULL; /* LCOV_EXCL_LINE */
    1923             : }
    1924             : 
    1925             : GEN
    1926          42 : QM_minors_coprime(GEN x, GEN D)
    1927             : {
    1928          42 :   pari_sp av = avma, av1;
    1929             :   long i, j, m, n, lP;
    1930             :   GEN P, y;
    1931             : 
    1932          42 :   n = lg(x)-1; if (!n) return gcopy(x);
    1933          42 :   m = nbrows(x);
    1934          42 :   if (n > m) pari_err_DOMAIN("QM_minors_coprime","n",">",strtoGENstr("m"),x);
    1935          35 :   y = x; x = cgetg(n+1,t_MAT);
    1936         112 :   for (j=1; j<=n; j++)
    1937             :   {
    1938          77 :     gel(x,j) = Q_primpart(gel(y,j));
    1939          77 :     RgV_check_ZV(gel(x,j), "QM_minors_coprime");
    1940             :   }
    1941             :   /* x now a ZM */
    1942          35 :   if (n==m)
    1943             :   {
    1944          21 :     if (gequal0(ZM_det(x)))
    1945          14 :       pari_err_DOMAIN("QM_minors_coprime", "rank(A)", "<",stoi(n),x);
    1946           7 :     set_avma(av); return matid(n);
    1947             :   }
    1948             :   /* m > n */
    1949          14 :   if (!D || gequal0(D))
    1950             :   {
    1951          14 :     pari_sp av2 = avma;
    1952          14 :     D = ZM_detmult(shallowtrans(x));
    1953          14 :     if (is_pm1(D)) { set_avma(av2); return ZM_copy(x); }
    1954             :   }
    1955          14 :   P = gel(Z_factor(D), 1); lP = lg(P);
    1956          14 :   av1 = avma;
    1957          56 :   for (i=1; i < lP; i++)
    1958             :   {
    1959          42 :     GEN p = gel(P,i), pov2 = shifti(p, -1);
    1960             :     for(;;)
    1961          42 :     {
    1962          84 :       GEN N, M = FpM_ker(x, p);
    1963          84 :       long lM = lg(M);
    1964          84 :       if (lM==1) break;
    1965             : 
    1966          42 :       FpM_center_inplace(M, p, pov2);
    1967          42 :       N = ZM_Z_divexact(ZM_mul(x,M), p);
    1968         126 :       for (j=1; j<lM; j++)
    1969             :       {
    1970         147 :         long k = n; while (!signe(gcoeff(M,k,j))) k--;
    1971          84 :         gel(x,k) = gel(N,j);
    1972             :       }
    1973          42 :       if (gc_needed(av1,1))
    1974             :       {
    1975           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"QM_minors_coprime, p = %Ps", p);
    1976           0 :         x = gc_GEN(av1, x); pov2 = shifti(p, -1);
    1977             :       }
    1978             :     }
    1979             :   }
    1980          14 :   return gc_GEN(av, x);
    1981             : }
    1982             : 
    1983             : static GEN
    1984        3479 : QM_ImZ_all_i(GEN A, GEN *U, long remove, long hnf, long linindep)
    1985             : {
    1986        3479 :   GEN V = NULL, D;
    1987        3479 :   A = Q_remove_denom(A,&D);
    1988        3479 :   if (D)
    1989             :   {
    1990             :     long l, lA;
    1991        1190 :     V = matkermod(A,D,NULL);
    1992        1190 :     l = lg(V); lA = lg(A);
    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        3479 :   if (!linindep && ZM_rank(A)==lg(A)-1) linindep = 1;
    2001        3479 :   if (hnf || !linindep) A = ZM_hnflll(A, U, remove);
    2002        3479 :   if (U && V)
    2003             :   {
    2004        1057 :     if (hnf)    *U = ZM_mul(V,*U);
    2005           0 :     else        *U = V;
    2006             :   }
    2007        3479 :   return A;
    2008             : }
    2009             : GEN
    2010          28 : QM_ImZ_all(GEN x, GEN *U, long remove, long hnf)
    2011             : {
    2012          28 :   pari_sp av = avma;
    2013          28 :   x = QM_ImZ_all_i(x, U, remove, hnf, 0);
    2014          28 :   return gc_all(av, U?2:1, &x, &U);
    2015             : }
    2016             : GEN
    2017           0 : QM_ImZ_hnfall(GEN x, GEN *U, long remove) { return QM_ImZ_all(x, U, remove, 1); }
    2018             : GEN
    2019           0 : QM_ImZ_hnf(GEN x) { return QM_ImZ_hnfall(x, NULL, 1); }
    2020             : GEN
    2021          28 : QM_ImZ(GEN x) { return QM_ImZ_all(x, NULL, 1, 0); }
    2022             : 
    2023             : GEN
    2024        3458 : QM_ImQ_all(GEN x, GEN *U, long remove, long hnf)
    2025             : {
    2026        3458 :   pari_sp av = avma;
    2027        3458 :   long i, n = lg(x), m;
    2028        3458 :   GEN ir, V, D, c, K = NULL;
    2029             : 
    2030        3458 :   if (U) *U = matid(n-1);
    2031        3458 :   if (n==1) return gcopy(x);
    2032        3451 :   m = lg(gel(x,1));
    2033             : 
    2034        3451 :   x = RgM_shallowcopy(x);
    2035       15029 :   for (i=1; i<n; i++)
    2036             :   {
    2037       11578 :     gel(x,i) = Q_primitive_part(gel(x,i), &c);
    2038       11578 :     if (U && c && signe(c)) gcoeff(*U,i,i) = ginv(c);
    2039             :   }
    2040             : 
    2041        3451 :   ir = ZM_indexrank(x);
    2042        3451 :   if (U)
    2043             :   {
    2044        2219 :     *U = vecpermute(*U, gel(ir,2));
    2045        2219 :     if (remove < 2) K = ZM_ker(x);
    2046             :   }
    2047        3451 :   x = vecpermute(x, gel(ir,2));
    2048             : 
    2049        3451 :   D = absi(ZM_det(rowpermute(x,gel(ir,1))));
    2050        3451 :   x = RgM_Rg_div(x, D);
    2051        3451 :   x = QM_ImZ_all_i(x, U? &V: NULL, remove, hnf, 1);
    2052             : 
    2053        3451 :   if (U)
    2054             :   {
    2055        2219 :     *U = RgM_Rg_div(RgM_mul(*U,V),D);
    2056        2219 :     if (remove < 2) *U = shallowconcat(K,*U);
    2057        2219 :     if (!remove) x = shallowconcat(zeromatcopy(m-1,lg(K)-1), x);
    2058        2219 :     (void)gc_all(av, 2, &x, U);
    2059             :   }
    2060        1232 :   else x = gc_GEN(av,x);
    2061        3451 :   return x;
    2062             : }
    2063             : GEN
    2064        3402 : QM_ImQ_hnfall(GEN x, GEN *U, long remove) { return QM_ImQ_all(x, U, remove, 1); }
    2065             : GEN
    2066        1183 : QM_ImQ_hnf(GEN x) { return QM_ImQ_hnfall(x, NULL, 1); }
    2067             : GEN
    2068          56 : QM_ImQ(GEN x) { return QM_ImQ_all(x, NULL, 1, 0); }
    2069             : 
    2070             : GEN
    2071        5691 : intersect(GEN x, GEN y)
    2072             : {
    2073        5691 :   long j, lx = lg(x);
    2074             :   pari_sp av;
    2075             :   GEN z;
    2076             : 
    2077        5691 :   if (typ(x)!=t_MAT) pari_err_TYPE("intersect",x);
    2078        5691 :   if (typ(y)!=t_MAT) pari_err_TYPE("intersect",y);
    2079        5691 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
    2080             : 
    2081        4207 :   av = avma; z = ker(shallowconcat(x,y));
    2082       18676 :   for (j=lg(z)-1; j; j--) setlg(z[j], lx);
    2083        4207 :   return gc_upto(av, image(RgM_mul(x,z)));
    2084             : }

Generated by: LCOV version 1.16