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

Generated by: LCOV version 1.11