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

Generated by: LCOV version 1.13