Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19623-dc26710) Lines: 1344 1447 92.9 %
Date: 2016-09-30 05:54:20 Functions: 107 122 87.7 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11