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

Generated by: LCOV version 1.13