Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19191-d7615ef) Lines: 1549 1665 93.0 %
Date: 2016-07-23 07:10:28 Functions: 117 132 88.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /**************************************************************/
      15             : /*                                                            */
      16             : /*                        NUMBER FIELDS                       */
      17             : /*                                                            */
      18             : /**************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : int new_galois_format = 0;
      23             : 
      24             : int
      25      211750 : checkrnf_i(GEN rnf)
      26      211750 : { return (typ(rnf)==t_VEC && lg(rnf)==13); }
      27             : 
      28             : void
      29      208656 : checkrnf(GEN rnf)
      30      208656 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
      31             : 
      32             : GEN
      33      686681 : checkbnf_i(GEN X)
      34             : {
      35      686681 :   if (typ(X) == t_VEC)
      36      686240 :     switch (lg(X))
      37             :     {
      38      684693 :       case 11: return X;
      39        1008 :       case 7:  return checkbnf_i(bnr_get_bnf(X));
      40             :     }
      41         980 :   return NULL;
      42             : }
      43             : 
      44             : GEN
      45     6452437 : checknf_i(GEN X)
      46             : {
      47     6452437 :   if (typ(X)==t_VEC)
      48     6451989 :     switch(lg(X))
      49             :     {
      50     6325607 :       case 10: return X;
      51      120110 :       case 11: return checknf_i(bnf_get_nf(X));
      52        5341 :       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        1372 :   return NULL;
      56             : }
      57             : 
      58             : GEN
      59      684546 : checkbnf(GEN x)
      60             : {
      61      684546 :   GEN bnf = checkbnf_i(x);
      62      684546 :   if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
      63      684546 :   return bnf;
      64             : }
      65             : 
      66             : GEN
      67     6324956 : checknf(GEN x)
      68             : {
      69     6324956 :   GEN nf = checknf_i(x);
      70     6324956 :   if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
      71     6324942 :   return nf;
      72             : }
      73             : 
      74             : void
      75      197230 : checkbnr(GEN bnr)
      76             : {
      77      197230 :   if (typ(bnr)!=t_VEC || lg(bnr)!=7)
      78           0 :     pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
      79      197230 :   (void)checkbnf(bnr_get_bnf(bnr));
      80      197230 : }
      81             : 
      82             : void
      83           0 : checkbnrgen(GEN bnr)
      84             : {
      85           0 :   checkbnr(bnr);
      86           0 :   if (lg(bnr_get_clgp(bnr))<=3)
      87           0 :     pari_err_TYPE("checkbnrgen [apply bnrinit(,,1), not bnrinit()]",bnr);
      88           0 : }
      89             : 
      90             : void
      91           0 : checksqmat(GEN x, long N)
      92             : {
      93           0 :   if (typ(x)!=t_MAT) pari_err_TYPE("checksqmat",x);
      94           0 :   if (lg(x) == 1 || lgcols(x) != N+1) pari_err_DIM("checksqmat");
      95           0 : }
      96             : 
      97             : GEN
      98      205987 : checkbid_i(GEN bid)
      99             : {
     100             :   GEN f;
     101      205987 :   if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(gel(bid,3)) != t_MAT) return NULL;
     102      199785 :   f = bid_get_mod(bid);
     103      199785 :   if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
     104      199785 :   return bid;
     105             : }
     106             : void
     107      199785 : checkbid(GEN bid)
     108             : {
     109      199785 :   if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
     110      199778 : }
     111             : void
     112       11403 : checkabgrp(GEN v)
     113             : {
     114       11403 :   if (typ(v) == t_VEC) switch(lg(v))
     115             :   {
     116       11312 :     case 4: if (typ(gel(v,3)) != t_VEC) break;
     117       11403 :     case 3: if (typ(gel(v,2)) != t_VEC) break;
     118       11375 :             if (typ(gel(v,1)) != t_INT) break;
     119       22750 :             return;/*OK*/
     120           0 :     default: break;
     121             :   }
     122          28 :   pari_err_TYPE("checkabgrp",v);
     123             : }
     124             : 
     125             : GEN
     126      369880 : checknfelt_mod(GEN nf, GEN x, const char *s)
     127             : {
     128      369880 :   GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
     129      369880 :   if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
     130      369810 :   return a;
     131             : }
     132             : 
     133             : void
     134        3234 : check_ZKmodule(GEN x, const char *s)
     135             : {
     136        3234 :   if (typ(x) != t_VEC || lg(x) < 3) pari_err_TYPE(s,x);
     137        3234 :   if (typ(gel(x,1)) != t_MAT) pari_err_TYPE(s,gel(x,1));
     138        3234 :   if (typ(gel(x,2)) != t_VEC) pari_err_TYPE(s,gel(x,2));
     139        3234 :   if (lg(gel(x,2)) != lgcols(x)) pari_err_DIM(s);
     140        3234 : }
     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         196 :   return typ_NULL;
     152             : }
     153             : 
     154             : GEN
     155       13041 : get_bnf(GEN x, long *t)
     156             : {
     157       13041 :   switch(typ(x))
     158             :   {
     159          56 :     case t_POL: *t = typ_POL;  return NULL;
     160          56 :     case t_QUAD: *t = typ_Q  ; return NULL;
     161             :     case t_VEC:
     162       12425 :       switch(lg(x))
     163             :       {
     164        4382 :         case 5: *t = typ_QUA; return NULL;
     165         322 :         case 6: *t = typv6(x); return NULL;
     166          77 :         case 7:  *t = typ_BNR;
     167          77 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     168          77 :           return x;
     169             :         case 9:
     170          56 :           x = gel(x,2);
     171          56 :           if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
     172          56 :           return NULL;
     173         217 :         case 10: *t = typ_NF; return NULL;
     174         189 :         case 11: *t = typ_BNF; return x;
     175          56 :         case 13: *t = typ_RNF; return NULL;
     176         266 :         case 17: *t = typ_ELL; return NULL;
     177             :       }
     178        6860 :       break;
     179             :     case t_COL:
     180         112 :       if (get_prid(x)) { *t = typ_MODPR; return NULL; }
     181          56 :       break;
     182             :   }
     183        7308 :   *t = typ_NULL; return NULL;
     184             : }
     185             : 
     186             : GEN
     187       62230 : get_nf(GEN x, long *t)
     188             : {
     189       62230 :   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       59815 :       switch(lg(x))
     195             :       {
     196             :         case 3:
     197         133 :           if (typ(gel(x,2)) != t_POLMOD) break;
     198         133 :           return get_nf(gel(x,1),t);
     199         133 :         case 5: *t = typ_QUA; return NULL;
     200       50302 :         case 6: *t = typv6(x); return NULL;
     201         140 :         case 7: *t = typ_BNR;
     202         140 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     203         140 :           x = bnf_get_nf(x);  if (typ(x)!=t_VEC || lg(x)!=10) break;
     204         140 :           return x;
     205             :         case 9:
     206         238 :           x = gel(x,2);
     207         238 :           if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
     208         238 :           return NULL;
     209         672 :         case 10: *t = typ_NF; return x;
     210        4788 :         case 11: *t = typ_BNF;
     211        4788 :           x = bnf_get_nf(x); if (typ(x)!=t_VEC || lg(x)!=10) break;
     212        4788 :           return x;
     213         336 :         case 13: *t = typ_RNF; return NULL;
     214        2940 :         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        2149 :   *t = typ_NULL; return NULL;
     222             : }
     223             : 
     224             : long
     225       38269 : nftyp(GEN x)
     226             : {
     227       38269 :   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       38241 :       switch(lg(x))
     233             :       {
     234         161 :         case 13: return typ_RNF;
     235             :         case 10:
     236       36456 :           if (typ(gel(x,1))!=t_POL) break;
     237       36449 :           return typ_NF;
     238             :         case 11:
     239          21 :           x = bnf_get_nf(x); if (typ(x)!=t_VEC || lg(x)!=10) break;
     240          21 :           return typ_BNF;
     241             :         case 7:
     242         896 :           x = bnr_get_bnf(x); if (typ(x)!=t_VEC || lg(x)!=11) break;
     243         889 :           x = bnf_get_nf(x);  if (typ(x)!=t_VEC || lg(x)!=10) break;
     244         889 :           return typ_BNR;
     245             :         case 6:
     246         672 :           return typv6(x);
     247             :         case 9:
     248           7 :           x = gel(x,2);
     249           7 :           if (typ(x) == t_VEC && lg(x) == 4) return typ_GAL;
     250          14 :         case 17: return typ_ELL;
     251             :       }
     252             :   }
     253          42 :   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        5950 : transroot(GEN x, int i, int j)
     265        5950 : { x = leafcopy(x); swap(gel(x,i), gel(x,j)); return x; }
     266             : 
     267             : GEN
     268        3157 : tschirnhaus(GEN x)
     269             : {
     270        3157 :   pari_sp av = avma, av2;
     271        3157 :   long a, v = varn(x);
     272        3157 :   GEN u, y = cgetg(5,t_POL);
     273             : 
     274        3157 :   if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
     275        3157 :   if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
     276        3157 :   if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
     277        3157 :   y[1] = evalsigne(1)|evalvarn(0);
     278             :   do
     279             :   {
     280        3283 :     a = random_bits(2); if (a==0) a  = 1; gel(y,4) = stoi(a);
     281        3283 :     a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
     282        3283 :     a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
     283        3283 :     u = RgXQ_charpoly(y,x,v); av2 = avma;
     284             :   }
     285        3283 :   while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
     286        3157 :   if (DEBUGLEVEL>1)
     287           0 :     err_printf("Tschirnhaus transform. New pol: %Ps",u);
     288        3157 :   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       12250 : ZX_Z_normalize(GEN pol, GEN *ptk)
     296             : {
     297       12250 :   long i,j, sk, n = degpol(pol); /* > 0 */
     298             :   GEN k, fa, P, E, a, POL;
     299             : 
     300       12250 :   a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
     301       46011 :   for (i = n-2; i >= 0; i--)
     302             :   {
     303       40775 :     k = gcdii(k, gel(a,i));
     304       40775 :     if (is_pm1(k)) { if (ptk) *ptk = gen_1; return pol; }
     305             :   }
     306        5236 :   sk = signe(k);
     307        5236 :   if (!sk) { if (ptk) *ptk = gen_1; return pol; /* monomial! */ }
     308        4256 :   fa = absi_factor_limit(k, 0); k = gen_1;
     309        4256 :   P = gel(fa,1);
     310        4256 :   E = gel(fa,2);
     311        4256 :   POL = leafcopy(pol); a = POL+2;
     312        9723 :   for (i = lg(P)-1; i > 0; i--)
     313             :   {
     314        5467 :     GEN p = gel(P,i), pv, pvj;
     315        5467 :     long vmin = itos(gel(E,i));
     316             :     /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
     317       38206 :     for (j=n-1; j>=0; j--)
     318             :     {
     319             :       long v;
     320       32739 :       if (!signe(gel(a,j))) continue;
     321       20804 :       v = Z_pval(gel(a,j), p) / (n - j);
     322       20804 :       if (v < vmin) vmin = v;
     323             :     }
     324        5467 :     if (!vmin) continue;
     325        1134 :     pvj = pv = powiu(p,vmin); k = mulii(k, pv);
     326             :     /* a[j] /= p^(v*(n-j)) */
     327        8085 :     for (j=n-1; j>=0; j--)
     328             :     {
     329        6951 :       if (j < n-1) pvj = mulii(pvj, pv);
     330        6951 :       gel(a,j) = diviiexact(gel(a,j), pvj);
     331             :     }
     332             :   }
     333        4256 :   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       12257 : ZX_primitive_to_monic(GEN pol, GEN *pL)
     341             : {
     342       12257 :   long i,j, n = degpol(pol);
     343       12257 :   GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
     344             : 
     345       12257 :   if (is_pm1(lc))
     346             :   {
     347       11998 :     if (pL) *pL = gen_1;
     348       11998 :     return signe(lc) < 0? ZX_neg(pol): pol;
     349             :   }
     350         259 :   if (signe(lc) < 0)
     351          35 :     POL = ZX_neg(pol);
     352             :   else
     353         224 :     POL = leafcopy(pol);
     354         259 :   a = POL+2; lc = gel(a,n);
     355         259 :   fa = Z_factor_limit(lc,0); L = gen_1;
     356         259 :   P = gel(fa,1);
     357         259 :   E = gel(fa,2);
     358         651 :   for (i = lg(P)-1; i > 0; i--)
     359             :   {
     360         392 :     GEN p = gel(P,i), pk, pku;
     361         392 :     long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
     362             : 
     363         392 :     if (d < 0) { k++; d += n; }
     364             :     /* k = ceil(e[i] / n); find d, k such that  p^d pol(x / p^k) monic */
     365        1519 :     for (j=n-1; j>0; j--)
     366             :     {
     367        1127 :       if (!signe(gel(a,j))) continue;
     368        1015 :       v = Z_pval(gel(a,j), p);
     369        1015 :       while (v + d < k * j) { k++; d += n; }
     370             :     }
     371         392 :     pk = powiu(p,k); j0 = d/k;
     372         392 :     L = mulii(L, pk);
     373             : 
     374         392 :     pku = powiu(p,d - k*j0);
     375             :     /* a[j] *= p^(d - kj) */
     376        1617 :     for (j=j0; j>=0; j--)
     377             :     {
     378        1225 :       if (j < j0) pku = mulii(pku, pk);
     379        1225 :       gel(a,j) = mulii(gel(a,j), pku);
     380             :     }
     381         392 :     j0++;
     382         392 :     pku = powiu(p,k*j0 - d);
     383             :     /* a[j] /= p^(kj - d) */
     384        1078 :     for (j=j0; j<=n; j++)
     385             :     {
     386         686 :       if (j > j0) pku = mulii(pku, pk);
     387         686 :       gel(a,j) = diviiexact(gel(a,j), pku);
     388             :     }
     389             :   }
     390         259 :   if (pL) *pL = L;
     391         259 :   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       11998 : ZX_Q_normalize(GEN pol, GEN *pL)
     399             : {
     400       11998 :   GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
     401       11998 :   POL = ZX_Z_normalize(POL, pL);
     402       11998 :   if (pL) *pL = gdiv(lc, *pL);
     403       11998 :   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        8736 :   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        1379 : roots_to_ZX(GEN z, long *e)
     428             : {
     429        1379 :   GEN a = roots_to_pol(z,0);
     430        1379 :   GEN b = grndtoi(real_i(a),e);
     431        1379 :   long e1 = gexpo(imag_i(a));
     432        1379 :   if (e1 > *e) *e = e1;
     433        1379 :   return b;
     434             : }
     435             : 
     436             : static GEN
     437        1057 : polgaloisnames(long a, long b)
     438             : {
     439        1057 :   const char * const t[]={"S1", "S2", "A3", "S3",
     440             :        "C(4) = 4", "E(4) = 2[x]2", "D(4)", "A4", "S4",
     441             :        "C(5) = 5", "D(5) = 5:2", "F(5) = 5:4", "A5", "S5",
     442             :        "C(6) = 6 = 3[x]2", "D_6(6) = [3]2", "D(6) = S(3)[x]2",
     443             :              "A_4(6) = [2^2]3", "F_18(6) = [3^2]2 = 3 wr 2",
     444             :              "2A_4(6) = [2^3]3 = 2 wr 3", "S_4(6d) = [2^2]S(3)",
     445             :              "S_4(6c) = 1/2[2^3]S(3)", "F_18(6):2 = [1/2.S(3)^2]2",
     446             :              "F_36(6) = 1/2[S(3)^2]2", "2S_4(6) = [2^3]S(3) = 2 wr S(3)",
     447             :              "L(6) = PSL(2,5) = A_5(6)", "F_36(6):2 = [S(3)^2]2 = S(3) wr 2",
     448             :              "L(6):2 = PGL(2,5) = S_5(6)", "A6", "S6",
     449             :        "C(7) = 7", "D(7) = 7:2", "F_21(7) = 7:3", "F_42(7) = 7:6",
     450             :              "L(7) = L(3,2)", "A7", "S7"};
     451             : 
     452        1057 :    const long idx[]={0,1,2,4,9,14,30};
     453        1057 :    return strtoGENstr(t[idx[a-1]+b-1]);
     454             : }
     455             : 
     456             : static GEN
     457        1057 : galois_res(long d, long n, long s, long k)
     458             : {
     459        1057 :   long kk = k;
     460        1057 :   GEN z = cgetg(5,t_VEC);
     461        1057 :   if (!new_galois_format)
     462             :   {
     463         553 :     switch (d) {
     464             :       case 6:
     465         238 :         kk = (k == 6 || k == 2)? 2: 1;
     466         238 :         break;
     467             :       default:
     468         315 :         kk = 1;
     469             :     }
     470             :   }
     471        1057 :   gel(z,1) = stoi(n);
     472        1057 :   gel(z,2) = stoi(s);
     473        1057 :   gel(z,3) = stoi(kk);
     474        1057 :   gel(z,4) = polgaloisnames(d,k);
     475        1057 :   return z;
     476             : }
     477             : 
     478             : GEN
     479        3059 : polgalois(GEN x, long prec)
     480             : {
     481        3059 :   pari_sp av = avma, av1;
     482             :   long i,j,k,n,f,l,l2,e,e1,pr,ind;
     483             :   GEN x1,p1,p2,p3,p4,p5,w,z,ee;
     484        3059 :   const int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
     485        3059 :   const int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
     486             :                       1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
     487             :                       1,2,4,6, 1,5,2,4, 1,3,2,6, 1,2,3,5, 1,4,2,3};
     488        3059 :   if (typ(x)!=t_POL) pari_err_TYPE("galois",x);
     489        3059 :   n=degpol(x);
     490        3059 :   if (n>11) pari_err_IMPL("galois of degree higher than 11");
     491        3059 :   x = Q_primpart(x);
     492        3059 :   RgX_check_ZX(x, "galois");
     493        3059 :   if (!ZX_is_irred(x)) pari_err_IRREDPOL("galois",x);
     494             : 
     495        3059 :   if (n<4)
     496             :   {
     497          91 :     if (n == 1) { avma = av; return galois_res(n,1, 1,1); }
     498          91 :     if (n == 2) { avma = av; return galois_res(n,2,-1,1); }
     499             :     /* n = 3 */
     500          63 :     f = Z_issquare(ZX_disc(x));
     501          63 :     avma = av;
     502          63 :     return f? galois_res(n,3,1,1):
     503             :               galois_res(n,6,-1,2);
     504             :   }
     505        2968 :   x1 = x = ZX_Q_normalize(x,NULL); av1=avma;
     506        2968 :   if (n > 7) return galoisbig(x, prec);
     507             :   for(;;)
     508             :   {
     509        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         168 :       case 5: z = cgetg(7,t_VEC);
     543         168 :         ee= cgetg(7,t_VECSMALL);
     544         168 :         w = cgetg(7,t_VECSMALL);
     545         168 :         prec = nbits2prec((long)(fb*21.) + 64);
     546             :         for(;;)
     547             :         {
     548             :           for(;;)
     549             :           {
     550         168 :             p1=QX_complex_roots(x,prec);
     551        1176 :             for (l=1; l<=6; l++)
     552             :             {
     553        1008 :               p2=(l==1)?p1: ((l<6)?transroot(p1,1,l): transroot(p1,2,5));
     554        1008 :               p3=gen_0;
     555        6048 :               for (k=0,i=1; i<=5; i++,k+=4)
     556             :               {
     557       10080 :                 p5 = gadd(gmul(gel(p2,ind5[k]),gel(p2,ind5[k+1])),
     558       10080 :                           gmul(gel(p2,ind5[k+2]),gel(p2,ind5[k+3])));
     559        5040 :                 p3 = gadd(p3, gmul(gsqr(gel(p2,i)),p5));
     560             :               }
     561        1008 :               gel(w,l) = grndtoi(real_i(p3),&e);
     562        1008 :               e1 = gexpo(imag_i(p3)); if (e1>e) e=e1;
     563        1008 :               ee[l]=e; gel(z,l) = p3;
     564             :             }
     565         168 :             p5 = roots_to_ZX(z, &e); if (e <= -10) break;
     566           0 :             prec = precdbl(prec);
     567           0 :           }
     568         168 :           if (!ZX_is_squarefree(p5)) goto tchi;
     569         168 :           p3=gel(ZX_factor(p5),1);
     570         168 :           f=Z_issquare(ZX_disc(x));
     571         168 :           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         105 :           if (!f) { avma = av; return galois_res(n,20,-1,3); }
     577             : 
     578          77 :           pr = - (prec2nbits(prec) >> 1);
     579         280 :           for (l=1; l<=6; l++)
     580         280 :             if (ee[l] <= pr && gequal0(poleval(p5,gel(w,l)))) break;
     581          77 :           if (l>6) pari_err_BUG("galois (bug4)");
     582          77 :           p2=(l==6)? transroot(p1,2,5):transroot(p1,1,l);
     583          77 :           p3=gen_0;
     584         462 :           for (i=1; i<=5; i++)
     585             :           {
     586         385 :             j = (i == 5)? 1: i+1;
     587         770 :             p3 = gadd(p3,gmul(gmul(gel(p2,i),gel(p2,j)),
     588         770 :                               gsub(gel(p2,j),gel(p2,i))));
     589             :           }
     590          77 :           p5=gsqr(p3); p4=grndtoi(real_i(p5),&e);
     591          77 :           e1 = gexpo(imag_i(p5)); if (e1>e) e=e1;
     592          77 :           if (e <= -10)
     593             :           {
     594          77 :             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       59808 :               for (i=1; i<=5; i++) for (j=i+1; j<=6; j++)
     613             :               {
     614      112140 :                 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         336 :               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        3024 :                   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         108 :                 case 3: f = Z_issquare(ZX_disc(gel(p2,1)))
     668          56 :                          || Z_issquare(ZX_disc(gel(p2,2)));
     669          56 :                   avma = av;
     670          56 :                   return f? galois_res(n,18,-1,5): galois_res(n,36,-1,9);
     671             :               }
     672             :             case 3:
     673         504 :               for (l2=1; l2<=3; l2++)
     674         378 :                 if (degpol(gel(p2,l2)) >= 3) p3 = gel(p2,l2);
     675         126 :               if (degpol(p3) == 3)
     676             :               {
     677          63 :                 f = Z_issquare(ZX_disc(p3)); avma = av;
     678          63 :                 return f? galois_res(n,6,-1,1): galois_res(n,12,-1,3);
     679             :               }
     680             :               else
     681             :               {
     682          63 :                 f = Z_issquare(ZX_disc(x)); avma = av;
     683          63 :                 return f? galois_res(n,12,1,4): galois_res(n,24,-1,8);
     684             :               }
     685          28 :             case 4: avma = av; return galois_res(n,6,-1,2);
     686           0 :             default: pari_err_BUG("galois (bug3)");
     687             :           }
     688           0 :         }
     689             : 
     690         203 :       case 7: z = cgetg(36,t_VEC);
     691         203 :         prec = nbits2prec((long)(fb*7.) + 64);
     692             :         for(;;)
     693             :         {
     694         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        4410 :               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          91 :           prec = precdbl(prec);
     703          91 :         }
     704         203 :         if (!ZX_is_squarefree(p5)) goto tchi;
     705         203 :         p2=gel(ZX_factor(p5),1);
     706         203 :         switch(lg(p2)-1)
     707             :         {
     708          56 :           case 1: f = Z_issquare(ZX_disc(x)); avma = av;
     709          56 :             return f? galois_res(n,2520,1,6): galois_res(n,5040,-1,7);
     710          63 :           case 2: f = degpol(gel(p2,1)); avma = av;
     711          63 :             return (f==7 || f==28)? galois_res(n,168,1,5): galois_res(n,42,-1,4);
     712          28 :           case 3: avma = av; return galois_res(n,21,1,3);
     713          28 :           case 4: avma = av; return galois_res(n,14,-1,2);
     714          28 :           case 5: avma = av; return galois_res(n,7,1,1);
     715           0 :           default: pari_err_BUG("galois (bug2)");
     716             :         }
     717             :     }
     718         210 :     tchi: avma = av1; x = tschirnhaus(x1);
     719         210 :   }
     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       19544 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
     740             : {
     741       19544 :   pari_sp av = avma;
     742       19544 :   long i = lg(pol)-1;
     743             :   GEN res, den;
     744       19544 :   if (i==1) return gen_0;
     745       19544 :   pol = Q_remove_denom(pol, &den);
     746       19544 :   res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
     747      104986 :   for (i-- ; i>=2; i--)
     748       85442 :     res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
     749       19544 :   if (den) res = RgC_Rg_div(res, den);
     750       19544 :   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        2807 :   return QX_table_nfpoleval(nf, x, m);
     777             : }
     778             : 
     779             : /* compute s(x), not stack clean */
     780             : static GEN
     781       10080 : ZC_galoisapply(GEN nf, GEN s, GEN x)
     782             : {
     783       10080 :   x = nf_to_scalar_or_alg(nf, x);
     784       10080 :   if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
     785       10017 :   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           7 :   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           0 :   return NULL; /* not reached */
     851             : }
     852             : 
     853             : GEN
     854        5047 : galoisapply(GEN nf, GEN aut, GEN x)
     855             : {
     856        5047 :   pari_sp av = avma;
     857             :   long lx, j;
     858             :   GEN y;
     859             : 
     860        5047 :   nf = checknf(nf);
     861        5047 :   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        3535 :       aut = algtobasis(nf, aut);
     886        3535 :       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         931 :       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           0 :   return NULL; /* not reached */
     898             : }
     899             : 
     900             : GEN
     901        1897 : nfgaloismatrix(GEN nf, GEN s)
     902             : {
     903             :   GEN zk, M, m;
     904             :   long k, l;
     905        1897 :   nf = checknf(nf);
     906        1897 :   zk = nf_get_zk(nf);
     907        1897 :   if (typ(s) != t_COL) s = algtobasis(nf, s); /* left on stack for efficiency */
     908        1897 :   m = zk_multable(nf, s);
     909        1897 :   l = lg(s); M = cgetg(l, t_MAT);
     910        1897 :   gel(M, 1) = col_ei(l-1, 1); /* s(1) = 1 */
     911        8617 :   for (k = 2; k < l; k++)
     912        6720 :     gel(M, k) = QX_table_nfpoleval(nf, gel(zk, k), m);
     913        1897 :   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           0 :   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        1400 :   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           0 :   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           0 :   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          21 :   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          14 :   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       14704 : get_nfpol(GEN x, GEN *nf)
    1159             : {
    1160       14704 :   if (typ(x) == t_POL) { *nf = NULL; return x; }
    1161        8628 :   *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          21 :   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        9063 : get_roots(GEN x, long r1, long prec)
    1280             : {
    1281             :   long i, ru;
    1282             :   GEN z;
    1283        9063 :   if (typ(x) != t_POL)
    1284             :   {
    1285           0 :     z = leafcopy(x);
    1286           0 :     ru = (lg(z)-1 + r1) >> 1;
    1287             :   }
    1288             :   else
    1289             :   {
    1290        9063 :     long n = degpol(x);
    1291        9063 :     z = (r1 == n)? realroots(x, NULL, prec): QX_complex_roots(x,prec);
    1292        9063 :     ru = (n+r1)>>1;
    1293             :   }
    1294        9063 :   for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
    1295        9063 :   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       58919 : quicktrace(GEN x, GEN sym)
    1307             : {
    1308       58919 :   GEN p1 = gen_0;
    1309             :   long i;
    1310             : 
    1311       58919 :   if (typ(x) != t_POL) return gmul(x, gel(sym,1));
    1312       58919 :   if (signe(x))
    1313             :   {
    1314       58919 :     sym--;
    1315      834127 :     for (i=lg(x)-1; i>1; i--)
    1316      775208 :       p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
    1317             :   }
    1318       58919 :   return p1;
    1319             : }
    1320             : 
    1321             : static GEN
    1322        5075 : get_Tr(GEN mul, GEN x, GEN basden)
    1323             : {
    1324        5075 :   GEN t, bas = gel(basden,1), den = gel(basden,2);
    1325        5075 :   long i, j, n = lg(bas)-1;
    1326        5075 :   GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
    1327             : 
    1328        5075 :   gel(TW,1) = utoipos(n);
    1329       19341 :   for (i=2; i<=n; i++)
    1330             :   {
    1331       14266 :     t = quicktrace(gel(bas,i), sym);
    1332       14266 :     if (den && den[i]) t = diviiexact(t,gel(den,i));
    1333       14266 :     gel(TW,i) = t; /* tr(w[i]) */
    1334             :   }
    1335        5075 :   gel(T,1) = TW;
    1336       19341 :   for (i=2; i<=n; i++)
    1337             :   {
    1338       14266 :     gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
    1339      117579 :     for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
    1340      103313 :       gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
    1341             :   }
    1342        5075 :   return T;
    1343             : }
    1344             : 
    1345             : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
    1346             : GEN
    1347       13016 : get_bas_den(GEN bas)
    1348             : {
    1349       13016 :   GEN b,d,den, dbas = leafcopy(bas);
    1350       13016 :   long i, l = lg(bas);
    1351       13016 :   int power = 1;
    1352       13016 :   den = cgetg(l,t_VEC);
    1353       67138 :   for (i=1; i<l; i++)
    1354             :   {
    1355       54122 :     b = Q_remove_denom(gel(bas,i), &d);
    1356       54122 :     gel(dbas,i) = b;
    1357       54122 :     gel(den,i) = d; if (d) power = 0;
    1358             :   }
    1359       13016 :   if (power) den = NULL; /* power basis */
    1360       13016 :   return mkvec2(dbas, den);
    1361             : }
    1362             : 
    1363             : /* Internal: nf partially filled. Require pol; fill zk, invzk, multable */
    1364             : void
    1365        5075 : nf_set_multable(GEN nf, GEN bas, GEN basden)
    1366             : {
    1367        5075 :   GEN T = nf_get_pol(nf), invbas, basM;
    1368        5075 :   long i,j, n = degpol(T);
    1369        5075 :   GEN w, den, mul = cgetg(n*n+1,t_MAT);
    1370             : 
    1371        5075 :   if (typ(bas) == t_MAT)
    1372           0 :   { basM = bas; bas = RgM_to_RgXV(basM, varn(T)); }
    1373             :   else
    1374        5075 :     basM = RgV_to_RgM(bas, n);
    1375        5075 :   gel(nf,7) = bas;
    1376        5075 :   gel(nf,8) = invbas = QM_inv(basM, gen_1);
    1377        5075 :   gel(nf,9) = mul;
    1378             : 
    1379        5075 :   if (!basden) basden = get_bas_den(nf_get_zk(nf)); /*integral basis*/
    1380        5075 :   w   = gel(basden,1);
    1381        5075 :   den = gel(basden,2);
    1382             :   /* i = 1 split for efficiency, assume w[1] = 1 */
    1383       24416 :   for (j=1; j<=n; j++)
    1384       19341 :     gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
    1385       19341 :   for (i=2; i<=n; i++)
    1386      117579 :     for (j=i; j<=n; j++)
    1387             :     {
    1388      103313 :       pari_sp av = avma;
    1389      103313 :       GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
    1390      103313 :       z = mulmat_pol(invbas, z); /* integral column */
    1391      103313 :       if (den)
    1392             :       {
    1393       71127 :         GEN d = mul_denom(gel(den,i), gel(den,j));
    1394       71127 :         if (d) z = ZC_Z_divexact(z, d);
    1395             :       }
    1396      103313 :       gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
    1397             :     }
    1398        5075 : }
    1399             : 
    1400             : /* as get_Tr, mul_table not precomputed */
    1401             : static GEN
    1402        2716 : make_Tr(GEN x, GEN basden)
    1403             : {
    1404        2716 :   long i,j, n = degpol(x);
    1405             :   GEN c, t, d;
    1406        2716 :   GEN T   = cgetg(n+1,t_MAT);
    1407        2716 :   GEN sym = polsym(x, n-1);
    1408        2716 :   GEN w   = gel(basden,1); /* W[i] = w[i]/den[i] */
    1409        2716 :   GEN den = gel(basden,2);
    1410             :   /* assume W[1] = 1, case i = 1 split for efficiency */
    1411        2716 :   c = cgetg(n+1,t_COL); gel(T,1) = c;
    1412        2716 :   gel(c, 1) = utoipos(n);
    1413        8183 :   for (j=2; j<=n; j++)
    1414             :   {
    1415        5467 :     pari_sp av = avma;
    1416        5467 :     t = quicktrace(gel(w,j), sym);
    1417        5467 :     if (den)
    1418             :     {
    1419        3962 :       d = gel(den,j);
    1420        3962 :       if (d) t = diviiexact(t, d);
    1421             :     }
    1422        5467 :     gel(c,j) = gerepileuptoint(av, t);
    1423             :   }
    1424        8183 :   for (i=2; i<=n; i++)
    1425             :   {
    1426        5467 :     c = cgetg(n+1,t_COL); gel(T,i) = c;
    1427        5467 :     for (j=1; j<i ; j++) gel(c,j) = gcoeff(T,i,j);
    1428       36512 :     for (   ; j<=n; j++)
    1429             :     {
    1430       31045 :       pari_sp av = avma;
    1431       31045 :       t = (i == j)? ZXQ_sqr(gel(w,i), x): ZXQ_mul(gel(w,i),gel(w,j), x);
    1432       31045 :       t = quicktrace(t, sym);
    1433       31045 :       if (den)
    1434             :       {
    1435       28441 :         d = mul_denom(gel(den,i),gel(den,j));
    1436       28441 :         if (d) t = diviiexact(t, d);
    1437             :       }
    1438       31045 :       gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
    1439             :     }
    1440             :   }
    1441        2716 :   return T;
    1442             : }
    1443             : 
    1444             : /* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */
    1445             : static void
    1446       10341 : get_roots_for_M(nffp_t *F)
    1447             : {
    1448             :   long n, eBD, PREC;
    1449             : 
    1450       10341 :   if (F->extraprec < 0)
    1451             :   { /* not initialized yet */
    1452             :     double er;
    1453       10341 :     n = degpol(F->x);
    1454       10341 :     eBD = 1 + gexpo(gel(F->basden,1));
    1455       10341 :     er  = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->x);
    1456       10341 :     if (er < 0) er = 0;
    1457       10341 :     F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
    1458             :   }
    1459             : 
    1460       10341 :   PREC = F->prec + F->extraprec;
    1461             :   /* make sure that default accuracy is the same on 32/64bit */
    1462             : #ifndef LONG_IS_64BIT
    1463        1515 :   if (odd(PREC)) PREC += EXTRAPRECWORD;
    1464             : #endif
    1465       20682 :   if (F->ro && gprecision(gel(F->ro,1)) >= PREC) return;
    1466        9063 :   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       10341 : make_M(nffp_t *F, int trunc)
    1472             : {
    1473       10341 :   GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
    1474             :   GEN m, d, M;
    1475       10341 :   long i, j, l = lg(ro), n = lg(bas);
    1476       10341 :   M = cgetg(n,t_MAT);
    1477       10341 :   gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
    1478       10341 :   for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
    1479       39351 :   for (i=1; i<l; i++)
    1480             :   {
    1481       29010 :     GEN r = gel(ro,i), ri;
    1482       29010 :     ri = (gexpo(r) > 1)? ginv(r): NULL;
    1483       29010 :     for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
    1484             :   }
    1485       10341 :   if (den)
    1486       28233 :     for (j=2; j<n; j++)
    1487             :     {
    1488       23870 :       d = gel(den,j); if (!d) continue;
    1489       19368 :       m = gel(M,j);
    1490       19368 :       for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
    1491             :     }
    1492             : 
    1493       10341 :   if (trunc && gprecision(M) > F->prec)
    1494             :   {
    1495        1857 :     M     = gprec_w(M, F->prec);
    1496        1857 :     F->ro = gprec_w(ro,F->prec);
    1497             :   }
    1498       10341 :   F->M = M;
    1499       10341 : }
    1500             : 
    1501             : /* return G real such that G~ * G = T_2 */
    1502             : static void
    1503       10341 : make_G(nffp_t *F)
    1504             : {
    1505       10341 :   GEN G, M = F->M;
    1506       10341 :   long i, j, k, r1 = F->r1, l = lg(M);
    1507             : 
    1508       10341 :   G = cgetg(l, t_MAT);
    1509       56390 :   for (j=1; j<l; j++)
    1510             :   {
    1511       46049 :     GEN g = cgetg(l, t_COL);
    1512       46049 :     GEN m = gel(M,j);
    1513       46049 :     gel(G,j) = g;
    1514       46049 :     for (k=i=1; i<=r1; i++) g[k++] = m[i];
    1515      266003 :     for (     ; k < l; i++)
    1516             :     {
    1517      219954 :       GEN r = gel(m,i);
    1518      219954 :       if (typ(r) == t_COMPLEX)
    1519             :       {
    1520      202915 :         gel(g,k++) = mpadd(gel(r,1), gel(r,2));
    1521      202915 :         gel(g,k++) = mpsub(gel(r,1), gel(r,2));
    1522             :       }
    1523             :       else
    1524             :       {
    1525       17039 :         gel(g,k++) = r;
    1526       17039 :         gel(g,k++) = r;
    1527             :       }
    1528             :     }
    1529             :   }
    1530       10341 :   F->G = G;
    1531       10341 : }
    1532             : 
    1533             : static void
    1534       10341 : make_M_G(nffp_t *F, int trunc)
    1535             : {
    1536       10341 :   get_roots_for_M(F);
    1537       10341 :   make_M(F, trunc);
    1538       10341 :   make_G(F);
    1539       10341 : }
    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        9589 : nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
    1554             : {
    1555        9589 :   F->x  = T->x;
    1556        9589 :   F->ro = ro;
    1557        9589 :   F->r1 = T->r1;
    1558        9589 :   if (!T->basden) T->basden = get_bas_den(T->bas);
    1559        9589 :   F->basden = T->basden;
    1560        9589 :   F->extraprec = -1;
    1561        9589 :   F->prec = prec;
    1562        9589 : }
    1563             : 
    1564             : static void
    1565        6208 : get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, int trunc, long prec)
    1566             : {
    1567        6208 :   nffp_init(F,T,ro,prec);
    1568        6208 :   make_M_G(F, trunc);
    1569        6208 : }
    1570             : 
    1571             : static GEN
    1572        5075 : get_sign(long r1, long n) { return mkvec2s(r1, (n-r1)>>1); }
    1573             : 
    1574             : GEN
    1575        5075 : nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec)
    1576             : {
    1577        5075 :   GEN nf = cgetg(10,t_VEC);
    1578        5075 :   GEN x = T->x, absdK, Tr, D, TI, A, dA, MDI, mat = cgetg(9,t_VEC);
    1579        5075 :   long n = degpol(T->x);
    1580             :   nffp_t F;
    1581        5075 :   get_nf_fp_compo(T, &F, ro, 0, prec);
    1582             : 
    1583        5075 :   gel(nf,1) = T->x;
    1584        5075 :   gel(nf,2) = get_sign(T->r1, n);
    1585        5075 :   gel(nf,3) = T->dK;
    1586        5075 :   gel(nf,4) = T->index;
    1587        5075 :   gel(nf,6) = F.ro;
    1588        5075 :   gel(nf,5) = mat;
    1589             : 
    1590        5075 :   gel(mat,1) = F.M;
    1591        5075 :   gel(mat,2) = F.G;
    1592             : 
    1593        5075 :   nf_set_multable(nf, T->bas, F.basden);
    1594             : 
    1595        5075 :   Tr = get_Tr(gel(nf,9), x, F.basden);
    1596        5075 :   absdK = T->dK; if (signe(absdK) < 0) absdK = negi(absdK);
    1597        5075 :   TI = ZM_inv(Tr, absdK); /* dK T^-1 */
    1598        5075 :   A = Q_primitive_part(TI, &dA);
    1599        5075 :   gel(mat,6) = A; /* primitive part of codifferent, dA its content */
    1600        5075 :   dA = dA? diviiexact(absdK, dA): absdK;
    1601        5075 :   A = ZM_hnfmodid(A, dA);
    1602        5075 :   MDI = idealtwoelt(nf, A);
    1603        5075 :   gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
    1604        5075 :   gel(mat,7) = MDI;
    1605        5075 :   if (is_pm1(T->index)) /* principal ideal (x'), whose norm is |dK| */
    1606             :   {
    1607        3528 :     D = zk_scalar_or_multable(nf, ZX_deriv(x));
    1608        3528 :     if (typ(D) == t_MAT) D = ZM_hnfmod(D, absdK);
    1609             :   }
    1610             :   else
    1611        1547 :     D = RgM_Rg_mul(idealinv(nf, A), dA);
    1612        5075 :   gel(mat,3) = RM_round_maxrank(F.G);
    1613        5075 :   gel(mat,4) = Tr;
    1614        5075 :   gel(mat,5) = D;
    1615        5075 :   gel(mat,8) = T->dKP? shallowtrans(T->dKP): cgetg(1,t_VEC);
    1616        5075 :   return nf;
    1617             : }
    1618             : 
    1619             : static GEN
    1620          91 : primes_certify(GEN dK, GEN dKP)
    1621             : {
    1622          91 :   long i, l = lg(dKP);
    1623          91 :   GEN v, w, D = dK;
    1624          91 :   v = vectrunc_init(l);
    1625          91 :   w = vectrunc_init(l);
    1626         427 :   for (i = 1; i < l; i++)
    1627             :   {
    1628         336 :     GEN p = gel(dKP,i);
    1629         336 :     vectrunc_append(isprime(p)? w: v, p);
    1630         336 :     (void)Z_pvalrem(D, p, &D);
    1631             :   }
    1632          91 :   if (!is_pm1(D))
    1633             :   {
    1634           0 :     if (signe(D) < 0) D = negi(D);
    1635           0 :     vectrunc_append(isprime(D)? w: v, D);
    1636             :   }
    1637          91 :   return mkvec2(v,w);
    1638             : }
    1639             : GEN
    1640           7 : nfcertify(GEN nf)
    1641             : {
    1642           7 :   pari_sp av = avma;
    1643             :   GEN vw;
    1644           7 :   nf = checknf(nf);
    1645           7 :   vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
    1646           7 :   return gerepilecopy(av, gel(vw,1));
    1647             : }
    1648             : 
    1649             : #if 0 /* used to check benches between HNF nf.zk and LLL-reduced nf.zk */
    1650             : static GEN
    1651             : hnffromLLL(GEN nf)
    1652             : {
    1653             :   GEN d, x;
    1654             :   x = RgV_to_RgM(nf_get_zk(nf), nf_get_degree(nf));
    1655             :   x = Q_remove_denom(x, &d);
    1656             :   if (!d) return x; /* power basis */
    1657             :   return RgM_solve(ZM_hnfmodid(x, d), x);
    1658             : }
    1659             : 
    1660             : static GEN
    1661             : nfbasechange(GEN u, GEN x)
    1662             : {
    1663             :   long i,lx;
    1664             :   GEN y;
    1665             :   switch(typ(x))
    1666             :   {
    1667             :     case t_COL: /* nfelt */
    1668             :       return RgM_RgC_mul(u, x);
    1669             : 
    1670             :     case t_MAT: /* ideal */
    1671             :       y = cgetg_copy(x, &lx);
    1672             :       for (i=1; i<lx; i++) gel(y,i) = RgM_RgC_mul(u, gel(x,i));
    1673             :       break;
    1674             : 
    1675             :     case t_VEC: /* pr */
    1676             :       checkprid(x); y = leafcopy(x);
    1677             :       gel(y,2) = RgM_RgC_mul(u, gel(y,2));
    1678             :       gel(y,5) = RgM_RgC_mul(u, gel(y,5));
    1679             :       break;
    1680             :     default: y = x;
    1681             :   }
    1682             :   return y;
    1683             : }
    1684             : 
    1685             : GEN
    1686             : nffromhnfbasis(GEN nf, GEN x)
    1687             : {
    1688             :   long tx = typ(x);
    1689             :   pari_sp av = avma;
    1690             :   GEN u;
    1691             :   if (!is_vec_t(tx)) return gcopy(x);
    1692             :   nf = checknf(nf);
    1693             :   u = hnffromLLL(nf);
    1694             :   return gerepilecopy(av, nfbasechange(u, x));
    1695             : }
    1696             : 
    1697             : GEN
    1698             : nftohnfbasis(GEN nf, GEN x)
    1699             : {
    1700             :   long tx = typ(x);
    1701             :   pari_sp av = avma;
    1702             :   GEN u;
    1703             :   if (!is_vec_t(tx)) return gcopy(x);
    1704             :   nf = checknf(nf);
    1705             :   u = ZM_inv(hnffromLLL(nf), gen_1);
    1706             :   return gerepilecopy(av, nfbasechange(u, x));
    1707             : }
    1708             : #endif
    1709             : 
    1710             : /* set *pro to roots of T->x */
    1711             : static GEN
    1712        3381 : get_red_G(nfbasic_t *T, GEN *pro)
    1713             : {
    1714        3381 :   GEN G, u, u0 = NULL;
    1715             :   pari_sp av;
    1716        3381 :   long i, prec, n = degpol(T->x);
    1717             :   nffp_t F;
    1718             : 
    1719        3381 :   prec = nbits2prec(n+32);
    1720        3381 :   nffp_init(&F, T, NULL, prec);
    1721        3381 :   av = avma;
    1722        3381 :   for (i=1; ; i++)
    1723             :   {
    1724        3381 :     F.prec = prec; make_M_G(&F, 0); G = F.G;
    1725        3381 :     if (u0) G = RgM_mul(G, u0);
    1726        3381 :     if (DEBUGLEVEL)
    1727           0 :       err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
    1728           0 :                   prec + F.extraprec, prec, F.extraprec);
    1729        3381 :     if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST|LLL_COMPATIBLE)))
    1730             :     {
    1731        3381 :       if (lg(u)-1 == n) break;
    1732             :       /* singular ==> loss of accuracy */
    1733           0 :       if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
    1734           0 :       else    u0 = gerepilecopy(av, u);
    1735             :     }
    1736           0 :     prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
    1737           0 :     F.ro = NULL;
    1738           0 :     if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
    1739           0 :   }
    1740        3381 :   if (u0) u = RgM_mul(u0,u);
    1741        3381 :   *pro = F.ro; return u;
    1742             : }
    1743             : 
    1744             : /* Compute an LLL-reduced basis for the integer basis of nf(T).
    1745             :  * set *pro = roots of x if computed [NULL if not computed] */
    1746             : static void
    1747        6097 : set_LLL_basis(nfbasic_t *T, GEN *pro, double DELTA)
    1748             : {
    1749        6097 :   GEN B = T->bas;
    1750        6097 :   if (T->r1 == degpol(T->x)) {
    1751        2716 :     pari_sp av = avma;
    1752        2716 :     GEN u, basden = T->basden;
    1753        2716 :     if (!basden) basden = get_bas_den(B);
    1754        2716 :     u = ZM_lll(make_Tr(T->x,basden), DELTA,
    1755             :                LLL_GRAM|LLL_KEEP_FIRST|LLL_IM|LLL_COMPATIBLE);
    1756        2716 :     B = gerepileupto(av, RgV_RgM_mul(B, u));
    1757        2716 :     *pro = NULL;
    1758             :   }
    1759             :   else
    1760        3381 :     B = RgV_RgM_mul(B, get_red_G(T, pro));
    1761        6097 :   T->bas = B;
    1762        6097 :   T->basden = get_bas_den(B);
    1763        6097 : }
    1764             : 
    1765             : static int
    1766        3172 : cmp_abs_ZX(GEN x, GEN y) { return gen_cmp_RgX((void*)&absi_cmp, x, y); }
    1767             : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
    1768             :  * (if so update *dx) */
    1769             : static int
    1770        4621 : ZX_is_better(GEN y, GEN x, GEN *dx)
    1771             : {
    1772        4621 :   GEN d = ZX_disc(y);
    1773             :   int cmp;
    1774        4621 :   if (!*dx) *dx = ZX_disc(x);
    1775        4621 :   cmp = absi_cmp(d, *dx);
    1776        4621 :   if (cmp < 0) { *dx = d; return 1; }
    1777        3781 :   if (cmp == 0) return cmp_abs_ZX(y, x) < 0;
    1778         609 :   return 0;
    1779             : }
    1780             : 
    1781             : static void polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
    1782             : /* Seek a simpler, polynomial pol defining the same number field as
    1783             :  * x (assumed to be monic at this point) */
    1784             : static GEN
    1785          84 : nfpolred(nfbasic_t *T, GEN *pro)
    1786             : {
    1787          84 :   GEN x = T->x, dx, b, rev;
    1788          84 :   long n = degpol(x), v = varn(x);
    1789             : 
    1790          84 :   if (n == 1) {
    1791           7 :     T->x = deg1pol_shallow(gen_1, gen_m1, v);
    1792           7 :     *pro = NULL; return pol_1(v);
    1793             :   }
    1794          77 :   polredbest_aux(T, pro, &x, &dx, &b);
    1795          77 :   if (x == T->x) return NULL; /* no improvement */
    1796          56 :   if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
    1797             : 
    1798             :   /* update T */
    1799          56 :   rev = QXQ_reverse(b, T->x);
    1800          56 :   T->bas = QXV_QXQ_eval(T->bas, rev, x);
    1801          56 :   (void)Z_issquareall(diviiexact(dx,T->dK), &(T->index));
    1802          56 :   T->basden = get_bas_den(T->bas);
    1803          56 :   T->dx = dx;
    1804          56 :   T->x = x;
    1805          56 :   *pro = NULL; /* reset */
    1806          56 :   return rev;
    1807             : }
    1808             : 
    1809             : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
    1810             :  * basis. Assume bas[1] is 1 and that the leading coefficient of elements
    1811             :  * of bas are of the form 1/b for a t_INT b */
    1812             : GEN
    1813         525 : get_nfindex(GEN bas)
    1814             : {
    1815         525 :   pari_sp av = avma;
    1816         525 :   long n = lg(bas)-1, i;
    1817             :   GEN D, d, mat;
    1818             : 
    1819         525 :   D = gen_1; /* assume bas[1] = 1 */
    1820        1197 :   for (i = 2; i <= n; i++)
    1821             :   { /* in most cases [e.g after nfbasis] basis is upper triangular! */
    1822         924 :     GEN B = gel(bas,i), lc;
    1823         924 :     if (degpol(B) != i-1) break;
    1824         672 :     lc = gel(B, i+1);
    1825         672 :     switch (typ(lc))
    1826             :     {
    1827         399 :       case t_INT: continue;
    1828         273 :       case t_FRAC: lc = gel(lc,2); break;
    1829           0 :       default: pari_err_TYPE("get_nfindex",lc);
    1830             :     }
    1831         273 :     D = mulii(D, lc);
    1832             :   }
    1833         525 :   if (i <= n)
    1834             :   { /* not triangular after all */
    1835         252 :     bas = Q_remove_denom(bas, &d);
    1836         252 :     if (!d) { avma = av; return D; }
    1837         238 :     mat = RgV_to_RgM(bas, n);
    1838         238 :     d = diviiexact(powiu(d, n), ZM_det(mat));
    1839         238 :     D = mulii(D,absi(d));
    1840             :   }
    1841         511 :   return gerepileuptoint(av, D);
    1842             : }
    1843             : 
    1844             : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
    1845             :  * basis (deprecated), or listP data (nfbasis input format) to specify
    1846             :  * a set of primes at with the basis order must be maximal.
    1847             :  * 1) nf type (or unrecognized): return t_VEC
    1848             :  * 2) ZX or [ZX, listP]: return t_POL
    1849             :  * 3) [ZX, order basis]: return 0 (deprecated)
    1850             :  * incorrect: return -1 */
    1851             : static long
    1852        6181 : nf_input_type(GEN x)
    1853             : {
    1854             :   GEN T, V;
    1855             :   long i, d, v;
    1856        6181 :   switch(typ(x))
    1857             :   {
    1858        5544 :     case t_POL: return t_POL;
    1859             :     case t_VEC:
    1860         637 :       if (lg(x) != 3) return t_VEC; /* nf or incorrect */
    1861         630 :       T = gel(x,1); V = gel(x,2);
    1862         630 :       if (typ(T) != t_POL) return -1;
    1863         630 :       switch(typ(V))
    1864             :       {
    1865          35 :         case t_INT: case t_MAT: return t_POL;
    1866             :         case t_VEC: case t_COL:
    1867         595 :           if (RgV_is_ZV(V)) return t_POL;
    1868         574 :           break;
    1869           0 :         default: return -1;
    1870             :       }
    1871         574 :       d = degpol(T); v = varn(T);
    1872         574 :       if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
    1873        3192 :       for (i = 1; i <= d; i++)
    1874             :       { /* check integer basis */
    1875        2639 :         GEN c = gel(V,i);
    1876        2639 :         switch(typ(c))
    1877             :         {
    1878          28 :           case t_INT: break;
    1879        2611 :           case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
    1880             :           /* fall through */
    1881          14 :           default: return -1;
    1882             :         }
    1883             :       }
    1884         553 :       return 0;
    1885             :   }
    1886           0 :   return t_VEC; /* nf or incorrect */
    1887             : }
    1888             : 
    1889             : static void
    1890        5068 : nfbasic_add_disc(nfbasic_t *T)
    1891             : {
    1892        5068 :   if (!T->index) T->index = get_nfindex(T->bas);
    1893        5068 :   if (!T->dx) T->dx = ZX_disc(T->x);
    1894        5068 :   if (!T->dK) T->dK = diviiexact(T->dx, sqri(T->index));
    1895        5068 : }
    1896             : 
    1897             : static void
    1898        6181 : nfbasic_init(GEN x, long flag, nfbasic_t *T)
    1899             : {
    1900        6181 :   GEN bas, dK, dx, index, unscale = gen_1;
    1901             :   long r1;
    1902             : 
    1903        6181 :   T->dKP = NULL;
    1904        6181 :   switch (nf_input_type(x))
    1905             :   {
    1906             :     case t_POL:
    1907             :     {
    1908             :       nfmaxord_t S;
    1909        5600 :       nfmaxord(&S, x, flag);
    1910        5593 :       x = S.T;
    1911        5593 :       T->x0 = S.T0;
    1912        5593 :       T->dKP = S.dKP;
    1913        5593 :       dK = S.dK;
    1914        5593 :       index = S.index;
    1915        5593 :       bas = S.basis;
    1916        5593 :       dx = S.dT;
    1917        5593 :       unscale = S.unscale;
    1918        5593 :       r1 = ZX_sturm(x);
    1919        5593 :       break;
    1920             :     }
    1921             :     case t_VEC:
    1922             :     { /* nf, bnf, bnr */
    1923           7 :       GEN nf = checknf(x);
    1924           7 :       x     = nf_get_pol(nf);
    1925           7 :       dK    = nf_get_disc(nf);
    1926           7 :       index = nf_get_index(nf);
    1927           7 :       bas   = nf_get_zk(nf);
    1928           7 :       T->x0 = x;
    1929           7 :       dx = NULL;
    1930           7 :       r1 = nf_get_r1(nf);
    1931           7 :       break;
    1932             :     }
    1933             :     case 0: /* monic integral polynomial + integer basis */
    1934         553 :       bas = gel(x,2); x = gel(x,1);
    1935         553 :       T->x0 = x;
    1936         553 :       index = NULL;
    1937         553 :       dx = NULL;
    1938         553 :       dK = NULL;
    1939         553 :       r1 = ZX_sturm(x);
    1940         553 :       break;
    1941             :     default: /* -1 */
    1942          21 :       pari_err_TYPE("nfbasic_init", x);
    1943        6153 :       return;
    1944             :   }
    1945        6153 :   T->x     = x;
    1946        6153 :   T->unscale = unscale;
    1947        6153 :   T->r1    = r1;
    1948        6153 :   T->dx    = dx;
    1949        6153 :   T->dK    = dK;
    1950        6153 :   T->bas   = bas;
    1951        6153 :   T->basden= NULL;
    1952        6153 :   T->index = index;
    1953             : }
    1954             : 
    1955             : void
    1956        5089 : nfinit_step1(nfbasic_t *T, GEN x, long flag)
    1957             : {
    1958        5089 :   nfbasic_init(x, flag, T);
    1959        5068 :   if (!ZX_is_irred(T->x)) pari_err_IRREDPOL("nfinit",x);
    1960        5068 : }
    1961             : 
    1962             : GEN
    1963        5068 : nfinit_step2(nfbasic_t *T, long flag, long prec)
    1964             : {
    1965             :   GEN nf, unscale;
    1966             : 
    1967        5068 :   if (!(flag & nf_RED) && !equali1(leading_coeff(T->x0)))
    1968             :   {
    1969          49 :     pari_warn(warner,"non-monic polynomial. Result of the form [nf,c]");
    1970          49 :     flag |= nf_RED | nf_ORIG;
    1971             :   }
    1972        5068 :   unscale = T->unscale;
    1973        5068 :   if (!(flag & nf_RED) && !isint1(unscale))
    1974             :   { /* implies lc(x0) = 1 and L := 1/unscale is integral */
    1975         105 :     long d = degpol(T->x0);
    1976         105 :     GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
    1977         105 :     GEN f= powiu(L, (d*(d-1)) >> 1);
    1978         105 :     T->x = T->x0; /* restore original user-supplied x0, unscale data */
    1979         105 :     T->unscale = gen_1;
    1980         105 :     T->dx    = gmul(T->dx, sqri(f));
    1981         105 :     T->bas   = RgXV_unscale(T->bas, unscale);
    1982         105 :     T->index = gmul(T->index, f);
    1983             :   }
    1984        5068 :   nfbasic_add_disc(T); /* more expensive after set_LLL_basis */
    1985        5068 :   if (flag & nf_RED)
    1986             :   {
    1987             :     GEN ro, rev;
    1988             :     /* lie to polred: more efficient to update *after* modreverse, than to
    1989             :      * unscale in the polred subsystem */
    1990          84 :     T->unscale = gen_1;
    1991          84 :     rev = nfpolred(T, &ro);
    1992          84 :     nf = nfbasic_to_nf(T, ro, prec);
    1993          84 :     if (flag & nf_ORIG)
    1994             :     {
    1995          56 :       if (!rev) rev = pol_x(varn(T->x)); /* no improvement */
    1996          56 :       if (!isint1(unscale)) rev = RgX_Rg_div(rev, unscale);
    1997          56 :       nf = mkvec2(nf, mkpolmod(rev, T->x));
    1998             :     }
    1999          84 :     T->unscale = unscale; /* restore */
    2000             :   } else {
    2001        4984 :     GEN ro; set_LLL_basis(T, &ro, 0.99);
    2002        4984 :     nf = nfbasic_to_nf(T, ro, prec);
    2003             :   }
    2004        5068 :   return nf;
    2005             : }
    2006             : /* Initialize the number field defined by the polynomial x (in variable v)
    2007             :  * flag & nf_RED:     try a polred first.
    2008             :  * flag & nf_ORIG
    2009             :  *    do a polred and return [nfinit(x), Mod(a,red)], where
    2010             :  *    Mod(a,red) = Mod(v,x) (i.e return the base change). */
    2011             : GEN
    2012        3094 : nfinitall(GEN x, long flag, long prec)
    2013             : {
    2014        3094 :   const pari_sp av = avma;
    2015             :   nfbasic_t T;
    2016             :   GEN nf;
    2017             : 
    2018        3094 :   if (checkrnf_i(x)) return check_and_build_nfabs(x, prec);
    2019        3087 :   nfinit_step1(&T, x, flag);
    2020        3066 :   nf = nfinit_step2(&T, flag, prec);
    2021        3066 :   return gerepilecopy(av, nf);
    2022             : }
    2023             : 
    2024             : GEN
    2025           0 : nfinitred(GEN x, long prec)  { return nfinitall(x, nf_RED, prec); }
    2026             : GEN
    2027           0 : nfinitred2(GEN x, long prec) { return nfinitall(x, nf_RED|nf_ORIG, prec); }
    2028             : GEN
    2029        1183 : nfinit(GEN x, long prec)     { return nfinitall(x, 0, prec); }
    2030             : 
    2031             : GEN
    2032        1694 : nfinit0(GEN x, long flag,long prec)
    2033             : {
    2034        1694 :   switch(flag)
    2035             :   {
    2036             :     case 0:
    2037        1673 :     case 1: return nfinitall(x,0,prec);
    2038          14 :     case 2: case 4: return nfinitall(x,nf_RED,prec);
    2039           7 :     case 3: case 5: return nfinitall(x,nf_RED|nf_ORIG,prec);
    2040           0 :     default: pari_err_FLAG("nfinit");
    2041             :   }
    2042           0 :   return NULL; /* not reached */
    2043             : }
    2044             : 
    2045             : /* assume x a bnr/bnf/nf */
    2046             : long
    2047       68350 : nf_get_prec(GEN x)
    2048             : {
    2049       68350 :   GEN nf = checknf(x), ro = nf_get_roots(nf);
    2050       68350 :   return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
    2051             : }
    2052             : 
    2053             : /* assume nf is an nf */
    2054             : GEN
    2055         752 : nfnewprec_shallow(GEN nf, long prec)
    2056             : {
    2057         752 :   GEN NF = leafcopy(nf);
    2058             :   nffp_t F;
    2059         752 :   gel(NF,5) = leafcopy(gel(NF,5));
    2060         752 :   remake_GM(NF, &F, prec);
    2061         752 :   gel(NF,6) = F.ro;
    2062         752 :   gmael(NF,5,1) = F.M;
    2063         752 :   gmael(NF,5,2) = F.G;
    2064         752 :   return NF;
    2065             : }
    2066             : 
    2067             : GEN
    2068          63 : nfnewprec(GEN nf, long prec)
    2069             : {
    2070             :   GEN z;
    2071          63 :   switch(nftyp(nf))
    2072             :   {
    2073          49 :     default: pari_err_TYPE("nfnewprec", nf);
    2074           7 :     case typ_BNF: z = bnfnewprec(nf,prec); break;
    2075           7 :     case typ_BNR: z = bnrnewprec(nf,prec); break;
    2076             :     case typ_NF: {
    2077           0 :       pari_sp av = avma;
    2078           0 :       z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
    2079           0 :       break;
    2080             :     }
    2081             :   }
    2082          14 :   return z;
    2083             : }
    2084             : 
    2085             : /********************************************************************/
    2086             : /**                                                                **/
    2087             : /**                           POLRED                               **/
    2088             : /**                                                                **/
    2089             : /********************************************************************/
    2090             : GEN
    2091           0 : embednorm_T2(GEN x, long r1)
    2092             : {
    2093           0 :   pari_sp av = avma;
    2094           0 :   GEN p = RgV_sumpart(x, r1);
    2095           0 :   GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
    2096           0 :   if (q != gen_0) p = gadd(p, gmul2n(q,1));
    2097           0 :   return avma == av? gcopy(p): gerepileupto(av, p);
    2098             : }
    2099             : 
    2100             : /* simplified version of gnorm for scalar, non-complex inputs, without GC */
    2101             : static GEN
    2102        6671 : real_norm(GEN x)
    2103             : {
    2104        6671 :   switch(typ(x))
    2105             :   {
    2106           0 :     case t_INT:  return sqri(x);
    2107        6671 :     case t_REAL: return sqrr(x);
    2108           0 :     case t_FRAC: return sqrfrac(x);
    2109             :   }
    2110           0 :   pari_err_TYPE("real_norm", x);
    2111           0 :   return NULL;
    2112             : }
    2113             : /* simplified version of gnorm, without GC */
    2114             : static GEN
    2115     1880831 : complex_norm(GEN x)
    2116             : {
    2117     1880831 :   return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
    2118             : }
    2119             : /* return T2(x), argument r1 needed in case x has components whose type
    2120             :  * is unexpected, e.g. all of them t_INT for embed(gen_1) */
    2121             : GEN
    2122        1583 : embed_T2(GEN x, long r1)
    2123             : {
    2124        1583 :   pari_sp av = avma;
    2125        1583 :   long i, l = lg(x);
    2126        1583 :   GEN c, s = NULL, t = NULL;
    2127        1583 :   if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
    2128        8254 :   for (i = 1; i <= r1; i++)
    2129             :   {
    2130        6671 :     c = real_norm(gel(x,i));
    2131        6671 :     s = s? gadd(s, c): c;
    2132             :   }
    2133        6338 :   for (; i < l; i++)
    2134             :   {
    2135        4755 :     c = complex_norm(gel(x,i));
    2136        4755 :     t = t? gadd(t, c): c;
    2137             :   }
    2138        1583 :   if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
    2139        1583 :   return gerepileupto(av, s);
    2140             : }
    2141             : /* return N(x) */
    2142             : GEN
    2143     1017715 : embed_norm(GEN x, long r1)
    2144             : {
    2145     1017715 :   pari_sp av = avma;
    2146     1017715 :   long i, l = lg(x);
    2147     1017715 :   GEN c, s = NULL, t = NULL;
    2148     1017715 :   if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
    2149     2691295 :   for (i = 1; i <= r1; i++)
    2150             :   {
    2151     1677050 :     c = gel(x,i);
    2152     1677050 :     s = s? gmul(s, c): c;
    2153             :   }
    2154     2890321 :   for (; i < l; i++)
    2155             :   {
    2156     1876076 :     c = complex_norm(gel(x,i));
    2157     1876076 :     t = t? gmul(t, c): c;
    2158             :   }
    2159     1014245 :   if (t) s = s? gmul(s,t): t;
    2160     1014245 :   return gerepileupto(av, s);
    2161             : }
    2162             : 
    2163             : typedef struct {
    2164             :   long r1, v, prec;
    2165             :   GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
    2166             :   GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
    2167             :   GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
    2168             :   GEN bound; /* T2 norm of the polynomial defining nf */
    2169             :   long expo_best_disc; /* expo(disc(x)), best generator so far */
    2170             : } CG_data;
    2171             : 
    2172             : /* characteristic pol of x (given by embeddings) */
    2173             : static GEN
    2174       24778 : get_pol(CG_data *d, GEN x)
    2175             : {
    2176             :   long e;
    2177       24778 :   GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
    2178       24778 :   return (e > -5)? NULL: g;
    2179             : }
    2180             : 
    2181             : /* characteristic pol of x (given as vector on (w_i)) */
    2182             : static GEN
    2183       10593 : get_polchar(CG_data *d, GEN x)
    2184       10593 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
    2185             : 
    2186             : /* Choose a canonical polynomial in the pair { z(X), (+/-)z(-X) }.
    2187             :  * z a ZX with lc(z) > 0. We want to keep that property, while
    2188             :  * ensuring that the leading coeff of the odd (resp. even) part of z is < 0
    2189             :  * if deg z is even (resp. odd).
    2190             :  * Either leave z alone (return 1) or set z <-- (-1)^deg(z) z(-X). In place. */
    2191             : static int
    2192       10733 : ZX_canon_neg(GEN z)
    2193             : {
    2194             :   long i,s;
    2195             : 
    2196      165778 :   for (i = lg(z)-2; i >= 2; i -= 2)
    2197             :   { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
    2198       77576 :     s = signe(gel(z,i));
    2199       77576 :     if (!s) continue;
    2200             :     /* non trivial */
    2201        5420 :     if (s < 0) break; /* the condition is already satisfied */
    2202             : 
    2203        1807 :     for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
    2204        1807 :     return 1;
    2205             :   }
    2206        8926 :   return 0;
    2207             : }
    2208             : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
    2209             :  * Return NULL on failure: discriminant too large or non primitive */
    2210             : static GEN
    2211       20272 : try_polmin(CG_data *d, nfbasic_t *T, GEN v, long flag, GEN *ai)
    2212             : {
    2213       20272 :   const long best = flag & nf_ABSOLUTE;
    2214             :   long ed;
    2215       20272 :   pari_sp av = avma;
    2216             :   GEN g;
    2217       20272 :   if (best)
    2218             :   {
    2219       19397 :     ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
    2220       19397 :     avma = av; if (d->expo_best_disc < ed) return NULL;
    2221             :   }
    2222             :   else
    2223         875 :     ed = 0;
    2224       11683 :   g = get_pol(d, v);
    2225             :   /* accuracy too low, compute algebraically */
    2226       11683 :   if (!g) { avma = av; g = ZXQ_charpoly(*ai, T->x, varn(T->x)); }
    2227       11683 :   (void)ZX_gcd_all(g, ZX_deriv(g), &g);
    2228       11683 :   if (best && degpol(g) != degpol(T->x)) { avma = av; return NULL; }
    2229        4718 :   g = gerepilecopy(av, g);
    2230        4718 :   d->expo_best_disc = ed;
    2231        4718 :   if (flag & nf_ORIG)
    2232             :   {
    2233        1001 :     if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
    2234        1001 :     if (!isint1(T->unscale)) *ai = RgX_unscale(*ai, T->unscale);
    2235             :   }
    2236             :   else
    2237        3717 :     (void)ZX_canon_neg(g);
    2238        4718 :   if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
    2239        4718 :   return g;
    2240             : }
    2241             : 
    2242             : /* does x generate the correct field ? */
    2243             : static GEN
    2244       10593 : chk_gen(void *data, GEN x)
    2245             : {
    2246       10593 :   pari_sp av = avma, av1;
    2247       10593 :   GEN h, g = get_polchar((CG_data*)data,x);
    2248       10593 :   if (!g) pari_err_PREC("chk_gen");
    2249       10593 :   av1 = avma;
    2250       10593 :   h = ZX_gcd(g, ZX_deriv(g));
    2251       10593 :   if (degpol(h)) { avma = av; return NULL; }
    2252        6050 :   if (DEBUGLEVEL>3) err_printf("  generator: %Ps\n",g);
    2253        6050 :   avma = av1; return gerepileupto(av, g);
    2254             : }
    2255             : 
    2256             : static long
    2257        1456 : chk_gen_prec(long N, long bit)
    2258        1456 : { return nbits2prec(10 + (long)log2((double)N) + bit); }
    2259             : 
    2260             : /* Remove duplicate polynomials in P, updating A (same indices), in place.
    2261             :  * Among elements having the same characteristic pol, choose the smallest
    2262             :  * according to ZV_abscmp */
    2263             : static void
    2264         343 : remove_duplicates(GEN P, GEN A)
    2265             : {
    2266         343 :   long k, i, l = lg(P);
    2267         343 :   pari_sp av = avma;
    2268             :   GEN x, a;
    2269             : 
    2270         686 :   if (l < 2) return;
    2271         343 :   (void)sort_factor_pol(mkmat2(P, A), cmpii);
    2272         343 :   x = gel(P,1); a = gel(A,1);
    2273        6015 :   for  (k=1,i=2; i<l; i++)
    2274        5672 :     if (ZX_equal(gel(P,i), x))
    2275             :     {
    2276        3480 :       if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
    2277             :     }
    2278             :     else
    2279             :     {
    2280        2192 :       gel(A,k) = a;
    2281        2192 :       gel(P,k) = x;
    2282        2192 :       k++;
    2283        2192 :       x = gel(P,i); a = gel(A,i);
    2284             :     }
    2285         343 :   l = k+1;
    2286         343 :   gel(A,k) = a; setlg(A,l);
    2287         343 :   gel(P,k) = x; setlg(P,l); avma = av;
    2288             : }
    2289             : 
    2290             : static long
    2291        1113 : polred_init(nfbasic_t *T, nffp_t *F, CG_data *d)
    2292             : {
    2293        1113 :   long e, prec, n = degpol(T->x);
    2294             :   double log2rho;
    2295             :   GEN ro;
    2296        1113 :   set_LLL_basis(T, &ro, 0.9999);
    2297             :   /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
    2298        1113 :   log2rho = ro ? (double)gexpo(ro): fujiwara_bound(T->x);
    2299        1113 :   e = n * (long)(log2rho + log2((double)n)) + 1;
    2300        1113 :   if (e < 0) e = 0; /* can occur if n = 1 */
    2301        1113 :   prec = chk_gen_prec(n, e);
    2302        1113 :   get_nf_fp_compo(T, F, ro, 1, prec);
    2303        1113 :   d->v = varn(T->x);
    2304        1113 :   d->expo_best_disc = -1;
    2305        1113 :   d->ZKembed = NULL;
    2306        1113 :   d->M = NULL;
    2307        1113 :   d->u = NULL;
    2308        1113 :   d->r1= T->r1; return prec;
    2309             : }
    2310             : static GEN
    2311         350 : findmindisc(GEN y, GEN *pa)
    2312             : {
    2313         350 :   GEN a = *pa, x = gel(y,1), b = gel(a,1), dx = NULL;
    2314         350 :   long i, l = lg(y);
    2315         421 :   for (i = 2; i < l; i++)
    2316             :   {
    2317          71 :     GEN yi = gel(y,i);
    2318          71 :     if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
    2319             :   }
    2320         350 :   *pa = b; return x;
    2321             : }
    2322             : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
    2323             :  * [ the best wrt discriminant ordering ], but keep all non-primitive
    2324             :  * polynomials */
    2325             : static void
    2326         770 : filter(GEN y, GEN b, long n)
    2327             : {
    2328             :   GEN x, a, dx;
    2329         770 :   long i, k = 1, l = lg(y);
    2330         770 :   a = x = dx = NULL;
    2331        5551 :   for (i = 1; i < l; i++)
    2332             :   {
    2333        4781 :     GEN yi = gel(y,i), ai = gel(b,i);
    2334        4781 :     if (degpol(yi) == n)
    2335             :     {
    2336        4613 :       pari_sp av = avma;
    2337        4613 :       if (dx && !ZX_is_better(yi,x,&dx)) { avma = av; continue; }
    2338        1085 :       if (!dx) dx = ZX_disc(yi);
    2339        1085 :       x = yi; a = ai; continue;
    2340             :     }
    2341         168 :     gel(y,k) = yi;
    2342         168 :     gel(b,k) = ai; k++;
    2343             :   }
    2344         770 :   if (dx)
    2345             :   {
    2346         770 :     gel(y,k) = x;
    2347         770 :     gel(b,k) = a; k++;
    2348             :   }
    2349         770 :   setlg(y, k);
    2350         770 :   setlg(b, k);
    2351         770 : }
    2352             : 
    2353             : static GEN
    2354         798 : polred_aux(nfbasic_t *T, GEN *pro, long flag)
    2355             : { /* only keep polynomials of max degree and best discriminant */
    2356         798 :   const long best = flag & nf_ABSOLUTE;
    2357         798 :   const long orig = flag & nf_ORIG;
    2358         798 :   GEN M, b, y, x = T->x;
    2359         798 :   long maxi, i, j, k, v = varn(x), n = lg(T->bas)-1;
    2360             :   nffp_t F;
    2361             :   CG_data d;
    2362             : 
    2363         798 :   if (n == 1)
    2364             :   {
    2365          28 :     if (!best)
    2366             :     {
    2367          14 :       GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2368          14 :       return orig? mkmat2(mkcol(ch),mkcol(gen_1)): mkvec(ch);
    2369             :     }
    2370             :     else
    2371          14 :       return orig? trivial_fact(): cgetg(1,t_VEC);
    2372             :   }
    2373             : 
    2374         770 :   (void)polred_init(T, &F, &d);
    2375         770 :   *pro = F.ro;
    2376         770 :   M = F.M;
    2377         770 :   if (best)
    2378             :   {
    2379         707 :     if (!T->dx) T->dx = ZX_disc(T->x);
    2380         707 :     d.expo_best_disc = expi(T->dx);
    2381             :   }
    2382             : 
    2383             :   /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
    2384         770 :   y = cgetg(n*n + 1, t_VEC);
    2385         770 :   b = cgetg(n*n + 1, t_COL);
    2386         770 :   k = 1;
    2387         770 :   if (!best)
    2388             :   {
    2389          63 :     GEN ch = deg1pol_shallow(gen_1, gen_m1, v);
    2390          63 :     gel(y,1) = ch; gel(b,1) = gen_1; k++;
    2391             :   }
    2392        4270 :   for (i = 2; i <= n; i++)
    2393             :   {
    2394             :     GEN ch, ai;
    2395        3500 :     ai = gel(T->bas,i);
    2396        3500 :     ch = try_polmin(&d, T, gel(M,i), flag, &ai);
    2397        3500 :     if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2398             :   }
    2399         770 :   maxi = minss(n, 3);
    2400        2884 :   for (i = 1; i <= maxi; i++)
    2401       10500 :     for (j = i+1; j <= n; j++)
    2402             :     {
    2403             :       GEN ch, ai, v;
    2404        8386 :       ai = gadd(gel(T->bas,i), gel(T->bas,j));
    2405        8386 :       v = RgV_add(gel(M,i), gel(M,j));
    2406             :       /* defining polynomial for Q(w_i+w_j) */
    2407        8386 :       ch = try_polmin(&d, T, v, flag, &ai);
    2408        8386 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2409             : 
    2410        8386 :       ai = gsub(gel(T->bas,i), gel(T->bas,j));
    2411        8386 :       v = RgV_sub(gel(M,i), gel(M,j));
    2412             :       /* defining polynomial for Q(w_i-w_j) */
    2413        8386 :       ch = try_polmin(&d, T, v, flag, &ai);
    2414        8386 :       if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
    2415             :     }
    2416         770 :   setlg(y, k);
    2417         770 :   setlg(b, k); filter(y, b, n);
    2418         770 :   if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
    2419         147 :   (void)sort_factor_pol(mkmat2(y, b), cmpii);
    2420         147 :   settyp(y, t_COL); return mkmat2(b, y);
    2421             : }
    2422             : 
    2423             : static GEN
    2424          84 : Polred(GEN x, long flag, GEN fa)
    2425             : {
    2426          84 :   pari_sp av = avma;
    2427             :   GEN ro;
    2428          84 :   nfbasic_t T; nfbasic_init(fa? mkvec2(x,fa): x, flag & nf_PARTIALFACT, &T);
    2429          77 :   return gerepilecopy(av, polred_aux(&T, &ro, flag));
    2430             : }
    2431             : 
    2432             : /* finds "best" polynomial in polred_aux list, defaulting to T->x if none of
    2433             :  * them is primitive. *px is the ZX, characteristic polynomial of Mod(*pb,T->x),
    2434             :  * *pdx its discriminant. Set *pro = polroots(T->x) [ NOT *px ]. */
    2435             : static void
    2436         721 : polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
    2437             : {
    2438         721 :   GEN y, x = T->x; /* default value */
    2439             :   long i, l;
    2440         721 :   y = polred_aux(T, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
    2441         721 :   *pdx = T->dx;
    2442         721 :   if (pb)
    2443             :   {
    2444         133 :     GEN a, b = deg1pol_shallow(T->unscale, gen_0, varn(x));
    2445         133 :     a = gel(y,1); l = lg(a);
    2446         133 :     y = gel(y,2);
    2447         259 :     for (i=1; i<l; i++)
    2448             :     {
    2449         126 :       GEN yi = gel(y,i);
    2450         126 :       pari_sp av = avma;
    2451         126 :       if (ZX_is_better(yi,x,pdx)) { x = yi; b = gel(a,i); } else avma = av;
    2452             :     }
    2453         133 :     *pb = b;
    2454             :   }
    2455             :   else
    2456             :   {
    2457         588 :     l = lg(y);
    2458        1169 :     for (i=1; i<l; i++)
    2459             :     {
    2460         581 :       GEN yi = gel(y,i);
    2461         581 :       pari_sp av = avma;
    2462         581 :       if (ZX_is_better(yi,x,pdx)) x = yi; else avma = av;
    2463             :     }
    2464             :   }
    2465         721 :   if (!*pdx) *pdx = ZX_disc(x);
    2466         721 :   *px = x;
    2467         721 : }
    2468             : GEN
    2469         616 : polredbest(GEN T0, long flag)
    2470             : {
    2471         616 :   pari_sp av = avma;
    2472             :   GEN T, dT, ro, a;
    2473             :   nfbasic_t S;
    2474         616 :   if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
    2475         616 :   T = T0; nfbasic_init(T, nf_PARTIALFACT, &S);
    2476         616 :   polredbest_aux(&S, &ro, &T, &dT, flag? &a: NULL);
    2477         616 :   if (flag)
    2478             :   { /* charpoly(Mod(a,T0)) = T */
    2479             :     GEN b;
    2480          28 :     if (T0 == T)
    2481           0 :       b = pol_x(varn(T)); /* no improvement */
    2482             :     else
    2483          28 :       b = QXQ_reverse(a, T0); /* charpoly(Mod(b,T)) = S.x */
    2484          28 :     b = (degpol(T) == 1)? gmodulo(b, T): mkpolmod(b,T);
    2485          28 :     T = mkvec2(T, b);
    2486             :   }
    2487         616 :   return gerepilecopy(av, T);
    2488             : }
    2489             : /* DEPRECATED: backward compatibility */
    2490             : GEN
    2491          70 : polred0(GEN x, long flag, GEN fa)
    2492             : {
    2493          70 :   long fl = 0;
    2494          70 :   if (flag & 1) fl |= nf_PARTIALFACT;
    2495          70 :   if (flag & 2) fl |= nf_ORIG;
    2496          70 :   return Polred(x, fl, fa);
    2497             : }
    2498             : 
    2499             : GEN
    2500          21 : polredord(GEN x)
    2501             : {
    2502          21 :   pari_sp av = avma;
    2503             :   GEN v, lt;
    2504             :   long i, n, vx;
    2505             : 
    2506          21 :   if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
    2507          21 :   x = Q_primpart(x); RgX_check_ZX(x,"polredord");
    2508          21 :   n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
    2509          21 :   if (n == 1) return gerepilecopy(av, mkvec(x));
    2510          14 :   lt = leading_coeff(x); vx = varn(x);
    2511          14 :   if (is_pm1(lt))
    2512             :   {
    2513           7 :     if (signe(lt) < 0) x = ZX_neg(x);
    2514           7 :     v = pol_x_powers(n, vx);
    2515             :   }
    2516             :   else
    2517             :   { GEN L;
    2518             :     /* basis for Dedekind order */
    2519           7 :     v = cgetg(n+1, t_VEC);
    2520           7 :     gel(v,1) = scalarpol_shallow(lt, vx);
    2521          14 :     for (i = 2; i <= n; i++)
    2522           7 :       gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
    2523           7 :     gel(v,1) = pol_1(vx);
    2524           7 :     x = ZX_Q_normalize(x, &L);
    2525           7 :     v = gsubst(v, vx, monomial(ginv(L),1,vx));
    2526          14 :     for (i=2; i <= n; i++)
    2527           7 :       if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = monomial(gen_1, i-1, vx);
    2528             :   }
    2529          14 :   return gerepileupto(av, polred(mkvec2(x, v)));
    2530             : }
    2531             : 
    2532             : GEN
    2533          14 : polred(GEN x) { return Polred(x, 0, NULL); }
    2534             : GEN
    2535           0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
    2536             : GEN
    2537           0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
    2538             : GEN
    2539           0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
    2540             : GEN
    2541           0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
    2542             : GEN
    2543           0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
    2544             : 
    2545             : /********************************************************************/
    2546             : /**                                                                **/
    2547             : /**                           POLREDABS                            **/
    2548             : /**                                                                **/
    2549             : /********************************************************************/
    2550             : /* set V[k] := matrix of multiplication by nk.zk[k] */
    2551             : static GEN
    2552        1302 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
    2553             : {
    2554        1302 :   GEN v, Mk = cgetg(N+1, t_MAT);
    2555             :   long i, e;
    2556        1302 :   for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
    2557       12320 :   for (     ; i <=N; i++)
    2558             :   {
    2559       11018 :     v = vecmul(gel(M,k), gel(M,i));
    2560       11018 :     v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
    2561       11018 :     gel(Mk,i) = grndtoi(v, &e);
    2562       11018 :     if (e > -5) return NULL;
    2563             :   }
    2564        1302 :   gel(V,k) = Mk; return Mk;
    2565             : }
    2566             : 
    2567             : static GEN
    2568        1427 : ZM_image_shallow(GEN M, long *pr)
    2569             : {
    2570             :   long j, k, r;
    2571        1427 :   GEN y, d = ZM_pivots(M, &k);
    2572        1427 :   r = lg(M)-1 - k;
    2573        1427 :   y = cgetg(r+1,t_MAT);
    2574       10550 :   for (j=k=1; j<=r; k++)
    2575        9123 :     if (d[k]) gel(y,j++) = gel(M,k);
    2576        1427 :   *pr = r; return y;
    2577             : }
    2578             : 
    2579             : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
    2580             :  * Q from algo 2.7.6] */
    2581             : static GEN
    2582         351 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
    2583             : {
    2584         351 :   CG_data *d = (CG_data*)chk->data;
    2585             :   GEN P, V, D, inv, bound, S, M;
    2586         351 :   long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
    2587         351 :   long i, j, prec, firstprim = 0, skipfirst = 0;
    2588             :   pari_sp av;
    2589             : 
    2590         351 :   d->u = U;
    2591         351 :   d->ZKembed = M = RgM_mul(d->M, U);
    2592             : 
    2593         351 :   av = avma; bound = d->bound;
    2594         351 :   D = cgetg(N+1, t_VECSMALL);
    2595        2831 :   for (i = 1; i <= N; i++)
    2596             :   {
    2597        2488 :     pari_sp av2 = avma;
    2598        2488 :     P = get_pol(d, gel(M,i));
    2599        2488 :     if (!P) pari_err_PREC("chk_gen_init");
    2600        2480 :     (void)ZX_gcd_all(P, ZX_deriv(P), &P);
    2601        2480 :     P = gerepilecopy(av2, P);
    2602        2480 :     D[i] = degpol(P);
    2603        2480 :     if (D[i] == N)
    2604             :     { /* primitive element */
    2605        1030 :       GEN B = embed_T2(gel(M,i), r1);
    2606        1030 :       if (!firstprim) firstprim = i; /* index of first primitive element */
    2607        1030 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2608        1030 :       if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
    2609             :     }
    2610             :     else
    2611             :     {
    2612        1450 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
    2613        1450 :       if (firstprim)
    2614             :       { /* cycle basis vectors so that primitive elements come last */
    2615         159 :         GEN u = d->u, e = M;
    2616         159 :         GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
    2617         159 :         long tS = D[i];
    2618         447 :         for (j = i; j > firstprim; j--)
    2619             :         {
    2620         288 :           u[j] = u[j-1];
    2621         288 :           e[j] = e[j-1];
    2622         288 :           R[j] = R[j-1];
    2623         288 :           D[j] = D[j-1];
    2624             :         }
    2625         159 :         gel(u,firstprim) = tu;
    2626         159 :         gel(e,firstprim) = te;
    2627         159 :         gel(R,firstprim) = tR;
    2628         159 :         D[firstprim] = tS; firstprim++;
    2629             :       }
    2630             :     }
    2631             :   }
    2632         343 :   if (!firstprim)
    2633             :   { /* try (a little) to find primitive elements to improve bound */
    2634          21 :     GEN x = cgetg(N+1, t_VECSMALL);
    2635          21 :     if (DEBUGLEVEL>1)
    2636           0 :       err_printf("chk_gen_init: difficult field, trying random elements\n");
    2637         231 :     for (i = 0; i < 10; i++)
    2638             :     {
    2639             :       GEN e, B;
    2640         210 :       for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
    2641         210 :       e = RgM_zc_mul(M, x);
    2642         210 :       B = embed_T2(e, r1);
    2643         210 :       if (gcmp(B,bound) >= 0) continue;
    2644          14 :       P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
    2645          14 :       if (!ZX_is_squarefree(P)) continue;
    2646          14 :       if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
    2647          14 :       bound = B ;
    2648             :     }
    2649             :   }
    2650             : 
    2651         343 :   if (firstprim != 1)
    2652             :   {
    2653         343 :     inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
    2654         343 :     V = gel(inv,1);
    2655         343 :     for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
    2656             :     /* V corresponds to 1_Z */
    2657         343 :     V = grndtoi(V, &j);
    2658         343 :     if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
    2659         343 :     S = mkmat(V); /* 1 */
    2660             : 
    2661         343 :     V = cgetg(N+1, t_VEC);
    2662        1456 :     for (i = 1; i <= N; i++,skipfirst++)
    2663             :     { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
    2664             :       GEN Mx, M2;
    2665        1456 :       long j, k, h, rkM, dP = D[i];
    2666             : 
    2667        1456 :       if (dP == N) break; /* primitive */
    2668        1302 :       Mx = set_mulid(V, M, inv, r1, r2, N, i);
    2669        1302 :       if (!Mx) break; /* prec. problem. Stop */
    2670        1302 :       if (dP == 1) continue;
    2671        1008 :       rkM = lg(S)-1;
    2672        1008 :       M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
    2673        1008 :       gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
    2674        1008 :       k = 2;
    2675        3428 :       for (h = 1; h < dP; h++)
    2676             :       {
    2677             :         long r; /* add to M2 the elts of S * nf.zk[i]  */
    2678        1427 :         for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
    2679        1427 :         setlg(M2, k); k = 1;
    2680        1427 :         S = ZM_image_shallow(shallowconcat(S,M2), &r);
    2681        2148 :         if (r == rkM) break;
    2682         895 :         if (r > rkM)
    2683             :         {
    2684         895 :           rkM = r;
    2685         895 :           if (rkM == N) break;
    2686             :         }
    2687             :       }
    2688        1008 :       if (rkM == N) break;
    2689             :       /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
    2690             :     }
    2691             :   }
    2692             :   /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
    2693         343 :   chk->skipfirst = skipfirst;
    2694         343 :   if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
    2695             : 
    2696             :   /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
    2697         343 :   bound = gerepileuptoleaf(av, bound);
    2698         343 :   prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
    2699         343 :   if (DEBUGLEVEL)
    2700           0 :     err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
    2701         343 :   if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
    2702         343 :   if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
    2703         343 :   return bound;
    2704             : }
    2705             : 
    2706             : /* z "small" minimal polynomial of Mod(a,x), deg z = deg x */
    2707             : static GEN
    2708        2478 : store(GEN x, GEN z, GEN a, nfbasic_t *T, long flag, GEN u)
    2709             : {
    2710             :   GEN y, b;
    2711             : 
    2712        2478 :   if (u) a = RgV_RgC_mul(T->bas, ZM_ZC_mul(u, a));
    2713        2478 :   if (flag & (nf_ORIG|nf_ADDZK))
    2714             :   {
    2715         245 :     b = QXQ_reverse(a, x);
    2716         245 :     if (!isint1(T->unscale)) b = gdiv(b, T->unscale); /* not RgX_Rg_div */
    2717             :   }
    2718             :   else
    2719        2233 :     b = NULL;
    2720             : 
    2721        2478 :   if (flag & nf_RAW)
    2722          28 :     y = mkvec2(z, a);
    2723        2450 :   else if (flag & nf_ORIG) /* store phi(b mod z). */
    2724         245 :     y = mkvec2(z, mkpolmod(b,z));
    2725             :   else
    2726        2205 :     y = z;
    2727        2478 :   if (flag & nf_ADDZK)
    2728             :   { /* append integral basis for number field Q[X]/(z) to result */
    2729           0 :     long n = degpol(x);
    2730           0 :     GEN t = RgV_RgM_mul(RgXQ_powers(b, n-1, z), RgV_to_RgM(T->bas,n));
    2731           0 :     y = mkvec2(y, t);
    2732             :   }
    2733        2478 :   return y;
    2734             : }
    2735             : static GEN
    2736         343 : polredabs_aux(nfbasic_t *T, GEN *u)
    2737             : {
    2738             :   long prec;
    2739             :   GEN v;
    2740         343 :   FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
    2741             :   nffp_t F;
    2742         343 :   CG_data d; chk.data = (void*)&d;
    2743             : 
    2744         343 :   prec = polred_init(T, &F, &d);
    2745         343 :   d.bound = embed_T2(F.ro, d.r1);
    2746         343 :   if (realprec(d.bound) > prec) d.bound = rtor(d.bound, prec);
    2747             :   for (;;)
    2748             :   {
    2749         363 :     GEN R = R_from_QR(F.G, prec);
    2750         363 :     if (R)
    2751             :     {
    2752         351 :       d.prec = prec;
    2753         351 :       d.M    = F.M;
    2754         351 :       v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
    2755         351 :       if (v) break;
    2756             :     }
    2757          20 :     prec = precdbl(prec);
    2758          20 :     get_nf_fp_compo(T, &F, NULL, 1, prec);
    2759          20 :     if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",prec);
    2760          20 :   }
    2761         343 :   *u = d.u; return v;
    2762             : }
    2763             : 
    2764             : GEN
    2765         357 : polredabs0(GEN x, long flag)
    2766             : {
    2767         357 :   pari_sp av = avma;
    2768             :   long i, l, vx;
    2769             :   GEN y, a, u;
    2770             :   nfbasic_t T;
    2771             : 
    2772         357 :   nfbasic_init(x, nf_PARTIALFACT, &T);
    2773         357 :   x = T.x; vx = varn(x);
    2774             : 
    2775         357 :   if (degpol(x) == 1)
    2776             :   {
    2777          14 :     u = NULL;
    2778          14 :     y = mkvec( pol_x(vx) );
    2779          14 :     a = mkvec( deg1pol_shallow(gen_1, negi(gel(x,2)), vx) );
    2780          14 :     l = 2;
    2781             :   }
    2782             :   else
    2783             :   {
    2784             :     GEN v;
    2785         343 :     if (!(flag & nf_PARTIALFACT) && T.dKP)
    2786             :     {
    2787          84 :       GEN vw = primes_certify(T.dK, T.dKP);
    2788          84 :       v = gel(vw,1); l = lg(v);
    2789          84 :       if (l != 1)
    2790             :       { /* fix integral basis */
    2791           7 :         GEN w = gel(vw,2);
    2792          14 :         for (i = 1; i < l; i++)
    2793           7 :           w = ZV_union_shallow(w, gel(Z_factor(gel(v,i)),1));
    2794           7 :         nfbasic_init(mkvec2(x,w), 0, &T);
    2795             :       }
    2796             :     }
    2797         343 :     v = polredabs_aux(&T, &u);
    2798         343 :     y = gel(v,1);
    2799         343 :     a = gel(v,2); l = lg(a);
    2800        6358 :     for (i=1; i<l; i++)
    2801        6015 :       if (ZX_canon_neg(gel(y,i))) gel(a,i) = ZC_neg(gel(a,i));
    2802         343 :     remove_duplicates(y,a);
    2803         343 :     l = lg(a);
    2804         343 :     if (l == 1)
    2805           0 :       pari_err_BUG("polredabs (missing vector)");
    2806             :   }
    2807         357 :   if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
    2808         357 :   if (flag & nf_ALL) {
    2809           7 :     for (i=1; i<l; i++) gel(y,i) = store(x, gel(y,i), gel(a,i), &T, flag, u);
    2810             :   } else {
    2811         350 :     GEN z = findmindisc(y, &a);
    2812         350 :     y = store(x, z, a, &T, flag, u);
    2813             :   }
    2814         357 :   return gerepilecopy(av, y);
    2815             : }
    2816             : 
    2817             : GEN
    2818           0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
    2819             : GEN
    2820           0 : polredabs(GEN x) { return polredabs0(x,0); }
    2821             : GEN
    2822           0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
    2823             : 
    2824             : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
    2825             :  * flag & nf_ORIG: + element (base change)
    2826             :  * flag & nf_ABSOLUTE: absolute polynomial */
    2827             : static GEN
    2828          63 : rnfpolred_i(GEN nf, GEN relpol, long flag, long best)
    2829             : {
    2830          63 :   const char *f = best? "rnfpolredbest": "rnfpolredabs";
    2831          63 :   const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
    2832             :   pari_timer ti;
    2833          63 :   GEN listP = NULL, red, bas, A, P, pol, T, rnfeq;
    2834          63 :   long ty = typ(relpol);
    2835          63 :   pari_sp av = avma;
    2836             : 
    2837          63 :   if (ty == t_VEC) {
    2838          14 :     if (lg(relpol) != 3) pari_err_TYPE(f,relpol);
    2839          14 :     listP = gel(relpol,2);
    2840          14 :     relpol = gel(relpol,1);
    2841             :   }
    2842          63 :   if (typ(relpol) != t_POL) pari_err_TYPE(f,relpol);
    2843          63 :   nf = checknf(nf);
    2844          63 :   if (DEBUGLEVEL>1) timer_start(&ti);
    2845          63 :   T = nf_get_pol(nf);
    2846          63 :   relpol = RgX_nffix(f, T, relpol, 0);
    2847          63 :   if (best || (flag & nf_PARTIALFACT))
    2848             :   {
    2849          42 :     if (abs)
    2850             :     {
    2851           7 :       rnfeq = nf_rnfeq(nf, relpol);
    2852           7 :       pol = gel(rnfeq,1);
    2853             :     }
    2854             :     else
    2855             :     {
    2856             :       long sa;
    2857          35 :       pol = rnfequationall(nf, relpol, &sa, NULL);
    2858          35 :       rnfeq = mkvec5(gen_0,gen_0,stoi(sa),T,liftpol_shallow(relpol));
    2859             :     }
    2860          42 :     bas = listP? mkvec2(pol, listP): pol;
    2861          84 :     if (best)
    2862             :     {
    2863          35 :       if (abs) red = polredbest(bas, 1);
    2864             :       else
    2865             :       {
    2866             :         GEN ro, x, dx, a;
    2867             :         nfbasic_t T;
    2868          28 :         nfbasic_init(bas, nf_PARTIALFACT, &T);
    2869          28 :         polredbest_aux(&T, &ro, &x, &dx, &a);
    2870          28 :         red = mkvec2(x, a);
    2871             :       }
    2872             :     }
    2873             :     else
    2874           7 :       red = polredabs0(bas, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
    2875             :   }
    2876             :   else
    2877             :   {
    2878          21 :     GEN rnf = rnfinit(nf, relpol), M = rnf_basM(rnf);
    2879          21 :     rnfeq = rnf_get_map(rnf);
    2880          21 :     pol = rnf_get_polabs(rnf);
    2881          21 :     bas = mkvec2(pol, RgM_to_RgXV(M, varn(pol)));
    2882          21 :     if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
    2883          21 :     red = polredabs0(bas, nf_RAW);
    2884             :   }
    2885          63 :   P = gel(red,1);
    2886          63 :   A = gel(red,2);
    2887          63 :   if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
    2888          63 :   if (flag & nf_ABSOLUTE)
    2889             :   {
    2890          14 :     if (flag & nf_ORIG)
    2891             :     {
    2892           7 :       GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
    2893           7 :       GEN k = gel(rnfeq,3); /* Mod(variable(relpol),relpol) + k*a root of pol */
    2894           7 :       a = RgX_RgXQ_eval(a, lift_intern(A), P); /* Mod(a, P) root of T */
    2895           7 :       P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
    2896             :     }
    2897          14 :     return gerepilecopy(av, P);
    2898             :   }
    2899          49 :   A = eltabstorel_lift(rnfeq, A);
    2900          49 :   P = RgXQ_charpoly(A, relpol, varn(relpol));
    2901          49 :   P = lift_if_rational(P);
    2902          49 :   if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,relpol),P));
    2903          49 :   return gerepilecopy(av, P);
    2904             : }
    2905             : GEN
    2906          28 : rnfpolredabs(GEN nf, GEN relpol, long flag)
    2907          28 : { return rnfpolred_i(nf,relpol,flag, 0); }
    2908             : GEN
    2909          35 : rnfpolredbest(GEN nf, GEN relpol, long flag)
    2910             : {
    2911          35 :   if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
    2912          35 :   return rnfpolred_i(nf,relpol,flag, 1);
    2913             : }

Generated by: LCOV version 1.11