Code coverage tests

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

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

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

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

Generated by: LCOV version 1.11