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

Generated by: LCOV version 1.13