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 22303-eb3e11d) Lines: 1534 1614 95.0 %
Date: 2018-04-20 06:16:30 Functions: 113 129 87.6 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11