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 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.10.0 lcov report (development 20459-9710128) Lines: 907 1018 89.1 %
Date: 2017-04-28 05:33:48 Functions: 67 72 93.1 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                         LINEAR ALGEBRA                         **/
      17             : /**                         (second part)                          **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : /*******************************************************************/
      23             : /*                                                                 */
      24             : /*                   CHARACTERISTIC POLYNOMIAL                     */
      25             : /*                                                                 */
      26             : /*******************************************************************/
      27             : 
      28             : static GEN
      29        8898 : Flm_charpoly_i(GEN x, ulong p)
      30             : {
      31        8898 :   long lx = lg(x), r, i;
      32        8898 :   GEN H, y = cgetg(lx+1, t_VEC);
      33        8898 :   gel(y,1) = pol1_Flx(0); H = Flm_hess(x, p);
      34       49572 :   for (r = 1; r < lx; r++)
      35             :   {
      36       40674 :     pari_sp av2 = avma;
      37       40674 :     ulong a = 1;
      38       40674 :     GEN z, b = zero_Flx(0);
      39      111450 :     for (i = r-1; i; i--)
      40             :     {
      41       84120 :       a = Fl_mul(a, ucoeff(H,i+1,i), p);
      42       84120 :       if (!a) break;
      43       70776 :       b = Flx_add(b, Flx_Fl_mul(gel(y,i), Fl_mul(a,ucoeff(H,i,r),p), p), p);
      44             :     }
      45       81348 :     z = Flx_sub(Flx_shift(gel(y,r), 1),
      46       81348 :                 Flx_Fl_mul(gel(y,r), ucoeff(H,r,r), p), p);
      47             :     /* (X - H[r,r])y[r] - b */
      48       40674 :     gel(y,r+1) = gerepileuptoleaf(av2, Flx_sub(z, b, p));
      49             :   }
      50        8898 :   return gel(y,lx);
      51             : }
      52             : 
      53             : GEN
      54           0 : Flm_charpoly(GEN x, ulong p)
      55             : {
      56           0 :   pari_sp av = avma;
      57           0 :   return gerepileuptoleaf(av, Flm_charpoly_i(x,p));
      58             : }
      59             : 
      60             : GEN
      61        6342 : FpM_charpoly(GEN x, GEN p)
      62             : {
      63        6342 :   pari_sp av = avma;
      64             :   long lx, r, i;
      65             :   GEN y, H;
      66             : 
      67        6342 :   if (lgefint(p) == 3)
      68             :   {
      69        6330 :     ulong pp = p[2];
      70        6330 :     y = Flx_to_ZX(Flm_charpoly_i(ZM_to_Flm(x,pp), pp));
      71        6330 :     return gerepileupto(av, y);
      72             :   }
      73          12 :   lx = lg(x); y = cgetg(lx+1, t_VEC);
      74          12 :   gel(y,1) = pol_1(0); H = FpM_hess(x, p);
      75          36 :   for (r = 1; r < lx; r++)
      76             :   {
      77          36 :     pari_sp av2 = avma;
      78          36 :     GEN z, a = gen_1, b = pol_0(0);
      79          72 :     for (i = r-1; i; i--)
      80             :     {
      81          42 :       a = Fp_mul(a, gcoeff(H,i+1,i), p);
      82          42 :       if (!signe(a)) break;
      83          36 :       b = ZX_add(b, ZX_Z_mul(gel(y,i), Fp_mul(a,gcoeff(H,i,r),p)));
      84             :     }
      85          36 :     b = FpX_red(b, p);
      86          72 :     z = FpX_sub(RgX_shift_shallow(gel(y,r), 1),
      87          72 :                 FpX_Fp_mul(gel(y,r), gcoeff(H,r,r), p), p);
      88          36 :     z = FpX_sub(z,b,p);
      89          36 :     if (r+1 == lx) { gel(y,lx) = z; break; }
      90          24 :     gel(y,r+1) = gerepileupto(av2, z); /* (X - H[r,r])y[r] - b */
      91             :   }
      92          12 :   return gerepileupto(av, gel(y,lx));
      93             : }
      94             : 
      95             : GEN
      96         246 : charpoly0(GEN x, long v, long flag)
      97             : {
      98         246 :   if (v<0) v = 0;
      99         246 :   switch(flag)
     100             :   {
     101          12 :     case 0: return caradj(x,v,NULL);
     102          12 :     case 1: return caract(x,v);
     103          12 :     case 2: return carhess(x,v);
     104          12 :     case 3: return carberkowitz(x,v);
     105             :     case 4:
     106           6 :       if (typ(x) != t_MAT) pari_err_TYPE("charpoly",x);
     107           6 :       RgM_check_ZM(x, "charpoly");
     108           6 :       x = ZM_charpoly(x); setvarn(x, v); return x;
     109             :     case 5:
     110         192 :       return charpoly(x, v);
     111             :   }
     112           0 :   pari_err_FLAG("charpoly");
     113             :   return NULL; /* LCOV_EXCL_LINE */
     114             : }
     115             : 
     116             : /* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
     117             : static GEN
     118        5382 : easychar(GEN x, long v)
     119             : {
     120             :   pari_sp av;
     121             :   long lx;
     122             :   GEN p1;
     123             : 
     124        5382 :   switch(typ(x))
     125             :   {
     126             :     case t_INT: case t_REAL: case t_INTMOD:
     127             :     case t_FRAC: case t_PADIC:
     128          30 :       p1=cgetg(4,t_POL);
     129          30 :       p1[1]=evalsigne(1) | evalvarn(v);
     130          30 :       gel(p1,2) = gneg(x); gel(p1,3) = gen_1;
     131          30 :       return p1;
     132             : 
     133             :     case t_COMPLEX: case t_QUAD:
     134          12 :       p1 = cgetg(5,t_POL);
     135          12 :       p1[1] = evalsigne(1) | evalvarn(v);
     136          12 :       gel(p1,2) = gnorm(x); av = avma;
     137          12 :       gel(p1,3) = gerepileupto(av, gneg(gtrace(x)));
     138          12 :       gel(p1,4) = gen_1; return p1;
     139             : 
     140             :     case t_FFELT: {
     141          24 :       pari_sp ltop=avma;
     142          24 :       p1 = FpX_to_mod(FF_charpoly(x), FF_p_i(x));
     143          24 :       setvarn(p1,v); return gerepileupto(ltop,p1);
     144             :     }
     145             : 
     146             :     case t_POLMOD:
     147             :     {
     148         714 :       GEN A = gel(x,2), T = gel(x,1);
     149         714 :       if (typ(A)==t_POL && RgX_is_QX(A) && RgX_is_ZX(T))
     150         702 :         return QXQ_charpoly(A, T, v);
     151             :       else
     152          12 :         return RgXQ_charpoly(A, T, v);
     153             :     }
     154             :     case t_MAT:
     155        4602 :       lx=lg(x);
     156        4602 :       if (lx==1) return pol_1(v);
     157        4554 :       if (lgcols(x) != lx) break;
     158        4548 :       return NULL;
     159             :   }
     160           6 :   pari_err_TYPE("easychar",x);
     161             :   return NULL; /* LCOV_EXCL_LINE */
     162             : }
     163             : /* compute charpoly by mapping to Fp first, return lift to Z */
     164             : static GEN
     165          30 : RgM_Fp_charpoly(GEN x, GEN p, long v)
     166             : {
     167             :   GEN T;
     168          30 :   if (lgefint(p) == 3)
     169             :   {
     170          18 :     ulong pp = itou(p);
     171          18 :     T = Flm_charpoly_i(RgM_to_Flm(x, pp), pp);
     172          18 :     T = Flx_to_ZX(T);
     173             :   }
     174             :   else
     175          12 :     T = FpM_charpoly(RgM_to_FpM(x, p), p);
     176          30 :   setvarn(T, v); return T;
     177             : }
     178             : GEN
     179        1752 : charpoly(GEN x, long v)
     180             : {
     181        1752 :   GEN T, p = NULL;
     182        1752 :   if ((T = easychar(x,v))) return T;
     183        1656 :   if (RgM_is_ZM(x))
     184             :   {
     185        1386 :     T = ZM_charpoly(x);
     186        1386 :     setvarn(T, v);
     187             :   }
     188         270 :   else if (RgM_is_FpM(x, &p) && BPSW_psp(p))
     189          30 :   {
     190          30 :     pari_sp av = avma;
     191          30 :     T = RgM_Fp_charpoly(x,p,v);
     192          30 :     T = gerepileupto(av, FpX_to_mod(T,p));
     193             :   }
     194         240 :   else if (isinexact(x))
     195           6 :     T = carhess(x, v);
     196             :   else
     197         234 :     T = carberkowitz(x, v);
     198        1650 :   return T;
     199             : }
     200             : 
     201             : /* We possibly worked with an "invalid" polynomial p, satisfying
     202             :  * varn(p) > gvar2(p). Fix this. */
     203             : static GEN
     204        2838 : fix_pol(pari_sp av, GEN p)
     205             : {
     206        2838 :   long w = gvar2(p), v = varn(p);
     207        2838 :   if (w == v) pari_err_PRIORITY("charpoly", p, "=", w);
     208        2832 :   if (varncmp(w,v) < 0) p = gerepileupto(av, poleval(p, pol_x(v)));
     209        2832 :   return p;
     210             : }
     211             : GEN
     212          12 : caract(GEN x, long v)
     213             : {
     214          12 :   pari_sp av = avma;
     215             :   GEN  T, C, x_k, Q;
     216             :   long k, n;
     217             : 
     218          12 :   if ((T = easychar(x,v))) return T;
     219             : 
     220          12 :   n = lg(x)-1;
     221          12 :   if (n == 1) return fix_pol(av, deg1pol(gen_1, gneg(gcoeff(x,1,1)), v));
     222             : 
     223          12 :   x_k = pol_x(v); /* to be modified in place */
     224          12 :   T = scalarpol(det(x), v); C = utoineg(n); Q = pol_x(v);
     225          24 :   for (k=1; k<=n; k++)
     226             :   {
     227          24 :     GEN mk = utoineg(k), d;
     228          24 :     gel(x_k,2) = mk;
     229          24 :     d = det(RgM_Rg_add_shallow(x, mk));
     230          24 :     T = RgX_add(RgX_mul(T, x_k), RgX_Rg_mul(Q, gmul(C, d)));
     231          24 :     if (k == n) break;
     232             : 
     233          12 :     Q = RgX_mul(Q, x_k);
     234          12 :     C = diviuexact(mulsi(k-n,C), k+1); /* (-1)^k binomial(n,k) */
     235             :   }
     236          12 :   return fix_pol(av, RgX_Rg_div(T, mpfact(n)));
     237             : }
     238             : 
     239             : /* C = charpoly(x, v) */
     240             : static GEN
     241          18 : RgM_adj_from_char(GEN x, long v, GEN C)
     242             : {
     243          18 :   if (varn(C) != v) /* problem with variable priorities */
     244             :   {
     245           6 :     C = gdiv(gsub(C, gsubst(C, v, gen_0)), pol_x(v));
     246           6 :     if (odd(lg(x))) C = RgX_neg(C); /* even dimension */
     247           6 :     return gsubst(C, v, x);
     248             :   }
     249             :   else
     250             :   {
     251          12 :     C = RgX_shift_shallow(C, -1);
     252          12 :     if (odd(lg(x))) C = RgX_neg(C); /* even dimension */
     253          12 :     return RgX_RgM_eval(C, x);
     254             :   }
     255             : }
     256             : /* assume x square matrice */
     257             : static GEN
     258        4194 : mattrace(GEN x)
     259             : {
     260        4194 :   long i, lx = lg(x);
     261             :   GEN t;
     262        4194 :   if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
     263        4164 :   t = gcoeff(x,1,1);
     264        4164 :   for (i = 2; i < lx; i++) t = gadd(t, gcoeff(x,i,i));
     265        4164 :   return t;
     266             : }
     267             : static int
     268        1992 : bad_char(GEN q, long n)
     269             : {
     270             :   forprime_t S;
     271             :   ulong p;
     272        1992 :   if (!signe(q)) return 0;
     273          36 :   (void)u_forprime_init(&S, 2, n);
     274          36 :   while ((p = u_forprime_next(&S)))
     275          60 :     if (!umodiu(q, p)) return 1;
     276          24 :   return 0;
     277             : }
     278             : /* Using traces: return the characteristic polynomial of x (in variable v).
     279             :  * If py != NULL, the adjoint matrix is put there. */
     280             : GEN
     281        2046 : caradj(GEN x, long v, GEN *py)
     282             : {
     283             :   pari_sp av, av0;
     284             :   long i, k, n;
     285             :   GEN T, y, t;
     286             : 
     287        2046 :   if ((T = easychar(x, v)))
     288             :   {
     289          36 :     if (py)
     290             :     {
     291          36 :       if (typ(x) != t_MAT) pari_err_TYPE("matadjoint",x);
     292          36 :       *py = cgetg(1,t_MAT);
     293             :     }
     294          36 :     return T;
     295             :   }
     296             : 
     297        2010 :   n = lg(x)-1; av0 = avma;
     298        2010 :   T = cgetg(n+3,t_POL); T[1] = evalsigne(1) | evalvarn(v);
     299        2010 :   gel(T,n+2) = gen_1;
     300        2010 :   if (!n) { if (py) *py = cgetg(1,t_MAT); return T; }
     301        2010 :   av = avma; t = gerepileupto(av, gneg(mattrace(x)));
     302        2010 :   gel(T,n+1) = t;
     303        2010 :   if (n == 1) {
     304           6 :     T = fix_pol(av0, T);
     305           6 :     if (py) *py = matid(1); return T;
     306             :   }
     307        2004 :   if (n == 2) {
     308          12 :     GEN a = gcoeff(x,1,1), b = gcoeff(x,1,2);
     309          12 :     GEN c = gcoeff(x,2,1), d = gcoeff(x,2,2);
     310          12 :     av = avma;
     311          12 :     gel(T,2) = gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
     312          12 :     T = fix_pol(av0, T);
     313          12 :     if (py) {
     314           6 :       y = cgetg(3, t_MAT);
     315           6 :       gel(y,1) = mkcol2(gcopy(d), gneg(c));
     316           6 :       gel(y,2) = mkcol2(gneg(b), gcopy(a));
     317           6 :       *py = y;
     318             :     }
     319          12 :     return T;
     320             :   }
     321             :   /* l > 3 */
     322        1992 :   if (bad_char(residual_characteristic(x), n))
     323             :   { /* n! not invertible in base ring */
     324          12 :     T = charpoly(x, v);
     325          12 :     if (!py) return gerepileupto(av, T);
     326          12 :     *py = RgM_adj_from_char(x, v, T);
     327          12 :     gerepileall(av, 2, &T,py);
     328          12 :     return T;
     329             :   }
     330        1980 :   av = avma; y = RgM_shallowcopy(x);
     331        1980 :   for (i = 1; i <= n; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
     332        3966 :   for (k = 2; k < n; k++)
     333             :   {
     334        1986 :     GEN y0 = y;
     335        1986 :     y = RgM_mul(y, x);
     336        1986 :     t = gdivgs(mattrace(y), -k);
     337        1986 :     for (i = 1; i <= n; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
     338        1986 :     y = gclone(y);
     339        1986 :     gel(T,n-k+2) = gerepilecopy(av, t); av = avma;
     340        1986 :     if (k > 2) gunclone(y0);
     341             :   }
     342        1980 :   t = gmul(gcoeff(x,1,1),gcoeff(y,1,1));
     343        1980 :   for (i=2; i<=n; i++) t = gadd(t, gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
     344        1980 :   gel(T,2) = gerepileupto(av, gneg(t));
     345        1980 :   T = fix_pol(av0, T);
     346        1980 :   if (py) *py = odd(n)? gcopy(y): RgM_neg(y);
     347        1980 :   gunclone(y); return T;
     348             : }
     349             : 
     350             : GEN
     351        2034 : adj(GEN x)
     352             : {
     353             :   GEN y;
     354        2034 :   (void)caradj(x, fetch_var(), &y);
     355        2034 :   (void)delete_var(); return y;
     356             : }
     357             : 
     358             : GEN
     359           6 : adjsafe(GEN x)
     360             : {
     361           6 :   const long v = fetch_var();
     362           6 :   pari_sp av = avma;
     363             :   GEN C, A;
     364           6 :   if (typ(x) != t_MAT) pari_err_TYPE("matadjoint",x);
     365           6 :   if (lg(x) < 3) return gcopy(x);
     366           6 :   C = charpoly(x,v);
     367           6 :   A = RgM_adj_from_char(x, v, C);
     368           6 :   (void)delete_var(); return gerepileupto(av, A);
     369             : }
     370             : 
     371             : GEN
     372          96 : matadjoint0(GEN x, long flag)
     373             : {
     374          96 :   switch(flag)
     375             :   {
     376          90 :     case 0: return adj(x);
     377           6 :     case 1: return adjsafe(x);
     378             :   }
     379           0 :   pari_err_FLAG("matadjoint");
     380             :   return NULL; /* LCOV_EXCL_LINE */
     381             : }
     382             : 
     383             : /*******************************************************************/
     384             : /*                                                                 */
     385             : /*                       Frobenius form                            */
     386             : /*                                                                 */
     387             : /*******************************************************************/
     388             : 
     389             : /* The following section implement a mix of Ozello and Storjohann algorithms
     390             : 
     391             : P. Ozello, doctoral thesis (in French):
     392             : Calcul exact des formes de Jordan et de Frobenius d'une matrice, Chapitre 2
     393             : http://tel.archives-ouvertes.fr/tel-00323705
     394             : 
     395             : A. Storjohann,  Diss. ETH No. 13922
     396             : Algorithms for Matrix Canonical Forms, Chapter 9
     397             : https://cs.uwaterloo.ca/~astorjoh/diss2up.pdf
     398             : 
     399             : We use Storjohann Lemma 9.14 (step1, step2, step3) Ozello theorem 4,
     400             : and Storjohann Lemma 9.18
     401             : */
     402             : 
     403             : /* Elementary transforms */
     404             : 
     405             : /* M <- U^(-1) M U, U = E_{i,j}(k) => U^(-1) = E{i,j}(-k)
     406             :  * P = U * P */
     407             : static void
     408       12522 : transL(GEN M, GEN P, GEN k, long i, long j)
     409             : {
     410       12522 :   long l, n = lg(M)-1;
     411      169206 :   for(l=1; l<=n; l++) /* M[,j]-=k*M[,i] */
     412      156684 :     gcoeff(M,l,j) = gsub(gcoeff(M,l,j), gmul(gcoeff(M,l,i), k));
     413      169206 :   for(l=1; l<=n; l++) /* M[i,]+=k*M[j,] */
     414      156684 :     gcoeff(M,i,l) = gadd(gcoeff(M,i,l), gmul(gcoeff(M,j,l), k));
     415       12522 :   if (P)
     416      160656 :     for(l=1; l<=n; l++)
     417      149076 :       gcoeff(P,i,l) = gadd(gcoeff(P,i,l), gmul(gcoeff(P,j,l), k));
     418       12522 : }
     419             : 
     420             : /* j = a or b */
     421             : static void
     422        1746 : transD(GEN M, GEN P, long a, long b, long j)
     423             : {
     424             :   long l, n;
     425        1746 :   GEN k = gcoeff(M,a,b), ki;
     426             : 
     427        3492 :   if (gequal1(k)) return;
     428         948 :   ki = ginv(k); n = lg(M)-1;
     429       11298 :   for(l=1; l<=n; l++)
     430       10350 :     if (l!=j)
     431             :     {
     432        9402 :       gcoeff(M,l,j) = gmul(gcoeff(M,l,j), k);
     433        9402 :       gcoeff(M,j,l) = (j==a && l==b)? gen_1: gmul(gcoeff(M,j,l), ki);
     434             :     }
     435         948 :   if (P)
     436        9288 :     for(l=1; l<=n; l++)
     437        8574 :       gcoeff(P,j,l) = gmul(gcoeff(P,j,l), ki);
     438             : }
     439             : 
     440             : static void
     441         174 : transS(GEN M, GEN P, long i, long j)
     442             : {
     443         174 :   long l, n = lg(M)-1;
     444         174 :   swap(gel(M,i), gel(M,j));
     445        1698 :   for (l=1; l<=n; l++)
     446        1524 :     swap(gcoeff(M,i,l), gcoeff(M,j,l));
     447         174 :   if (P)
     448         468 :     for (l=1; l<=n; l++)
     449         432 :       swap(gcoeff(P,i,l), gcoeff(P,j,l));
     450         174 : }
     451             : 
     452             : /* Convert companion matrix to polynomial*/
     453             : static GEN
     454         168 : minpoly_polslice(GEN M, long i, long j, long v)
     455             : {
     456         168 :   long k, d = j+1-i;
     457         168 :   GEN P = cgetg(d+3,t_POL);
     458         168 :   P[1] = evalsigne(1)|evalvarn(v);
     459         846 :   for (k=0; k<d; k++)
     460         678 :     gel(P,k+2) = gneg(gcoeff(M,i+k, j));
     461         168 :   gel(P,d+2) = gen_1;
     462         168 :   return P;
     463             : }
     464             : 
     465             : static GEN
     466          42 : minpoly_listpolslice(GEN M, GEN V, long v)
     467             : {
     468          42 :   long i, n = lg(M)-1, nb = lg(V)-1;
     469          42 :   GEN W = cgetg(nb+1, t_VEC);
     470         126 :   for (i=1; i<=nb; i++)
     471          84 :     gel(W,i) = minpoly_polslice(M, V[i], i < nb? V[i+1]-1: n, v);
     472          42 :   return W;
     473             : }
     474             : 
     475             : static int
     476          42 : minpoly_dvdslice(GEN M, long i, long j, long k)
     477             : {
     478          42 :   pari_sp av = avma;
     479          42 :   long r = signe(RgX_rem(minpoly_polslice(M, i, j-1, 0),
     480             :                         minpoly_polslice(M, j, k, 0)));
     481          42 :   avma = av; return r==0;
     482             : }
     483             : 
     484             : static void
     485           0 : RgM_replace(GEN M, GEN M2)
     486             : {
     487           0 :   long n = lg(M)-1, m = nbrows(M), i, j;
     488           0 :   for(i=1; i<=n; i++)
     489           0 :     for(j=1; j<=m; j++)
     490           0 :       gcoeff(M, i, j) = gcoeff(M2, i, j);
     491           0 : }
     492             : 
     493             : static void
     494           0 : gerepilemat2_inplace(pari_sp av, GEN M, GEN P)
     495             : {
     496           0 :   GEN M2 = M, P2 = P;
     497           0 :   gerepileall(av, P ? 2: 1, &M2, &P2);
     498           0 :   RgM_replace(M, M2);
     499           0 :   if (P) RgM_replace(P, P2);
     500           0 : }
     501             : 
     502             : /* Lemma 9.14 */
     503             : static long
     504         408 : weakfrobenius_step1(GEN M, GEN P, long j0)
     505             : {
     506         408 :   pari_sp av = avma;
     507         408 :   long n = lg(M)-1, k, j;
     508        2136 :   for (j = j0; j < n; ++j)
     509             :   {
     510        1830 :     if (gequal0(gcoeff(M, j+1, j)))
     511             :     {
     512         816 :       for (k = j+2; k <= n; ++k)
     513         714 :         if (!gequal0(gcoeff(M,k,j))) break;
     514         258 :       if (k > n) return j;
     515         156 :       transS(M, P, k, j+1);
     516             :     }
     517        1728 :     transD(M, P, j+1, j, j+1);
     518             :     /* Now M[j+1,j] = 1 */
     519       20136 :     for (k = 1; k <= n; ++k)
     520       18408 :       if (k != j+1 && !gequal0(gcoeff(M,k,j))) /* zero M[k,j] */
     521             :       {
     522       12228 :         transL(M, P, gneg(gcoeff(M,k,j)), k, j+1);
     523       12228 :         gcoeff(M,k,j) = gen_0; /* avoid approximate 0 */
     524             :       }
     525        1728 :     if (gc_needed(av,1))
     526             :     {
     527           0 :       if (DEBUGMEM > 1)
     528           0 :         pari_warn(warnmem,"RgM_minpoly stage 1: j0=%ld, j=%ld", j0, j);
     529           0 :       gerepilemat2_inplace(av, M, P);
     530             :     }
     531             :   }
     532         306 :   return n;
     533             : }
     534             : 
     535             : static void
     536         408 : weakfrobenius_step2(GEN M, GEN P, long j)
     537             : {
     538         408 :   pari_sp av = avma;
     539         408 :   long i, k, n = lg(M)-1;
     540        2532 :   for(i=j; i>=2; i--)
     541             :   {
     542        4650 :     for(k=j+1; k<=n; k++)
     543        2526 :       if (!gequal0(gcoeff(M,i,k)))
     544         294 :         transL(M, P, gcoeff(M,i,k), i-1, k);
     545        2124 :     if (gc_needed(av,1))
     546             :     {
     547           0 :       if (DEBUGMEM > 1)
     548           0 :         pari_warn(warnmem,"RgM_minpoly stage 2: j=%ld, i=%ld", j, i);
     549           0 :       gerepilemat2_inplace(av, M, P);
     550             :     }
     551             :   }
     552         408 : }
     553             : 
     554             : static long
     555         408 : weakfrobenius_step3(GEN M, GEN P, long j0, long j)
     556             : {
     557         408 :   long i, k, n = lg(M)-1;
     558         408 :   if (j == n) return 0;
     559         102 :   if (gequal0(gcoeff(M, j0, j+1)))
     560             :   {
     561         372 :     for (k=j+2; k<=n; k++)
     562         288 :       if (!gequal0(gcoeff(M, j0, k))) break;
     563          84 :     if (k > n) return 0;
     564           0 :     transS(M, P, k, j+1);
     565             :   }
     566          18 :   transD(M, P, j0, j+1, j+1);
     567          18 :   for (i=j+2; i<=n; i++)
     568           0 :     if (!gequal0(gcoeff(M, j0, i)))
     569           0 :       transL(M, P, gcoeff(M, j0, i),j+1, i);
     570          18 :   return 1;
     571             : }
     572             : 
     573             : /* flag: 0 -> full Frobenius from , 1 -> weak Frobenius form */
     574             : static GEN
     575         306 : RgM_Frobenius(GEN M, long flag, GEN *pt_P, GEN *pt_v)
     576             : {
     577         306 :   pari_sp av = avma, av2, ltop;
     578         306 :   long n = lg(M)-1, eps, j0 = 1, nb = 0;
     579             :   GEN v, P;
     580         306 :   v = cgetg(n+1, t_VECSMALL);
     581         306 :   ltop = avma;
     582         306 :   P = pt_P ? matid(n): NULL;
     583         306 :   M = RgM_shallowcopy(M);
     584         306 :   av2 = avma;
     585        1020 :   while (j0 <= n)
     586             :   {
     587         408 :     long j = weakfrobenius_step1(M, P, j0);
     588         408 :     weakfrobenius_step2(M, P, j);
     589         408 :     eps = weakfrobenius_step3(M, P, j0, j);
     590         408 :     if (eps == 0)
     591             :     {
     592         390 :       v[++nb] = j0;
     593         390 :       if (flag == 0 && nb > 1 && !minpoly_dvdslice(M, v[nb-1], j0, j))
     594             :       {
     595           0 :         j = j0; j0 = v[nb-1]; nb -= 2;
     596           0 :         transL(M, P, gen_1, j, j0); /*lemma 9.18*/
     597             :       } else
     598         390 :         j0 = j+1;
     599             :     }
     600             :     else
     601          18 :       transS(M, P, j0, j+1); /*theorem 4*/
     602         408 :     if (gc_needed(av,1))
     603             :     {
     604           0 :       if (DEBUGMEM > 1)
     605           0 :         pari_warn(warnmem,"weakfrobenius j0=%ld",j0);
     606           0 :       gerepilemat2_inplace(av2, M, P);
     607             :     }
     608             :   }
     609         306 :   fixlg(v, nb+1);
     610         306 :   if (pt_v) *pt_v = v;
     611         306 :   gerepileall(pt_v ? ltop: av, P? 2: 1, &M, &P);
     612         306 :   if (pt_P) *pt_P = P;
     613         306 :   return M;
     614             : }
     615             : 
     616             : static GEN
     617          42 : RgM_minpoly(GEN M, long v)
     618             : {
     619          42 :   pari_sp av = avma;
     620             :   GEN V, W;
     621          42 :   M = RgM_Frobenius(M, 1, NULL, &V);
     622          42 :   W = minpoly_listpolslice(M, V, v);
     623          42 :   if (varncmp(v,gvar2(W)) >= 0)
     624           0 :     pari_err_PRIORITY("matfrobenius", M, "<=", v);
     625          42 :   return gerepileupto(av, RgX_normalize(glcm0(W, NULL)));
     626             : }
     627             : 
     628             : GEN
     629           0 : Frobeniusform(GEN V, long n)
     630             : {
     631             :   long i, j, k;
     632           0 :   GEN M = zeromatcopy(n,n);
     633           0 :   for (k=1,i=1;i<lg(V);i++,k++)
     634             :   {
     635           0 :     GEN  P = gel(V,i);
     636           0 :     long d = degpol(P);
     637           0 :     if (k+d-1 > n) pari_err_PREC("matfrobenius");
     638           0 :     for (j=0; j<d-1; j++, k++) gcoeff(M,k+1,k) = gen_1;
     639           0 :     for (j=0; j<d; j++) gcoeff(M,k-j,k) = gneg(gel(P, 1+d-j));
     640             :   }
     641           0 :   return M;
     642             : }
     643             : 
     644             : GEN
     645         264 : matfrobenius(GEN M, long flag, long v)
     646             : {
     647             :   long n;
     648         264 :   if (typ(M)!=t_MAT) pari_err_TYPE("matfrobenius",M);
     649         264 :   if (v < 0) v = 0;
     650         264 :   n = lg(M)-1;
     651         264 :   if (n && lgcols(M)!=n+1) pari_err_DIM("matfrobenius");
     652         264 :   if (flag > 2) pari_err_FLAG("matfrobenius");
     653         264 :   switch (flag)
     654             :   {
     655             :   case 0:
     656           6 :     return RgM_Frobenius(M, 0, NULL, NULL);
     657             :   case 1:
     658             :     {
     659           0 :       pari_sp av = avma;
     660             :       GEN V, W, F;
     661           0 :       F = RgM_Frobenius(M, 0, NULL, &V);
     662           0 :       W = minpoly_listpolslice(F, V, v);
     663           0 :       if (varncmp(v, gvar2(W)) >= 0)
     664           0 :         pari_err_PRIORITY("matfrobenius", M, "<=", v);
     665           0 :       return gerepileupto(av, W);
     666             :     }
     667             :   case 2:
     668             :     {
     669         258 :       GEN P, F, R = cgetg(3, t_VEC);
     670         258 :       F = RgM_Frobenius(M, 0, &P, NULL);
     671         258 :       gel(R,1) = F; gel(R,2) = P;
     672         258 :       return R;
     673             :     }
     674             :   default:
     675           0 :     pari_err_FLAG("matfrobenius");
     676             :   }
     677             :   return NULL; /*LCOV_EXCL_LINE*/
     678             : }
     679             : 
     680             : /*******************************************************************/
     681             : /*                                                                 */
     682             : /*                       MINIMAL POLYNOMIAL                        */
     683             : /*                                                                 */
     684             : /*******************************************************************/
     685             : 
     686             : static GEN
     687         792 : easymin(GEN x, long v)
     688             : {
     689         792 :   pari_sp ltop=avma;
     690             :   GEN G, R, dR;
     691         792 :   if (typ(x)==t_POLMOD && !issquarefree(gel(x,1)))
     692          48 :     return NULL;
     693         744 :   R = easychar(x, v);
     694         744 :   if (!R) return R;
     695         702 :   dR=RgX_deriv(R);
     696         702 :   if (!lgpol(dR)) {avma=ltop; return NULL;}
     697         702 :   G=RgX_gcd(R,dR);
     698         702 :   G=RgX_Rg_div(G,leading_coeff(G));
     699         702 :   return gerepileupto(ltop, RgX_div(R,G));
     700             : }
     701             : 
     702             : GEN
     703         816 : minpoly(GEN x, long v)
     704             : {
     705         816 :   pari_sp ltop=avma;
     706             :   GEN P;
     707         816 :   if (v<0) v = 0;
     708         816 :   if (typ(x)==t_FFELT)
     709             :   {
     710          24 :       GEN p1 = FpX_to_mod(FF_minpoly(x), FF_p_i(x));
     711          24 :       setvarn(p1,v); return gerepileupto(ltop,p1);
     712             :   }
     713             : 
     714         792 :   P = easymin(x,v);
     715         792 :   if (P) return P;
     716          90 :   if (typ(x)==t_POLMOD)
     717             :   {
     718          48 :     P = gcopy(RgXQ_minpoly_naive(gel(x,2), gel(x,1)));
     719          48 :     setvarn(P,v);
     720          48 :     return gerepileupto(ltop,P);
     721             :   }
     722          42 :   if (typ(x)!=t_MAT) pari_err_TYPE("minpoly",x);
     723          42 :   if (lg(x) == 1) return pol_1(v);
     724          42 :   return RgM_minpoly(x,v);
     725             : }
     726             : 
     727             : /*******************************************************************/
     728             : /*                                                                 */
     729             : /*                       HESSENBERG FORM                           */
     730             : /*                                                                 */
     731             : /*******************************************************************/
     732             : GEN
     733          24 : hess(GEN x)
     734             : {
     735          24 :   pari_sp av = avma;
     736          24 :   long lx = lg(x), m, i, j;
     737             : 
     738          24 :   if (typ(x) != t_MAT) pari_err_TYPE("hess",x);
     739          24 :   if (lx == 1) return cgetg(1,t_MAT);
     740          24 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     741             : 
     742          24 :   x = RgM_shallowcopy(x);
     743          54 :   for (m=2; m<lx-1; m++)
     744             :   {
     745          30 :     GEN t = NULL;
     746          30 :     for (i=m+1; i<lx; i++) { t = gcoeff(x,i,m-1); if (!gequal0(t)) break; }
     747          30 :     if (i == lx) continue;
     748          30 :     for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
     749          30 :     swap(gel(x,i), gel(x,m)); t = ginv(t);
     750             : 
     751         120 :     for (i=m+1; i<lx; i++)
     752             :     {
     753          90 :       GEN c = gcoeff(x,i,m-1);
     754          90 :       if (gequal0(c)) continue;
     755             : 
     756          90 :       c = gmul(c,t); gcoeff(x,i,m-1) = gen_0;
     757         510 :       for (j=m; j<lx; j++)
     758         420 :         gcoeff(x,i,j) = gsub(gcoeff(x,i,j), gmul(c,gcoeff(x,m,j)));
     759         720 :       for (j=1; j<lx; j++)
     760         630 :         gcoeff(x,j,m) = gadd(gcoeff(x,j,m), gmul(c,gcoeff(x,j,i)));
     761          90 :       if (gc_needed(av,2))
     762             :       {
     763           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
     764           0 :         gerepileall(av,2, &x, &t);
     765             :       }
     766             :     }
     767             :   }
     768          24 :   return gerepilecopy(av,x);
     769             : }
     770             : 
     771             : GEN
     772        8898 : Flm_hess(GEN x, ulong p)
     773             : {
     774        8898 :   long lx = lg(x), m, i, j;
     775        8898 :   if (lx == 1) return cgetg(1,t_MAT);
     776        8898 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     777             : 
     778        8898 :   x = Flm_copy(x);
     779       32076 :   for (m=2; m<lx-1; m++)
     780             :   {
     781       23178 :     ulong t = 0;
     782       23178 :     for (i=m+1; i<lx; i++) { t = ucoeff(x,i,m-1); if (t) break; }
     783       23178 :     if (i == lx) continue;
     784       14604 :     for (j=m-1; j<lx; j++) lswap(ucoeff(x,i,j), ucoeff(x,m,j));
     785       14604 :     swap(gel(x,i), gel(x,m)); t = Fl_inv(t, p);
     786             : 
     787      110526 :     for (i=m+1; i<lx; i++)
     788             :     {
     789       95922 :       ulong c = ucoeff(x,i,m-1);
     790       95922 :       if (!c) continue;
     791             : 
     792       56328 :       c = Fl_mul(c,t,p); ucoeff(x,i,m-1) = 0;
     793      646728 :       for (j=m; j<lx; j++)
     794      590400 :         ucoeff(x,i,j) = Fl_sub(ucoeff(x,i,j), Fl_mul(c,ucoeff(x,m,j), p), p);
     795      931800 :       for (j=1; j<lx; j++)
     796      875472 :         ucoeff(x,j,m) = Fl_add(ucoeff(x,j,m), Fl_mul(c,ucoeff(x,j,i), p), p);
     797             :     }
     798             :   }
     799        8898 :   return x;
     800             : }
     801             : GEN
     802          12 : FpM_hess(GEN x, GEN p)
     803             : {
     804          12 :   pari_sp av = avma;
     805          12 :   long lx = lg(x), m, i, j;
     806          12 :   if (lx == 1) return cgetg(1,t_MAT);
     807          12 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     808          12 :   if (lgefint(p) == 3)
     809             :   {
     810           0 :     ulong pp = p[2];
     811           0 :     x = Flm_hess(ZM_to_Flm(x, pp), pp);
     812           0 :     return gerepileupto(av, Flm_to_ZM(x));
     813             :   }
     814          12 :   x = RgM_shallowcopy(x);
     815          24 :   for (m=2; m<lx-1; m++)
     816             :   {
     817          12 :     GEN t = NULL;
     818          12 :     for (i=m+1; i<lx; i++) { t = gcoeff(x,i,m-1); if (signe(t)) break; }
     819          12 :     if (i == lx) continue;
     820          12 :     for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
     821          12 :     swap(gel(x,i), gel(x,m)); t = Fp_inv(t, p);
     822             : 
     823          30 :     for (i=m+1; i<lx; i++)
     824             :     {
     825          18 :       GEN c = gcoeff(x,i,m-1);
     826          18 :       if (!signe(c)) continue;
     827             : 
     828          18 :       c = Fp_mul(c,t, p); gcoeff(x,i,m-1) = gen_0;
     829          66 :       for (j=m; j<lx; j++)
     830          48 :         gcoeff(x,i,j) = Fp_sub(gcoeff(x,i,j), Fp_mul(c,gcoeff(x,m,j),p), p);
     831          90 :       for (j=1; j<lx; j++)
     832          72 :         gcoeff(x,j,m) = Fp_add(gcoeff(x,j,m), Fp_mul(c,gcoeff(x,j,i),p), p);
     833          18 :       if (gc_needed(av,2))
     834             :       {
     835           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
     836           0 :         gerepileall(av,2, &x, &t);
     837             :       }
     838             :     }
     839             :   }
     840          12 :   return gerepilecopy(av,x);
     841             : }
     842             : GEN
     843          18 : carhess(GEN x, long v)
     844             : {
     845             :   pari_sp av;
     846             :   long lx, r, i;
     847             :   GEN y, H;
     848             : 
     849          18 :   if ((H = easychar(x,v))) return H;
     850             : 
     851          18 :   lx = lg(x); av = avma; y = cgetg(lx+1, t_VEC);
     852          18 :   gel(y,1) = pol_1(v); H = hess(x);
     853          54 :   for (r = 1; r < lx; r++)
     854             :   {
     855          36 :     pari_sp av2 = avma;
     856          36 :     GEN z, a = gen_1, b = pol_0(v);
     857          48 :     for (i = r-1; i; i--)
     858             :     {
     859          18 :       a = gmul(a, gcoeff(H,i+1,i));
     860          18 :       if (gequal0(a)) break;
     861          12 :       b = RgX_add(b, RgX_Rg_mul(gel(y,i), gmul(a,gcoeff(H,i,r))));
     862             :     }
     863          72 :     z = RgX_sub(RgX_shift_shallow(gel(y,r), 1),
     864          72 :                 RgX_Rg_mul(gel(y,r), gcoeff(H,r,r)));
     865          36 :     gel(y,r+1) = gerepileupto(av2, RgX_sub(z, b)); /* (X - H[r,r])y[r] - b */
     866             :   }
     867          18 :   return fix_pol(av, gel(y,lx));
     868             : }
     869             : 
     870             : /* Bound for sup norm of charpoly(M/dM), M integral: let B = |M|oo / |dM|,
     871             :  *   s = max_k binomial(n,k) (kB^2)^(k/2),
     872             :  * return ceil(log2(s)) */
     873             : static double
     874        1890 : charpoly_bound(GEN M, GEN dM)
     875             : {
     876        1890 :   pari_sp av = avma;
     877        1890 :   GEN B = itor(ZM_supnorm(M), LOWDEFAULTPREC);
     878        1890 :   GEN s = real_0(LOWDEFAULTPREC), bin, B2;
     879        1890 :   long n = lg(M)-1, k;
     880             :   double d;
     881        1890 :   bin = gen_1;
     882        1890 :   if (dM) B = divri(B, dM);
     883        1890 :   B2 = sqrr(B);
     884        9372 :   for (k = n; k >= (n+1)>>1; k--)
     885             :   {
     886        7482 :     GEN t = mulri(powruhalf(mulur(k, B2), k), bin);
     887        7482 :     if (abscmprr(t, s) > 0) s = t;
     888        7482 :     bin = diviuexact(muliu(bin, k), n-k+1);
     889             :   }
     890        1890 :   d = dbllog2(s); avma = av; return ceil(d);
     891             : }
     892             : 
     893             : /* Return char_{M/d}(X) = d^(-n) char_M(dX) modulo p. Assume dp = d mod p. */
     894             : static GEN
     895        2550 : QM_charpoly_Flx(GEN M, ulong dp, ulong p)
     896             : {
     897        2550 :   pari_sp av = avma;
     898        2550 :   GEN H = Flm_charpoly_i(ZM_to_Flm(M,p), p);
     899        2550 :   if (dp) H = Flx_rescale(H, Fl_inv(dp,p), p);
     900        2550 :   return gerepileuptoleaf(av, H);
     901             : }
     902             : 
     903             : static int
     904        2550 : ZX_CRT(GEN *H, GEN Hp, GEN *q, ulong p, long bit)
     905             : {
     906        2550 :   if (!*H)
     907             :   {
     908        2094 :     *H = ZX_init_CRT(Hp, p, 0);
     909        2094 :     if (DEBUGLEVEL>5)
     910           0 :       err_printf("charpoly mod %lu, bound = 2^%ld\n", p, expu(p));
     911        2094 :     if (expu(p) > bit) return 1;
     912         288 :     *q = utoipos(p);
     913             :   }
     914             :   else
     915             :   {
     916         456 :     int stable = ZX_incremental_CRT(H, Hp, q,p);
     917         456 :     if (DEBUGLEVEL>5)
     918           0 :       err_printf("charpoly mod %lu (stable=%ld), bound = 2^%ld\n",
     919             :                  p, stable, expi(*q));
     920         456 :     if (stable && expi(*q) > bit) return 1;
     921             :   }
     922         456 :   return 0;
     923             : }
     924             : 
     925             : /* Assume M a square ZM, dM integer. Return charpoly(M / dM) in Z[X] */
     926             : static GEN
     927        2136 : QM_charpoly_ZX_i(GEN M, GEN dM, long bit)
     928             : {
     929        2136 :   long n = lg(M)-1;
     930        2136 :   GEN q = NULL, H = NULL;
     931             :   forprime_t S;
     932             :   ulong p;
     933        2136 :   if (!n) return pol_1(0);
     934             : 
     935        2094 :   if (bit < 0) bit = (long)charpoly_bound(M, dM) + 1;
     936        2094 :   if (DEBUGLEVEL>5) err_printf("ZM_charpoly: bit-bound 2^%ld\n", bit);
     937        2094 :   init_modular_big(&S);
     938        2094 :   while ((p = u_forprime_next(&S)))
     939             :   {
     940        2550 :     ulong dMp = 0;
     941             :     GEN Hp;
     942        2550 :     if (dM && !(dMp = umodiu(dM, p))) continue;
     943        2550 :     Hp = QM_charpoly_Flx(M, dMp, p);
     944        2550 :     if (ZX_CRT(&H, Hp, &q,p, bit)) break;
     945             :   }
     946        2094 :   if (!p) pari_err_OVERFLOW("charpoly [ran out of primes]");
     947        2094 :   return H;
     948             : }
     949             : GEN
     950         246 : QM_charpoly_ZX_bound(GEN M, long bit)
     951             : {
     952         246 :   pari_sp av = avma;
     953         246 :   GEN dM; M = Q_remove_denom(M, &dM);
     954         246 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, bit));
     955             : }
     956             : GEN
     957         498 : QM_charpoly_ZX(GEN M)
     958             : {
     959         498 :   pari_sp av = avma;
     960         498 :   GEN dM; M = Q_remove_denom(M, &dM);
     961         498 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, -1));
     962             : }
     963             : GEN
     964        1392 : ZM_charpoly(GEN M)
     965             : {
     966        1392 :   pari_sp av = avma;
     967        1392 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, NULL, -1));
     968             : }
     969             : 
     970             : /*******************************************************************/
     971             : /*                                                                 */
     972             : /*        CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)        */
     973             : /*                                                                 */
     974             : /*******************************************************************/
     975             : GEN
     976         810 : carberkowitz(GEN x, long v)
     977             : {
     978             :   long lx, i, j, k, r;
     979             :   GEN V, S, C, Q;
     980             :   pari_sp av0, av;
     981         810 :   if ((V = easychar(x,v))) return V;
     982         810 :   lx = lg(x); av0 = avma;
     983         810 :   V = cgetg(lx+1, t_VEC);
     984         810 :   S = cgetg(lx+1, t_VEC);
     985         810 :   C = cgetg(lx+1, t_VEC);
     986         810 :   Q = cgetg(lx+1, t_VEC);
     987         810 :   av = avma;
     988         810 :   gel(C,1) = gen_m1;
     989         810 :   gel(V,1) = gen_m1;
     990         810 :   for (i=2;i<=lx; i++) gel(C,i) = gel(Q,i) = gel(S,i) = gel(V,i) = gen_0;
     991         810 :   gel(V,2) = gcoeff(x,1,1);
     992        5136 :   for (r = 2; r < lx; r++)
     993             :   {
     994             :     pari_sp av2;
     995             :     GEN t;
     996             : 
     997        4326 :     for (i = 1; i < r; i++) gel(S,i) = gcoeff(x,i,r);
     998        4326 :     gel(C,2) = gcoeff(x,r,r);
     999       19212 :     for (i = 1; i < r-1; i++)
    1000             :     {
    1001       14886 :       av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1002       14886 :       for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1003       14886 :       gel(C,i+2) = gerepileupto(av2, t);
    1004      114138 :       for (j = 1; j < r; j++)
    1005             :       {
    1006       99252 :         av2 = avma; t = gmul(gcoeff(x,j,1), gel(S,1));
    1007       99252 :         for (k = 2; k < r; k++) t = gadd(t, gmul(gcoeff(x,j,k), gel(S,k)));
    1008       99252 :         gel(Q,j) = gerepileupto(av2, t);
    1009             :       }
    1010       14886 :       for (j = 1; j < r; j++) gel(S,j) = gel(Q,j);
    1011             :     }
    1012        4326 :     av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1013        4326 :     for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1014        4326 :     gel(C,r+1) = gerepileupto(av2, t);
    1015        4326 :     if (gc_needed(av0,1))
    1016             :     {
    1017           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"carberkowitz");
    1018           0 :       gerepileall(av, 2, &C, &V);
    1019             :     }
    1020       32190 :     for (i = 1; i <= r+1; i++)
    1021             :     {
    1022       27864 :       av2 = avma; t = gmul(gel(C,i), gel(V,1));
    1023      115914 :       for (j = 2; j <= minss(r,i); j++)
    1024       88050 :         t = gadd(t, gmul(gel(C,i+1-j), gel(V,j)));
    1025       27864 :       gel(Q,i) = gerepileupto(av2, t);
    1026             :     }
    1027        4326 :     for (i = 1; i <= r+1; i++) gel(V,i) = gel(Q,i);
    1028             :   }
    1029         810 :   V = RgV_to_RgX(vecreverse(V), v); /* not gtopoly: fail if v > gvar(V) */
    1030         810 :   V = odd(lx)? gcopy(V): RgX_neg(V);
    1031         810 :   return fix_pol(av0, V);
    1032             : }
    1033             : 
    1034             : /*******************************************************************/
    1035             : /*                                                                 */
    1036             : /*                            NORMS                                */
    1037             : /*                                                                 */
    1038             : /*******************************************************************/
    1039             : GEN
    1040      155988 : gnorm(GEN x)
    1041             : {
    1042             :   pari_sp av;
    1043             :   long lx, i;
    1044             :   GEN y;
    1045             : 
    1046      155988 :   switch(typ(x))
    1047             :   {
    1048         624 :     case t_INT:  return sqri(x);
    1049       16374 :     case t_REAL: return sqrr(x);
    1050          90 :     case t_FRAC: return sqrfrac(x);
    1051      135126 :     case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
    1052         204 :     case t_QUAD:    av = avma; return gerepileupto(av, quadnorm(x));
    1053             : 
    1054          12 :     case t_POL: case t_SER: case t_RFRAC: av = avma;
    1055          12 :       return gerepileupto(av, greal(gmul(gconj(x),x)));
    1056             : 
    1057             :     case t_FFELT:
    1058          24 :       y = cgetg(3, t_INTMOD);
    1059          24 :       gel(y,1) = FF_p(x);
    1060          24 :       gel(y,2) = FF_norm(x); return y;
    1061             : 
    1062             :     case t_POLMOD:
    1063             :     {
    1064        3534 :       GEN T = gel(x,1), a = gel(x,2);
    1065        3534 :       if (typ(a) != t_POL || varn(a) != varn(T)) return gpowgs(a, degpol(T));
    1066        3348 :       return RgXQ_norm(a, T);
    1067             :     }
    1068             :     case t_VEC: case t_COL: case t_MAT:
    1069           0 :       y = cgetg_copy(x, &lx);
    1070           0 :       for (i=1; i<lx; i++) gel(y,i) = gnorm(gel(x,i));
    1071           0 :       return y;
    1072             :   }
    1073           0 :   pari_err_TYPE("gnorm",x);
    1074             :   return NULL; /* LCOV_EXCL_LINE */
    1075             : }
    1076             : 
    1077             : /* return |q|^2, complex modulus */
    1078             : static GEN
    1079          24 : cxquadnorm(GEN q, long prec)
    1080             : {
    1081          24 :   GEN X = gel(q,1), c = gel(X,2); /* (1-D)/4, -D/4 */
    1082          24 :   if (signe(c) > 0) return quadnorm(q); /* imaginary */
    1083          18 :   if (!prec) pari_err_TYPE("gnorml2", q);
    1084           6 :   return sqrr(quadtofp(q, prec));
    1085             : }
    1086             : 
    1087             : static GEN
    1088    10233540 : gnorml2_i(GEN x, long prec)
    1089             : {
    1090             :   pari_sp av;
    1091             :   long i, lx;
    1092             :   GEN s;
    1093             : 
    1094    10233540 :   switch(typ(x))
    1095             :   {
    1096     5785662 :     case t_INT:  return sqri(x);
    1097     2321736 :     case t_REAL: return sqrr(x);
    1098           6 :     case t_FRAC: return sqrfrac(x);
    1099      352560 :     case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
    1100          18 :     case t_QUAD:    av = avma; return gerepileupto(av, cxquadnorm(x,prec));
    1101             : 
    1102        6126 :     case t_POL: lx = lg(x)-1; x++; break;
    1103             : 
    1104             :     case t_VEC:
    1105             :     case t_COL:
    1106     1767432 :     case t_MAT: lx = lg(x); break;
    1107             : 
    1108           0 :     default: pari_err_TYPE("gnorml2",x);
    1109             :       return NULL; /* LCOV_EXCL_LINE */
    1110             :   }
    1111     1773558 :   if (lx == 1) return gen_0;
    1112     1773558 :   av = avma;
    1113     1773558 :   s = gnorml2(gel(x,1));
    1114     8436180 :   for (i=2; i<lx; i++)
    1115             :   {
    1116     6662622 :     s = gadd(s, gnorml2(gel(x,i)));
    1117     6662622 :     if (gc_needed(av,1))
    1118             :     {
    1119           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnorml2");
    1120           0 :       s = gerepileupto(av, s);
    1121             :     }
    1122             :   }
    1123     1773558 :   return gerepileupto(av,s);
    1124             : }
    1125             : GEN
    1126    10233516 : gnorml2(GEN x) { return gnorml2_i(x, 0); }
    1127             : 
    1128             : static GEN pnormlp(GEN,GEN,long);
    1129             : static GEN
    1130          54 : pnormlpvec(long i0, GEN x, GEN p, long prec)
    1131             : {
    1132          54 :   pari_sp av = avma;
    1133          54 :   long i, lx = lg(x);
    1134          54 :   GEN s = gen_0;
    1135         192 :   for (i=i0; i<lx; i++)
    1136             :   {
    1137         138 :     s = gadd(s, pnormlp(gel(x,i),p,prec));
    1138         138 :     if (gc_needed(av,1))
    1139             :     {
    1140           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnormlp, i = %ld", i);
    1141           0 :       s = gerepileupto(av, s);
    1142             :     }
    1143             :   }
    1144          54 :   return s;
    1145             : }
    1146             : /* (||x||_p)^p */
    1147             : static GEN
    1148         168 : pnormlp(GEN x, GEN p, long prec)
    1149             : {
    1150         168 :   switch(typ(x))
    1151             :   {
    1152         102 :     case t_INT: case t_REAL: x = mpabs(x); break;
    1153           0 :     case t_FRAC: x = absfrac(x); break;
    1154          12 :     case t_COMPLEX: case t_QUAD: x = gabs(x,prec); break;
    1155           6 :     case t_POL: return pnormlpvec(2, x, p, prec);
    1156          48 :     case t_VEC: case t_COL: case t_MAT: return pnormlpvec(1, x, p, prec);
    1157           0 :     default: pari_err_TYPE("gnormlp",x);
    1158             :   }
    1159         114 :   return gpow(x, p, prec);
    1160             : }
    1161             : 
    1162             : GEN
    1163         282 : gnormlp(GEN x, GEN p, long prec)
    1164             : {
    1165         282 :   pari_sp av = avma;
    1166         282 :   if (!p || (typ(p) == t_INFINITY && inf_get_sign(p) > 0))
    1167         144 :     return gsupnorm(x, prec);
    1168         138 :   if (gsigne(p) <= 0) pari_err_DOMAIN("normlp", "p", "<=", gen_0, p);
    1169         132 :   if (is_scalar_t(typ(x))) return gabs(x, prec);
    1170          78 :   if (typ(p) == t_INT)
    1171             :   {
    1172          54 :     ulong pp = itou_or_0(p);
    1173          54 :     switch(pp)
    1174             :     {
    1175          24 :       case 1: return gnorml1(x, prec);
    1176          24 :       case 2: x = gnorml2_i(x, prec); break;
    1177           6 :       default: x = pnormlp(x, p, prec); break;
    1178             :     }
    1179          30 :     if (pp && typ(x) == t_INT && Z_ispowerall(x, pp, &x))
    1180           6 :       return gerepileuptoleaf(av, x);
    1181          24 :     if (pp == 2) return gerepileupto(av, gsqrt(x, prec));
    1182             :   }
    1183             :   else
    1184          24 :     x = pnormlp(x, p, prec);
    1185          24 :   x = gpow(x, ginv(p), prec);
    1186          24 :   return gerepileupto(av, x);
    1187             : }
    1188             : 
    1189             : GEN
    1190         144 : gnorml1(GEN x,long prec)
    1191             : {
    1192         144 :   pari_sp av = avma;
    1193             :   long lx,i;
    1194             :   GEN s;
    1195         144 :   switch(typ(x))
    1196             :   {
    1197          84 :     case t_INT: case t_REAL: return mpabs(x);
    1198           0 :     case t_FRAC: return absfrac(x);
    1199             : 
    1200             :     case t_COMPLEX: case t_QUAD:
    1201          12 :       return gabs(x,prec);
    1202             : 
    1203             :     case t_POL:
    1204           6 :       lx = lg(x); s = gen_0;
    1205           6 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1206           6 :       break;
    1207             : 
    1208             :     case t_VEC: case t_COL: case t_MAT:
    1209          42 :       lx = lg(x); s = gen_0;
    1210          42 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1211          42 :       break;
    1212             : 
    1213           0 :     default: pari_err_TYPE("gnorml1",x);
    1214             :       return NULL; /* LCOV_EXCL_LINE */
    1215             :   }
    1216          48 :   return gerepileupto(av, s);
    1217             : }
    1218             : /* As gnorml1, except for t_QUAD and t_COMPLEX: |x + wy| := |x| + |y|
    1219             :  * Still a norm of R-vector spaces, and can be cheaply computed without
    1220             :  * square roots */
    1221             : GEN
    1222           0 : gnorml1_fake(GEN x)
    1223             : {
    1224           0 :   pari_sp av = avma;
    1225             :   long lx, i;
    1226             :   GEN s;
    1227           0 :   switch(typ(x))
    1228             :   {
    1229           0 :     case t_INT: case t_REAL: return mpabs(x);
    1230           0 :     case t_FRAC: return absfrac(x);
    1231             : 
    1232             :     case t_COMPLEX:
    1233           0 :       s = gadd(gnorml1_fake(gel(x,1)), gnorml1_fake(gel(x,2)));
    1234           0 :       break;
    1235             :     case t_QUAD:
    1236           0 :       s = gadd(gnorml1_fake(gel(x,2)), gnorml1_fake(gel(x,3)));
    1237           0 :       break;
    1238             : 
    1239             :     case t_POL:
    1240           0 :       lx = lg(x); s = gen_0;
    1241           0 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1242           0 :       break;
    1243             : 
    1244             :     case t_VEC: case t_COL: case t_MAT:
    1245           0 :       lx = lg(x); s = gen_0;
    1246           0 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1247           0 :       break;
    1248             : 
    1249           0 :     default: pari_err_TYPE("gnorml1_fake",x);
    1250             :       return NULL; /* LCOV_EXCL_LINE */
    1251             :   }
    1252           0 :   return gerepileupto(av, s);
    1253             : }
    1254             : 
    1255             : static void
    1256       43548 : store(GEN z, GEN *m) { if (!*m || gcmp(z, *m) > 0) *m = z; }
    1257             : /* compare |x| to *m or |x|^2 to *msq, whichever is easiest, and update
    1258             :  * the pointed value if x is larger */
    1259             : void
    1260       48084 : gsupnorm_aux(GEN x, GEN *m, GEN *msq, long prec)
    1261             : {
    1262             :   long i, lx;
    1263             :   GEN z;
    1264       48084 :   switch(typ(x))
    1265             :   {
    1266        3990 :     case t_COMPLEX: z = cxnorm(x); store(z, msq); return;
    1267           6 :     case t_QUAD:  z = cxquadnorm(x,prec); store(z, msq); return;
    1268       39498 :     case t_INT: case t_REAL: z = mpabs(x); store(z,m); return;
    1269          54 :     case t_FRAC: z = absfrac(x); store(z,m); return;
    1270             : 
    1271        2688 :     case t_POL: lx = lg(x)-1; x++; break;
    1272             : 
    1273             :     case t_VEC:
    1274             :     case t_COL:
    1275        1848 :     case t_MAT: lx = lg(x); break;
    1276             : 
    1277           0 :     default: pari_err_TYPE("gsupnorm",x);
    1278             :       return; /* LCOV_EXCL_LINE */
    1279             :   }
    1280        4536 :   for (i=1; i<lx; i++) gsupnorm_aux(gel(x,i), m, msq, prec);
    1281             : }
    1282             : GEN
    1283        4530 : gsupnorm(GEN x, long prec)
    1284             : {
    1285        4530 :   GEN m = NULL, msq = NULL;
    1286        4530 :   pari_sp av = avma;
    1287        4530 :   gsupnorm_aux(x, &m, &msq, prec);
    1288             :   /* now set m = max (m, sqrt(msq)) */
    1289        4530 :   if (msq) {
    1290         516 :     msq = gsqrt(msq, prec);
    1291         516 :     if (!m || gcmp(m, msq) < 0) m = msq;
    1292        4014 :   } else if (!m) m = gen_0;
    1293        4530 :   return gerepilecopy(av, m);
    1294             : }
    1295             : 
    1296             : /*******************************************************************/
    1297             : /*                                                                 */
    1298             : /*                            TRACES                               */
    1299             : /*                                                                 */
    1300             : /*******************************************************************/
    1301             : GEN
    1302          24 : matcompanion(GEN x)
    1303             : {
    1304          24 :   long n = degpol(x), j;
    1305             :   GEN y, c;
    1306             : 
    1307          24 :   if (typ(x)!=t_POL) pari_err_TYPE("matcompanion",x);
    1308          24 :   if (!signe(x)) pari_err_DOMAIN("matcompanion","polynomial","=",gen_0,x);
    1309          18 :   if (n == 0) return cgetg(1, t_MAT);
    1310             : 
    1311          18 :   y = cgetg(n+1,t_MAT);
    1312          18 :   for (j=1; j < n; j++) gel(y,j) = col_ei(n, j+1);
    1313          18 :   c = cgetg(n+1,t_COL); gel(y,n) = c;
    1314          18 :   if (gequal1(gel(x, n+2)))
    1315          12 :     for (j=1; j<=n; j++) gel(c,j) = gneg(gel(x,j+1));
    1316             :   else
    1317             :   { /* not monic. Hardly ever used */
    1318           6 :     pari_sp av = avma;
    1319           6 :     GEN d = gclone(gneg(gel(x,n+2)));
    1320           6 :     avma = av;
    1321           6 :     for (j=1; j<=n; j++) gel(c,j) = gdiv(gel(x,j+1), d);
    1322           6 :     gunclone(d);
    1323             :   }
    1324          18 :   return y;
    1325             : }
    1326             : 
    1327             : GEN
    1328       90666 : gtrace(GEN x)
    1329             : {
    1330             :   pari_sp av;
    1331       90666 :   long i, lx, tx = typ(x);
    1332             :   GEN y, z;
    1333             : 
    1334       90666 :   switch(tx)
    1335             :   {
    1336             :     case t_INT: case t_REAL: case t_FRAC:
    1337         636 :       return gmul2n(x,1);
    1338             : 
    1339             :     case t_COMPLEX:
    1340       89406 :       return gmul2n(gel(x,1),1);
    1341             : 
    1342             :     case t_QUAD:
    1343          24 :       y = gel(x,1);
    1344          24 :       if (!gequal0(gel(y,3)))
    1345             :       { /* assume quad. polynomial is either x^2 + d or x^2 - x + d */
    1346          24 :         av = avma;
    1347          24 :         return gerepileupto(av, gadd(gel(x,3), gmul2n(gel(x,2),1)));
    1348             :       }
    1349           0 :       return gmul2n(gel(x,2),1);
    1350             : 
    1351             :     case t_POL:
    1352           6 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1353           6 :       for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1354           6 :       return normalizepol_lg(y, lx);
    1355             : 
    1356             :     case t_SER:
    1357          12 :       if (ser_isexactzero(x)) return gcopy(x);
    1358           6 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1359           6 :       for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1360           6 :       return normalize(y);
    1361             : 
    1362             :     case t_POLMOD:
    1363         342 :       y = gel(x,1); z = gel(x,2);
    1364         342 :       if (typ(z) != t_POL || varn(y) != varn(z)) return gmulsg(degpol(y), z);
    1365         222 :       av = avma;
    1366         222 :       return gerepileupto(av, quicktrace(z, polsym(y, degpol(y)-1)));
    1367             : 
    1368             :     case t_FFELT:
    1369          24 :       y=cgetg(3, t_INTMOD);
    1370          24 :       gel(y,1) = FF_p(x);
    1371          24 :       gel(y,2) = FF_trace(x);
    1372          24 :       return y;
    1373             : 
    1374             : 
    1375             :     case t_RFRAC:
    1376           6 :       return gadd(x, gconj(x));
    1377             : 
    1378             :     case t_VEC: case t_COL:
    1379           0 :       y = cgetg_copy(x, &lx);
    1380           0 :       for (i=1; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1381           0 :       return y;
    1382             : 
    1383             :     case t_MAT:
    1384         210 :       lx = lg(x); if (lx == 1) return gen_0;
    1385             :       /*now lx >= 2*/
    1386         204 :       if (lx != lgcols(x)) pari_err_DIM("gtrace");
    1387         198 :       av = avma; return gerepileupto(av, mattrace(x));
    1388             :   }
    1389           0 :   pari_err_TYPE("gtrace",x);
    1390             :   return NULL; /* LCOV_EXCL_LINE */
    1391             : }
    1392             : 
    1393             : /* Cholesky decomposition for positive definite matrix a
    1394             :  * [GTM138, Algo 2.7.6, matrix Q]
    1395             :  * If a is not positive definite return NULL. */
    1396             : GEN
    1397         858 : qfgaussred_positive(GEN a)
    1398             : {
    1399         858 :   pari_sp av = avma;
    1400             :   GEN b;
    1401         858 :   long i,j,k, n = lg(a);
    1402             : 
    1403         858 :   if (typ(a)!=t_MAT) pari_err_TYPE("qfgaussred_positive",a);
    1404         858 :   if (n == 1) return cgetg(1, t_MAT);
    1405         852 :   if (lgcols(a)!=n) pari_err_DIM("qfgaussred_positive");
    1406         852 :   b = cgetg(n,t_MAT);
    1407        5790 :   for (j=1; j<n; j++)
    1408             :   {
    1409        4938 :     GEN p1=cgetg(n,t_COL), p2=gel(a,j);
    1410             : 
    1411        4938 :     gel(b,j) = p1;
    1412        4938 :     for (i=1; i<=j; i++) gel(p1,i) = gel(p2,i);
    1413        4938 :     for (   ; i<n ; i++) gel(p1,i) = gen_0;
    1414             :   }
    1415        5790 :   for (k=1; k<n; k++)
    1416             :   {
    1417        4938 :     GEN bk, p = gcoeff(b,k,k), invp;
    1418        4938 :     if (gsigne(p)<=0) { avma = av; return NULL; } /* not positive definite */
    1419        4938 :     invp = ginv(p);
    1420        4938 :     bk = row(b, k);
    1421        4938 :     for (i=k+1; i<n; i++) gcoeff(b,k,i) = gmul(gel(bk,i), invp);
    1422       25398 :     for (i=k+1; i<n; i++)
    1423             :     {
    1424       20460 :       GEN c = gel(bk, i);
    1425      115818 :       for (j=i; j<n; j++)
    1426       95358 :         gcoeff(b,i,j) = gsub(gcoeff(b,i,j), gmul(c,gcoeff(b,k,j)));
    1427             :     }
    1428        4938 :     if (gc_needed(av,1))
    1429             :     {
    1430           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfgaussred_positive");
    1431           0 :       b=gerepilecopy(av,b);
    1432             :     }
    1433             :   }
    1434         852 :   return gerepilecopy(av,b);
    1435             : }
    1436             : 
    1437             : /* Maximal pivot strategy: x is a suitable pivot if it is non zero and either
    1438             :  * - an exact type, or
    1439             :  * - it is maximal among remaining non-zero (t_REAL) pivots */
    1440             : static int
    1441       20172 : suitable(GEN x, long k, GEN *pp, long *pi)
    1442             : {
    1443       20172 :   long t = typ(x);
    1444       20172 :   switch(t)
    1445             :   {
    1446        3876 :     case t_INT: return signe(x) != 0;
    1447       16164 :     case t_FRAC: return 1;
    1448             :     case t_REAL: {
    1449         132 :       GEN p = *pp;
    1450         132 :       if (signe(x) && (!p || abscmprr(p, x) < 0)) { *pp = x; *pi = k; }
    1451         132 :       return 0;
    1452             :     }
    1453           0 :     default: return !gequal0(x);
    1454             :   }
    1455             : }
    1456             : 
    1457             : /* Gauss reduction (arbitrary symetric matrix, only the part above the
    1458             :  * diagonal is considered). If signature is non-zero, return only the
    1459             :  * signature, in which case gsigne() should be defined for elements of a. */
    1460             : static GEN
    1461        3654 : gaussred(GEN a, long signature)
    1462             : {
    1463             :   GEN r, ak, al;
    1464             :   pari_sp av, av1;
    1465        3654 :   long n = lg(a), i, j, k, l, sp, sn, t;
    1466             : 
    1467        3654 :   if (typ(a) != t_MAT) pari_err_TYPE("gaussred",a);
    1468        3654 :   if (n == 1) return signature? mkvec2(gen_0, gen_0): cgetg(1, t_MAT);
    1469        3654 :   if (lgcols(a) != n) pari_err_DIM("gaussred");
    1470        3654 :   n--;
    1471             : 
    1472        3654 :   av = avma;
    1473        3654 :   r = const_vecsmall(n, 1);
    1474        3654 :   av1= avma;
    1475        3654 :   a = RgM_shallowcopy(a);
    1476        3654 :   t = n; sp = sn = 0;
    1477       27360 :   while (t)
    1478             :   {
    1479       20052 :     long pind = 0;
    1480       20052 :     GEN invp, p = NULL;
    1481       20052 :     k=1; while (k<=n && (!r[k] || !suitable(gcoeff(a,k,k), k, &p, &pind))) k++;
    1482       20052 :     if (k > n && p) k = pind;
    1483       20052 :     if (k <= n)
    1484             :     {
    1485       20046 :       p = gcoeff(a,k,k); invp = ginv(p); /* != 0 */
    1486       20046 :       if (signature) { /* skip if (!signature): gsigne may fail ! */
    1487       20004 :         if (gsigne(p) > 0) sp++; else sn++;
    1488             :       }
    1489       20046 :       r[k] = 0; t--;
    1490       20046 :       ak = row(a, k);
    1491      141156 :       for (i=1; i<=n; i++)
    1492      121110 :         gcoeff(a,k,i) = r[i]? gmul(gcoeff(a,k,i), invp): gen_0;
    1493             : 
    1494      141156 :       for (i=1; i<=n; i++) if (r[i])
    1495             :       {
    1496       50520 :         GEN c = gel(ak,i); /* - p * a[k,i] */
    1497       50520 :         if (gequal0(c)) continue;
    1498      381540 :         for (j=1; j<=n; j++) if (r[j])
    1499      203970 :           gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
    1500             :       }
    1501       20046 :       gcoeff(a,k,k) = p;
    1502       20046 :       if (gc_needed(av1,1))
    1503             :       {
    1504           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"gaussred (t = %ld)", t);
    1505           0 :         a = gerepilecopy(av1, a);
    1506             :       }
    1507             :     }
    1508             :     else
    1509             :     { /* all remaining diagonal coeffs are currently 0 */
    1510           6 :       for (k=1; k<=n; k++) if (r[k])
    1511             :       {
    1512           6 :         l=k+1; while (l<=n && (!r[l] || !suitable(gcoeff(a,k,l), l, &p, &pind))) l++;
    1513           6 :         if (l > n && p) l = pind;
    1514           6 :         if (l > n) continue;
    1515             : 
    1516           6 :         p = gcoeff(a,k,l); invp = ginv(p);
    1517           6 :         sp++; sn++;
    1518           6 :         r[k] = r[l] = 0; t -= 2;
    1519           6 :         ak = row(a, k);
    1520           6 :         al = row(a, l);
    1521          30 :         for (i=1; i<=n; i++) if (r[i])
    1522             :         {
    1523          12 :           gcoeff(a,k,i) = gmul(gcoeff(a,k,i), invp);
    1524          12 :           gcoeff(a,l,i) = gmul(gcoeff(a,l,i), invp);
    1525             :         } else {
    1526          12 :           gcoeff(a,k,i) = gen_0;
    1527          12 :           gcoeff(a,l,i) = gen_0;
    1528             :         }
    1529             : 
    1530          30 :         for (i=1; i<=n; i++) if (r[i])
    1531             :         { /* c = a[k,i] * p, d = a[l,i] * p; */
    1532          12 :           GEN c = gel(ak,i), d = gel(al,i);
    1533          60 :           for (j=1; j<=n; j++) if (r[j])
    1534          72 :             gcoeff(a,i,j) = gsub(gcoeff(a,i,j),
    1535          24 :                                  gadd(gmul(gcoeff(a,l,j), c),
    1536          24 :                                       gmul(gcoeff(a,k,j), d)));
    1537             :         }
    1538          30 :         for (i=1; i<=n; i++) if (r[i])
    1539             :         {
    1540          12 :           GEN c = gcoeff(a,k,i), d = gcoeff(a,l,i);
    1541          12 :           gcoeff(a,k,i) = gadd(c, d);
    1542          12 :           gcoeff(a,l,i) = gsub(c, d);
    1543             :         }
    1544           6 :         gcoeff(a,k,l) = gen_1;
    1545           6 :         gcoeff(a,l,k) = gen_m1;
    1546           6 :         gcoeff(a,k,k) = gmul2n(p,-1);
    1547           6 :         gcoeff(a,l,l) = gneg(gcoeff(a,k,k));
    1548           6 :         if (gc_needed(av1,1))
    1549             :         {
    1550           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gaussred");
    1551           0 :           a = gerepilecopy(av1, a);
    1552             :         }
    1553           6 :         break;
    1554             :       }
    1555           6 :       if (k > n) break;
    1556             :     }
    1557             :   }
    1558        3654 :   if (!signature) return gerepilecopy(av, a);
    1559        3642 :   avma = av; return mkvec2s(sp, sn);
    1560             : }
    1561             : 
    1562             : GEN
    1563          12 : qfgaussred(GEN a) { return gaussred(a,0); }
    1564             : 
    1565             : GEN
    1566        3642 : qfsign(GEN a) { return gaussred(a,1); }
    1567             : 
    1568             : /* x -= s(y+u*x) */
    1569             : /* y += s(x-u*y), simultaneously */
    1570             : static void
    1571        4080 : rot(GEN x, GEN y, GEN s, GEN u) {
    1572        4080 :   GEN x1 = subrr(x, mulrr(s,addrr(y,mulrr(u,x))));
    1573        4080 :   GEN y1 = addrr(y, mulrr(s,subrr(x,mulrr(u,y))));
    1574        4080 :   affrr(x1,x);
    1575        4080 :   affrr(y1,y);
    1576        4080 : }
    1577             : 
    1578             : /* Diagonalization of a REAL symetric matrix. Return a vector [L, r]:
    1579             :  * L = vector of eigenvalues
    1580             :  * r = matrix of eigenvectors */
    1581             : GEN
    1582           6 : jacobi(GEN a, long prec)
    1583             : {
    1584             :   pari_sp av1;
    1585           6 :   long de, e, e1, e2, i, j, p, q, l = lg(a);
    1586             :   GEN c, ja, L, r, L2, r2, unr;
    1587             : 
    1588           6 :   if (typ(a) != t_MAT) pari_err_TYPE("jacobi",a);
    1589           6 :   ja = cgetg(3,t_VEC);
    1590           6 :   L = cgetg(l,t_COL); gel(ja,1) = L;
    1591           6 :   r = cgetg(l,t_MAT); gel(ja,2) = r;
    1592           6 :   if (l == 1) return ja;
    1593           6 :   if (lgcols(a) != l) pari_err_DIM("jacobi");
    1594             : 
    1595           6 :   e1 = HIGHEXPOBIT-1;
    1596          42 :   for (j=1; j<l; j++)
    1597             :   {
    1598          36 :     GEN z = gtofp(gcoeff(a,j,j), prec);
    1599          36 :     gel(L,j) = z;
    1600          36 :     e = expo(z); if (e < e1) e1 = e;
    1601             :   }
    1602          42 :   for (j=1; j<l; j++)
    1603             :   {
    1604          36 :     gel(r,j) = cgetg(l,t_COL);
    1605          36 :     for (i=1; i<l; i++) gcoeff(r,i,j) = utor(i==j? 1: 0, prec);
    1606             :   }
    1607           6 :   av1 = avma;
    1608             : 
    1609           6 :   e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1610           6 :   c = cgetg(l,t_MAT);
    1611          42 :   for (j=1; j<l; j++)
    1612             :   {
    1613          36 :     gel(c,j) = cgetg(j,t_COL);
    1614         126 :     for (i=1; i<j; i++)
    1615             :     {
    1616          90 :       GEN z = gtofp(gcoeff(a,i,j), prec);
    1617          90 :       gcoeff(c,i,j) = z;
    1618          90 :       if (!signe(z)) continue;
    1619          90 :       e = expo(z); if (e > e2) { e2 = e; p = i; q = j; }
    1620             :     }
    1621             :   }
    1622           6 :   a = c; unr = real_1(prec);
    1623           6 :   de = prec2nbits(prec);
    1624             : 
    1625             :  /* e1 = min expo(a[i,i])
    1626             :   * e2 = max expo(a[i,j]), i != j */
    1627         420 :   while (e1-e2 < de)
    1628             :   {
    1629         408 :     pari_sp av2 = avma;
    1630             :     GEN x, y, t, c, s, u;
    1631             :     /* compute attached rotation in the plane formed by basis vectors number
    1632             :      * p and q */
    1633         408 :     x = subrr(gel(L,q),gel(L,p));
    1634         408 :     if (signe(x))
    1635             :     {
    1636         408 :       x = divrr(x, shiftr(gcoeff(a,p,q),1));
    1637         408 :       y = sqrtr(addrr(unr, sqrr(x)));
    1638         408 :       t = invr((signe(x)>0)? addrr(x,y): subrr(x,y));
    1639             :     }
    1640             :     else
    1641           0 :       y = t = unr;
    1642         408 :     c = sqrtr(addrr(unr,sqrr(t)));
    1643         408 :     s = divrr(t,c);
    1644         408 :     u = divrr(t,addrr(unr,c));
    1645             : 
    1646             :     /* compute successive transforms of a and the matrix of accumulated
    1647             :      * rotations (r) */
    1648         408 :     for (i=1;   i<p; i++) rot(gcoeff(a,i,p), gcoeff(a,i,q), s,u);
    1649         408 :     for (i=p+1; i<q; i++) rot(gcoeff(a,p,i), gcoeff(a,i,q), s,u);
    1650         408 :     for (i=q+1; i<l; i++) rot(gcoeff(a,p,i), gcoeff(a,q,i), s,u);
    1651         408 :     y = gcoeff(a,p,q);
    1652         408 :     t = mulrr(t, y); shiftr_inplace(y, -de - 1);
    1653         408 :     x = gel(L,p); subrrz(x,t, x);
    1654         408 :     y = gel(L,q); addrrz(y,t, y);
    1655         408 :     for (i=1; i<l; i++) rot(gcoeff(r,i,p), gcoeff(r,i,q), s,u);
    1656             : 
    1657         408 :     e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1658        2856 :     for (j=1; j<l; j++)
    1659             :     {
    1660        8568 :       for (i=1; i<j; i++)
    1661             :       {
    1662        6120 :         GEN z = gcoeff(a,i,j);
    1663        6120 :         if (!signe(z)) continue;
    1664        6120 :         e = expo(z); if (e > e2) { e2=e; p=i; q=j; }
    1665             :       }
    1666        8568 :       for (i=j+1; i<l; i++)
    1667             :       {
    1668        6120 :         GEN z = gcoeff(a,j,i);
    1669        6120 :         if (!signe(z)) continue;
    1670        6120 :         e = expo(z); if (e > e2) { e2=e; p=j; q=i; }
    1671             :       }
    1672             :     }
    1673         408 :     avma = av2;
    1674             :   }
    1675             :   /* sort eigenvalues from smallest to largest */
    1676           6 :   c = indexsort(L);
    1677           6 :   r2 = vecpermute(r, c); for (i=1; i<l; i++) gel(r,i) = gel(r2,i);
    1678           6 :   L2 = vecpermute(L, c); for (i=1; i<l; i++) gel(L,i) = gel(L2,i);
    1679           6 :   avma = av1; return ja;
    1680             : }
    1681             : 
    1682             : /*************************************************************************/
    1683             : /**                                                                     **/
    1684             : /**                   Q-vector space -> Z-modules                       **/
    1685             : /**                                                                     **/
    1686             : /*************************************************************************/
    1687             : 
    1688             : GEN
    1689          48 : matrixqz0(GEN x,GEN p)
    1690             : {
    1691          48 :   if (typ(x) != t_MAT) pari_err_TYPE("QM_minors_coprime",x);
    1692          48 :   if (!p) return QM_minors_coprime(x,NULL);
    1693          30 :   if (typ(p) != t_INT) pari_err_TYPE("QM_minors_coprime",p);
    1694          30 :   if (signe(p)>=0) return QM_minors_coprime(x,p);
    1695          24 :   if (absequaliu(p,1)) return QM_ImZ_hnf(x); /* p = -1 */
    1696          12 :   if (absequaliu(p,2)) return QM_ImQ_hnf(x); /* p = -2 */
    1697           0 :   pari_err_FLAG("QM_minors_coprime");
    1698             :   return NULL; /* LCOV_EXCL_LINE */
    1699             : }
    1700             : 
    1701             : GEN
    1702          24 : QM_minors_coprime(GEN x, GEN D)
    1703             : {
    1704          24 :   pari_sp av = avma, av1;
    1705             :   long i, j, m, n, lP;
    1706             :   GEN P, y;
    1707             : 
    1708          24 :   n = lg(x)-1; if (!n) return gcopy(x);
    1709          24 :   m = nbrows(x);
    1710          24 :   if (n > m) pari_err_DOMAIN("QM_minors_coprime","n",">",strtoGENstr("m"),x);
    1711          18 :   y = x; x = cgetg(n+1,t_MAT);
    1712          48 :   for (j=1; j<=n; j++)
    1713             :   {
    1714          30 :     gel(x,j) = Q_primpart(gel(y,j));
    1715          30 :     RgV_check_ZV(gel(x,j), "QM_minors_coprime");
    1716             :   }
    1717             :   /* x now a ZM */
    1718          18 :   if (n==m)
    1719             :   {
    1720          12 :     if (gequal0(ZM_det(x)))
    1721           6 :       pari_err_DOMAIN("QM_minors_coprime", "rank(A)", "<",stoi(n),x);
    1722           6 :     avma = av; return matid(n);
    1723             :   }
    1724             :   /* m > n */
    1725           6 :   if (!D || gequal0(D))
    1726             :   {
    1727           6 :     pari_sp av2 = avma;
    1728           6 :     D = ZM_detmult(shallowtrans(x));
    1729           6 :     if (is_pm1(D)) { avma = av2; return ZM_copy(x); }
    1730             :   }
    1731           6 :   P = gel(Z_factor(D), 1); lP = lg(P);
    1732           6 :   av1 = avma;
    1733          12 :   for (i=1; i < lP; i++)
    1734             :   {
    1735           6 :     GEN p = gel(P,i), pov2 = shifti(p, -1);
    1736             :     for(;;)
    1737             :     {
    1738          18 :       GEN N, M = FpM_ker(x, p);
    1739          18 :       long lM = lg(M);
    1740          18 :       if (lM==1) break;
    1741             : 
    1742          12 :       FpM_center_inplace(M, p, pov2);
    1743          12 :       N = ZM_Z_divexact(ZM_mul(x,M), p);
    1744          24 :       for (j=1; j<lM; j++)
    1745             :       {
    1746          12 :         long k = n; while (!signe(gcoeff(M,k,j))) k--;
    1747          12 :         gel(x,k) = gel(N,j);
    1748             :       }
    1749          12 :       if (gc_needed(av1,1))
    1750             :       {
    1751           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"QM_minors_coprime, p = %Ps", p);
    1752           0 :         x = gerepilecopy(av1, x); pov2 = shifti(p, -1);
    1753             :       }
    1754          12 :     }
    1755             :   }
    1756           6 :   return gerepilecopy(av, x);
    1757             : }
    1758             : 
    1759             : static GEN
    1760         834 : QM_imZ_hnf_aux(GEN A)
    1761             : {
    1762             :   GEN D;
    1763         834 :   long l = lg(A);
    1764         834 :   if (l == 1) return cgetg(1,t_MAT);
    1765         834 :   if (l == 2) {
    1766             :     GEN c;
    1767           0 :     A = Q_primitive_part(A, &c);
    1768           0 :     if (!c) A = ZM_copy(A); else if ( isintzero(c) ) A = cgetg(1,t_MAT);
    1769           0 :     return A;
    1770             :   }
    1771         834 :   A = Q_remove_denom(A,&D);
    1772         834 :   if (D)
    1773             :   {
    1774             :     GEN B, U;
    1775          30 :     A = shallowconcat(A, scalarmat_shallow(D, nbrows(A)));
    1776          30 :     B = ZM_hnfall(A, &U, 1);
    1777          30 :     U = vecslice(U, 1, lg(A) - lg(B));
    1778          30 :     A = rowslice(U, l, nbrows(U));
    1779             :   }
    1780         834 :   return ZM_hnflll(A,NULL,1);
    1781             : }
    1782             : GEN
    1783          12 : QM_ImZ_hnf(GEN x)
    1784             : {
    1785          12 :   pari_sp av = avma;
    1786          12 :   return gerepileupto(av, QM_imZ_hnf_aux(x));
    1787             : }
    1788             : 
    1789             : GEN
    1790         822 : QM_ImQ_hnf(GEN x)
    1791             : {
    1792         822 :   pari_sp av = avma, av1;
    1793             :   long j,j1,k,m,n;
    1794             :   GEN c;
    1795             : 
    1796         822 :   n = lg(x); if (n==1) return gcopy(x);
    1797         822 :   m = lgcols(x); x = RgM_shallowcopy(x);
    1798         822 :   c = zero_zv(n-1);
    1799         822 :   av1 = avma;
    1800        9474 :   for (k=1; k<m; k++)
    1801             :   {
    1802        8652 :     j=1; while (j<n && (c[j] || gequal0(gcoeff(x,k,j)))) j++;
    1803        8652 :     if (j==n) continue;
    1804             : 
    1805        3990 :     c[j]=k; gel(x,j) = RgC_Rg_div(gel(x,j),gcoeff(x,k,j));
    1806       39714 :     for (j1=1; j1<n; j1++)
    1807       35724 :       if (j1!=j)
    1808             :       {
    1809       31734 :         GEN t = gcoeff(x,k,j1);
    1810       31734 :         if (!gequal0(t)) gel(x,j1) = RgC_sub(gel(x,j1), RgC_Rg_mul(gel(x,j),t));
    1811             :       }
    1812        3990 :     if (gc_needed(av1,1))
    1813             :     {
    1814           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QM_ImQ_hnf");
    1815           0 :       x = gerepilecopy(av1,x);
    1816             :     }
    1817             :   }
    1818         822 :   return gerepileupto(av, QM_imZ_hnf_aux(x));
    1819             : }
    1820             : 
    1821             : GEN
    1822           6 : intersect(GEN x, GEN y)
    1823             : {
    1824           6 :   long j, lx = lg(x);
    1825             :   pari_sp av;
    1826             :   GEN z;
    1827             : 
    1828           6 :   if (typ(x)!=t_MAT) pari_err_TYPE("intersect",x);
    1829           6 :   if (typ(y)!=t_MAT) pari_err_TYPE("intersect",y);
    1830           6 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
    1831             : 
    1832           6 :   av = avma; z = ker(shallowconcat(x,y));
    1833           6 :   for (j=lg(z)-1; j; j--) setlg(z[j], lx);
    1834           6 :   return gerepileupto(av, RgM_mul(x,z));
    1835             : }

Generated by: LCOV version 1.11