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-bordeaux1.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 18583-be9c204) Lines: 951 1083 87.8 %
Date: 2016-02-07 Functions: 70 76 92.1 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 674 908 74.2 %

           Branch data     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                 :       6027 : Flm_charpoly_i(GEN x, ulong p)
      30                 :            : {
      31                 :       6027 :   long lx = lg(x), r, i;
      32                 :       6027 :   GEN H, y = cgetg(lx+1, t_VEC);
      33                 :       6027 :   gel(y,1) = pol1_Flx(0); H = Flm_hess(x, p);
      34         [ +  + ]:      57995 :   for (r = 1; r < lx; r++)
      35                 :            :   {
      36                 :      51968 :     pari_sp av2 = avma;
      37                 :      51968 :     ulong a = 1;
      38                 :      51968 :     GEN z, b = zero_Flx(0);
      39         [ +  + ]:     141680 :     for (i = r-1; i; i--)
      40                 :            :     {
      41                 :     120064 :       a = Fl_mul(a, ucoeff(H,i+1,i), p);
      42         [ +  + ]:     120064 :       if (!a) break;
      43                 :      89712 :       b = Flx_add(b, Flx_Fl_mul(gel(y,i), Fl_mul(a,ucoeff(H,i,r),p), p), p);
      44                 :            :     }
      45                 :      51968 :     z = Flx_sub(Flx_shift(gel(y,r), 1),
      46                 :     103936 :                 Flx_Fl_mul(gel(y,r), ucoeff(H,r,r), p), p);
      47                 :            :     /* (X - H[r,r])y[r] - b */
      48                 :      51968 :     gel(y,r+1) = gerepileuptoleaf(av2, Flx_sub(z, b, p));
      49                 :            :   }
      50                 :       6027 :   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                 :       2198 : FpM_charpoly(GEN x, GEN p)
      62                 :            : {
      63                 :       2198 :   pari_sp av = avma;
      64                 :            :   long lx, r, i;
      65                 :            :   GEN y, H;
      66                 :            : 
      67         [ +  + ]:       2198 :   if (lgefint(p) == 3)
      68                 :            :   {
      69                 :       2184 :     ulong pp = p[2];
      70                 :       2184 :     y = Flx_to_ZX(Flm_charpoly_i(ZM_to_Flm(x,pp), pp));
      71                 :       2184 :     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                 :         42 :     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                 :       2198 :   return gerepileupto(av, gel(y,lx));
      93                 :            : }
      94                 :            : 
      95                 :            : GEN
      96                 :        280 : charpoly0(GEN x, long v, long flag)
      97                 :            : {
      98         [ +  + ]:        280 :   if (v<0) v = 0;
      99   [ +  +  +  +  :        280 :   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                 :        217 :       return charpoly(x, v);
     111                 :            :   }
     112                 :        273 :   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                 :      50036 : easychar(GEN x, long v)
     118                 :            : {
     119                 :            :   pari_sp av;
     120                 :            :   long lx;
     121                 :            :   GEN p1;
     122                 :            : 
     123   [ +  +  +  +  :      50036 :   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                 :      45367 :       GEN A = gel(x,2), T = gel(x,1);
     148 [ +  + ][ +  - ]:      45367 :       if (RgX_is_QX(A) && RgX_is_ZX(T))
     149                 :      45311 :         return QXQ_charpoly(A, T, v);
     150                 :            :       else
     151                 :         56 :         return RgXQ_charpoly(A, T, v);
     152                 :            :     }
     153                 :            :     case t_MAT:
     154                 :       4592 :       lx=lg(x);
     155         [ +  + ]:       4592 :       if (lx==1) return pol_1(v);
     156         [ +  + ]:       4536 :       if (lgcols(x) != lx) break;
     157                 :       4529 :       return NULL;
     158                 :            :   }
     159                 :          7 :   pari_err_TYPE("easychar",x);
     160                 :      50029 :   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                 :       1218 : charpoly(GEN x, long v)
     179                 :            : {
     180                 :       1218 :   GEN T, p = NULL;
     181         [ +  + ]:       1218 :   if ((T = easychar(x,v))) return T;
     182         [ +  + ]:       1106 :   if (RgM_is_ZM(x))
     183                 :            :   {
     184                 :        665 :     T = ZM_charpoly(x);
     185                 :        665 :     setvarn(T, v);
     186                 :            :   }
     187 [ +  + ][ +  - ]:        441 :   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         [ +  + ]:        406 :   else if (isinexact(x))
     194                 :          7 :     T = carhess(x, v);
     195                 :            :   else
     196                 :        399 :     T = carberkowitz(x, v);
     197                 :       1204 :   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                 :       3360 : fix_pol(pari_sp av, GEN p)
     204                 :            : {
     205                 :       3360 :   long w = gvar2(p), v = varn(p);
     206         [ +  + ]:       3360 :   if (w == v) pari_err_PRIORITY("charpoly", p, "=", w);
     207         [ +  + ]:       3353 :   if (varncmp(w,v) < 0) p = gerepileupto(av, poleval(p, pol_x(v)));
     208                 :       3353 :   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                 :         21 :     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         [ +  + ]:      14357 :   for (i = 2; i < lx; i++) t = gadd(t, gcoeff(x,i,i));
     264                 :       4795 :   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         [ +  + ]:         98 :   while ((p = u_forprime_next(&S)))
     274         [ +  + ]:         70 :     if (!umodiu(q, p)) return 1;
     275                 :       2324 :   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         [ +  + ]:       9247 :   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         [ +  + ]:       9282 :     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         [ +  + ]:       6937 :   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                 :       2387 :   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                 :        112 :   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         [ +  + ]:       3143 :   if (gcmp1(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                 :        476 :   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                 :        476 :   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         [ +  + ]:        833 :   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                 :        308 :   return NULL; /*NOT REACHED*/
     676                 :            : }
     677                 :            : 
     678                 :            : /*******************************************************************/
     679                 :            : /*                                                                 */
     680                 :            : /*                       MINIMAL POLYNOMIAL                        */
     681                 :            : /*                                                                 */
     682                 :            : /*******************************************************************/
     683                 :            : 
     684                 :            : static GEN
     685                 :      45458 : easymin(GEN x, long v)
     686                 :            : {
     687                 :      45458 :   pari_sp ltop=avma;
     688                 :            :   GEN G, R, dR;
     689 [ +  + ][ +  + ]:      45458 :   if (typ(x)==t_POLMOD && !issquarefree(gel(x,1)))
     690                 :         56 :     return NULL;
     691                 :      45402 :   R = easychar(x, v);
     692         [ +  + ]:      45402 :   if (!R) return R;
     693                 :      45353 :   dR=RgX_deriv(R);
     694         [ -  + ]:      45353 :   if (!lgpol(dR)) {avma=ltop; return NULL;}
     695                 :      45353 :   G=RgX_gcd(R,dR);
     696                 :      45353 :   G=RgX_Rg_div(G,leading_term(G));
     697                 :      45458 :   return gerepileupto(ltop, RgX_div(R,G));
     698                 :            : }
     699                 :            : 
     700                 :            : GEN
     701                 :      45486 : minpoly(GEN x, long v)
     702                 :            : {
     703                 :      45486 :   pari_sp ltop=avma;
     704                 :            :   GEN P;
     705         [ +  - ]:      45486 :   if (v<0) v = 0;
     706         [ +  + ]:      45486 :   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                 :      45458 :   P = easymin(x,v);
     713         [ +  + ]:      45458 :   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                 :      45486 :   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         [ +  + ]:        210 :     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                 :       6027 : Flm_hess(GEN x, ulong p)
     771                 :            : {
     772                 :       6027 :   long lx = lg(x), m, i, j;
     773         [ -  + ]:       6027 :   if (lx == 1) return cgetg(1,t_MAT);
     774         [ -  + ]:       6027 :   if (lgcols(x) != lx) pari_err_DIM("hess");
     775                 :            : 
     776                 :       6027 :   x = Flm_copy(x);
     777         [ +  + ]:      46186 :   for (m=2; m<lx-1; m++)
     778                 :            :   {
     779                 :      40159 :     ulong t = 0;
     780 [ +  + ][ +  + ]:     276381 :     for (i=m+1; i<lx; i++) { t = ucoeff(x,i,m-1); if (t) break; }
     781         [ +  + ]:      40159 :     if (i == lx) continue;
     782         [ +  + ]:     245112 :     for (j=m-1; j<lx; j++) lswap(ucoeff(x,i,j), ucoeff(x,m,j));
     783                 :      19649 :     swap(gel(x,i), gel(x,m)); t = Fl_inv(t, p);
     784                 :            : 
     785         [ +  + ]:     205814 :     for (i=m+1; i<lx; i++)
     786                 :            :     {
     787                 :     186165 :       ulong c = ucoeff(x,i,m-1);
     788         [ +  + ]:     186165 :       if (!c) continue;
     789                 :            : 
     790                 :      96362 :       c = Fl_mul(c,t,p); ucoeff(x,i,m-1) = 0;
     791         [ +  + ]:    1374408 :       for (j=m; j<lx; j++)
     792                 :    1278046 :         ucoeff(x,i,j) = Fl_sub(ucoeff(x,i,j), Fl_mul(c,ucoeff(x,m,j), p), p);
     793         [ +  + ]:    1980342 :       for (j=1; j<lx; j++)
     794                 :    1883980 :         ucoeff(x,j,m) = Fl_add(ucoeff(x,j,m), Fl_mul(c,ucoeff(x,j,i), p), p);
     795                 :            :     }
     796                 :            :   }
     797                 :       6027 :   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         [ +  + ]:         63 :     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                 :         42 :     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                 :       1085 : charpoly_bound(GEN M, GEN dM)
     873                 :            : {
     874                 :       1085 :   pari_sp av = avma;
     875                 :       1085 :   GEN B = itor(ZM_supnorm(M), LOWDEFAULTPREC);
     876                 :       1085 :   GEN s = real_0(LOWDEFAULTPREC), bin, B2;
     877                 :       1085 :   long n = lg(M)-1, k;
     878                 :            :   double d;
     879                 :       1085 :   bin = gen_1;
     880         [ +  + ]:       1085 :   if (dM) B = divri(B, dM);
     881                 :       1085 :   B2 = sqrr(B);
     882         [ +  + ]:       6307 :   for (k = n; k >= (n+1)>>1; k--)
     883                 :            :   {
     884                 :       5222 :     GEN t = mulri(powruhalf(mulur(k, B2), k), bin);
     885         [ +  + ]:       5222 :     if (absr_cmp(t, s) > 0) s = t;
     886                 :       5222 :     bin = diviuexact(muliu(bin, k), n-k+1);
     887                 :            :   }
     888                 :       1085 :   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                 :       3822 : QM_charpoly_Flx(GEN M, ulong dp, ulong p)
     894                 :            : {
     895                 :       3822 :   pari_sp av = avma;
     896                 :       3822 :   GEN H = Flm_charpoly_i(ZM_to_Flm(M,p), p);
     897         [ +  + ]:       3822 :   if (dp) H = Flx_rescale(H, Fl_inv(dp,p), p);
     898                 :       3822 :   return gerepileuptoleaf(av, H);
     899                 :            : }
     900                 :            : 
     901                 :            : static int
     902                 :       3822 : ZX_CRT(GEN *H, GEN Hp, GEN *q, ulong p, long bit)
     903                 :            : {
     904         [ +  + ]:       3822 :   if (!*H)
     905                 :            :   {
     906                 :       1190 :     *H = ZX_init_CRT(Hp, p, 0);
     907         [ -  + ]:       1190 :     if (DEBUGLEVEL>5)
     908                 :          0 :       err_printf("charpoly mod %lu, bound = 2^%ld\n", p, expu(p));
     909         [ +  + ]:       1190 :     if (expu(p) > bit) return 1;
     910                 :        616 :     *q = utoipos(p);
     911                 :            :   }
     912                 :            :   else
     913                 :            :   {
     914                 :       2632 :     int stable = ZX_incremental_CRT(H, Hp, q,p);
     915         [ -  + ]:       2632 :     if (DEBUGLEVEL>5)
     916                 :          0 :       err_printf("charpoly mod %lu (stable=%ld), bound = 2^%ld\n",
     917                 :            :                  p, stable, expi(*q));
     918 [ +  + ][ +  + ]:       2632 :     if (stable && expi(*q) > bit) return 1;
     919                 :            :   }
     920                 :       3822 :   return 0;
     921                 :            : }
     922                 :            : /* Let V = V1 \oplus V2 Q-vector spaces and f in End(V) stabilizing the Vi
     923                 :            :  * Let M, M2 be square QM representing f, f|V2.
     924                 :            :  * Return H := char(f|V1) = char(f) / char(f|V2) assuming H is in Z[X]
     925                 :            :  * and log 2 |H|oo <= bit */
     926                 :            : GEN
     927                 :          0 : QM_charpoly_ZX2_bound(GEN M, GEN M2, long bit)
     928                 :            : {
     929                 :          0 :   long n = lg(M)-1, n2 = lg(M2)-1;
     930                 :          0 :   GEN q = NULL, H1 = NULL, dM, dM2;
     931                 :            :   forprime_t S;
     932                 :            :   ulong p;
     933         [ #  # ]:          0 :   if (n == n2) return pol_1(0);
     934         [ #  # ]:          0 :   if (!n2) return QM_charpoly_ZX_bound(M, bit);
     935                 :          0 :   M = Q_remove_denom(M,&dM);
     936                 :          0 :   M2= Q_remove_denom(M2,&dM2);
     937                 :            : 
     938         [ #  # ]:          0 :   if (DEBUGLEVEL>5) err_printf("QM_charpoly_ZX2_bound: bit-bound 2^%ld\n", bit);
     939                 :          0 :   init_modular(&S);
     940         [ #  # ]:          0 :   while ((p = u_forprime_next(&S)))
     941                 :            :   {
     942                 :          0 :     ulong dMp = 0, dM2p = 0;
     943                 :            :     GEN Hp, H1p, H2p;
     944 [ #  # ][ #  # ]:          0 :     if (dM && !(dMp = umodiu(dM, p))) continue;
     945 [ #  # ][ #  # ]:          0 :     if (dM2 && !(dM2p = umodiu(dM2, p))) continue;
     946                 :          0 :     Hp  = QM_charpoly_Flx(M,  dMp,  p);
     947                 :          0 :     H2p = QM_charpoly_Flx(M2, dM2p, p);
     948                 :          0 :     H1p = Flx_div(Hp, H2p, p);
     949         [ #  # ]:          0 :     if (ZX_CRT(&H1, H1p, &q,p, bit)) break;
     950                 :            :   }
     951         [ #  # ]:          0 :   if (!p) pari_err_OVERFLOW("charpoly [ran out of primes]");
     952                 :          0 :   return H1;
     953                 :            : }
     954                 :            : 
     955                 :            : /* Assume M a square ZM, dM integer. Return charpoly(M / dM) in Z[X] */
     956                 :            : static GEN
     957                 :       1190 : QM_charpoly_ZX_i(GEN M, GEN dM, long bit)
     958                 :            : {
     959                 :       1190 :   long n = lg(M)-1;
     960                 :       1190 :   GEN q = NULL, H = NULL;
     961                 :            :   forprime_t S;
     962                 :            :   ulong p;
     963         [ -  + ]:       1190 :   if (!n) return pol_1(0);
     964                 :            : 
     965         [ +  + ]:       1190 :   if (bit < 0) bit = (long)charpoly_bound(M, dM) + 1;
     966         [ -  + ]:       1190 :   if (DEBUGLEVEL>5) err_printf("ZM_charpoly: bit-bound 2^%ld\n", bit);
     967                 :       1190 :   init_modular(&S);
     968         [ +  - ]:       3822 :   while ((p = u_forprime_next(&S)))
     969                 :            :   {
     970                 :       3822 :     ulong dMp = 0;
     971                 :            :     GEN Hp;
     972 [ +  + ][ -  + ]:       3822 :     if (dM && !(dMp = umodiu(dM, p))) continue;
     973                 :       3822 :     Hp = QM_charpoly_Flx(M, dMp, p);
     974         [ +  + ]:       3822 :     if (ZX_CRT(&H, Hp, &q,p, bit)) break;
     975                 :            :   }
     976         [ -  + ]:       1190 :   if (!p) pari_err_OVERFLOW("charpoly [ran out of primes]");
     977                 :       1190 :   return H;
     978                 :            : }
     979                 :            : GEN
     980                 :        105 : QM_charpoly_ZX_bound(GEN M, long bit)
     981                 :            : {
     982                 :        105 :   pari_sp av = avma;
     983                 :        105 :   GEN dM; M = Q_remove_denom(M, &dM);
     984                 :        105 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, bit));
     985                 :            : }
     986                 :            : GEN
     987                 :        413 : QM_charpoly_ZX(GEN M)
     988                 :            : {
     989                 :        413 :   pari_sp av = avma;
     990                 :        413 :   GEN dM; M = Q_remove_denom(M, &dM);
     991                 :        413 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, dM, -1));
     992                 :            : }
     993                 :            : GEN
     994                 :        672 : ZM_charpoly(GEN M)
     995                 :            : {
     996                 :        672 :   pari_sp av = avma;
     997                 :        672 :   return gerepilecopy(av, QM_charpoly_ZX_i(M, NULL, -1));
     998                 :            : }
     999                 :            : 
    1000                 :            : /*******************************************************************/
    1001                 :            : /*                                                                 */
    1002                 :            : /*        CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)        */
    1003                 :            : /*                                                                 */
    1004                 :            : /*******************************************************************/
    1005                 :            : GEN
    1006                 :        994 : carberkowitz(GEN x, long v)
    1007                 :            : {
    1008                 :            :   long lx, i, j, k, r;
    1009                 :            :   GEN V, S, C, Q;
    1010                 :            :   pari_sp av0, av;
    1011         [ -  + ]:        994 :   if ((V = easychar(x,v))) return V;
    1012                 :        994 :   lx = lg(x); av0 = avma;
    1013                 :        994 :   V = cgetg(lx+1, t_VEC);
    1014                 :        994 :   S = cgetg(lx+1, t_VEC);
    1015                 :        994 :   C = cgetg(lx+1, t_VEC);
    1016                 :        994 :   Q = cgetg(lx+1, t_VEC);
    1017                 :        994 :   av = avma;
    1018                 :        994 :   gel(C,1) = gen_m1;
    1019                 :        994 :   gel(V,1) = gen_m1;
    1020         [ +  + ]:       7056 :   for (i=2;i<=lx; i++) gel(C,i) = gel(Q,i) = gel(S,i) = gel(V,i) = gen_0;
    1021                 :        994 :   gel(V,2) = gcoeff(x,1,1);
    1022         [ +  + ]:       6062 :   for (r = 2; r < lx; r++)
    1023                 :            :   {
    1024                 :            :     pari_sp av2;
    1025                 :            :     GEN t;
    1026                 :            : 
    1027         [ +  + ]:      27377 :     for (i = 1; i < r; i++) gel(S,i) = gcoeff(x,i,r);
    1028                 :       5068 :     gel(C,2) = gcoeff(x,r,r);
    1029         [ +  + ]:      22309 :     for (i = 1; i < r-1; i++)
    1030                 :            :     {
    1031                 :      17241 :       av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1032         [ +  + ]:     115066 :       for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1033                 :      17241 :       gel(C,i+2) = gerepileupto(av2, t);
    1034         [ +  + ]:     132307 :       for (j = 1; j < r; j++)
    1035                 :            :       {
    1036                 :     115066 :         av2 = avma; t = gmul(gcoeff(x,j,1), gel(S,1));
    1037         [ +  + ]:     866222 :         for (k = 2; k < r; k++) t = gadd(t, gmul(gcoeff(x,j,k), gel(S,k)));
    1038                 :     115066 :         gel(Q,j) = gerepileupto(av2, t);
    1039                 :            :       }
    1040         [ +  + ]:     132307 :       for (j = 1; j < r; j++) gel(S,j) = gel(Q,j);
    1041                 :            :     }
    1042                 :       5068 :     av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    1043         [ +  + ]:      22309 :     for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    1044                 :       5068 :     gel(C,r+1) = gerepileupto(av2, t);
    1045         [ -  + ]:       5068 :     if (gc_needed(av0,1))
    1046                 :            :     {
    1047         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"carberkowitz");
    1048                 :          0 :       gerepileall(av, 2, &C, &V);
    1049                 :            :     }
    1050         [ +  + ]:      37513 :     for (i = 1; i <= r+1; i++)
    1051                 :            :     {
    1052                 :      32445 :       av2 = avma; t = gmul(gel(C,i), gel(V,1));
    1053         [ +  + ]:     134596 :       for (j = 2; j <= minss(r,i); j++)
    1054                 :     102151 :         t = gadd(t, gmul(gel(C,i+1-j), gel(V,j)));
    1055                 :      32445 :       gel(Q,i) = gerepileupto(av2, t);
    1056                 :            :     }
    1057         [ +  + ]:      37513 :     for (i = 1; i <= r+1; i++) gel(V,i) = gel(Q,i);
    1058                 :            :   }
    1059                 :        994 :   V = RgV_to_RgX(vecreverse(V), v); /* not gtopoly: fail if v > gvar(V) */
    1060         [ +  + ]:        994 :   V = odd(lx)? gcopy(V): RgX_neg(V);
    1061                 :        994 :   return fix_pol(av0, V);
    1062                 :            : }
    1063                 :            : 
    1064                 :            : /*******************************************************************/
    1065                 :            : /*                                                                 */
    1066                 :            : /*                            NORMS                                */
    1067                 :            : /*                                                                 */
    1068                 :            : /*******************************************************************/
    1069                 :            : GEN
    1070                 :     185569 : gnorm(GEN x)
    1071                 :            : {
    1072                 :            :   pari_sp av;
    1073                 :            :   long lx, i;
    1074                 :            :   GEN y;
    1075                 :            : 
    1076   [ +  +  +  +  :     185569 :   switch(typ(x))
          +  +  +  +  -  
                      - ]
    1077                 :            :   {
    1078                 :        716 :     case t_INT:  return sqri(x);
    1079                 :      18774 :     case t_REAL: return sqrr(x);
    1080                 :        105 :     case t_FRAC: return sqrfrac(x);
    1081                 :     161648 :     case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
    1082                 :        238 :     case t_QUAD:    av = avma; return gerepileupto(av, quadnorm(x));
    1083                 :            : 
    1084                 :         14 :     case t_POL: case t_SER: case t_RFRAC: av = avma;
    1085                 :         14 :       return gerepileupto(av, greal(gmul(gconj(x),x)));
    1086                 :            : 
    1087                 :            :     case t_FFELT:
    1088                 :         28 :       y = cgetg(3, t_INTMOD);
    1089                 :         28 :       gel(y,1) = FF_p(x);
    1090                 :         28 :       gel(y,2) = FF_norm(x); return y;
    1091                 :            : 
    1092                 :            :     case t_POLMOD:
    1093                 :            :     {
    1094                 :       4046 :       GEN T = gel(x,1), a = gel(x,2);
    1095 [ +  + ][ +  + ]:       4046 :       if (typ(a) != t_POL || varn(a) != varn(T)) return gpowgs(a, degpol(T));
    1096                 :       3913 :       return RgXQ_norm(a, T);
    1097                 :            :     }
    1098                 :            :     case t_VEC: case t_COL: case t_MAT:
    1099                 :          0 :       y = cgetg_copy(x, &lx);
    1100         [ #  # ]:          0 :       for (i=1; i<lx; i++) gel(y,i) = gnorm(gel(x,i));
    1101                 :          0 :       return y;
    1102                 :            :   }
    1103                 :          0 :   pari_err_TYPE("gnorm",x);
    1104                 :     185569 :   return NULL; /* not reached */
    1105                 :            : }
    1106                 :            : 
    1107                 :            : /* return |q|^2, complex modulus */
    1108                 :            : static GEN
    1109                 :         28 : cxquadnorm(GEN q, long prec)
    1110                 :            : {
    1111                 :         28 :   GEN X = gel(q,1), c = gel(X,2); /* (1-D)/4, -D/4 */
    1112         [ +  + ]:         28 :   if (signe(c) > 0) return quadnorm(q); /* imaginary */
    1113         [ +  + ]:         21 :   if (!prec) pari_err_TYPE("gnorml2", q);
    1114                 :         14 :   return sqrr(quadtofp(q, prec));
    1115                 :            : }
    1116                 :            : 
    1117                 :            : static GEN
    1118                 :    5585886 : gnorml2_i(GEN x, long prec)
    1119                 :            : {
    1120                 :            :   pari_sp av;
    1121                 :            :   long i, lx;
    1122                 :            :   GEN s;
    1123                 :            : 
    1124   [ +  +  +  +  :    5585886 :   switch(typ(x))
             +  +  +  - ]
    1125                 :            :   {
    1126                 :     448013 :     case t_INT:  return sqri(x);
    1127                 :    4268989 :     case t_REAL: return sqrr(x);
    1128                 :          7 :     case t_FRAC: return sqrfrac(x);
    1129                 :     391916 :     case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
    1130                 :         21 :     case t_QUAD:    av = avma; return gerepileupto(av, cxquadnorm(x,prec));
    1131                 :            : 
    1132                 :      10736 :     case t_POL: lx = lg(x)-1; x++; break;
    1133                 :            : 
    1134                 :            :     case t_VEC:
    1135                 :            :     case t_COL:
    1136                 :     466204 :     case t_MAT: lx = lg(x); break;
    1137                 :            : 
    1138                 :          0 :     default: pari_err_TYPE("gnorml2",x);
    1139                 :          0 :       return NULL; /* not reached */
    1140                 :            :   }
    1141         [ -  + ]:     476940 :   if (lx == 1) return gen_0;
    1142                 :     476940 :   av = avma;
    1143                 :     476940 :   s = gnorml2(gel(x,1));
    1144         [ +  + ]:    5076417 :   for (i=2; i<lx; i++)
    1145                 :            :   {
    1146                 :    4599477 :     s = gadd(s, gnorml2(gel(x,i)));
    1147         [ -  + ]:    4599477 :     if (gc_needed(av,1))
    1148                 :            :     {
    1149         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnorml2");
    1150                 :          0 :       s = gerepileupto(av, s);
    1151                 :            :     }
    1152                 :            :   }
    1153                 :    5585872 :   return gerepileupto(av,s);
    1154                 :            : }
    1155                 :            : GEN
    1156                 :    5585858 : gnorml2(GEN x) { return gnorml2_i(x, 0); }
    1157                 :            : 
    1158                 :            : static GEN pnormlp(GEN,GEN,long);
    1159                 :            : static GEN
    1160                 :         63 : pnormlpvec(long i0, GEN x, GEN p, long prec)
    1161                 :            : {
    1162                 :         63 :   pari_sp av = avma;
    1163                 :         63 :   long i, lx = lg(x);
    1164                 :         63 :   GEN s = gen_0;
    1165         [ +  + ]:        224 :   for (i=i0; i<lx; i++)
    1166                 :            :   {
    1167                 :        161 :     s = gadd(s, pnormlp(gel(x,i),p,prec));
    1168         [ -  + ]:        161 :     if (gc_needed(av,1))
    1169                 :            :     {
    1170         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gnormlp, i = %ld", i);
    1171                 :          0 :       s = gerepileupto(av, s);
    1172                 :            :     }
    1173                 :            :   }
    1174                 :         63 :   return s;
    1175                 :            : }
    1176                 :            : /* (||x||_p)^p */
    1177                 :            : static GEN
    1178                 :        196 : pnormlp(GEN x, GEN p, long prec)
    1179                 :            : {
    1180   [ +  -  +  +  :        196 :   switch(typ(x))
                   +  - ]
    1181                 :            :   {
    1182                 :        119 :     case t_INT: case t_REAL: x = mpabs(x); break;
    1183                 :          0 :     case t_FRAC: x = absfrac(x); break;
    1184                 :         14 :     case t_COMPLEX: case t_QUAD: x = gabs(x,prec); break;
    1185                 :          7 :     case t_POL: return pnormlpvec(2, x, p, prec);
    1186                 :         56 :     case t_VEC: case t_COL: case t_MAT: return pnormlpvec(1, x, p, prec);
    1187                 :          0 :     default: pari_err_TYPE("gnormlp",x);
    1188                 :            :   }
    1189                 :        196 :   return gpow(x, p, prec);
    1190                 :            : }
    1191                 :            : 
    1192                 :            : GEN
    1193                 :        315 : gnormlp(GEN x, GEN p, long prec)
    1194                 :            : {
    1195                 :        315 :   pari_sp av = avma;
    1196         [ +  + ]:        315 :   if (!p) return gsupnorm(x, prec);
    1197         [ -  + ]:        154 :   if (gsigne(p) <= 0) pari_err_DOMAIN("normlp", "p", "<=", gen_0, p);
    1198         [ +  + ]:        154 :   if (is_scalar_t(typ(x))) return gabs(x, prec);
    1199         [ +  + ]:         91 :   if (typ(p) == t_INT)
    1200                 :            :   {
    1201                 :         63 :     ulong pp = itou_or_0(p);
    1202      [ +  +  + ]:         63 :     switch(pp)
    1203                 :            :     {
    1204                 :         28 :       case 1: return gnorml1(x, prec);
    1205                 :         28 :       case 2: x = gnorml2_i(x, prec); break;
    1206                 :          7 :       default: x = pnormlp(x, p, prec); break;
    1207                 :            :     }
    1208 [ +  - ][ +  - ]:         35 :     if (pp && typ(x) == t_INT && Z_ispowerall(x, pp, &x))
                 [ +  + ]
    1209                 :          7 :       return gerepileuptoleaf(av, x);
    1210         [ +  - ]:         28 :     if (pp == 2) return gerepileupto(av, gsqrt(x, prec));
    1211                 :            :   }
    1212                 :            :   else
    1213                 :         28 :     x = pnormlp(x, p, prec);
    1214                 :         28 :   x = gpow(x, ginv(p), prec);
    1215                 :        315 :   return gerepileupto(av, x);
    1216                 :            : }
    1217                 :            : 
    1218                 :            : GEN
    1219                 :        168 : gnorml1(GEN x,long prec)
    1220                 :            : {
    1221                 :        168 :   pari_sp av = avma;
    1222                 :            :   long lx,i;
    1223                 :            :   GEN s;
    1224   [ +  -  +  +  :        168 :   switch(typ(x))
                   +  - ]
    1225                 :            :   {
    1226                 :         98 :     case t_INT: case t_REAL: return mpabs(x);
    1227                 :          0 :     case t_FRAC: return absfrac(x);
    1228                 :            : 
    1229                 :            :     case t_COMPLEX: case t_QUAD:
    1230                 :         14 :       return gabs(x,prec);
    1231                 :            : 
    1232                 :            :     case t_POL:
    1233                 :          7 :       lx = lg(x); s = gen_0;
    1234         [ +  + ]:         28 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1235                 :          7 :       break;
    1236                 :            : 
    1237                 :            :     case t_VEC: case t_COL: case t_MAT:
    1238                 :         49 :       lx = lg(x); s = gen_0;
    1239         [ +  + ]:        168 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1(gel(x,i),prec));
    1240                 :         49 :       break;
    1241                 :            : 
    1242                 :          0 :     default: pari_err_TYPE("gnorml1",x);
    1243                 :          0 :       return NULL; /* not reached */
    1244                 :            :   }
    1245                 :        168 :   return gerepileupto(av, s);
    1246                 :            : }
    1247                 :            : /* As gnorml1, except for t_QUAD and t_COMPLEX: |x + wy| := |x| + |y|
    1248                 :            :  * Still a norm of R-vector spaces, and can be cheaply computed without
    1249                 :            :  * square roots */
    1250                 :            : GEN
    1251                 :          0 : gnorml1_fake(GEN x)
    1252                 :            : {
    1253                 :          0 :   pari_sp av = avma;
    1254                 :            :   long lx, i;
    1255                 :            :   GEN s;
    1256   [ #  #  #  #  :          0 :   switch(typ(x))
                #  #  # ]
    1257                 :            :   {
    1258                 :          0 :     case t_INT: case t_REAL: return mpabs(x);
    1259                 :          0 :     case t_FRAC: return absfrac(x);
    1260                 :            : 
    1261                 :            :     case t_COMPLEX:
    1262                 :          0 :       s = gadd(gnorml1_fake(gel(x,1)), gnorml1_fake(gel(x,2)));
    1263                 :          0 :       break;
    1264                 :            :     case t_QUAD:
    1265                 :          0 :       s = gadd(gnorml1_fake(gel(x,2)), gnorml1_fake(gel(x,3)));
    1266                 :          0 :       break;
    1267                 :            : 
    1268                 :            :     case t_POL:
    1269                 :          0 :       lx = lg(x); s = gen_0;
    1270         [ #  # ]:          0 :       for (i=2; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1271                 :          0 :       break;
    1272                 :            : 
    1273                 :            :     case t_VEC: case t_COL: case t_MAT:
    1274                 :          0 :       lx = lg(x); s = gen_0;
    1275         [ #  # ]:          0 :       for (i=1; i<lx; i++) s = gadd(s, gnorml1_fake(gel(x,i)));
    1276                 :          0 :       break;
    1277                 :            : 
    1278                 :          0 :     default: pari_err_TYPE("gnorml1_fake",x);
    1279                 :          0 :       return NULL; /* not reached */
    1280                 :            :   }
    1281                 :          0 :   return gerepileupto(av, s);
    1282                 :            : }
    1283                 :            : 
    1284                 :            : static void
    1285 [ +  + ][ +  + ]:      45129 : store(GEN z, GEN *m) { if (!*m || gcmp(z, *m) > 0) *m = z; }
    1286                 :            : /* compare |x| to *m or |x|^2 to *msq, whichever is easiest, and update
    1287                 :            :  * the pointed value if x is larger */
    1288                 :            : void
    1289                 :      50099 : gsupnorm_aux(GEN x, GEN *m, GEN *msq, long prec)
    1290                 :            : {
    1291                 :            :   long i, lx;
    1292                 :            :   GEN z;
    1293   [ +  +  +  +  :      50099 :   switch(typ(x))
                +  +  - ]
    1294                 :            :   {
    1295                 :       4543 :     case t_COMPLEX: z = cxnorm(x); store(z, msq); return;
    1296                 :          7 :     case t_QUAD:  z = cxquadnorm(x,prec); store(z, msq); return;
    1297                 :      40516 :     case t_INT: case t_REAL: z = mpabs(x); store(z,m); return;
    1298                 :         63 :     case t_FRAC: z = absfrac(x); store(z,m); return;
    1299                 :            : 
    1300                 :       2919 :     case t_POL: lx = lg(x)-1; x++; break;
    1301                 :            : 
    1302                 :            :     case t_VEC:
    1303                 :            :     case t_COL:
    1304                 :       2051 :     case t_MAT: lx = lg(x); break;
    1305                 :            : 
    1306                 :          0 :     default: pari_err_TYPE("gsupnorm",x);
    1307                 :          0 :       return; /* not reached */
    1308                 :            :   }
    1309         [ +  + ]:      95242 :   for (i=1; i<lx; i++) gsupnorm_aux(gel(x,i), m, msq, prec);
    1310                 :            : }
    1311                 :            : GEN
    1312                 :       4956 : gsupnorm(GEN x, long prec)
    1313                 :            : {
    1314                 :       4956 :   GEN m = NULL, msq = NULL;
    1315                 :       4956 :   pari_sp av = avma;
    1316                 :       4956 :   gsupnorm_aux(x, &m, &msq, prec);
    1317                 :            :   /* now set m = max (m, sqrt(msq)) */
    1318         [ +  + ]:       4956 :   if (msq) {
    1319                 :        609 :     msq = gsqrt(msq, prec);
    1320 [ +  + ][ +  - ]:        609 :     if (!m || gcmp(m, msq) < 0) m = msq;
    1321         [ +  + ]:       4347 :   } else if (!m) m = gen_0;
    1322                 :       4956 :   return gerepilecopy(av, m);
    1323                 :            : }
    1324                 :            : 
    1325                 :            : /*******************************************************************/
    1326                 :            : /*                                                                 */
    1327                 :            : /*                            TRACES                               */
    1328                 :            : /*                                                                 */
    1329                 :            : /*******************************************************************/
    1330                 :            : GEN
    1331                 :         21 : matcompanion(GEN x)
    1332                 :            : {
    1333                 :         21 :   long n = degpol(x), j;
    1334                 :            :   GEN y, c;
    1335                 :            : 
    1336         [ -  + ]:         21 :   if (typ(x)!=t_POL) pari_err_TYPE("matcompanion",x);
    1337         [ +  + ]:         21 :   if (!signe(x)) pari_err_DOMAIN("matcompanion","polynomial","=",gen_0,x);
    1338         [ -  + ]:         14 :   if (n == 0) return cgetg(1, t_MAT);
    1339                 :            : 
    1340                 :         14 :   y = cgetg(n+1,t_MAT);
    1341         [ +  + ]:         49 :   for (j=1; j < n; j++) gel(y,j) = col_ei(n, j+1);
    1342                 :         14 :   c = cgetg(n+1,t_COL); gel(y,n) = c;
    1343         [ +  + ]:         14 :   if (gequal1(gel(x, n+2)))
    1344         [ +  + ]:         42 :     for (j=1; j<=n; j++) gel(c,j) = gneg(gel(x,j+1));
    1345                 :            :   else
    1346                 :            :   { /* not monic. Hardly ever used */
    1347                 :          7 :     pari_sp av = avma;
    1348                 :          7 :     GEN d = gclone(gneg(gel(x,n+2)));
    1349                 :          7 :     avma = av;
    1350         [ +  + ]:         21 :     for (j=1; j<=n; j++) gel(c,j) = gdiv(gel(x,j+1), d);
    1351                 :          7 :     gunclone(d);
    1352                 :            :   }
    1353                 :         14 :   return y;
    1354                 :            : }
    1355                 :            : 
    1356                 :            : GEN
    1357                 :     105803 : gtrace(GEN x)
    1358                 :            : {
    1359                 :            :   pari_sp av;
    1360                 :     105803 :   long i, lx, tx = typ(x);
    1361                 :            :   GEN y, z;
    1362                 :            : 
    1363   [ +  +  +  +  :     105803 :   switch(tx)
          +  +  +  +  -  
                   +  - ]
    1364                 :            :   {
    1365                 :            :     case t_INT: case t_REAL: case t_FRAC:
    1366                 :        730 :       return gmul2n(x,1);
    1367                 :            : 
    1368                 :            :     case t_COMPLEX:
    1369                 :     104443 :       return gmul2n(gel(x,1),1);
    1370                 :            : 
    1371                 :            :     case t_QUAD:
    1372                 :         28 :       y = gel(x,1);
    1373         [ +  - ]:         28 :       if (!gequal0(gel(y,3)))
    1374                 :            :       { /* assume quad. polynomial is either x^2 + d or x^2 - x + d */
    1375                 :         28 :         av = avma;
    1376                 :         28 :         return gerepileupto(av, gadd(gel(x,3), gmul2n(gel(x,2),1)));
    1377                 :            :       }
    1378                 :          0 :       return gmul2n(gel(x,2),1);
    1379                 :            : 
    1380                 :            :     case t_POL:
    1381                 :          7 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1382         [ +  + ]:         21 :       for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1383                 :          7 :       return normalizepol_lg(y, lx);
    1384                 :            : 
    1385                 :            :     case t_SER:
    1386         [ +  + ]:         14 :       if (ser_isexactzero(x)) return gcopy(x);
    1387                 :          7 :       y = cgetg_copy(x, &lx); y[1] = x[1];
    1388         [ +  + ]:         21 :       for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1389                 :          7 :       return normalize(y);
    1390                 :            : 
    1391                 :            :     case t_POLMOD:
    1392                 :        399 :       y = gel(x,1); z = gel(x,2);
    1393 [ +  + ][ -  + ]:        399 :       if (typ(z) != t_POL || varn(y) != varn(z)) return gmulsg(degpol(y), z);
    1394                 :        259 :       av = avma;
    1395                 :        259 :       return gerepileupto(av, quicktrace(z, polsym(y, degpol(y)-1)));
    1396                 :            : 
    1397                 :            :     case t_FFELT:
    1398                 :         28 :       y=cgetg(3, t_INTMOD);
    1399                 :         28 :       gel(y,1) = FF_p(x);
    1400                 :         28 :       gel(y,2) = FF_trace(x);
    1401                 :         28 :       return y;
    1402                 :            : 
    1403                 :            : 
    1404                 :            :     case t_RFRAC:
    1405                 :          7 :       return gadd(x, gconj(x));
    1406                 :            : 
    1407                 :            :     case t_VEC: case t_COL:
    1408                 :          0 :       y = cgetg_copy(x, &lx);
    1409         [ #  # ]:          0 :       for (i=1; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
    1410                 :          0 :       return y;
    1411                 :            : 
    1412                 :            :     case t_MAT:
    1413         [ +  + ]:        147 :       lx = lg(x); if (lx == 1) return gen_0;
    1414                 :            :       /*now lx >= 2*/
    1415         [ +  + ]:        140 :       if (lx != lgcols(x)) pari_err_DIM("gtrace");
    1416                 :        133 :       av = avma; return gerepileupto(av, mattrace(x));
    1417                 :            :   }
    1418                 :          0 :   pari_err_TYPE("gtrace",x);
    1419                 :     105796 :   return NULL; /* not reached */
    1420                 :            : }
    1421                 :            : 
    1422                 :            : /* Cholesky decomposition for positive definite matrix a
    1423                 :            :  * [GTM138, Algo 2.7.6, matrix Q]
    1424                 :            :  * If a is not positive definite return NULL. */
    1425                 :            : GEN
    1426                 :        867 : qfgaussred_positive(GEN a)
    1427                 :            : {
    1428                 :        867 :   pari_sp av = avma;
    1429                 :            :   GEN b;
    1430                 :        867 :   long i,j,k, n = lg(a);
    1431                 :            : 
    1432         [ -  + ]:        867 :   if (typ(a)!=t_MAT) pari_err_TYPE("qfgaussred_positive",a);
    1433         [ +  + ]:        867 :   if (n == 1) return cgetg(1, t_MAT);
    1434         [ -  + ]:        860 :   if (lgcols(a)!=n) pari_err_DIM("qfgaussred_positive");
    1435                 :        860 :   b = cgetg(n,t_MAT);
    1436         [ +  + ]:       5765 :   for (j=1; j<n; j++)
    1437                 :            :   {
    1438                 :       4905 :     GEN p1=cgetg(n,t_COL), p2=gel(a,j);
    1439                 :            : 
    1440                 :       4905 :     gel(b,j) = p1;
    1441         [ +  + ]:      30319 :     for (i=1; i<=j; i++) gel(p1,i) = gel(p2,i);
    1442         [ +  + ]:      25414 :     for (   ; i<n ; i++) gel(p1,i) = gen_0;
    1443                 :            :   }
    1444         [ +  + ]:       5765 :   for (k=1; k<n; k++)
    1445                 :            :   {
    1446                 :       4905 :     GEN bk, p = gcoeff(b,k,k), invp;
    1447         [ -  + ]:       4905 :     if (gsigne(p)<=0) { avma = av; return NULL; } /* not positive definite */
    1448                 :       4905 :     invp = ginv(p);
    1449                 :       4905 :     bk = row(b, k);
    1450         [ +  + ]:      25414 :     for (i=k+1; i<n; i++) gcoeff(b,k,i) = gmul(gel(bk,i), invp);
    1451         [ +  + ]:      25414 :     for (i=k+1; i<n; i++)
    1452                 :            :     {
    1453                 :      20509 :       GEN c = gel(bk, i);
    1454         [ +  + ]:     116933 :       for (j=i; j<n; j++)
    1455                 :      96424 :         gcoeff(b,i,j) = gsub(gcoeff(b,i,j), gmul(c,gcoeff(b,k,j)));
    1456                 :            :     }
    1457         [ -  + ]:       4905 :     if (gc_needed(av,1))
    1458                 :            :     {
    1459         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfgaussred_positive");
    1460                 :          0 :       b=gerepilecopy(av,b);
    1461                 :            :     }
    1462                 :            :   }
    1463                 :        867 :   return gerepilecopy(av,b);
    1464                 :            : }
    1465                 :            : 
    1466                 :            : /* Maximal pivot strategy: x is a suitable pivot if it is non zero and either
    1467                 :            :  * - an exact type, or
    1468                 :            :  * - it is maximal among remaining non-zero (t_REAL) pivots */
    1469                 :            : static int
    1470                 :      23534 : suitable(GEN x, long k, GEN *pp, long *pi)
    1471                 :            : {
    1472                 :      23534 :   long t = typ(x);
    1473   [ +  +  +  - ]:      23534 :   switch(t)
    1474                 :            :   {
    1475                 :       4522 :     case t_INT: return signe(x) != 0;
    1476                 :      18858 :     case t_FRAC: return 1;
    1477                 :            :     case t_REAL: {
    1478                 :        154 :       GEN p = *pp;
    1479 [ +  - ][ +  + ]:        154 :       if (signe(x) && (!p || absr_cmp(p, x) < 0)) { *pp = x; *pi = k; }
                 [ -  + ]
    1480                 :        154 :       return 0;
    1481                 :            :     }
    1482                 :      23534 :     default: return !gequal0(x);
    1483                 :            :   }
    1484                 :            : }
    1485                 :            : 
    1486                 :            : /* Gauss reduction (arbitrary symetric matrix, only the part above the
    1487                 :            :  * diagonal is considered). If signature is non-zero, return only the
    1488                 :            :  * signature, in which case gsigne() should be defined for elements of a. */
    1489                 :            : static GEN
    1490                 :       4263 : gaussred(GEN a, long signature)
    1491                 :            : {
    1492                 :            :   GEN r, ak, al;
    1493                 :            :   pari_sp av, av1;
    1494                 :       4263 :   long n = lg(a), i, j, k, l, sp, sn, t;
    1495                 :            : 
    1496         [ -  + ]:       4263 :   if (typ(a) != t_MAT) pari_err_TYPE("gaussred",a);
    1497 [ -  + ][ #  # ]:       4263 :   if (n == 1) return signature? mkvec2(gen_0, gen_0): cgetg(1, t_MAT);
    1498         [ -  + ]:       4263 :   if (lgcols(a) != n) pari_err_DIM("gaussred");
    1499                 :       4263 :   n--;
    1500                 :            : 
    1501                 :       4263 :   av = avma;
    1502                 :       4263 :   r = const_vecsmall(n, 1);
    1503                 :       4263 :   av1= avma;
    1504                 :       4263 :   a = RgM_shallowcopy(a);
    1505                 :       4263 :   t = n; sp = sn = 0;
    1506         [ +  + ]:      27657 :   while (t)
    1507                 :            :   {
    1508                 :      23394 :     long pind = 0;
    1509                 :      23394 :     GEN invp, p = NULL;
    1510 [ +  + ][ +  + ]:      82544 :     k=1; while (k<=n && (!r[k] || !suitable(gcoeff(a,k,k), k, &p, &pind))) k++;
                 [ +  + ]
    1511 [ +  + ][ +  + ]:      23394 :     if (k > n && p) k = pind;
    1512         [ +  + ]:      23394 :     if (k <= n)
    1513                 :            :     {
    1514                 :      23387 :       p = gcoeff(a,k,k); invp = ginv(p); /* != 0 */
    1515         [ +  + ]:      23387 :       if (signature) { /* skip if (!signature): gsigne may fail ! */
    1516         [ +  + ]:      23338 :         if (gsigne(p) > 0) sp++; else sn++;
    1517                 :            :       }
    1518                 :      23387 :       r[k] = 0; t--;
    1519                 :      23387 :       ak = row(a, k);
    1520         [ +  + ]:     164682 :       for (i=1; i<=n; i++)
    1521         [ +  + ]:     141295 :         gcoeff(a,k,i) = r[i]? gmul(gcoeff(a,k,i), invp): gen_0;
    1522                 :            : 
    1523 [ +  + ][ +  + ]:     164682 :       for (i=1; i<=n; i++) if (r[i])
    1524                 :            :       {
    1525                 :      58940 :         GEN c = gel(ak,i); /* - p * a[k,i] */
    1526         [ +  + ]:      58940 :         if (gequal0(c)) continue;
    1527 [ +  + ][ +  + ]:     445130 :         for (j=1; j<=n; j++) if (r[j])
    1528                 :     237965 :           gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
    1529                 :            :       }
    1530                 :      23387 :       gcoeff(a,k,k) = p;
    1531         [ -  + ]:      23387 :       if (gc_needed(av1,1))
    1532                 :            :       {
    1533         [ #  # ]:          0 :         if(DEBUGMEM>1) pari_warn(warnmem,"gaussred (t = %ld)", t);
    1534                 :          0 :         a = gerepilecopy(av1, a);
    1535                 :            :       }
    1536                 :            :     }
    1537                 :            :     else
    1538                 :            :     { /* all remaining diagonal coeffs are currently 0 */
    1539 [ +  - ][ +  - ]:          7 :       for (k=1; k<=n; k++) if (r[k])
    1540                 :            :       {
    1541 [ +  - ][ -  + ]:          7 :         l=k+1; while (l<=n && (!r[l] || !suitable(gcoeff(a,k,l), l, &p, &pind))) l++;
                 [ -  + ]
    1542 [ -  + ][ #  # ]:          7 :         if (l > n && p) l = pind;
    1543         [ -  + ]:          7 :         if (l > n) continue;
    1544                 :            : 
    1545                 :          7 :         p = gcoeff(a,k,l); invp = ginv(p);
    1546                 :          7 :         sp++; sn++;
    1547                 :          7 :         r[k] = r[l] = 0; t -= 2;
    1548                 :          7 :         ak = row(a, k);
    1549                 :          7 :         al = row(a, l);
    1550 [ +  + ][ +  + ]:         35 :         for (i=1; i<=n; i++) if (r[i])
    1551                 :            :         {
    1552                 :         14 :           gcoeff(a,k,i) = gmul(gcoeff(a,k,i), invp);
    1553                 :         14 :           gcoeff(a,l,i) = gmul(gcoeff(a,l,i), invp);
    1554                 :            :         } else {
    1555                 :         14 :           gcoeff(a,k,i) = gen_0;
    1556                 :         14 :           gcoeff(a,l,i) = gen_0;
    1557                 :            :         }
    1558                 :            : 
    1559 [ +  + ][ +  + ]:         35 :         for (i=1; i<=n; i++) if (r[i])
    1560                 :            :         { /* c = a[k,i] * p, d = a[l,i] * p; */
    1561                 :         14 :           GEN c = gel(ak,i), d = gel(al,i);
    1562 [ +  + ][ +  + ]:         70 :           for (j=1; j<=n; j++) if (r[j])
    1563                 :         28 :             gcoeff(a,i,j) = gsub(gcoeff(a,i,j),
    1564                 :         28 :                                  gadd(gmul(gcoeff(a,l,j), c),
    1565                 :         28 :                                       gmul(gcoeff(a,k,j), d)));
    1566                 :            :         }
    1567 [ +  + ][ +  + ]:         35 :         for (i=1; i<=n; i++) if (r[i])
    1568                 :            :         {
    1569                 :         14 :           GEN c = gcoeff(a,k,i), d = gcoeff(a,l,i);
    1570                 :         14 :           gcoeff(a,k,i) = gadd(c, d);
    1571                 :         14 :           gcoeff(a,l,i) = gsub(c, d);
    1572                 :            :         }
    1573                 :          7 :         gcoeff(a,k,l) = gen_1;
    1574                 :          7 :         gcoeff(a,l,k) = gen_m1;
    1575                 :          7 :         gcoeff(a,k,k) = gmul2n(p,-1);
    1576                 :          7 :         gcoeff(a,l,l) = gneg(gcoeff(a,k,k));
    1577         [ -  + ]:          7 :         if (gc_needed(av1,1))
    1578                 :            :         {
    1579         [ #  # ]:          0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gaussred");
    1580                 :          0 :           a = gerepilecopy(av1, a);
    1581                 :            :         }
    1582                 :          7 :         break;
    1583                 :            :       }
    1584         [ +  - ]:      23394 :       if (k > n) break;
    1585                 :            :     }
    1586                 :            :   }
    1587         [ +  + ]:       4263 :   if (!signature) return gerepilecopy(av, a);
    1588                 :       4263 :   avma = av; return mkvec2s(sp, sn);
    1589                 :            : }
    1590                 :            : 
    1591                 :            : GEN
    1592                 :         14 : qfgaussred(GEN a) { return gaussred(a,0); }
    1593                 :            : 
    1594                 :            : GEN
    1595                 :       4249 : qfsign(GEN a) { return gaussred(a,1); }
    1596                 :            : 
    1597                 :            : /* x -= s(y+u*x) */
    1598                 :            : /* y += s(x-u*y), simultaneously */
    1599                 :            : static void
    1600                 :       4760 : rot(GEN x, GEN y, GEN s, GEN u) {
    1601                 :       4760 :   GEN x1 = subrr(x, mulrr(s,addrr(y,mulrr(u,x))));
    1602                 :       4760 :   GEN y1 = addrr(y, mulrr(s,subrr(x,mulrr(u,y))));
    1603                 :       4760 :   affrr(x1,x);
    1604                 :       4760 :   affrr(y1,y);
    1605                 :       4760 : }
    1606                 :            : 
    1607                 :            : /* Diagonalization of a REAL symetric matrix. Return a vector [L, r]:
    1608                 :            :  * L = vector of eigenvalues
    1609                 :            :  * r = matrix of eigenvectors */
    1610                 :            : GEN
    1611                 :          7 : jacobi(GEN a, long prec)
    1612                 :            : {
    1613                 :            :   pari_sp av1;
    1614                 :          7 :   long de, e, e1, e2, i, j, p, q, l = lg(a);
    1615                 :            :   GEN c, ja, L, r, L2, r2, unr;
    1616                 :            : 
    1617         [ -  + ]:          7 :   if (typ(a) != t_MAT) pari_err_TYPE("jacobi",a);
    1618                 :          7 :   ja = cgetg(3,t_VEC);
    1619                 :          7 :   L = cgetg(l,t_COL); gel(ja,1) = L;
    1620                 :          7 :   r = cgetg(l,t_MAT); gel(ja,2) = r;
    1621         [ -  + ]:          7 :   if (l == 1) return ja;
    1622         [ -  + ]:          7 :   if (lgcols(a) != l) pari_err_DIM("jacobi");
    1623                 :            : 
    1624                 :          7 :   e1 = HIGHEXPOBIT-1;
    1625         [ +  + ]:         49 :   for (j=1; j<l; j++)
    1626                 :            :   {
    1627                 :         42 :     GEN z = gtofp(gcoeff(a,j,j), prec);
    1628                 :         42 :     gel(L,j) = z;
    1629         [ +  + ]:         42 :     e = expo(z); if (e < e1) e1 = e;
    1630                 :            :   }
    1631         [ +  + ]:         49 :   for (j=1; j<l; j++)
    1632                 :            :   {
    1633                 :         42 :     gel(r,j) = cgetg(l,t_COL);
    1634         [ +  + ]:        294 :     for (i=1; i<l; i++) gcoeff(r,i,j) = utor(i==j? 1: 0, prec);
    1635                 :            :   }
    1636                 :          7 :   av1 = avma;
    1637                 :            : 
    1638                 :          7 :   e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1639                 :          7 :   c = cgetg(l,t_MAT);
    1640         [ +  + ]:         49 :   for (j=1; j<l; j++)
    1641                 :            :   {
    1642                 :         42 :     gel(c,j) = cgetg(j,t_COL);
    1643         [ +  + ]:        147 :     for (i=1; i<j; i++)
    1644                 :            :     {
    1645                 :        105 :       GEN z = gtofp(gcoeff(a,i,j), prec);
    1646                 :        105 :       gcoeff(c,i,j) = z;
    1647         [ -  + ]:        105 :       if (!signe(z)) continue;
    1648         [ +  + ]:        105 :       e = expo(z); if (e > e2) { e2 = e; p = i; q = j; }
    1649                 :            :     }
    1650                 :            :   }
    1651                 :          7 :   a = c; unr = real_1(prec);
    1652                 :          7 :   de = prec2nbits(prec);
    1653                 :            : 
    1654                 :            :  /* e1 = min expo(a[i,i])
    1655                 :            :   * e2 = max expo(a[i,j]), i != j */
    1656         [ +  + ]:        483 :   while (e1-e2 < de)
    1657                 :            :   {
    1658                 :        476 :     pari_sp av2 = avma;
    1659                 :            :     GEN x, y, t, c, s, u;
    1660                 :            :     /* compute associated rotation in the plane formed by basis vectors number
    1661                 :            :      * p and q */
    1662                 :        476 :     x = subrr(gel(L,q),gel(L,p));
    1663         [ +  - ]:        476 :     if (signe(x))
    1664                 :            :     {
    1665                 :        476 :       x = divrr(x, shiftr(gcoeff(a,p,q),1));
    1666                 :        476 :       y = sqrtr(addrr(unr, sqrr(x)));
    1667         [ +  + ]:        476 :       t = invr((signe(x)>0)? addrr(x,y): subrr(x,y));
    1668                 :            :     }
    1669                 :            :     else
    1670                 :          0 :       y = t = unr;
    1671                 :        476 :     c = sqrtr(addrr(unr,sqrr(t)));
    1672                 :        476 :     s = divrr(t,c);
    1673                 :        476 :     u = divrr(t,addrr(unr,c));
    1674                 :            : 
    1675                 :            :     /* compute successive transforms of a and the matrix of accumulated
    1676                 :            :      * rotations (r) */
    1677         [ +  + ]:       1029 :     for (i=1;   i<p; i++) rot(gcoeff(a,i,p), gcoeff(a,i,q), s,u);
    1678         [ +  + ]:       1113 :     for (i=p+1; i<q; i++) rot(gcoeff(a,p,i), gcoeff(a,i,q), s,u);
    1679         [ +  + ]:       1190 :     for (i=q+1; i<l; i++) rot(gcoeff(a,p,i), gcoeff(a,q,i), s,u);
    1680                 :        476 :     y = gcoeff(a,p,q);
    1681                 :        476 :     t = mulrr(t, y); shiftr_inplace(y, -de - 1);
    1682                 :        476 :     x = gel(L,p); subrrz(x,t, x);
    1683                 :        476 :     y = gel(L,q); addrrz(y,t, y);
    1684         [ +  + ]:       3332 :     for (i=1; i<l; i++) rot(gcoeff(r,i,p), gcoeff(r,i,q), s,u);
    1685                 :            : 
    1686                 :        476 :     e2 = -(long)HIGHEXPOBIT; p = q = 1;
    1687         [ +  + ]:       3332 :     for (j=1; j<l; j++)
    1688                 :            :     {
    1689         [ +  + ]:       9996 :       for (i=1; i<j; i++)
    1690                 :            :       {
    1691                 :       7140 :         GEN z = gcoeff(a,i,j);
    1692         [ -  + ]:       7140 :         if (!signe(z)) continue;
    1693         [ -  + ]:       7140 :         e = expo(z); if (e > e2) { e2=e; p=i; q=j; }
    1694                 :            :       }
    1695         [ +  + ]:       9996 :       for (i=j+1; i<l; i++)
    1696                 :            :       {
    1697                 :       7140 :         GEN z = gcoeff(a,j,i);
    1698         [ -  + ]:       7140 :         if (!signe(z)) continue;
    1699         [ +  + ]:       7140 :         e = expo(z); if (e > e2) { e2=e; p=j; q=i; }
    1700                 :            :       }
    1701                 :            :     }
    1702                 :        476 :     avma = av2;
    1703                 :            :   }
    1704                 :            :   /* sort eigenvalues from smallest to largest */
    1705                 :          7 :   c = indexsort(L);
    1706         [ +  + ]:         49 :   r2 = vecpermute(r, c); for (i=1; i<l; i++) gel(r,i) = gel(r2,i);
    1707         [ +  + ]:         49 :   L2 = vecpermute(L, c); for (i=1; i<l; i++) gel(L,i) = gel(L2,i);
    1708                 :          7 :   avma = av1; return ja;
    1709                 :            : }
    1710                 :            : 
    1711                 :            : /*************************************************************************/
    1712                 :            : /**                                                                     **/
    1713                 :            : /**              MATRICE RATIONNELLE --> ENTIERE                        **/
    1714                 :            : /**                                                                     **/
    1715                 :            : /*************************************************************************/
    1716                 :            : 
    1717                 :            : GEN
    1718                 :         49 : matrixqz0(GEN x,GEN p)
    1719                 :            : {
    1720         [ -  + ]:         49 :   if (typ(x) != t_MAT) pari_err_TYPE("QM_minors_coprime",x);
    1721         [ +  + ]:         49 :   if (!p) return QM_minors_coprime(x,NULL);
    1722         [ -  + ]:         28 :   if (typ(p) != t_INT) pari_err_TYPE("QM_minors_coprime",p);
    1723         [ +  + ]:         28 :   if (signe(p)>=0) return QM_minors_coprime(x,p);
    1724         [ +  + ]:         21 :   if (equaliu(p,1)) return QM_ImZ_hnf(x); /* p = -1 */
    1725         [ +  - ]:          7 :   if (equaliu(p,2)) return QM_ImQ_hnf(x); /* p = -2 */
    1726                 :         35 :   pari_err_FLAG("QM_minors_coprime"); return NULL; /* not reached */
    1727                 :            : }
    1728                 :            : 
    1729                 :            : GEN
    1730                 :         28 : QM_minors_coprime(GEN x, GEN D)
    1731                 :            : {
    1732                 :         28 :   pari_sp av = avma, av1;
    1733                 :            :   long i, j, m, n, lP;
    1734                 :            :   GEN P, y;
    1735                 :            : 
    1736         [ -  + ]:         28 :   n = lg(x)-1; if (!n) return gcopy(x);
    1737                 :         28 :   m = nbrows(x);
    1738         [ +  + ]:         28 :   if (n > m) pari_err_DOMAIN("QM_minors_coprime","n",">",strtoGENstr("m"),x);
    1739                 :         21 :   y = x; x = cgetg(n+1,t_MAT);
    1740         [ +  + ]:         56 :   for (j=1; j<=n; j++)
    1741                 :            :   {
    1742                 :         35 :     gel(x,j) = Q_primpart(gel(y,j));
    1743                 :         35 :     RgV_check_ZV(gel(x,j), "QM_minors_coprime");
    1744                 :            :   }
    1745                 :            :   /* x now a ZM */
    1746         [ +  + ]:         21 :   if (n==m)
    1747                 :            :   {
    1748         [ +  + ]:         14 :     if (gequal0(ZM_det(x)))
    1749                 :          7 :       pari_err_DOMAIN("QM_minors_coprime", "rank(A)", "<",stoi(n),x);
    1750                 :          7 :     avma = av; return matid(n);
    1751                 :            :   }
    1752                 :            :   /* m > n */
    1753 [ +  - ][ +  - ]:          7 :   if (!D || gequal0(D))
    1754                 :            :   {
    1755                 :          7 :     pari_sp av2 = avma;
    1756                 :          7 :     D = ZM_detmult(shallowtrans(x));
    1757         [ -  + ]:          7 :     if (is_pm1(D)) { avma = av2; return ZM_copy(x); }
    1758                 :            :   }
    1759                 :          7 :   P = gel(Z_factor(D), 1); lP = lg(P);
    1760                 :          7 :   av1 = avma;
    1761         [ +  + ]:         14 :   for (i=1; i < lP; i++)
    1762                 :            :   {
    1763                 :          7 :     GEN p = gel(P,i), pov2 = shifti(p, -1);
    1764                 :            :     for(;;)
    1765                 :            :     {
    1766                 :         21 :       GEN N, M = FpM_ker(x, p);
    1767                 :         21 :       long lM = lg(M);
    1768         [ +  + ]:         21 :       if (lM==1) break;
    1769                 :            : 
    1770                 :         14 :       FpM_center_inplace(M, p, pov2);
    1771                 :         14 :       N = ZM_Z_divexact(ZM_mul(x,M), p);
    1772         [ +  + ]:         28 :       for (j=1; j<lM; j++)
    1773                 :            :       {
    1774         [ -  + ]:         14 :         long k = n; while (!signe(gcoeff(M,k,j))) k--;
    1775                 :         14 :         gel(x,k) = gel(N,j);
    1776                 :            :       }
    1777         [ -  + ]:         14 :       if (gc_needed(av1,1))
    1778                 :            :       {
    1779         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"QM_minors_coprime, p = %Ps", p);
    1780                 :          0 :         x = gerepilecopy(av1, x); pov2 = shifti(p, -1);
    1781                 :            :       }
    1782                 :         14 :     }
    1783                 :            :   }
    1784                 :         14 :   return gerepilecopy(av, x);
    1785                 :            : }
    1786                 :            : 
    1787                 :            : static GEN
    1788                 :        140 : RgC_Z_mul(GEN A, GEN u)
    1789                 :            : {
    1790                 :        140 :   long s = signe(u);
    1791 [ +  + ][ +  + ]:        140 :   if (is_pm1(u)) return s > 0? A: RgC_neg(A);
    1792         [ +  - ]:        140 :   return s? gmul(u,A): zerocol(lg(A)-1);
    1793                 :            : }
    1794                 :            : 
    1795                 :            : /* u,v integral, A,B RgC */
    1796                 :            : static GEN
    1797                 :         70 : RgC_lincomb(GEN u, GEN v, GEN A, GEN B)
    1798                 :            : {
    1799         [ -  + ]:         70 :   if (!signe(u)) return RgC_Z_mul(B,v);
    1800         [ -  + ]:         70 :   if (!signe(v)) return RgC_Z_mul(A,u);
    1801                 :         70 :   return RgC_add(RgC_Z_mul(A,u), RgC_Z_mul(B,v));
    1802                 :            : }
    1803                 :            : 
    1804                 :            : /* cf ZC_elem */
    1805                 :            : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
    1806                 :            :  * A[j] and A[k] of determinant 1. */
    1807                 :            : static void
    1808                 :         98 : QC_elem(GEN aj, GEN ak, GEN A, long j, long k)
    1809                 :            : {
    1810                 :            :   GEN p1, u, v, d;
    1811                 :            : 
    1812         [ +  + ]:         98 :   if (gequal0(ak)) { swap(gel(A,j), gel(A,k)); return; }
    1813         [ +  + ]:         49 :   if (typ(aj) == t_INT) {
    1814         [ +  + ]:         14 :     if (typ(ak) != t_INT) { aj = mulii(aj, gel(ak,2)); ak = gel(ak,1); }
    1815                 :            :   } else {
    1816         [ -  + ]:         35 :     if (typ(ak) == t_INT) { ak = mulii(ak, gel(aj,2)); aj = gel(aj,1); }
    1817                 :            :     else {
    1818                 :         35 :       GEN daj = gel(aj,2), dak = gel(ak,2), D = gcdii(daj, dak);
    1819                 :         35 :       aj = gel(aj,1); ak = gel(ak,1);
    1820         [ +  + ]:         35 :       if (!is_pm1(D)) { daj = diviiexact(daj, D); dak = diviiexact(dak, D); }
    1821         [ +  + ]:         35 :       if (!is_pm1(dak)) aj = mulii(aj, dak);
    1822         [ +  + ]:         35 :       if (!is_pm1(daj)) ak = mulii(ak, daj);
    1823                 :            :     }
    1824                 :            :   }
    1825                 :            :   /* aj,ak were multiplied by their least common denominator */
    1826                 :            : 
    1827                 :         49 :   d = bezout(aj,ak,&u,&v);
    1828                 :            :   /* frequent special case (u,v) = (1,0) or (0,1) */
    1829         [ +  + ]:         49 :   if (!signe(u))
    1830                 :            :   { /* ak | aj */
    1831                 :         14 :     GEN c = negi(diviiexact(aj,ak));
    1832                 :         14 :     gel(A,j) = RgC_lincomb(gen_1, c, gel(A,j), gel(A,k));
    1833                 :            :     return;
    1834                 :            :   }
    1835         [ +  + ]:         35 :   if (!signe(v))
    1836                 :            :   { /* aj | ak */
    1837                 :         14 :     GEN c = negi(diviiexact(ak,aj));
    1838                 :         14 :     gel(A,k) = RgC_lincomb(gen_1, c, gel(A,k), gel(A,j));
    1839                 :         14 :     swap(gel(A,j), gel(A,k));
    1840                 :            :     return;
    1841                 :            :   }
    1842                 :            : 
    1843         [ -  + ]:         21 :   if (!is_pm1(d)) { aj = diviiexact(aj,d); ak = diviiexact(ak,d); }
    1844                 :         21 :   p1 = gel(A,k);
    1845                 :         21 :   gel(A,k) = RgC_lincomb(u,v, gel(A,j),p1);
    1846                 :         98 :   gel(A,j) = RgC_lincomb(negi(aj),ak, p1,gel(A,j));
    1847                 :            : }
    1848                 :            : 
    1849                 :            : static GEN
    1850                 :         21 : QM_imZ_hnf_aux(GEN A)
    1851                 :            : {
    1852                 :         21 :   pari_sp av = avma;
    1853                 :            :   long i,j,k,n,m;
    1854                 :            : 
    1855                 :         21 :   n = lg(A);
    1856         [ -  + ]:         21 :   if (n == 1) return cgetg(1,t_MAT);
    1857         [ -  + ]:         21 :   if (n == 2) {
    1858                 :            :     GEN c;
    1859                 :          0 :     A = Q_primitive_part(A, &c);
    1860 [ #  # ][ #  # ]:          0 :     if (!c) A = ZM_copy(A); else if ( isintzero(c) ) A = cgetg(1,t_MAT);
    1861                 :          0 :     return A;
    1862                 :            :   }
    1863                 :         21 :   m = lgcols(A);
    1864         [ +  + ]:         77 :   for (i=1; i<m; i++)
    1865                 :            :   {
    1866                 :            :     GEN b;
    1867         [ +  + ]:        189 :     for (j = k = 1; j<n; j++)
    1868                 :            :     {
    1869                 :        133 :       GEN a = gcoeff(A,i,j);
    1870         [ +  + ]:        133 :       if (gequal0(a)) continue;
    1871                 :            : 
    1872         [ +  + ]:         98 :       k = j+1; if (k == n) k = 1;
    1873                 :            :       /* zero a = Aij  using  b = Aik */
    1874                 :         98 :       QC_elem(a, gcoeff(A,i,k), A, j,k);
    1875                 :            :     }
    1876                 :         56 :     b = gcoeff(A,i,k);
    1877         [ +  + ]:         56 :     if (!gequal0(b))
    1878                 :            :     {
    1879                 :         42 :       b = Q_denom(b);
    1880         [ +  + ]:         42 :       if (!is_pm1(b)) gel(A,k) = RgC_Rg_mul(gel(A,k), b);
    1881                 :            :     }
    1882         [ -  + ]:         56 :     if (gc_needed(av,1))
    1883                 :            :     {
    1884         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QM_imZ_hnf_aux");
    1885                 :          0 :       A = gerepilecopy(av,A);
    1886                 :            :     }
    1887                 :            :   }
    1888                 :         21 :   return ZM_hnf(A);
    1889                 :            : }
    1890                 :            : 
    1891                 :            : GEN
    1892                 :         14 : QM_ImZ_hnf(GEN x)
    1893                 :            : {
    1894                 :         14 :   pari_sp av = avma;
    1895                 :         14 :   return gerepileupto(av, QM_imZ_hnf_aux( RgM_shallowcopy(x) ));
    1896                 :            : }
    1897                 :            : 
    1898                 :            : GEN
    1899                 :          7 : QM_ImQ_hnf(GEN x)
    1900                 :            : {
    1901                 :          7 :   pari_sp av = avma, av1;
    1902                 :            :   long j,j1,k,m,n;
    1903                 :            :   GEN c;
    1904                 :            : 
    1905         [ -  + ]:          7 :   n = lg(x); if (n==1) return gcopy(x);
    1906                 :          7 :   m = lgcols(x); x = RgM_shallowcopy(x);
    1907                 :          7 :   c = zero_zv(n-1);
    1908                 :          7 :   av1 = avma;
    1909         [ +  + ]:         28 :   for (k=1; k<m; k++)
    1910                 :            :   {
    1911 [ +  + ][ +  + ]:         42 :     j=1; while (j<n && (c[j] || gequal0(gcoeff(x,k,j)))) j++;
                 [ -  + ]
    1912         [ +  + ]:         21 :     if (j==n) continue;
    1913                 :            : 
    1914                 :         14 :     c[j]=k; gel(x,j) = RgC_Rg_div(gel(x,j),gcoeff(x,k,j));
    1915         [ +  + ]:         42 :     for (j1=1; j1<n; j1++)
    1916         [ +  + ]:         28 :       if (j1!=j)
    1917                 :            :       {
    1918                 :         14 :         GEN t = gcoeff(x,k,j1);
    1919         [ +  - ]:         14 :         if (!gequal0(t)) gel(x,j1) = RgC_sub(gel(x,j1), RgC_Rg_mul(gel(x,j),t));
    1920                 :            :       }
    1921         [ -  + ]:         14 :     if (gc_needed(av1,1))
    1922                 :            :     {
    1923         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QM_ImQ_hnf");
    1924                 :          0 :       x = gerepilecopy(av1,x);
    1925                 :            :     }
    1926                 :            :   }
    1927                 :          7 :   return gerepileupto(av, QM_imZ_hnf_aux(x));
    1928                 :            : }
    1929                 :            : 
    1930                 :            : GEN
    1931                 :          7 : intersect(GEN x, GEN y)
    1932                 :            : {
    1933                 :          7 :   long j, lx = lg(x);
    1934                 :            :   pari_sp av;
    1935                 :            :   GEN z;
    1936                 :            : 
    1937         [ -  + ]:          7 :   if (typ(x)!=t_MAT) pari_err_TYPE("intersect",x);
    1938         [ -  + ]:          7 :   if (typ(y)!=t_MAT) pari_err_TYPE("intersect",y);
    1939 [ +  - ][ -  + ]:          7 :   if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
    1940                 :            : 
    1941                 :          7 :   av = avma; z = ker(shallowconcat(x,y));
    1942         [ +  + ]:         14 :   for (j=lg(z)-1; j; j--) setlg(z[j], lx);
    1943                 :          7 :   return gerepileupto(av, RgM_mul(x,z));
    1944                 :            : }

Generated by: LCOV version 1.9