Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 18616-fe48336) Lines: 1541 1651 93.3 %
Date: 2016-02-14 Functions: 115 130 88.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 984 1219 80.7 %

           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                 :     213703 : checkrnf_i(GEN rnf)
      26 [ +  + ][ +  + ]:     213703 : { return (typ(rnf)==t_VEC && lg(rnf)==13); }
      27                 :            : 
      28                 :            : void
      29                 :     208194 : checkrnf(GEN rnf)
      30         [ +  + ]:     208194 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
      31                 :            : 
      32                 :            : GEN
      33                 :     680773 : checkbnf_i(GEN X)
      34                 :            : {
      35         [ +  + ]:     680773 :   if (typ(X) == t_VEC)
      36      [ +  +  + ]:     680360 :     switch (lg(X))
      37                 :            :     {
      38                 :     678890 :       case 11: return X;
      39                 :        994 :       case 7:  return checkbnf_i(bnr_get_bnf(X));
      40                 :            :     }
      41                 :     680773 :   return NULL;
      42                 :            : }
      43                 :            : 
      44                 :            : GEN
      45                 :    5786318 : checknf_i(GEN X)
      46                 :            : {
      47         [ +  + ]:    5786318 :   if (typ(X)==t_VEC)
      48   [ +  +  +  +  :    5785898 :     switch(lg(X))
                      + ]
      49                 :            :     {
      50                 :    5571197 :       case 10: return X;
      51                 :     208674 :       case 11: return checknf_i(bnf_get_nf(X));
      52                 :       5565 :       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                 :    5786318 :   return NULL;
      56                 :            : }
      57                 :            : 
      58                 :            : GEN
      59                 :     678757 : checkbnf(GEN x)
      60                 :            : {
      61                 :     678757 :   GEN bnf = checkbnf_i(x);
      62         [ -  + ]:     678757 :   if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
      63                 :     678757 :   return bnf;
      64                 :            : }
      65                 :            : 
      66                 :            : GEN
      67                 :    5570609 : checknf(GEN x)
      68                 :            : {
      69                 :    5570609 :   GEN nf = checknf_i(x);
      70         [ +  + ]:    5570609 :   if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
      71                 :    5570595 :   return nf;
      72                 :            : }
      73                 :            : 
      74                 :            : void
      75                 :     195368 : checkbnr(GEN bnr)
      76                 :            : {
      77 [ +  - ][ -  + ]:     195368 :   if (typ(bnr)!=t_VEC || lg(bnr)!=7)
      78                 :          0 :     pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
      79                 :     195368 :   (void)checkbnf(bnr_get_bnf(bnr));
      80                 :     195368 : }
      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                 :     204321 : checkbid_i(GEN bid)
      99                 :            : {
     100                 :            :   GEN f;
     101 [ +  + ][ +  + ]:     204321 :   if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(gel(bid,3)) != t_MAT) return NULL;
                 [ +  + ]
     102                 :     198343 :   f = bid_get_mod(bid);
     103 [ +  - ][ -  + ]:     198343 :   if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
     104                 :     204321 :   return bid;
     105                 :            : }
     106                 :            : void
     107                 :     198343 : checkbid(GEN bid)
     108                 :            : {
     109         [ +  + ]:     198343 :   if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
     110                 :     198336 : }
     111                 :            : void
     112                 :      11389 : checkabgrp(GEN v)
     113                 :            : {
     114         [ +  - ]:      11389 :   if (typ(v) == t_VEC) switch(lg(v))
              [ +  +  - ]
     115                 :            :   {
     116         [ -  + ]:      11298 :     case 4: if (typ(gel(v,3)) != t_VEC) break;
     117         [ +  + ]:      11389 :     case 3: if (typ(gel(v,2)) != t_VEC) break;
     118         [ -  + ]:      11361 :             if (typ(gel(v,1)) != t_INT) break;
     119                 :      11361 :             return;/*OK*/
     120                 :          0 :     default: break;
     121                 :            :   }
     122                 :         28 :   pari_err_TYPE("checkabgrp",v);
     123                 :            : }
     124                 :            : 
     125                 :            : GEN
     126                 :      82726 : checknfelt_mod(GEN nf, GEN x, const char *s)
     127                 :            : {
     128                 :      82726 :   GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
     129         [ +  + ]:      82726 :   if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
     130                 :      82656 :   return a;
     131                 :            : }
     132                 :            : 
     133                 :            : void
     134                 :       3038 : check_ZKmodule(GEN x, const char *s)
     135                 :            : {
     136 [ +  - ][ -  + ]:       3038 :   if (typ(x) != t_VEC || lg(x) < 3) pari_err_TYPE(s,x);
     137         [ -  + ]:       3038 :   if (typ(gel(x,1)) != t_MAT) pari_err_TYPE(s,gel(x,1));
     138         [ -  + ]:       3038 :   if (typ(gel(x,2)) != t_VEC) pari_err_TYPE(s,gel(x,2));
     139         [ -  + ]:       3038 :   if (lg(gel(x,2)) != lgcols(x)) pari_err_DIM(s);
     140                 :       3038 : }
     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                 :      13027 : get_bnf(GEN x, long *t)
     156                 :            : {
     157   [ +  +  +  +  :      13027 :   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   [ +  +  +  +  :      12411 :       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                 :        175 :         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                 :      13027 :   *t = typ_NULL; return NULL;
     184                 :            : }
     185                 :            : 
     186                 :            : GEN
     187                 :      62083 : get_nf(GEN x, long *t)
     188                 :            : {
     189   [ +  +  +  +  :      62083 :   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   [ +  +  +  +  :      59668 :       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                 :       4641 :         case 11: *t = typ_BNF;
     211 [ +  - ][ +  - ]:       4641 :           x = bnf_get_nf(x); if (typ(x)!=t_VEC || lg(x)!=10) break;
     212                 :       4641 :           return x;
     213                 :        343 :         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                 :      62083 :   *t = typ_NULL; return NULL;
     222                 :            : }
     223                 :            : 
     224                 :            : long
     225                 :      38129 : nftyp(GEN x)
     226                 :            : {
     227   [ +  +  +  + ]:      38129 :   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   [ +  +  +  +  :      38101 :       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 [ +  + ][ +  - ]:        770 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     243 [ +  - ][ +  - ]:        763 :           x = bnf_get_nf(x);  if (typ(x)!=t_VEC || lg(x)!=10) break;
     244                 :        763 :           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                 :      38129 :   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                 :      12537 : ZX_Z_normalize(GEN pol, GEN *ptk)
     296                 :            : {
     297                 :      12537 :   long i,j, sk, n = degpol(pol); /* > 0 */
     298                 :            :   GEN k, fa, P, E, a, POL;
     299                 :            : 
     300                 :      12537 :   a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
     301         [ +  + ]:      46354 :   for (i = n-2; i >= 0; i--)
     302                 :            :   {
     303                 :      41230 :     k = gcdii(k, gel(a,i));
     304 [ +  + ][ +  + ]:      41230 :     if (is_pm1(k)) { if (ptk) *ptk = gen_1; return pol; }
     305                 :            :   }
     306                 :       5124 :   sk = signe(k);
     307 [ +  + ][ +  - ]:       5124 :   if (!sk) { if (ptk) *ptk = gen_1; return pol; /* monomial! */ }
     308                 :       4298 :   fa = absi_factor_limit(k, 0); k = gen_1;
     309                 :       4298 :   P = gel(fa,1);
     310                 :       4298 :   E = gel(fa,2);
     311                 :       4298 :   POL = leafcopy(pol); a = POL+2;
     312         [ +  + ]:       9779 :   for (i = lg(P)-1; i > 0; i--)
     313                 :            :   {
     314                 :       5481 :     GEN p = gel(P,i), pv, pvj;
     315                 :       5481 :     long vmin = itos(gel(E,i));
     316                 :            :     /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
     317         [ +  + ]:      37646 :     for (j=n-1; j>=0; j--)
     318                 :            :     {
     319                 :            :       long v;
     320         [ +  + ]:      32165 :       if (!signe(gel(a,j))) continue;
     321                 :      20146 :       v = Z_pval(gel(a,j), p) / (n - j);
     322         [ +  + ]:      20146 :       if (v < vmin) vmin = v;
     323                 :            :     }
     324         [ +  + ]:       5481 :     if (!vmin) continue;
     325                 :       1015 :     pvj = pv = powiu(p,vmin); k = mulii(k, pv);
     326                 :            :     /* a[j] /= p^(v*(n-j)) */
     327         [ +  + ]:       7294 :     for (j=n-1; j>=0; j--)
     328                 :            :     {
     329         [ +  + ]:       6279 :       if (j < n-1) pvj = mulii(pvj, pv);
     330                 :       6279 :       gel(a,j) = diviiexact(gel(a,j), pvj);
     331                 :            :     }
     332                 :            :   }
     333         [ +  + ]:      12537 :   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                 :      12544 : ZX_primitive_to_monic(GEN pol, GEN *pL)
     341                 :            : {
     342                 :      12544 :   long i,j, n = degpol(pol);
     343                 :      12544 :   GEN lc = leading_term(pol), L, fa, P, E, a, POL;
     344                 :            : 
     345         [ +  + ]:      12544 :   if (is_pm1(lc))
     346                 :            :   {
     347         [ +  - ]:      12299 :     if (pL) *pL = gen_1;
     348         [ +  + ]:      12299 :     return signe(lc) < 0? ZX_neg(pol): pol;
     349                 :            :   }
     350         [ +  + ]:        245 :   if (signe(lc) < 0)
     351                 :         35 :     POL = ZX_neg(pol);
     352                 :            :   else
     353                 :        210 :     POL = leafcopy(pol);
     354                 :        245 :   a = POL+2; lc = gel(a,n);
     355                 :        245 :   fa = Z_factor_limit(lc,0); L = gen_1;
     356                 :        245 :   P = gel(fa,1);
     357                 :        245 :   E = gel(fa,2);
     358         [ +  + ]:        623 :   for (i = lg(P)-1; i > 0; i--)
     359                 :            :   {
     360                 :        378 :     GEN p = gel(P,i), pk, pku;
     361                 :        378 :     long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
     362                 :            : 
     363         [ +  + ]:        378 :     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         [ +  + ]:       1498 :     for (j=n-1; j>0; j--)
     366                 :            :     {
     367         [ +  + ]:       1120 :       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                 :        378 :     pk = powiu(p,k); j0 = d/k;
     372                 :        378 :     L = mulii(L, pk);
     373                 :            : 
     374                 :        378 :     pku = powiu(p,d - k*j0);
     375                 :            :     /* a[j] *= p^(d - kj) */
     376         [ +  + ]:       1582 :     for (j=j0; j>=0; j--)
     377                 :            :     {
     378         [ +  + ]:       1204 :       if (j < j0) pku = mulii(pku, pk);
     379                 :       1204 :       gel(a,j) = mulii(gel(a,j), pku);
     380                 :            :     }
     381                 :        378 :     j0++;
     382                 :        378 :     pku = powiu(p,k*j0 - d);
     383                 :            :     /* a[j] /= p^(kj - d) */
     384         [ +  + ]:       1050 :     for (j=j0; j<=n; j++)
     385                 :            :     {
     386         [ +  + ]:        672 :       if (j > j0) pku = mulii(pku, pk);
     387                 :        672 :       gel(a,j) = diviiexact(gel(a,j), pku);
     388                 :            :     }
     389                 :            :   }
     390         [ +  - ]:        245 :   if (pL) *pL = L;
     391                 :      12544 :   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                 :      12285 : ZX_Q_normalize(GEN pol, GEN *pL)
     399                 :            : {
     400                 :      12285 :   GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
     401                 :      12285 :   POL = ZX_Z_normalize(POL, pL);
     402         [ +  + ]:      12285 :   if (pL) *pL = gdiv(lc, *pL);
     403                 :      12285 :   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                 :       1386 : roots_to_ZX(GEN z, long *e)
     428                 :            : {
     429                 :       1386 :   GEN a = roots_to_pol(z,0);
     430                 :       1386 :   GEN b = grndtoi(real_i(a),e);
     431                 :       1386 :   long e1 = gexpo(imag_i(a));
     432         [ +  + ]:       1386 :   if (e1 > *e) *e = e1;
     433                 :       1386 :   return b;
     434                 :            : }
     435                 :            : 
     436                 :            : static GEN
     437                 :       1050 : polgaloisnames(long a, long b)
     438                 :            : {
     439                 :       1050 :   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                 :       1050 :    const long idx[]={0,1,2,4,9,14,30};
     453                 :       1050 :    return strtoGENstr(t[idx[a-1]+b-1]);
     454                 :            : }
     455                 :            : 
     456                 :            : static GEN
     457                 :       1050 : galois_res(long d, long n, long s, long k)
     458                 :            : {
     459                 :       1050 :   long kk = k;
     460                 :       1050 :   GEN z = cgetg(5,t_VEC);
     461         [ +  + ]:       1050 :   if (!new_galois_format)
     462                 :            :   {
     463         [ +  + ]:        546 :     switch (d) {
     464                 :            :       case 6:
     465 [ +  + ][ +  + ]:        238 :         kk = (k == 6 || k == 2)? 2: 1;
     466                 :        238 :         break;
     467                 :            :       default:
     468                 :        308 :         kk = 1;
     469                 :            :     }
     470                 :            :   }
     471                 :       1050 :   gel(z,1) = stoi(n);
     472                 :       1050 :   gel(z,2) = stoi(s);
     473                 :       1050 :   gel(z,3) = stoi(kk);
     474                 :       1050 :   gel(z,4) = polgaloisnames(d,k);
     475                 :       1050 :   return z;
     476                 :            : }
     477                 :            : 
     478                 :            : GEN
     479                 :       3052 : polgalois(GEN x, long prec)
     480                 :            : {
     481                 :       3052 :   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                 :       3052 :   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                 :       3052 :   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         [ -  + ]:       3052 :   if (typ(x)!=t_POL) pari_err_TYPE("galois",x);
     489                 :       3052 :   n=degpol(x);
     490         [ -  + ]:       3052 :   if (n>11) pari_err_IMPL("galois of degree higher than 11");
     491                 :       3052 :   x = Q_primpart(x);
     492                 :       3052 :   RgX_check_ZX(x, "galois");
     493         [ -  + ]:       3052 :   if (!ZX_is_irred(x)) pari_err_IRREDPOL("galois",x);
     494                 :            : 
     495         [ +  + ]:       3052 :   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                 :       2961 :   x1 = x = ZX_Q_normalize(x,NULL); av1=avma;
     506         [ +  + ]:       2961 :   if (n > 7) return galoisbig(x, prec);
     507                 :            :   for(;;)
     508                 :            :   {
     509                 :       1176 :     double fb = fujiwara_bound(x);
     510   [ +  +  +  +  :       1176 :     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                 :        196 :       case 7: z = cgetg(36,t_VEC);
     691                 :        196 :         prec = nbits2prec((long)(fb*7.) + 64);
     692                 :            :         for(;;)
     693                 :            :         {
     694                 :        294 :           ind = 0; p1=QX_complex_roots(x,prec);
     695         [ +  + ]:       1764 :           for (i=1; i<=5; i++)
     696         [ +  + ]:       5880 :             for (j=i+1; j<=6; j++)
     697                 :            :             {
     698                 :       4410 :               GEN t = gadd(gel(p1,i),gel(p1,j));
     699         [ +  + ]:      14700 :               for (k=j+1; k<=7; k++) gel(z,++ind) = gadd(t, gel(p1,k));
     700                 :            :             }
     701         [ +  + ]:        294 :           p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     702                 :         98 :           prec = precdbl(prec);
     703                 :         98 :         }
     704         [ -  + ]:        196 :         if (!ZX_is_squarefree(p5)) goto tchi;
     705                 :        196 :         p2=gel(ZX_factor(p5),1);
     706   [ +  +  +  +  :        196 :         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                 :         56 :           case 2: f = degpol(gel(p2,1)); avma = av;
     711 [ +  + ][ -  + ]:         56 :             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                 :       3269 :     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                 :      19502 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
     740                 :            : {
     741                 :      19502 :   pari_sp av = avma;
     742                 :      19502 :   long i = lg(pol)-1;
     743                 :            :   GEN res, den;
     744         [ -  + ]:      19502 :   if (i==1) return gen_0;
     745                 :      19502 :   pol = Q_remove_denom(pol, &den);
     746                 :      19502 :   res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
     747         [ +  + ]:     104902 :   for (i-- ; i>=2; i--)
     748                 :      85400 :     res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
     749         [ +  + ]:      19502 :   if (den) res = RgC_Rg_div(res, den);
     750                 :      19502 :   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                 :      10059 : ZC_galoisapply(GEN nf, GEN s, GEN x)
     782                 :            : {
     783                 :      10059 :   x = nf_to_scalar_or_alg(nf, x);
     784         [ +  + ]:      10059 :   if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
     785                 :      10059 :   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                 :       5026 : galoisapply(GEN nf, GEN aut, GEN x)
     855                 :            : {
     856                 :       5026 :   pari_sp av = avma;
     857                 :            :   long lx, j;
     858                 :            :   GEN y;
     859                 :            : 
     860                 :       5026 :   nf = checknf(nf);
     861   [ +  +  +  +  :       5026 :   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                 :       3514 :       aut = algtobasis(nf, aut);
     886                 :       3514 :       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                 :       5026 :   return NULL; /* not reached */
     898                 :            : }
     899                 :            : 
     900                 :            : GEN
     901                 :       1876 : nfgaloismatrix(GEN nf, GEN s)
     902                 :            : {
     903                 :            :   GEN zk, M, m;
     904                 :            :   long k, l;
     905                 :       1876 :   nf = checknf(nf);
     906                 :       1876 :   zk = nf_get_zk(nf);
     907         [ -  + ]:       1876 :   if (typ(s) != t_COL) s = algtobasis(nf, s); /* left on stack for efficiency */
     908                 :       1876 :   m = zk_multable(nf, s);
     909                 :       1876 :   l = lg(s); M = cgetg(l, t_MAT);
     910                 :       1876 :   gel(M, 1) = col_ei(l-1, 1); /* s(1) = 1 */
     911         [ +  + ]:       8575 :   for (k = 2; k < l; k++)
     912                 :       6699 :     gel(M, k) = QX_table_nfpoleval(nf, gel(zk, k), m);
     913                 :       1876 :   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                 :      15257 : get_nfpol(GEN x, GEN *nf)
    1159                 :            : {
    1160         [ +  + ]:      15257 :   if (typ(x) == t_POL) { *nf = NULL; return x; }
    1161                 :      15257 :   *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                 :       9560 : get_roots(GEN x, long r1, long prec)
    1280                 :            : {
    1281                 :            :   long i, ru;
    1282                 :            :   GEN z;
    1283         [ -  + ]:       9560 :   if (typ(x) != t_POL)
    1284                 :            :   {
    1285                 :          0 :     z = leafcopy(x);
    1286                 :          0 :     ru = (lg(z)-1 + r1) >> 1;
    1287                 :            :   }
    1288                 :            :   else
    1289                 :            :   {
    1290                 :       9560 :     long n = degpol(x);
    1291         [ +  + ]:       9560 :     z = (r1 == n)? realroots(x, NULL, prec): QX_complex_roots(x,prec);
    1292                 :       9560 :     ru = (n+r1)>>1;
    1293                 :            :   }
    1294         [ +  + ]:      23803 :   for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
    1295                 :       9560 :   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                 :      67039 : quicktrace(GEN x, GEN sym)
    1307                 :            : {
    1308                 :      67039 :   GEN p1 = gen_0;
    1309                 :            :   long i;
    1310                 :            : 
    1311         [ -  + ]:      67039 :   if (typ(x) != t_POL) return gmul(x, gel(sym,1));
    1312         [ +  - ]:      67039 :   if (signe(x))
    1313                 :            :   {
    1314                 :      67039 :     sym--;
    1315         [ +  + ]:    1005123 :     for (i=lg(x)-1; i>1; i--)
    1316                 :     938084 :       p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
    1317                 :            :   }
    1318                 :      67039 :   return p1;
    1319                 :            : }
    1320                 :            : 
    1321                 :            : static GEN
    1322                 :       5495 : get_Tr(GEN mul, GEN x, GEN basden)
    1323                 :            : {
    1324                 :       5495 :   GEN t, bas = gel(basden,1), den = gel(basden,2);
    1325                 :       5495 :   long i, j, n = lg(bas)-1;
    1326                 :       5495 :   GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
    1327                 :            : 
    1328                 :       5495 :   gel(TW,1) = utoipos(n);
    1329         [ +  + ]:      21910 :   for (i=2; i<=n; i++)
    1330                 :            :   {
    1331                 :      16415 :     t = quicktrace(gel(bas,i), sym);
    1332 [ +  + ][ +  + ]:      16415 :     if (den && den[i]) t = diviiexact(t,gel(den,i));
    1333                 :      16415 :     gel(TW,i) = t; /* tr(w[i]) */
    1334                 :            :   }
    1335                 :       5495 :   gel(T,1) = TW;
    1336         [ +  + ]:      21910 :   for (i=2; i<=n; i++)
    1337                 :            :   {
    1338                 :      16415 :     gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
    1339         [ +  + ]:     129206 :     for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
    1340                 :     112791 :       gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
    1341                 :            :   }
    1342                 :       5495 :   return T;
    1343                 :            : }
    1344                 :            : 
    1345                 :            : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
    1346                 :            : GEN
    1347                 :      13849 : get_bas_den(GEN bas)
    1348                 :            : {
    1349                 :      13849 :   GEN b,d,den, dbas = leafcopy(bas);
    1350                 :      13849 :   long i, l = lg(bas);
    1351                 :      13849 :   int power = 1;
    1352                 :      13849 :   den = cgetg(l,t_VEC);
    1353         [ +  + ]:      73046 :   for (i=1; i<l; i++)
    1354                 :            :   {
    1355                 :      59197 :     b = Q_remove_denom(gel(bas,i), &d);
    1356                 :      59197 :     gel(dbas,i) = b;
    1357         [ +  + ]:      59197 :     gel(den,i) = d; if (d) power = 0;
    1358                 :            :   }
    1359         [ +  + ]:      13849 :   if (power) den = NULL; /* power basis */
    1360                 :      13849 :   return mkvec2(dbas, den);
    1361                 :            : }
    1362                 :            : 
    1363                 :            : /* Internal: nf partially filled. Require pol; fill zk, invzk, multable */
    1364                 :            : void
    1365                 :       5495 : nf_set_multable(GEN nf, GEN bas, GEN basden)
    1366                 :            : {
    1367                 :       5495 :   GEN T = nf_get_pol(nf), invbas, basM;
    1368                 :       5495 :   long i,j, n = degpol(T);
    1369                 :       5495 :   GEN w, den, mul = cgetg(n*n+1,t_MAT);
    1370                 :            : 
    1371         [ -  + ]:       5495 :   if (typ(bas) == t_MAT)
    1372                 :          0 :   { basM = bas; bas = RgM_to_RgXV(basM, varn(T)); }
    1373                 :            :   else
    1374                 :       5495 :     basM = RgV_to_RgM(bas, n);
    1375                 :       5495 :   gel(nf,7) = bas;
    1376                 :       5495 :   gel(nf,8) = invbas = QM_inv(basM, gen_1);
    1377                 :       5495 :   gel(nf,9) = mul;
    1378                 :            : 
    1379         [ -  + ]:       5495 :   if (!basden) basden = get_bas_den(nf_get_zk(nf)); /*integral basis*/
    1380                 :       5495 :   w   = gel(basden,1);
    1381                 :       5495 :   den = gel(basden,2);
    1382                 :            :   /* i = 1 split for efficiency, assume w[1] = 1 */
    1383         [ +  + ]:      27405 :   for (j=1; j<=n; j++)
    1384                 :      21910 :     gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
    1385         [ +  + ]:      21910 :   for (i=2; i<=n; i++)
    1386         [ +  + ]:     129206 :     for (j=i; j<=n; j++)
    1387                 :            :     {
    1388                 :     112791 :       pari_sp av = avma;
    1389         [ +  + ]:     112791 :       GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
    1390                 :     112791 :       z = mulmat_pol(invbas, z); /* integral column */
    1391         [ +  + ]:     112791 :       if (den)
    1392                 :            :       {
    1393                 :      79121 :         GEN d = mul_denom(gel(den,i), gel(den,j));
    1394         [ +  + ]:      79121 :         if (d) z = ZC_Z_divexact(z, d);
    1395                 :            :       }
    1396                 :     112791 :       gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
    1397                 :            :     }
    1398                 :       5495 : }
    1399                 :            : 
    1400                 :            : /* as get_Tr, mul_table not precomputed */
    1401                 :            : static GEN
    1402                 :       2723 : make_Tr(GEN x, GEN basden)
    1403                 :            : {
    1404                 :       2723 :   long i,j, n = degpol(x);
    1405                 :            :   GEN c, t, d;
    1406                 :       2723 :   GEN T   = cgetg(n+1,t_MAT);
    1407                 :       2723 :   GEN sym = polsym(x, n-1);
    1408                 :       2723 :   GEN w   = gel(basden,1); /* W[i] = w[i]/den[i] */
    1409                 :       2723 :   GEN den = gel(basden,2);
    1410                 :            :   /* assume W[1] = 1, case i = 1 split for efficiency */
    1411                 :       2723 :   c = cgetg(n+1,t_COL); gel(T,1) = c;
    1412                 :       2723 :   gel(c, 1) = utoipos(n);
    1413         [ +  + ]:       8883 :   for (j=2; j<=n; j++)
    1414                 :            :   {
    1415                 :       6160 :     pari_sp av = avma;
    1416                 :       6160 :     t = quicktrace(gel(w,j), sym);
    1417         [ +  + ]:       6160 :     if (den)
    1418                 :            :     {
    1419                 :       4249 :       d = gel(den,j);
    1420         [ +  + ]:       4249 :       if (d) t = diviiexact(t, d);
    1421                 :            :     }
    1422                 :       6160 :     gel(c,j) = gerepileuptoint(av, t);
    1423                 :            :   }
    1424         [ +  + ]:       8883 :   for (i=2; i<=n; i++)
    1425                 :            :   {
    1426                 :       6160 :     c = cgetg(n+1,t_COL); gel(T,i) = c;
    1427         [ +  + ]:      42714 :     for (j=1; j<i ; j++) gel(c,j) = gcoeff(T,i,j);
    1428         [ +  + ]:      42714 :     for (   ; j<=n; j++)
    1429                 :            :     {
    1430                 :      36554 :       pari_sp av = avma;
    1431         [ +  + ]:      36554 :       t = (i == j)? ZXQ_sqr(gel(w,i), x): ZXQ_mul(gel(w,i),gel(w,j), x);
    1432                 :      36554 :       t = quicktrace(t, sym);
    1433         [ +  + ]:      36554 :       if (den)
    1434                 :            :       {
    1435                 :      33040 :         d = mul_denom(gel(den,i),gel(den,j));
    1436         [ +  + ]:      33040 :         if (d) t = diviiexact(t, d);
    1437                 :            :       }
    1438                 :      36554 :       gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
    1439                 :            :     }
    1440                 :            :   }
    1441                 :       2723 :   return T;
    1442                 :            : }
    1443                 :            : 
    1444                 :            : /* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */
    1445                 :            : static void
    1446                 :      11160 : get_roots_for_M(nffp_t *F)
    1447                 :            : {
    1448                 :            :   long n, eBD, PREC;
    1449                 :            : 
    1450         [ +  - ]:      11160 :   if (F->extraprec < 0)
    1451                 :            :   { /* not initialized yet */
    1452                 :            :     double er;
    1453                 :      11160 :     n = degpol(F->x);
    1454                 :      11160 :     eBD = 1 + gexpo(gel(F->basden,1));
    1455         [ +  + ]:      11160 :     er  = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->x);
    1456         [ +  + ]:      11160 :     if (er < 0) er = 0;
    1457                 :      11160 :     F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
    1458                 :            :   }
    1459                 :            : 
    1460                 :      11160 :   PREC = F->prec + F->extraprec;
    1461                 :            :   /* make sure that default accuracy is the same on 32/64bit */
    1462                 :            : #ifndef LONG_IS_64BIT
    1463         [ +  + ]:       1632 :   if (odd(PREC)) PREC += EXTRAPRECWORD;
    1464                 :            : #endif
    1465 [ +  + ][ +  + ]:      20720 :   if (F->ro && gprecision(gel(F->ro,1)) >= PREC) return;
    1466                 :       9560 :   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                 :      11160 : make_M(nffp_t *F, int trunc)
    1472                 :            : {
    1473                 :      11160 :   GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
    1474                 :            :   GEN m, d, M;
    1475                 :      11160 :   long i, j, l = lg(ro), n = lg(bas);
    1476                 :      11160 :   M = cgetg(n,t_MAT);
    1477                 :      11160 :   gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
    1478         [ +  + ]:      50410 :   for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
    1479         [ +  + ]:      43187 :   for (i=1; i<l; i++)
    1480                 :            :   {
    1481                 :      32027 :     GEN r = gel(ro,i), ri;
    1482         [ +  + ]:      32027 :     ri = (gexpo(r) > 1)? ginv(r): NULL;
    1483         [ +  + ]:     330402 :     for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
    1484                 :            :   }
    1485         [ +  + ]:      11160 :   if (den)
    1486         [ +  + ]:      31173 :     for (j=2; j<n; j++)
    1487                 :            :     {
    1488         [ +  + ]:      26285 :       d = gel(den,j); if (!d) continue;
    1489                 :      21055 :       m = gel(M,j);
    1490         [ +  + ]:     206258 :       for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
    1491                 :            :     }
    1492                 :            : 
    1493 [ +  + ][ +  + ]:      11160 :   if (trunc && gprecision(M) > F->prec)
    1494                 :            :   {
    1495                 :       1843 :     M     = gprec_w(M, F->prec);
    1496                 :       1843 :     F->ro = gprec_w(ro,F->prec);
    1497                 :            :   }
    1498                 :      11160 :   F->M = M;
    1499                 :      11160 : }
    1500                 :            : 
    1501                 :            : /* return G real such that G~ * G = T_2 */
    1502                 :            : static void
    1503                 :      11160 : make_G(nffp_t *F)
    1504                 :            : {
    1505                 :      11160 :   GEN G, M = F->M;
    1506                 :      11160 :   long i, j, k, r1 = F->r1, l = lg(M);
    1507                 :            : 
    1508                 :      11160 :   G = cgetg(l, t_MAT);
    1509         [ +  + ]:      61570 :   for (j=1; j<l; j++)
    1510                 :            :   {
    1511                 :      50410 :     GEN g = cgetg(l, t_COL);
    1512                 :      50410 :     GEN m = gel(M,j);
    1513                 :      50410 :     gel(G,j) = g;
    1514         [ +  + ]:     153718 :     for (k=i=1; i<=r1; i++) g[k++] = m[i];
    1515         [ +  + ]:     277504 :     for (     ; k < l; i++)
    1516                 :            :     {
    1517                 :     227094 :       GEN r = gel(m,i);
    1518         [ +  + ]:     227094 :       if (typ(r) == t_COMPLEX)
    1519                 :            :       {
    1520                 :     208711 :         gel(g,k++) = mpadd(gel(r,1), gel(r,2));
    1521                 :     208711 :         gel(g,k++) = mpsub(gel(r,1), gel(r,2));
    1522                 :            :       }
    1523                 :            :       else
    1524                 :            :       {
    1525                 :      18383 :         gel(g,k++) = r;
    1526                 :      18383 :         gel(g,k++) = r;
    1527                 :            :       }
    1528                 :            :     }
    1529                 :            :   }
    1530                 :      11160 :   F->G = G;
    1531                 :      11160 : }
    1532                 :            : 
    1533                 :            : static void
    1534                 :      11160 : make_M_G(nffp_t *F, int trunc)
    1535                 :            : {
    1536                 :      11160 :   get_roots_for_M(F);
    1537                 :      11160 :   make_M(F, trunc);
    1538                 :      11160 :   make_G(F);
    1539                 :      11160 : }
    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                 :      10408 : nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
    1554                 :            : {
    1555                 :      10408 :   F->x  = T->x;
    1556                 :      10408 :   F->ro = ro;
    1557                 :      10408 :   F->r1 = T->r1;
    1558         [ +  + ]:      10408 :   if (!T->basden) T->basden = get_bas_den(T->bas);
    1559                 :      10408 :   F->basden = T->basden;
    1560                 :      10408 :   F->extraprec = -1;
    1561                 :      10408 :   F->prec = prec;
    1562                 :      10408 : }
    1563                 :            : 
    1564                 :            : static void
    1565                 :       6614 : get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, int trunc, long prec)
    1566                 :            : {
    1567                 :       6614 :   nffp_init(F,T,ro,prec);
    1568                 :       6614 :   make_M_G(F, trunc);
    1569                 :       6614 : }
    1570                 :            : 
    1571                 :            : static GEN
    1572                 :       5495 : get_sign(long r1, long n) { return mkvec2s(r1, (n-r1)>>1); }
    1573                 :            : 
    1574                 :            : GEN
    1575                 :       5495 : nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec)
    1576                 :            : {
    1577                 :       5495 :   GEN nf = cgetg(10,t_VEC);
    1578                 :       5495 :   GEN x = T->x, absdK, Tr, D, TI, A, dA, MDI, mat = cgetg(9,t_VEC);
    1579                 :       5495 :   long n = degpol(T->x);
    1580                 :            :   nffp_t F;
    1581                 :       5495 :   get_nf_fp_compo(T, &F, ro, 0, prec);
    1582                 :            : 
    1583                 :       5495 :   gel(nf,1) = T->x;
    1584                 :       5495 :   gel(nf,2) = get_sign(T->r1, n);
    1585                 :       5495 :   gel(nf,3) = T->dK;
    1586                 :       5495 :   gel(nf,4) = T->index;
    1587                 :       5495 :   gel(nf,6) = F.ro;
    1588                 :       5495 :   gel(nf,5) = mat;
    1589                 :            : 
    1590                 :       5495 :   gel(mat,1) = F.M;
    1591                 :       5495 :   gel(mat,2) = F.G;
    1592                 :            : 
    1593                 :       5495 :   nf_set_multable(nf, T->bas, F.basden);
    1594                 :            : 
    1595                 :       5495 :   Tr = get_Tr(gel(nf,9), x, F.basden);
    1596         [ +  + ]:       5495 :   absdK = T->dK; if (signe(absdK) < 0) absdK = negi(absdK);
    1597                 :       5495 :   TI = ZM_inv(Tr, absdK); /* dK T^-1 */
    1598                 :       5495 :   A = Q_primitive_part(TI, &dA);
    1599                 :       5495 :   gel(mat,6) = A; /* primitive part of codifferent, dA its content */
    1600         [ +  + ]:       5495 :   dA = dA? diviiexact(absdK, dA): absdK;
    1601                 :       5495 :   A = ZM_hnfmodid(A, dA);
    1602                 :       5495 :   MDI = idealtwoelt(nf, A);
    1603                 :       5495 :   gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
    1604                 :       5495 :   gel(mat,7) = MDI;
    1605         [ +  + ]:       5495 :   if (is_pm1(T->index)) /* principal ideal (x'), whose norm is |dK| */
    1606                 :            :   {
    1607                 :       3675 :     D = zk_scalar_or_multable(nf, ZX_deriv(x));
    1608         [ +  + ]:       3675 :     if (typ(D) == t_MAT) D = ZM_hnfmod(D, absdK);
    1609                 :            :   }
    1610                 :            :   else
    1611                 :       1820 :     D = RgM_Rg_mul(idealinv(nf, A), dA);
    1612                 :       5495 :   gel(mat,3) = RM_round_maxrank(F.G);
    1613                 :       5495 :   gel(mat,4) = Tr;
    1614                 :       5495 :   gel(mat,5) = D;
    1615         [ +  + ]:       5495 :   gel(mat,8) = T->dKP? shallowtrans(T->dKP): cgetg(1,t_VEC);
    1616                 :       5495 :   return nf;
    1617                 :            : }
    1618                 :            : 
    1619                 :            : static GEN
    1620                 :         84 : primes_certify(GEN dK, GEN dKP)
    1621                 :            : {
    1622                 :         84 :   pari_sp av = avma;
    1623                 :         84 :   long i, l = lg(dKP);
    1624                 :         84 :   GEN v, D = dK;
    1625                 :         84 :   v = vectrunc_init(l);
    1626         [ +  + ]:        413 :   for (i = 1; i < l; i++)
    1627                 :            :   {
    1628                 :        329 :     GEN p = gel(dKP,i);
    1629         [ +  + ]:        329 :     if (!isprime(p)) vectrunc_append(v, p);
    1630                 :        329 :     (void)Z_pvalrem(D, p, &D);
    1631                 :            :   }
    1632         [ -  + ]:         84 :   if (!is_pm1(D))
    1633                 :            :   {
    1634         [ #  # ]:          0 :     if (signe(D) < 0) D = negi(D);
    1635         [ #  # ]:          0 :     if (!isprime(D)) vectrunc_append(v, D);
    1636                 :            :   }
    1637                 :         84 :   fixlg(v, lg(v)); return gerepilecopy(av, v);
    1638                 :            : }
    1639                 :            : GEN
    1640                 :          7 : nfcertify(GEN nf)
    1641                 :            : {
    1642                 :          7 :   nf = checknf(nf);
    1643                 :          7 :   return primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
    1644                 :            : }
    1645                 :            : 
    1646                 :            : #if 0 /* used to check benches between HNF nf.zk and LLL-reduced nf.zk */
    1647                 :            : static GEN
    1648                 :            : hnffromLLL(GEN nf)
    1649                 :            : {
    1650                 :            :   GEN d, x;
    1651                 :            :   x = RgV_to_RgM(nf_get_zk(nf), nf_get_degree(nf));
    1652                 :            :   x = Q_remove_denom(x, &d);
    1653                 :            :   if (!d) return x; /* power basis */
    1654                 :            :   return RgM_solve(ZM_hnfmodid(x, d), x);
    1655                 :            : }
    1656                 :            : 
    1657                 :            : static GEN
    1658                 :            : nfbasechange(GEN u, GEN x)
    1659                 :            : {
    1660                 :            :   long i,lx;
    1661                 :            :   GEN y;
    1662                 :            :   switch(typ(x))
    1663                 :            :   {
    1664                 :            :     case t_COL: /* nfelt */
    1665                 :            :       return RgM_RgC_mul(u, x);
    1666                 :            : 
    1667                 :            :     case t_MAT: /* ideal */
    1668                 :            :       y = cgetg_copy(x, &lx);
    1669                 :            :       for (i=1; i<lx; i++) gel(y,i) = RgM_RgC_mul(u, gel(x,i));
    1670                 :            :       break;
    1671                 :            : 
    1672                 :            :     case t_VEC: /* pr */
    1673                 :            :       checkprid(x); y = leafcopy(x);
    1674                 :            :       gel(y,2) = RgM_RgC_mul(u, gel(y,2));
    1675                 :            :       gel(y,5) = RgM_RgC_mul(u, gel(y,5));
    1676                 :            :       break;
    1677                 :            :     default: y = x;
    1678                 :            :   }
    1679                 :            :   return y;
    1680                 :            : }
    1681                 :            : 
    1682                 :            : GEN
    1683                 :            : nffromhnfbasis(GEN nf, GEN x)
    1684                 :            : {
    1685                 :            :   long tx = typ(x);
    1686                 :            :   pari_sp av = avma;
    1687                 :            :   GEN u;
    1688                 :            :   if (!is_vec_t(tx)) return gcopy(x);
    1689                 :            :   nf = checknf(nf);
    1690                 :            :   u = hnffromLLL(nf);
    1691                 :            :   return gerepilecopy(av, nfbasechange(u, x));
    1692                 :            : }
    1693                 :            : 
    1694                 :            : GEN
    1695                 :            : nftohnfbasis(GEN nf, GEN x)
    1696                 :            : {
    1697                 :            :   long tx = typ(x);
    1698                 :            :   pari_sp av = avma;
    1699                 :            :   GEN u;
    1700                 :            :   if (!is_vec_t(tx)) return gcopy(x);
    1701                 :            :   nf = checknf(nf);
    1702                 :            :   u = ZM_inv(hnffromLLL(nf), gen_1);
    1703                 :            :   return gerepilecopy(av, nfbasechange(u, x));
    1704                 :            : }
    1705                 :            : #endif
    1706                 :            : 
    1707                 :            : /* set *pro to roots of T->x */
    1708                 :            : static GEN
    1709                 :       3794 : get_red_G(nfbasic_t *T, GEN *pro)
    1710                 :            : {
    1711                 :       3794 :   GEN G, u, u0 = NULL;
    1712                 :            :   pari_sp av;
    1713                 :       3794 :   long i, prec, n = degpol(T->x);
    1714                 :            :   nffp_t F;
    1715                 :            : 
    1716                 :       3794 :   prec = nbits2prec(n+32);
    1717                 :       3794 :   nffp_init(&F, T, NULL, prec);
    1718                 :       3794 :   av = avma;
    1719                 :       3794 :   for (i=1; ; i++)
    1720                 :            :   {
    1721                 :       3794 :     F.prec = prec; make_M_G(&F, 0); G = F.G;
    1722         [ -  + ]:       3794 :     if (u0) G = RgM_mul(G, u0);
    1723         [ -  + ]:       3794 :     if (DEBUGLEVEL)
    1724                 :          0 :       err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
    1725                 :          0 :                   prec + F.extraprec, prec, F.extraprec);
    1726         [ +  - ]:       3794 :     if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST|LLL_COMPATIBLE)))
    1727                 :            :     {
    1728         [ +  - ]:       3794 :       if (lg(u)-1 == n) break;
    1729                 :            :       /* singular ==> loss of accuracy */
    1730         [ #  # ]:          0 :       if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
    1731                 :          0 :       else    u0 = gerepilecopy(av, u);
    1732                 :            :     }
    1733                 :          0 :     prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
    1734                 :          0 :     F.ro = NULL;
    1735         [ #  # ]:          0 :     if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
    1736                 :          0 :   }
    1737         [ -  + ]:       3794 :   if (u0) u = RgM_mul(u0,u);
    1738                 :       3794 :   *pro = F.ro; return u;
    1739                 :            : }
    1740                 :            : 
    1741                 :            : /* Compute an LLL-reduced basis for the integer basis of nf(T).
    1742                 :            :  * set *pro = roots of x if computed [NULL if not computed] */
    1743                 :            : static void
    1744                 :       6517 : set_LLL_basis(nfbasic_t *T, GEN *pro, double DELTA)
    1745                 :            : {
    1746                 :       6517 :   GEN B = T->bas;
    1747         [ +  + ]:       6517 :   if (T->r1 == degpol(T->x)) {
    1748                 :       2723 :     pari_sp av = avma;
    1749                 :       2723 :     GEN u, basden = T->basden;
    1750         [ +  - ]:       2723 :     if (!basden) basden = get_bas_den(B);
    1751                 :       2723 :     u = ZM_lll(make_Tr(T->x,basden), DELTA,
    1752                 :            :                LLL_GRAM|LLL_KEEP_FIRST|LLL_IM|LLL_COMPATIBLE);
    1753                 :       2723 :     B = gerepileupto(av, RgV_RgM_mul(B, u));
    1754                 :       2723 :     *pro = NULL;
    1755                 :            :   }
    1756                 :            :   else
    1757                 :       3794 :     B = RgV_RgM_mul(B, get_red_G(T, pro));
    1758                 :       6517 :   T->bas = B;
    1759                 :       6517 :   T->basden = get_bas_den(B);
    1760                 :       6517 : }
    1761                 :            : 
    1762                 :            : static int
    1763                 :       3151 : cmp_abs_ZX(GEN x, GEN y) { return gen_cmp_RgX((void*)&absi_cmp, x, y); }
    1764                 :            : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
    1765                 :            :  * (if so update *dx) */
    1766                 :            : static int
    1767                 :       4600 : ZX_is_better(GEN y, GEN x, GEN *dx)
    1768                 :            : {
    1769                 :       4600 :   GEN d = ZX_disc(y);
    1770                 :            :   int cmp;
    1771         [ +  + ]:       4600 :   if (!*dx) *dx = ZX_disc(x);
    1772                 :       4600 :   cmp = absi_cmp(d, *dx);
    1773         [ +  + ]:       4600 :   if (cmp < 0) { *dx = d; return 1; }
    1774         [ +  + ]:       3760 :   if (cmp == 0) return cmp_abs_ZX(y, x) < 0;
    1775                 :       4600 :   return 0;
    1776                 :            : }
    1777                 :            : 
    1778                 :            : static void polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
    1779                 :            : /* Seek a simpler, polynomial pol defining the same number field as
    1780                 :            :  * x (assumed to be monic at this point) */
    1781                 :            : static GEN
    1782                 :         70 : nfpolred(nfbasic_t *T, GEN *pro)
    1783                 :            : {
    1784                 :         70 :   GEN x = T->x, dx, b, rev;
    1785                 :         70 :   long n = degpol(x), v = varn(x);
    1786                 :            : 
    1787         [ -  + ]:         70 :   if (n == 1) {
    1788                 :          0 :     T->x = deg1pol_shallow(gen_1, gen_m1, v);
    1789                 :          0 :     *pro = NULL; return pol_1(v);
    1790                 :            :   }
    1791                 :         70 :   polredbest_aux(T, pro, &x, &dx, &b);
    1792         [ +  + ]:         70 :   if (x == T->x) return NULL; /* no improvement */
    1793         [ -  + ]:         56 :   if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
    1794                 :            : 
    1795                 :            :   /* update T */
    1796                 :         56 :   rev = QXQ_reverse(b, T->x);
    1797                 :         56 :   T->bas = QXV_QXQ_eval(T->bas, rev, x);
    1798                 :         56 :   (void)Z_issquareall(diviiexact(dx,T->dK), &(T->index));
    1799                 :         56 :   T->basden = get_bas_den(T->bas);
    1800                 :         56 :   T->dx = dx;
    1801                 :         56 :   T->x = x;
    1802                 :         56 :   *pro = NULL; /* reset */
    1803                 :         70 :   return rev;
    1804                 :            : }
    1805                 :            : 
    1806                 :            : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
    1807                 :            :  * basis. Assume bas[1] is 1 and that the leading coefficient of elements
    1808                 :            :  * of bas are of the form 1/b for a t_INT b */
    1809                 :            : GEN
    1810                 :        637 : get_nfindex(GEN bas)
    1811                 :            : {
    1812                 :        637 :   pari_sp av = avma;
    1813                 :        637 :   long n = lg(bas)-1, i;
    1814                 :            :   GEN D, d, mat;
    1815                 :            : 
    1816                 :        637 :   D = gen_1; /* assume bas[1] = 1 */
    1817         [ +  + ]:       1309 :   for (i = 2; i <= n; i++)
    1818                 :            :   { /* in most cases [e.g after nfbasis] basis is upper triangular! */
    1819                 :       1036 :     GEN B = gel(bas,i), lc;
    1820         [ +  + ]:       1036 :     if (degpol(B) != i-1) break;
    1821                 :        672 :     lc = gel(B, i+1);
    1822      [ +  +  - ]:        672 :     switch (typ(lc))
    1823                 :            :     {
    1824                 :        399 :       case t_INT: continue;
    1825                 :        273 :       case t_FRAC: lc = gel(lc,2); break;
    1826                 :          0 :       default: pari_err_TYPE("get_nfindex",lc);
    1827                 :            :     }
    1828                 :        273 :     D = mulii(D, lc);
    1829                 :            :   }
    1830         [ +  + ]:        637 :   if (i <= n)
    1831                 :            :   { /* not triangular after all */
    1832                 :        364 :     bas = Q_remove_denom(bas, &d);
    1833         [ +  + ]:        364 :     if (!d) { avma = av; return D; }
    1834                 :        350 :     mat = RgV_to_RgM(bas, n);
    1835                 :        350 :     d = diviiexact(powiu(d, n), ZM_det(mat));
    1836                 :        350 :     D = mulii(D,absi(d));
    1837                 :            :   }
    1838                 :        637 :   return gerepileuptoint(av, D);
    1839                 :            : }
    1840                 :            : 
    1841                 :            : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
    1842                 :            :  * basis (deprecated), or listP data (nfbasis input format) to specify
    1843                 :            :  * a set of primes at with the basis order must be maximal.
    1844                 :            :  * 1) nf type (or unrecognized): return t_VEC
    1845                 :            :  * 2) ZX or [ZX, listP]: return t_POL
    1846                 :            :  * 3) [ZX, order basis]: return 0 (deprecated)
    1847                 :            :  * incorrect: return -1 */
    1848                 :            : static long
    1849                 :       6587 : nf_input_type(GEN x)
    1850                 :            : {
    1851                 :            :   GEN T, V;
    1852                 :            :   long i, d, v;
    1853      [ +  +  - ]:       6587 :   switch(typ(x))
    1854                 :            :   {
    1855                 :       5845 :     case t_POL: return t_POL;
    1856                 :            :     case t_VEC:
    1857         [ +  + ]:        742 :       if (lg(x) != 3) return t_VEC; /* nf or incorrect */
    1858                 :        735 :       T = gel(x,1); V = gel(x,2);
    1859         [ -  + ]:        735 :       if (typ(T) != t_POL) return -1;
    1860      [ +  +  - ]:        735 :       switch(typ(V))
    1861                 :            :       {
    1862                 :         35 :         case t_INT: case t_MAT: return t_POL;
    1863                 :            :         case t_VEC: case t_COL:
    1864         [ +  + ]:        700 :           if (RgV_is_ZV(V)) return t_POL;
    1865                 :        686 :           break;
    1866                 :          0 :         default: return -1;
    1867                 :            :       }
    1868                 :        686 :       d = degpol(T); v = varn(T);
    1869 [ +  - ][ +  - ]:        686 :       if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
         [ +  - ][ +  + ]
    1870         [ +  + ]:       3976 :       for (i = 1; i <= d; i++)
    1871                 :            :       { /* check integer basis */
    1872                 :       3311 :         GEN c = gel(V,i);
    1873      [ +  +  - ]:       3311 :         switch(typ(c))
    1874                 :            :         {
    1875                 :         28 :           case t_INT: break;
    1876 [ +  + ][ +  - ]:       3283 :           case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
                 [ +  + ]
    1877                 :            :           /* fall through */
    1878                 :         14 :           default: return -1;
    1879                 :            :         }
    1880                 :            :       }
    1881                 :        665 :       return 0;
    1882                 :            :   }
    1883                 :       6587 :   return t_VEC; /* nf or incorrect */
    1884                 :            : }
    1885                 :            : 
    1886                 :            : static void
    1887                 :       5488 : nfbasic_add_disc(nfbasic_t *T)
    1888                 :            : {
    1889         [ +  + ]:       5488 :   if (!T->index) T->index = get_nfindex(T->bas);
    1890         [ +  + ]:       5488 :   if (!T->dx) T->dx = ZX_disc(T->x);
    1891         [ +  + ]:       5488 :   if (!T->dK) T->dK = diviiexact(T->dx, sqri(T->index));
    1892                 :       5488 : }
    1893                 :            : 
    1894                 :            : static void
    1895                 :       6587 : nfbasic_init(GEN x, long flag, nfbasic_t *T)
    1896                 :            : {
    1897                 :       6587 :   GEN bas, dK, dx, index, unscale = gen_1;
    1898                 :            :   long r1;
    1899                 :            : 
    1900                 :       6587 :   T->dKP = NULL;
    1901   [ +  +  +  + ]:       6587 :   switch (nf_input_type(x))
    1902                 :            :   {
    1903                 :            :     case t_POL:
    1904                 :            :     {
    1905                 :            :       nfmaxord_t S;
    1906                 :       5894 :       nfmaxord(&S, x, flag);
    1907                 :       5887 :       x = S.T;
    1908                 :       5887 :       T->x0 = S.T0;
    1909                 :       5887 :       T->dKP = S.dKP;
    1910                 :       5887 :       dK = S.dK;
    1911                 :       5887 :       index = S.index;
    1912                 :       5887 :       bas = S.basis;
    1913                 :       5887 :       dx = S.dT;
    1914                 :       5887 :       unscale = S.unscale;
    1915                 :       5887 :       r1 = ZX_sturm(x);
    1916                 :            :       break;
    1917                 :            :     }
    1918                 :            :     case t_VEC:
    1919                 :            :     { /* nf, bnf, bnr */
    1920                 :          7 :       GEN nf = checknf(x);
    1921                 :          7 :       x     = nf_get_pol(nf);
    1922                 :          7 :       dK    = nf_get_disc(nf);
    1923                 :          7 :       index = nf_get_index(nf);
    1924                 :          7 :       bas   = nf_get_zk(nf);
    1925                 :          7 :       T->x0 = x;
    1926                 :          7 :       dx = NULL;
    1927                 :          7 :       r1 = nf_get_r1(nf);
    1928                 :          7 :       break;
    1929                 :            :     }
    1930                 :            :     case 0: /* monic integral polynomial + integer basis */
    1931                 :        665 :       bas = gel(x,2); x = gel(x,1);
    1932                 :        665 :       T->x0 = x;
    1933                 :        665 :       index = NULL;
    1934                 :        665 :       dx = NULL;
    1935                 :        665 :       dK = NULL;
    1936                 :        665 :       r1 = ZX_sturm(x);
    1937                 :        665 :       break;
    1938                 :            :     default: /* -1 */
    1939                 :         21 :       pari_err_TYPE("nfbasic_init", x);
    1940                 :       6559 :       return;
    1941                 :            :   }
    1942                 :       6559 :   T->x     = x;
    1943                 :       6559 :   T->unscale = unscale;
    1944                 :       6559 :   T->r1    = r1;
    1945                 :       6559 :   T->dx    = dx;
    1946                 :       6559 :   T->dK    = dK;
    1947                 :       6559 :   T->bas   = bas;
    1948                 :       6559 :   T->basden= NULL;
    1949                 :       6559 :   T->index = index;
    1950                 :            : }
    1951                 :            : 
    1952                 :            : /* Initialize the number field defined by the polynomial x (in variable v)
    1953                 :            :  * flag & nf_RED:     try a polred first.
    1954                 :            :  * flag & nf_ORIG
    1955                 :            :  *    do a polred and return [nfinit(x), Mod(a,red)], where
    1956                 :            :  *    Mod(a,red) = Mod(v,x) (i.e return the base change). */
    1957                 :            : GEN
    1958                 :       5509 : nfinitall(GEN x, long flag, long prec)
    1959                 :            : {
    1960                 :       5509 :   const pari_sp av = avma;
    1961                 :            :   GEN nf, unscale;
    1962                 :            :   nfbasic_t T;
    1963                 :            : 
    1964         [ -  + ]:       5509 :   if (checkrnf_i(x)) return check_and_build_nfabs(x, prec);
    1965                 :       5509 :   nfbasic_init(x, flag, &T);
    1966         [ -  + ]:       5488 :   if (!ZX_is_irred(T.x)) pari_err_IRREDPOL("nfinit",x);
    1967 [ +  + ][ +  + ]:       5488 :   if (!equali1(leading_term(T.x0)) && !(flag & nf_RED))
    1968                 :            :   {
    1969                 :         49 :     pari_warn(warner,"non-monic polynomial. Result of the form [nf,c]");
    1970                 :         49 :     flag |= nf_RED | nf_ORIG;
    1971                 :            :   }
    1972                 :       5488 :   unscale = T.unscale;
    1973 [ +  + ][ +  + ]:       5488 :   if (!(flag & nf_RED) && !isint1(unscale))
    1974                 :            :   { /* implies lc(x0) = 1 and L := 1/unscale is integral */
    1975                 :         91 :     long d = degpol(T.x0);
    1976                 :         91 :     GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
    1977                 :         91 :     GEN f= powiu(L, (d*(d-1)) >> 1);
    1978                 :         91 :     T.x = T.x0; /* restore original user-supplied x0, unscale data */
    1979                 :         91 :     T.unscale = gen_1;
    1980                 :         91 :     T.dx    = gmul(T.dx, sqri(f));
    1981                 :         91 :     T.bas   = RgXV_unscale(T.bas, unscale);
    1982                 :         91 :     T.index = gmul(T.index, f);
    1983                 :            :   }
    1984                 :       5488 :   nfbasic_add_disc(&T); /* more expensive after set_LLL_basis */
    1985         [ +  + ]:       5488 :   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                 :         70 :     T.unscale = gen_1;
    1991                 :         70 :     rev = nfpolred(&T, &ro);
    1992                 :         70 :     nf = nfbasic_to_nf(&T, ro, prec);
    1993         [ +  + ]:         70 :     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                 :         70 :     T.unscale = unscale; /* restore */
    2000                 :            :   } else {
    2001                 :       5418 :     GEN ro; set_LLL_basis(&T, &ro, 0.99);
    2002                 :       5418 :     nf = nfbasic_to_nf(&T, ro, prec);
    2003                 :            :   }
    2004                 :       5488 :   return gerepilecopy(av, nf);
    2005                 :            : }
    2006                 :            : 
    2007                 :            : GEN
    2008                 :          0 : nfinitred(GEN x, long prec)  { return nfinitall(x, nf_RED, prec); }
    2009                 :            : GEN
    2010                 :          0 : nfinitred2(GEN x, long prec) { return nfinitall(x, nf_RED|nf_ORIG, prec); }
    2011                 :            : GEN
    2012                 :       3087 : nfinit(GEN x, long prec)     { return nfinitall(x, 0, prec); }
    2013                 :            : 
    2014                 :            : GEN
    2015                 :       1680 : nfinit0(GEN x, long flag,long prec)
    2016                 :            : {
    2017   [ +  +  +  - ]:       1680 :   switch(flag)
    2018                 :            :   {
    2019                 :            :     case 0:
    2020                 :       1659 :     case 1: return nfinitall(x,0,prec);
    2021                 :         14 :     case 2: case 4: return nfinitall(x,nf_RED,prec);
    2022                 :          7 :     case 3: case 5: return nfinitall(x,nf_RED|nf_ORIG,prec);
    2023                 :          0 :     default: pari_err_FLAG("nfinit");
    2024                 :            :   }
    2025                 :       1659 :   return NULL; /* not reached */
    2026                 :            : }
    2027                 :            : 
    2028                 :            : /* assume x a bnr/bnf/nf */
    2029                 :            : long
    2030                 :      66391 : nf_get_prec(GEN x)
    2031                 :            : {
    2032                 :      66391 :   GEN nf = checknf(x), ro = nf_get_roots(nf);
    2033         [ +  - ]:      66391 :   return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
    2034                 :            : }
    2035                 :            : 
    2036                 :            : /* assume nf is an nf */
    2037                 :            : GEN
    2038                 :        752 : nfnewprec_shallow(GEN nf, long prec)
    2039                 :            : {
    2040                 :        752 :   GEN NF = leafcopy(nf);
    2041                 :            :   nffp_t F;
    2042                 :        752 :   gel(NF,5) = leafcopy(gel(NF,5));
    2043                 :        752 :   remake_GM(NF, &F, prec);
    2044                 :        752 :   gel(NF,6) = F.ro;
    2045                 :        752 :   gmael(NF,5,1) = F.M;
    2046                 :        752 :   gmael(NF,5,2) = F.G;
    2047                 :        752 :   return NF;
    2048                 :            : }
    2049                 :            : 
    2050                 :            : GEN
    2051                 :         63 : nfnewprec(GEN nf, long prec)
    2052                 :            : {
    2053                 :            :   GEN z;
    2054   [ +  +  +  - ]:         63 :   switch(nftyp(nf))
    2055                 :            :   {
    2056                 :         49 :     default: pari_err_TYPE("nfnewprec", nf);
    2057                 :          7 :     case typ_BNF: z = bnfnewprec(nf,prec); break;
    2058                 :          7 :     case typ_BNR: z = bnrnewprec(nf,prec); break;
    2059                 :            :     case typ_NF: {
    2060                 :          0 :       pari_sp av = avma;
    2061                 :          0 :       z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
    2062                 :          0 :       break;
    2063                 :            :     }
    2064                 :            :   }
    2065                 :         14 :   return z;
    2066                 :            : }
    2067                 :            : 
    2068                 :            : /********************************************************************/
    2069                 :            : /**                                                                **/
    2070                 :            : /**                           POLRED                               **/
    2071                 :            : /**                                                                **/
    2072                 :            : /********************************************************************/
    2073                 :            : GEN
    2074                 :          0 : embednorm_T2(GEN x, long r1)
    2075                 :            : {
    2076                 :          0 :   pari_sp av = avma;
    2077                 :          0 :   GEN p = RgV_sumpart(x, r1);
    2078                 :          0 :   GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
    2079         [ #  # ]:          0 :   if (q != gen_0) p = gadd(p, gmul2n(q,1));
    2080         [ #  # ]:          0 :   return avma == av? gcopy(p): gerepileupto(av, p);
    2081                 :            : }
    2082                 :            : 
    2083                 :            : /* simplified version of gnorm for scalar, non-complex inputs, without GC */
    2084                 :            : static GEN
    2085                 :       6636 : real_norm(GEN x)
    2086                 :            : {
    2087   [ -  +  -  - ]:       6636 :   switch(typ(x))
    2088                 :            :   {
    2089                 :          0 :     case t_INT:  return sqri(x);
    2090                 :       6636 :     case t_REAL: return sqrr(x);
    2091                 :          0 :     case t_FRAC: return sqrfrac(x);
    2092                 :            :   }
    2093                 :          0 :   pari_err_TYPE("real_norm", x);
    2094                 :       6636 :   return NULL;
    2095                 :            : }
    2096                 :            : /* simplified version of gnorm, without GC */
    2097                 :            : static GEN
    2098                 :    1867783 : complex_norm(GEN x)
    2099                 :            : {
    2100         [ +  - ]:    1867783 :   return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
    2101                 :            : }
    2102                 :            : /* return T2(x), argument r1 needed in case x has components whose type
    2103                 :            :  * is unexpected, e.g. all of them t_INT for embed(gen_1) */
    2104                 :            : GEN
    2105                 :       1548 : embed_T2(GEN x, long r1)
    2106                 :            : {
    2107                 :       1548 :   pari_sp av = avma;
    2108                 :       1548 :   long i, l = lg(x);
    2109                 :       1548 :   GEN c, s = NULL, t = NULL;
    2110         [ -  + ]:       1548 :   if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
    2111         [ +  + ]:       8184 :   for (i = 1; i <= r1; i++)
    2112                 :            :   {
    2113                 :       6636 :     c = real_norm(gel(x,i));
    2114         [ +  + ]:       6636 :     s = s? gadd(s, c): c;
    2115                 :            :   }
    2116         [ +  + ]:       6233 :   for (; i < l; i++)
    2117                 :            :   {
    2118                 :       4685 :     c = complex_norm(gel(x,i));
    2119         [ +  + ]:       4685 :     t = t? gadd(t, c): c;
    2120                 :            :   }
    2121 [ +  + ][ +  + ]:       1548 :   if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
    2122                 :       1548 :   return gerepileupto(av, s);
    2123                 :            : }
    2124                 :            : /* return N(x) */
    2125                 :            : GEN
    2126                 :    1012493 : embed_norm(GEN x, long r1)
    2127                 :            : {
    2128                 :    1012493 :   pari_sp av = avma;
    2129                 :    1012493 :   long i, l = lg(x);
    2130                 :    1012493 :   GEN c, s = NULL, t = NULL;
    2131         [ +  + ]:    1012493 :   if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
    2132         [ +  + ]:    2691729 :   for (i = 1; i <= r1; i++)
    2133                 :            :   {
    2134                 :    1682370 :     c = gel(x,i);
    2135         [ +  + ]:    1682370 :     s = s? gmul(s, c): c;
    2136                 :            :   }
    2137         [ +  + ]:    2872457 :   for (; i < l; i++)
    2138                 :            :   {
    2139                 :    1863098 :     c = complex_norm(gel(x,i));
    2140         [ +  + ]:    1863098 :     t = t? gmul(t, c): c;
    2141                 :            :   }
    2142 [ +  + ][ +  + ]:    1009359 :   if (t) s = s? gmul(s,t): t;
    2143                 :    1012493 :   return gerepileupto(av, s);
    2144                 :            : }
    2145                 :            : 
    2146                 :            : typedef struct {
    2147                 :            :   long r1, v, prec;
    2148                 :            :   GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
    2149                 :            :   GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
    2150                 :            :   GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
    2151                 :            :   GEN bound; /* T2 norm of the polynomial defining nf */
    2152                 :            :   long expo_best_disc; /* expo(disc(x)), best generator so far */
    2153                 :            : } CG_data;
    2154                 :            : 
    2155                 :            : /* characteristic pol of x (given by embeddings) */
    2156                 :            : static GEN
    2157                 :      24715 : get_pol(CG_data *d, GEN x)
    2158                 :            : {
    2159                 :            :   long e;
    2160                 :      24715 :   GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
    2161         [ +  + ]:      24715 :   return (e > -5)? NULL: g;
    2162                 :            : }
    2163                 :            : 
    2164                 :            : /* characteristic pol of x (given as vector on (w_i)) */
    2165                 :            : static GEN
    2166                 :      10586 : get_polchar(CG_data *d, GEN x)
    2167                 :      10586 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
    2168                 :            : 
    2169                 :            : /* Choose a canonical polynomial in the pair { z(X), (+/-)z(-X) }.
    2170                 :            :  * z a ZX with lc(z) > 0. We want to keep that property, while
    2171                 :            :  * ensuring that the leading coeff of the odd (resp. even) part of z is < 0
    2172                 :            :  * if deg z is even (resp. odd).
    2173                 :            :  * Either leave z alone (return 1) or set z <-- (-1)^deg(z) z(-X). In place. */
    2174                 :            : static int
    2175                 :      10705 : ZX_canon_neg(GEN z)
    2176                 :            : {
    2177                 :            :   long i,s;
    2178                 :            : 
    2179         [ +  + ]:      82840 :   for (i = lg(z)-2; i >= 2; i -= 2)
    2180                 :            :   { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
    2181                 :      77534 :     s = signe(gel(z,i));
    2182         [ +  + ]:      77534 :     if (!s) continue;
    2183                 :            :     /* non trivial */
    2184         [ +  + ]:       5399 :     if (s < 0) break; /* the condition is already satisfied */
    2185                 :            : 
    2186         [ +  + ]:       8997 :     for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
    2187                 :       1807 :     return 1;
    2188                 :            :   }
    2189                 :      10705 :   return 0;
    2190                 :            : }
    2191                 :            : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
    2192                 :            :  * Return NULL on failure: discriminant too large or non primitive */
    2193                 :            : static GEN
    2194                 :      20251 : try_polmin(CG_data *d, nfbasic_t *T, GEN v, long flag, GEN *ai)
    2195                 :            : {
    2196                 :      20251 :   const long best = flag & nf_ABSOLUTE;
    2197                 :            :   long ed;
    2198                 :      20251 :   pari_sp av = avma;
    2199                 :            :   GEN g;
    2200         [ +  + ]:      20251 :   if (best)
    2201                 :            :   {
    2202                 :      19376 :     ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
    2203         [ +  + ]:      19376 :     avma = av; if (d->expo_best_disc < ed) return NULL;
    2204                 :            :   }
    2205                 :            :   else
    2206                 :        875 :     ed = 0;
    2207                 :      11662 :   g = get_pol(d, v);
    2208                 :            :   /* accuracy too low, compute algebraically */
    2209         [ +  + ]:      11662 :   if (!g) { avma = av; g = ZXQ_charpoly(*ai, T->x, varn(T->x)); }
    2210                 :      11662 :   (void)ZX_gcd_all(g, ZX_deriv(g), &g);
    2211 [ +  + ][ +  + ]:      11662 :   if (best && degpol(g) != degpol(T->x)) { avma = av; return NULL; }
    2212                 :       4697 :   g = gerepilecopy(av, g);
    2213                 :       4697 :   d->expo_best_disc = ed;
    2214         [ +  + ]:       4697 :   if (flag & nf_ORIG)
    2215                 :            :   {
    2216         [ +  + ]:        980 :     if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
    2217         [ +  + ]:        980 :     if (!isint1(T->unscale)) *ai = RgX_unscale(*ai, T->unscale);
    2218                 :            :   }
    2219                 :            :   else
    2220                 :       3717 :     (void)ZX_canon_neg(g);
    2221         [ -  + ]:       4697 :   if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
    2222                 :      20251 :   return g;
    2223                 :            : }
    2224                 :            : 
    2225                 :            : /* does x generate the correct field ? */
    2226                 :            : static GEN
    2227                 :      10586 : chk_gen(void *data, GEN x)
    2228                 :            : {
    2229                 :      10586 :   pari_sp av = avma, av1;
    2230                 :      10586 :   GEN h, g = get_polchar((CG_data*)data,x);
    2231         [ -  + ]:      10586 :   if (!g) pari_err_PREC("chk_gen");
    2232                 :      10586 :   av1 = avma;
    2233                 :      10586 :   h = ZX_gcd(g, ZX_deriv(g));
    2234         [ +  + ]:      10586 :   if (degpol(h)) { avma = av; return NULL; }
    2235         [ -  + ]:       6043 :   if (DEBUGLEVEL>3) err_printf("  generator: %Ps\n",g);
    2236                 :      10586 :   avma = av1; return gerepileupto(av, g);
    2237                 :            : }
    2238                 :            : 
    2239                 :            : static long
    2240                 :       1435 : chk_gen_prec(long N, long bit)
    2241                 :       1435 : { return nbits2prec(10 + (long)log2((double)N) + bit); }
    2242                 :            : 
    2243                 :            : /* Remove duplicate polynomials in P, updating A (same indices), in place.
    2244                 :            :  * Among elements having the same characteristic pol, choose the smallest
    2245                 :            :  * according to ZV_abscmp */
    2246                 :            : static void
    2247                 :        336 : remove_duplicates(GEN P, GEN A)
    2248                 :            : {
    2249                 :        336 :   long k, i, l = lg(P);
    2250                 :        336 :   pari_sp av = avma;
    2251                 :            :   GEN x, a;
    2252                 :            : 
    2253         [ -  + ]:        672 :   if (l < 2) return;
    2254                 :        336 :   (void)sort_factor_pol(mkmat2(P, A), cmpii);
    2255                 :        336 :   x = gel(P,1); a = gel(A,1);
    2256         [ +  + ]:       6008 :   for  (k=1,i=2; i<l; i++)
    2257         [ +  + ]:       5672 :     if (ZX_equal(gel(P,i), x))
    2258                 :            :     {
    2259         [ +  + ]:       3480 :       if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
    2260                 :            :     }
    2261                 :            :     else
    2262                 :            :     {
    2263                 :       2192 :       gel(A,k) = a;
    2264                 :       2192 :       gel(P,k) = x;
    2265                 :       2192 :       k++;
    2266                 :       2192 :       x = gel(P,i); a = gel(A,i);
    2267                 :            :     }
    2268                 :        336 :   gel(A,k) = a;
    2269                 :        336 :   gel(P,k) = x;
    2270                 :        336 :   l = k+1; setlg(A,l); setlg(P,l);
    2271                 :        336 :   avma = av;
    2272                 :            : }
    2273                 :            : 
    2274                 :            : static long
    2275                 :       1099 : polred_init(nfbasic_t *T, nffp_t *F, CG_data *d)
    2276                 :            : {
    2277                 :       1099 :   long e, prec, n = degpol(T->x);
    2278                 :            :   double log2rho;
    2279                 :            :   GEN ro;
    2280                 :       1099 :   set_LLL_basis(T, &ro, 0.9999);
    2281                 :            :   /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
    2282         [ +  + ]:       1099 :   log2rho = ro ? (double)gexpo(ro): fujiwara_bound(T->x);
    2283                 :       1099 :   e = n * (long)(log2rho + log2((double)n)) + 1;
    2284         [ -  + ]:       1099 :   if (e < 0) e = 0; /* can occur if n = 1 */
    2285                 :       1099 :   prec = chk_gen_prec(n, e);
    2286                 :       1099 :   get_nf_fp_compo(T, F, ro, 1, prec);
    2287                 :       1099 :   d->v = varn(T->x);
    2288                 :       1099 :   d->expo_best_disc = -1;
    2289                 :       1099 :   d->ZKembed = NULL;
    2290                 :       1099 :   d->M = NULL;
    2291                 :       1099 :   d->u = NULL;
    2292                 :       1099 :   d->r1= T->r1; return prec;
    2293                 :            : }
    2294                 :            : static GEN
    2295                 :        343 : findmindisc(GEN y, GEN *pa)
    2296                 :            : {
    2297                 :        343 :   GEN a = *pa, x = gel(y,1), b = gel(a,1), dx = NULL;
    2298                 :        343 :   long i, l = lg(y);
    2299         [ +  + ]:        414 :   for (i = 2; i < l; i++)
    2300                 :            :   {
    2301                 :         71 :     GEN yi = gel(y,i);
    2302         [ +  + ]:         71 :     if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
    2303                 :            :   }
    2304                 :        343 :   *pa = b; return x;
    2305                 :            : }
    2306                 :            : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
    2307                 :            :  * [ the best wrt discriminant ordering ], but keep all non-primitive
    2308                 :            :  * polynomials */
    2309                 :            : static void
    2310                 :        763 : filter(GEN y, GEN b, long n)
    2311                 :            : {
    2312                 :            :   GEN x, a, dx;
    2313                 :        763 :   long i, k = 1, l = lg(y);
    2314                 :        763 :   a = x = dx = NULL;
    2315         [ +  + ]:       5523 :   for (i = 1; i < l; i++)
    2316                 :            :   {
    2317                 :       4760 :     GEN yi = gel(y,i), ai = gel(b,i);
    2318         [ +  + ]:       4760 :     if (degpol(yi) == n)
    2319                 :            :     {
    2320                 :       4592 :       pari_sp av = avma;
    2321 [ +  + ][ +  + ]:       4592 :       if (dx && !ZX_is_better(yi,x,&dx)) { avma = av; continue; }
    2322         [ +  + ]:       1078 :       if (!dx) dx = ZX_disc(yi);
    2323                 :       1078 :       x = yi; a = ai; continue;
    2324                 :            :     }
    2325                 :        168 :     gel(y,k) = yi;
    2326                 :        168 :     gel(b,k) = ai; k++;
    2327                 :            :   }
    2328         [ +  - ]:        763 :   if (dx)
    2329                 :            :   {
    2330                 :        763 :     gel(y,k) = x;
    2331                 :        763 :     gel(b,k) = a; k++;
    2332                 :            :   }
    2333                 :        763 :   setlg(y, k);
    2334                 :        763 :   setlg(b, k);
    2335                 :        763 : }
    2336                 :            : 
    2337                 :            : static GEN
    2338                 :        791 : polred_aux(nfbasic_t *T, GEN *pro, long flag)
    2339                 :            : { /* only keep polynomials of max degree and best discriminant */
    2340                 :        791 :   const long best = flag & nf_ABSOLUTE;
    2341                 :        791 :   const long orig = flag & nf_ORIG;
    2342                 :        791 :   GEN M, b, y, x = T->x;
    2343                 :        791 :   long maxi, i, j, k, v = varn(x), n = lg(T->bas)-1;
    2344                 :            :   nffp_t F;
    2345                 :            :   CG_data d;
    2346                 :            : 
    2347         [ +  + ]:        791 :   if (n == 1)
    2348                 :            :   {
    2349         [ +  + ]:         28 :     if (!best)
    2350                 :            :     {
    2351                 :         14 :       GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2352         [ -  + ]:         14 :       return orig? mkmat2(mkcol(ch),mkcol(gen_1)): mkvec(ch);
    2353                 :            :     }
    2354                 :            :     else
    2355         [ +  + ]:         14 :       return orig? trivial_fact(): cgetg(1,t_VEC);
    2356                 :            :   }
    2357                 :            : 
    2358                 :        763 :   (void)polred_init(T, &F, &d);
    2359                 :        763 :   *pro = F.ro;
    2360                 :        763 :   M = F.M;
    2361         [ +  + ]:        763 :   if (best)
    2362                 :            :   {
    2363         [ -  + ]:        700 :     if (!T->dx) T->dx = ZX_disc(T->x);
    2364                 :        700 :     d.expo_best_disc = expi(T->dx);
    2365                 :            :   }
    2366                 :            : 
    2367                 :            :   /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
    2368                 :        763 :   y = cgetg(n*n + 1, t_VEC);
    2369                 :        763 :   b = cgetg(n*n + 1, t_COL);
    2370                 :        763 :   k = 1;
    2371         [ +  + ]:        763 :   if (!best)
    2372                 :            :   {
    2373                 :         63 :     GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2374                 :         63 :     gel(y,1) = ch; gel(b,1) = gen_1; k++;
    2375                 :            :   }
    2376         [ +  + ]:       4256 :   for (i = 2; i <= n; i++)
    2377                 :            :   {
    2378                 :            :     GEN ch, ai;
    2379                 :       3493 :     ai = gel(T->bas,i);
    2380                 :       3493 :     ch = try_polmin(&d, T, gel(M,i), flag, &ai);
    2381         [ +  + ]:       3493 :     if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2382                 :            :   }
    2383                 :        763 :   maxi = minss(n, 3);
    2384         [ +  + ]:       2863 :   for (i = 1; i <= maxi; i++)
    2385         [ +  + ]:      10479 :     for (j = i+1; j <= n; j++)
    2386                 :            :     {
    2387                 :            :       GEN ch, ai, v;
    2388                 :       8379 :       ai = gadd(gel(T->bas,i), gel(T->bas,j));
    2389                 :       8379 :       v = RgV_add(gel(M,i), gel(M,j));
    2390                 :            :       /* defining polynomial for Q(w_i+w_j) */
    2391                 :       8379 :       ch = try_polmin(&d, T, v, flag, &ai);
    2392         [ +  + ]:       8379 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2393                 :            : 
    2394                 :       8379 :       ai = gsub(gel(T->bas,i), gel(T->bas,j));
    2395                 :       8379 :       v = RgV_sub(gel(M,i), gel(M,j));
    2396                 :            :       /* defining polynomial for Q(w_i-w_j) */
    2397                 :       8379 :       ch = try_polmin(&d, T, v, flag, &ai);
    2398         [ +  + ]:       8379 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2399                 :            :     }
    2400                 :        763 :   setlg(y, k);
    2401                 :        763 :   setlg(b, k); filter(y, b, n);
    2402         [ +  + ]:        763 :   if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
    2403                 :        140 :   (void)sort_factor_pol(mkmat2(y, b), cmpii);
    2404                 :        791 :   settyp(y, t_COL); return mkmat2(b, y);
    2405                 :            : }
    2406                 :            : 
    2407                 :            : static GEN
    2408                 :         84 : Polred(GEN x, long flag, GEN fa)
    2409                 :            : {
    2410                 :         84 :   pari_sp av = avma;
    2411                 :            :   GEN ro;
    2412         [ +  + ]:         84 :   nfbasic_t T; nfbasic_init(fa? mkvec2(x,fa): x, flag & nf_PARTIALFACT, &T);
    2413                 :         77 :   return gerepilecopy(av, polred_aux(&T, &ro, flag));
    2414                 :            : }
    2415                 :            : 
    2416                 :            : /* finds "best" polynomial in polred_aux list, defaulting to T->x if none of
    2417                 :            :  * them is primitive. *px is the ZX, characteristic polynomial of Mod(*pb,T->x),
    2418                 :            :  * *pdx its discriminant. Set *pro = polroots(T->x) [ NOT *px ]. */
    2419                 :            : static void
    2420                 :        714 : polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
    2421                 :            : {
    2422                 :        714 :   GEN y, x = T->x; /* default value */
    2423                 :            :   long i, l;
    2424         [ +  + ]:        714 :   y = polred_aux(T, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
    2425                 :        714 :   *pdx = T->dx;
    2426         [ +  + ]:        714 :   if (pb)
    2427                 :            :   {
    2428                 :        126 :     GEN a, b = deg1pol_shallow(T->unscale, gen_0, varn(x));
    2429                 :        126 :     a = gel(y,1); l = lg(a);
    2430                 :        126 :     y = gel(y,2);
    2431         [ +  + ]:        245 :     for (i=1; i<l; i++)
    2432                 :            :     {
    2433                 :        119 :       GEN yi = gel(y,i);
    2434                 :        119 :       pari_sp av = avma;
    2435         [ +  + ]:        119 :       if (ZX_is_better(yi,x,pdx)) { x = yi; b = gel(a,i); } else avma = av;
    2436                 :            :     }
    2437                 :        126 :     *pb = b;
    2438                 :            :   }
    2439                 :            :   else
    2440                 :            :   {
    2441                 :        588 :     l = lg(y);
    2442         [ +  + ]:       1169 :     for (i=1; i<l; i++)
    2443                 :            :     {
    2444                 :        581 :       GEN yi = gel(y,i);
    2445                 :        581 :       pari_sp av = avma;
    2446         [ +  + ]:        581 :       if (ZX_is_better(yi,x,pdx)) x = yi; else avma = av;
    2447                 :            :     }
    2448                 :            :   }
    2449         [ -  + ]:        714 :   if (!*pdx) *pdx = ZX_disc(x);
    2450                 :        714 :   *px = x;
    2451                 :        714 : }
    2452                 :            : GEN
    2453                 :        616 : polredbest(GEN T0, long flag)
    2454                 :            : {
    2455                 :        616 :   pari_sp av = avma;
    2456                 :            :   GEN T, dT, ro, a;
    2457                 :            :   nfbasic_t S;
    2458 [ +  - ][ -  + ]:        616 :   if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
    2459                 :        616 :   T = T0; nfbasic_init(T, nf_PARTIALFACT, &S);
    2460         [ +  + ]:        616 :   polredbest_aux(&S, &ro, &T, &dT, flag? &a: NULL);
    2461         [ +  + ]:        616 :   if (flag)
    2462                 :            :   { /* charpoly(Mod(a,T0)) = T */
    2463                 :            :     GEN b;
    2464         [ -  + ]:         28 :     if (T0 == T)
    2465                 :          0 :       b = pol_x(varn(T)); /* no improvement */
    2466                 :            :     else
    2467                 :         28 :       b = QXQ_reverse(a, T0); /* charpoly(Mod(b,T)) = S.x */
    2468         [ +  + ]:         28 :     b = (degpol(T) == 1)? gmodulo(b, T): mkpolmod(b,T);
    2469                 :         28 :     T = mkvec2(T, b);
    2470                 :            :   }
    2471                 :        616 :   return gerepilecopy(av, T);
    2472                 :            : }
    2473                 :            : /* DEPRECATED: backward compatibility */
    2474                 :            : GEN
    2475                 :         70 : polred0(GEN x, long flag, GEN fa)
    2476                 :            : {
    2477                 :         70 :   long fl = 0;
    2478         [ +  + ]:         70 :   if (flag & 1) fl |= nf_PARTIALFACT;
    2479         [ +  + ]:         70 :   if (flag & 2) fl |= nf_ORIG;
    2480                 :         70 :   return Polred(x, fl, fa);
    2481                 :            : }
    2482                 :            : 
    2483                 :            : GEN
    2484                 :         21 : polredord(GEN x)
    2485                 :            : {
    2486                 :         21 :   pari_sp av = avma;
    2487                 :            :   GEN v, lt;
    2488                 :            :   long i, n, vx;
    2489                 :            : 
    2490         [ -  + ]:         21 :   if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
    2491                 :         21 :   x = Q_primpart(x); RgX_check_ZX(x,"polredord");
    2492         [ -  + ]:         21 :   n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
    2493         [ +  + ]:         21 :   if (n == 1) return gerepilecopy(av, mkvec(x));
    2494                 :         14 :   lt = leading_term(x); vx = varn(x);
    2495         [ +  + ]:         14 :   if (is_pm1(lt))
    2496                 :            :   {
    2497         [ -  + ]:          7 :     if (signe(lt) < 0) x = ZX_neg(x);
    2498                 :          7 :     v = pol_x_powers(n, vx);
    2499                 :            :   }
    2500                 :            :   else
    2501                 :            :   { GEN L;
    2502                 :            :     /* basis for Dedekind order */
    2503                 :          7 :     v = cgetg(n+1, t_VEC);
    2504                 :          7 :     gel(v,1) = scalarpol_shallow(lt, vx);
    2505         [ +  + ]:         14 :     for (i = 2; i <= n; i++)
    2506                 :          7 :       gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
    2507                 :          7 :     gel(v,1) = pol_1(vx);
    2508                 :          7 :     x = ZX_Q_normalize(x, &L);
    2509                 :          7 :     v = gsubst(v, vx, monomial(ginv(L),1,vx));
    2510         [ +  + ]:         14 :     for (i=2; i <= n; i++)
    2511         [ +  - ]:          7 :       if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = monomial(gen_1, i-1, vx);
    2512                 :            :   }
    2513                 :         21 :   return gerepileupto(av, polred(mkvec2(x, v)));
    2514                 :            : }
    2515                 :            : 
    2516                 :            : GEN
    2517                 :         14 : polred(GEN x) { return Polred(x, 0, NULL); }
    2518                 :            : GEN
    2519                 :          0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
    2520                 :            : GEN
    2521                 :          0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
    2522                 :            : GEN
    2523                 :          0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
    2524                 :            : GEN
    2525                 :          0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
    2526                 :            : GEN
    2527                 :          0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
    2528                 :            : 
    2529                 :            : /********************************************************************/
    2530                 :            : /**                                                                **/
    2531                 :            : /**                           POLREDABS                            **/
    2532                 :            : /**                                                                **/
    2533                 :            : /********************************************************************/
    2534                 :            : /* set V[k] := matrix of multiplication by nk.zk[k] */
    2535                 :            : static GEN
    2536                 :       1295 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
    2537                 :            : {
    2538                 :       1295 :   GEN v, Mk = cgetg(N+1, t_MAT);
    2539                 :            :   long i, e;
    2540         [ +  + ]:       4452 :   for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
    2541         [ +  + ]:      12278 :   for (     ; i <=N; i++)
    2542                 :            :   {
    2543                 :      10983 :     v = vecmul(gel(M,k), gel(M,i));
    2544                 :      10983 :     v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
    2545                 :      10983 :     gel(Mk,i) = grndtoi(v, &e);
    2546         [ -  + ]:      10983 :     if (e > -5) return NULL;
    2547                 :            :   }
    2548                 :       1295 :   gel(V,k) = Mk; return Mk;
    2549                 :            : }
    2550                 :            : 
    2551                 :            : static GEN
    2552                 :       1427 : ZM_image_shallow(GEN M, long *pr)
    2553                 :            : {
    2554                 :            :   long j, k, r;
    2555                 :       1427 :   GEN y, d = ZM_pivots(M, &k);
    2556                 :       1427 :   r = lg(M)-1 - k;
    2557                 :       1427 :   y = cgetg(r+1,t_MAT);
    2558         [ +  + ]:      10550 :   for (j=k=1; j<=r; k++)
    2559         [ +  + ]:       9123 :     if (d[k]) gel(y,j++) = gel(M,k);
    2560                 :       1427 :   *pr = r; return y;
    2561                 :            : }
    2562                 :            : 
    2563                 :            : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
    2564                 :            :  * Q from algo 2.7.6] */
    2565                 :            : static GEN
    2566                 :        344 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
    2567                 :            : {
    2568                 :        344 :   CG_data *d = (CG_data*)chk->data;
    2569                 :            :   GEN P, V, D, inv, bound, S, M;
    2570                 :        344 :   long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
    2571                 :        344 :   long i, j, prec, firstprim = 0, skipfirst = 0;
    2572                 :            :   pari_sp av;
    2573                 :            : 
    2574                 :        344 :   d->u = U;
    2575                 :        344 :   d->ZKembed = M = RgM_mul(d->M, U);
    2576                 :            : 
    2577                 :        344 :   av = avma; bound = d->bound;
    2578                 :        344 :   D = cgetg(N+1, t_VECSMALL);
    2579         [ +  + ]:       2789 :   for (i = 1; i <= N; i++)
    2580                 :            :   {
    2581                 :       2453 :     pari_sp av2 = avma;
    2582                 :       2453 :     P = get_pol(d, gel(M,i));
    2583         [ +  + ]:       2453 :     if (!P) pari_err_PREC("chk_gen_init");
    2584                 :       2445 :     (void)ZX_gcd_all(P, ZX_deriv(P), &P);
    2585                 :       2445 :     P = gerepilecopy(av2, P);
    2586                 :       2445 :     D[i] = degpol(P);
    2587         [ +  + ]:       2445 :     if (D[i] == N)
    2588                 :            :     { /* primitive element */
    2589                 :       1002 :       GEN B = embed_T2(gel(M,i), r1);
    2590         [ +  + ]:       1002 :       if (!firstprim) firstprim = i; /* index of first primitive element */
    2591         [ -  + ]:       1002 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2592         [ +  + ]:       1002 :       if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
    2593                 :            :     }
    2594                 :            :     else
    2595                 :            :     {
    2596         [ -  + ]:       1443 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
    2597         [ +  + ]:       1443 :       if (firstprim)
    2598                 :            :       { /* cycle basis vectors so that primitive elements come last */
    2599                 :        159 :         GEN u = d->u, e = M;
    2600                 :        159 :         GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
    2601                 :        159 :         long tS = D[i];
    2602         [ +  + ]:        447 :         for (j = i; j > firstprim; j--)
    2603                 :            :         {
    2604                 :        288 :           u[j] = u[j-1];
    2605                 :        288 :           e[j] = e[j-1];
    2606                 :        288 :           R[j] = R[j-1];
    2607                 :        288 :           D[j] = D[j-1];
    2608                 :            :         }
    2609                 :        159 :         gel(u,firstprim) = tu;
    2610                 :        159 :         gel(e,firstprim) = te;
    2611                 :        159 :         gel(R,firstprim) = tR;
    2612                 :        159 :         D[firstprim] = tS; firstprim++;
    2613                 :            :       }
    2614                 :            :     }
    2615                 :            :   }
    2616         [ +  + ]:        336 :   if (!firstprim)
    2617                 :            :   { /* try (a little) to find primitive elements to improve bound */
    2618                 :         21 :     GEN x = cgetg(N+1, t_VECSMALL);
    2619         [ -  + ]:         21 :     if (DEBUGLEVEL>1)
    2620                 :          0 :       err_printf("chk_gen_init: difficult field, trying random elements\n");
    2621         [ +  + ]:        231 :     for (i = 0; i < 10; i++)
    2622                 :            :     {
    2623                 :            :       GEN e, B;
    2624         [ +  + ]:       2310 :       for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
    2625                 :        210 :       e = RgM_zc_mul(M, x);
    2626                 :        210 :       B = embed_T2(e, r1);
    2627         [ +  + ]:        210 :       if (gcmp(B,bound) >= 0) continue;
    2628         [ -  + ]:         14 :       P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
    2629         [ -  + ]:         14 :       if (!ZX_is_squarefree(P)) continue;
    2630         [ -  + ]:         14 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2631                 :         14 :       bound = B ;
    2632                 :            :     }
    2633                 :            :   }
    2634                 :            : 
    2635         [ +  - ]:        336 :   if (firstprim != 1)
    2636                 :            :   {
    2637                 :        336 :     inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
    2638                 :        336 :     V = gel(inv,1);
    2639         [ +  + ]:       1533 :     for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
    2640                 :            :     /* V corresponds to 1_Z */
    2641                 :        336 :     V = grndtoi(V, &j);
    2642         [ -  + ]:        336 :     if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
    2643                 :        336 :     S = mkmat(V); /* 1 */
    2644                 :            : 
    2645                 :        336 :     V = cgetg(N+1, t_VEC);
    2646         [ +  - ]:       1442 :     for (i = 1; i <= N; i++,skipfirst++)
    2647                 :            :     { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
    2648                 :            :       GEN Mx, M2;
    2649                 :       1442 :       long j, k, h, rkM, dP = D[i];
    2650                 :            : 
    2651         [ +  + ]:       1442 :       if (dP == N) break; /* primitive */
    2652                 :       1295 :       Mx = set_mulid(V, M, inv, r1, r2, N, i);
    2653         [ -  + ]:       1295 :       if (!Mx) break; /* prec. problem. Stop */
    2654         [ +  + ]:       1295 :       if (dP == 1) continue;
    2655                 :       1008 :       rkM = lg(S)-1;
    2656                 :       1008 :       M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
    2657                 :       1008 :       gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
    2658                 :       1008 :       k = 2;
    2659         [ +  + ]:       2435 :       for (h = 1; h < dP; h++)
    2660                 :            :       {
    2661                 :            :         long r; /* add to M2 the elts of S * nf.zk[i]  */
    2662         [ +  + ]:       7801 :         for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
    2663                 :       1427 :         setlg(M2, k); k = 1;
    2664                 :       1427 :         S = ZM_image_shallow(shallowconcat(S,M2), &r);
    2665         [ +  + ]:       1427 :         if (r == rkM) break;
    2666         [ +  - ]:        895 :         if (r > rkM)
    2667                 :            :         {
    2668                 :        895 :           rkM = r;
    2669         [ +  + ]:        895 :           if (rkM == N) break;
    2670                 :            :         }
    2671                 :            :       }
    2672         [ +  + ]:       1008 :       if (rkM == N) break;
    2673                 :            :       /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
    2674                 :            :     }
    2675                 :            :   }
    2676                 :            :   /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
    2677                 :        336 :   chk->skipfirst = skipfirst;
    2678         [ -  + ]:        336 :   if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
    2679                 :            : 
    2680                 :            :   /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
    2681                 :        336 :   bound = gerepileuptoleaf(av, bound);
    2682                 :        336 :   prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
    2683         [ -  + ]:        336 :   if (DEBUGLEVEL)
    2684                 :          0 :     err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
    2685         [ -  + ]:        336 :   if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
    2686         [ +  + ]:        336 :   if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
    2687                 :        336 :   return bound;
    2688                 :            : }
    2689                 :            : 
    2690                 :            : /* z "small" minimal polynomial of Mod(a,x), deg z = deg x */
    2691                 :            : static GEN
    2692                 :       2471 : store(GEN x, GEN z, GEN a, nfbasic_t *T, long flag, GEN u)
    2693                 :            : {
    2694                 :            :   GEN y, b;
    2695                 :            : 
    2696         [ +  + ]:       2471 :   if (u) a = RgV_RgC_mul(T->bas, ZM_ZC_mul(u, a));
    2697         [ +  + ]:       2471 :   if (flag & (nf_ORIG|nf_ADDZK))
    2698                 :            :   {
    2699                 :        245 :     b = QXQ_reverse(a, x);
    2700         [ +  + ]:        245 :     if (!isint1(T->unscale)) b = gdiv(b, T->unscale); /* not RgX_Rg_div */
    2701                 :            :   }
    2702                 :            :   else
    2703                 :       2226 :     b = NULL;
    2704                 :            : 
    2705         [ +  + ]:       2471 :   if (flag & nf_RAW)
    2706                 :         28 :     y = mkvec2(z, a);
    2707         [ +  + ]:       2443 :   else if (flag & nf_ORIG) /* store phi(b mod z). */
    2708                 :        245 :     y = mkvec2(z, mkpolmod(b,z));
    2709                 :            :   else
    2710                 :       2198 :     y = z;
    2711         [ -  + ]:       2471 :   if (flag & nf_ADDZK)
    2712                 :            :   { /* append integral basis for number field Q[X]/(z) to result */
    2713                 :          0 :     long n = degpol(x);
    2714                 :          0 :     GEN t = RgV_RgM_mul(RgXQ_powers(b, n-1, z), RgV_to_RgM(T->bas,n));
    2715                 :          0 :     y = mkvec2(y, t);
    2716                 :            :   }
    2717                 :       2471 :   return y;
    2718                 :            : }
    2719                 :            : static GEN
    2720                 :        336 : polredabs_aux(nfbasic_t *T, GEN *u)
    2721                 :            : {
    2722                 :            :   long prec;
    2723                 :            :   GEN v;
    2724                 :        336 :   FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
    2725                 :            :   nffp_t F;
    2726                 :        336 :   CG_data d; chk.data = (void*)&d;
    2727                 :            : 
    2728                 :        336 :   prec = polred_init(T, &F, &d);
    2729                 :        336 :   d.bound = embed_T2(F.ro, d.r1);
    2730         [ +  + ]:        336 :   if (realprec(d.bound) > prec) d.bound = rtor(d.bound, prec);
    2731                 :            :   for (;;)
    2732                 :            :   {
    2733                 :        356 :     GEN R = R_from_QR(F.G, prec);
    2734         [ +  + ]:        356 :     if (R)
    2735                 :            :     {
    2736                 :        344 :       d.prec = prec;
    2737                 :        344 :       d.M    = F.M;
    2738                 :        344 :       v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
    2739         [ +  + ]:        344 :       if (v) break;
    2740                 :            :     }
    2741                 :         20 :     prec = precdbl(prec);
    2742                 :         20 :     get_nf_fp_compo(T, &F, NULL, 1, prec);
    2743         [ -  + ]:         20 :     if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",prec);
    2744                 :         20 :   }
    2745                 :        336 :   *u = d.u; return v;
    2746                 :            : }
    2747                 :            : 
    2748                 :            : GEN
    2749                 :        350 : polredabs0(GEN x, long flag)
    2750                 :            : {
    2751                 :        350 :   pari_sp av = avma;
    2752                 :            :   long i, l, vx;
    2753                 :            :   GEN y, a, u;
    2754                 :            :   nfbasic_t T;
    2755                 :            : 
    2756                 :        350 :   nfbasic_init(x, flag & nf_PARTIALFACT, &T);
    2757                 :        350 :   x = T.x; vx = varn(x);
    2758                 :            : 
    2759         [ +  + ]:        350 :   if (degpol(x) == 1)
    2760                 :            :   {
    2761                 :         14 :     u = NULL;
    2762                 :         14 :     y = mkvec( pol_x(vx) );
    2763                 :         14 :     a = mkvec( deg1pol_shallow(gen_1, negi(gel(x,2)), vx) );
    2764                 :         14 :     l = 2;
    2765                 :            :   }
    2766                 :            :   else
    2767                 :            :   {
    2768                 :            :     GEN v;
    2769         [ +  + ]:        336 :     if (!(flag & nf_PARTIALFACT)
    2770 [ +  + ][ -  + ]:         98 :         && T.dKP && lg(primes_certify(T.dK, T.dKP)) != 1) return gen_0;
    2771                 :        336 :     v = polredabs_aux(&T, &u);
    2772                 :        336 :     y = gel(v,1);
    2773                 :        336 :     a = gel(v,2); l = lg(a);
    2774         [ +  + ]:       6344 :     for (i=1; i<l; i++)
    2775         [ +  + ]:       6008 :       if (ZX_canon_neg(gel(y,i))) gel(a,i) = ZC_neg(gel(a,i));
    2776                 :        336 :     remove_duplicates(y,a);
    2777                 :        336 :     l = lg(a);
    2778         [ -  + ]:        336 :     if (l == 1)
    2779                 :          0 :       pari_err_BUG("polredabs (missing vector)");
    2780                 :            :   }
    2781         [ -  + ]:        350 :   if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
    2782         [ +  + ]:        350 :   if (flag & nf_ALL) {
    2783         [ +  + ]:       2135 :     for (i=1; i<l; i++) gel(y,i) = store(x, gel(y,i), gel(a,i), &T, flag, u);
    2784                 :            :   } else {
    2785                 :        343 :     GEN z = findmindisc(y, &a);
    2786                 :        343 :     y = store(x, z, a, &T, flag, u);
    2787                 :            :   }
    2788                 :        350 :   return gerepilecopy(av, y);
    2789                 :            : }
    2790                 :            : 
    2791                 :            : GEN
    2792                 :          0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
    2793                 :            : GEN
    2794                 :          0 : polredabs(GEN x) { return polredabs0(x,0); }
    2795                 :            : GEN
    2796                 :          0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
    2797                 :            : 
    2798                 :            : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
    2799                 :            :  * flag & nf_ORIG: + element (base change)
    2800                 :            :  * flag & nf_ABSOLUTE: absolute polynomial */
    2801                 :            : static GEN
    2802                 :         63 : rnfpolred_i(GEN nf, GEN relpol, long flag, long best)
    2803                 :            : {
    2804         [ +  + ]:         63 :   const char *f = best? "rnfpolredbest": "rnfpolredabs";
    2805 [ +  + ][ +  + ]:         63 :   const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
    2806                 :            :   pari_timer ti;
    2807                 :         63 :   GEN listP = NULL, red, bas, A, P, pol, T, rnfeq;
    2808                 :         63 :   long ty = typ(relpol);
    2809                 :         63 :   pari_sp av = avma;
    2810                 :            : 
    2811         [ +  + ]:         63 :   if (ty == t_VEC) {
    2812         [ -  + ]:         14 :     if (lg(relpol) != 3) pari_err_TYPE(f,relpol);
    2813                 :         14 :     listP = gel(relpol,2);
    2814                 :         14 :     relpol = gel(relpol,1);
    2815                 :            :   }
    2816         [ -  + ]:         63 :   if (typ(relpol) != t_POL) pari_err_TYPE(f,relpol);
    2817                 :         63 :   nf = checknf(nf);
    2818         [ -  + ]:         63 :   if (DEBUGLEVEL>1) timer_start(&ti);
    2819                 :         63 :   T = nf_get_pol(nf);
    2820                 :         63 :   relpol = RgX_nffix(f, T, relpol, 0);
    2821 [ +  + ][ +  + ]:         63 :   if (best || (flag & nf_PARTIALFACT))
    2822                 :            :   {
    2823         [ +  + ]:         42 :     if (abs)
    2824                 :            :     {
    2825                 :          7 :       rnfeq = nf_rnfeq(nf, relpol);
    2826                 :          7 :       pol = gel(rnfeq,1);
    2827                 :            :     }
    2828                 :            :     else
    2829                 :            :     {
    2830                 :            :       long sa;
    2831                 :         35 :       pol = rnfequationall(nf, relpol, &sa, NULL);
    2832                 :         35 :       rnfeq = mkvec5(gen_0,gen_0,stoi(sa),T,liftpol_shallow(relpol));
    2833                 :            :     }
    2834         [ +  + ]:         42 :     bas = listP? mkvec2(pol, listP): pol;
    2835         [ +  + ]:         84 :     if (best)
    2836                 :            :     {
    2837         [ +  + ]:         35 :       if (abs) red = polredbest(bas, 1);
    2838                 :            :       else
    2839                 :            :       {
    2840                 :            :         GEN ro, x, dx, a;
    2841                 :            :         nfbasic_t T;
    2842                 :         28 :         nfbasic_init(bas, nf_PARTIALFACT, &T);
    2843                 :         28 :         polredbest_aux(&T, &ro, &x, &dx, &a);
    2844                 :         28 :         red = mkvec2(x, a);
    2845                 :            :       }
    2846                 :            :     }
    2847                 :            :     else
    2848         [ -  + ]:          7 :       red = polredabs0(bas, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
    2849                 :            :   }
    2850                 :            :   else
    2851                 :            :   {
    2852                 :         21 :     GEN rnf = rnfinit(nf, relpol), M = rnf_basM(rnf);
    2853                 :         21 :     rnfeq = rnf_get_map(rnf);
    2854                 :         21 :     pol = rnf_get_polabs(rnf);
    2855                 :         21 :     bas = mkvec2(pol, RgM_to_RgXV(M, varn(pol)));
    2856         [ -  + ]:         21 :     if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
    2857                 :         21 :     red = polredabs0(bas, nf_RAW);
    2858                 :            :   }
    2859                 :         63 :   P = gel(red,1);
    2860                 :         63 :   A = gel(red,2);
    2861         [ -  + ]:         63 :   if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
    2862         [ +  + ]:         63 :   if (flag & nf_ABSOLUTE)
    2863                 :            :   {
    2864         [ +  + ]:         14 :     if (flag & nf_ORIG)
    2865                 :            :     {
    2866                 :          7 :       GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
    2867                 :          7 :       GEN k = gel(rnfeq,3); /* Mod(variable(relpol),relpol) + k*a root of pol */
    2868                 :          7 :       a = RgX_RgXQ_eval(a, lift_intern(A), P); /* Mod(a, P) root of T */
    2869                 :          7 :       P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
    2870                 :            :     }
    2871                 :         14 :     return gerepilecopy(av, P);
    2872                 :            :   }
    2873                 :         49 :   A = eltabstorel_lift(rnfeq, A);
    2874                 :         49 :   P = RgXQ_charpoly(A, relpol, varn(relpol));
    2875                 :         49 :   P = lift_if_rational(P);
    2876         [ +  + ]:         49 :   if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,relpol),P));
    2877                 :         63 :   return gerepilecopy(av, P);
    2878                 :            : }
    2879                 :            : GEN
    2880                 :         28 : rnfpolredabs(GEN nf, GEN relpol, long flag)
    2881                 :         28 : { return rnfpolred_i(nf,relpol,flag, 0); }
    2882                 :            : GEN
    2883                 :         35 : rnfpolredbest(GEN nf, GEN relpol, long flag)
    2884                 :            : {
    2885 [ +  - ][ -  + ]:         35 :   if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
    2886                 :         35 :   return rnfpolred_i(nf,relpol,flag, 1);
    2887                 :            : }

Generated by: LCOV version 1.9