Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 18948-dd2a64a) Lines: 1556 1664 93.5 %
Date: 2016-05-28 Functions: 117 132 88.6 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 989 1221 81.0 %

           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                 :            : /*                        NUMBER FIELDS                       */
      17                 :            : /*                                                            */
      18                 :            : /**************************************************************/
      19                 :            : #include "pari.h"
      20                 :            : #include "paripriv.h"
      21                 :            : 
      22                 :            : int new_galois_format = 0;
      23                 :            : 
      24                 :            : int
      25                 :     211498 : checkrnf_i(GEN rnf)
      26 [ +  + ][ +  + ]:     211498 : { return (typ(rnf)==t_VEC && lg(rnf)==13); }
      27                 :            : 
      28                 :            : void
      29                 :     208425 : checkrnf(GEN rnf)
      30         [ +  + ]:     208425 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
      31                 :            : 
      32                 :            : GEN
      33                 :     686681 : checkbnf_i(GEN X)
      34                 :            : {
      35         [ +  + ]:     686681 :   if (typ(X) == t_VEC)
      36      [ +  +  + ]:     686240 :     switch (lg(X))
      37                 :            :     {
      38                 :     684693 :       case 11: return X;
      39                 :       1008 :       case 7:  return checkbnf_i(bnr_get_bnf(X));
      40                 :            :     }
      41                 :     686681 :   return NULL;
      42                 :            : }
      43                 :            : 
      44                 :            : GEN
      45                 :    5834995 : checknf_i(GEN X)
      46                 :            : {
      47         [ +  + ]:    5834995 :   if (typ(X)==t_VEC)
      48   [ +  +  +  +  :    5834547 :     switch(lg(X))
                      + ]
      49                 :            :     {
      50                 :    5614568 :       case 10: return X;
      51                 :     213917 :       case 11: return checknf_i(bnf_get_nf(X));
      52                 :       5593 :       case 7:  return checknf_i(bnr_get_bnf(X));
      53         [ +  + ]:        175 :       case 3: if (typ(gel(X,2)) == t_POLMOD) return checknf_i(gel(X,1));
      54                 :            :     }
      55                 :    5834995 :   return NULL;
      56                 :            : }
      57                 :            : 
      58                 :            : GEN
      59                 :     684546 : checkbnf(GEN x)
      60                 :            : {
      61                 :     684546 :   GEN bnf = checkbnf_i(x);
      62         [ -  + ]:     684546 :   if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
      63                 :     684546 :   return bnf;
      64                 :            : }
      65                 :            : 
      66                 :            : GEN
      67                 :    5613917 : checknf(GEN x)
      68                 :            : {
      69                 :    5613917 :   GEN nf = checknf_i(x);
      70         [ +  + ]:    5613917 :   if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
      71                 :    5613903 :   return nf;
      72                 :            : }
      73                 :            : 
      74                 :            : void
      75                 :     197230 : checkbnr(GEN bnr)
      76                 :            : {
      77 [ +  - ][ -  + ]:     197230 :   if (typ(bnr)!=t_VEC || lg(bnr)!=7)
      78                 :          0 :     pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
      79                 :     197230 :   (void)checkbnf(bnr_get_bnf(bnr));
      80                 :     197230 : }
      81                 :            : 
      82                 :            : void
      83                 :          0 : checkbnrgen(GEN bnr)
      84                 :            : {
      85                 :          0 :   checkbnr(bnr);
      86         [ #  # ]:          0 :   if (lg(bnr_get_clgp(bnr))<=3)
      87                 :          0 :     pari_err_TYPE("checkbnrgen [apply bnrinit(,,1), not bnrinit()]",bnr);
      88                 :          0 : }
      89                 :            : 
      90                 :            : void
      91                 :          0 : checksqmat(GEN x, long N)
      92                 :            : {
      93         [ #  # ]:          0 :   if (typ(x)!=t_MAT) pari_err_TYPE("checksqmat",x);
      94 [ #  # ][ #  # ]:          0 :   if (lg(x) == 1 || lgcols(x) != N+1) pari_err_DIM("checksqmat");
      95                 :          0 : }
      96                 :            : 
      97                 :            : GEN
      98                 :     205987 : checkbid_i(GEN bid)
      99                 :            : {
     100                 :            :   GEN f;
     101 [ +  + ][ +  + ]:     205987 :   if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(gel(bid,3)) != t_MAT) return NULL;
                 [ +  + ]
     102                 :     199785 :   f = bid_get_mod(bid);
     103 [ +  - ][ -  + ]:     199785 :   if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
     104                 :     205987 :   return bid;
     105                 :            : }
     106                 :            : void
     107                 :     199785 : checkbid(GEN bid)
     108                 :            : {
     109         [ +  + ]:     199785 :   if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
     110                 :     199778 : }
     111                 :            : void
     112                 :      11403 : checkabgrp(GEN v)
     113                 :            : {
     114         [ +  - ]:      11403 :   if (typ(v) == t_VEC) switch(lg(v))
              [ +  +  - ]
     115                 :            :   {
     116         [ -  + ]:      11312 :     case 4: if (typ(gel(v,3)) != t_VEC) break;
     117         [ +  + ]:      11403 :     case 3: if (typ(gel(v,2)) != t_VEC) break;
     118         [ -  + ]:      11375 :             if (typ(gel(v,1)) != t_INT) break;
     119                 :      11375 :             return;/*OK*/
     120                 :          0 :     default: break;
     121                 :            :   }
     122                 :         28 :   pari_err_TYPE("checkabgrp",v);
     123                 :            : }
     124                 :            : 
     125                 :            : GEN
     126                 :      84217 : checknfelt_mod(GEN nf, GEN x, const char *s)
     127                 :            : {
     128                 :      84217 :   GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
     129         [ +  + ]:      84217 :   if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
     130                 :      84147 :   return a;
     131                 :            : }
     132                 :            : 
     133                 :            : void
     134                 :       3206 : check_ZKmodule(GEN x, const char *s)
     135                 :            : {
     136 [ +  - ][ -  + ]:       3206 :   if (typ(x) != t_VEC || lg(x) < 3) pari_err_TYPE(s,x);
     137         [ -  + ]:       3206 :   if (typ(gel(x,1)) != t_MAT) pari_err_TYPE(s,gel(x,1));
     138         [ -  + ]:       3206 :   if (typ(gel(x,2)) != t_VEC) pari_err_TYPE(s,gel(x,2));
     139         [ -  + ]:       3206 :   if (lg(gel(x,2)) != lgcols(x)) pari_err_DIM(s);
     140                 :       3206 : }
     141                 :            : 
     142                 :            : static long
     143                 :      51296 : typv6(GEN x)
     144                 :            : {
     145 [ +  + ][ +  - ]:      51296 :   if (typ(gel(x,1)) == t_VEC && lg(gel(x,3)) == 3)
     146                 :            :   {
     147                 :       1806 :     long t = typ(gel(x,3));
     148 [ +  + ][ +  - ]:       1806 :     return (t == t_MAT || t == t_VEC)? typ_BID: typ_NULL;
     149                 :            :   }
     150 [ +  + ][ +  - ]:      49490 :   if (typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT) return typ_PRID;
     151                 :      51296 :   return typ_NULL;
     152                 :            : }
     153                 :            : 
     154                 :            : GEN
     155                 :      13041 : get_bnf(GEN x, long *t)
     156                 :            : {
     157   [ +  +  +  +  :      13041 :   switch(typ(x))
                      + ]
     158                 :            :   {
     159                 :         56 :     case t_POL: *t = typ_POL;  return NULL;
     160                 :         56 :     case t_QUAD: *t = typ_Q  ; return NULL;
     161                 :            :     case t_VEC:
     162   [ +  +  +  +  :      12425 :       switch(lg(x))
             +  +  +  +  
                      + ]
     163                 :            :       {
     164                 :       4382 :         case 5: *t = typ_QUA; return NULL;
     165                 :        322 :         case 6: *t = typv6(x); return NULL;
     166                 :         77 :         case 7:  *t = typ_BNR;
     167 [ +  - ][ +  - ]:         77 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     168                 :         77 :           return x;
     169                 :            :         case 9:
     170                 :         56 :           x = gel(x,2);
     171 [ +  - ][ +  - ]:         56 :           if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
     172                 :         56 :           return NULL;
     173                 :        217 :         case 10: *t = typ_NF; return NULL;
     174                 :        189 :         case 11: *t = typ_BNF; return x;
     175                 :         56 :         case 13: *t = typ_RNF; return NULL;
     176                 :        266 :         case 17: *t = typ_ELL; return NULL;
     177                 :            :       }
     178                 :       6860 :       break;
     179                 :            :     case t_COL:
     180         [ +  + ]:        112 :       if (get_prid(x)) { *t = typ_MODPR; return NULL; }
     181                 :         56 :       break;
     182                 :            :   }
     183                 :      13041 :   *t = typ_NULL; return NULL;
     184                 :            : }
     185                 :            : 
     186                 :            : GEN
     187                 :      62223 : get_nf(GEN x, long *t)
     188                 :            : {
     189   [ +  +  +  +  :      62223 :   switch(typ(x))
                      + ]
     190                 :            :   {
     191                 :        133 :     case t_POL : *t = typ_POL; return NULL;
     192                 :        133 :     case t_QUAD: *t = typ_Q  ; return NULL;
     193                 :            :     case t_VEC:
     194   [ +  +  +  +  :      59808 :       switch(lg(x))
          +  +  +  +  +  
                      + ]
     195                 :            :       {
     196                 :            :         case 3:
     197         [ -  + ]:        133 :           if (typ(gel(x,2)) != t_POLMOD) break;
     198                 :        133 :           return get_nf(gel(x,1),t);
     199                 :        133 :         case 5: *t = typ_QUA; return NULL;
     200                 :      50302 :         case 6: *t = typv6(x); return NULL;
     201                 :        140 :         case 7: *t = typ_BNR;
     202 [ +  - ][ +  - ]:        140 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     203 [ +  - ][ +  - ]:        140 :           x = bnf_get_nf(x);  if (typ(x)!=t_VEC || lg(x)!=10) break;
     204                 :        140 :           return x;
     205                 :            :         case 9:
     206                 :        238 :           x = gel(x,2);
     207 [ +  - ][ +  - ]:        238 :           if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
     208                 :        238 :           return NULL;
     209                 :        672 :         case 10: *t = typ_NF; return x;
     210                 :       4788 :         case 11: *t = typ_BNF;
     211 [ +  - ][ +  - ]:       4788 :           x = bnf_get_nf(x); if (typ(x)!=t_VEC || lg(x)!=10) break;
     212                 :       4788 :           return x;
     213                 :        336 :         case 13: *t = typ_RNF; return NULL;
     214                 :       2933 :         case 17: *t = typ_ELL; return NULL;
     215                 :            :       }
     216                 :        133 :       break;
     217                 :            :     case t_COL:
     218         [ +  + ]:        266 :       if (get_prid(x)) { *t = typ_MODPR; return NULL; }
     219                 :        133 :       break;
     220                 :            :   }
     221                 :      62223 :   *t = typ_NULL; return NULL;
     222                 :            : }
     223                 :            : 
     224                 :            : long
     225                 :      38255 : nftyp(GEN x)
     226                 :            : {
     227   [ +  +  +  + ]:      38255 :   switch(typ(x))
     228                 :            :   {
     229                 :         14 :     case t_POL : return typ_POL;
     230                 :          7 :     case t_QUAD: return typ_Q;
     231                 :            :     case t_VEC:
     232   [ +  +  +  +  :      38227 :       switch(lg(x))
             +  +  +  + ]
     233                 :            :       {
     234                 :        161 :         case 13: return typ_RNF;
     235                 :            :         case 10:
     236         [ +  + ]:      36435 :           if (typ(gel(x,1))!=t_POL) break;
     237                 :      36428 :           return typ_NF;
     238                 :            :         case 11:
     239 [ +  - ][ +  - ]:         28 :           x = bnf_get_nf(x); if (typ(x)!=t_VEC || lg(x)!=10) break;
     240                 :         28 :           return typ_BNF;
     241                 :            :         case 7:
     242 [ +  + ][ +  - ]:        896 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     243 [ +  - ][ +  - ]:        889 :           x = bnf_get_nf(x);  if (typ(x)!=t_VEC || lg(x)!=10) break;
     244                 :        889 :           return typ_BNR;
     245                 :            :         case 6:
     246                 :        672 :           return typv6(x);
     247                 :            :         case 9:
     248                 :          7 :           x = gel(x,2);
     249 [ -  + ][ #  # ]:          7 :           if (typ(x) == t_VEC && lg(x) == 4) return typ_GAL;
     250                 :         14 :         case 17: return typ_ELL;
     251                 :            :       }
     252                 :            :   }
     253                 :      38255 :   return typ_NULL;
     254                 :            : }
     255                 :            : 
     256                 :            : /*************************************************************************/
     257                 :            : /**                                                                     **/
     258                 :            : /**                           GALOIS GROUP                              **/
     259                 :            : /**                                                                     **/
     260                 :            : /*************************************************************************/
     261                 :            : 
     262                 :            : /* exchange elements i and j in vector x */
     263                 :            : static GEN
     264                 :       5992 : transroot(GEN x, int i, int j)
     265                 :       5992 : { x = leafcopy(x); swap(gel(x,i), gel(x,j)); return x; }
     266                 :            : 
     267                 :            : GEN
     268                 :       3164 : tschirnhaus(GEN x)
     269                 :            : {
     270                 :       3164 :   pari_sp av = avma, av2;
     271                 :       3164 :   long a, v = varn(x);
     272                 :       3164 :   GEN u, y = cgetg(5,t_POL);
     273                 :            : 
     274         [ -  + ]:       3164 :   if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
     275         [ -  + ]:       3164 :   if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
     276         [ -  + ]:       3164 :   if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
     277                 :       3164 :   y[1] = evalsigne(1)|evalvarn(0);
     278                 :            :   do
     279                 :            :   {
     280         [ +  + ]:       3234 :     a = random_bits(2); if (a==0) a  = 1; gel(y,4) = stoi(a);
     281         [ +  + ]:       3234 :     a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
     282         [ +  + ]:       3234 :     a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
     283                 :       3234 :     u = RgXQ_charpoly(y,x,v); av2 = avma;
     284                 :            :   }
     285         [ +  + ]:       3234 :   while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
     286         [ -  + ]:       3164 :   if (DEBUGLEVEL>1)
     287                 :          0 :     err_printf("Tschirnhaus transform. New pol: %Ps",u);
     288                 :       3164 :   avma=av2; return gerepileupto(av,u);
     289                 :            : }
     290                 :            : 
     291                 :            : /* Assume pol in Z[X], monic of degree n. Find L in Z such that
     292                 :            :  * POL = L^(-n) pol(L x) is monic in Z[X]. Return POL and set *ptk = L.
     293                 :            :  * No GC. */
     294                 :            : GEN
     295                 :      12236 : ZX_Z_normalize(GEN pol, GEN *ptk)
     296                 :            : {
     297                 :      12236 :   long i,j, sk, n = degpol(pol); /* > 0 */
     298                 :            :   GEN k, fa, P, E, a, POL;
     299                 :            : 
     300                 :      12236 :   a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
     301         [ +  + ]:      45675 :   for (i = n-2; i >= 0; i--)
     302                 :            :   {
     303                 :      40509 :     k = gcdii(k, gel(a,i));
     304 [ +  + ][ +  + ]:      40509 :     if (is_pm1(k)) { if (ptk) *ptk = gen_1; return pol; }
     305                 :            :   }
     306                 :       5166 :   sk = signe(k);
     307 [ +  + ][ +  - ]:       5166 :   if (!sk) { if (ptk) *ptk = gen_1; return pol; /* monomial! */ }
     308                 :       4186 :   fa = absi_factor_limit(k, 0); k = gen_1;
     309                 :       4186 :   P = gel(fa,1);
     310                 :       4186 :   E = gel(fa,2);
     311                 :       4186 :   POL = leafcopy(pol); a = POL+2;
     312         [ +  + ]:       9562 :   for (i = lg(P)-1; i > 0; i--)
     313                 :            :   {
     314                 :       5376 :     GEN p = gel(P,i), pv, pvj;
     315                 :       5376 :     long vmin = itos(gel(E,i));
     316                 :            :     /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
     317         [ +  + ]:      37275 :     for (j=n-1; j>=0; j--)
     318                 :            :     {
     319                 :            :       long v;
     320         [ +  + ]:      31899 :       if (!signe(gel(a,j))) continue;
     321                 :      19964 :       v = Z_pval(gel(a,j), p) / (n - j);
     322         [ +  + ]:      19964 :       if (v < vmin) vmin = v;
     323                 :            :     }
     324         [ +  + ]:       5376 :     if (!vmin) continue;
     325                 :       1050 :     pvj = pv = powiu(p,vmin); k = mulii(k, pv);
     326                 :            :     /* a[j] /= p^(v*(n-j)) */
     327         [ +  + ]:       7420 :     for (j=n-1; j>=0; j--)
     328                 :            :     {
     329         [ +  + ]:       6370 :       if (j < n-1) pvj = mulii(pvj, pv);
     330                 :       6370 :       gel(a,j) = diviiexact(gel(a,j), pvj);
     331                 :            :     }
     332                 :            :   }
     333         [ +  + ]:      12236 :   if (ptk) *ptk = k; return POL;
     334                 :            : }
     335                 :            : 
     336                 :            : /* Assume pol != 0 in Z[X]. Find C in Q, L in Z such that POL = C pol(x/L) monic
     337                 :            :  * in Z[X]. Return POL and set *pL = L. Wasteful (but correct) if pol is not
     338                 :            :  * primitive: better if caller used Q_primpart already. No GC. */
     339                 :            : GEN
     340                 :      12243 : ZX_primitive_to_monic(GEN pol, GEN *pL)
     341                 :            : {
     342                 :      12243 :   long i,j, n = degpol(pol);
     343                 :      12243 :   GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
     344                 :            : 
     345         [ +  + ]:      12243 :   if (is_pm1(lc))
     346                 :            :   {
     347         [ +  - ]:      11984 :     if (pL) *pL = gen_1;
     348         [ +  + ]:      11984 :     return signe(lc) < 0? ZX_neg(pol): pol;
     349                 :            :   }
     350         [ +  + ]:        259 :   if (signe(lc) < 0)
     351                 :         35 :     POL = ZX_neg(pol);
     352                 :            :   else
     353                 :        224 :     POL = leafcopy(pol);
     354                 :        259 :   a = POL+2; lc = gel(a,n);
     355                 :        259 :   fa = Z_factor_limit(lc,0); L = gen_1;
     356                 :        259 :   P = gel(fa,1);
     357                 :        259 :   E = gel(fa,2);
     358         [ +  + ]:        651 :   for (i = lg(P)-1; i > 0; i--)
     359                 :            :   {
     360                 :        392 :     GEN p = gel(P,i), pk, pku;
     361                 :        392 :     long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
     362                 :            : 
     363         [ +  + ]:        392 :     if (d < 0) { k++; d += n; }
     364                 :            :     /* k = ceil(e[i] / n); find d, k such that  p^d pol(x / p^k) monic */
     365         [ +  + ]:       1519 :     for (j=n-1; j>0; j--)
     366                 :            :     {
     367         [ +  + ]:       1127 :       if (!signe(gel(a,j))) continue;
     368                 :       1015 :       v = Z_pval(gel(a,j), p);
     369         [ +  + ]:       1085 :       while (v + d < k * j) { k++; d += n; }
     370                 :            :     }
     371                 :        392 :     pk = powiu(p,k); j0 = d/k;
     372                 :        392 :     L = mulii(L, pk);
     373                 :            : 
     374                 :        392 :     pku = powiu(p,d - k*j0);
     375                 :            :     /* a[j] *= p^(d - kj) */
     376         [ +  + ]:       1617 :     for (j=j0; j>=0; j--)
     377                 :            :     {
     378         [ +  + ]:       1225 :       if (j < j0) pku = mulii(pku, pk);
     379                 :       1225 :       gel(a,j) = mulii(gel(a,j), pku);
     380                 :            :     }
     381                 :        392 :     j0++;
     382                 :        392 :     pku = powiu(p,k*j0 - d);
     383                 :            :     /* a[j] /= p^(kj - d) */
     384         [ +  + ]:       1078 :     for (j=j0; j<=n; j++)
     385                 :            :     {
     386         [ +  + ]:        686 :       if (j > j0) pku = mulii(pku, pk);
     387                 :        686 :       gel(a,j) = diviiexact(gel(a,j), pku);
     388                 :            :     }
     389                 :            :   }
     390         [ +  - ]:        259 :   if (pL) *pL = L;
     391                 :      12243 :   return POL;
     392                 :            : }
     393                 :            : /* Assume pol != 0 in Z[X]. Find C,L in Q such that POL = C pol(x/L)
     394                 :            :  * monic in Z[X]. Return POL and set *pL = L.
     395                 :            :  * Wasteful (but correct) if pol is not primitive: better if caller used
     396                 :            :  * Q_primpart already. No GC. */
     397                 :            : GEN
     398                 :      11984 : ZX_Q_normalize(GEN pol, GEN *pL)
     399                 :            : {
     400                 :      11984 :   GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
     401                 :      11984 :   POL = ZX_Z_normalize(POL, pL);
     402         [ +  + ]:      11984 :   if (pL) *pL = gdiv(lc, *pL);
     403                 :      11984 :   return POL;
     404                 :            : }
     405                 :            : /* pol != 0 in Z[x], returns a monic polynomial POL in Z[x] generating the
     406                 :            :  * same field: there exist C in Q, L in Z such that POL(x) = C pol(x/L).
     407                 :            :  * Set *L = NULL if L = 1, and to L otherwise. No garbage collecting. */
     408                 :            : GEN
     409                 :          0 : ZX_to_monic(GEN pol, GEN *L)
     410                 :            : {
     411                 :          0 :   long n = lg(pol)-1;
     412                 :          0 :   GEN lc = gel(pol,n);
     413 [ #  # ][ #  # ]:          0 :   if (is_pm1(lc)) { *L = gen_1; return signe(lc) > 0? pol: ZX_neg(pol); }
     414                 :          0 :   return ZX_primitive_to_monic(Q_primpart(pol), L);
     415                 :            : }
     416                 :            : 
     417                 :            : /* x1*x2^2 + x2*x3^2 + x3*x4^2 + x4*x1^2 */
     418                 :            : static GEN
     419                 :       1092 : F4(GEN x)
     420                 :            : {
     421                 :       1092 :   return gadd(
     422                 :       4368 :     gmul(gel(x,1), gadd(gsqr(gel(x,2)), gmul(gel(x,4),gel(x,1)))),
     423                 :       4368 :     gmul(gel(x,3), gadd(gsqr(gel(x,4)), gmul(gel(x,2),gel(x,3)))));
     424                 :            : }
     425                 :            : 
     426                 :            : static GEN
     427                 :       1393 : roots_to_ZX(GEN z, long *e)
     428                 :            : {
     429                 :       1393 :   GEN a = roots_to_pol(z,0);
     430                 :       1393 :   GEN b = grndtoi(real_i(a),e);
     431                 :       1393 :   long e1 = gexpo(imag_i(a));
     432         [ +  + ]:       1393 :   if (e1 > *e) *e = e1;
     433                 :       1393 :   return b;
     434                 :            : }
     435                 :            : 
     436                 :            : static GEN
     437                 :       1057 : polgaloisnames(long a, long b)
     438                 :            : {
     439                 :       1057 :   const char * const t[]={"S1", "S2", "A3", "S3",
     440                 :            :        "C(4) = 4", "E(4) = 2[x]2", "D(4)", "A4", "S4",
     441                 :            :        "C(5) = 5", "D(5) = 5:2", "F(5) = 5:4", "A5", "S5",
     442                 :            :        "C(6) = 6 = 3[x]2", "D_6(6) = [3]2", "D(6) = S(3)[x]2",
     443                 :            :              "A_4(6) = [2^2]3", "F_18(6) = [3^2]2 = 3 wr 2",
     444                 :            :              "2A_4(6) = [2^3]3 = 2 wr 3", "S_4(6d) = [2^2]S(3)",
     445                 :            :              "S_4(6c) = 1/2[2^3]S(3)", "F_18(6):2 = [1/2.S(3)^2]2",
     446                 :            :              "F_36(6) = 1/2[S(3)^2]2", "2S_4(6) = [2^3]S(3) = 2 wr S(3)",
     447                 :            :              "L(6) = PSL(2,5) = A_5(6)", "F_36(6):2 = [S(3)^2]2 = S(3) wr 2",
     448                 :            :              "L(6):2 = PGL(2,5) = S_5(6)", "A6", "S6",
     449                 :            :        "C(7) = 7", "D(7) = 7:2", "F_21(7) = 7:3", "F_42(7) = 7:6",
     450                 :            :              "L(7) = L(3,2)", "A7", "S7"};
     451                 :            : 
     452                 :       1057 :    const long idx[]={0,1,2,4,9,14,30};
     453                 :       1057 :    return strtoGENstr(t[idx[a-1]+b-1]);
     454                 :            : }
     455                 :            : 
     456                 :            : static GEN
     457                 :       1057 : galois_res(long d, long n, long s, long k)
     458                 :            : {
     459                 :       1057 :   long kk = k;
     460                 :       1057 :   GEN z = cgetg(5,t_VEC);
     461         [ +  + ]:       1057 :   if (!new_galois_format)
     462                 :            :   {
     463         [ +  + ]:        553 :     switch (d) {
     464                 :            :       case 6:
     465 [ +  + ][ +  + ]:        238 :         kk = (k == 6 || k == 2)? 2: 1;
     466                 :        238 :         break;
     467                 :            :       default:
     468                 :        315 :         kk = 1;
     469                 :            :     }
     470                 :            :   }
     471                 :       1057 :   gel(z,1) = stoi(n);
     472                 :       1057 :   gel(z,2) = stoi(s);
     473                 :       1057 :   gel(z,3) = stoi(kk);
     474                 :       1057 :   gel(z,4) = polgaloisnames(d,k);
     475                 :       1057 :   return z;
     476                 :            : }
     477                 :            : 
     478                 :            : GEN
     479                 :       3059 : polgalois(GEN x, long prec)
     480                 :            : {
     481                 :       3059 :   pari_sp av = avma, av1;
     482                 :            :   long i,j,k,n,f,l,l2,e,e1,pr,ind;
     483                 :            :   GEN x1,p1,p2,p3,p4,p5,w,z,ee;
     484                 :       3059 :   const int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
     485                 :       3059 :   const int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
     486                 :            :                       1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
     487                 :            :                       1,2,4,6, 1,5,2,4, 1,3,2,6, 1,2,3,5, 1,4,2,3};
     488         [ -  + ]:       3059 :   if (typ(x)!=t_POL) pari_err_TYPE("galois",x);
     489                 :       3059 :   n=degpol(x);
     490         [ -  + ]:       3059 :   if (n>11) pari_err_IMPL("galois of degree higher than 11");
     491                 :       3059 :   x = Q_primpart(x);
     492                 :       3059 :   RgX_check_ZX(x, "galois");
     493         [ -  + ]:       3059 :   if (!ZX_is_irred(x)) pari_err_IRREDPOL("galois",x);
     494                 :            : 
     495         [ +  + ]:       3059 :   if (n<4)
     496                 :            :   {
     497         [ -  + ]:         91 :     if (n == 1) { avma = av; return galois_res(n,1, 1,1); }
     498         [ +  + ]:         91 :     if (n == 2) { avma = av; return galois_res(n,2,-1,1); }
     499                 :            :     /* n = 3 */
     500                 :         63 :     f = Z_issquare(ZX_disc(x));
     501                 :         63 :     avma = av;
     502         [ +  + ]:         63 :     return f? galois_res(n,3,1,1):
     503                 :            :               galois_res(n,6,-1,2);
     504                 :            :   }
     505                 :       2968 :   x1 = x = ZX_Q_normalize(x,NULL); av1=avma;
     506         [ +  + ]:       2968 :   if (n > 7) return galoisbig(x, prec);
     507                 :            :   for(;;)
     508                 :            :   {
     509                 :       1183 :     double fb = fujiwara_bound(x);
     510   [ +  +  +  +  :       1183 :     switch(n)
                      - ]
     511                 :            :     {
     512                 :        182 :       case 4: z = cgetg(7,t_VEC);
     513                 :        182 :         prec = nbits2prec((long)(fb*18.) + 64);
     514                 :            :         for(;;)
     515                 :            :         {
     516                 :        182 :           p1=QX_complex_roots(x,prec);
     517                 :        182 :           gel(z,1) = F4(p1);
     518                 :        182 :           gel(z,2) = F4(transroot(p1,1,2));
     519                 :        182 :           gel(z,3) = F4(transroot(p1,1,3));
     520                 :        182 :           gel(z,4) = F4(transroot(p1,1,4));
     521                 :        182 :           gel(z,5) = F4(transroot(p1,2,3));
     522                 :        182 :           gel(z,6) = F4(transroot(p1,3,4));
     523         [ +  - ]:        182 :           p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     524                 :          0 :           prec = precdbl(prec);
     525                 :          0 :         }
     526         [ +  + ]:        182 :         if (!ZX_is_squarefree(p5)) goto tchi;
     527                 :        147 :         p2 = gel(ZX_factor(p5),1);
     528   [ +  +  +  - ]:        147 :         switch(lg(p2)-1)
     529                 :            :         {
     530                 :         56 :           case 1: f = Z_issquare(ZX_disc(x)); avma = av;
     531         [ +  + ]:         56 :             return f? galois_res(n,12,1,4): galois_res(n,24,-1,5);
     532                 :            : 
     533                 :         35 :           case 2: avma = av; return galois_res(n,8,-1,3);
     534                 :            : 
     535                 :         56 :           case 3: avma = av;
     536                 :        112 :             return (degpol(gel(p2,1))==2)? galois_res(n,4,1,2)
     537         [ +  + ]:         56 :                                          : galois_res(n,4,-1,1);
     538                 :            : 
     539                 :          0 :           default: pari_err_BUG("galois (bug1)");
     540                 :            :         }
     541                 :            : 
     542                 :        175 :       case 5: z = cgetg(7,t_VEC);
     543                 :        175 :         ee= cgetg(7,t_VECSMALL);
     544                 :        175 :         w = cgetg(7,t_VECSMALL);
     545                 :        175 :         prec = nbits2prec((long)(fb*21.) + 64);
     546                 :            :         for(;;)
     547                 :            :         {
     548                 :            :           for(;;)
     549                 :            :           {
     550                 :        175 :             p1=QX_complex_roots(x,prec);
     551         [ +  + ]:       1225 :             for (l=1; l<=6; l++)
     552                 :            :             {
     553 [ +  + ][ +  + ]:       1050 :               p2=(l==1)?p1: ((l<6)?transroot(p1,1,l): transroot(p1,2,5));
     554                 :       1050 :               p3=gen_0;
     555         [ +  + ]:       6300 :               for (k=0,i=1; i<=5; i++,k+=4)
     556                 :            :               {
     557                 :       5250 :                 p5 = gadd(gmul(gel(p2,ind5[k]),gel(p2,ind5[k+1])),
     558                 :      10500 :                           gmul(gel(p2,ind5[k+2]),gel(p2,ind5[k+3])));
     559                 :       5250 :                 p3 = gadd(p3, gmul(gsqr(gel(p2,i)),p5));
     560                 :            :               }
     561                 :       1050 :               gel(w,l) = grndtoi(real_i(p3),&e);
     562         [ +  + ]:       1050 :               e1 = gexpo(imag_i(p3)); if (e1>e) e=e1;
     563                 :       1050 :               ee[l]=e; gel(z,l) = p3;
     564                 :            :             }
     565         [ +  - ]:        175 :             p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     566                 :          0 :             prec = precdbl(prec);
     567                 :          0 :           }
     568         [ -  + ]:        175 :           if (!ZX_is_squarefree(p5)) goto tchi;
     569                 :        175 :           p3=gel(ZX_factor(p5),1);
     570                 :        175 :           f=Z_issquare(ZX_disc(x));
     571         [ +  + ]:        175 :           if (lg(p3)-1==1)
     572                 :            :           {
     573                 :         63 :             avma = av;
     574         [ +  + ]:         63 :             return f? galois_res(n,60,1,4): galois_res(n,120,-1,5);
     575                 :            :           }
     576         [ +  + ]:        112 :           if (!f) { avma = av; return galois_res(n,20,-1,3); }
     577                 :            : 
     578                 :         84 :           pr = - (prec2nbits(prec) >> 1);
     579         [ +  - ]:        315 :           for (l=1; l<=6; l++)
     580 [ +  + ][ +  - ]:        315 :             if (ee[l] <= pr && gequal0(poleval(p5,gel(w,l)))) break;
     581         [ -  + ]:         84 :           if (l>6) pari_err_BUG("galois (bug4)");
     582         [ -  + ]:         84 :           p2=(l==6)? transroot(p1,2,5):transroot(p1,1,l);
     583                 :         84 :           p3=gen_0;
     584         [ +  + ]:        504 :           for (i=1; i<=5; i++)
     585                 :            :           {
     586         [ +  + ]:        420 :             j = (i == 5)? 1: i+1;
     587                 :        420 :             p3 = gadd(p3,gmul(gmul(gel(p2,i),gel(p2,j)),
     588                 :        840 :                               gsub(gel(p2,j),gel(p2,i))));
     589                 :            :           }
     590                 :         84 :           p5=gsqr(p3); p4=grndtoi(real_i(p5),&e);
     591         [ +  + ]:         84 :           e1 = gexpo(imag_i(p5)); if (e1>e) e=e1;
     592         [ +  - ]:         84 :           if (e <= -10)
     593                 :            :           {
     594         [ +  + ]:         84 :             if (gequal0(p4)) goto tchi;
     595                 :         63 :             f = Z_issquare(p4); avma = av;
     596         [ +  + ]:         63 :             return f? galois_res(n,5,1,1): galois_res(n,10,1,2);
     597                 :            :           }
     598                 :          0 :           prec = precdbl(prec);
     599                 :          0 :         }
     600                 :            : 
     601                 :        623 :       case 6: z = cgetg(7, t_VEC);
     602                 :        623 :         prec = nbits2prec((long) (fb * 42) + 64);
     603                 :            :         for(;;)
     604                 :            :         {
     605                 :            :           for(;;)
     606                 :            :           {
     607                 :        623 :             p1=QX_complex_roots(x,prec);
     608         [ +  + ]:       4361 :             for (l=1; l<=6; l++)
     609                 :            :             {
     610         [ +  + ]:       3738 :               p2=(l==1)?p1:transroot(p1,1,l);
     611                 :       3738 :               p3=gen_0; k=0;
     612 [ +  + ][ +  + ]:      78498 :               for (i=1; i<=5; i++) for (j=i+1; j<=6; j++)
     613                 :            :               {
     614                 :      56070 :                 p5=gadd(gmul(gel(p2,ind6[k]),gel(p2,ind6[k+1])),
     615                 :     112140 :                         gmul(gel(p2,ind6[k+2]),gel(p2,ind6[k+3])));
     616                 :      56070 :                 p3=gadd(p3,gmul(gsqr(gmul(gel(p2,i),gel(p2,j))),p5));
     617                 :      56070 :                 k += 4;
     618                 :            :               }
     619                 :       3738 :               gel(z,l) = p3;
     620                 :            :             }
     621         [ +  - ]:        623 :             p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     622                 :          0 :             prec = precdbl(prec);
     623                 :          0 :           }
     624         [ +  + ]:        623 :           if (!ZX_is_squarefree(p5)) goto tchi;
     625                 :        462 :           p2=gel(ZX_factor(p5),1);
     626   [ +  +  +  +  :        462 :           switch(lg(p2)-1)
                      - ]
     627                 :            :           {
     628                 :            :             case 1:
     629                 :        112 :               z = cgetg(11,t_VEC); ind=0;
     630                 :        112 :               p3=gadd(gmul(gmul(gel(p1,1),gel(p1,2)),gel(p1,3)),
     631                 :        336 :                       gmul(gmul(gel(p1,4),gel(p1,5)),gel(p1,6)));
     632                 :        112 :               gel(z,++ind) = p3;
     633         [ +  + ]:        448 :               for (i=1; i<=3; i++)
     634         [ +  + ]:       1344 :                 for (j=4; j<=6; j++)
     635                 :            :                 {
     636                 :       1008 :                   p2=transroot(p1,i,j);
     637                 :       1008 :                   p3=gadd(gmul(gmul(gel(p2,1),gel(p2,2)),gel(p2,3)),
     638                 :       3024 :                           gmul(gmul(gel(p2,4),gel(p2,5)),gel(p2,6)));
     639                 :       1008 :                   gel(z,++ind) = p3;
     640                 :            :                 }
     641                 :        112 :               p5 = roots_to_ZX(z, &e);
     642         [ +  - ]:        112 :               if (e <= -10)
     643                 :            :               {
     644         [ -  + ]:        112 :                 if (!ZX_is_squarefree(p5)) goto tchi;
     645                 :        112 :                 p2 = gel(ZX_factor(p5),1);
     646                 :        112 :                 f = Z_issquare(ZX_disc(x));
     647                 :        112 :                 avma = av;
     648         [ +  + ]:        112 :                 if (lg(p2)-1==1)
     649         [ +  + ]:         56 :                   return f? galois_res(n,360,1,15): galois_res(n,720,-1,16);
     650                 :            :                 else
     651         [ +  + ]:         56 :                   return f? galois_res(n,36,1,10): galois_res(n,72,-1,13);
     652                 :            :               }
     653                 :          0 :               prec = precdbl(prec); break;
     654                 :            : 
     655         [ -  + ]:        196 :             case 2: l2=degpol(gel(p2,1)); if (l2>3) l2=6-l2;
     656   [ +  +  +  - ]:        196 :               switch(l2)
     657                 :            :               {
     658                 :         56 :                 case 1: f = Z_issquare(ZX_disc(x));
     659                 :         56 :                   avma = av;
     660         [ +  + ]:         56 :                   return f? galois_res(n,60,1,12): galois_res(n,120,-1,14);
     661                 :         84 :                 case 2: f = Z_issquare(ZX_disc(x));
     662         [ +  + ]:         84 :                   if (f) { avma = av; return galois_res(n,24,1,7); }
     663         [ +  - ]:         56 :                   p3 = (degpol(gel(p2,1))==2)? gel(p2,2): gel(p2,1);
     664                 :         56 :                   f = Z_issquare(ZX_disc(p3));
     665                 :         56 :                   avma = av;
     666         [ +  + ]:         56 :                   return f? galois_res(n,24,-1,6): galois_res(n,48,-1,11);
     667         [ +  - ]:         64 :                 case 3: f = Z_issquare(ZX_disc(gel(p2,1)))
     668         [ +  - ]:         56 :                          || Z_issquare(ZX_disc(gel(p2,2)));
           [ +  +  +  + ]
     669                 :         56 :                   avma = av;
     670         [ +  + ]:         56 :                   return f? galois_res(n,18,-1,5): galois_res(n,36,-1,9);
     671                 :            :               }
     672                 :            :             case 3:
     673         [ +  + ]:        504 :               for (l2=1; l2<=3; l2++)
     674         [ +  + ]:        378 :                 if (degpol(gel(p2,l2)) >= 3) p3 = gel(p2,l2);
     675         [ +  + ]:        126 :               if (degpol(p3) == 3)
     676                 :            :               {
     677                 :         63 :                 f = Z_issquare(ZX_disc(p3)); avma = av;
     678         [ +  + ]:         63 :                 return f? galois_res(n,6,-1,1): galois_res(n,12,-1,3);
     679                 :            :               }
     680                 :            :               else
     681                 :            :               {
     682                 :         63 :                 f = Z_issquare(ZX_disc(x)); avma = av;
     683         [ +  + ]:         63 :                 return f? galois_res(n,12,1,4): galois_res(n,24,-1,8);
     684                 :            :               }
     685                 :         28 :             case 4: avma = av; return galois_res(n,6,-1,2);
     686                 :          0 :             default: pari_err_BUG("galois (bug3)");
     687                 :            :           }
     688                 :          0 :         }
     689                 :            : 
     690                 :        203 :       case 7: z = cgetg(36,t_VEC);
     691                 :        203 :         prec = nbits2prec((long)(fb*7.) + 64);
     692                 :            :         for(;;)
     693                 :            :         {
     694                 :        301 :           ind = 0; p1=QX_complex_roots(x,prec);
     695         [ +  + ]:       1806 :           for (i=1; i<=5; i++)
     696         [ +  + ]:       6020 :             for (j=i+1; j<=6; j++)
     697                 :            :             {
     698                 :       4515 :               GEN t = gadd(gel(p1,i),gel(p1,j));
     699         [ +  + ]:      15050 :               for (k=j+1; k<=7; k++) gel(z,++ind) = gadd(t, gel(p1,k));
     700                 :            :             }
     701         [ +  + ]:        301 :           p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     702                 :         98 :           prec = precdbl(prec);
     703                 :         98 :         }
     704         [ -  + ]:        203 :         if (!ZX_is_squarefree(p5)) goto tchi;
     705                 :        203 :         p2=gel(ZX_factor(p5),1);
     706   [ +  +  +  +  :        203 :         switch(lg(p2)-1)
                   +  - ]
     707                 :            :         {
     708                 :         56 :           case 1: f = Z_issquare(ZX_disc(x)); avma = av;
     709         [ +  + ]:         56 :             return f? galois_res(n,2520,1,6): galois_res(n,5040,-1,7);
     710                 :         63 :           case 2: f = degpol(gel(p2,1)); avma = av;
     711 [ +  + ][ -  + ]:         63 :             return (f==7 || f==28)? galois_res(n,168,1,5): galois_res(n,42,-1,4);
     712                 :         28 :           case 3: avma = av; return galois_res(n,21,1,3);
     713                 :         28 :           case 4: avma = av; return galois_res(n,14,-1,2);
     714                 :         28 :           case 5: avma = av; return galois_res(n,7,1,1);
     715                 :          0 :           default: pari_err_BUG("galois (bug2)");
     716                 :            :         }
     717                 :            :     }
     718                 :       3276 :     tchi: avma = av1; x = tschirnhaus(x1);
     719                 :        217 :   }
     720                 :            : }
     721                 :            : 
     722                 :            : #undef _res
     723                 :            : 
     724                 :            : /* Evaluate pol in s using nfelt arithmetic and Horner rule */
     725                 :            : GEN
     726                 :      11417 : nfpoleval(GEN nf, GEN pol, GEN s)
     727                 :            : {
     728                 :      11417 :   pari_sp av=avma;
     729                 :      11417 :   long i=lg(pol)-1;
     730                 :            :   GEN res;
     731         [ -  + ]:      11417 :   if (i==1) return gen_0;
     732                 :      11417 :   res = nf_to_scalar_or_basis(nf, gel(pol,i));
     733         [ +  + ]:      28609 :   for (i-- ; i>=2; i--)
     734                 :      17192 :     res = nfadd(nf, nfmul(nf, s, res), gel(pol,i));
     735                 :      11417 :   return gerepileupto(av, res);
     736                 :            : }
     737                 :            : 
     738                 :            : static GEN
     739                 :      19558 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
     740                 :            : {
     741                 :      19558 :   pari_sp av = avma;
     742                 :      19558 :   long i = lg(pol)-1;
     743                 :            :   GEN res, den;
     744         [ -  + ]:      19558 :   if (i==1) return gen_0;
     745                 :      19558 :   pol = Q_remove_denom(pol, &den);
     746                 :      19558 :   res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
     747         [ +  + ]:     105014 :   for (i-- ; i>=2; i--)
     748                 :      85456 :     res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
     749         [ +  + ]:      19558 :   if (den) res = RgC_Rg_div(res, den);
     750                 :      19558 :   return gerepileupto(av, res);
     751                 :            : }
     752                 :            : 
     753                 :            : GEN
     754                 :       1456 : FpX_FpC_nfpoleval(GEN nf, GEN pol, GEN a, GEN p)
     755                 :            : {
     756                 :       1456 :   pari_sp av=avma;
     757                 :       1456 :   long i=lg(pol)-1, n=nf_get_degree(nf);
     758                 :            :   GEN res, Ma;
     759         [ -  + ]:       1456 :   if (i==1) return zerocol(n);
     760                 :       1456 :   Ma = FpM_red(zk_multable(nf, a), p);
     761                 :       1456 :   res = scalarcol(gel(pol,i),n);
     762         [ +  + ]:       3612 :   for (i-- ; i>=2; i--)
     763                 :            :   {
     764                 :       2156 :     res = FpM_FpC_mul(Ma, res, p);
     765                 :       2156 :     gel(res,1) = Fp_add(gel(res,1), gel(pol,i), p);
     766                 :            :   }
     767                 :       1456 :   return gerepileupto(av, res);
     768                 :            : }
     769                 :            : 
     770                 :            : /* compute s(x), not stack clean */
     771                 :            : static GEN
     772                 :       3738 : table_galoisapply(GEN nf, GEN m, GEN x)
     773                 :            : {
     774                 :       3738 :   x = nf_to_scalar_or_alg(nf, x);
     775         [ +  + ]:       3738 :   if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
     776                 :       3738 :   return QX_table_nfpoleval(nf, x, m);
     777                 :            : }
     778                 :            : 
     779                 :            : /* compute s(x), not stack clean */
     780                 :            : static GEN
     781                 :      10087 : ZC_galoisapply(GEN nf, GEN s, GEN x)
     782                 :            : {
     783                 :      10087 :   x = nf_to_scalar_or_alg(nf, x);
     784         [ +  + ]:      10087 :   if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
     785                 :      10087 :   return QX_table_nfpoleval(nf, x, zk_multable(nf, s));
     786                 :            : }
     787                 :            : 
     788                 :            : static GEN
     789                 :       1456 : QX_galoisapplymod(GEN nf, GEN pol, GEN S, GEN p)
     790                 :            : {
     791                 :       1456 :   GEN den, P = Q_remove_denom(pol,&den);
     792                 :            :   GEN pe, pe1, denpe, R;
     793         [ +  + ]:       1456 :   if (den)
     794                 :            :   {
     795                 :         98 :     ulong e = Z_pval(den, p);
     796                 :         98 :     pe = powiu(p, e); pe1 = mulii(pe, p);
     797                 :         98 :     denpe = Fp_inv(diviiexact(den, pe), pe1);
     798                 :            :   } else {
     799                 :       1358 :     pe = gen_1; pe1 = p; denpe = gen_1;
     800                 :            :   }
     801                 :       1456 :   R = FpX_FpC_nfpoleval(nf, FpX_red(P, pe1), FpC_red(S, pe1), pe1);
     802                 :       1456 :   return gdivexact(FpC_Fp_mul(R, denpe, pe1), pe);
     803                 :            : }
     804                 :            : 
     805                 :            : static GEN
     806                 :          7 : pr_galoisapply(GEN nf, GEN pr, GEN aut)
     807                 :            : {
     808                 :            :   GEN p, t, u;
     809         [ -  + ]:          7 :   if (typ(pr_get_tau(pr)) == t_INT) return pr; /* inert */
     810                 :          7 :   p = pr_get_p(pr);
     811                 :          7 :   u = QX_galoisapplymod(nf, coltoliftalg(nf, pr_get_gen(pr)), aut, p);
     812                 :          7 :   t = FpM_deplin(zk_multable(nf, u), p);
     813                 :          7 :   t = zk_scalar_or_multable(nf, t);
     814                 :          7 :   return mkvec5(p, u, gel(pr,3), gel(pr,4), t);
     815                 :            : }
     816                 :            : 
     817                 :            : static GEN
     818                 :          7 : vecgaloisapply(GEN nf, GEN aut, GEN v)
     819                 :            : {
     820                 :            :   long i, l;
     821                 :          7 :   GEN V = cgetg_copy(v, &l);
     822         [ +  + ]:         21 :   for (i = 1; i < l; i++) gel(V,i) = galoisapply(nf, aut, gel(v,i));
     823                 :          7 :   return V;
     824                 :            : }
     825                 :            : 
     826                 :            : /* x: famat or standard algebraic number, aut automorphism in ZC form
     827                 :            :  * simplified from general galoisapply */
     828                 :            : static GEN
     829                 :         49 : elt_galoisapply(GEN nf, GEN aut, GEN x)
     830                 :            : {
     831                 :         49 :   pari_sp av = avma;
     832   [ +  +  +  +  :         49 :   switch(typ(x))
                +  +  - ]
     833                 :            :   {
     834                 :          7 :     case t_INT:  return icopy(x);
     835                 :          7 :     case t_FRAC: return gcopy(x);
     836                 :          7 :     case t_POLMOD: x = gel(x,2); /* fall through */
     837                 :            :     case t_POL: {
     838                 :         14 :       GEN y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
     839                 :         14 :       return gerepileupto(av,y);
     840                 :            :     }
     841                 :            :     case t_COL:
     842                 :          7 :       return gerepileupto(av, ZC_galoisapply(nf, aut, x));
     843                 :            :     case t_MAT:
     844      [ +  +  - ]:         14 :       switch(lg(x)) {
     845                 :          7 :         case 1: return cgetg(1, t_MAT);
     846                 :          7 :         case 3: retmkmat2(vecgaloisapply(nf,aut,gel(x,1)), ZC_copy(gel(x,2)));
     847                 :            :       }
     848                 :            :   }
     849                 :          0 :   pari_err_TYPE("galoisapply",x);
     850                 :         49 :   return NULL; /* not reached */
     851                 :            : }
     852                 :            : 
     853                 :            : GEN
     854                 :       5054 : galoisapply(GEN nf, GEN aut, GEN x)
     855                 :            : {
     856                 :       5054 :   pari_sp av = avma;
     857                 :            :   long lx, j;
     858                 :            :   GEN y;
     859                 :            : 
     860                 :       5054 :   nf = checknf(nf);
     861   [ +  +  +  +  :       5054 :   switch(typ(x))
             +  +  +  - ]
     862                 :            :   {
     863                 :         70 :     case t_INT:  return icopy(x);
     864                 :          7 :     case t_FRAC: return gcopy(x);
     865                 :            : 
     866                 :         35 :     case t_POLMOD: x = gel(x,2); /* fall through */
     867                 :            :     case t_POL:
     868                 :        448 :       aut = algtobasis(nf, aut);
     869                 :        448 :       y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
     870                 :        448 :       return gerepileupto(av,y);
     871                 :            : 
     872                 :            :     case t_VEC:
     873                 :         56 :       aut = algtobasis(nf, aut);
     874      [ +  +  - ]:         56 :       switch(lg(x))
     875                 :            :       {
     876                 :          7 :         case 6: return gerepilecopy(av, pr_galoisapply(nf, x, aut));
     877                 :         49 :         case 3: y = cgetg(3,t_VEC);
     878                 :         49 :           gel(y,1) = galoisapply(nf, aut, gel(x,1));
     879                 :         49 :           gel(y,2) = elt_galoisapply(nf, aut, gel(x,2));
     880                 :         49 :           return gerepileupto(av, y);
     881                 :            :       }
     882                 :          0 :       break;
     883                 :            : 
     884                 :            :     case t_COL:
     885                 :       3542 :       aut = algtobasis(nf, aut);
     886                 :       3542 :       return gerepileupto(av, ZC_galoisapply(nf, aut, x));
     887                 :            : 
     888                 :            :     case t_MAT: /* ideal */
     889         [ -  + ]:        931 :       lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
     890         [ -  + ]:        931 :       if (nbrows(x) != nf_get_degree(nf)) break;
     891                 :        931 :       aut = zk_multable(nf, algtobasis(nf, aut));
     892                 :        931 :       y = cgetg(lx,t_MAT);
     893         [ +  + ]:       4669 :       for (j=1; j<lx; j++) gel(y,j) = table_galoisapply(nf, aut, gel(x,j));
     894                 :        931 :       return gerepileupto(av, idealhnf_shallow(nf,y));
     895                 :            :   }
     896                 :          0 :   pari_err_TYPE("galoisapply",x);
     897                 :       5054 :   return NULL; /* not reached */
     898                 :            : }
     899                 :            : 
     900                 :            : GEN
     901                 :       1904 : nfgaloismatrix(GEN nf, GEN s)
     902                 :            : {
     903                 :            :   GEN zk, M, m;
     904                 :            :   long k, l;
     905                 :       1904 :   nf = checknf(nf);
     906                 :       1904 :   zk = nf_get_zk(nf);
     907         [ -  + ]:       1904 :   if (typ(s) != t_COL) s = algtobasis(nf, s); /* left on stack for efficiency */
     908                 :       1904 :   m = zk_multable(nf, s);
     909                 :       1904 :   l = lg(s); M = cgetg(l, t_MAT);
     910                 :       1904 :   gel(M, 1) = col_ei(l-1, 1); /* s(1) = 1 */
     911         [ +  + ]:       8631 :   for (k = 2; k < l; k++)
     912                 :       6727 :     gel(M, k) = QX_table_nfpoleval(nf, gel(zk, k), m);
     913                 :       1904 :   return M;
     914                 :            : }
     915                 :            : 
     916                 :            : static GEN
     917                 :       3017 : idealquasifrob(GEN nf, GEN gal, GEN grp, GEN pr, GEN subg, GEN *S, GEN aut)
     918                 :            : {
     919                 :       3017 :   pari_sp av = avma;
     920                 :       3017 :   long i, n = nf_get_degree(nf), f = pr_get_f(pr);
     921                 :       3017 :   GEN pi = pr_get_gen(pr);
     922         [ +  - ]:      15274 :   for (i=1; i<=n; i++)
     923                 :            :   {
     924                 :      15274 :     GEN g = gel(grp,i);
     925 [ +  + ][ +  + ]:      15274 :     if ((!subg && perm_order(g)==f)
     926 [ +  + ][ +  + ]:       9478 :       || (subg && perm_relorder(g, subg)==f))
     927                 :            :     {
     928         [ +  + ]:       5824 :       *S = aut ? gel(aut, i): poltobasis(nf, galoispermtopol(gal, g));
     929         [ +  + ]:       5824 :       if (ZC_prdvd(nf, ZC_galoisapply(nf, *S, pi), pr)) return g;
     930                 :       2807 :       avma = av;
     931                 :            :     }
     932                 :            :   }
     933                 :          0 :   pari_err_BUG("idealquasifrob [Frobenius not found]");
     934                 :       3017 :   return NULL; /*NOT REACHED*/
     935                 :            : }
     936                 :            : 
     937                 :            : GEN
     938                 :         14 : nfgaloispermtobasis(GEN nf, GEN gal)
     939                 :            : {
     940                 :         14 :   GEN grp = gal_get_group(gal);
     941                 :         14 :   long i, n = lg(grp)-1;
     942                 :         14 :   GEN aut = cgetg(n+1, t_VEC);
     943         [ +  + ]:        126 :   for(i=1; i<=n; i++)
     944                 :        112 :     gel(aut, i) = poltobasis(nf, galoispermtopol(gal, gel(grp, i)));
     945                 :         14 :   return aut;
     946                 :            : }
     947                 :            : 
     948                 :            : static void
     949                 :        182 : gal_check_pol(const char *f, GEN x, GEN y)
     950         [ -  + ]:        182 : { if (!RgX_equal_var(x,y)) pari_err_MODULUS(f,x,y); }
     951                 :            : 
     952                 :            : GEN
     953                 :       3024 : idealfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN aut)
     954                 :            : {
     955                 :       3024 :   pari_sp av = avma;
     956                 :       3024 :   GEN S=NULL, g=NULL; /*-Wall*/
     957                 :            :   GEN T, p, a, b, modpr;
     958                 :            :   long f, n, s;
     959                 :       3024 :   f = pr_get_f(pr); n = nf_get_degree(nf);
     960         [ +  + ]:       3024 :   if (f==1) { avma = av; return identity_perm(n); }
     961                 :       2933 :   g = idealquasifrob(nf, gal, gal_get_group(gal), pr, NULL, &S, aut);
     962         [ +  + ]:       2933 :   if (f==2) return gerepileupto(av, g);
     963                 :       1400 :   modpr = zk_to_Fq_init(nf,&pr,&T,&p);
     964                 :       1400 :   a = pol_x(nf_get_varn(nf));
     965                 :       1400 :   b = nf_to_Fq(nf, QX_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
     966         [ +  + ]:       3003 :   for (s = 1; s < f-1; s++)
     967                 :            :   {
     968                 :       2709 :     a = Fq_pow(a, p, T, p);
     969         [ +  + ]:       2709 :     if (ZX_equal(a, b)) break;
     970                 :            :   }
     971                 :       1400 :   g = perm_pow(g, Fl_inv(s, f));
     972                 :       3024 :   return gerepileupto(av, g);
     973                 :            : }
     974                 :            : 
     975                 :            : GEN
     976                 :         63 : idealfrobenius(GEN nf, GEN gal, GEN pr)
     977                 :            : {
     978                 :         63 :   nf = checknf(nf);
     979                 :         63 :   checkgal(gal);
     980                 :         63 :   checkprid(pr);
     981                 :         63 :   gal_check_pol("idealfrobenius",nf_get_pol(nf),gal_get_pol(gal));
     982         [ +  + ]:         63 :   if (pr_get_e(pr)>1) pari_err_DOMAIN("idealfrobenius","pr.e", ">", gen_1,pr);
     983                 :         56 :   return idealfrobenius_aut(nf, gal, pr, NULL);
     984                 :            : }
     985                 :            : 
     986                 :            : GEN
     987                 :         14 : idealramfrobenius(GEN nf, GEN gal, GEN pr, GEN ram)
     988                 :            : {
     989                 :         14 :   pari_sp av = avma;
     990                 :         14 :   GEN S=NULL, g=NULL; /*-Wall*/
     991                 :            :   GEN T, p, a, b, modpr;
     992                 :            :   GEN isog, deco;
     993                 :            :   long f, n, s;
     994                 :         14 :   f = pr_get_f(pr); n = nf_get_degree(nf);
     995         [ +  - ]:         14 :   if (f==1) { avma = av; return identity_perm(n); }
     996                 :          0 :   modpr = zk_to_Fq_init(nf,&pr,&T,&p);
     997                 :          0 :   deco = group_elts(gel(ram,1), nf_get_degree(nf));
     998                 :          0 :   isog = group_set(gel(ram,2),  nf_get_degree(nf));
     999                 :          0 :   g = idealquasifrob(nf, gal, deco, pr, isog, &S, NULL);
    1000                 :          0 :   a = pol_x(nf_get_varn(nf));
    1001                 :          0 :   b = nf_to_Fq(nf, QX_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
    1002         [ #  # ]:          0 :   for (s=0; !ZX_equal(a, b); s++)
    1003                 :          0 :     a = Fq_pow(a, p, T, p);
    1004                 :          0 :   g = perm_pow(g, Fl_inv(s, f));
    1005                 :         14 :   return gerepileupto(av, g);
    1006                 :            : }
    1007                 :            : 
    1008                 :            : static GEN
    1009                 :         42 : idealinertiagroup(GEN nf, GEN gal, GEN pr)
    1010                 :            : {
    1011                 :         42 :   long i, n = nf_get_degree(nf);
    1012                 :         42 :   GEN p, T, modpr = zk_to_Fq_init(nf,&pr,&T,&p);
    1013                 :         42 :   GEN b = modpr_genFq(modpr);
    1014                 :         42 :   long e = pr_get_e(pr), coprime = cgcd(e, pr_get_f(pr)) == 1;
    1015                 :         42 :   GEN grp = gal_get_group(gal), pi = pr_get_gen(pr);
    1016                 :         42 :   pari_sp ltop = avma;
    1017         [ +  - ]:        322 :   for (i=1; i<=n; i++)
    1018                 :            :   {
    1019                 :        322 :     GEN iso = gel(grp,i);
    1020         [ +  + ]:        322 :     if (perm_order(iso) == e)
    1021                 :            :     {
    1022                 :         98 :       GEN S = poltobasis(nf, galoispermtopol(gal, iso));
    1023         [ +  + ]:         98 :       if (ZC_prdvd(nf, ZC_galoisapply(nf, S, pi), pr)
    1024 [ +  + ][ +  - ]:         42 :           && (coprime || gequalX(nf_to_Fq(nf, galoisapply(nf,S,b), modpr))))
    1025                 :         42 :           return iso;
    1026                 :         56 :       avma = ltop;
    1027                 :            :     }
    1028                 :            :   }
    1029                 :          0 :   pari_err_BUG("idealinertiagroup [no isotropic element]");
    1030                 :         42 :   return NULL;
    1031                 :            : }
    1032                 :            : 
    1033                 :            : static GEN
    1034                 :        105 : idealramgroupstame(GEN nf, GEN gal, GEN pr)
    1035                 :            : {
    1036                 :        105 :   pari_sp av = avma;
    1037                 :            :   GEN iso, frob, giso, isog, S, res;
    1038                 :        105 :   long e = pr_get_e(pr), f = pr_get_f(pr);
    1039         [ +  + ]:        105 :   if (e == 1)
    1040                 :            :   {
    1041         [ -  + ]:         63 :     if (f==1)
    1042                 :          0 :       return cgetg(1,t_VEC);
    1043                 :         63 :     frob = idealquasifrob(nf, gal, gal_get_group(gal), pr, NULL, &S, NULL);
    1044                 :         63 :     avma = av;
    1045                 :         63 :     res = cgetg(2, t_VEC);
    1046                 :         63 :     gel(res, 1) = cyclicgroup(frob, f);
    1047                 :         63 :     return res;
    1048                 :            :   }
    1049                 :         42 :   res = cgetg(3, t_VEC);
    1050                 :         42 :   av = avma;
    1051                 :         42 :   iso = idealinertiagroup(nf, gal, pr);
    1052                 :         42 :   avma = av;
    1053                 :         42 :   giso = cyclicgroup(iso, e);
    1054                 :         42 :   gel(res, 2) = giso;
    1055         [ +  + ]:         42 :   if (f==1)
    1056                 :            :   {
    1057                 :         21 :     gel(res, 1) = giso;
    1058                 :         21 :     return res;
    1059                 :            :   }
    1060                 :         21 :   av = avma;
    1061                 :         21 :   isog = group_set(giso, nf_get_degree(nf));
    1062                 :         21 :   frob = idealquasifrob(nf, gal, gal_get_group(gal), pr, isog, &S, NULL);
    1063                 :         21 :   avma = av;
    1064                 :         21 :   gel(res, 1) = dicyclicgroup(iso,frob,e,f);
    1065                 :        105 :   return res;
    1066                 :            : }
    1067                 :            : 
    1068                 :            : static GEN
    1069                 :         14 : idealramgroupindex(GEN nf, GEN gal, GEN pr)
    1070                 :            : {
    1071                 :         14 :   pari_sp av = avma;
    1072                 :            :   GEN p, T, g, idx, modpr;
    1073                 :            :   long i, e, f, n;
    1074                 :            :   ulong nt,rorder;
    1075                 :         14 :   GEN grp = vecvecsmall_sort(gal_get_group(gal));
    1076                 :         14 :   e = pr_get_e(pr); f = pr_get_f(pr); n = nf_get_degree(nf);
    1077                 :         14 :   modpr = zk_to_Fq_init(nf,&pr,&T,&p);
    1078                 :         14 :   (void) u_pvalrem(n,p,&nt);
    1079                 :         14 :   rorder = e*f*(n/nt);
    1080                 :         14 :   idx = const_vecsmall(n,-1);
    1081                 :         14 :   g = modpr_genFq(modpr);
    1082         [ +  + ]:        266 :   for (i=2; i<=n; i++)
    1083                 :            :   {
    1084                 :            :     GEN iso;
    1085                 :            :     long o;
    1086         [ -  + ]:        252 :     if (idx[i]>=0) continue;
    1087                 :        252 :     iso = gel(grp,i); o = perm_order(iso);
    1088         [ +  + ]:        252 :     if (rorder%o == 0)
    1089                 :            :     {
    1090                 :        154 :       GEN piso = iso;
    1091                 :        154 :       GEN S = poltobasis(nf, galoispermtopol(gal, iso));
    1092                 :        154 :       GEN pi = pr_get_gen(pr);
    1093                 :        154 :       GEN spi = ZC_galoisapply(nf, S, pi);
    1094                 :            :       long j;
    1095                 :        154 :       idx[i] = idealval(nf, gsub(spi,pi), pr);
    1096         [ +  + ]:        154 :       if (idx[i] >=1)
    1097                 :            :       {
    1098         [ +  + ]:         56 :         if (f>1)
    1099                 :            :         {
    1100                 :         49 :           GEN b = nf_to_Fq(nf, QX_galoisapplymod(nf, g, S, p), modpr);
    1101         [ +  + ]:         49 :           if (!gequalX(b)) idx[i] = 0;
    1102                 :            :         }
    1103                 :            :       }
    1104                 :         98 :       else idx[i] = -1;
    1105         [ -  + ]:        154 :       for(j=2;j<o;j++)
    1106                 :            :       {
    1107                 :          0 :         piso = perm_mul(piso,iso);
    1108         [ #  # ]:          0 :         if(cgcd(j,o)==1) idx[piso[1]] = idx[i];
    1109                 :            :       }
    1110                 :            :     }
    1111                 :            :   }
    1112                 :         14 :   return gerepileuptoleaf(av, idx);
    1113                 :            : }
    1114                 :            : 
    1115                 :            : GEN
    1116                 :        119 : idealramgroups(GEN nf, GEN gal, GEN pr)
    1117                 :            : {
    1118                 :        119 :   pari_sp av = avma;
    1119                 :            :   GEN tbl, idx, res, set, sub;
    1120                 :            :   long i, j, e, n, maxm, p;
    1121                 :            :   ulong et;
    1122                 :        119 :   nf = checknf(nf);
    1123                 :        119 :   checkgal(gal);
    1124                 :        119 :   checkprid(pr);
    1125                 :        119 :   gal_check_pol("idealramgroups",nf_get_pol(nf),gal_get_pol(gal));
    1126                 :        119 :   e = pr_get_e(pr); n = nf_get_degree(nf);
    1127                 :        119 :   p = itos(pr_get_p(pr));
    1128         [ +  + ]:        119 :   if (e%p) return idealramgroupstame(nf, gal, pr);
    1129                 :         14 :   (void) u_lvalrem(e,p,&et);
    1130                 :         14 :   idx = idealramgroupindex(nf, gal, pr);
    1131                 :         14 :   sub = group_subgroups(galois_group(gal));
    1132                 :         14 :   tbl = subgroups_tableset(sub, n);
    1133                 :         14 :   maxm = vecsmall_max(idx)+1;
    1134                 :         14 :   res = cgetg(maxm+1,t_VEC);
    1135                 :         14 :   set = zero_F2v(n); F2v_set(set,1);
    1136         [ +  + ]:         77 :   for(i=maxm; i>0; i--)
    1137                 :            :   {
    1138         [ +  + ]:       1183 :     for(j=1;j<=n;j++)
    1139         [ +  + ]:       1120 :       if (idx[j]==i-1)
    1140                 :         56 :         F2v_set(set,j);
    1141                 :         63 :     gel(res,i) = gel(sub, tableset_find_index(tbl, set));
    1142                 :            :   }
    1143                 :        119 :   return gerepilecopy(av, res);
    1144                 :            : }
    1145                 :            : 
    1146                 :            : /* x = relative polynomial nf = absolute nf, bnf = absolute bnf */
    1147                 :            : GEN
    1148                 :        112 : get_bnfpol(GEN x, GEN *bnf, GEN *nf)
    1149                 :            : {
    1150                 :        112 :   *bnf = checkbnf_i(x);
    1151                 :        112 :   *nf  = checknf_i(x);
    1152         [ +  + ]:        112 :   if (*nf) x = nf_get_pol(*nf);
    1153         [ -  + ]:        112 :   if (typ(x) != t_POL) pari_err_TYPE("get_bnfpol",x);
    1154                 :        112 :   return x;
    1155                 :            : }
    1156                 :            : 
    1157                 :            : GEN
    1158                 :      14725 : get_nfpol(GEN x, GEN *nf)
    1159                 :            : {
    1160         [ +  + ]:      14725 :   if (typ(x) == t_POL) { *nf = NULL; return x; }
    1161                 :      14725 :   *nf = checknf(x); return nf_get_pol(*nf);
    1162                 :            : }
    1163                 :            : 
    1164                 :            : /* is isomorphism / inclusion (a \subset b) compatible with what we know about
    1165                 :            :  * basic invariants ? (degree, signature, discriminant) */
    1166                 :            : static int
    1167                 :         49 : tests_OK(GEN a, GEN nfa, GEN b, GEN nfb, long fliso)
    1168                 :            : {
    1169                 :            :   GEN da, db, fa, P, E, U;
    1170                 :         49 :   long i, nP, m = degpol(a), n = degpol(b), q = m / n; /* relative degree */
    1171                 :            : 
    1172         [ -  + ]:         49 :   if (m <= 0) pari_err_IRREDPOL("nfisincl",a);
    1173         [ -  + ]:         49 :   if (n <= 0) pari_err_IRREDPOL("nfisincl",b);
    1174 [ +  + ][ -  + ]:         49 :   if (fliso) { if (n != m) return 0; } else { if (n % m) return 0; }
                 [ -  + ]
    1175         [ +  + ]:         49 :   if (m == 1) return 1;
    1176                 :            : 
    1177 [ -  + ][ #  # ]:         42 :   if (nfa && nfb) /* both nf structures available */
    1178                 :            :   {
    1179                 :          0 :     long r1a = nf_get_r1(nfa), r1b = nf_get_r1(nfb) ;
    1180         [ #  # ]:          0 :     if (fliso)
    1181 [ #  # ][ #  # ]:          0 :       return (r1a == r1b && equalii(nf_get_disc(nfa), nf_get_disc(nfb)));
    1182                 :            :     else
    1183   [ #  #  #  # ]:          0 :       return (r1b <= r1a * q &&
    1184                 :          0 :               dvdii(nf_get_disc(nfb), powiu(nf_get_disc(nfa), q)));
    1185                 :            :   }
    1186         [ -  + ]:         42 :   da = nfa? nf_get_disc(nfa): ZX_disc(a);
    1187         [ +  + ]:         42 :   if (!signe(da)) pari_err_IRREDPOL("nfisincl",a);
    1188         [ +  + ]:         35 :   db = nfb? nf_get_disc(nfb): ZX_disc(b);
    1189         [ -  + ]:         35 :   if (!signe(db)) pari_err_IRREDPOL("nfisincl",a);
    1190         [ +  + ]:         35 :   if (fliso) return issquare(gdiv(da,db));
    1191                 :            : 
    1192 [ +  + ][ -  + ]:         21 :   if (odd(q) && signe(da) != signe(db)) return 0;
    1193                 :         21 :   fa = absi_factor_limit(da, 0);
    1194                 :         21 :   P = gel(fa,1);
    1195                 :         21 :   E = gel(fa,2); nP = lg(P) - 1;
    1196         [ +  + ]:         77 :   for (i=1; i<nP; i++)
    1197 [ +  + ][ -  + ]:         56 :     if (mod2(gel(E,i)) && !dvdii(db, powiu(gel(P,i),q))) return 0;
    1198                 :         21 :   U = gel(P,nP);
    1199 [ -  + ][ #  # ]:         21 :   if (mod2(gel(E,i)) && expi(U) < 150)
    1200                 :            :   { /* "unfactored" cofactor is small, finish */
    1201         [ #  # ]:          0 :     if (cmpiu(U, maxprime()) > 0)
    1202                 :            :     {
    1203                 :          0 :       fa = Z_factor(U);
    1204                 :          0 :       P = gel(fa,1);
    1205                 :          0 :       E = gel(fa,2);
    1206                 :            :     }
    1207                 :            :     else
    1208                 :            :     {
    1209                 :          0 :       P = mkvec(U);
    1210                 :          0 :       E = mkvec(gen_1);
    1211                 :            :     }
    1212                 :          0 :     nP = lg(P) - 1;
    1213         [ #  # ]:          0 :     for (i=1; i<=nP; i++)
    1214 [ #  # ][ #  # ]:          0 :       if (mod2(gel(E,i)) && !dvdii(db, powiu(gel(P,i),q))) return 0;
    1215                 :            :   }
    1216                 :         42 :   return 1;
    1217                 :            : }
    1218                 :            : 
    1219                 :            : /* if fliso test for isomorphism, for inclusion otherwise. */
    1220                 :            : static GEN
    1221                 :         49 : nfiso0(GEN a, GEN b, long fliso)
    1222                 :            : {
    1223                 :         49 :   pari_sp av = avma;
    1224                 :            :   long i, vb, lx;
    1225                 :            :   GEN nfa, nfb, y, la, lb;
    1226                 :            :   int newvar;
    1227                 :            : 
    1228                 :         49 :   a = get_nfpol(a, &nfa);
    1229                 :         49 :   b = get_nfpol(b, &nfb);
    1230         [ +  + ]:         49 :   if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nsiso0"); }
    1231         [ +  + ]:         49 :   if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nsiso0"); }
    1232 [ +  + ][ -  + ]:         49 :   if (fliso && nfa && !nfb) { swap(a,b); nfb = nfa; nfa = NULL; }
                 [ #  # ]
    1233         [ -  + ]:         49 :   if (!tests_OK(a, nfa, b, nfb, fliso)) { avma = av; return gen_0; }
    1234                 :            : 
    1235         [ +  + ]:         42 :   if (nfb) lb = gen_1; else b = ZX_Q_normalize(b,&lb);
    1236         [ +  + ]:         42 :   if (nfa) la = gen_1; else a = ZX_Q_normalize(a,&la);
    1237                 :         42 :   vb = varn(b); newvar = (varncmp(vb,varn(a)) <= 0);
    1238         [ +  + ]:         42 :   if (newvar) { a = leafcopy(a); setvarn(a, fetch_var_higher()); }
    1239         [ +  + ]:         42 :   if (nfb)
    1240                 :         14 :     y = lift_intern(nfroots(nfb,a));
    1241                 :            :   else
    1242                 :            :   {
    1243                 :         28 :     y = gel(polfnf(a,b),1); lx = lg(y);
    1244         [ +  + ]:        154 :     for (i=1; i<lx; i++)
    1245                 :            :     {
    1246                 :        126 :       GEN t = gel(y,i);
    1247         [ -  + ]:        126 :       if (degpol(t) != 1) { setlg(y,i); break; }
    1248                 :        126 :       gel(y,i) = gneg_i(lift_intern(gel(t,2)));
    1249                 :            :     }
    1250                 :         28 :     settyp(y, t_VEC);
    1251                 :         28 :     gen_sort_inplace(y, (void*)&cmp_RgX, &cmp_nodata, NULL);
    1252                 :            :   }
    1253         [ +  + ]:         42 :   if (newvar) (void)delete_var();
    1254         [ -  + ]:         42 :   lx = lg(y); if (lx==1) { avma=av; return gen_0; }
    1255         [ +  + ]:        189 :   for (i=1; i<lx; i++)
    1256                 :            :   {
    1257                 :        147 :     GEN t = gel(y,i);
    1258         [ +  + ]:        147 :     if (typ(t) == t_POL) setvarn(t, vb); else t = scalarpol(t, vb);
    1259         [ +  + ]:        147 :     if (lb != gen_1) t = RgX_unscale(t, lb);
    1260         [ +  + ]:        147 :     if (la != gen_1) t = RgX_Rg_div(t, la);
    1261                 :        147 :     gel(y,i) = t;
    1262                 :            :   }
    1263                 :         42 :   return gerepilecopy(av,y);
    1264                 :            : }
    1265                 :            : 
    1266                 :            : GEN
    1267                 :         14 : nfisisom(GEN a, GEN b) { return nfiso0(a,b,1); }
    1268                 :            : 
    1269                 :            : GEN
    1270                 :         35 : nfisincl(GEN a, GEN b) { return nfiso0(a,b,0); }
    1271                 :            : 
    1272                 :            : /*************************************************************************/
    1273                 :            : /**                                                                     **/
    1274                 :            : /**                               INITALG                               **/
    1275                 :            : /**                                                                     **/
    1276                 :            : /*************************************************************************/
    1277                 :            : 
    1278                 :            : GEN
    1279                 :       9043 : get_roots(GEN x, long r1, long prec)
    1280                 :            : {
    1281                 :            :   long i, ru;
    1282                 :            :   GEN z;
    1283         [ -  + ]:       9043 :   if (typ(x) != t_POL)
    1284                 :            :   {
    1285                 :          0 :     z = leafcopy(x);
    1286                 :          0 :     ru = (lg(z)-1 + r1) >> 1;
    1287                 :            :   }
    1288                 :            :   else
    1289                 :            :   {
    1290                 :       9043 :     long n = degpol(x);
    1291         [ +  + ]:       9043 :     z = (r1 == n)? realroots(x, NULL, prec): QX_complex_roots(x,prec);
    1292                 :       9043 :     ru = (n+r1)>>1;
    1293                 :            :   }
    1294         [ +  + ]:      22399 :   for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
    1295                 :       9043 :   z[0]=evaltyp(t_VEC)|evallg(ru+1); return z;
    1296                 :            : }
    1297                 :            : 
    1298                 :            : GEN
    1299                 :          0 : nf_get_allroots(GEN nf)
    1300                 :            : {
    1301                 :          0 :   return embed_roots(nf_get_roots(nf), nf_get_r1(nf));
    1302                 :            : }
    1303                 :            : 
    1304                 :            : /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
    1305                 :            : GEN
    1306                 :      58891 : quicktrace(GEN x, GEN sym)
    1307                 :            : {
    1308                 :      58891 :   GEN p1 = gen_0;
    1309                 :            :   long i;
    1310                 :            : 
    1311         [ -  + ]:      58891 :   if (typ(x) != t_POL) return gmul(x, gel(sym,1));
    1312         [ +  - ]:      58891 :   if (signe(x))
    1313                 :            :   {
    1314                 :      58891 :     sym--;
    1315         [ +  + ]:     834043 :     for (i=lg(x)-1; i>1; i--)
    1316                 :     775152 :       p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
    1317                 :            :   }
    1318                 :      58891 :   return p1;
    1319                 :            : }
    1320                 :            : 
    1321                 :            : static GEN
    1322                 :       5061 : get_Tr(GEN mul, GEN x, GEN basden)
    1323                 :            : {
    1324                 :       5061 :   GEN t, bas = gel(basden,1), den = gel(basden,2);
    1325                 :       5061 :   long i, j, n = lg(bas)-1;
    1326                 :       5061 :   GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
    1327                 :            : 
    1328                 :       5061 :   gel(TW,1) = utoipos(n);
    1329         [ +  + ]:      19313 :   for (i=2; i<=n; i++)
    1330                 :            :   {
    1331                 :      14252 :     t = quicktrace(gel(bas,i), sym);
    1332 [ +  + ][ +  + ]:      14252 :     if (den && den[i]) t = diviiexact(t,gel(den,i));
    1333                 :      14252 :     gel(TW,i) = t; /* tr(w[i]) */
    1334                 :            :   }
    1335                 :       5061 :   gel(T,1) = TW;
    1336         [ +  + ]:      19313 :   for (i=2; i<=n; i++)
    1337                 :            :   {
    1338                 :      14252 :     gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
    1339         [ +  + ]:     117551 :     for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
    1340                 :     103299 :       gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
    1341                 :            :   }
    1342                 :       5061 :   return T;
    1343                 :            : }
    1344                 :            : 
    1345                 :            : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
    1346                 :            : GEN
    1347                 :      12988 : get_bas_den(GEN bas)
    1348                 :            : {
    1349                 :      12988 :   GEN b,d,den, dbas = leafcopy(bas);
    1350                 :      12988 :   long i, l = lg(bas);
    1351                 :      12988 :   int power = 1;
    1352                 :      12988 :   den = cgetg(l,t_VEC);
    1353         [ +  + ]:      67054 :   for (i=1; i<l; i++)
    1354                 :            :   {
    1355                 :      54066 :     b = Q_remove_denom(gel(bas,i), &d);
    1356                 :      54066 :     gel(dbas,i) = b;
    1357         [ +  + ]:      54066 :     gel(den,i) = d; if (d) power = 0;
    1358                 :            :   }
    1359         [ +  + ]:      12988 :   if (power) den = NULL; /* power basis */
    1360                 :      12988 :   return mkvec2(dbas, den);
    1361                 :            : }
    1362                 :            : 
    1363                 :            : /* Internal: nf partially filled. Require pol; fill zk, invzk, multable */
    1364                 :            : void
    1365                 :       5061 : nf_set_multable(GEN nf, GEN bas, GEN basden)
    1366                 :            : {
    1367                 :       5061 :   GEN T = nf_get_pol(nf), invbas, basM;
    1368                 :       5061 :   long i,j, n = degpol(T);
    1369                 :       5061 :   GEN w, den, mul = cgetg(n*n+1,t_MAT);
    1370                 :            : 
    1371         [ -  + ]:       5061 :   if (typ(bas) == t_MAT)
    1372                 :          0 :   { basM = bas; bas = RgM_to_RgXV(basM, varn(T)); }
    1373                 :            :   else
    1374                 :       5061 :     basM = RgV_to_RgM(bas, n);
    1375                 :       5061 :   gel(nf,7) = bas;
    1376                 :       5061 :   gel(nf,8) = invbas = QM_inv(basM, gen_1);
    1377                 :       5061 :   gel(nf,9) = mul;
    1378                 :            : 
    1379         [ -  + ]:       5061 :   if (!basden) basden = get_bas_den(nf_get_zk(nf)); /*integral basis*/
    1380                 :       5061 :   w   = gel(basden,1);
    1381                 :       5061 :   den = gel(basden,2);
    1382                 :            :   /* i = 1 split for efficiency, assume w[1] = 1 */
    1383         [ +  + ]:      24374 :   for (j=1; j<=n; j++)
    1384                 :      19313 :     gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
    1385         [ +  + ]:      19313 :   for (i=2; i<=n; i++)
    1386         [ +  + ]:     117551 :     for (j=i; j<=n; j++)
    1387                 :            :     {
    1388                 :     103299 :       pari_sp av = avma;
    1389         [ +  + ]:     103299 :       GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
    1390                 :     103299 :       z = mulmat_pol(invbas, z); /* integral column */
    1391         [ +  + ]:     103299 :       if (den)
    1392                 :            :       {
    1393                 :      71127 :         GEN d = mul_denom(gel(den,i), gel(den,j));
    1394         [ +  + ]:      71127 :         if (d) z = ZC_Z_divexact(z, d);
    1395                 :            :       }
    1396                 :     103299 :       gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
    1397                 :            :     }
    1398                 :       5061 : }
    1399                 :            : 
    1400                 :            : /* as get_Tr, mul_table not precomputed */
    1401                 :            : static GEN
    1402                 :       2709 : make_Tr(GEN x, GEN basden)
    1403                 :            : {
    1404                 :       2709 :   long i,j, n = degpol(x);
    1405                 :            :   GEN c, t, d;
    1406                 :       2709 :   GEN T   = cgetg(n+1,t_MAT);
    1407                 :       2709 :   GEN sym = polsym(x, n-1);
    1408                 :       2709 :   GEN w   = gel(basden,1); /* W[i] = w[i]/den[i] */
    1409                 :       2709 :   GEN den = gel(basden,2);
    1410                 :            :   /* assume W[1] = 1, case i = 1 split for efficiency */
    1411                 :       2709 :   c = cgetg(n+1,t_COL); gel(T,1) = c;
    1412                 :       2709 :   gel(c, 1) = utoipos(n);
    1413         [ +  + ]:       8169 :   for (j=2; j<=n; j++)
    1414                 :            :   {
    1415                 :       5460 :     pari_sp av = avma;
    1416                 :       5460 :     t = quicktrace(gel(w,j), sym);
    1417         [ +  + ]:       5460 :     if (den)
    1418                 :            :     {
    1419                 :       3962 :       d = gel(den,j);
    1420         [ +  + ]:       3962 :       if (d) t = diviiexact(t, d);
    1421                 :            :     }
    1422                 :       5460 :     gel(c,j) = gerepileuptoint(av, t);
    1423                 :            :   }
    1424         [ +  + ]:       8169 :   for (i=2; i<=n; i++)
    1425                 :            :   {
    1426                 :       5460 :     c = cgetg(n+1,t_COL); gel(T,i) = c;
    1427         [ +  + ]:      36498 :     for (j=1; j<i ; j++) gel(c,j) = gcoeff(T,i,j);
    1428         [ +  + ]:      36498 :     for (   ; j<=n; j++)
    1429                 :            :     {
    1430                 :      31038 :       pari_sp av = avma;
    1431         [ +  + ]:      31038 :       t = (i == j)? ZXQ_sqr(gel(w,i), x): ZXQ_mul(gel(w,i),gel(w,j), x);
    1432                 :      31038 :       t = quicktrace(t, sym);
    1433         [ +  + ]:      31038 :       if (den)
    1434                 :            :       {
    1435                 :      28441 :         d = mul_denom(gel(den,i),gel(den,j));
    1436         [ +  + ]:      28441 :         if (d) t = diviiexact(t, d);
    1437                 :            :       }
    1438                 :      31038 :       gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
    1439                 :            :     }
    1440                 :            :   }
    1441                 :       2709 :   return T;
    1442                 :            : }
    1443                 :            : 
    1444                 :            : /* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */
    1445                 :            : static void
    1446                 :      10320 : get_roots_for_M(nffp_t *F)
    1447                 :            : {
    1448                 :            :   long n, eBD, PREC;
    1449                 :            : 
    1450         [ +  - ]:      10320 :   if (F->extraprec < 0)
    1451                 :            :   { /* not initialized yet */
    1452                 :            :     double er;
    1453                 :      10320 :     n = degpol(F->x);
    1454                 :      10320 :     eBD = 1 + gexpo(gel(F->basden,1));
    1455         [ +  + ]:      10320 :     er  = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->x);
    1456         [ +  + ]:      10320 :     if (er < 0) er = 0;
    1457                 :      10320 :     F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
    1458                 :            :   }
    1459                 :            : 
    1460                 :      10320 :   PREC = F->prec + F->extraprec;
    1461                 :            :   /* make sure that default accuracy is the same on 32/64bit */
    1462                 :            : #ifndef LONG_IS_64BIT
    1463         [ +  + ]:       1512 :   if (odd(PREC)) PREC += EXTRAPRECWORD;
    1464                 :            : #endif
    1465 [ +  + ][ +  + ]:      19363 :   if (F->ro && gprecision(gel(F->ro,1)) >= PREC) return;
    1466                 :       9043 :   F->ro = get_roots(F->x, F->r1, PREC);
    1467                 :            : }
    1468                 :            : 
    1469                 :            : /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
    1470                 :            : static void
    1471                 :      10320 : make_M(nffp_t *F, int trunc)
    1472                 :            : {
    1473                 :      10320 :   GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
    1474                 :            :   GEN m, d, M;
    1475                 :      10320 :   long i, j, l = lg(ro), n = lg(bas);
    1476                 :      10320 :   M = cgetg(n,t_MAT);
    1477                 :      10320 :   gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
    1478         [ +  + ]:      46007 :   for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
    1479         [ +  + ]:      39302 :   for (i=1; i<l; i++)
    1480                 :            :   {
    1481                 :      28982 :     GEN r = gel(ro,i), ri;
    1482         [ +  + ]:      28982 :     ri = (gexpo(r) > 1)? ginv(r): NULL;
    1483         [ +  + ]:     306525 :     for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
    1484                 :            :   }
    1485         [ +  + ]:      10320 :   if (den)
    1486         [ +  + ]:      28233 :     for (j=2; j<n; j++)
    1487                 :            :     {
    1488         [ +  + ]:      23870 :       d = gel(den,j); if (!d) continue;
    1489                 :      19403 :       m = gel(M,j);
    1490         [ +  + ]:     190501 :       for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
    1491                 :            :     }
    1492                 :            : 
    1493 [ +  + ][ +  + ]:      10320 :   if (trunc && gprecision(M) > F->prec)
    1494                 :            :   {
    1495                 :       1857 :     M     = gprec_w(M, F->prec);
    1496                 :       1857 :     F->ro = gprec_w(ro,F->prec);
    1497                 :            :   }
    1498                 :      10320 :   F->M = M;
    1499                 :      10320 : }
    1500                 :            : 
    1501                 :            : /* return G real such that G~ * G = T_2 */
    1502                 :            : static void
    1503                 :      10320 : make_G(nffp_t *F)
    1504                 :            : {
    1505                 :      10320 :   GEN G, M = F->M;
    1506                 :      10320 :   long i, j, k, r1 = F->r1, l = lg(M);
    1507                 :            : 
    1508                 :      10320 :   G = cgetg(l, t_MAT);
    1509         [ +  + ]:      56327 :   for (j=1; j<l; j++)
    1510                 :            :   {
    1511                 :      46007 :     GEN g = cgetg(l, t_COL);
    1512                 :      46007 :     GEN m = gel(M,j);
    1513                 :      46007 :     gel(G,j) = g;
    1514         [ +  + ]:     132606 :     for (k=i=1; i<=r1; i++) g[k++] = m[i];
    1515         [ +  + ]:     265933 :     for (     ; k < l; i++)
    1516                 :            :     {
    1517                 :     219926 :       GEN r = gel(m,i);
    1518         [ +  + ]:     219926 :       if (typ(r) == t_COMPLEX)
    1519                 :            :       {
    1520                 :     202901 :         gel(g,k++) = mpadd(gel(r,1), gel(r,2));
    1521                 :     202901 :         gel(g,k++) = mpsub(gel(r,1), gel(r,2));
    1522                 :            :       }
    1523                 :            :       else
    1524                 :            :       {
    1525                 :      17025 :         gel(g,k++) = r;
    1526                 :      17025 :         gel(g,k++) = r;
    1527                 :            :       }
    1528                 :            :     }
    1529                 :            :   }
    1530                 :      10320 :   F->G = G;
    1531                 :      10320 : }
    1532                 :            : 
    1533                 :            : static void
    1534                 :      10320 : make_M_G(nffp_t *F, int trunc)
    1535                 :            : {
    1536                 :      10320 :   get_roots_for_M(F);
    1537                 :      10320 :   make_M(F, trunc);
    1538                 :      10320 :   make_G(F);
    1539                 :      10320 : }
    1540                 :            : 
    1541                 :            : void
    1542                 :        752 : remake_GM(GEN nf, nffp_t *F, long prec)
    1543                 :            : {
    1544                 :        752 :   F->x  = nf_get_pol(nf);
    1545                 :        752 :   F->ro = NULL;
    1546                 :        752 :   F->r1 = nf_get_r1(nf);
    1547                 :        752 :   F->basden = get_bas_den(nf_get_zk(nf));
    1548                 :        752 :   F->extraprec = -1;
    1549                 :        752 :   F->prec = prec; make_M_G(F, 1);
    1550                 :        752 : }
    1551                 :            : 
    1552                 :            : static void
    1553                 :       9568 : nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
    1554                 :            : {
    1555                 :       9568 :   F->x  = T->x;
    1556                 :       9568 :   F->ro = ro;
    1557                 :       9568 :   F->r1 = T->r1;
    1558         [ +  + ]:       9568 :   if (!T->basden) T->basden = get_bas_den(T->bas);
    1559                 :       9568 :   F->basden = T->basden;
    1560                 :       9568 :   F->extraprec = -1;
    1561                 :       9568 :   F->prec = prec;
    1562                 :       9568 : }
    1563                 :            : 
    1564                 :            : static void
    1565                 :       6194 : get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, int trunc, long prec)
    1566                 :            : {
    1567                 :       6194 :   nffp_init(F,T,ro,prec);
    1568                 :       6194 :   make_M_G(F, trunc);
    1569                 :       6194 : }
    1570                 :            : 
    1571                 :            : static GEN
    1572                 :       5061 : get_sign(long r1, long n) { return mkvec2s(r1, (n-r1)>>1); }
    1573                 :            : 
    1574                 :            : GEN
    1575                 :       5061 : nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec)
    1576                 :            : {
    1577                 :       5061 :   GEN nf = cgetg(10,t_VEC);
    1578                 :       5061 :   GEN x = T->x, absdK, Tr, D, TI, A, dA, MDI, mat = cgetg(9,t_VEC);
    1579                 :       5061 :   long n = degpol(T->x);
    1580                 :            :   nffp_t F;
    1581                 :       5061 :   get_nf_fp_compo(T, &F, ro, 0, prec);
    1582                 :            : 
    1583                 :       5061 :   gel(nf,1) = T->x;
    1584                 :       5061 :   gel(nf,2) = get_sign(T->r1, n);
    1585                 :       5061 :   gel(nf,3) = T->dK;
    1586                 :       5061 :   gel(nf,4) = T->index;
    1587                 :       5061 :   gel(nf,6) = F.ro;
    1588                 :       5061 :   gel(nf,5) = mat;
    1589                 :            : 
    1590                 :       5061 :   gel(mat,1) = F.M;
    1591                 :       5061 :   gel(mat,2) = F.G;
    1592                 :            : 
    1593                 :       5061 :   nf_set_multable(nf, T->bas, F.basden);
    1594                 :            : 
    1595                 :       5061 :   Tr = get_Tr(gel(nf,9), x, F.basden);
    1596         [ +  + ]:       5061 :   absdK = T->dK; if (signe(absdK) < 0) absdK = negi(absdK);
    1597                 :       5061 :   TI = ZM_inv(Tr, absdK); /* dK T^-1 */
    1598                 :       5061 :   A = Q_primitive_part(TI, &dA);
    1599                 :       5061 :   gel(mat,6) = A; /* primitive part of codifferent, dA its content */
    1600         [ +  + ]:       5061 :   dA = dA? diviiexact(absdK, dA): absdK;
    1601                 :       5061 :   A = ZM_hnfmodid(A, dA);
    1602                 :       5061 :   MDI = idealtwoelt(nf, A);
    1603                 :       5061 :   gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
    1604                 :       5061 :   gel(mat,7) = MDI;
    1605         [ +  + ]:       5061 :   if (is_pm1(T->index)) /* principal ideal (x'), whose norm is |dK| */
    1606                 :            :   {
    1607                 :       3514 :     D = zk_scalar_or_multable(nf, ZX_deriv(x));
    1608         [ +  + ]:       3514 :     if (typ(D) == t_MAT) D = ZM_hnfmod(D, absdK);
    1609                 :            :   }
    1610                 :            :   else
    1611                 :       1547 :     D = RgM_Rg_mul(idealinv(nf, A), dA);
    1612                 :       5061 :   gel(mat,3) = RM_round_maxrank(F.G);
    1613                 :       5061 :   gel(mat,4) = Tr;
    1614                 :       5061 :   gel(mat,5) = D;
    1615         [ +  + ]:       5061 :   gel(mat,8) = T->dKP? shallowtrans(T->dKP): cgetg(1,t_VEC);
    1616                 :       5061 :   return nf;
    1617                 :            : }
    1618                 :            : 
    1619                 :            : static GEN
    1620                 :         91 : primes_certify(GEN dK, GEN dKP)
    1621                 :            : {
    1622                 :         91 :   long i, l = lg(dKP);
    1623                 :         91 :   GEN v, w, D = dK;
    1624                 :         91 :   v = vectrunc_init(l);
    1625                 :         91 :   w = vectrunc_init(l);
    1626         [ +  + ]:        427 :   for (i = 1; i < l; i++)
    1627                 :            :   {
    1628                 :        336 :     GEN p = gel(dKP,i);
    1629         [ +  + ]:        336 :     vectrunc_append(isprime(p)? w: v, p);
    1630                 :        336 :     (void)Z_pvalrem(D, p, &D);
    1631                 :            :   }
    1632         [ -  + ]:         91 :   if (!is_pm1(D))
    1633                 :            :   {
    1634         [ #  # ]:          0 :     if (signe(D) < 0) D = negi(D);
    1635         [ #  # ]:          0 :     vectrunc_append(isprime(D)? w: v, D);
    1636                 :            :   }
    1637                 :         91 :   return mkvec2(v,w);
    1638                 :            : }
    1639                 :            : GEN
    1640                 :          7 : nfcertify(GEN nf)
    1641                 :            : {
    1642                 :          7 :   pari_sp av = avma;
    1643                 :            :   GEN vw;
    1644                 :          7 :   nf = checknf(nf);
    1645                 :          7 :   vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
    1646                 :          7 :   return gerepilecopy(av, gel(vw,1));
    1647                 :            : }
    1648                 :            : 
    1649                 :            : #if 0 /* used to check benches between HNF nf.zk and LLL-reduced nf.zk */
    1650                 :            : static GEN
    1651                 :            : hnffromLLL(GEN nf)
    1652                 :            : {
    1653                 :            :   GEN d, x;
    1654                 :            :   x = RgV_to_RgM(nf_get_zk(nf), nf_get_degree(nf));
    1655                 :            :   x = Q_remove_denom(x, &d);
    1656                 :            :   if (!d) return x; /* power basis */
    1657                 :            :   return RgM_solve(ZM_hnfmodid(x, d), x);
    1658                 :            : }
    1659                 :            : 
    1660                 :            : static GEN
    1661                 :            : nfbasechange(GEN u, GEN x)
    1662                 :            : {
    1663                 :            :   long i,lx;
    1664                 :            :   GEN y;
    1665                 :            :   switch(typ(x))
    1666                 :            :   {
    1667                 :            :     case t_COL: /* nfelt */
    1668                 :            :       return RgM_RgC_mul(u, x);
    1669                 :            : 
    1670                 :            :     case t_MAT: /* ideal */
    1671                 :            :       y = cgetg_copy(x, &lx);
    1672                 :            :       for (i=1; i<lx; i++) gel(y,i) = RgM_RgC_mul(u, gel(x,i));
    1673                 :            :       break;
    1674                 :            : 
    1675                 :            :     case t_VEC: /* pr */
    1676                 :            :       checkprid(x); y = leafcopy(x);
    1677                 :            :       gel(y,2) = RgM_RgC_mul(u, gel(y,2));
    1678                 :            :       gel(y,5) = RgM_RgC_mul(u, gel(y,5));
    1679                 :            :       break;
    1680                 :            :     default: y = x;
    1681                 :            :   }
    1682                 :            :   return y;
    1683                 :            : }
    1684                 :            : 
    1685                 :            : GEN
    1686                 :            : nffromhnfbasis(GEN nf, GEN x)
    1687                 :            : {
    1688                 :            :   long tx = typ(x);
    1689                 :            :   pari_sp av = avma;
    1690                 :            :   GEN u;
    1691                 :            :   if (!is_vec_t(tx)) return gcopy(x);
    1692                 :            :   nf = checknf(nf);
    1693                 :            :   u = hnffromLLL(nf);
    1694                 :            :   return gerepilecopy(av, nfbasechange(u, x));
    1695                 :            : }
    1696                 :            : 
    1697                 :            : GEN
    1698                 :            : nftohnfbasis(GEN nf, GEN x)
    1699                 :            : {
    1700                 :            :   long tx = typ(x);
    1701                 :            :   pari_sp av = avma;
    1702                 :            :   GEN u;
    1703                 :            :   if (!is_vec_t(tx)) return gcopy(x);
    1704                 :            :   nf = checknf(nf);
    1705                 :            :   u = ZM_inv(hnffromLLL(nf), gen_1);
    1706                 :            :   return gerepilecopy(av, nfbasechange(u, x));
    1707                 :            : }
    1708                 :            : #endif
    1709                 :            : 
    1710                 :            : /* set *pro to roots of T->x */
    1711                 :            : static GEN
    1712                 :       3374 : get_red_G(nfbasic_t *T, GEN *pro)
    1713                 :            : {
    1714                 :       3374 :   GEN G, u, u0 = NULL;
    1715                 :            :   pari_sp av;
    1716                 :       3374 :   long i, prec, n = degpol(T->x);
    1717                 :            :   nffp_t F;
    1718                 :            : 
    1719                 :       3374 :   prec = nbits2prec(n+32);
    1720                 :       3374 :   nffp_init(&F, T, NULL, prec);
    1721                 :       3374 :   av = avma;
    1722                 :       3374 :   for (i=1; ; i++)
    1723                 :            :   {
    1724                 :       3374 :     F.prec = prec; make_M_G(&F, 0); G = F.G;
    1725         [ -  + ]:       3374 :     if (u0) G = RgM_mul(G, u0);
    1726         [ -  + ]:       3374 :     if (DEBUGLEVEL)
    1727                 :          0 :       err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
    1728                 :          0 :                   prec + F.extraprec, prec, F.extraprec);
    1729         [ +  - ]:       3374 :     if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST|LLL_COMPATIBLE)))
    1730                 :            :     {
    1731         [ +  - ]:       3374 :       if (lg(u)-1 == n) break;
    1732                 :            :       /* singular ==> loss of accuracy */
    1733         [ #  # ]:          0 :       if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
    1734                 :          0 :       else    u0 = gerepilecopy(av, u);
    1735                 :            :     }
    1736                 :          0 :     prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
    1737                 :          0 :     F.ro = NULL;
    1738         [ #  # ]:          0 :     if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
    1739                 :          0 :   }
    1740         [ -  + ]:       3374 :   if (u0) u = RgM_mul(u0,u);
    1741                 :       3374 :   *pro = F.ro; return u;
    1742                 :            : }
    1743                 :            : 
    1744                 :            : /* Compute an LLL-reduced basis for the integer basis of nf(T).
    1745                 :            :  * set *pro = roots of x if computed [NULL if not computed] */
    1746                 :            : static void
    1747                 :       6083 : set_LLL_basis(nfbasic_t *T, GEN *pro, double DELTA)
    1748                 :            : {
    1749                 :       6083 :   GEN B = T->bas;
    1750         [ +  + ]:       6083 :   if (T->r1 == degpol(T->x)) {
    1751                 :       2709 :     pari_sp av = avma;
    1752                 :       2709 :     GEN u, basden = T->basden;
    1753         [ +  - ]:       2709 :     if (!basden) basden = get_bas_den(B);
    1754                 :       2709 :     u = ZM_lll(make_Tr(T->x,basden), DELTA,
    1755                 :            :                LLL_GRAM|LLL_KEEP_FIRST|LLL_IM|LLL_COMPATIBLE);
    1756                 :       2709 :     B = gerepileupto(av, RgV_RgM_mul(B, u));
    1757                 :       2709 :     *pro = NULL;
    1758                 :            :   }
    1759                 :            :   else
    1760                 :       3374 :     B = RgV_RgM_mul(B, get_red_G(T, pro));
    1761                 :       6083 :   T->bas = B;
    1762                 :       6083 :   T->basden = get_bas_den(B);
    1763                 :       6083 : }
    1764                 :            : 
    1765                 :            : static int
    1766                 :       3172 : cmp_abs_ZX(GEN x, GEN y) { return gen_cmp_RgX((void*)&absi_cmp, x, y); }
    1767                 :            : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
    1768                 :            :  * (if so update *dx) */
    1769                 :            : static int
    1770                 :       4621 : ZX_is_better(GEN y, GEN x, GEN *dx)
    1771                 :            : {
    1772                 :       4621 :   GEN d = ZX_disc(y);
    1773                 :            :   int cmp;
    1774         [ +  + ]:       4621 :   if (!*dx) *dx = ZX_disc(x);
    1775                 :       4621 :   cmp = absi_cmp(d, *dx);
    1776         [ +  + ]:       4621 :   if (cmp < 0) { *dx = d; return 1; }
    1777         [ +  + ]:       3781 :   if (cmp == 0) return cmp_abs_ZX(y, x) < 0;
    1778                 :       4621 :   return 0;
    1779                 :            : }
    1780                 :            : 
    1781                 :            : static void polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
    1782                 :            : /* Seek a simpler, polynomial pol defining the same number field as
    1783                 :            :  * x (assumed to be monic at this point) */
    1784                 :            : static GEN
    1785                 :         84 : nfpolred(nfbasic_t *T, GEN *pro)
    1786                 :            : {
    1787                 :         84 :   GEN x = T->x, dx, b, rev;
    1788                 :         84 :   long n = degpol(x), v = varn(x);
    1789                 :            : 
    1790         [ +  + ]:         84 :   if (n == 1) {
    1791                 :          7 :     T->x = deg1pol_shallow(gen_1, gen_m1, v);
    1792                 :          7 :     *pro = NULL; return pol_1(v);
    1793                 :            :   }
    1794                 :         77 :   polredbest_aux(T, pro, &x, &dx, &b);
    1795         [ +  + ]:         77 :   if (x == T->x) return NULL; /* no improvement */
    1796         [ -  + ]:         56 :   if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
    1797                 :            : 
    1798                 :            :   /* update T */
    1799                 :         56 :   rev = QXQ_reverse(b, T->x);
    1800                 :         56 :   T->bas = QXV_QXQ_eval(T->bas, rev, x);
    1801                 :         56 :   (void)Z_issquareall(diviiexact(dx,T->dK), &(T->index));
    1802                 :         56 :   T->basden = get_bas_den(T->bas);
    1803                 :         56 :   T->dx = dx;
    1804                 :         56 :   T->x = x;
    1805                 :         56 :   *pro = NULL; /* reset */
    1806                 :         84 :   return rev;
    1807                 :            : }
    1808                 :            : 
    1809                 :            : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
    1810                 :            :  * basis. Assume bas[1] is 1 and that the leading coefficient of elements
    1811                 :            :  * of bas are of the form 1/b for a t_INT b */
    1812                 :            : GEN
    1813                 :        525 : get_nfindex(GEN bas)
    1814                 :            : {
    1815                 :        525 :   pari_sp av = avma;
    1816                 :        525 :   long n = lg(bas)-1, i;
    1817                 :            :   GEN D, d, mat;
    1818                 :            : 
    1819                 :        525 :   D = gen_1; /* assume bas[1] = 1 */
    1820         [ +  + ]:       1197 :   for (i = 2; i <= n; i++)
    1821                 :            :   { /* in most cases [e.g after nfbasis] basis is upper triangular! */
    1822                 :        924 :     GEN B = gel(bas,i), lc;
    1823         [ +  + ]:        924 :     if (degpol(B) != i-1) break;
    1824                 :        672 :     lc = gel(B, i+1);
    1825      [ +  +  - ]:        672 :     switch (typ(lc))
    1826                 :            :     {
    1827                 :        399 :       case t_INT: continue;
    1828                 :        273 :       case t_FRAC: lc = gel(lc,2); break;
    1829                 :          0 :       default: pari_err_TYPE("get_nfindex",lc);
    1830                 :            :     }
    1831                 :        273 :     D = mulii(D, lc);
    1832                 :            :   }
    1833         [ +  + ]:        525 :   if (i <= n)
    1834                 :            :   { /* not triangular after all */
    1835                 :        252 :     bas = Q_remove_denom(bas, &d);
    1836         [ +  + ]:        252 :     if (!d) { avma = av; return D; }
    1837                 :        238 :     mat = RgV_to_RgM(bas, n);
    1838                 :        238 :     d = diviiexact(powiu(d, n), ZM_det(mat));
    1839                 :        238 :     D = mulii(D,absi(d));
    1840                 :            :   }
    1841                 :        525 :   return gerepileuptoint(av, D);
    1842                 :            : }
    1843                 :            : 
    1844                 :            : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
    1845                 :            :  * basis (deprecated), or listP data (nfbasis input format) to specify
    1846                 :            :  * a set of primes at with the basis order must be maximal.
    1847                 :            :  * 1) nf type (or unrecognized): return t_VEC
    1848                 :            :  * 2) ZX or [ZX, listP]: return t_POL
    1849                 :            :  * 3) [ZX, order basis]: return 0 (deprecated)
    1850                 :            :  * incorrect: return -1 */
    1851                 :            : static long
    1852                 :       6167 : nf_input_type(GEN x)
    1853                 :            : {
    1854                 :            :   GEN T, V;
    1855                 :            :   long i, d, v;
    1856      [ +  +  - ]:       6167 :   switch(typ(x))
    1857                 :            :   {
    1858                 :       5530 :     case t_POL: return t_POL;
    1859                 :            :     case t_VEC:
    1860         [ +  + ]:        637 :       if (lg(x) != 3) return t_VEC; /* nf or incorrect */
    1861                 :        630 :       T = gel(x,1); V = gel(x,2);
    1862         [ -  + ]:        630 :       if (typ(T) != t_POL) return -1;
    1863      [ +  +  - ]:        630 :       switch(typ(V))
    1864                 :            :       {
    1865                 :         35 :         case t_INT: case t_MAT: return t_POL;
    1866                 :            :         case t_VEC: case t_COL:
    1867         [ +  + ]:        595 :           if (RgV_is_ZV(V)) return t_POL;
    1868                 :        574 :           break;
    1869                 :          0 :         default: return -1;
    1870                 :            :       }
    1871                 :        574 :       d = degpol(T); v = varn(T);
    1872 [ +  - ][ +  - ]:        574 :       if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
         [ +  - ][ +  + ]
    1873         [ +  + ]:       3192 :       for (i = 1; i <= d; i++)
    1874                 :            :       { /* check integer basis */
    1875                 :       2639 :         GEN c = gel(V,i);
    1876      [ +  +  - ]:       2639 :         switch(typ(c))
    1877                 :            :         {
    1878                 :         28 :           case t_INT: break;
    1879 [ +  + ][ +  - ]:       2611 :           case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
                 [ +  + ]
    1880                 :            :           /* fall through */
    1881                 :         14 :           default: return -1;
    1882                 :            :         }
    1883                 :            :       }
    1884                 :        553 :       return 0;
    1885                 :            :   }
    1886                 :       6167 :   return t_VEC; /* nf or incorrect */
    1887                 :            : }
    1888                 :            : 
    1889                 :            : static void
    1890                 :       5054 : nfbasic_add_disc(nfbasic_t *T)
    1891                 :            : {
    1892         [ +  + ]:       5054 :   if (!T->index) T->index = get_nfindex(T->bas);
    1893         [ +  + ]:       5054 :   if (!T->dx) T->dx = ZX_disc(T->x);
    1894         [ +  + ]:       5054 :   if (!T->dK) T->dK = diviiexact(T->dx, sqri(T->index));
    1895                 :       5054 : }
    1896                 :            : 
    1897                 :            : static void
    1898                 :       6167 : nfbasic_init(GEN x, long flag, nfbasic_t *T)
    1899                 :            : {
    1900                 :       6167 :   GEN bas, dK, dx, index, unscale = gen_1;
    1901                 :            :   long r1;
    1902                 :            : 
    1903                 :       6167 :   T->dKP = NULL;
    1904   [ +  +  +  + ]:       6167 :   switch (nf_input_type(x))
    1905                 :            :   {
    1906                 :            :     case t_POL:
    1907                 :            :     {
    1908                 :            :       nfmaxord_t S;
    1909                 :       5586 :       nfmaxord(&S, x, flag);
    1910                 :       5579 :       x = S.T;
    1911                 :       5579 :       T->x0 = S.T0;
    1912                 :       5579 :       T->dKP = S.dKP;
    1913                 :       5579 :       dK = S.dK;
    1914                 :       5579 :       index = S.index;
    1915                 :       5579 :       bas = S.basis;
    1916                 :       5579 :       dx = S.dT;
    1917                 :       5579 :       unscale = S.unscale;
    1918                 :       5579 :       r1 = ZX_sturm(x);
    1919                 :            :       break;
    1920                 :            :     }
    1921                 :            :     case t_VEC:
    1922                 :            :     { /* nf, bnf, bnr */
    1923                 :          7 :       GEN nf = checknf(x);
    1924                 :          7 :       x     = nf_get_pol(nf);
    1925                 :          7 :       dK    = nf_get_disc(nf);
    1926                 :          7 :       index = nf_get_index(nf);
    1927                 :          7 :       bas   = nf_get_zk(nf);
    1928                 :          7 :       T->x0 = x;
    1929                 :          7 :       dx = NULL;
    1930                 :          7 :       r1 = nf_get_r1(nf);
    1931                 :          7 :       break;
    1932                 :            :     }
    1933                 :            :     case 0: /* monic integral polynomial + integer basis */
    1934                 :        553 :       bas = gel(x,2); x = gel(x,1);
    1935                 :        553 :       T->x0 = x;
    1936                 :        553 :       index = NULL;
    1937                 :        553 :       dx = NULL;
    1938                 :        553 :       dK = NULL;
    1939                 :        553 :       r1 = ZX_sturm(x);
    1940                 :        553 :       break;
    1941                 :            :     default: /* -1 */
    1942                 :         21 :       pari_err_TYPE("nfbasic_init", x);
    1943                 :       6139 :       return;
    1944                 :            :   }
    1945                 :       6139 :   T->x     = x;
    1946                 :       6139 :   T->unscale = unscale;
    1947                 :       6139 :   T->r1    = r1;
    1948                 :       6139 :   T->dx    = dx;
    1949                 :       6139 :   T->dK    = dK;
    1950                 :       6139 :   T->bas   = bas;
    1951                 :       6139 :   T->basden= NULL;
    1952                 :       6139 :   T->index = index;
    1953                 :            : }
    1954                 :            : 
    1955                 :            : void
    1956                 :       5075 : nfinit_step1(nfbasic_t *T, GEN x, long flag)
    1957                 :            : {
    1958                 :       5075 :   nfbasic_init(x, flag, T);
    1959         [ -  + ]:       5054 :   if (!ZX_is_irred(T->x)) pari_err_IRREDPOL("nfinit",x);
    1960                 :       5054 : }
    1961                 :            : 
    1962                 :            : GEN
    1963                 :       5054 : nfinit_step2(nfbasic_t *T, long flag, long prec)
    1964                 :            : {
    1965                 :            :   GEN nf, unscale;
    1966                 :            : 
    1967 [ +  + ][ +  + ]:       5054 :   if (!(flag & nf_RED) && !equali1(leading_coeff(T->x0)))
    1968                 :            :   {
    1969                 :         49 :     pari_warn(warner,"non-monic polynomial. Result of the form [nf,c]");
    1970                 :         49 :     flag |= nf_RED | nf_ORIG;
    1971                 :            :   }
    1972                 :       5054 :   unscale = T->unscale;
    1973 [ +  + ][ +  + ]:       5054 :   if (!(flag & nf_RED) && !isint1(unscale))
    1974                 :            :   { /* implies lc(x0) = 1 and L := 1/unscale is integral */
    1975                 :        105 :     long d = degpol(T->x0);
    1976                 :        105 :     GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
    1977                 :        105 :     GEN f= powiu(L, (d*(d-1)) >> 1);
    1978                 :        105 :     T->x = T->x0; /* restore original user-supplied x0, unscale data */
    1979                 :        105 :     T->unscale = gen_1;
    1980                 :        105 :     T->dx    = gmul(T->dx, sqri(f));
    1981                 :        105 :     T->bas   = RgXV_unscale(T->bas, unscale);
    1982                 :        105 :     T->index = gmul(T->index, f);
    1983                 :            :   }
    1984                 :       5054 :   nfbasic_add_disc(T); /* more expensive after set_LLL_basis */
    1985         [ +  + ]:       5054 :   if (flag & nf_RED)
    1986                 :            :   {
    1987                 :            :     GEN ro, rev;
    1988                 :            :     /* lie to polred: more efficient to update *after* modreverse, than to
    1989                 :            :      * unscale in the polred subsystem */
    1990                 :         84 :     T->unscale = gen_1;
    1991                 :         84 :     rev = nfpolred(T, &ro);
    1992                 :         84 :     nf = nfbasic_to_nf(T, ro, prec);
    1993         [ +  + ]:         84 :     if (flag & nf_ORIG)
    1994                 :            :     {
    1995         [ +  + ]:         56 :       if (!rev) rev = pol_x(varn(T->x)); /* no improvement */
    1996         [ +  + ]:         56 :       if (!isint1(unscale)) rev = RgX_Rg_div(rev, unscale);
    1997                 :         56 :       nf = mkvec2(nf, mkpolmod(rev, T->x));
    1998                 :            :     }
    1999                 :         84 :     T->unscale = unscale; /* restore */
    2000                 :            :   } else {
    2001                 :       4970 :     GEN ro; set_LLL_basis(T, &ro, 0.99);
    2002                 :       4970 :     nf = nfbasic_to_nf(T, ro, prec);
    2003                 :            :   }
    2004                 :       5054 :   return nf;
    2005                 :            : }
    2006                 :            : /* Initialize the number field defined by the polynomial x (in variable v)
    2007                 :            :  * flag & nf_RED:     try a polred first.
    2008                 :            :  * flag & nf_ORIG
    2009                 :            :  *    do a polred and return [nfinit(x), Mod(a,red)], where
    2010                 :            :  *    Mod(a,red) = Mod(v,x) (i.e return the base change). */
    2011                 :            : GEN
    2012                 :       3073 : nfinitall(GEN x, long flag, long prec)
    2013                 :            : {
    2014                 :       3073 :   const pari_sp av = avma;
    2015                 :            :   nfbasic_t T;
    2016                 :            :   GEN nf;
    2017                 :            : 
    2018         [ +  + ]:       3073 :   if (checkrnf_i(x)) return check_and_build_nfabs(x, prec);
    2019                 :       3066 :   nfinit_step1(&T, x, flag);
    2020                 :       3045 :   nf = nfinit_step2(&T, flag, prec);
    2021                 :       3052 :   return gerepilecopy(av, nf);
    2022                 :            : }
    2023                 :            : 
    2024                 :            : GEN
    2025                 :          0 : nfinitred(GEN x, long prec)  { return nfinitall(x, nf_RED, prec); }
    2026                 :            : GEN
    2027                 :          0 : nfinitred2(GEN x, long prec) { return nfinitall(x, nf_RED|nf_ORIG, prec); }
    2028                 :            : GEN
    2029                 :       1183 : nfinit(GEN x, long prec)     { return nfinitall(x, 0, prec); }
    2030                 :            : 
    2031                 :            : GEN
    2032                 :       1673 : nfinit0(GEN x, long flag,long prec)
    2033                 :            : {
    2034   [ +  +  +  - ]:       1673 :   switch(flag)
    2035                 :            :   {
    2036                 :            :     case 0:
    2037                 :       1652 :     case 1: return nfinitall(x,0,prec);
    2038                 :         14 :     case 2: case 4: return nfinitall(x,nf_RED,prec);
    2039                 :          7 :     case 3: case 5: return nfinitall(x,nf_RED|nf_ORIG,prec);
    2040                 :          0 :     default: pari_err_FLAG("nfinit");
    2041                 :            :   }
    2042                 :       1652 :   return NULL; /* not reached */
    2043                 :            : }
    2044                 :            : 
    2045                 :            : /* assume x a bnr/bnf/nf */
    2046                 :            : long
    2047                 :      68217 : nf_get_prec(GEN x)
    2048                 :            : {
    2049                 :      68217 :   GEN nf = checknf(x), ro = nf_get_roots(nf);
    2050         [ +  - ]:      68217 :   return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
    2051                 :            : }
    2052                 :            : 
    2053                 :            : /* assume nf is an nf */
    2054                 :            : GEN
    2055                 :        752 : nfnewprec_shallow(GEN nf, long prec)
    2056                 :            : {
    2057                 :        752 :   GEN NF = leafcopy(nf);
    2058                 :            :   nffp_t F;
    2059                 :        752 :   gel(NF,5) = leafcopy(gel(NF,5));
    2060                 :        752 :   remake_GM(NF, &F, prec);
    2061                 :        752 :   gel(NF,6) = F.ro;
    2062                 :        752 :   gmael(NF,5,1) = F.M;
    2063                 :        752 :   gmael(NF,5,2) = F.G;
    2064                 :        752 :   return NF;
    2065                 :            : }
    2066                 :            : 
    2067                 :            : GEN
    2068                 :         63 : nfnewprec(GEN nf, long prec)
    2069                 :            : {
    2070                 :            :   GEN z;
    2071   [ +  +  +  - ]:         63 :   switch(nftyp(nf))
    2072                 :            :   {
    2073                 :         49 :     default: pari_err_TYPE("nfnewprec", nf);
    2074                 :          7 :     case typ_BNF: z = bnfnewprec(nf,prec); break;
    2075                 :          7 :     case typ_BNR: z = bnrnewprec(nf,prec); break;
    2076                 :            :     case typ_NF: {
    2077                 :          0 :       pari_sp av = avma;
    2078                 :          0 :       z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
    2079                 :          0 :       break;
    2080                 :            :     }
    2081                 :            :   }
    2082                 :         14 :   return z;
    2083                 :            : }
    2084                 :            : 
    2085                 :            : /********************************************************************/
    2086                 :            : /**                                                                **/
    2087                 :            : /**                           POLRED                               **/
    2088                 :            : /**                                                                **/
    2089                 :            : /********************************************************************/
    2090                 :            : GEN
    2091                 :          0 : embednorm_T2(GEN x, long r1)
    2092                 :            : {
    2093                 :          0 :   pari_sp av = avma;
    2094                 :          0 :   GEN p = RgV_sumpart(x, r1);
    2095                 :          0 :   GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
    2096         [ #  # ]:          0 :   if (q != gen_0) p = gadd(p, gmul2n(q,1));
    2097         [ #  # ]:          0 :   return avma == av? gcopy(p): gerepileupto(av, p);
    2098                 :            : }
    2099                 :            : 
    2100                 :            : /* simplified version of gnorm for scalar, non-complex inputs, without GC */
    2101                 :            : static GEN
    2102                 :       6671 : real_norm(GEN x)
    2103                 :            : {
    2104   [ -  +  -  - ]:       6671 :   switch(typ(x))
    2105                 :            :   {
    2106                 :          0 :     case t_INT:  return sqri(x);
    2107                 :       6671 :     case t_REAL: return sqrr(x);
    2108                 :          0 :     case t_FRAC: return sqrfrac(x);
    2109                 :            :   }
    2110                 :          0 :   pari_err_TYPE("real_norm", x);
    2111                 :       6671 :   return NULL;
    2112                 :            : }
    2113                 :            : /* simplified version of gnorm, without GC */
    2114                 :            : static GEN
    2115                 :    1870303 : complex_norm(GEN x)
    2116                 :            : {
    2117         [ +  - ]:    1870303 :   return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
    2118                 :            : }
    2119                 :            : /* return T2(x), argument r1 needed in case x has components whose type
    2120                 :            :  * is unexpected, e.g. all of them t_INT for embed(gen_1) */
    2121                 :            : GEN
    2122                 :       1583 : embed_T2(GEN x, long r1)
    2123                 :            : {
    2124                 :       1583 :   pari_sp av = avma;
    2125                 :       1583 :   long i, l = lg(x);
    2126                 :       1583 :   GEN c, s = NULL, t = NULL;
    2127         [ -  + ]:       1583 :   if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
    2128         [ +  + ]:       8254 :   for (i = 1; i <= r1; i++)
    2129                 :            :   {
    2130                 :       6671 :     c = real_norm(gel(x,i));
    2131         [ +  + ]:       6671 :     s = s? gadd(s, c): c;
    2132                 :            :   }
    2133         [ +  + ]:       6338 :   for (; i < l; i++)
    2134                 :            :   {
    2135                 :       4755 :     c = complex_norm(gel(x,i));
    2136         [ +  + ]:       4755 :     t = t? gadd(t, c): c;
    2137                 :            :   }
    2138 [ +  + ][ +  + ]:       1583 :   if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
    2139                 :       1583 :   return gerepileupto(av, s);
    2140                 :            : }
    2141                 :            : /* return N(x) */
    2142                 :            : GEN
    2143                 :    1015083 : embed_norm(GEN x, long r1)
    2144                 :            : {
    2145                 :    1015083 :   pari_sp av = avma;
    2146                 :    1015083 :   long i, l = lg(x);
    2147                 :    1015083 :   GEN c, s = NULL, t = NULL;
    2148         [ +  + ]:    1015083 :   if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
    2149         [ +  + ]:    2688579 :   for (i = 1; i <= r1; i++)
    2150                 :            :   {
    2151                 :    1676966 :     c = gel(x,i);
    2152         [ +  + ]:    1676966 :     s = s? gmul(s, c): c;
    2153                 :            :   }
    2154         [ +  + ]:    2877161 :   for (; i < l; i++)
    2155                 :            :   {
    2156                 :    1865548 :     c = complex_norm(gel(x,i));
    2157         [ +  + ]:    1865548 :     t = t? gmul(t, c): c;
    2158                 :            :   }
    2159 [ +  + ][ +  + ]:    1011613 :   if (t) s = s? gmul(s,t): t;
    2160                 :    1015083 :   return gerepileupto(av, s);
    2161                 :            : }
    2162                 :            : 
    2163                 :            : typedef struct {
    2164                 :            :   long r1, v, prec;
    2165                 :            :   GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
    2166                 :            :   GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
    2167                 :            :   GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
    2168                 :            :   GEN bound; /* T2 norm of the polynomial defining nf */
    2169                 :            :   long expo_best_disc; /* expo(disc(x)), best generator so far */
    2170                 :            : } CG_data;
    2171                 :            : 
    2172                 :            : /* characteristic pol of x (given by embeddings) */
    2173                 :            : static GEN
    2174                 :      24778 : get_pol(CG_data *d, GEN x)
    2175                 :            : {
    2176                 :            :   long e;
    2177                 :      24778 :   GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
    2178         [ +  + ]:      24778 :   return (e > -5)? NULL: g;
    2179                 :            : }
    2180                 :            : 
    2181                 :            : /* characteristic pol of x (given as vector on (w_i)) */
    2182                 :            : static GEN
    2183                 :      10593 : get_polchar(CG_data *d, GEN x)
    2184                 :      10593 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
    2185                 :            : 
    2186                 :            : /* Choose a canonical polynomial in the pair { z(X), (+/-)z(-X) }.
    2187                 :            :  * z a ZX with lc(z) > 0. We want to keep that property, while
    2188                 :            :  * ensuring that the leading coeff of the odd (resp. even) part of z is < 0
    2189                 :            :  * if deg z is even (resp. odd).
    2190                 :            :  * Either leave z alone (return 1) or set z <-- (-1)^deg(z) z(-X). In place. */
    2191                 :            : static int
    2192                 :      10733 : ZX_canon_neg(GEN z)
    2193                 :            : {
    2194                 :            :   long i,s;
    2195                 :            : 
    2196         [ +  + ]:      82889 :   for (i = lg(z)-2; i >= 2; i -= 2)
    2197                 :            :   { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
    2198                 :      77576 :     s = signe(gel(z,i));
    2199         [ +  + ]:      77576 :     if (!s) continue;
    2200                 :            :     /* non trivial */
    2201         [ +  + ]:       5420 :     if (s < 0) break; /* the condition is already satisfied */
    2202                 :            : 
    2203         [ +  + ]:       8997 :     for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
    2204                 :       1807 :     return 1;
    2205                 :            :   }
    2206                 :      10733 :   return 0;
    2207                 :            : }
    2208                 :            : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
    2209                 :            :  * Return NULL on failure: discriminant too large or non primitive */
    2210                 :            : static GEN
    2211                 :      20272 : try_polmin(CG_data *d, nfbasic_t *T, GEN v, long flag, GEN *ai)
    2212                 :            : {
    2213                 :      20272 :   const long best = flag & nf_ABSOLUTE;
    2214                 :            :   long ed;
    2215                 :      20272 :   pari_sp av = avma;
    2216                 :            :   GEN g;
    2217         [ +  + ]:      20272 :   if (best)
    2218                 :            :   {
    2219                 :      19397 :     ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
    2220         [ +  + ]:      19397 :     avma = av; if (d->expo_best_disc < ed) return NULL;
    2221                 :            :   }
    2222                 :            :   else
    2223                 :        875 :     ed = 0;
    2224                 :      11683 :   g = get_pol(d, v);
    2225                 :            :   /* accuracy too low, compute algebraically */
    2226         [ +  + ]:      11683 :   if (!g) { avma = av; g = ZXQ_charpoly(*ai, T->x, varn(T->x)); }
    2227                 :      11683 :   (void)ZX_gcd_all(g, ZX_deriv(g), &g);
    2228 [ +  + ][ +  + ]:      11683 :   if (best && degpol(g) != degpol(T->x)) { avma = av; return NULL; }
    2229                 :       4718 :   g = gerepilecopy(av, g);
    2230                 :       4718 :   d->expo_best_disc = ed;
    2231         [ +  + ]:       4718 :   if (flag & nf_ORIG)
    2232                 :            :   {
    2233         [ +  + ]:       1001 :     if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
    2234         [ +  + ]:       1001 :     if (!isint1(T->unscale)) *ai = RgX_unscale(*ai, T->unscale);
    2235                 :            :   }
    2236                 :            :   else
    2237                 :       3717 :     (void)ZX_canon_neg(g);
    2238         [ -  + ]:       4718 :   if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
    2239                 :      20272 :   return g;
    2240                 :            : }
    2241                 :            : 
    2242                 :            : /* does x generate the correct field ? */
    2243                 :            : static GEN
    2244                 :      10593 : chk_gen(void *data, GEN x)
    2245                 :            : {
    2246                 :      10593 :   pari_sp av = avma, av1;
    2247                 :      10593 :   GEN h, g = get_polchar((CG_data*)data,x);
    2248         [ -  + ]:      10593 :   if (!g) pari_err_PREC("chk_gen");
    2249                 :      10593 :   av1 = avma;
    2250                 :      10593 :   h = ZX_gcd(g, ZX_deriv(g));
    2251         [ +  + ]:      10593 :   if (degpol(h)) { avma = av; return NULL; }
    2252         [ -  + ]:       6050 :   if (DEBUGLEVEL>3) err_printf("  generator: %Ps\n",g);
    2253                 :      10593 :   avma = av1; return gerepileupto(av, g);
    2254                 :            : }
    2255                 :            : 
    2256                 :            : static long
    2257                 :       1456 : chk_gen_prec(long N, long bit)
    2258                 :       1456 : { return nbits2prec(10 + (long)log2((double)N) + bit); }
    2259                 :            : 
    2260                 :            : /* Remove duplicate polynomials in P, updating A (same indices), in place.
    2261                 :            :  * Among elements having the same characteristic pol, choose the smallest
    2262                 :            :  * according to ZV_abscmp */
    2263                 :            : static void
    2264                 :        343 : remove_duplicates(GEN P, GEN A)
    2265                 :            : {
    2266                 :        343 :   long k, i, l = lg(P);
    2267                 :        343 :   pari_sp av = avma;
    2268                 :            :   GEN x, a;
    2269                 :            : 
    2270         [ -  + ]:        686 :   if (l < 2) return;
    2271                 :        343 :   (void)sort_factor_pol(mkmat2(P, A), cmpii);
    2272                 :        343 :   x = gel(P,1); a = gel(A,1);
    2273         [ +  + ]:       6015 :   for  (k=1,i=2; i<l; i++)
    2274         [ +  + ]:       5672 :     if (ZX_equal(gel(P,i), x))
    2275                 :            :     {
    2276         [ +  + ]:       3480 :       if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
    2277                 :            :     }
    2278                 :            :     else
    2279                 :            :     {
    2280                 :       2192 :       gel(A,k) = a;
    2281                 :       2192 :       gel(P,k) = x;
    2282                 :       2192 :       k++;
    2283                 :       2192 :       x = gel(P,i); a = gel(A,i);
    2284                 :            :     }
    2285                 :        343 :   l = k+1;
    2286                 :        343 :   gel(A,k) = a; setlg(A,l);
    2287                 :        343 :   gel(P,k) = x; setlg(P,l); avma = av;
    2288                 :            : }
    2289                 :            : 
    2290                 :            : static long
    2291                 :       1113 : polred_init(nfbasic_t *T, nffp_t *F, CG_data *d)
    2292                 :            : {
    2293                 :       1113 :   long e, prec, n = degpol(T->x);
    2294                 :            :   double log2rho;
    2295                 :            :   GEN ro;
    2296                 :       1113 :   set_LLL_basis(T, &ro, 0.9999);
    2297                 :            :   /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
    2298         [ +  + ]:       1113 :   log2rho = ro ? (double)gexpo(ro): fujiwara_bound(T->x);
    2299                 :       1113 :   e = n * (long)(log2rho + log2((double)n)) + 1;
    2300         [ -  + ]:       1113 :   if (e < 0) e = 0; /* can occur if n = 1 */
    2301                 :       1113 :   prec = chk_gen_prec(n, e);
    2302                 :       1113 :   get_nf_fp_compo(T, F, ro, 1, prec);
    2303                 :       1113 :   d->v = varn(T->x);
    2304                 :       1113 :   d->expo_best_disc = -1;
    2305                 :       1113 :   d->ZKembed = NULL;
    2306                 :       1113 :   d->M = NULL;
    2307                 :       1113 :   d->u = NULL;
    2308                 :       1113 :   d->r1= T->r1; return prec;
    2309                 :            : }
    2310                 :            : static GEN
    2311                 :        350 : findmindisc(GEN y, GEN *pa)
    2312                 :            : {
    2313                 :        350 :   GEN a = *pa, x = gel(y,1), b = gel(a,1), dx = NULL;
    2314                 :        350 :   long i, l = lg(y);
    2315         [ +  + ]:        421 :   for (i = 2; i < l; i++)
    2316                 :            :   {
    2317                 :         71 :     GEN yi = gel(y,i);
    2318         [ +  + ]:         71 :     if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
    2319                 :            :   }
    2320                 :        350 :   *pa = b; return x;
    2321                 :            : }
    2322                 :            : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
    2323                 :            :  * [ the best wrt discriminant ordering ], but keep all non-primitive
    2324                 :            :  * polynomials */
    2325                 :            : static void
    2326                 :        770 : filter(GEN y, GEN b, long n)
    2327                 :            : {
    2328                 :            :   GEN x, a, dx;
    2329                 :        770 :   long i, k = 1, l = lg(y);
    2330                 :        770 :   a = x = dx = NULL;
    2331         [ +  + ]:       5551 :   for (i = 1; i < l; i++)
    2332                 :            :   {
    2333                 :       4781 :     GEN yi = gel(y,i), ai = gel(b,i);
    2334         [ +  + ]:       4781 :     if (degpol(yi) == n)
    2335                 :            :     {
    2336                 :       4613 :       pari_sp av = avma;
    2337 [ +  + ][ +  + ]:       4613 :       if (dx && !ZX_is_better(yi,x,&dx)) { avma = av; continue; }
    2338         [ +  + ]:       1085 :       if (!dx) dx = ZX_disc(yi);
    2339                 :       1085 :       x = yi; a = ai; continue;
    2340                 :            :     }
    2341                 :        168 :     gel(y,k) = yi;
    2342                 :        168 :     gel(b,k) = ai; k++;
    2343                 :            :   }
    2344         [ +  - ]:        770 :   if (dx)
    2345                 :            :   {
    2346                 :        770 :     gel(y,k) = x;
    2347                 :        770 :     gel(b,k) = a; k++;
    2348                 :            :   }
    2349                 :        770 :   setlg(y, k);
    2350                 :        770 :   setlg(b, k);
    2351                 :        770 : }
    2352                 :            : 
    2353                 :            : static GEN
    2354                 :        798 : polred_aux(nfbasic_t *T, GEN *pro, long flag)
    2355                 :            : { /* only keep polynomials of max degree and best discriminant */
    2356                 :        798 :   const long best = flag & nf_ABSOLUTE;
    2357                 :        798 :   const long orig = flag & nf_ORIG;
    2358                 :        798 :   GEN M, b, y, x = T->x;
    2359                 :        798 :   long maxi, i, j, k, v = varn(x), n = lg(T->bas)-1;
    2360                 :            :   nffp_t F;
    2361                 :            :   CG_data d;
    2362                 :            : 
    2363         [ +  + ]:        798 :   if (n == 1)
    2364                 :            :   {
    2365         [ +  + ]:         28 :     if (!best)
    2366                 :            :     {
    2367                 :         14 :       GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2368         [ -  + ]:         14 :       return orig? mkmat2(mkcol(ch),mkcol(gen_1)): mkvec(ch);
    2369                 :            :     }
    2370                 :            :     else
    2371         [ +  + ]:         14 :       return orig? trivial_fact(): cgetg(1,t_VEC);
    2372                 :            :   }
    2373                 :            : 
    2374                 :        770 :   (void)polred_init(T, &F, &d);
    2375                 :        770 :   *pro = F.ro;
    2376                 :        770 :   M = F.M;
    2377         [ +  + ]:        770 :   if (best)
    2378                 :            :   {
    2379         [ -  + ]:        707 :     if (!T->dx) T->dx = ZX_disc(T->x);
    2380                 :        707 :     d.expo_best_disc = expi(T->dx);
    2381                 :            :   }
    2382                 :            : 
    2383                 :            :   /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
    2384                 :        770 :   y = cgetg(n*n + 1, t_VEC);
    2385                 :        770 :   b = cgetg(n*n + 1, t_COL);
    2386                 :        770 :   k = 1;
    2387         [ +  + ]:        770 :   if (!best)
    2388                 :            :   {
    2389                 :         63 :     GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2390                 :         63 :     gel(y,1) = ch; gel(b,1) = gen_1; k++;
    2391                 :            :   }
    2392         [ +  + ]:       4270 :   for (i = 2; i <= n; i++)
    2393                 :            :   {
    2394                 :            :     GEN ch, ai;
    2395                 :       3500 :     ai = gel(T->bas,i);
    2396                 :       3500 :     ch = try_polmin(&d, T, gel(M,i), flag, &ai);
    2397         [ +  + ]:       3500 :     if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2398                 :            :   }
    2399                 :        770 :   maxi = minss(n, 3);
    2400         [ +  + ]:       2884 :   for (i = 1; i <= maxi; i++)
    2401         [ +  + ]:      10500 :     for (j = i+1; j <= n; j++)
    2402                 :            :     {
    2403                 :            :       GEN ch, ai, v;
    2404                 :       8386 :       ai = gadd(gel(T->bas,i), gel(T->bas,j));
    2405                 :       8386 :       v = RgV_add(gel(M,i), gel(M,j));
    2406                 :            :       /* defining polynomial for Q(w_i+w_j) */
    2407                 :       8386 :       ch = try_polmin(&d, T, v, flag, &ai);
    2408         [ +  + ]:       8386 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2409                 :            : 
    2410                 :       8386 :       ai = gsub(gel(T->bas,i), gel(T->bas,j));
    2411                 :       8386 :       v = RgV_sub(gel(M,i), gel(M,j));
    2412                 :            :       /* defining polynomial for Q(w_i-w_j) */
    2413                 :       8386 :       ch = try_polmin(&d, T, v, flag, &ai);
    2414         [ +  + ]:       8386 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2415                 :            :     }
    2416                 :        770 :   setlg(y, k);
    2417                 :        770 :   setlg(b, k); filter(y, b, n);
    2418         [ +  + ]:        770 :   if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
    2419                 :        147 :   (void)sort_factor_pol(mkmat2(y, b), cmpii);
    2420                 :        798 :   settyp(y, t_COL); return mkmat2(b, y);
    2421                 :            : }
    2422                 :            : 
    2423                 :            : static GEN
    2424                 :         84 : Polred(GEN x, long flag, GEN fa)
    2425                 :            : {
    2426                 :         84 :   pari_sp av = avma;
    2427                 :            :   GEN ro;
    2428         [ +  + ]:         84 :   nfbasic_t T; nfbasic_init(fa? mkvec2(x,fa): x, flag & nf_PARTIALFACT, &T);
    2429                 :         77 :   return gerepilecopy(av, polred_aux(&T, &ro, flag));
    2430                 :            : }
    2431                 :            : 
    2432                 :            : /* finds "best" polynomial in polred_aux list, defaulting to T->x if none of
    2433                 :            :  * them is primitive. *px is the ZX, characteristic polynomial of Mod(*pb,T->x),
    2434                 :            :  * *pdx its discriminant. Set *pro = polroots(T->x) [ NOT *px ]. */
    2435                 :            : static void
    2436                 :        721 : polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
    2437                 :            : {
    2438                 :        721 :   GEN y, x = T->x; /* default value */
    2439                 :            :   long i, l;
    2440         [ +  + ]:        721 :   y = polred_aux(T, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
    2441                 :        721 :   *pdx = T->dx;
    2442         [ +  + ]:        721 :   if (pb)
    2443                 :            :   {
    2444                 :        133 :     GEN a, b = deg1pol_shallow(T->unscale, gen_0, varn(x));
    2445                 :        133 :     a = gel(y,1); l = lg(a);
    2446                 :        133 :     y = gel(y,2);
    2447         [ +  + ]:        259 :     for (i=1; i<l; i++)
    2448                 :            :     {
    2449                 :        126 :       GEN yi = gel(y,i);
    2450                 :        126 :       pari_sp av = avma;
    2451         [ +  + ]:        126 :       if (ZX_is_better(yi,x,pdx)) { x = yi; b = gel(a,i); } else avma = av;
    2452                 :            :     }
    2453                 :        133 :     *pb = b;
    2454                 :            :   }
    2455                 :            :   else
    2456                 :            :   {
    2457                 :        588 :     l = lg(y);
    2458         [ +  + ]:       1169 :     for (i=1; i<l; i++)
    2459                 :            :     {
    2460                 :        581 :       GEN yi = gel(y,i);
    2461                 :        581 :       pari_sp av = avma;
    2462         [ +  + ]:        581 :       if (ZX_is_better(yi,x,pdx)) x = yi; else avma = av;
    2463                 :            :     }
    2464                 :            :   }
    2465         [ -  + ]:        721 :   if (!*pdx) *pdx = ZX_disc(x);
    2466                 :        721 :   *px = x;
    2467                 :        721 : }
    2468                 :            : GEN
    2469                 :        616 : polredbest(GEN T0, long flag)
    2470                 :            : {
    2471                 :        616 :   pari_sp av = avma;
    2472                 :            :   GEN T, dT, ro, a;
    2473                 :            :   nfbasic_t S;
    2474 [ +  - ][ -  + ]:        616 :   if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
    2475                 :        616 :   T = T0; nfbasic_init(T, nf_PARTIALFACT, &S);
    2476         [ +  + ]:        616 :   polredbest_aux(&S, &ro, &T, &dT, flag? &a: NULL);
    2477         [ +  + ]:        616 :   if (flag)
    2478                 :            :   { /* charpoly(Mod(a,T0)) = T */
    2479                 :            :     GEN b;
    2480         [ -  + ]:         28 :     if (T0 == T)
    2481                 :          0 :       b = pol_x(varn(T)); /* no improvement */
    2482                 :            :     else
    2483                 :         28 :       b = QXQ_reverse(a, T0); /* charpoly(Mod(b,T)) = S.x */
    2484         [ +  + ]:         28 :     b = (degpol(T) == 1)? gmodulo(b, T): mkpolmod(b,T);
    2485                 :         28 :     T = mkvec2(T, b);
    2486                 :            :   }
    2487                 :        616 :   return gerepilecopy(av, T);
    2488                 :            : }
    2489                 :            : /* DEPRECATED: backward compatibility */
    2490                 :            : GEN
    2491                 :         70 : polred0(GEN x, long flag, GEN fa)
    2492                 :            : {
    2493                 :         70 :   long fl = 0;
    2494         [ +  + ]:         70 :   if (flag & 1) fl |= nf_PARTIALFACT;
    2495         [ +  + ]:         70 :   if (flag & 2) fl |= nf_ORIG;
    2496                 :         70 :   return Polred(x, fl, fa);
    2497                 :            : }
    2498                 :            : 
    2499                 :            : GEN
    2500                 :         21 : polredord(GEN x)
    2501                 :            : {
    2502                 :         21 :   pari_sp av = avma;
    2503                 :            :   GEN v, lt;
    2504                 :            :   long i, n, vx;
    2505                 :            : 
    2506         [ -  + ]:         21 :   if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
    2507                 :         21 :   x = Q_primpart(x); RgX_check_ZX(x,"polredord");
    2508         [ -  + ]:         21 :   n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
    2509         [ +  + ]:         21 :   if (n == 1) return gerepilecopy(av, mkvec(x));
    2510                 :         14 :   lt = leading_coeff(x); vx = varn(x);
    2511         [ +  + ]:         14 :   if (is_pm1(lt))
    2512                 :            :   {
    2513         [ -  + ]:          7 :     if (signe(lt) < 0) x = ZX_neg(x);
    2514                 :          7 :     v = pol_x_powers(n, vx);
    2515                 :            :   }
    2516                 :            :   else
    2517                 :            :   { GEN L;
    2518                 :            :     /* basis for Dedekind order */
    2519                 :          7 :     v = cgetg(n+1, t_VEC);
    2520                 :          7 :     gel(v,1) = scalarpol_shallow(lt, vx);
    2521         [ +  + ]:         14 :     for (i = 2; i <= n; i++)
    2522                 :          7 :       gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
    2523                 :          7 :     gel(v,1) = pol_1(vx);
    2524                 :          7 :     x = ZX_Q_normalize(x, &L);
    2525                 :          7 :     v = gsubst(v, vx, monomial(ginv(L),1,vx));
    2526         [ +  + ]:         14 :     for (i=2; i <= n; i++)
    2527         [ +  - ]:          7 :       if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = monomial(gen_1, i-1, vx);
    2528                 :            :   }
    2529                 :         21 :   return gerepileupto(av, polred(mkvec2(x, v)));
    2530                 :            : }
    2531                 :            : 
    2532                 :            : GEN
    2533                 :         14 : polred(GEN x) { return Polred(x, 0, NULL); }
    2534                 :            : GEN
    2535                 :          0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
    2536                 :            : GEN
    2537                 :          0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
    2538                 :            : GEN
    2539                 :          0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
    2540                 :            : GEN
    2541                 :          0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
    2542                 :            : GEN
    2543                 :          0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
    2544                 :            : 
    2545                 :            : /********************************************************************/
    2546                 :            : /**                                                                **/
    2547                 :            : /**                           POLREDABS                            **/
    2548                 :            : /**                                                                **/
    2549                 :            : /********************************************************************/
    2550                 :            : /* set V[k] := matrix of multiplication by nk.zk[k] */
    2551                 :            : static GEN
    2552                 :       1302 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
    2553                 :            : {
    2554                 :       1302 :   GEN v, Mk = cgetg(N+1, t_MAT);
    2555                 :            :   long i, e;
    2556         [ +  + ]:       4459 :   for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
    2557         [ +  + ]:      12320 :   for (     ; i <=N; i++)
    2558                 :            :   {
    2559                 :      11018 :     v = vecmul(gel(M,k), gel(M,i));
    2560                 :      11018 :     v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
    2561                 :      11018 :     gel(Mk,i) = grndtoi(v, &e);
    2562         [ -  + ]:      11018 :     if (e > -5) return NULL;
    2563                 :            :   }
    2564                 :       1302 :   gel(V,k) = Mk; return Mk;
    2565                 :            : }
    2566                 :            : 
    2567                 :            : static GEN
    2568                 :       1427 : ZM_image_shallow(GEN M, long *pr)
    2569                 :            : {
    2570                 :            :   long j, k, r;
    2571                 :       1427 :   GEN y, d = ZM_pivots(M, &k);
    2572                 :       1427 :   r = lg(M)-1 - k;
    2573                 :       1427 :   y = cgetg(r+1,t_MAT);
    2574         [ +  + ]:      10550 :   for (j=k=1; j<=r; k++)
    2575         [ +  + ]:       9123 :     if (d[k]) gel(y,j++) = gel(M,k);
    2576                 :       1427 :   *pr = r; return y;
    2577                 :            : }
    2578                 :            : 
    2579                 :            : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
    2580                 :            :  * Q from algo 2.7.6] */
    2581                 :            : static GEN
    2582                 :        351 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
    2583                 :            : {
    2584                 :        351 :   CG_data *d = (CG_data*)chk->data;
    2585                 :            :   GEN P, V, D, inv, bound, S, M;
    2586                 :        351 :   long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
    2587                 :        351 :   long i, j, prec, firstprim = 0, skipfirst = 0;
    2588                 :            :   pari_sp av;
    2589                 :            : 
    2590                 :        351 :   d->u = U;
    2591                 :        351 :   d->ZKembed = M = RgM_mul(d->M, U);
    2592                 :            : 
    2593                 :        351 :   av = avma; bound = d->bound;
    2594                 :        351 :   D = cgetg(N+1, t_VECSMALL);
    2595         [ +  + ]:       2831 :   for (i = 1; i <= N; i++)
    2596                 :            :   {
    2597                 :       2488 :     pari_sp av2 = avma;
    2598                 :       2488 :     P = get_pol(d, gel(M,i));
    2599         [ +  + ]:       2488 :     if (!P) pari_err_PREC("chk_gen_init");
    2600                 :       2480 :     (void)ZX_gcd_all(P, ZX_deriv(P), &P);
    2601                 :       2480 :     P = gerepilecopy(av2, P);
    2602                 :       2480 :     D[i] = degpol(P);
    2603         [ +  + ]:       2480 :     if (D[i] == N)
    2604                 :            :     { /* primitive element */
    2605                 :       1030 :       GEN B = embed_T2(gel(M,i), r1);
    2606         [ +  + ]:       1030 :       if (!firstprim) firstprim = i; /* index of first primitive element */
    2607         [ -  + ]:       1030 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2608         [ +  + ]:       1030 :       if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
    2609                 :            :     }
    2610                 :            :     else
    2611                 :            :     {
    2612         [ -  + ]:       1450 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
    2613         [ +  + ]:       1450 :       if (firstprim)
    2614                 :            :       { /* cycle basis vectors so that primitive elements come last */
    2615                 :        159 :         GEN u = d->u, e = M;
    2616                 :        159 :         GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
    2617                 :        159 :         long tS = D[i];
    2618         [ +  + ]:        447 :         for (j = i; j > firstprim; j--)
    2619                 :            :         {
    2620                 :        288 :           u[j] = u[j-1];
    2621                 :        288 :           e[j] = e[j-1];
    2622                 :        288 :           R[j] = R[j-1];
    2623                 :        288 :           D[j] = D[j-1];
    2624                 :            :         }
    2625                 :        159 :         gel(u,firstprim) = tu;
    2626                 :        159 :         gel(e,firstprim) = te;
    2627                 :        159 :         gel(R,firstprim) = tR;
    2628                 :        159 :         D[firstprim] = tS; firstprim++;
    2629                 :            :       }
    2630                 :            :     }
    2631                 :            :   }
    2632         [ +  + ]:        343 :   if (!firstprim)
    2633                 :            :   { /* try (a little) to find primitive elements to improve bound */
    2634                 :         21 :     GEN x = cgetg(N+1, t_VECSMALL);
    2635         [ -  + ]:         21 :     if (DEBUGLEVEL>1)
    2636                 :          0 :       err_printf("chk_gen_init: difficult field, trying random elements\n");
    2637         [ +  + ]:        231 :     for (i = 0; i < 10; i++)
    2638                 :            :     {
    2639                 :            :       GEN e, B;
    2640         [ +  + ]:       2310 :       for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
    2641                 :        210 :       e = RgM_zc_mul(M, x);
    2642                 :        210 :       B = embed_T2(e, r1);
    2643         [ +  + ]:        210 :       if (gcmp(B,bound) >= 0) continue;
    2644         [ -  + ]:         14 :       P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
    2645         [ -  + ]:         14 :       if (!ZX_is_squarefree(P)) continue;
    2646         [ -  + ]:         14 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2647                 :         14 :       bound = B ;
    2648                 :            :     }
    2649                 :            :   }
    2650                 :            : 
    2651         [ +  - ]:        343 :   if (firstprim != 1)
    2652                 :            :   {
    2653                 :        343 :     inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
    2654                 :        343 :     V = gel(inv,1);
    2655         [ +  + ]:       1554 :     for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
    2656                 :            :     /* V corresponds to 1_Z */
    2657                 :        343 :     V = grndtoi(V, &j);
    2658         [ -  + ]:        343 :     if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
    2659                 :        343 :     S = mkmat(V); /* 1 */
    2660                 :            : 
    2661                 :        343 :     V = cgetg(N+1, t_VEC);
    2662         [ +  - ]:       1456 :     for (i = 1; i <= N; i++,skipfirst++)
    2663                 :            :     { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
    2664                 :            :       GEN Mx, M2;
    2665                 :       1456 :       long j, k, h, rkM, dP = D[i];
    2666                 :            : 
    2667         [ +  + ]:       1456 :       if (dP == N) break; /* primitive */
    2668                 :       1302 :       Mx = set_mulid(V, M, inv, r1, r2, N, i);
    2669         [ -  + ]:       1302 :       if (!Mx) break; /* prec. problem. Stop */
    2670         [ +  + ]:       1302 :       if (dP == 1) continue;
    2671                 :       1008 :       rkM = lg(S)-1;
    2672                 :       1008 :       M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
    2673                 :       1008 :       gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
    2674                 :       1008 :       k = 2;
    2675         [ +  + ]:       2435 :       for (h = 1; h < dP; h++)
    2676                 :            :       {
    2677                 :            :         long r; /* add to M2 the elts of S * nf.zk[i]  */
    2678         [ +  + ]:       7801 :         for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
    2679                 :       1427 :         setlg(M2, k); k = 1;
    2680                 :       1427 :         S = ZM_image_shallow(shallowconcat(S,M2), &r);
    2681         [ +  + ]:       1427 :         if (r == rkM) break;
    2682         [ +  - ]:        895 :         if (r > rkM)
    2683                 :            :         {
    2684                 :        895 :           rkM = r;
    2685         [ +  + ]:        895 :           if (rkM == N) break;
    2686                 :            :         }
    2687                 :            :       }
    2688         [ +  + ]:       1008 :       if (rkM == N) break;
    2689                 :            :       /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
    2690                 :            :     }
    2691                 :            :   }
    2692                 :            :   /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
    2693                 :        343 :   chk->skipfirst = skipfirst;
    2694         [ -  + ]:        343 :   if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
    2695                 :            : 
    2696                 :            :   /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
    2697                 :        343 :   bound = gerepileuptoleaf(av, bound);
    2698                 :        343 :   prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
    2699         [ -  + ]:        343 :   if (DEBUGLEVEL)
    2700                 :          0 :     err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
    2701         [ -  + ]:        343 :   if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
    2702         [ +  + ]:        343 :   if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
    2703                 :        343 :   return bound;
    2704                 :            : }
    2705                 :            : 
    2706                 :            : /* z "small" minimal polynomial of Mod(a,x), deg z = deg x */
    2707                 :            : static GEN
    2708                 :       2478 : store(GEN x, GEN z, GEN a, nfbasic_t *T, long flag, GEN u)
    2709                 :            : {
    2710                 :            :   GEN y, b;
    2711                 :            : 
    2712         [ +  + ]:       2478 :   if (u) a = RgV_RgC_mul(T->bas, ZM_ZC_mul(u, a));
    2713         [ +  + ]:       2478 :   if (flag & (nf_ORIG|nf_ADDZK))
    2714                 :            :   {
    2715                 :        245 :     b = QXQ_reverse(a, x);
    2716         [ +  + ]:        245 :     if (!isint1(T->unscale)) b = gdiv(b, T->unscale); /* not RgX_Rg_div */
    2717                 :            :   }
    2718                 :            :   else
    2719                 :       2233 :     b = NULL;
    2720                 :            : 
    2721         [ +  + ]:       2478 :   if (flag & nf_RAW)
    2722                 :         28 :     y = mkvec2(z, a);
    2723         [ +  + ]:       2450 :   else if (flag & nf_ORIG) /* store phi(b mod z). */
    2724                 :        245 :     y = mkvec2(z, mkpolmod(b,z));
    2725                 :            :   else
    2726                 :       2205 :     y = z;
    2727         [ -  + ]:       2478 :   if (flag & nf_ADDZK)
    2728                 :            :   { /* append integral basis for number field Q[X]/(z) to result */
    2729                 :          0 :     long n = degpol(x);
    2730                 :          0 :     GEN t = RgV_RgM_mul(RgXQ_powers(b, n-1, z), RgV_to_RgM(T->bas,n));
    2731                 :          0 :     y = mkvec2(y, t);
    2732                 :            :   }
    2733                 :       2478 :   return y;
    2734                 :            : }
    2735                 :            : static GEN
    2736                 :        343 : polredabs_aux(nfbasic_t *T, GEN *u)
    2737                 :            : {
    2738                 :            :   long prec;
    2739                 :            :   GEN v;
    2740                 :        343 :   FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
    2741                 :            :   nffp_t F;
    2742                 :        343 :   CG_data d; chk.data = (void*)&d;
    2743                 :            : 
    2744                 :        343 :   prec = polred_init(T, &F, &d);
    2745                 :        343 :   d.bound = embed_T2(F.ro, d.r1);
    2746         [ +  + ]:        343 :   if (realprec(d.bound) > prec) d.bound = rtor(d.bound, prec);
    2747                 :            :   for (;;)
    2748                 :            :   {
    2749                 :        363 :     GEN R = R_from_QR(F.G, prec);
    2750         [ +  + ]:        363 :     if (R)
    2751                 :            :     {
    2752                 :        351 :       d.prec = prec;
    2753                 :        351 :       d.M    = F.M;
    2754                 :        351 :       v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
    2755         [ +  + ]:        351 :       if (v) break;
    2756                 :            :     }
    2757                 :         20 :     prec = precdbl(prec);
    2758                 :         20 :     get_nf_fp_compo(T, &F, NULL, 1, prec);
    2759         [ -  + ]:         20 :     if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",prec);
    2760                 :         20 :   }
    2761                 :        343 :   *u = d.u; return v;
    2762                 :            : }
    2763                 :            : 
    2764                 :            : GEN
    2765                 :        357 : polredabs0(GEN x, long flag)
    2766                 :            : {
    2767                 :        357 :   pari_sp av = avma;
    2768                 :            :   long i, l, vx;
    2769                 :            :   GEN y, a, u;
    2770                 :            :   nfbasic_t T;
    2771                 :            : 
    2772                 :        357 :   nfbasic_init(x, nf_PARTIALFACT, &T);
    2773                 :        357 :   x = T.x; vx = varn(x);
    2774                 :            : 
    2775         [ +  + ]:        357 :   if (degpol(x) == 1)
    2776                 :            :   {
    2777                 :         14 :     u = NULL;
    2778                 :         14 :     y = mkvec( pol_x(vx) );
    2779                 :         14 :     a = mkvec( deg1pol_shallow(gen_1, negi(gel(x,2)), vx) );
    2780                 :         14 :     l = 2;
    2781                 :            :   }
    2782                 :            :   else
    2783                 :            :   {
    2784                 :            :     GEN v;
    2785 [ +  + ][ +  + ]:        343 :     if (!(flag & nf_PARTIALFACT) && T.dKP)
    2786                 :            :     {
    2787                 :         84 :       GEN vw = primes_certify(T.dK, T.dKP);
    2788                 :         84 :       v = gel(vw,1); l = lg(v);
    2789         [ +  + ]:         84 :       if (l != 1)
    2790                 :            :       { /* fix integral basis */
    2791                 :          7 :         GEN w = gel(vw,2);
    2792         [ +  + ]:         14 :         for (i = 1; i < l; i++)
    2793                 :          7 :           w = ZV_union_shallow(w, gel(Z_factor(gel(v,i)),1));
    2794                 :          7 :         nfbasic_init(mkvec2(x,w), 0, &T);
    2795                 :            :       }
    2796                 :            :     }
    2797                 :        343 :     v = polredabs_aux(&T, &u);
    2798                 :        343 :     y = gel(v,1);
    2799                 :        343 :     a = gel(v,2); l = lg(a);
    2800         [ +  + ]:       6358 :     for (i=1; i<l; i++)
    2801         [ +  + ]:       6015 :       if (ZX_canon_neg(gel(y,i))) gel(a,i) = ZC_neg(gel(a,i));
    2802                 :        343 :     remove_duplicates(y,a);
    2803                 :        343 :     l = lg(a);
    2804         [ -  + ]:        343 :     if (l == 1)
    2805                 :          0 :       pari_err_BUG("polredabs (missing vector)");
    2806                 :            :   }
    2807         [ -  + ]:        357 :   if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
    2808         [ +  + ]:        357 :   if (flag & nf_ALL) {
    2809         [ +  + ]:       2135 :     for (i=1; i<l; i++) gel(y,i) = store(x, gel(y,i), gel(a,i), &T, flag, u);
    2810                 :            :   } else {
    2811                 :        350 :     GEN z = findmindisc(y, &a);
    2812                 :        350 :     y = store(x, z, a, &T, flag, u);
    2813                 :            :   }
    2814                 :        357 :   return gerepilecopy(av, y);
    2815                 :            : }
    2816                 :            : 
    2817                 :            : GEN
    2818                 :          0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
    2819                 :            : GEN
    2820                 :          0 : polredabs(GEN x) { return polredabs0(x,0); }
    2821                 :            : GEN
    2822                 :          0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
    2823                 :            : 
    2824                 :            : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
    2825                 :            :  * flag & nf_ORIG: + element (base change)
    2826                 :            :  * flag & nf_ABSOLUTE: absolute polynomial */
    2827                 :            : static GEN
    2828                 :         63 : rnfpolred_i(GEN nf, GEN relpol, long flag, long best)
    2829                 :            : {
    2830         [ +  + ]:         63 :   const char *f = best? "rnfpolredbest": "rnfpolredabs";
    2831 [ +  + ][ +  + ]:         63 :   const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
    2832                 :            :   pari_timer ti;
    2833                 :         63 :   GEN listP = NULL, red, bas, A, P, pol, T, rnfeq;
    2834                 :         63 :   long ty = typ(relpol);
    2835                 :         63 :   pari_sp av = avma;
    2836                 :            : 
    2837         [ +  + ]:         63 :   if (ty == t_VEC) {
    2838         [ -  + ]:         14 :     if (lg(relpol) != 3) pari_err_TYPE(f,relpol);
    2839                 :         14 :     listP = gel(relpol,2);
    2840                 :         14 :     relpol = gel(relpol,1);
    2841                 :            :   }
    2842         [ -  + ]:         63 :   if (typ(relpol) != t_POL) pari_err_TYPE(f,relpol);
    2843                 :         63 :   nf = checknf(nf);
    2844         [ -  + ]:         63 :   if (DEBUGLEVEL>1) timer_start(&ti);
    2845                 :         63 :   T = nf_get_pol(nf);
    2846                 :         63 :   relpol = RgX_nffix(f, T, relpol, 0);
    2847 [ +  + ][ +  + ]:         63 :   if (best || (flag & nf_PARTIALFACT))
    2848                 :            :   {
    2849         [ +  + ]:         42 :     if (abs)
    2850                 :            :     {
    2851                 :          7 :       rnfeq = nf_rnfeq(nf, relpol);
    2852                 :          7 :       pol = gel(rnfeq,1);
    2853                 :            :     }
    2854                 :            :     else
    2855                 :            :     {
    2856                 :            :       long sa;
    2857                 :         35 :       pol = rnfequationall(nf, relpol, &sa, NULL);
    2858                 :         35 :       rnfeq = mkvec5(gen_0,gen_0,stoi(sa),T,liftpol_shallow(relpol));
    2859                 :            :     }
    2860         [ +  + ]:         42 :     bas = listP? mkvec2(pol, listP): pol;
    2861         [ +  + ]:         84 :     if (best)
    2862                 :            :     {
    2863         [ +  + ]:         35 :       if (abs) red = polredbest(bas, 1);
    2864                 :            :       else
    2865                 :            :       {
    2866                 :            :         GEN ro, x, dx, a;
    2867                 :            :         nfbasic_t T;
    2868                 :         28 :         nfbasic_init(bas, nf_PARTIALFACT, &T);
    2869                 :         28 :         polredbest_aux(&T, &ro, &x, &dx, &a);
    2870                 :         28 :         red = mkvec2(x, a);
    2871                 :            :       }
    2872                 :            :     }
    2873                 :            :     else
    2874         [ -  + ]:          7 :       red = polredabs0(bas, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
    2875                 :            :   }
    2876                 :            :   else
    2877                 :            :   {
    2878                 :         21 :     GEN rnf = rnfinit(nf, relpol), M = rnf_basM(rnf);
    2879                 :         21 :     rnfeq = rnf_get_map(rnf);
    2880                 :         21 :     pol = rnf_get_polabs(rnf);
    2881                 :         21 :     bas = mkvec2(pol, RgM_to_RgXV(M, varn(pol)));
    2882         [ -  + ]:         21 :     if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
    2883                 :         21 :     red = polredabs0(bas, nf_RAW);
    2884                 :            :   }
    2885                 :         63 :   P = gel(red,1);
    2886                 :         63 :   A = gel(red,2);
    2887         [ -  + ]:         63 :   if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
    2888         [ +  + ]:         63 :   if (flag & nf_ABSOLUTE)
    2889                 :            :   {
    2890         [ +  + ]:         14 :     if (flag & nf_ORIG)
    2891                 :            :     {
    2892                 :          7 :       GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
    2893                 :          7 :       GEN k = gel(rnfeq,3); /* Mod(variable(relpol),relpol) + k*a root of pol */
    2894                 :          7 :       a = RgX_RgXQ_eval(a, lift_intern(A), P); /* Mod(a, P) root of T */
    2895                 :          7 :       P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
    2896                 :            :     }
    2897                 :         14 :     return gerepilecopy(av, P);
    2898                 :            :   }
    2899                 :         49 :   A = eltabstorel_lift(rnfeq, A);
    2900                 :         49 :   P = RgXQ_charpoly(A, relpol, varn(relpol));
    2901                 :         49 :   P = lift_if_rational(P);
    2902         [ +  + ]:         49 :   if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,relpol),P));
    2903                 :         63 :   return gerepilecopy(av, P);
    2904                 :            : }
    2905                 :            : GEN
    2906                 :         28 : rnfpolredabs(GEN nf, GEN relpol, long flag)
    2907                 :         28 : { return rnfpolred_i(nf,relpol,flag, 0); }
    2908                 :            : GEN
    2909                 :         35 : rnfpolredbest(GEN nf, GEN relpol, long flag)
    2910                 :            : {
    2911 [ +  - ][ -  + ]:         35 :   if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
    2912                 :         35 :   return rnfpolred_i(nf,relpol,flag, 1);
    2913                 :            : }

Generated by: LCOV version 1.9