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-bordeaux1.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 - polarit3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17097-9391e68) Lines: 1453 1684 86.3 %
Date: 2014-11-21 Functions: 130 147 88.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 817 1183 69.1 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000-2005  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                 :            : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17                 :            : /**                         (third part)                              **/
      18                 :            : /**                                                                   **/
      19                 :            : /***********************************************************************/
      20                 :            : #include "pari.h"
      21                 :            : #include "paripriv.h"
      22                 :            : 
      23                 :            : /************************************************************************
      24                 :            :  **                                                                    **
      25                 :            :  **                      Ring membership                               **
      26                 :            :  **                                                                    **
      27                 :            :  ************************************************************************/
      28                 :            : struct charact {
      29                 :            :   GEN q;
      30                 :            :   int isprime;
      31                 :            : };
      32                 :            : static void
      33                 :        465 : char_update_prime(struct charact *S, GEN p)
      34                 :            : {
      35         [ +  + ]:        465 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36         [ -  + ]:        465 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37                 :        465 : }
      38                 :            : static void
      39                 :        415 : char_update_int(struct charact *S, GEN n)
      40                 :            : {
      41         [ -  + ]:        415 :   if (S->isprime)
      42                 :            :   {
      43         [ #  # ]:        415 :     if (dvdii(n, S->q)) return;
      44                 :          0 :     pari_err_MODULUS("characteristic", S->q, n);
      45                 :            :   }
      46                 :        415 :   S->q = gcdii(S->q, n);
      47                 :            : }
      48                 :            : static void
      49                 :     384870 : charact(struct charact *S, GEN x)
      50                 :            : {
      51                 :     384870 :   const long tx = typ(x);
      52                 :            :   long i, l;
      53   [ +  +  +  +  :     384870 :   switch(tx)
                      + ]
      54                 :            :   {
      55                 :        225 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56                 :        420 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      57                 :            :     case t_COMPLEX: case t_QUAD:
      58                 :            :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      59                 :            :     case t_VEC: case t_COL: case t_MAT:
      60                 :       1175 :       l = lg(x);
      61         [ +  + ]:     384380 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62                 :       1175 :       break;
      63                 :            :     case t_LIST:
      64                 :          5 :       x = list_data(x);
      65         [ -  + ]:          5 :       if (x) charact(S, x);
      66                 :          5 :       break;
      67                 :            :   }
      68                 :     384870 : }
      69                 :            : static void
      70                 :      21575 : charact_res(struct charact *S, GEN x)
      71                 :            : {
      72                 :      21575 :   const long tx = typ(x);
      73                 :            :   long i, l;
      74   [ +  -  +  +  :      21575 :   switch(tx)
                   -  + ]
      75                 :            :   {
      76                 :        190 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      77                 :          0 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      78                 :         45 :     case t_PADIC: char_update_prime(S, gel(x,2)); break;
      79                 :            :     case t_COMPLEX: case t_QUAD:
      80                 :            :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      81                 :            :     case t_VEC: case t_COL: case t_MAT:
      82                 :       6635 :       l = lg(x);
      83         [ +  + ]:      26555 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84                 :       6635 :       break;
      85                 :            :     case t_LIST:
      86                 :          0 :       x = list_data(x);
      87         [ #  # ]:          0 :       if (x) charact_res(S, x);
      88                 :          0 :       break;
      89                 :            :   }
      90                 :      21575 : }
      91                 :            : GEN
      92                 :       1665 : characteristic(GEN x)
      93                 :            : {
      94                 :            :   struct charact S;
      95                 :       1665 :   S.q = gen_0; S.isprime = 0;
      96                 :       1665 :   charact(&S, x); return S.q;
      97                 :            : }
      98                 :            : GEN
      99                 :       1655 : residual_characteristic(GEN x)
     100                 :            : {
     101                 :            :   struct charact S;
     102                 :       1655 :   S.q = gen_0; S.isprime = 0;
     103                 :       1655 :   charact_res(&S, x); return S.q;
     104                 :            : }
     105                 :            : 
     106                 :            : int
     107                 :   12821676 : Rg_is_Fp(GEN x, GEN *pp)
     108                 :            : {
     109                 :            :   GEN mod;
     110      [ +  +  + ]:   12821676 :   switch(typ(x))
     111                 :            :   {
     112                 :            :   case t_INTMOD:
     113                 :    2724695 :     mod = gel(x,1);
     114         [ +  + ]:    2724695 :     if (!*pp) *pp = mod;
     115 [ +  + ][ -  + ]:    2689550 :     else if (mod != *pp && !equalii(mod, *pp))
     116                 :            :     {
     117         [ #  # ]:          0 :       if (DEBUGMEM) pari_warn(warner,"different moduli in Rg_is_Fp");
     118                 :          0 :       return 0;
     119                 :            :     }
     120                 :    2724695 :     return 1;
     121                 :            :   case t_INT:
     122                 :    9625934 :     return 1;
     123                 :   12821676 :   default: return 0;
     124                 :            :   }
     125                 :            : }
     126                 :            : 
     127                 :            : int
     128                 :     400111 : RgX_is_FpX(GEN x, GEN *pp)
     129                 :            : {
     130                 :     400111 :   long i, lx = lg(x);
     131         [ +  + ]:    4418104 :   for (i=2; i<lx; i++)
     132         [ +  + ]:    4122520 :     if (!Rg_is_Fp(gel(x, i), pp))
     133                 :     104527 :       return 0;
     134                 :     400111 :   return 1;
     135                 :            : }
     136                 :            : 
     137                 :            : int
     138                 :    1528996 : RgV_is_FpV(GEN x, GEN *pp)
     139                 :            : {
     140                 :    1528996 :   long i, lx = lg(x);
     141         [ +  + ]:    9861127 :   for (i=1; i<lx; i++)
     142         [ +  + ]:    8698646 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     143                 :    1528996 :   return 1;
     144                 :            : }
     145                 :            : 
     146                 :            : int
     147                 :     507135 : RgM_is_FpM(GEN x, GEN *pp)
     148                 :            : {
     149                 :     507135 :   long i, lx = lg(x);
     150         [ +  + ]:    1665616 :   for (i=1; i<lx; i++)
     151         [ +  + ]:    1524996 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     152                 :     507135 :   return 1;
     153                 :            : }
     154                 :            : 
     155                 :            : int
     156                 :      37985 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157                 :            : {
     158                 :            :   GEN pol, mod, p;
     159   [ +  +  +  +  :      37985 :   switch(typ(x))
                   +  + ]
     160                 :            :   {
     161                 :            :   case t_INTMOD:
     162                 :        500 :     return Rg_is_Fp(x, pp);
     163                 :            :   case t_INT:
     164                 :      11660 :     return 1;
     165                 :            :   case t_POL:
     166                 :       5870 :     return RgX_is_FpX(x, pp);
     167                 :            :   case t_FFELT:
     168                 :       3160 :     mod = FF_1(x); p = FF_p_i(x);
     169         [ +  + ]:       3160 :     if (!*pp) *pp = p;
     170         [ +  + ]:       3160 :     if (!*pT) *pT = mod;
     171 [ +  + ][ +  - ]:       3160 :     if ((p != *pp && !equalii(p, *pp)) || (mod != *pT && !gequal(mod, *pT)))
         [ +  + ][ -  + ]
     172                 :            :     {
     173         [ #  # ]:          0 :       if (DEBUGMEM) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     174                 :          0 :       return 0;
     175                 :            :     }
     176                 :       3160 :     return 1;
     177                 :            :   case t_POLMOD:
     178                 :       1140 :     mod = gel(x,1); pol = gel(x, 2);
     179         [ -  + ]:       1140 :     if (!RgX_is_FpX(mod, pp)) return 0;
     180         [ +  + ]:       1140 :     if (typ(pol)==t_POL)
     181                 :            :     {
     182         [ +  + ]:       1130 :       if (!RgX_is_FpX(pol, pp)) return 0;
     183                 :            :     }
     184         [ +  + ]:         10 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     185         [ +  + ]:       1125 :     if (!*pT) *pT = mod;
     186 [ +  - ][ -  + ]:        910 :     else if (mod != *pT && !gequal(mod, *pT))
     187                 :            :     {
     188         [ #  # ]:          0 :       if (DEBUGMEM) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     189                 :          0 :       return 0;
     190                 :            :     }
     191                 :       1125 :     return 1;
     192                 :            : 
     193                 :      37985 :   default: return 0;
     194                 :            :   }
     195                 :            : }
     196                 :            : 
     197                 :            : int
     198                 :      19270 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199                 :            : {
     200                 :      19270 :   long i, lx = lg(x);
     201         [ +  + ]:      40630 :   for (i = 2; i < lx; i++)
     202         [ +  + ]:      37810 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203                 :      19270 :   return 1;
     204                 :            : }
     205                 :            : 
     206                 :            : /************************************************************************
     207                 :            :  **                                                                    **
     208                 :            :  **                      Ring conversion                               **
     209                 :            :  **                                                                    **
     210                 :            :  ************************************************************************/
     211                 :            : 
     212                 :            : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
     213                 :            :  * If x is an INTMOD, assume modulus is a multiple of p. */
     214                 :            : GEN
     215                 :    8748504 : Rg_to_Fp(GEN x, GEN p)
     216                 :            : {
     217         [ +  + ]:    8748504 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218   [ +  +  -  +  :     263540 :   switch(typ(x))
                      - ]
     219                 :            :   {
     220                 :      10996 :     case t_INT: return modii(x, p);
     221                 :            :     case t_FRAC: {
     222                 :         52 :       pari_sp av = avma;
     223                 :         52 :       GEN z = modii(gel(x,1), p);
     224         [ -  + ]:         52 :       if (z == gen_0) return gen_0;
     225                 :         52 :       return gerepileuptoint(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
     226                 :            :     }
     227                 :          0 :     case t_PADIC: return padic_to_Fp(x, p);
     228                 :            :     case t_INTMOD: {
     229                 :     252492 :       GEN q = gel(x,1), a = gel(x,2);
     230         [ +  - ]:     252492 :       if (equalii(q, p)) return icopy(a);
     231         [ #  # ]:          0 :       if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
     232                 :          0 :       return remii(a, p);
     233                 :            :     }
     234                 :          0 :     default: pari_err_TYPE("Rg_to_Fp",x);
     235                 :    8748424 :       return NULL; /* not reached */
     236                 :            :   }
     237                 :            : }
     238                 :            : /* If x is a POLMOD, assume modulus is a multiple of T. */
     239                 :            : GEN
     240                 :      17156 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241                 :            : {
     242                 :      17156 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243                 :            :   GEN a, b;
     244         [ +  + ]:      17156 :   if (is_const_t(tx))
     245                 :            :   {
     246         [ +  + ]:      17031 :     if (tx == t_FFELT)
     247                 :            :     {
     248                 :       6759 :       GEN z = FF_to_FpXQ(x);
     249                 :       6759 :       setvarn(z, v);
     250                 :       6759 :       return z;
     251                 :            :     }
     252                 :      10272 :     return scalar_ZX(Rg_to_Fp(x, p), v);
     253                 :            :   }
     254   [ +  +  -  - ]:        125 :   switch(tx)
     255                 :            :   {
     256                 :            :     case t_POLMOD:
     257                 :         35 :       b = gel(x,1);
     258                 :         35 :       a = gel(x,2); ta = typ(a);
     259         [ -  + ]:         35 :       if (is_const_t(ta)) return scalar_ZX(Rg_to_Fp(a, p), v);
     260         [ -  + ]:         35 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     261         [ +  - ]:         35 :       a = RgX_to_FpX(a, p); if (ZX_equal(b,get_FpX_mod(T))) return a;
     262                 :          0 :       return FpX_rem(a, T, p);
     263                 :            :     case t_POL:
     264         [ -  + ]:         90 :       if (varn(x) != v) break;
     265                 :         90 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     266                 :            :     case t_RFRAC:
     267                 :          0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     268                 :          0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     269                 :          0 :       return FpXQ_div(a,b, T,p);
     270                 :            :   }
     271                 :          0 :   pari_err_TYPE("Rg_to_FpXQ",x);
     272                 :      17156 :   return NULL; /* not reached */
     273                 :            : }
     274                 :            : GEN
     275                 :      92670 : RgX_to_FpX(GEN x, GEN p)
     276                 :            : {
     277                 :            :   long i, l;
     278                 :      92670 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     279         [ +  + ]:    2715535 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     280                 :      92670 :   return FpX_renormalize(z, l);
     281                 :            : }
     282                 :            : 
     283                 :            : GEN
     284                 :        765 : RgV_to_FpV(GEN x, GEN p)
     285                 :            : {
     286                 :        765 :   long i, l = lg(x);
     287                 :        765 :   GEN z = cgetg(l, t_VEC);
     288         [ +  + ]:       2300 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     289                 :        765 :   return z;
     290                 :            : }
     291                 :            : 
     292                 :            : GEN
     293                 :       1556 : RgC_to_FpC(GEN x, GEN p)
     294                 :            : {
     295                 :       1556 :   long i, l = lg(x);
     296                 :       1556 :   GEN z = cgetg(l, t_COL);
     297         [ +  + ]:      55173 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     298                 :       1556 :   return z;
     299                 :            : }
     300                 :            : 
     301                 :            : GEN
     302                 :        203 : RgM_to_FpM(GEN x, GEN p)
     303                 :            : {
     304                 :        203 :   long i, l = lg(x);
     305                 :        203 :   GEN z = cgetg(l, t_MAT);
     306         [ +  + ]:       1744 :   for (i = 1; i < l; i++) gel(z,i) = RgC_to_FpC(gel(x,i), p);
     307                 :        203 :   return z;
     308                 :            : }
     309                 :            : GEN
     310                 :       8964 : RgC_to_Flc(GEN x, ulong p)
     311                 :            : {
     312                 :       8964 :   long l = lg(x), i;
     313                 :       8964 :   GEN a = cgetg(l, t_VECSMALL);
     314         [ +  + ]:      75182 :   for (i=1; i<l; i++) a[i] = Rg_to_Fl(gel(x,i), p);
     315                 :       8964 :   return a;
     316                 :            : }
     317                 :            : GEN
     318                 :       1242 : RgM_to_Flm(GEN x, ulong p)
     319                 :            : {
     320                 :            :   long l, i;
     321                 :       1242 :   GEN a = cgetg_copy(x, &l);
     322         [ +  + ]:       9886 :   for (i=1; i<l; i++) gel(a,i) = RgC_to_Flc(gel(x,i), p);
     323                 :       1242 :   return a;
     324                 :            : }
     325                 :            : 
     326                 :            : GEN
     327                 :         15 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     328                 :            : {
     329                 :         15 :   long i, l = lg(x);
     330                 :         15 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     331         [ +  + ]:        105 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     332                 :         15 :   return FpXQX_renormalize(z, l);
     333                 :            : }
     334                 :            : GEN
     335                 :        345 : RgX_to_FqX(GEN x, GEN T, GEN p)
     336                 :            : {
     337                 :        345 :   long i, l = lg(x);
     338                 :        345 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     339         [ +  - ]:        345 :   if (T)
     340         [ +  + ]:       4975 :     for (i = 2; i < l; i++)
     341                 :       4630 :       gel(z,i) = simplify_shallow(Rg_to_FpXQ(gel(x,i), T, p));
     342                 :            :   else
     343         [ #  # ]:          0 :     for (i = 2; i < l; i++)
     344                 :          0 :       gel(z,i) = Rg_to_Fp(gel(x,i), p);
     345                 :        345 :   return FpXQX_renormalize(z, l);
     346                 :            : }
     347                 :            : 
     348                 :            : /* lg(V) > 1 */
     349                 :            : GEN
     350                 :     303890 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     351                 :            : {
     352                 :     303890 :   pari_sp av = avma;
     353                 :     303890 :   long i, l = lg(V);
     354                 :     303890 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     355         [ +  + ]:    1505340 :   for(i=2; i<l; i++)
     356                 :            :   {
     357                 :    1201450 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     358         [ +  + ]:    1201450 :     if ((i & 7) == 0) z = gerepileupto(av, z);
     359                 :            :   }
     360                 :     303890 :   return gerepileupto(av, FpX_red(z,p));
     361                 :            : }
     362                 :            : 
     363                 :            : GEN
     364                 :       4742 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     365                 :            : {
     366                 :       4742 :   long i, lz = lg(y);
     367                 :            :   GEN z;
     368         [ +  + ]:       4742 :   if (!T) return FpX_Fp_add(y, x, p);
     369         [ +  + ]:        877 :   if (lz == 2) return scalarpol(x, varn(y));
     370                 :        871 :   z = cgetg(lz,t_POL); z[1] = y[1];
     371                 :        871 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     372         [ +  + ]:        871 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     373                 :            :   else
     374         [ +  + ]:       1221 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     375                 :       4742 :   return z;
     376                 :            : }
     377                 :            : 
     378                 :            : GEN
     379                 :      11435 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     380                 :            : {
     381                 :            :   long i, lP;
     382                 :      11435 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     383         [ +  + ]:      33075 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     384                 :      11435 :   gel(res,lP-1) = gen_1; return res;
     385                 :            : }
     386                 :            : 
     387                 :            : GEN
     388                 :      30480 : FqX_normalize(GEN z, GEN T, GEN p)
     389                 :            : {
     390                 :            :   GEN lc;
     391         [ +  + ]:      30480 :   if (!T) return FpX_normalize(z,p);
     392         [ -  + ]:      24985 :   if (lg(z) == 2) return z;
     393                 :      24985 :   lc = leading_term(z);
     394         [ +  + ]:      24985 :   if (typ(lc) == t_POL)
     395                 :            :   {
     396         [ +  + ]:      11225 :     if (lg(lc) > 3) /* non-constant */
     397                 :      10775 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     398                 :            :     /* constant */
     399                 :        450 :     lc = gel(lc,2);
     400                 :        450 :     z = shallowcopy(z);
     401                 :        450 :     gel(z, lg(z)-1) = lc;
     402                 :            :   }
     403                 :            :   /* lc a t_INT */
     404         [ +  + ]:      14210 :   if (equali1(lc)) return z;
     405                 :      30480 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     406                 :            : }
     407                 :            : 
     408                 :            : GEN
     409                 :      73665 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     410                 :            : {
     411                 :            :   pari_sp av;
     412                 :            :   GEN p1, r;
     413                 :      73665 :   long j, i=lg(x)-1;
     414         [ +  + ]:      73665 :   if (i<=2)
     415         [ +  + ]:      12227 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     416                 :      61438 :   av=avma; p1=gel(x,i);
     417                 :            :   /* specific attention to sparse polynomials (see poleval)*/
     418                 :            :   /*You've guessed it! It's a copy-paste(tm)*/
     419         [ +  + ]:     186796 :   for (i--; i>=2; i=j-1)
     420                 :            :   {
     421         [ +  + ]:     125543 :     for (j=i; !signe(gel(x,j)); j--)
     422         [ +  + ]:        185 :       if (j==2)
     423                 :            :       {
     424         [ +  + ]:        125 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     425                 :        125 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     426                 :            :       }
     427         [ +  + ]:     125358 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     428                 :     125358 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     429                 :            :   }
     430                 :      73665 :   return gerepileupto(av, p1);
     431                 :            : }
     432                 :            : 
     433                 :            : GEN
     434                 :      19305 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     435                 :            : {
     436                 :      19305 :   long i, lb = lg(Q);
     437                 :            :   GEN z;
     438         [ +  + ]:      19305 :   if (!T) return FpXY_evalx(Q, x, p);
     439                 :      14490 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     440         [ +  + ]:      79685 :   for (i=2; i<lb; i++)
     441                 :            :   {
     442                 :      65195 :     GEN q = gel(Q,i);
     443         [ +  + ]:      65195 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     444                 :            :   }
     445                 :      19305 :   return FpXQX_renormalize(z, lb);
     446                 :            : }
     447                 :            : 
     448                 :            : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     449                 :            : GEN
     450                 :       6560 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     451                 :            : {
     452                 :       6560 :   pari_sp av = avma;
     453         [ +  + ]:       6560 :   if (!T) return FpXY_eval(Q, y, x, p);
     454                 :       6560 :   return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     455                 :            : }
     456                 :            : 
     457                 :            : /* a X^d */
     458                 :            : GEN
     459                 :    1582395 : monomial(GEN a, long d, long v)
     460                 :            : {
     461                 :    1582395 :   long i, lP = d+3;
     462                 :            :   GEN P;
     463         [ +  + ]:    1582395 :   if (d < 0) {
     464         [ -  + ]:         40 :     if (isrationalzero(a)) return pol_0(v);
     465                 :         40 :     P = cgetg(3, t_RFRAC);
     466                 :         40 :     gel(P,1) = a;
     467                 :         40 :     gel(P,2) = monomial(gen_1, -d, v);
     468                 :            :   } else {
     469                 :    1582355 :     P = cgetg(lP, t_POL);
     470         [ +  + ]:    1582355 :     if (gequal0(a))
     471                 :            :     {
     472         [ +  - ]:        190 :       if (isexactzero(a)) return scalarpol_shallow(a,v);
     473                 :          0 :       P[1] = evalsigne(0) | evalvarn(v);
     474                 :            :     }
     475                 :            :     else
     476                 :    1582165 :       P[1] = evalsigne(1) | evalvarn(v);
     477                 :    1582165 :     lP--; gel(P,lP) = a;
     478         [ +  + ]:    3048298 :     for (i=2; i<lP; i++) gel(P,i) = gen_0;
     479                 :            :   }
     480                 :    1582395 :   return P;
     481                 :            : }
     482                 :            : GEN
     483                 :    5295739 : monomialcopy(GEN a, long d, long v)
     484                 :            : {
     485                 :    5295739 :   long i, lP = d+3;
     486                 :            :   GEN P;
     487         [ +  + ]:    5295739 :   if (d < 0) {
     488         [ -  + ]:          5 :     if (isrationalzero(a)) return pol_0(v);
     489                 :          5 :     P = cgetg(3, t_RFRAC);
     490                 :          5 :     gel(P,1) = gcopy(a);
     491                 :          5 :     gel(P,2) = monomial(gen_1, -d, v);
     492                 :            :   } else {
     493                 :    5295734 :     P = cgetg(lP, t_POL);
     494         [ +  + ]:    5295734 :     if (gequal0(a))
     495                 :            :     {
     496         [ +  - ]:          5 :       if (isexactzero(a)) return scalarpol(a,v);
     497                 :          0 :       P[1] = evalsigne(0) | evalvarn(v);
     498                 :            :     }
     499                 :            :     else
     500                 :    5295729 :       P[1] = evalsigne(1) | evalvarn(v);
     501                 :    5295729 :     lP--; gel(P,lP) = gcopy(a);
     502         [ +  + ]:    9737506 :     for (i=2; i<lP; i++) gel(P,i) = gen_0;
     503                 :            :   }
     504                 :    5295739 :   return P;
     505                 :            : }
     506                 :            : GEN
     507                 :      13750 : pol_x_powers(long N, long v)
     508                 :            : {
     509                 :      13750 :   GEN L = cgetg(N+1,t_VEC);
     510                 :            :   long i;
     511         [ +  + ]:      46670 :   for (i=1; i<=N; i++) gel(L,i) = monomial(gen_1, i-1, v);
     512                 :      13750 :   return L;
     513                 :            : }
     514                 :            : 
     515                 :            : GEN
     516                 :          0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     517                 :            : {
     518         [ #  # ]:          0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     519                 :            : }
     520                 :            : 
     521                 :            : GEN
     522                 :          0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     523                 :            : {
     524         [ #  # ]:          0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     525                 :            : }
     526                 :            : 
     527                 :            : /*******************************************************************/
     528                 :            : /*                                                                 */
     529                 :            : /*                             Fq                                  */
     530                 :            : /*                                                                 */
     531                 :            : /*******************************************************************/
     532                 :            : 
     533                 :            : GEN
     534                 :    3050983 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     535                 :            : {
     536                 :            :   (void)T;
     537 [ +  + ][ +  +  :    3050983 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
                +  +  - ]
     538                 :            :   {
     539                 :    1362548 :     case 0: return Fp_add(x,y,p);
     540                 :     131584 :     case 1: return FpX_Fp_add(x,y,p);
     541                 :     209604 :     case 2: return FpX_Fp_add(y,x,p);
     542                 :    1347247 :     case 3: return FpX_add(x,y,p);
     543                 :            :   }
     544                 :    3050983 :   return NULL;
     545                 :            : }
     546                 :            : 
     547                 :            : GEN
     548                 :    8228551 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     549                 :            : {
     550                 :            :   (void)T;
     551 [ +  + ][ +  +  :    8228551 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
                +  +  - ]
     552                 :            :   {
     553                 :      70857 :     case 0: return Fp_sub(x,y,p);
     554                 :      18791 :     case 1: return FpX_Fp_sub(x,y,p);
     555                 :      11565 :     case 2: return Fp_FpX_sub(x,y,p);
     556                 :    8127338 :     case 3: return FpX_sub(x,y,p);
     557                 :            :   }
     558                 :    8228551 :   return NULL;
     559                 :            : }
     560                 :            : 
     561                 :            : GEN
     562                 :     158180 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     563                 :            : {
     564                 :            :   (void)T;
     565         [ +  + ]:     158180 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     566                 :            : }
     567                 :            : 
     568                 :            : GEN
     569                 :       5910 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     570                 :            : {
     571                 :            :   (void)T;
     572         [ +  + ]:       5910 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     573                 :            : }
     574                 :            : 
     575                 :            : /* If T==NULL do not reduce*/
     576                 :            : GEN
     577                 :   35733458 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     578                 :            : {
     579 [ +  + ][ +  +  :   35733458 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
                +  +  - ]
     580                 :            :   {
     581                 :    1431779 :     case 0: return Fp_mul(x,y,p);
     582                 :     495118 :     case 1: return FpX_Fp_mul(x,y,p);
     583                 :     208609 :     case 2: return FpX_Fp_mul(y,x,p);
     584         [ +  + ]:   33597952 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     585                 :    7651383 :             else return FpX_mul(x,y,p);
     586                 :            :   }
     587                 :   35733458 :   return NULL;
     588                 :            : }
     589                 :            : 
     590                 :            : /* If T==NULL do not reduce*/
     591                 :            : GEN
     592                 :     337781 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     593                 :            : {
     594                 :            :   (void) T;
     595         [ +  + ]:     337781 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     596                 :            : }
     597                 :            : 
     598                 :            : /* y t_INT */
     599                 :            : GEN
     600                 :      16250 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     601                 :            : {
     602                 :            :   (void)T;
     603                 :      16250 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     604         [ +  + ]:      16250 :                           : Fp_mul(x,y,p);
     605                 :            : }
     606                 :            : /* If T==NULL do not reduce*/
     607                 :            : GEN
     608                 :      40444 : Fq_sqr(GEN x, GEN T, GEN p)
     609                 :            : {
     610         [ +  + ]:      40444 :   if (typ(x) == t_POL)
     611                 :            :   {
     612         [ +  - ]:       7755 :     if (T) return FpXQ_sqr(x,T,p);
     613                 :          0 :     else return FpX_sqr(x,p);
     614                 :            :   }
     615                 :            :   else
     616                 :      40444 :     return Fp_sqr(x,p);
     617                 :            :   return NULL;
     618                 :            : }
     619                 :            : 
     620                 :            : GEN
     621                 :          0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     622                 :            : {
     623         [ #  # ]:          0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     624                 :          0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     625                 :            : }
     626                 :            : 
     627                 :            : GEN
     628                 :          0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     629                 :            : {
     630         [ #  # ]:          0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     631                 :          0 :   return FpXQ_invsafe(x,pol,p);
     632                 :            : }
     633                 :            : 
     634                 :            : GEN
     635                 :      25177 : Fq_inv(GEN x, GEN pol, GEN p)
     636                 :            : {
     637         [ +  + ]:      25177 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     638                 :      25177 :   return FpXQ_inv(x,pol,p);
     639                 :            : }
     640                 :            : 
     641                 :            : GEN
     642                 :     228419 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     643                 :            : {
     644 [ +  + ][ +  +  :     228419 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
                +  +  - ]
     645                 :            :   {
     646                 :     210934 :     case 0: return Fp_div(x,y,p);
     647                 :      14080 :     case 1: return FpX_Fp_mul(x,Fp_inv(y,p),p);
     648                 :         75 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     649                 :       3330 :     case 3: return FpXQ_div(x,y,pol,p);
     650                 :            :   }
     651                 :     228419 :   return NULL;
     652                 :            : }
     653                 :            : 
     654                 :            : GEN
     655                 :       4010 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     656                 :            : {
     657         [ +  + ]:       4010 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     658                 :       3980 :   return FpXQ_pow(x,n,pol,p);
     659                 :            : }
     660                 :            : 
     661                 :            : GEN
     662                 :       9540 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     663                 :            : {
     664         [ +  + ]:       9540 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     665                 :       9540 :   return FpXQ_powu(x,n,pol,p);
     666                 :            : }
     667                 :            : 
     668                 :            : GEN
     669                 :       4265 : Fq_sqrt(GEN x, GEN T, GEN p)
     670                 :            : {
     671         [ +  + ]:       4265 :   if (typ(x) == t_INT)
     672                 :            :   {
     673 [ +  - ][ -  + ]:        240 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     674                 :        240 :     x = scalarpol_shallow(x, get_FpX_var(T));
     675                 :            :   }
     676                 :       4265 :   return FpXQ_sqrt(x,T,p);
     677                 :            : }
     678                 :            : GEN
     679                 :        510 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     680                 :            : {
     681         [ -  + ]:        510 :   if (typ(x) == t_INT)
     682                 :            :   {
     683                 :            :     long d;
     684         [ #  # ]:          0 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     685                 :          0 :     d = get_FpX_degree(T);
     686         [ #  # ]:          0 :     if (ugcd(umodiu(n,d),d) == 1)
     687                 :            :     {
     688         [ #  # ]:          0 :       if (!zeta)
     689                 :          0 :         return Fp_sqrtn(x,n,p,NULL);
     690                 :            :       else
     691                 :            :       {
     692                 :            :         /* gcd(n,p-1)=gcd(n,p^d-1) <=> same number of solutions if Fp and F_{p^d} */
     693         [ #  # ]:          0 :         if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     694                 :          0 :           return Fp_sqrtn(x,n,p,zeta);
     695                 :            :       }
     696                 :            :     }
     697                 :          0 :     x = scalarpol_shallow(x, get_FpX_var(T));
     698                 :            :   }
     699                 :        510 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     700                 :            : }
     701                 :            : 
     702                 :            : struct _Fq_field
     703                 :            : {
     704                 :            :   GEN T, p;
     705                 :            : };
     706                 :            : 
     707                 :            : static GEN
     708                 :      32022 : _Fq_red(void *E, GEN x)
     709                 :      32022 : { struct _Fq_field *s = (struct _Fq_field *)E;
     710                 :      32022 :   return Fq_red(x, s->T, s->p);
     711                 :            : }
     712                 :            : 
     713                 :            : static GEN
     714                 :     421704 : _Fq_add(void *E, GEN x, GEN y)
     715                 :            : {
     716                 :            :   (void) E;
     717         [ +  - ]:     421704 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
           [ -  -  +  + ]
     718                 :            :   {
     719                 :          0 :     case 0: return addii(x,y);
     720                 :          0 :     case 1: return ZX_Z_add(x,y);
     721                 :         51 :     case 2: return ZX_Z_add(y,x);
     722                 :     421704 :     default: return ZX_add(x,y);
     723                 :            :   }
     724                 :            : }
     725                 :            : 
     726                 :            : static GEN
     727         [ +  - ]:        868 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     728                 :            : 
     729                 :            : static GEN
     730                 :     436468 : _Fq_mul(void *E, GEN x, GEN y)
     731                 :            : {
     732                 :            :   (void) E;
     733         [ +  + ]:     436468 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
           [ -  +  +  + ]
     734                 :            :   {
     735                 :          0 :     case 0: return mulii(x,y);
     736                 :         30 :     case 1: return ZX_Z_mul(x,y);
     737                 :          5 :     case 2: return ZX_Z_mul(y,x);
     738                 :     436468 :     default: return ZX_mul(x,y);
     739                 :            :   }
     740                 :            : }
     741                 :            : 
     742                 :            : static GEN
     743                 :        818 : _Fq_inv(void *E, GEN x)
     744                 :        818 : { struct _Fq_field *s = (struct _Fq_field *)E;
     745                 :        818 :   return Fq_inv(x,s->T,s->p);
     746                 :            : }
     747                 :            : 
     748                 :            : static int
     749                 :      17157 : _Fq_equal0(GEN x) { return signe(x)==0; }
     750                 :            : 
     751                 :            : static GEN
     752                 :        832 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     753                 :            : 
     754                 :            : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     755                 :            :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     756                 :            : 
     757                 :         91 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     758                 :            : {
     759                 :         91 :   GEN z = new_chunk(sizeof(struct _Fq_field));
     760                 :         91 :   struct _Fq_field *e = (struct _Fq_field *) z;
     761                 :         91 :   e->T = T; e->p  = p; *E = (void*)e;
     762                 :         91 :   return &Fq_field;
     763                 :            : }
     764                 :            : 
     765                 :            : /*******************************************************************/
     766                 :            : /*                                                                 */
     767                 :            : /*                             Fq[X]                               */
     768                 :            : /*                                                                 */
     769                 :            : /*******************************************************************/
     770                 :            : /* P(X + c) */
     771                 :            : GEN
     772                 :          0 : FpX_translate(GEN P, GEN c, GEN p)
     773                 :            : {
     774                 :          0 :   pari_sp av = avma;
     775                 :            :   GEN Q, *R;
     776                 :            :   long i, k, n;
     777                 :            : 
     778 [ #  # ][ #  # ]:          0 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     779                 :          0 :   Q = leafcopy(P);
     780                 :          0 :   R = (GEN*)(Q+2); n = degpol(P);
     781         [ #  # ]:          0 :   for (i=1; i<=n; i++)
     782                 :            :   {
     783         [ #  # ]:          0 :     for (k=n-i; k<n; k++)
     784                 :          0 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     785                 :            : 
     786         [ #  # ]:          0 :     if (gc_needed(av,2))
     787                 :            :     {
     788         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     789                 :          0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     790                 :            :     }
     791                 :            :   }
     792                 :          0 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     793                 :            : }
     794                 :            : /* P(X + c), c an Fq */
     795                 :            : GEN
     796                 :      26700 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     797                 :            : {
     798                 :      26700 :   pari_sp av = avma;
     799                 :            :   GEN Q, *R;
     800                 :            :   long i, k, n;
     801                 :            : 
     802                 :            :   /* signe works for t_(INT|POL) */
     803 [ +  - ][ +  + ]:      26700 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     804                 :      26675 :   Q = leafcopy(P);
     805                 :      26675 :   R = (GEN*)(Q+2); n = degpol(P);
     806         [ +  + ]:     124230 :   for (i=1; i<=n; i++)
     807                 :            :   {
     808         [ +  + ]:     551335 :     for (k=n-i; k<n; k++)
     809                 :     453780 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     810                 :            : 
     811         [ -  + ]:      97555 :     if (gc_needed(av,2))
     812                 :            :     {
     813         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     814                 :          0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     815                 :            :     }
     816                 :            :   }
     817                 :      26700 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     818                 :            : }
     819                 :            : 
     820                 :            : GEN
     821                 :        280 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     822                 :            : {
     823                 :        280 :   pari_sp ltop = avma;
     824                 :            :   long k;
     825                 :            :   GEN W;
     826         [ +  + ]:        280 :   if (lgefint(p) == 3)
     827                 :            :   {
     828                 :        234 :     ulong pp = p[2];
     829                 :        234 :     GEN Tl = ZX_to_Flx(T, pp);
     830                 :        234 :     GEN Vl = FqV_to_FlxV(V, T, p);
     831                 :        234 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     832                 :        234 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     833                 :            :   }
     834                 :         46 :   W = cgetg(lg(V),t_VEC);
     835         [ +  + ]:        294 :   for(k=1; k < lg(V); k++)
     836                 :        248 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     837                 :        280 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     838                 :            : }
     839                 :            : 
     840                 :            : GEN
     841                 :       3430 : FqV_red(GEN z, GEN T, GEN p)
     842                 :            : {
     843                 :       3430 :   long i, l = lg(z);
     844                 :       3430 :   GEN res = cgetg(l, typ(z));
     845         [ +  + ]:      29280 :   for(i=1;i<l;i++) gel(res,i) = Fq_red(gel(z,i),T,p);
     846                 :       3430 :   return res;
     847                 :            : }
     848                 :            : 
     849                 :            : GEN
     850                 :          0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     851                 :            : {
     852                 :          0 :   long i, lx = lg(x);
     853                 :            :   GEN z;
     854         [ #  # ]:          0 :   if (!T) return FpC_add(x, y, p);
     855                 :          0 :   z = cgetg(lx, t_COL);
     856         [ #  # ]:          0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_add(gel(x, i), gel(y, i), T, p);
     857                 :          0 :   return z;
     858                 :            : }
     859                 :            : 
     860                 :            : GEN
     861                 :          0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     862                 :            : {
     863                 :          0 :   long i, lx = lg(x);
     864                 :            :   GEN z;
     865         [ #  # ]:          0 :   if (!T) return FpC_sub(x, y, p);
     866                 :          0 :   z = cgetg(lx, t_COL);
     867         [ #  # ]:          0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_sub(gel(x, i), gel(y, i), T, p);
     868                 :          0 :   return z;
     869                 :            : }
     870                 :            : 
     871                 :            : GEN
     872                 :          0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     873                 :            : {
     874                 :          0 :   long i, l = lg(x);
     875                 :            :   GEN z;
     876         [ #  # ]:          0 :   if (!T) return FpC_Fp_mul(x, y, p);
     877                 :          0 :   z = cgetg(l, t_COL);
     878         [ #  # ]:          0 :   for (i=1;i<l;i++) gel(z,i) = Fq_mul(gel(x,i),y,T,p);
     879                 :          0 :   return z;
     880                 :            : }
     881                 :            : 
     882                 :            : GEN
     883                 :        234 : FqV_to_FlxV(GEN v, GEN T, GEN pp)
     884                 :            : {
     885                 :        234 :   long j, N = lg(v);
     886                 :        234 :   long vT = varn(T);
     887                 :        234 :   ulong p = pp[2];
     888                 :        234 :   GEN y = cgetg(N, t_VEC);
     889         [ +  + ]:       1351 :   for (j=1; j<N; j++)
     890                 :       2234 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     891         [ -  + ]:       1117 :                                     : ZX_to_Flx(gel(v,j), p));
     892                 :        234 :   return y;
     893                 :            : }
     894                 :            : 
     895                 :            : GEN
     896                 :        515 : FqC_to_FlxC(GEN v, GEN T, GEN pp)
     897                 :            : {
     898                 :        515 :   long j, N = lg(v);
     899                 :        515 :   long vT = get_FpX_var(T);
     900                 :        515 :   ulong p = pp[2];
     901                 :        515 :   GEN y = cgetg(N, t_COL);
     902         [ +  + ]:      15835 :   for (j=1; j<N; j++)
     903                 :      33300 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     904         [ +  + ]:      15320 :                                     : ZX_to_Flx(gel(v,j), p));
     905                 :        515 :   return y;
     906                 :            : }
     907                 :            : 
     908                 :            : GEN
     909                 :         80 : FqM_to_FlxM(GEN x, GEN T, GEN pp)
     910                 :            : {
     911                 :         80 :   long j, n = lg(x);
     912                 :         80 :   GEN y = cgetg(n,t_MAT);
     913         [ -  + ]:         80 :   if (n == 1) return y;
     914         [ +  + ]:        595 :   for (j=1; j<n; j++)
     915                 :        515 :     gel(y,j) = FqC_to_FlxC(gel(x,j), T, pp);
     916                 :         80 :   return y;
     917                 :            : }
     918                 :            : 
     919                 :            : 
     920                 :            : /*******************************************************************/
     921                 :            : /*  Isomorphisms between finite fields                             */
     922                 :            : /*******************************************************************/
     923                 :            : 
     924                 :            : /* compute the reciprocical isomorphism of S mod T,p, i.e. V such that
     925                 :            :    V(S)=X  mod T,p*/
     926                 :            : 
     927                 :            : GEN
     928                 :         20 : Flxq_ffisom_inv(GEN S,GEN T, ulong p)
     929                 :            : {
     930                 :         20 :   pari_sp ltop = avma;
     931                 :         20 :   long n = degpol(T);
     932                 :         20 :   GEN M = Flxq_matrix_pow(S,n,n,T,p);
     933                 :         20 :   GEN V = Flm_Flc_invimage(M, vecsmall_ei(n, 2), p);
     934                 :         20 :   return gerepileupto(ltop, Flv_to_Flx(V, T[1]));
     935                 :            : }
     936                 :            : 
     937                 :            : GEN
     938                 :        905 : FpXQ_ffisom_inv(GEN S,GEN T, GEN p)
     939                 :            : {
     940                 :        905 :   pari_sp ltop = avma;
     941                 :        905 :   long n = degpol(T);
     942                 :        905 :   GEN V, M = FpXQ_matrix_pow(S,n,n,T,p);
     943                 :        905 :   V = FpM_FpC_invimage(M, col_ei(n, 2), p);
     944                 :        905 :   return gerepilecopy(ltop, RgV_to_RgX(V, varn(T)));
     945                 :            : }
     946                 :            : 
     947                 :            : /* Let M the matrix of the x^p Frobenius automorphism.
     948                 :            :  * Compute x^(p^i) for i=0..r */
     949                 :            : static GEN
     950                 :         20 : FpM_Frobenius(GEN M, long r, GEN p, long v)
     951                 :            : {
     952                 :         20 :   GEN W, V = cgetg(r+2,t_VEC);
     953                 :            :   long i;
     954         [ -  + ]:         20 :   gel(V,1) = pol_x(v); if (!r) return V;
     955                 :         20 :   gel(V,2) = RgV_to_RgX(gel(M,2),v);
     956                 :         20 :   W = gel(M,2);
     957         [ +  + ]:        200 :   for (i = 3; i <= r+1; ++i)
     958                 :            :   {
     959                 :        180 :     W = FpM_FpC_mul(M,W,p);
     960                 :        180 :     gel(V,i) = RgV_to_RgX(W,v);
     961                 :            :   }
     962                 :         20 :   return V;
     963                 :            : }
     964                 :            : 
     965                 :            : /* Let M the matrix of the x^p Frobenius automorphism.
     966                 :            :  * Compute x^(p^i) for i=0..r */
     967                 :            : static GEN
     968                 :       1390 : Flm_Frobenius(GEN M, long r, ulong p, long v)
     969                 :            : {
     970                 :       1390 :   GEN W, V = cgetg(r+2,t_VEC);
     971                 :            :   long i;
     972         [ -  + ]:       1390 :   gel(V,1) = polx_Flx(v); if (!r) return V;
     973                 :       1390 :   gel(V,2) = Flv_to_Flx(gel(M,2),v);
     974                 :       1390 :   W = gel(M,2);
     975         [ +  + ]:       4930 :   for (i = 3; i <= r+1; ++i)
     976                 :            :   {
     977                 :       3540 :     W = Flm_Flc_mul(M,W,p);
     978                 :       3540 :     gel(V,i) = Flv_to_Flx(W,v);
     979                 :            :   }
     980                 :       1390 :   return V;
     981                 :            : }
     982                 :            : 
     983                 :            : /* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in
     984                 :            :  * FFp[X]/T. Compute P(M)
     985                 :            :  * V=FpM_Frobenius(M, p, degpol(P), v)
     986                 :            :  * not stack clean
     987                 :            :  */
     988                 :            : 
     989                 :            : static GEN
     990                 :         20 : FpXQV_FpX_Frobenius(GEN V, GEN P, GEN T, GEN p)
     991                 :            : {
     992                 :            :   pari_sp btop;
     993                 :            :   long i;
     994                 :         20 :   long l = get_FpX_degree(T);
     995                 :         20 :   long v = get_FpX_var(T);
     996                 :            :   GEN M,W,Mi;
     997                 :            :   GEN *gptr[2];
     998                 :         20 :   long lV=lg(V);
     999                 :         20 :   GEN  PV=RgX_to_RgC(P, lgpol(P));
    1000                 :         20 :   M=cgetg(l+1,t_VEC);
    1001                 :         20 :   gel(M,1) = scalar_ZX_shallow(FpX_eval(P,gen_1,p),v);
    1002                 :         20 :   gel(M,2) = FpXV_FpC_mul(V,PV,p);
    1003                 :         20 :   btop=avma;
    1004                 :         20 :   gptr[0]=&Mi;
    1005                 :         20 :   gptr[1]=&W;
    1006                 :         20 :   W = leafcopy(V);
    1007         [ +  + ]:        200 :   for(i=3;i<=l;i++)
    1008                 :            :   {
    1009                 :            :     long j;
    1010                 :            :     pari_sp bbot;
    1011                 :        180 :     GEN W2=cgetg(lV,t_VEC);
    1012         [ +  + ]:       3440 :     for(j=1;j<lV;j++)
    1013                 :       3260 :       gel(W2,j) = FpXQ_mul(gel(W,j),gel(V,j),T,p);
    1014                 :        180 :     bbot=avma;
    1015                 :        180 :     Mi=FpXV_FpC_mul(W2,PV,p);
    1016                 :        180 :     W=gcopy(W2);
    1017                 :        180 :     gerepilemanysp(btop,bbot,gptr,2);
    1018                 :        180 :     btop=(pari_sp)W;
    1019                 :        180 :     gel(M,i) = Mi;
    1020                 :            :   }
    1021                 :         20 :   return RgXV_to_RgM(M,l);
    1022                 :            : }
    1023                 :            : 
    1024                 :            : static GEN
    1025                 :       1390 : FlxqV_Flx_Frobenius(GEN V, GEN P, GEN T, ulong p)
    1026                 :            : {
    1027                 :            :   pari_sp btop;
    1028                 :            :   long i;
    1029                 :       1390 :   long l = get_Flx_degree(T);
    1030                 :       1390 :   long v = get_Flx_var(T);
    1031                 :            :   GEN M,W,Mi;
    1032                 :       1390 :   GEN PV=Flx_to_Flv(P, lgpol(P));
    1033                 :            :   GEN *gptr[2];
    1034                 :       1390 :   long lV=lg(V);
    1035                 :       1390 :   M=cgetg(l+1,t_VEC);
    1036                 :       1390 :   gel(M,1) = Fl_to_Flx(Flx_eval(P,1,p),v);
    1037                 :       1390 :   gel(M,2) = FlxV_Flc_mul(V,PV,p);
    1038                 :       1390 :   btop=avma;
    1039                 :       1390 :   gptr[0]=&Mi;
    1040                 :       1390 :   gptr[1]=&W;
    1041                 :       1390 :   W=gcopy(V);
    1042         [ +  + ]:       8660 :   for(i=3;i<=l;i++)
    1043                 :            :   {
    1044                 :            :     long j;
    1045                 :            :     pari_sp bbot;
    1046                 :       7270 :     GEN W2=cgetg(lV,t_VEC);
    1047         [ +  + ]:     138890 :     for(j=1;j<lV;j++)
    1048                 :     131620 :       gel(W2,j) = Flxq_mul(gel(W,j),gel(V,j),T,p);
    1049                 :       7270 :     bbot=avma;
    1050                 :       7270 :     Mi=FlxV_Flc_mul(W2,PV,p);
    1051                 :       7270 :     W=gcopy(W2);
    1052                 :       7270 :     gerepilemanysp(btop,bbot,gptr,2);
    1053                 :       7270 :     btop=(pari_sp)W;
    1054                 :       7270 :     gel(M,i) = Mi;
    1055                 :            :   }
    1056                 :       1390 :   return FlxV_to_Flm(M,l);
    1057                 :            : }
    1058                 :            : 
    1059                 :            : /* Let M the matrix of the Frobenius automorphism of Fp[X]/(T).
    1060                 :            :  * Compute M^d
    1061                 :            :  * TODO: use left-right binary (tricky!)
    1062                 :            :  */
    1063                 :            : GEN
    1064                 :          0 : Flm_Frobenius_pow(GEN M, long d, GEN T, ulong p)
    1065                 :            : {
    1066                 :          0 :   pari_sp ltop=avma;
    1067                 :          0 :   long i,l=degpol(T);
    1068                 :          0 :   GEN R, W = gel(M,2);
    1069         [ #  # ]:          0 :   for (i = 2; i <= d; ++i) W = Flm_Flc_mul(M,W,p);
    1070                 :          0 :   R=Flxq_matrix_pow(Flv_to_Flx(W,T[2]),l,l,T,p);
    1071                 :          0 :   return gerepileupto(ltop,R);
    1072                 :            : }
    1073                 :            : 
    1074                 :            : GEN
    1075                 :          0 : FpM_Frobenius_pow(GEN M, long d, GEN T, GEN p)
    1076                 :            : {
    1077                 :          0 :   pari_sp ltop=avma;
    1078                 :          0 :   long i,l=degpol(T);
    1079                 :          0 :   GEN R, W = gel(M,2);
    1080         [ #  # ]:          0 :   for (i = 2; i <= d; ++i) W = FpM_FpC_mul(M,W,p);
    1081                 :          0 :   R=FpXQ_matrix_pow(RgV_to_RgX(W,varn(T)),l,l,T,p);
    1082                 :          0 :   return gerepilecopy(ltop,R);
    1083                 :            : }
    1084                 :            : 
    1085                 :            : /* Essentially we want to compute
    1086                 :            :  * FqM_ker(MA-pol_x(v),U,l)
    1087                 :            :  * To avoid use of matrix in Fq we procede as follows:
    1088                 :            :  * We compute FpM_ker(U(MA),l) and then we recover
    1089                 :            :  * the eigen value by Galois action, see formula.
    1090                 :            :  */
    1091                 :            : 
    1092                 :            : static GEN
    1093                 :       1390 : Flx_intersect_ker(GEN P, GEN MA, GEN U, ulong p)
    1094                 :            : {
    1095                 :       1390 :   pari_sp ltop = avma;
    1096                 :       1390 :   long i, vp = P[1], vu = U[1], r = degpol(U);
    1097                 :            :   GEN A, R;
    1098                 :            :   ulong ib0;
    1099                 :            :   pari_timer T;
    1100                 :            :   GEN M, V;
    1101         [ -  + ]:       1390 :   if (DEBUGLEVEL>=4) timer_start(&T);
    1102                 :       1390 :   V = Flm_Frobenius(MA, r, p, U[1]);
    1103         [ -  + ]:       1390 :   if (DEBUGLEVEL>=4) timer_printf(&T,"pol[Frobenius]");
    1104                 :       1390 :   M = FlxqV_Flx_Frobenius(V, U, P, p);
    1105         [ +  + ]:       1390 :   if (p==2)
    1106                 :        880 :     A = F2m_to_Flm(F2m_ker(Flm_to_F2m(M)));
    1107                 :            :   else
    1108                 :        510 :     A = Flm_ker(M,p);
    1109         [ -  + ]:       1390 :   if (DEBUGLEVEL>=4) timer_printf(&T,"matrix polcyclo");
    1110         [ -  + ]:       1390 :   if (lg(A)!=r+1) pari_err_IRREDPOL("FpX_ffintersect", Flx_to_ZX(P));
    1111                 :       1390 :   A = gerepileupto(ltop,A);
    1112                 :            :   /*The formula is
    1113                 :            :    * a_{r-1} = -\phi(a_0)/b_0
    1114                 :            :    * a_{i-1} = \phi(a_i)+b_ia_{r-1}  i=r-1 to 1
    1115                 :            :    * Where a_0=A[1] and b_i=U[i+2] */
    1116                 :       1390 :   ib0 = Fl_inv(Fl_neg(U[2], p), p);
    1117                 :       1390 :   R = cgetg(r+1,t_MAT);
    1118                 :       1390 :   gel(R,1) = gel(A,1);
    1119                 :       1390 :   gel(R,r) = Flm_Flc_mul(MA, Flc_Fl_mul(gel(A,1),ib0, p), p);
    1120         [ +  + ]:       3540 :   for(i=r-1; i>1; i--)
    1121                 :       2150 :     gel(R,i) = Flv_add(Flm_Flc_mul(MA,gel(R,i+1),p),
    1122                 :       4300 :                        Flc_Fl_mul(gel(R,r), U[i+2], p), p);
    1123                 :       1390 :   return gerepileupto(ltop, Flm_to_FlxX(Flm_transpose(R),vp,vu));
    1124                 :            : }
    1125                 :            : 
    1126                 :            : static GEN
    1127                 :       1020 : FpX_intersect_ker(GEN P, GEN MA, GEN U, GEN l)
    1128                 :            : {
    1129                 :       1020 :   pari_sp ltop = avma;
    1130                 :       1020 :   long i, vp = varn(P), vu = varn(U), r = degpol(U);
    1131                 :            :   GEN V, A, R, ib0;
    1132                 :            :   pari_timer T;
    1133         [ +  + ]:       1020 :   if (lgefint(l)==3)
    1134                 :            :   {
    1135                 :       1000 :     ulong p = l[2];
    1136                 :       1000 :     GEN res = Flx_intersect_ker(ZX_to_Flx(P,p), ZM_to_Flm(MA,p), ZX_to_Flx(U,p), p);
    1137                 :       1000 :     return gerepileupto(ltop, FlxX_to_ZXX(res));
    1138                 :            :   }
    1139         [ -  + ]:         20 :   if (DEBUGLEVEL>=4) timer_start(&T);
    1140                 :         20 :   V = FpM_Frobenius(MA,r,l,vu);
    1141         [ -  + ]:         20 :   if (DEBUGLEVEL>=4) timer_printf(&T,"pol[Frobenius]");
    1142                 :         20 :   A = FpM_ker(FpXQV_FpX_Frobenius(V, U, P, l), l);
    1143         [ -  + ]:         20 :   if (DEBUGLEVEL>=4) timer_printf(&T,"matrix polcyclo");
    1144         [ -  + ]:         20 :   if (lg(A)!=r+1) pari_err_IRREDPOL("FpX_ffintersect", P);
    1145                 :         20 :   A = gerepileupto(ltop,A);
    1146                 :            :   /*The formula is
    1147                 :            :    * a_{r-1} = -\phi(a_0)/b_0
    1148                 :            :    * a_{i-1} = \phi(a_i)+b_ia_{r-1}  i=r-1 to 1
    1149                 :            :    * Where a_0=A[1] and b_i=U[i+2] */
    1150                 :         20 :   ib0 = Fp_inv(negi(gel(U,2)),l);
    1151                 :         20 :   R = cgetg(r+1,t_MAT);
    1152                 :         20 :   gel(R,1) = gel(A,1);
    1153                 :         20 :   gel(R,r) = FpM_FpC_mul(MA, FpC_Fp_mul(gel(A,1),ib0,l), l);
    1154         [ +  + ]:        180 :   for(i=r-1;i>1;i--)
    1155                 :        160 :     gel(R,i) = FpC_add(FpM_FpC_mul(MA,gel(R,i+1),l),
    1156                 :        320 :         FpC_Fp_mul(gel(R,r), gel(U,i+2), l),l);
    1157                 :       1020 :   return gerepilecopy(ltop,RgM_to_RgXX(shallowtrans(R),vp,vu));
    1158                 :            : }
    1159                 :            : 
    1160                 :            : /* n must divide both the degree of P and Q.  Compute SP and SQ such
    1161                 :            :   that the subfield of FF_l[X]/(P) generated by SP and the subfield of
    1162                 :            :   FF_l[X]/(Q) generated by SQ are isomorphic of degree n.  P and Q do
    1163                 :            :   not need to be of the same variable.  if MA (resp. MB) is not NULL,
    1164                 :            :   must be the matrix of the Frobenius map in FF_l[X]/(P) (resp.
    1165                 :            :   FF_l[X]/(Q) ).  */
    1166                 :            : /* Note on the implementation choice:
    1167                 :            :  * We assume the prime p is very large
    1168                 :            :  * so we handle Frobenius as matrices.
    1169                 :            :  */
    1170                 :            : 
    1171                 :            : void
    1172                 :        650 : Flx_ffintersect(GEN P, GEN Q, long n, ulong l,GEN *SP, GEN *SQ, GEN MA, GEN MB)
    1173                 :            : {
    1174                 :        650 :   pari_sp ltop = avma;
    1175                 :        650 :   long vp = P[1], vq = Q[1], np = degpol(P), nq = degpol(Q), e;
    1176                 :            :   ulong pg;
    1177                 :            :   GEN A, B, Ap, Bp;
    1178         [ -  + ]:        650 :   if (np<=0) pari_err_IRREDPOL("FpX_ffintersect", P);
    1179         [ -  + ]:        650 :   if (nq<=0) pari_err_IRREDPOL("FpX_ffintersect", Q);
    1180 [ +  - ][ +  - ]:        650 :   if (n<=0 || np%n || nq%n)
                 [ -  + ]
    1181                 :          0 :     pari_err_TYPE("FpX_ffintersect [bad degrees]",stoi(n));
    1182                 :        650 :   e = u_lvalrem(n, l, &pg);
    1183         [ +  + ]:        650 :   if(!MA) MA = Flx_matFrobenius(P,l);
    1184         [ +  + ]:        650 :   if(!MB) MB = Flx_matFrobenius(Q,l);
    1185                 :        650 :   A = Ap = pol0_Flx(vp);
    1186                 :        650 :   B = Bp = pol0_Flx(vq);
    1187         [ +  + ]:        650 :   if (pg > 1)
    1188                 :            :   {
    1189                 :            :     pari_timer T;
    1190                 :        560 :     GEN ipg = utoipos(pg);
    1191         [ +  + ]:        560 :     if (l%pg == 1)
    1192                 :            :     /* No need to use relative extension, so don't. (Well, now we don't
    1193                 :            :      * in the other case either, but this special case is more efficient) */
    1194                 :            :     {
    1195                 :            :       GEN L;
    1196                 :            :       ulong z, An, Bn;
    1197                 :        365 :       z = Fl_neg(rootsof1_Fl(pg, l), l);
    1198         [ -  + ]:        365 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1199                 :        365 :       A = Flm_ker(Flm_Fl_add(MA, z, l),l);
    1200         [ -  + ]:        365 :       if (lg(A)!=2) pari_err_IRREDPOL("FpX_ffintersect",P);
    1201                 :        365 :       A = Flv_to_Flx(gel(A,1),vp);
    1202                 :            : 
    1203                 :        365 :       B = Flm_ker(Flm_Fl_add(MB, z, l),l);
    1204         [ -  + ]:        365 :       if (lg(B)!=2) pari_err_IRREDPOL("FpX_ffintersect",Q);
    1205                 :        365 :       B = Flv_to_Flx(gel(B,1),vq);
    1206                 :            : 
    1207         [ -  + ]:        365 :       if (DEBUGLEVEL>=4) timer_printf(&T, "FpM_ker");
    1208                 :        365 :       An = Flxq_powu(A,pg,P,l)[2];
    1209                 :        365 :       Bn = Flxq_powu(B,pg,Q,l)[2];
    1210         [ -  + ]:        365 :       if (!Bn) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1211                 :        365 :       z = Fl_div(An,Bn,l);
    1212                 :        365 :       L = Fp_sqrtn(utoi(z),ipg,utoipos(l),NULL);
    1213         [ -  + ]:        365 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1214         [ -  + ]:        365 :       if (DEBUGLEVEL>=4) timer_printf(&T, "Fp_sqrtn");
    1215                 :        365 :       B = Flx_Fl_mul(B,itou(L),l);
    1216                 :            :     }
    1217                 :            :     else
    1218                 :            :     {
    1219                 :            :       GEN L, An, Bn, z, U;
    1220                 :        195 :       U = gmael(Flx_factor(ZX_to_Flx(polcyclo(pg, fetch_var()),l),l),1,1);
    1221                 :        195 :       A = Flx_intersect_ker(P, MA, U, l);
    1222                 :        195 :       B = Flx_intersect_ker(Q, MB, U, l);
    1223         [ -  + ]:        195 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1224                 :        195 :       An = gel(FlxYqq_pow(A,ipg,P,U,l),2);
    1225                 :        195 :       Bn = gel(FlxYqq_pow(B,ipg,Q,U,l),2);
    1226         [ -  + ]:        195 :       if (DEBUGLEVEL>=4) timer_printf(&T,"pows [P,Q]");
    1227                 :        195 :       z = Flxq_div(An,Bn,U,l);
    1228                 :        195 :       L = Flxq_sqrtn(z,ipg,U,l,NULL);
    1229         [ -  + ]:        195 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1230         [ -  + ]:        195 :       if (DEBUGLEVEL>=4) timer_printf(&T,"FpXQ_sqrtn");
    1231                 :        195 :       B = FlxqX_Flxq_mul(B,L,U,l);
    1232                 :        195 :       A = FlxY_evalx(A,0,l);
    1233                 :        195 :       B = FlxY_evalx(B,0,l);
    1234                 :        560 :       (void)delete_var();
    1235                 :            :     }
    1236                 :            :   }
    1237         [ +  + ]:        650 :   if (e)
    1238                 :            :   {
    1239                 :            :     GEN VP, VQ, Ay, By;
    1240                 :        125 :     ulong lmun = l-1;
    1241                 :            :     long j;
    1242                 :        125 :     MA = Flm_Fl_add(MA,lmun,l);
    1243                 :        125 :     MB = Flm_Fl_add(MB,lmun,l);
    1244                 :        125 :     Ay = pol1_Flx(vp);
    1245                 :        125 :     By = pol1_Flx(vq);
    1246                 :        125 :     VP = vecsmall_ei(np, 1);
    1247         [ +  + ]:        125 :     VQ = np == nq? VP: vecsmall_ei(nq, 1); /* save memory */
    1248         [ +  + ]:        375 :     for(j=0;j<e;j++)
    1249                 :            :     {
    1250         [ +  + ]:        250 :       if (j)
    1251                 :            :       {
    1252                 :        125 :         Ay = Flxq_mul(Ay,Flxq_powu(Ap,lmun,P,l),P,l);
    1253                 :        125 :         VP = Flx_to_Flv(Ay,np);
    1254                 :            :       }
    1255                 :        250 :       Ap = Flm_Flc_invimage(MA,VP,l);
    1256                 :        250 :       Ap = Flv_to_Flx(Ap,vp);
    1257                 :            : 
    1258         [ +  + ]:        250 :       if (j)
    1259                 :            :       {
    1260                 :        125 :         By = Flxq_mul(By,Flxq_powu(Bp,lmun,Q,l),Q,l);
    1261                 :        125 :         VQ = Flx_to_Flv(By,nq);
    1262                 :            :       }
    1263                 :        250 :       Bp = Flm_Flc_invimage(MB,VQ,l);
    1264                 :        250 :       Bp = Flv_to_Flx(Bp,vq);
    1265                 :            :     }
    1266                 :            :   }
    1267                 :        650 :   *SP = Flx_add(A,Ap,l);
    1268                 :        650 :   *SQ = Flx_add(B,Bp,l);
    1269                 :        650 :   gerepileall(ltop,2,SP,SQ);
    1270                 :        650 : }
    1271                 :            : 
    1272                 :            : /* Let l be a prime number, P, Q in ZZ[X].  P and Q are both
    1273                 :            :  * irreducible modulo l and degree(P) divides degree(Q).  Output a
    1274                 :            :  * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
    1275                 :            :  * that Q | P(R) mod l.  If P and Q have the same degree, it is of course an
    1276                 :            :  * isomorphism.  */
    1277                 :            : GEN
    1278                 :         20 : Flx_ffisom(GEN P,GEN Q,ulong l)
    1279                 :            : {
    1280                 :         20 :   pari_sp av = avma;
    1281                 :            :   GEN SP, SQ, R;
    1282                 :         20 :   Flx_ffintersect(P,Q,degpol(P),l,&SP,&SQ,NULL,NULL);
    1283                 :         20 :   R = Flxq_ffisom_inv(SP,P,l);
    1284                 :         20 :   return gerepileupto(av, Flx_Flxq_eval(R,SQ,Q,l));
    1285                 :            : }
    1286                 :            : 
    1287                 :            : void
    1288                 :        830 : FpX_ffintersect(GEN P, GEN Q, long n, GEN l, GEN *SP, GEN *SQ, GEN MA, GEN MB)
    1289                 :            : {
    1290                 :        830 :   pari_sp ltop = avma;
    1291                 :            :   long vp, vq, np, nq, e;
    1292                 :            :   ulong pg;
    1293                 :            :   GEN A, B, Ap, Bp;
    1294                 :        830 :   vp = varn(P); np = degpol(P);
    1295                 :        830 :   vq = varn(Q); nq = degpol(Q);
    1296         [ -  + ]:        830 :   if (np<=0) pari_err_IRREDPOL("FpX_ffintersect", P);
    1297         [ -  + ]:        830 :   if (nq<=0) pari_err_IRREDPOL("FpX_ffintersect", Q);
    1298 [ +  - ][ +  - ]:        830 :   if (n<=0 || np%n || nq%n)
                 [ -  + ]
    1299                 :          0 :     pari_err_TYPE("FpX_ffintersect [bad degrees]",stoi(n));
    1300                 :        830 :   e = u_pvalrem(n, l, &pg);
    1301         [ +  + ]:        830 :   if(!MA) MA = FpX_matFrobenius(P, l);
    1302         [ +  + ]:        830 :   if(!MB) MB = FpX_matFrobenius(Q, l);
    1303                 :        830 :   A = Ap = pol_0(vp);
    1304                 :        830 :   B = Bp = pol_0(vq);
    1305         [ +  + ]:        830 :   if (pg > 1)
    1306                 :            :   {
    1307                 :        740 :     GEN ipg = utoipos(pg);
    1308                 :            :     pari_timer T;
    1309         [ +  + ]:        740 :     if (umodiu(l,pg) == 1)
    1310                 :            :     /* No need to use relative extension, so don't. (Well, now we don't
    1311                 :            :      * in the other case either, but this special case is more efficient) */
    1312                 :            :     {
    1313                 :            :       GEN L, An, Bn, z;
    1314                 :        230 :       z = negi( rootsof1u_Fp(pg, l) );
    1315         [ -  + ]:        230 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1316                 :        230 :       A = FpM_ker(RgM_Rg_add_shallow(MA, z),l);
    1317         [ -  + ]:        230 :       if (lg(A)!=2) pari_err_IRREDPOL("FpX_ffintersect",P);
    1318                 :        230 :       A = RgV_to_RgX(gel(A,1),vp);
    1319                 :            : 
    1320                 :        230 :       B = FpM_ker(RgM_Rg_add_shallow(MB, z),l);
    1321         [ -  + ]:        230 :       if (lg(B)!=2) pari_err_IRREDPOL("FpX_ffintersect",Q);
    1322                 :        230 :       B = RgV_to_RgX(gel(B,1),vq);
    1323                 :            : 
    1324         [ -  + ]:        230 :       if (DEBUGLEVEL>=4) timer_printf(&T, "FpM_ker");
    1325                 :        230 :       An = gel(FpXQ_pow(A,ipg,P,l),2);
    1326                 :        230 :       Bn = gel(FpXQ_pow(B,ipg,Q,l),2);
    1327         [ -  + ]:        230 :       if (!signe(Bn)) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1328                 :        230 :       z = Fp_div(An,Bn,l);
    1329                 :        230 :       L = Fp_sqrtn(z,ipg,l,NULL);
    1330         [ -  + ]:        230 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1331         [ -  + ]:        230 :       if (DEBUGLEVEL>=4) timer_printf(&T, "Fp_sqrtn");
    1332                 :        230 :       B = FpX_Fp_mul(B,L,l);
    1333                 :            :     }
    1334                 :            :     else
    1335                 :            :     {
    1336                 :            :       GEN L, An, Bn, z, U;
    1337                 :        510 :       U = gmael(FpX_factor(polcyclo(pg,fetch_var()),l),1,1);
    1338                 :        510 :       A = FpX_intersect_ker(P, MA, U, l);
    1339                 :        510 :       B = FpX_intersect_ker(Q, MB, U, l);
    1340         [ -  + ]:        510 :       if (DEBUGLEVEL>=4) timer_start(&T);
    1341                 :        510 :       An = gel(FpXYQQ_pow(A,ipg,P,U,l),2);
    1342                 :        510 :       Bn = gel(FpXYQQ_pow(B,ipg,Q,U,l),2);
    1343         [ -  + ]:        510 :       if (DEBUGLEVEL>=4) timer_printf(&T,"pows [P,Q]");
    1344         [ -  + ]:        510 :       if (!signe(Bn)) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1345                 :        510 :       z = Fq_div(An,Bn,U,l);
    1346                 :        510 :       L = Fq_sqrtn(z,ipg,U,l,NULL);
    1347         [ -  + ]:        510 :       if (!L) pari_err_IRREDPOL("FpX_ffintersect", mkvec2(P,Q));
    1348         [ -  + ]:        510 :       if (DEBUGLEVEL>=4) timer_printf(&T,"FpXQ_sqrtn");
    1349                 :        510 :       B = FqX_Fq_mul(B,L,U,l);
    1350                 :        510 :       A = FpXY_evalx(A,gen_0,l);
    1351                 :        510 :       B = FpXY_evalx(B,gen_0,l);
    1352                 :        740 :       (void)delete_var();
    1353                 :            :     }
    1354                 :            :   }
    1355         [ +  + ]:        830 :   if (e)
    1356                 :            :   {
    1357                 :         95 :     GEN VP, VQ, Ay, By, lmun = addis(l,-1);
    1358                 :            :     long j;
    1359                 :         95 :     MA = RgM_Rg_add_shallow(MA,gen_m1);
    1360                 :         95 :     MB = RgM_Rg_add_shallow(MB,gen_m1);
    1361                 :         95 :     Ay = pol_1(vp);
    1362                 :         95 :     By = pol_1(vq);
    1363                 :         95 :     VP = col_ei(np, 1);
    1364         [ +  + ]:         95 :     VQ = np == nq? VP: col_ei(nq, 1); /* save memory */
    1365         [ +  + ]:        245 :     for(j=0;j<e;j++)
    1366                 :            :     {
    1367         [ +  + ]:        150 :       if (j)
    1368                 :            :       {
    1369                 :         55 :         Ay = FpXQ_mul(Ay,FpXQ_pow(Ap,lmun,P,l),P,l);
    1370                 :         55 :         VP = RgX_to_RgC(Ay,np);
    1371                 :            :       }
    1372                 :        150 :       Ap = FpM_FpC_invimage(MA,VP,l);
    1373                 :        150 :       Ap = RgV_to_RgX(Ap,vp);
    1374                 :            : 
    1375         [ +  + ]:        150 :       if (j)
    1376                 :            :       {
    1377                 :         55 :         By = FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l);
    1378                 :         55 :         VQ = RgX_to_RgC(By,nq);
    1379                 :            :       }
    1380                 :        150 :       Bp = FpM_FpC_invimage(MB,VQ,l);
    1381                 :        150 :       Bp = RgV_to_RgX(Bp,vq);
    1382                 :            :     }
    1383                 :            :   }
    1384                 :        830 :   *SP = FpX_add(A,Ap,l);
    1385                 :        830 :   *SQ = FpX_add(B,Bp,l);
    1386                 :        830 :   gerepileall(ltop,2,SP,SQ);
    1387                 :        830 : }
    1388                 :            : /* Let l be a prime number, P, Q in ZZ[X].  P and Q are both
    1389                 :            :  * irreducible modulo l and degree(P) divides degree(Q).  Output a
    1390                 :            :  * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
    1391                 :            :  * that Q | P(R) mod l.  If P and Q have the same degree, it is of course an
    1392                 :            :  * isomorphism.  */
    1393                 :            : GEN
    1394                 :        815 : FpX_ffisom(GEN P,GEN Q,GEN l)
    1395                 :            : {
    1396                 :        815 :   pari_sp av = avma;
    1397                 :            :   GEN SP, SQ, R;
    1398                 :        815 :   FpX_ffintersect(P,Q,degpol(P),l,&SP,&SQ,NULL,NULL);
    1399                 :        815 :   R = FpXQ_ffisom_inv(SP,P,l);
    1400                 :        815 :   return gerepileupto(av, FpX_FpXQ_eval(R,SQ,Q,l));
    1401                 :            : }
    1402                 :            : 
    1403                 :            : /* Let l be a prime number, P a ZX irreducible modulo l, MP the matrix of the
    1404                 :            :  * Frobenius automorphism of F_l[X]/(P).
    1405                 :            :  * Factor P over the subfield of F_l[X]/(P) of index d. */
    1406                 :            : static GEN
    1407                 :         15 : FpX_factorgalois(GEN P, GEN l, long d, long w, GEN MP)
    1408                 :            : {
    1409                 :         15 :   pari_sp ltop = avma;
    1410                 :            :   GEN R, V, Tl, z, M;
    1411                 :         15 :   long k, n = degpol(P), m = n/d;
    1412                 :         15 :   long v = varn(P);
    1413                 :            : 
    1414                 :            :   /* x - y */
    1415         [ +  - ]:         15 :   if (m == 1) return deg1pol_shallow(gen_1, deg1pol_shallow(subis(l,1), gen_0, w), v);
    1416                 :          0 :   M = FpM_Frobenius_pow(MP,d,P,l);
    1417                 :            : 
    1418                 :          0 :   Tl = leafcopy(P); setvarn(Tl,w);
    1419                 :          0 :   V = cgetg(m+1,t_VEC);
    1420                 :          0 :   gel(V,1) = pol_x(w);
    1421                 :          0 :   z = gel(M,2);
    1422                 :          0 :   gel(V,2) = RgV_to_RgX(z,w);
    1423         [ #  # ]:          0 :   for(k=3;k<=m;k++)
    1424                 :            :   {
    1425                 :          0 :     z = FpM_FpC_mul(M,z,l);
    1426                 :          0 :     gel(V,k) = RgV_to_RgX(z,w);
    1427                 :            :   }
    1428                 :          0 :   R = FqV_roots_to_pol(V,Tl,l,v);
    1429                 :         15 :   return gerepileupto(ltop,R);
    1430                 :            : }
    1431                 :            : /* same: P is an Flx, MP an Flm */
    1432                 :            : static GEN
    1433                 :        630 : Flx_factorgalois(GEN P, ulong l, long d, long w, GEN MP)
    1434                 :            : {
    1435                 :        630 :   pari_sp ltop = avma;
    1436                 :            :   GEN R, V, Tl, z, M;
    1437                 :        630 :   long k, n = degpol(P), m = n/d;
    1438                 :        630 :   long v = P[1];
    1439                 :            : 
    1440         [ +  - ]:        630 :   if (m == 1) {
    1441                 :        630 :     R = polx_Flx(v);
    1442                 :        630 :     gel(R,2) = z = polx_Flx(w); z[3] = l - 1; /* - y */
    1443                 :        630 :     gel(R,3) = pol1_Flx(w);
    1444                 :        630 :     return R; /* x - y */
    1445                 :            :   }
    1446                 :          0 :   M = Flm_Frobenius_pow(MP,d,P,l);
    1447                 :            : 
    1448                 :          0 :   Tl = leafcopy(P); setvarn(Tl,w);
    1449                 :          0 :   V = cgetg(m+1,t_VEC);
    1450                 :          0 :   gel(V,1) = polx_Flx(w);
    1451                 :          0 :   z = gel(M,2);
    1452                 :          0 :   gel(V,2) = Flv_to_Flx(z,w);
    1453         [ #  # ]:          0 :   for(k=3;k<=m;k++)
    1454                 :            :   {
    1455                 :          0 :     z = Flm_Flc_mul(M,z,l);
    1456                 :          0 :     gel(V,k) = Flv_to_Flx(z,w);
    1457                 :            :   }
    1458                 :          0 :   R = FlxqV_roots_to_pol(V,Tl,l,v);
    1459                 :        630 :   return gerepileupto(ltop,R);
    1460                 :            : }
    1461                 :            : 
    1462                 :            : GEN
    1463                 :        630 : Flx_factorff_irred(GEN P, GEN Q, ulong p)
    1464                 :            : {
    1465                 :        630 :   pari_sp ltop = avma, av;
    1466                 :            :   GEN SP, SQ, MP, MQ, M, FP, FQ, E, V, IR, res;
    1467                 :        630 :   long np = degpol(P), nq = degpol(Q), d = cgcd(np,nq);
    1468                 :        630 :   long i, vp = P[1], vq = Q[1];
    1469         [ -  + ]:        630 :   if (d==1) retmkcol(Flx_to_FlxX(P, vq));
    1470                 :        630 :   FQ = Flx_matFrobenius(Q,p);
    1471                 :        630 :   av = avma;
    1472                 :        630 :   FP = Flx_matFrobenius(P,p);
    1473                 :        630 :   Flx_ffintersect(P,Q,d,p,&SP,&SQ, FP, FQ);
    1474                 :        630 :   E = Flx_factorgalois(P,p,d,vq, FP);
    1475                 :        630 :   E = FlxX_to_Flm(E,np);
    1476                 :        630 :   MP= Flxq_matrix_pow(SP,np,d,P,p);
    1477                 :        630 :   IR= gel(Flm_indexrank(MP,p),1);
    1478                 :        630 :   E = rowpermute(E, IR);
    1479                 :        630 :   M = rowpermute(MP,IR);
    1480                 :        630 :   M = Flm_inv(M,p);
    1481                 :        630 :   MQ= Flxq_matrix_pow(SQ,nq,d,Q,p);
    1482                 :        630 :   M = Flm_mul(MQ,M,p);
    1483                 :        630 :   M = Flm_mul(M,E,p);
    1484                 :        630 :   M = gerepileupto(av,M);
    1485                 :        630 :   V = cgetg(d+1,t_VEC);
    1486                 :        630 :   gel(V,1) = M;
    1487         [ +  + ]:       4930 :   for(i=2;i<=d;i++)
    1488                 :       4300 :     gel(V,i) = Flm_mul(FQ,gel(V,i-1),p);
    1489                 :        630 :   res = cgetg(d+1,t_COL);
    1490         [ +  + ]:       5560 :   for(i=1;i<=d;i++)
    1491                 :       4930 :     gel(res,i) = Flm_to_FlxX(gel(V,i),vp,vq);
    1492                 :        630 :   return gerepileupto(ltop,res);
    1493                 :            : }
    1494                 :            : 
    1495                 :            : /* P,Q irreducible over F_p. Factor P over FF_p[X] / Q  [factors are ordered as
    1496                 :            :  * a Frobenius cycle] */
    1497                 :            : GEN
    1498                 :       1440 : FpX_factorff_irred(GEN P, GEN Q, GEN p)
    1499                 :            : {
    1500                 :       1440 :   pari_sp ltop = avma, av;
    1501                 :            :   GEN res;
    1502                 :       1440 :   long np = degpol(P), nq = degpol(Q), d = cgcd(np,nq);
    1503         [ +  + ]:       1440 :   if (d==1) return mkcolcopy(P);
    1504                 :            : 
    1505         [ +  + ]:        645 :   if (lgefint(p)==3)
    1506                 :            :   {
    1507                 :        630 :     ulong pp = p[2];
    1508                 :        630 :     GEN F = Flx_factorff_irred(ZX_to_Flx(P,pp), ZX_to_Flx(Q,pp), pp);
    1509                 :        630 :     long i, lF = lg(F);
    1510                 :        630 :     res = cgetg(lF, t_COL);
    1511         [ +  + ]:       5560 :     for(i=1; i<lF; i++)
    1512                 :       4930 :       gel(res,i) = FlxX_to_ZXX(gel(F,i));
    1513                 :            :   }
    1514                 :            :   else
    1515                 :            :   {
    1516                 :            :     GEN SP, SQ, MP, MQ, M, FP, FQ, E, V, IR;
    1517                 :         15 :     long i, vp = varn(P), vq = varn(Q);
    1518                 :         15 :     FQ = FpX_matFrobenius(Q,p);
    1519                 :         15 :     av = avma;
    1520                 :         15 :     FP = FpX_matFrobenius(P,p);
    1521                 :         15 :     FpX_ffintersect(P,Q,d,p,&SP,&SQ,FP,FQ);
    1522                 :            : 
    1523                 :         15 :     E = FpX_factorgalois(P,p,d,vq,FP);
    1524                 :         15 :     E = RgXX_to_RgM(E,np);
    1525                 :         15 :     MP= FpXQ_matrix_pow(SP,np,d,P,p);
    1526                 :         15 :     IR= gel(FpM_indexrank(MP,p),1);
    1527                 :         15 :     E = rowpermute(E, IR);
    1528                 :         15 :     M = rowpermute(MP,IR);
    1529                 :         15 :     M = FpM_inv(M,p);
    1530                 :         15 :     MQ= FpXQ_matrix_pow(SQ,nq,d,Q,p);
    1531                 :         15 :     M = FpM_mul(MQ,M,p);
    1532                 :         15 :     M = FpM_mul(M,E,p);
    1533                 :         15 :     M = gerepileupto(av,M);
    1534                 :         15 :     V = cgetg(d+1,t_VEC);
    1535                 :         15 :     gel(V,1) = M;
    1536         [ +  + ]:        205 :     for(i=2;i<=d;i++)
    1537                 :        190 :       gel(V,i) = FpM_mul(FQ,gel(V,i-1),p);
    1538                 :         15 :     res = cgetg(d+1,t_COL);
    1539         [ +  + ]:        220 :     for(i=1;i<=d;i++)
    1540                 :        205 :       gel(res,i) = RgM_to_RgXX(gel(V,i),vp,vq);
    1541                 :            :   }
    1542                 :       1440 :   return gerepilecopy(ltop,res);
    1543                 :            : }
    1544                 :            : 
    1545                 :            : /*******************************************************************/
    1546                 :            : /*                                                                 */
    1547                 :            : /*                          MODULAR GCD                            */
    1548                 :            : /*                                                                 */
    1549                 :            : /*******************************************************************/
    1550                 :            : /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
    1551                 :            : static GEN
    1552                 :   21375055 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
    1553                 :            : {
    1554                 :   21375055 :   ulong d, amod = umodiu(a, p);
    1555                 :   21375055 :   pari_sp av = avma;
    1556                 :            :   GEN ax;
    1557                 :            : 
    1558         [ +  + ]:   21375055 :   if (b == amod) return NULL;
    1559         [ +  + ]:    8129836 :   d = (b > amod)? b - amod: p - (amod - b); /* (b - a) mod p */
    1560                 :    8129836 :   (void)new_chunk(lgefint(pq)<<1); /* HACK */
    1561                 :    8129836 :   ax = mului(Fl_mul(d,qinv,p), q); /* d mod p, 0 mod q */
    1562                 :   21375055 :   avma = av; return addii(a, ax); /* in ]-q, pq[ assuming a in -]-q,q[ */
    1563                 :            : }
    1564                 :            : GEN
    1565                 :      11895 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1566                 :            : GEN
    1567                 :    2224235 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1568                 :            : {
    1569                 :    2224235 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1570                 :    2224235 :   GEN H = cgetg(l, t_POL);
    1571                 :    2224235 :   H[1] = evalsigne(1) | evalvarn(v);
    1572         [ +  + ]:    7868712 :   for (i=2; i<l; i++)
    1573                 :    5644477 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1574                 :    2224235 :   return H;
    1575                 :            : }
    1576                 :            : 
    1577                 :            : /* assume lg(Hp) > 1 */
    1578                 :            : GEN
    1579                 :      65562 : ZM_init_CRT(GEN Hp, ulong p)
    1580                 :            : {
    1581                 :      65562 :   long i,j, m = lgcols(Hp), l = lg(Hp), lim = (long)(p>>1);
    1582                 :      65562 :   GEN c,cp,H = cgetg(l, t_MAT);
    1583         [ +  + ]:     550455 :   for (j=1; j<l; j++)
    1584                 :            :   {
    1585                 :     484893 :     cp = gel(Hp,j);
    1586                 :     484893 :     c = cgetg(m, t_COL);
    1587                 :     484893 :     gel(H,j) = c;
    1588         [ +  + ]:    4843310 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1589                 :            :   }
    1590                 :      65562 :   return H;
    1591                 :            : }
    1592                 :            : 
    1593                 :            : int
    1594                 :     156123 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1595                 :            : {
    1596                 :     156123 :   GEN h, q = *ptq, qp = muliu(q,p), lim = shifti(qp,-1);
    1597                 :     156123 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1598                 :     156123 :   int stable = 1;
    1599                 :     156123 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp);
    1600         [ +  + ]:     156123 :   if (h)
    1601                 :            :   {
    1602         [ +  + ]:      43405 :     if (cmpii(h,lim) > 0) h = subii(h,qp);
    1603                 :      43405 :     *H = h; stable = 0;
    1604                 :            :   }
    1605                 :     156123 :   *ptq = qp; return stable;
    1606                 :            : }
    1607                 :            : 
    1608                 :            : static int
    1609                 :    2450653 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1610                 :            : {
    1611                 :    2450653 :   GEN H = *ptH, h, lim = shifti(qp,-1);
    1612                 :    2450653 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1613                 :    2450653 :   long i, l = lg(H), lp = lg(Hp);
    1614                 :    2450653 :   int stable = 1;
    1615                 :            : 
    1616         [ +  + ]:    2450653 :   if (l < lp)
    1617                 :            :   { /* degree increases */
    1618                 :        300 :     GEN x = cgetg(lp, t_POL);
    1619         [ +  + ]:        600 :     for (i=1; i<l; i++)  x[i] = H[i];
    1620         [ +  + ]:       4810 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1621                 :        300 :     *ptH = H = x;
    1622                 :        300 :     stable = 0;
    1623         [ +  + ]:    2450353 :   } else if (l > lp)
    1624                 :            :   { /* degree decreases */
    1625                 :          5 :     GEN x = cgetg(l, t_VECSMALL);
    1626         [ +  + ]:         80 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1627         [ +  + ]:         10 :     for (   ; i<l; i++) x[i] = 0;
    1628                 :          5 :     Hp = x; lp = l;
    1629                 :            :   }
    1630         [ +  + ]:   12086926 :   for (i=2; i<lp; i++)
    1631                 :            :   {
    1632                 :    9636273 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp);
    1633         [ +  + ]:    9636273 :     if (h)
    1634                 :            :     {
    1635         [ +  + ]:    3490680 :       if (cmpii(h,lim) > 0) h = subii(h,qp);
    1636                 :    3490680 :       gel(H,i) = h; stable = 0;
    1637                 :            :     }
    1638                 :            :   }
    1639                 :    2450653 :   return stable;
    1640                 :            : }
    1641                 :            : 
    1642                 :            : int
    1643                 :    2245729 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1644                 :            : {
    1645                 :    2245729 :   GEN q = *ptq, qp = muliu(q,p);
    1646                 :    2245729 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1647                 :    2245729 :   *ptq = qp; return stable;
    1648                 :            : }
    1649                 :            : 
    1650                 :            : int
    1651                 :     115754 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1652                 :            : {
    1653                 :     115754 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), lim = shifti(qp,-1);
    1654                 :     115754 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1655                 :     115754 :   long i,j, l = lg(H), m = lgcols(H);
    1656                 :     115754 :   int stable = 1;
    1657         [ +  + ]:    1093000 :   for (j=1; j<l; j++)
    1658         [ +  + ]:   12559905 :     for (i=1; i<m; i++)
    1659                 :            :     {
    1660                 :   11582659 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
    1661         [ +  + ]:   11582659 :       if (h)
    1662                 :            :       {
    1663         [ +  + ]:    4595751 :         if (cmpii(h,lim) > 0) h = subii(h,qp);
    1664                 :    4595751 :         gcoeff(H,i,j) = h; stable = 0;
    1665                 :            :       }
    1666                 :            :     }
    1667                 :     115754 :   *ptq = qp; return stable;
    1668                 :            : }
    1669                 :            : 
    1670                 :            : /* record the degrees of Euclidean remainders (make them as large as
    1671                 :            :  * possible : smaller values correspond to a degenerate sequence) */
    1672                 :            : static void
    1673                 :       1095 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1674                 :            : {
    1675                 :            :   long da,db,dc, ind;
    1676                 :       1095 :   pari_sp av = avma;
    1677                 :            : 
    1678 [ +  - ][ -  + ]:       1095 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1679                 :       1095 :   da = degpol(a);
    1680                 :       1095 :   db = degpol(b);
    1681         [ -  + ]:       1095 :   if (db > da)
    1682                 :          0 :   { swapspec(a,b, da,db); }
    1683         [ -  + ]:       1095 :   else if (!da) return;
    1684                 :       1095 :   ind = 0;
    1685         [ +  + ]:       5815 :   while (db)
    1686                 :            :   {
    1687                 :       4720 :     GEN c = Flx_rem(a,b, p);
    1688                 :       4720 :     a = b; b = c; dc = degpol(c);
    1689         [ -  + ]:       4720 :     if (dc < 0) break;
    1690                 :            : 
    1691                 :       4720 :     ind++;
    1692         [ +  + ]:       4720 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1693         [ -  + ]:       4720 :     if (gc_needed(av,2))
    1694                 :            :     {
    1695         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1696                 :          0 :       gerepileall(av, 2, &a,&b);
    1697                 :            :     }
    1698                 :       4720 :     db = dc; /* = degpol(b) */
    1699                 :            :   }
    1700         [ +  + ]:       1095 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1701                 :       1095 :   avma = av; return;
    1702                 :            : }
    1703                 :            : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1704                 :            :  * "generic" degree sequence as given by dglist, set *Ci and return
    1705                 :            :  * resultant(a,b). Modular version of Collins's subresultant */
    1706                 :            : static ulong
    1707                 :       6293 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1708                 :            : {
    1709                 :            :   long da,db,dc, ind;
    1710                 :       6293 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1711                 :       6293 :   int s = 1;
    1712                 :       6293 :   pari_sp av = avma;
    1713                 :            : 
    1714                 :       6293 :   *C0 = 1; *C1 = 0;
    1715 [ +  - ][ -  + ]:       6293 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1716                 :       6293 :   da = degpol(a);
    1717                 :       6293 :   db = degpol(b);
    1718         [ -  + ]:       6293 :   if (db > da)
    1719                 :            :   {
    1720                 :          0 :     swapspec(a,b, da,db);
    1721         [ #  # ]:          0 :     if (both_odd(da,db)) s = -s;
    1722                 :            :   }
    1723         [ -  + ]:       6293 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1724                 :       6293 :   ind = 0;
    1725         [ +  + ]:      30538 :   while (db)
    1726                 :            :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1727                 :            :      * da = deg a, db = deg b */
    1728                 :      24642 :     GEN c = Flx_rem(a,b, p);
    1729                 :      24642 :     long delta = da - db;
    1730                 :            : 
    1731         [ +  + ]:      24642 :     if (both_odd(da,db)) s = -s;
    1732                 :      24642 :     lb = Fl_mul(b[db+2], cb, p);
    1733                 :      24642 :     a = b; b = c; dc = degpol(c);
    1734                 :      24642 :     ind++;
    1735         [ +  + ]:      24642 :     if (dc != dglist[ind]) { avma = av; return 0; } /* degenerates */
    1736         [ +  + ]:      24245 :     if (g == h)
    1737                 :            :     { /* frequent */
    1738                 :      21726 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1739                 :      21726 :       ca = cb;
    1740                 :      21726 :       cb = cc;
    1741                 :            :     }
    1742                 :            :     else
    1743                 :            :     {
    1744                 :       2519 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1745                 :       2519 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1746                 :       2519 :       ca = cb;
    1747                 :       2519 :       cb = Fl_div(cc, ghdelta, p);
    1748                 :            :     }
    1749                 :      24245 :     da = db; /* = degpol(a) */
    1750                 :      24245 :     db = dc; /* = degpol(b) */
    1751                 :            : 
    1752                 :      24245 :     g = lb;
    1753         [ +  + ]:      24245 :     if (delta == 1)
    1754                 :      14117 :       h = g; /* frequent */
    1755                 :            :     else
    1756                 :      10128 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1757                 :            : 
    1758         [ -  + ]:      24245 :     if (gc_needed(av,2))
    1759                 :            :     {
    1760         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1761                 :          0 :       gerepileall(av, 2, &a,&b);
    1762                 :            :     }
    1763                 :            :   }
    1764         [ -  + ]:       5896 :   if (da > 1) return 0; /* Failure */
    1765                 :            :   /* last non-constant polynomial has degree 1 */
    1766                 :       5896 :   *C0 = Fl_mul(ca, a[2], p);
    1767                 :       5896 :   *C1 = Fl_mul(ca, a[3], p);
    1768                 :       5896 :   res = Fl_mul(cb, b[2], p);
    1769         [ +  + ]:       5896 :   if (s == -1) res = p - res;
    1770                 :       6293 :   avma = av; return res;
    1771                 :            : }
    1772                 :            : 
    1773                 :            : /* u P(X) + v P(-X) */
    1774                 :            : static GEN
    1775                 :          4 : pol_comp(GEN P, GEN u, GEN v)
    1776                 :            : {
    1777                 :          4 :   long i, l = lg(P);
    1778                 :          4 :   GEN y = cgetg(l, t_POL);
    1779         [ +  + ]:         20 :   for (i=2; i<l; i++)
    1780                 :            :   {
    1781                 :         16 :     GEN t = gel(P,i);
    1782         [ +  - ]:         32 :     gel(y,i) = gequal0(t)? gen_0:
    1783                 :          8 :                          (i&1)? gmul(t, gsub(u,v)) /*  odd degree */
    1784         [ +  + ]:         24 :                               : gmul(t, gadd(u,v));/* even degree */
    1785                 :            :   }
    1786                 :          4 :   y[1] = P[1]; return normalizepol_lg(y,l);
    1787                 :            : }
    1788                 :            : 
    1789                 :            : GEN
    1790                 :          0 : polint_triv(GEN xa, GEN ya)
    1791                 :            : {
    1792                 :          0 :   GEN P = NULL, Q = roots_to_pol(xa,0);
    1793                 :          0 :   long i, n = lg(xa);
    1794                 :          0 :   pari_sp av = avma;
    1795         [ #  # ]:          0 :   for (i=1; i<n; i++)
    1796                 :            :   {
    1797                 :            :     GEN T, dP, r;
    1798         [ #  # ]:          0 :     if (gequal0(gel(ya,i))) continue;
    1799                 :          0 :     T = RgX_div_by_X_x(Q, gel(xa,i), NULL);
    1800                 :          0 :     r = poleval(T, gel(xa,i));
    1801 [ #  # ][ #  # ]:          0 :     if (i < n-1 && absi_equal(gel(xa,i), gel(xa,i+1)))
    1802                 :            :     { /* x_i = -x_{i+1} */
    1803                 :          0 :       dP = pol_comp(gdiv(T, r), gel(ya,i), gel(ya,i+1));
    1804                 :          0 :       i++;
    1805                 :            :     }
    1806                 :            :     else
    1807                 :          0 :       dP = gdiv(gmul(gel(ya,i), T), r);
    1808         [ #  # ]:          0 :     P = P? gadd(P, dP): dP;
    1809         [ #  # ]:          0 :     if (gc_needed(av,2))
    1810                 :            :     {
    1811         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"polint_triv2 (i = %ld)",i);
    1812                 :          0 :       P = gerepileupto(av, P);
    1813                 :            :     }
    1814                 :            :   }
    1815         [ #  # ]:          0 :   return P? P: pol_0(0);
    1816                 :            : }
    1817                 :            : 
    1818                 :            : GEN
    1819                 :          2 : FpV_polint(GEN xa, GEN ya, GEN p, long v)
    1820                 :            : {
    1821                 :          2 :   GEN inv,T,dP, P = NULL, Q = FpV_roots_to_pol(xa, p, v);
    1822                 :          2 :   long i, n = lg(xa);
    1823                 :          2 :   pari_sp av = avma;
    1824         [ +  + ]:          6 :   for (i=1; i<n; i++)
    1825                 :            :   {
    1826         [ -  + ]:          4 :     if (!signe(gel(ya,i))) continue;
    1827                 :          4 :     T = FpX_div_by_X_x(Q, gel(xa,i), p, NULL);
    1828                 :          4 :     inv = Fp_inv(FpX_eval(T,gel(xa,i), p), p);
    1829 [ +  - ][ +  - ]:          4 :     if (i < n-1 && equalii(addii(gel(xa,i), gel(xa,i+1)), p))
    1830                 :            :     {
    1831                 :          4 :       dP = pol_comp(T, Fp_mul(gel(ya,i),  inv,p),
    1832                 :          4 :                        Fp_mul(gel(ya,i+1),inv,p));
    1833                 :          4 :       i++; /* x_i = -x_{i+1} */
    1834                 :            :     }
    1835                 :            :     else
    1836                 :          0 :       dP = FpX_Fp_mul(T, Fp_mul(gel(ya,i),inv,p), p);
    1837         [ +  + ]:          4 :     P = P? FpX_add(P, dP, p): dP;
    1838         [ -  + ]:          4 :     if (gc_needed(av,2))
    1839                 :            :     {
    1840         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpV_polint");
    1841                 :          0 :       P = gerepileupto(av, P);
    1842                 :            :     }
    1843                 :            :   }
    1844         [ -  + ]:          2 :   return P? P: pol_0(v);
    1845                 :            : }
    1846                 :            : 
    1847                 :            : static void
    1848                 :        368 : Flv_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p,
    1849                 :            :                GEN *pHp, GEN *pH0p, GEN *pH1p)
    1850                 :            : {
    1851                 :        368 :   GEN T,Q = Flv_roots_to_pol(xa, p, 0);
    1852                 :        368 :   GEN dP  = NULL,  P = NULL;
    1853                 :        368 :   GEN dP0 = NULL, P0= NULL;
    1854                 :        368 :   GEN dP1 = NULL, P1= NULL;
    1855                 :        368 :   long i, n = lg(xa);
    1856                 :            :   ulong inv;
    1857         [ +  + ]:       6264 :   for (i=1; i<n; i++)
    1858                 :            :   {
    1859                 :       5896 :     T = Flx_div_by_X_x(Q, xa[i], p, NULL);
    1860                 :       5896 :     inv = Fl_inv(Flx_eval(T,xa[i], p), p);
    1861                 :            : 
    1862         [ +  - ]:       5896 :     if (ya[i])
    1863                 :            :     {
    1864                 :       5896 :       dP = Flx_Fl_mul(T, Fl_mul(ya[i],inv,p), p);
    1865         [ +  + ]:       5896 :       P = P ? Flx_add(P , dP , p): dP;
    1866                 :            :     }
    1867         [ +  + ]:       5896 :     if (C0[i])
    1868                 :            :     {
    1869                 :       5855 :       dP0= Flx_Fl_mul(T, Fl_mul(C0[i],inv,p), p);
    1870         [ +  + ]:       5855 :       P0= P0? Flx_add(P0, dP0, p): dP0;
    1871                 :            :     }
    1872         [ +  - ]:       5896 :     if (C1[i])
    1873                 :            :     {
    1874                 :       5896 :       dP1= Flx_Fl_mul(T, Fl_mul(C1[i],inv,p), p);
    1875         [ +  + ]:       5896 :       P1= P1? Flx_add(P1, dP1, p): dP1;
    1876                 :            :     }
    1877                 :            :   }
    1878         [ -  + ]:        368 :   *pHp  = (P ? P : zero_Flx(0));
    1879         [ -  + ]:        368 :   *pH0p = (P0? P0: zero_Flx(0));
    1880         [ -  + ]:        368 :   *pH1p = (P1? P1: zero_Flx(0));
    1881                 :        368 : }
    1882                 :            : 
    1883                 :            : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1884                 :            :  * Return 0 in case of degree drop. */
    1885                 :            : static GEN
    1886                 :       7388 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1887                 :            : {
    1888                 :            :   GEN z;
    1889                 :       7388 :   long i, lb = lg(Q);
    1890                 :       7388 :   ulong leadz = Flx_eval(leading_term(Q), x, p);
    1891                 :       7388 :   long vs=mael(Q,2,1);
    1892         [ -  + ]:       7388 :   if (!leadz) return zero_Flx(vs);
    1893                 :            : 
    1894                 :       7388 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1895         [ +  + ]:      61267 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1896                 :       7388 :   z[i] = leadz; return z;
    1897                 :            : }
    1898                 :            : 
    1899                 :            : GEN
    1900                 :       9170 : FpXY_Fq_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1901                 :            : {
    1902                 :       9170 :   pari_sp av = avma;
    1903                 :       9170 :   long i, lb = lg(Q);
    1904                 :            :   GEN z;
    1905         [ +  + ]:       9170 :   if (!T) return FpXY_evaly(Q, y, p, vx);
    1906         [ -  + ]:        730 :   if (lb == 2) return pol_0(vx);
    1907                 :        730 :   z = gel(Q, lb-1);
    1908 [ +  - ][ -  + ]:        730 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
                 [ #  # ]
    1909                 :            : 
    1910         [ +  + ]:        730 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1911         [ +  + ]:      17400 :   for (i=lb-2; i>=2; i--)
    1912                 :            :   {
    1913                 :      16670 :     GEN c = gel(Q,i);
    1914                 :      16670 :     z = FqX_Fq_mul(z, y, T, p);
    1915         [ +  + ]:      16670 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1916                 :            :   }
    1917                 :       9170 :   return gerepileupto(av, z);
    1918                 :            : }
    1919                 :            : 
    1920                 :            : static GEN
    1921                 :       3890 : ZX_norml1(GEN x)
    1922                 :            : {
    1923                 :       3890 :   long i, l = lg(x);
    1924                 :            :   GEN s;
    1925                 :            : 
    1926         [ -  + ]:       3890 :   if (l == 2) return gen_0;
    1927                 :       3890 :   s = gel(x, l-1); /* != 0 */
    1928         [ +  + ]:      21005 :   for (i = l-2; i > 1; i--) {
    1929                 :      17115 :     GEN xi = gel(x,i);
    1930         [ -  + ]:      17115 :     if (!signe(x)) continue;
    1931                 :      17115 :     s = addii_sign(s,1, xi,1);
    1932                 :            :   }
    1933                 :       3890 :   return s;
    1934                 :            : }
    1935                 :            : 
    1936                 :            : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1937                 :            :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1938                 :            :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1939                 :            :  * Return e such that Res(A, B) < 2^e */
    1940                 :            : ulong
    1941                 :     266119 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1942                 :            : {
    1943                 :     266119 :   pari_sp av = avma;
    1944                 :     266119 :   GEN a = gen_0, b = gen_0;
    1945                 :     266119 :   long i , lA = lg(A), lB = lg(B);
    1946                 :            :   double loga, logb;
    1947         [ +  + ]:    1594502 :   for (i=2; i<lA; i++) a = addii(a, sqri(gel(A,i)));
    1948         [ +  + ]:    1014951 :   for (i=2; i<lB; i++)
    1949                 :            :   {
    1950                 :     748832 :     GEN t = gel(B,i);
    1951         [ +  + ]:     748832 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1952                 :     748832 :     b = addii(b, sqri(t));
    1953                 :            :   }
    1954                 :     266119 :   loga = dbllog2(a);
    1955         [ +  + ]:     266119 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1956                 :     266119 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1957         [ +  + ]:     266119 :   avma = av; return (i <= 0)? 1: 1 + (ulong)i;
    1958                 :            : }
    1959                 :            : 
    1960                 :            : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_term(a) [for efficiency] */
    1961                 :            : static ulong
    1962                 :     714701 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1963                 :            : {
    1964                 :     714701 :   GEN ev = FlxY_evalx(b, n, p);
    1965                 :     714701 :   long drop = lg(b) - lg(ev);
    1966                 :     714701 :   ulong r = Flx_resultant(a, ev, p);
    1967 [ +  + ][ -  + ]:     714701 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1968                 :     714701 :   return r;
    1969                 :            : }
    1970                 :            : static GEN
    1971                 :          8 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la)
    1972                 :            : {
    1973                 :          8 :   GEN ev = FpXY_evalx(b, n, p);
    1974                 :          8 :   long drop=lg(b)-lg(ev);
    1975                 :          8 :   GEN r = FpX_resultant(a, ev, p);
    1976 [ -  + ][ #  # ]:          8 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1977                 :          8 :   return r;
    1978                 :            : }
    1979                 :            : 
    1980                 :            : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1981                 :            : /* Return a Fly */
    1982                 :            : static GEN
    1983                 :      17822 : Flx_FlyX_resultant_polint(GEN a, GEN b, ulong p, ulong dres, long sx)
    1984                 :            : {
    1985                 :      17822 :   ulong i, n, la = Flx_lead(a);
    1986                 :      17822 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1987                 :      17822 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1988                 :            :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1989                 :            :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1990         [ +  + ]:     369183 :   for (i=0,n = 1; i < dres; n++)
    1991                 :            :   {
    1992                 :     351361 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1993                 :     351361 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1994                 :            :   }
    1995         [ +  + ]:      17822 :   if (i == dres)
    1996                 :            :   {
    1997                 :      11979 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1998                 :            :   }
    1999                 :      17822 :   return Flv_polint(x,y, p, sx);
    2000                 :            : }
    2001                 :            : 
    2002                 :            : static GEN
    2003                 :       1040 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    2004                 :            : {
    2005                 :       1040 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    2006                 :       1040 :   pari_sp av = avma, av2;
    2007                 :            : 
    2008         [ -  + ]:       1040 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    2009                 :       1040 :   (void)new_chunk(2);
    2010                 :       1040 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    2011                 :       1040 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    2012                 :       1040 :   av2 = avma;
    2013                 :            :   for (;;)
    2014                 :            :   {
    2015                 :       3880 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    2016         [ +  + ]:      13675 :     for (i=1; i<=dy; i++)
    2017                 :       9795 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    2018                 :       9795 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    2019         [ +  + ]:      21510 :     for (   ; i<=dx; i++)
    2020                 :      17630 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    2021 [ +  - ][ +  + ]:       4045 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    2022         [ +  + ]:       3880 :     if (dx < dy) break;
    2023         [ -  + ]:       2840 :     if (gc_needed(av2,1))
    2024                 :            :     {
    2025         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    2026                 :          0 :       gerepilecoeffs(av2,x,dx+1);
    2027                 :            :     }
    2028                 :       2840 :   }
    2029         [ -  + ]:       1040 :   if (dx < 0) return zero_Flx(0);
    2030                 :       1040 :   lx = dx+3; x -= 2;
    2031                 :       1040 :   x[0]=evaltyp(t_POL) | evallg(lx);
    2032                 :       1040 :   x[1]=evalsigne(1) | evalvarn(vx);
    2033                 :       1040 :   x = RgX_recip_shallow(x);
    2034         [ +  + ]:       1040 :   if (dp)
    2035                 :            :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    2036                 :        105 :     GEN t = Flx_powu(gel(y,0), dp, p);
    2037         [ +  + ]:        460 :     for (i=2; i<lx; i++)
    2038                 :        355 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    2039                 :            :   }
    2040                 :       1040 :   return gerepilecopy(av, x);
    2041                 :            : }
    2042                 :            : 
    2043                 :            : /* return a Flx */
    2044                 :            : GEN
    2045                 :        410 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    2046                 :            : {
    2047                 :        410 :   pari_sp av = avma, av2;
    2048                 :            :   long degq,dx,dy,du,dv,dr,signh;
    2049                 :            :   GEN z,g,h,r,p1;
    2050                 :            : 
    2051                 :        410 :   dx=degpol(u); dy=degpol(v); signh=1;
    2052         [ +  + ]:        410 :   if (dx < dy)
    2053                 :            :   {
    2054                 :         15 :     swap(u,v); lswap(dx,dy);
    2055         [ +  + ]:         15 :     if (both_odd(dx, dy)) signh = -signh;
    2056                 :            :   }
    2057         [ -  + ]:        410 :   if (dy < 0) return zero_Flx(sx);
    2058         [ +  + ]:        410 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    2059                 :            : 
    2060                 :        360 :   g = h = pol1_Flx(sx); av2 = avma;
    2061                 :            :   for(;;)
    2062                 :            :   {
    2063                 :       1040 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    2064         [ -  + ]:       1040 :     if (dr == 2) { avma = av; return zero_Flx(sx); }
    2065                 :       1040 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    2066                 :       1040 :     u = v; p1 = g; g = leading_term(u);
    2067      [ -  +  + ]:       1040 :     switch(degq)
    2068                 :            :     {
    2069                 :          0 :       case 0: break;
    2070                 :            :       case 1:
    2071                 :        775 :         p1 = Flx_mul(h,p1, p); h = g; break;
    2072                 :            :       default:
    2073                 :        265 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    2074                 :        265 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    2075                 :            :     }
    2076         [ +  + ]:       1040 :     if (both_odd(du,dv)) signh = -signh;
    2077                 :       1040 :     v = FlxY_Flx_div(r, p1, p);
    2078         [ +  + ]:       1040 :     if (dr==3) break;
    2079         [ -  + ]:        680 :     if (gc_needed(av2,1))
    2080                 :            :     {
    2081         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, dr = %ld",dr);
    2082                 :          0 :       gerepileall(av2,4, &u, &v, &g, &h);
    2083                 :            :     }
    2084                 :        680 :   }
    2085                 :        360 :   z = gel(v,2);
    2086         [ +  + ]:        360 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    2087         [ +  + ]:        360 :   if (signh < 0) z = Flx_neg(z,p);
    2088                 :        410 :   return gerepileupto(av, z);
    2089                 :            : }
    2090                 :            : 
    2091                 :            : /* Warning:
    2092                 :            :  * This function switches between valid and invalid variable ordering*/
    2093                 :            : 
    2094                 :            : static GEN
    2095                 :       2538 : FlxY_to_FlyX(GEN b, long sv)
    2096                 :            : {
    2097                 :       2538 :   long i, n=-1;
    2098                 :       2538 :   long sw = b[1]&VARNBITS;
    2099         [ +  + ]:       8399 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    2100                 :       2538 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    2101                 :            : }
    2102                 :            : 
    2103                 :            : /* Return a Fly*/
    2104                 :            : GEN
    2105                 :       2538 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    2106                 :            : {
    2107                 :       2538 :   pari_sp ltop=avma;
    2108                 :       2538 :   long dres = degpol(a)*degpol(b);
    2109                 :       2538 :   long sx=a[1], sy=b[1]&VARNBITS;
    2110                 :            :   GEN z;
    2111                 :       2538 :   b = FlxY_to_FlyX(b,sx);
    2112         [ +  + ]:       2538 :   if ((ulong)dres >= pp)
    2113                 :        410 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    2114                 :            :   else
    2115                 :       2128 :     z = Flx_FlyX_resultant_polint(a, b, pp, (ulong)dres, sy);
    2116                 :       2538 :   return gerepileupto(ltop,z);
    2117                 :            : }
    2118                 :            : 
    2119                 :            : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    2120                 :            :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    2121                 :            :  * We could return a vector of coeffs, but it is convenient to have degpol()
    2122                 :            :  * and friends available. Even in that case, it will behave nicely with all
    2123                 :            :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    2124                 :            :  * FOR INTERNAL USE! */
    2125                 :            : GEN
    2126                 :       2442 : swap_vars(GEN b0, long v)
    2127                 :            : {
    2128                 :       2442 :   long i, n = RgX_degree(b0, v);
    2129                 :            :   GEN b, x;
    2130         [ -  + ]:       2442 :   if (n < 0) return pol_0(v);
    2131                 :       2442 :   b = cgetg(n+3, t_POL); x = b + 2;
    2132                 :       2442 :   b[1] = evalsigne(1) | evalvarn(v);
    2133         [ +  + ]:      27531 :   for (i=0; i<=n; i++) gel(x,i) = polcoeff_i(b0, i, v);
    2134                 :       2442 :   return b;
    2135                 :            : }
    2136                 :            : 
    2137                 :            : /* assume varn(b) << varn(a) */
    2138                 :            : /* return a FpY*/
    2139                 :            : GEN
    2140                 :       2397 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    2141                 :            : {
    2142                 :       2397 :   long i,n,dres, vX = varn(b), vY = varn(a);
    2143                 :            :   GEN la,x,y;
    2144                 :            : 
    2145         [ +  + ]:       2397 :   if (lgefint(p) == 3)
    2146                 :            :   {
    2147                 :       2395 :     ulong pp = uel(p,2);
    2148                 :       2395 :     b = ZXX_to_FlxX(b, pp, vY);
    2149                 :       2395 :     a = ZX_to_Flx(a, pp);
    2150                 :       2395 :     x = Flx_FlxY_resultant(a, b, pp);
    2151                 :       2395 :     return Flx_to_ZX(x);
    2152                 :            :   }
    2153                 :          2 :   dres = degpol(a)*degpol(b);
    2154                 :          2 :   b = swap_vars(b, vY);
    2155                 :          2 :   la = leading_term(a);
    2156                 :          2 :   x = cgetg(dres+2, t_VEC);
    2157                 :          2 :   y = cgetg(dres+2, t_VEC);
    2158                 :            :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    2159                 :            :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    2160         [ +  + ]:          6 :   for (i=0,n = 1; i < dres; n++)
    2161                 :            :   {
    2162                 :          4 :     gel(x,++i) = utoipos(n);
    2163                 :          4 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la);
    2164                 :          4 :     gel(x,++i) = subis(p,n);
    2165                 :          4 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la);
    2166                 :            :   }
    2167         [ -  + ]:          2 :   if (i == dres)
    2168                 :            :   {
    2169                 :          0 :     gel(x,++i) = gen_0;
    2170                 :          0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la);
    2171                 :            :   }
    2172                 :       2397 :   return FpV_polint(x,y, p, vX);
    2173                 :            : }
    2174                 :            : 
    2175                 :            : GEN
    2176                 :        285 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    2177                 :            : {
    2178                 :        285 :   GEN x = deg1pol_shallow(gen_1, pol_x(varn(a)), fetch_var_higher()); /* x+y */
    2179                 :        285 :   x = FpX_FpXY_resultant(a, poleval(b,x),p);
    2180                 :        285 :   (void)delete_var(); return x;
    2181                 :            : }
    2182                 :            : 
    2183                 :            : /* 0, 1, -1, 2, -2, ... */
    2184                 :            : #define next_lambda(a) (a>0 ? -a : 1-a)
    2185                 :            : GEN
    2186                 :          0 : FpX_compositum(GEN a, GEN b, GEN p)
    2187                 :            : {
    2188                 :          0 :   long k, v = fetch_var_higher();
    2189         [ #  # ]:          0 :   for (k = 1;; k = next_lambda(k))
    2190                 :            :   {
    2191                 :          0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    2192                 :          0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    2193         [ #  # ]:          0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    2194                 :          0 :   }
    2195                 :            : }
    2196                 :            : 
    2197                 :            : #if 0
    2198                 :            : /* Return x such that theta(x) - theta(27448) >= ln(2)*bound */
    2199                 :            : /* NB: theta(27449) ~ 27225.387, theta(x) > 0.98 x for x>7481
    2200                 :            :  * (Schoenfeld, 1976 for x > 1155901 + direct calculations) */
    2201                 :            : static ulong
    2202                 :            : get_theta_x(ulong bound)
    2203                 :            : { return (ulong)ceil((bound * LOG2 + 27225.388) / 0.98); }
    2204                 :            : #endif
    2205                 :            : /* 27449 = prime(3000) */
    2206                 :            : void
    2207                 :    2958758 : init_modular(forprime_t *S) { u_forprime_init(S, 27449, ULONG_MAX); }
    2208                 :            : 
    2209                 :            : /* Assume A in Z[Y], B in Q[Y][X], and Res_Y(A, B) in Z[X].
    2210                 :            :  * If lambda = NULL, return Res_Y(A,B).
    2211                 :            :  * Otherwise, find a small lambda (start from *lambda, use the sequence above)
    2212                 :            :  * such that R(X) = Res_Y(A(Y), B(X + lambda Y)) is squarefree, reset *lambda
    2213                 :            :  * to the chosen value and return R
    2214                 :            :  *
    2215                 :            :  * If LERS is non-NULL, set it to the Last non-constant polynomial in the
    2216                 :            :  * Euclidean Remainder Sequence */
    2217                 :            : GEN
    2218                 :       2470 : ZX_ZXY_resultant_all(GEN A, GEN B0, long *plambda, GEN *LERS)
    2219                 :            : {
    2220                 :       2470 :   int checksqfree = plambda? 1: 0, stable;
    2221         [ +  + ]:       2470 :   long lambda = plambda? *plambda: 0, cnt = 0;
    2222                 :            :   ulong bound, p, dp;
    2223                 :       2470 :   pari_sp av = avma, av2 = 0;
    2224                 :       2470 :   long i,n, lb, degA = degpol(A), dres = degA*degpol(B0);
    2225                 :       2470 :   long v = fetch_var_higher();
    2226                 :       2470 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    2227                 :       2470 :   long sX = evalvarn(vX);
    2228                 :            :   GEN x, y, dglist, dB, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    2229                 :            :   forprime_t S;
    2230                 :            : 
    2231                 :       2470 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    2232         [ +  + ]:       2470 :   if (LERS)
    2233                 :            :   {
    2234         [ -  + ]:        320 :     if (!checksqfree)
    2235                 :          0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    2236                 :        320 :     C0 = cgetg(dres+2, t_VECSMALL);
    2237                 :        320 :     C1 = cgetg(dres+2, t_VECSMALL);
    2238                 :        320 :     dglist = cgetg(dres+1, t_VECSMALL);
    2239                 :            :   }
    2240                 :       2470 :   x = cgetg(dres+2, t_VECSMALL);
    2241                 :       2470 :   y = cgetg(dres+2, t_VECSMALL);
    2242                 :       2470 :   B0 = Q_remove_denom(B0, &dB);
    2243         [ +  - ]:       2470 :   if (!dB) B0 = leafcopy(B0);
    2244                 :       2470 :   A = leafcopy(A);
    2245                 :       2470 :   B = B0;
    2246                 :       2470 :   setvarn(A,v);
    2247         [ +  + ]:       2470 :   if (vX == vY) setvarn(B0, v);
    2248                 :            : 
    2249                 :            :   /* make sure p large enough */
    2250                 :       2470 :   u_forprime_init(&S, maxuu(dres << 1, 27499), ULONG_MAX);
    2251                 :            : INIT:
    2252                 :            :   /* allways except the first time */
    2253 [ +  + ][ +  + ]:       2875 :   if (av2) { avma = av2; lambda = next_lambda(lambda); }
    2254         [ +  + ]:       2875 :   if (vX == vY)
    2255                 :            :   {
    2256                 :        460 :     B = RgX_translate(B0, pol_x(vX));
    2257         [ +  - ]:        460 :     if (lambda) B = RgX_unscale(B, stoi(lambda));
    2258                 :            :   }
    2259                 :            :   else
    2260                 :            :   {
    2261         [ +  + ]:       2415 :     if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    2262                 :       2415 :     B = swap_vars(B, vY); setvarn(B,v);
    2263                 :            :   }
    2264                 :            :   /* B0(lambda v + x, v) */
    2265 [ -  + ][ #  # ]:       2875 :   if (DEBUGLEVEL>4 && checksqfree) err_printf("Trying lambda = %ld\n", lambda);
    2266                 :       2875 :   av2 = avma;
    2267                 :            : 
    2268         [ +  + ]:       2875 :   if (degA <= 3)
    2269                 :            :   { /* sub-resultant faster for small degrees */
    2270         [ +  + ]:       1555 :     if (LERS)
    2271                 :            :     { /* implies checksqfree */
    2272                 :        300 :       H = resultant_all(A,B,&q);
    2273 [ +  + ][ -  + ]:        300 :       if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    2274                 :        230 :       H0 = gel(q,2);
    2275         [ +  + ]:        230 :       if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    2276                 :        230 :       H1 = gel(q,3);
    2277         [ +  + ]:        230 :       if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    2278                 :            :     }
    2279                 :            :     else
    2280                 :       1255 :       H = resultant(A,B);
    2281 [ +  - ][ +  + ]:       1485 :     if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
    2282                 :       1375 :     goto END;
    2283                 :            :   }
    2284                 :            : 
    2285                 :       1320 :   H = H0 = H1 = NULL;
    2286                 :       1320 :   lb = lg(B);
    2287                 :       1320 :   bound = ZX_ZXY_ResBound(A, B, dB);
    2288         [ -  + ]:       1320 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2289                 :       1320 :   dp = 1;
    2290         [ +  - ]:      16062 :   while ((p = u_forprime_next(&S)))
    2291                 :            :   {
    2292 [ -  + ][ #  # ]:      16062 :     if (dB) { dp = smodis(dB, p); if (!dp) continue; }
    2293                 :            : 
    2294                 :      16062 :     a = ZX_to_Flx(A, p);
    2295                 :      16062 :     b = ZXX_to_FlxX(B, p, varn(A));
    2296         [ +  + ]:      16062 :     if (LERS)
    2297                 :            :     {
    2298 [ +  - ][ -  + ]:        368 :       if (degpol(a) < degA || lg(b) < lb) continue; /* p | lc(A)lc(B) */
    2299         [ +  + ]:        368 :       if (checksqfree)
    2300                 :            :       { /* find degree list for generic Euclidean Remainder Sequence */
    2301                 :        120 :         long goal = minss(degpol(a), degpol(b)); /* longest possible */
    2302         [ +  + ]:        855 :         for (n=1; n <= goal; n++) dglist[n] = 0;
    2303                 :        120 :         setlg(dglist, 1);
    2304         [ +  + ]:       1155 :         for (n=0; n <= dres; n++)
    2305                 :            :         {
    2306                 :       1095 :           ev = FlxY_evalx_drop(b, n, p);
    2307                 :       1095 :           Flx_resultant_set_dglist(a, ev, dglist, p);
    2308         [ +  + ]:       1095 :           if (lg(dglist)-1 == goal) break;
    2309                 :            :         }
    2310                 :            :         /* last pol in ERS has degree > 1 ? */
    2311                 :        120 :         goal = lg(dglist)-1;
    2312 [ +  + ][ -  + ]:        120 :         if (degpol(B) == 1) { if (!goal) goto INIT; }
    2313                 :            :         else
    2314                 :            :         {
    2315         [ -  + ]:        115 :           if (goal <= 1) goto INIT;
    2316 [ -  + ][ -  + ]:        115 :           if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    2317                 :            :         }
    2318         [ -  + ]:        120 :         if (DEBUGLEVEL>4)
    2319                 :          0 :           err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    2320                 :            :       }
    2321                 :            : 
    2322         [ +  + ]:       6661 :       for (i=0,n = 0; i <= dres; n++)
    2323                 :            :       {
    2324                 :       6293 :         ev = FlxY_evalx_drop(b, n, p);
    2325                 :       6293 :         x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    2326         [ +  + ]:       6293 :         if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    2327                 :            :       }
    2328                 :        368 :       Flv_polint_all(x,y,C0,C1, p, &Hp, &H0p, &H1p);
    2329                 :            :     }
    2330                 :            :     else
    2331                 :            :     {
    2332                 :      15694 :       long dropa = degA - degpol(a), dropb = lb - lg(b);
    2333                 :      15694 :       Hp = Flx_FlyX_resultant_polint(a, b, p, (ulong)dres, sX);
    2334 [ -  + ][ #  # ]:      15694 :       if (dropa && dropb)
    2335                 :          0 :         Hp = zero_Flx(sX);
    2336                 :            :       else {
    2337         [ -  + ]:      15694 :         if (dropa)
    2338                 :            :         { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2339                 :          0 :           GEN c = gel(b,lb-1); /* lc(B) */
    2340         [ #  # ]:          0 :           if (!odd(lb)) c = Flx_neg(c, p); /* deg B = lb - 3 */
    2341         [ #  # ]:          0 :           if (!Flx_equal1(c)) {
    2342                 :          0 :             c = Flx_powu(c, dropa, p);
    2343         [ #  # ]:          0 :             if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    2344                 :            :           }
    2345                 :            :         }
    2346         [ -  + ]:      15694 :         else if (dropb)
    2347                 :            :         { /* multiply by lc(A)^(deg B - deg b) */
    2348                 :          0 :           ulong c = a[degA+2]; /* lc(A) */
    2349                 :          0 :           c = Fl_powu(c, dropb, p);
    2350         [ #  # ]:          0 :           if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    2351                 :            :         }
    2352                 :            :       }
    2353                 :            :     }
    2354 [ +  + ][ -  + ]:      16062 :     if (!H && degpol(Hp) != dres) continue;
    2355         [ -  + ]:      16062 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2356         [ +  + ]:      16062 :     if (checksqfree) {
    2357         [ +  + ]:        635 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    2358         [ -  + ]:        410 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2359                 :        410 :       checksqfree = 0;
    2360                 :            :     }
    2361                 :            : 
    2362         [ +  + ]:      15837 :     if (!H)
    2363                 :            :     { /* initialize */
    2364                 :       1095 :       q = utoipos(p); stable = 0;
    2365                 :       1095 :       H = ZX_init_CRT(Hp, p,vX);
    2366         [ +  + ]:       1095 :       if (LERS) {
    2367                 :        120 :         H0= ZX_init_CRT(H0p, p,vX);
    2368                 :        120 :         H1= ZX_init_CRT(H1p, p,vX);
    2369                 :            :       }
    2370                 :            :     }
    2371                 :            :     else
    2372                 :            :     {
    2373         [ +  + ]:      14742 :       if (LERS) {
    2374                 :        248 :         GEN qp = muliu(q,p);
    2375                 :        248 :         stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    2376                 :        248 :                 & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    2377                 :        248 :                 & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    2378                 :        248 :         q = qp;
    2379                 :            :       }
    2380                 :            :       else
    2381                 :      14494 :         stable = ZX_incremental_CRT(&H, Hp, &q, p);
    2382                 :            :     }
    2383                 :            :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    2384                 :            :      * Probabilistic anyway for H0, H1 */
    2385 [ -  + ][ #  # ]:      15837 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
                 [ #  # ]
    2386         [ #  # ]:          0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    2387 [ +  + ][ +  + ]:      15837 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    2388         [ +  + ]:      14742 :     if (gc_needed(av,2))
    2389                 :            :     {
    2390         [ -  + ]:         10 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    2391         [ -  + ]:         10 :       gerepileall(av2, LERS? 4: 2, &H, &q, &H0, &H1);
    2392                 :            :     }
    2393                 :            :   }
    2394         [ -  + ]:       1095 :   if (!p) pari_err_OVERFLOW("ZX_ZXY_rnfequation [ran out of primes]");
    2395                 :            : END:
    2396         [ -  + ]:       2470 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    2397                 :       2470 :   setvarn(H, vX); (void)delete_var();
    2398         [ +  + ]:       2470 :   if (plambda) *plambda = lambda;
    2399         [ +  + ]:       2470 :   if (LERS)
    2400                 :            :   {
    2401                 :        320 :     *LERS = mkvec2(H0,H1);
    2402                 :        320 :     gerepileall(av, 2, &H, LERS);
    2403                 :        320 :     return H;
    2404                 :            :   }
    2405                 :       2470 :   return gerepilecopy(av, H);
    2406                 :            : }
    2407                 :            : 
    2408                 :            : GEN
    2409                 :       1775 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    2410                 :            : {
    2411                 :       1775 :   return ZX_ZXY_resultant_all(A, B, lambda, NULL);
    2412                 :            : }
    2413                 :            : 
    2414                 :            : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    2415                 :            :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    2416                 :            :  * squarefree */
    2417                 :            : GEN
    2418                 :       1640 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    2419                 :            : {
    2420                 :       1640 :   pari_sp av = avma;
    2421                 :            :   GEN R, a;
    2422                 :            :   long dA;
    2423                 :            :   int delvar;
    2424                 :            : 
    2425         [ -  + ]:       1640 :   if (v < 0) v = 0;
    2426         [ +  - ]:       1640 :   switch (typ(A))
    2427                 :            :   {
    2428         [ +  - ]:       1640 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    2429         [ #  # ]:          0 :       A = dA? gel(A,2): gen_0; /* fall through */
    2430                 :            :     default:
    2431         [ #  # ]:          0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    2432                 :          0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    2433                 :            :   }
    2434                 :       1640 :   delvar = 0;
    2435         [ +  + ]:       1640 :   if (varn(T) == 0)
    2436                 :            :   {
    2437                 :       1320 :     long v0 = fetch_var(); delvar = 1;
    2438                 :       1320 :     T = leafcopy(T); setvarn(T,v0);
    2439                 :       1320 :     A = leafcopy(A); setvarn(A,v0);
    2440                 :            :   }
    2441                 :       1640 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    2442         [ +  + ]:       1640 :   if (delvar) (void)delete_var();
    2443                 :       1640 :   setvarn(R, v); a = leading_term(T);
    2444         [ +  + ]:       1640 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    2445                 :       1640 :   return gerepileupto(av, R);
    2446                 :            : }
    2447                 :            : 
    2448                 :            : 
    2449                 :            : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    2450                 :            : GEN
    2451                 :       7105 : ZXQ_charpoly(GEN A, GEN T, long v)
    2452                 :            : {
    2453         [ +  + ]:       7105 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    2454                 :            : }
    2455                 :            : 
    2456                 :            : static GEN
    2457                 :     540528 : trivial_case(GEN A, GEN B)
    2458                 :            : {
    2459                 :            :   long d;
    2460         [ +  + ]:     540528 :   if (typ(A) == t_INT) return powiu(A, degpol(B));
    2461                 :     537182 :   d = degpol(A);
    2462         [ +  + ]:     537182 :   if (d == 0) return trivial_case(gel(A,2),B);
    2463         [ +  + ]:     534696 :   if (d < 0) return gen_0;
    2464                 :     540528 :   return NULL;
    2465                 :            : }
    2466                 :            : 
    2467                 :            : /* floating point resultant */
    2468                 :            : static GEN
    2469                 :         21 : fp_resultant(GEN a, GEN b)
    2470                 :            : {
    2471                 :            :   long da, db, dc;
    2472                 :         21 :   GEN res = gen_1;
    2473                 :            :   pari_sp av;
    2474                 :            : 
    2475 [ +  - ][ -  + ]:         21 :   if (lgpol(a)==0 || lgpol(b)==0) return gen_0;
    2476                 :         21 :   da = degpol(a);
    2477                 :         21 :   db = degpol(b);
    2478         [ +  + ]:         21 :   if (db > da)
    2479                 :            :   {
    2480                 :          6 :     swapspec(a,b, da,db);
    2481         [ -  + ]:          6 :     if (both_odd(da,db)) res = gneg(res);
    2482                 :            :   }
    2483         [ -  + ]:         15 :   else if (!da) return gen_1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
    2484                 :         21 :   av = avma;
    2485         [ +  + ]:       2034 :   while (db)
    2486                 :            :   {
    2487                 :       2014 :     GEN lb = gel(b,db+2), c = RgX_rem(a,b);
    2488                 :       2014 :     c = normalizepol_approx(c, lg(c)); /* kill leading zeroes without warning */
    2489                 :       2014 :     a = b; b = c; dc = degpol(c);
    2490         [ +  + ]:       2014 :     if (dc < 0) { avma = av; return gen_0; }
    2491                 :            : 
    2492         [ +  + ]:       2013 :     if (both_odd(da,db)) res = gneg(res);
    2493                 :       2013 :     res = gmul(res, gpowgs(lb, da - dc));
    2494         [ -  + ]:       2013 :     if (gc_needed(av,1)) {
    2495         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"fp_resultant");
    2496                 :          0 :       gerepileall(av, 3, &a,&b,&res);
    2497                 :            :     }
    2498                 :       2013 :     da = db; /* = degpol(a) */
    2499                 :       2013 :     db = dc; /* = degpol(b) */
    2500                 :            :   }
    2501                 :         21 :   return gerepileupto(av, gmul(res, gpowgs(gel(b,2), da)));
    2502                 :            : }
    2503                 :            : 
    2504                 :            : #ifdef LONG_IS_64BIT
    2505                 :            : #define  ZXRES_PRIME  4611686018427388039UL
    2506                 :            : #else
    2507                 :            : #define  ZXRES_PRIME  27449UL
    2508                 :            : #endif
    2509                 :            : 
    2510                 :            : static long
    2511                 :     266109 : get_nbprimes(ulong bound, GEN dB, ulong *pt_start)
    2512                 :            : {
    2513                 :     266109 :   pari_sp av = avma;
    2514                 :            :   ulong p;
    2515                 :     266109 :   long i = 0, e;
    2516                 :     266109 :   ulong pstart = ZXRES_PRIME;
    2517                 :     266109 :   ulong s = 0;
    2518                 :            :   forprime_t S;
    2519                 :     266109 :   *pt_start = pstart;
    2520 [ +  + ][ +  + ]:     266109 :   if ((ulong)expu(pstart) > bound && (!dB || umodiu(dB, pstart)!=0)) return 1;
                 [ +  - ]
    2521                 :      81506 :   u_forprime_init(&S, pstart, ULONG_MAX);
    2522         [ +  - ]:     502192 :   while ((p = u_forprime_next(&S)))
    2523                 :            :   {
    2524                 :     502192 :     i++;
    2525 [ +  + ][ -  + ]:     502192 :     if (dB && umodiu(dB, p)==0) continue;
    2526                 :     502192 :     e = expu(p);
    2527                 :     502192 :     s += e;
    2528         [ +  + ]:     502192 :     if (s > bound) break;
    2529                 :            :   }
    2530                 :     266109 :   avma = av; return i;
    2531                 :            : }
    2532                 :            : 
    2533                 :            : static ulong
    2534                 :     686795 : ZX_resultant_prime(GEN a, GEN b, ulong dp, long degA, long degB, ulong p)
    2535                 :            : {
    2536                 :     686795 :   pari_sp av = avma;
    2537                 :            :   ulong H;
    2538                 :            :   long dropa, dropb;
    2539         [ +  + ]:     686795 :   if (!b) b = Flx_deriv(a, p);
    2540                 :     686795 :   dropa = degA - degpol(a);
    2541                 :     686795 :   dropb = degB - degpol(b);
    2542 [ +  + ][ -  + ]:     686795 :   if (dropa && dropb) /* p | lc(A), p | lc(B) */
    2543                 :          0 :   { avma = av; return 0; }
    2544                 :     686795 :   H = Flx_resultant(a, b, p);
    2545         [ +  + ]:     686795 :   if (dropa)
    2546                 :            :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2547                 :          7 :     ulong c = b[degB+2]; /* lc(B) */
    2548         [ -  + ]:          7 :     if (odd(degB)) c = p - c;
    2549                 :          7 :     c = Fl_powu(c, dropa, p);
    2550         [ +  + ]:          7 :     if (c != 1) H = Fl_mul(H, c, p);
    2551                 :            :   }
    2552         [ -  + ]:     686788 :   else if (dropb)
    2553                 :            :   { /* multiply by lc(A)^(deg B - deg b) */
    2554                 :          0 :     ulong c = a[degA+2]; /* lc(A) */
    2555                 :          0 :     c = Fl_powu(c, dropb, p);
    2556         [ #  # ]:          0 :     if (c != 1) H = Fl_mul(H, c, p);
    2557                 :            :   }
    2558         [ +  + ]:     686795 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2559                 :     686795 :   avma = av; return H;
    2560                 :            : }
    2561                 :            : 
    2562                 :            : /* If B=NULL, assume B=A' */
    2563                 :            : static GEN
    2564                 :     514254 : ZX_resultant_slice(GEN A, GEN B, GEN dB, ulong p, long n, ulong *plast, GEN *mod)
    2565                 :            : {
    2566                 :     514254 :   ulong dp = 1;
    2567                 :     514254 :   pari_sp av = avma;
    2568                 :            :   long degA, degB, i;
    2569                 :            :   GEN H, P, T;
    2570                 :            :   forprime_t S;
    2571                 :            : 
    2572                 :     514254 :   degA = degpol(A);
    2573         [ +  + ]:     514254 :   degB = B ? degpol(B): degA - 1;
    2574         [ +  + ]:     514254 :   if (n == 1)
    2575                 :            :   {
    2576                 :            :     ulong Hp;
    2577                 :            :     GEN a, b;
    2578                 :     467208 :     p = unextprime(p); *plast = p+1;
    2579         [ +  + ]:     467208 :     a = ZX_to_Flx(A, p), b = B ? ZX_to_Flx(B, p): NULL;
    2580         [ +  + ]:     467208 :     if (dB) dp = umodiu(dB, p);
    2581                 :     467208 :     Hp = ZX_resultant_prime(a, b, dp, degA, degB, p);
    2582                 :     467208 :     avma = av;
    2583                 :     467208 :     *mod = utoi(p); return utoi(Hp);
    2584                 :            :   }
    2585                 :      47046 :   u_forprime_init(&S, p, ULONG_MAX);
    2586                 :      47046 :   P = cgetg(n+1, t_VECSMALL);
    2587         [ +  + ]:     266633 :   for (i=1; i <= n; i++)
    2588                 :     219587 :     P[i] = u_forprime_next(&S);
    2589                 :      47046 :   T = ZV_producttree(P);
    2590                 :      47046 :   A = ZX_nv_mod_tree(A, P, T);
    2591         [ +  + ]:      47046 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    2592                 :      47046 :   H = cgetg(n+1, t_VECSMALL);
    2593         [ +  + ]:     266633 :   for(i=1; i <= n; i++)
    2594                 :            :   {
    2595                 :     219587 :     ulong p = P[i];
    2596                 :            :     GEN a, b;
    2597 [ +  + ][ -  + ]:     219587 :     if (dB) { dp = umodiu(dB, p); if (!dp) continue; }
    2598         [ +  + ]:     219587 :     a = gel(A, i); b = B ? gel(B, i): NULL;
    2599                 :     219587 :     H[i] = ZX_resultant_prime(a, b, dp, degA, degB, p);
    2600                 :            :   }
    2601                 :      47046 :   H = ZV_chinese_tree(H, P, T, mod); *plast=P[n]+1;
    2602                 :      47046 :   gerepileall(av, 2, &H, mod);
    2603                 :     514254 :   return H;
    2604                 :            : }
    2605                 :            : 
    2606                 :            : /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
    2607                 :            : /* if B=NULL, take B = A' */
    2608                 :            : GEN
    2609                 :     269466 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    2610                 :            : {
    2611                 :            :   ulong p;
    2612                 :     269466 :   pari_sp av = avma;
    2613                 :            :   long degA, n, m;
    2614                 :            :   GEN  H, P, mod;
    2615                 :     269466 :   int is_disc = !B;
    2616         [ +  + ]:     269466 :   if (is_disc) B = ZX_deriv(A);
    2617                 :            : 
    2618 [ +  + ][ +  + ]:     269466 :   if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
    2619                 :     266109 :   degA = degpol(A);
    2620         [ +  + ]:     266109 :   if (!bound)
    2621                 :            :   {
    2622                 :     264799 :     bound = ZX_ZXY_ResBound(A, B, dB);
    2623         [ +  + ]:     264799 :     if (bound > 10000)
    2624                 :            :     {
    2625                 :         20 :       const long CNTMAX = 5; /* to avoid oo loops if R = 0 */
    2626                 :         20 :       long bnd = 0, cnt;
    2627                 :         20 :       long prec = nbits2prec( maxss(gexpo(A), gexpo(B)) + 1 );
    2628         [ -  + ]:         20 :       long bndden = dB? (long)(dbllog2(dB)*degA): 0;
    2629         [ +  - ]:         21 :       for(cnt = 1; cnt < CNTMAX; cnt++, prec = precdbl(prec))
    2630                 :            :       {
    2631                 :         21 :         GEN R = fp_resultant(RgX_gtofp(A, prec), RgX_gtofp(B, prec));
    2632                 :         21 :         bnd = gexpo(R) - bndden + 1;
    2633 [ +  + ][ +  - ]:         21 :         if (bnd >= 0 && bnd <= (long)bound && !gequal0(R))
                 [ +  - ]
    2634                 :            :         {
    2635                 :         20 :           bound = bnd; break;
    2636                 :            :         }
    2637                 :            :       }
    2638                 :            :     }
    2639                 :            :   }
    2640                 :     266109 :   n = get_nbprimes(bound+1, dB, &p);/* +1 to account for sign */
    2641         [ +  + ]:     266109 :   if (is_disc)
    2642                 :      16047 :     B = NULL;
    2643         [ +  + ]:     266109 :   m = minss(degpol(A)+(B ? degpol(B): 0), n);
    2644         [ +  + ]:     266109 :   if (m == 1)
    2645                 :     184603 :     H = ZX_resultant_slice(A, B, dB, p, n, &p, &mod);
    2646                 :            :   else
    2647                 :            :   {
    2648                 :      81506 :     long i, s = n/m, r = n - m*s;
    2649         [ -  + ]:      81506 :     if (DEBUGLEVEL > 4)
    2650                 :          0 :       err_printf("ZX_resultant: bound 2^%ld, nb primes: %ld\n",bound, n);
    2651                 :      81506 :     H = cgetg(m+1+!!r, t_VEC); P = cgetg(m+1+!!r, t_VEC);
    2652         [ +  + ]:     398646 :     for (i=1; i<=m; i++)
    2653                 :            :     {
    2654                 :     317140 :       gel(H, i) = ZX_resultant_slice(A, B, dB, p, s, &p, &gel(P, i));
    2655         [ -  + ]:     317140 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",100*i/m);
    2656                 :            :     }
    2657         [ +  + ]:      81506 :     if (r)
    2658                 :      12511 :       gel(H, i) = ZX_resultant_slice(A, B, dB, p, r, &p, &gel(P, i));
    2659                 :      81506 :     H = ZV_chinese(H, P, &mod);
    2660         [ -  + ]:      81506 :     if (DEBUGLEVEL>5) err_printf("done\n");
    2661                 :            :   }
    2662                 :     266109 :   H = Fp_center(H, mod, shifti(mod,-1));
    2663                 :     269466 :   return gerepileuptoint(av, H);
    2664                 :            : }
    2665                 :            : 
    2666                 :            : /* A0 and B0 in Q[X] */
    2667                 :            : GEN
    2668                 :       7122 : QX_resultant(GEN A0, GEN B0)
    2669                 :            : {
    2670                 :            :   GEN s, a, b, A, B;
    2671                 :       7122 :   pari_sp av = avma;
    2672                 :            : 
    2673                 :       7122 :   A = Q_primitive_part(A0, &a);
    2674                 :       7122 :   B = Q_primitive_part(B0, &b);
    2675                 :       7122 :   s = ZX_resultant(A, B);
    2676         [ -  + ]:       7122 :   if (!signe(s)) { avma = av; return gen_0; }
    2677         [ -  + ]:       7122 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    2678         [ +  - ]:       7122 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    2679                 :       7122 :   return gerepileupto(av, s);
    2680                 :            : }
    2681                 :            : 
    2682                 :            : GEN
    2683                 :      25060 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    2684                 :            : 
    2685                 :            : GEN
    2686                 :          0 : QXQ_intnorm(GEN A, GEN B)
    2687                 :            : {
    2688                 :            :   GEN c, n, R, lB;
    2689                 :          0 :   long dA = degpol(A), dB = degpol(B);
    2690                 :          0 :   pari_sp av = avma;
    2691         [ #  # ]:          0 :   if (dA < 0) return gen_0;
    2692                 :          0 :   A = Q_primitive_part(A, &c);
    2693 [ #  # ][ #  # ]:          0 :   if (!c || typ(c) == t_INT) {
    2694                 :          0 :     n = c;
    2695                 :          0 :     R = ZX_resultant(B, A);
    2696                 :            :   } else {
    2697                 :          0 :     n = gel(c,1);
    2698                 :          0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2699                 :            :   }
    2700 [ #  # ][ #  # ]:          0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2701                 :          0 :   lB = leading_term(B);
    2702         [ #  # ]:          0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2703                 :          0 :   return gerepileuptoint(av, R);
    2704                 :            : }
    2705                 :            : 
    2706                 :            : GEN
    2707                 :          0 : QXQ_norm(GEN A, GEN B)
    2708                 :            : {
    2709                 :            :   GEN c, R, lB;
    2710                 :          0 :   long dA = degpol(A), dB = degpol(B);
    2711                 :          0 :   pari_sp av = avma;
    2712         [ #  # ]:          0 :   if (dA < 0) return gen_0;
    2713                 :          0 :   A = Q_primitive_part(A, &c);
    2714                 :          0 :   R = ZX_resultant(B, A);
    2715         [ #  # ]:          0 :   if (c) R = gmul(R, gpowgs(c, dB));
    2716                 :          0 :   lB = leading_term(B);
    2717         [ #  # ]:          0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2718                 :          0 :   return gerepileupto(av, R);
    2719                 :            : }
    2720                 :            : 
    2721                 :            : /* assume x has integral coefficients */
    2722                 :            : GEN
    2723                 :      16197 : ZX_disc_all(GEN x, ulong bound)
    2724                 :            : {
    2725                 :      16197 :   pari_sp av = avma;
    2726                 :            :   GEN l, R;
    2727                 :      16197 :   long s, d = degpol(x);
    2728 [ +  + ][ +  - ]:      16197 :   if (d <= 1) return d ? gen_1: gen_0;
    2729         [ +  + ]:      16047 :   s = (d & 2) ? -1: 1;
    2730                 :      16047 :   l = leading_term(x);
    2731                 :      16047 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2732         [ +  + ]:      16047 :   if (is_pm1(l))
    2733         [ -  + ]:      14117 :   { if (signe(l) < 0) s = -s; }
    2734                 :            :   else
    2735                 :       1930 :     R = diviiexact(R,l);
    2736         [ +  + ]:      16047 :   if (s == -1) togglesign_safe(&R);
    2737                 :      16197 :   return gerepileuptoint(av,R);
    2738                 :            : }
    2739                 :      14882 : GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2740                 :            : 
    2741                 :            : GEN
    2742                 :          0 : QX_disc(GEN x)
    2743                 :            : {
    2744                 :          0 :   pari_sp av = avma;
    2745                 :          0 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2746         [ #  # ]:          0 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2747                 :          0 :   return gerepileupto(av, d);
    2748                 :            : }
    2749                 :            : 
    2750                 :            : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2751                 :            : GEN
    2752                 :      76033 : QXQ_inv(GEN A, GEN B)
    2753                 :            : {
    2754                 :            :   GEN D, cU, q, U, V;
    2755                 :            :   ulong p;
    2756                 :      76033 :   pari_sp av2, av = avma;
    2757                 :            :   forprime_t S;
    2758                 :            : 
    2759         [ -  + ]:      76033 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2760                 :            :   /* A a QX, B a ZX */
    2761         [ +  + ]:      76033 :   if (degpol(A) < 15) return RgXQ_inv(A,B);
    2762                 :       1155 :   A = Q_primitive_part(A, &D);
    2763                 :            :   /* A, B in Z[X] */
    2764                 :       1155 :   init_modular(&S);
    2765                 :       1155 :   av2 = avma; U = NULL;
    2766         [ +  - ]:     103245 :   while ((p = u_forprime_next(&S)))
    2767                 :            :   {
    2768                 :            :     GEN a, b, qp, Up, Vp;
    2769                 :            :     int stable;
    2770                 :            : 
    2771                 :     103245 :     a = ZX_to_Flx(A, p);
    2772                 :     103245 :     b = ZX_to_Flx(B, p);
    2773                 :            :     /* if p | Res(A/G, B/G), discard */
    2774         [ -  + ]:     103245 :     if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
    2775                 :            : 
    2776         [ +  + ]:     103245 :     if (!U)
    2777                 :            :     { /* First time */
    2778                 :       1155 :       U = ZX_init_CRT(Up,p,varn(A));
    2779                 :       1155 :       V = ZX_init_CRT(Vp,p,varn(A));
    2780                 :       1155 :       q = utoipos(p); continue;
    2781                 :            :     }
    2782         [ -  + ]:     102090 :     if (DEBUGLEVEL>5) err_printf("QXQ_inv: mod %ld (bound 2^%ld)", p,expi(q));
    2783                 :     102090 :     qp = muliu(q,p);
    2784                 :     102090 :     stable = ZX_incremental_CRT_raw(&U, Up, q,qp, p)
    2785                 :     102090 :            & ZX_incremental_CRT_raw(&V, Vp, q,qp, p);
    2786         [ +  + ]:     102090 :     if (stable)
    2787                 :            :     { /* all stable: check divisibility */
    2788                 :       1155 :       GEN res = ZX_add(ZX_mul(A,U), ZX_mul(B,V));
    2789         [ +  - ]:       1155 :       if (degpol(res) == 0) {
    2790                 :       1155 :         res = gel(res,2);
    2791         [ +  - ]:       1155 :         D = D? gmul(D, res): res;
    2792                 :            :         break;
    2793                 :            :       } /* DONE */
    2794         [ #  # ]:          0 :       if (DEBUGLEVEL) err_printf("QXQ_inv: char 0 check failed");
    2795                 :            :     }
    2796                 :     100935 :     q = qp;
    2797         [ -  + ]:     100935 :     if (gc_needed(av,1))
    2798                 :            :     {
    2799         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_inv");
    2800                 :     102090 :       gerepileall(av2, 3, &q,&U,&V);
    2801                 :            :     }
    2802                 :            :   }
    2803         [ -  + ]:       1155 :   if (!p) pari_err_OVERFLOW("QXQ_inv [ran out of primes]");
    2804                 :       1155 :   cU = ZX_content(U);
    2805         [ +  - ]:       1155 :   if (!is_pm1(cU)) { U = Q_div_to_int(U, cU); D = gdiv(D, cU); }
    2806                 :      76033 :   return gerepileupto(av, RgX_Rg_div(U, D));
    2807                 :            : }
    2808                 :            : 
    2809                 :            : /************************************************************************
    2810                 :            :  *                                                                      *
    2811                 :            :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2812                 :            :  *                                                                      *
    2813                 :            :  ************************************************************************/
    2814                 :            : 
    2815                 :            : /* irreducible (unitary) polynomial of degree n over Fp */
    2816                 :            : GEN
    2817                 :          0 : ffinit_rand(GEN p,long n)
    2818                 :            : {
    2819                 :            :   for(;;) {
    2820                 :          0 :     pari_sp av = avma;
    2821                 :          0 :     GEN pol = ZX_add(monomial(gen_1, n, 0), random_FpX(n-1,0, p));
    2822         [ #  # ]:          0 :     if (FpX_is_irred(pol, p)) return pol;
    2823                 :          0 :     avma = av;
    2824                 :          0 :   }
    2825                 :            : }
    2826                 :            : 
    2827                 :            : /* return an extension of degree 2^l of F_2, assume l > 0
    2828                 :            :  * Not stack clean. */
    2829                 :            : static GEN
    2830                 :         35 : f2init(long l)
    2831                 :            : {
    2832                 :            :   GEN Q, T, S;
    2833                 :            :   long i, v;
    2834                 :            : 
    2835         [ +  + ]:         35 :   if (l == 1) return polcyclo(3, 0);
    2836                 :         20 :   v = fetch_var_higher();
    2837                 :         20 :   S = mkpoln(4, gen_1,gen_1,gen_0,gen_0); /* y(y^2 + y) */
    2838                 :         20 :   Q = mkpoln(3, gen_1,gen_1, S); /* x^2 + x + y(y^2+y) */
    2839                 :         20 :   setvarn(Q, v);
    2840                 :            : 
    2841                 :            :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2842                 :         20 :   T = mkpoln(5, gen_1,gen_0,gen_0,gen_1,gen_1);
    2843                 :         20 :   setvarn(T, v);
    2844                 :            :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2845                 :            :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2846                 :            :    * ==> x^2 + x + (b^2+b)b */
    2847         [ +  + ]:         50 :   for (i=2; i<l; i++) T = FpX_FpXY_resultant(T, Q, gen_2); /* minpoly(b) / F2*/
    2848                 :         35 :   (void)delete_var(); setvarn(T,0); return T;
    2849                 :            : }
    2850                 :            : 
    2851                 :            : /* return an extension of degree p^l of F_p, assume l > 0
    2852                 :            :  * Not stack clean. */
    2853                 :            : GEN
    2854                 :          0 : ffinit_Artin_Shreier(GEN ip, long l)
    2855                 :            : {
    2856                 :          0 :   long i, v, p = itos(ip);
    2857                 :          0 :   GEN T, Q, xp = monomial(gen_1,p,0); /* x^p */
    2858                 :          0 :   T = ZX_sub(xp, deg1pol_shallow(gen_1,gen_1,0)); /* x^p - x - 1 */
    2859         [ #  # ]:          0 :   if (l == 1) return T;
    2860                 :            : 
    2861                 :          0 :   v = fetch_var_higher();
    2862                 :          0 :   setvarn(xp, v);
    2863                 :          0 :   Q = ZX_sub(monomial(gen_1,2*p-1,0), monomial(gen_1,p,0));
    2864                 :          0 :   Q = gsub(xp, deg1pol_shallow(gen_1, Q, v)); /* x^p - x - (y^(2p-1)-y^p) */
    2865         [ #  # ]:          0 :   for (i = 2; i <= l; ++i) T = FpX_FpXY_resultant(T, Q, ip);
    2866                 :          0 :   (void)delete_var(); setvarn(T,0); return T;
    2867                 :            : }
    2868                 :            : 
    2869                 :            : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    2870                 :            : static long
    2871                 :       7250 : fpinit_check(GEN p, long n, long l)
    2872                 :            : {
    2873                 :            :   ulong q;
    2874         [ +  + ]:       7250 :   if (!uisprime(n)) return 0;
    2875         [ +  + ]:       3655 :   q = umodiu(p,n); if (!q) return 0;
    2876                 :       7250 :   return cgcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2877                 :            : }
    2878                 :            : 
    2879                 :            : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    2880                 :            :  * Return an irreducible polynomial of degree l over F_p.
    2881                 :            :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    2882                 :            :  * finite fields", ACM, 1986 (5) 350--355.
    2883                 :            :  * Not stack clean */
    2884                 :            : static GEN
    2885                 :       1915 : fpinit(GEN p, long l)
    2886                 :            : {
    2887                 :       1915 :   ulong n = 1+l;
    2888         [ +  + ]:       5315 :   while (!fpinit_check(p,n,l)) n += l;
    2889         [ -  + ]:       1915 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    2890                 :       1915 :   return FpX_red(polsubcyclo(n,l,0),p);
    2891                 :            : }
    2892                 :            : 
    2893                 :            : static GEN
    2894                 :       1600 : ffinit_fact(GEN p, long n)
    2895                 :            : {
    2896                 :       1600 :   GEN P, F = gel(factoru_pow(n),3);
    2897                 :            :   long i;
    2898 [ +  + ][ +  + ]:       1600 :   if (!odd(n) && equaliu(p, 2))
    2899                 :         35 :     P = f2init(vals(n)); /* if n is even, F[1] = 2^vals(n)*/
    2900                 :            :   else
    2901                 :       1565 :     P = fpinit(p, F[1]);
    2902         [ +  + ]:       1885 :   for (i = 2; i < lg(F); ++i)
    2903                 :        285 :     P = FpX_direct_compositum(fpinit(p, F[i]), P, p);
    2904                 :       1600 :   return P;
    2905                 :            : }
    2906                 :            : 
    2907                 :            : static GEN
    2908                 :         65 : ffinit_nofact(GEN p, long n)
    2909                 :            : {
    2910                 :         65 :   GEN P, Q = NULL;
    2911         [ -  + ]:         65 :   if (lgefint(p)==3)
    2912                 :            :   {
    2913                 :          0 :     ulong pp = p[2], q;
    2914                 :          0 :     long v = u_lvalrem(n,pp,&q);
    2915         [ #  # ]:          0 :     if (v>0)
    2916                 :            :     {
    2917         [ #  # ]:          0 :       Q = (pp == 2)? f2init(v): fpinit(p,n/q);
    2918                 :          0 :       n = q;
    2919                 :            :     }
    2920                 :            :   }
    2921                 :            :   /* n coprime to p */
    2922         [ -  + ]:         65 :   if (n==1) P = Q;
    2923                 :            :   else
    2924                 :            :   {
    2925                 :         65 :     P = fpinit(p, n);
    2926         [ -  + ]:         65 :     if (Q) P = FpX_direct_compositum(P, Q, p);
    2927                 :            :   }
    2928                 :         65 :   return P;
    2929                 :            : }
    2930                 :            : 
    2931                 :            : static GEN
    2932                 :       2035 : init_Fq_i(GEN p, long n, long v)
    2933                 :            : {
    2934                 :            :   GEN P;
    2935         [ -  + ]:       2035 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    2936         [ -  + ]:       2035 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    2937         [ -  + ]:       2035 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    2938         [ +  + ]:       2035 :   if (v < 0) v = 0;
    2939         [ +  + ]:       2035 :   if (n == 1) return pol_x(v);
    2940         [ +  + ]:       1935 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    2941         [ +  + ]:       1665 :   if (lgefint(p)-2 <= expu(n))
    2942                 :       1600 :     P = ffinit_fact(p,n);
    2943                 :            :   else
    2944                 :         65 :     P = ffinit_nofact(p,n);
    2945                 :       2035 :   setvarn(P, v); return P;
    2946                 :            : }
    2947                 :            : GEN
    2948                 :       1995 : init_Fq(GEN p, long n, long v)
    2949                 :            : {
    2950                 :       1995 :   pari_sp av = avma;
    2951                 :       1995 :   return gerepileupto(av, init_Fq_i(p, n, v));
    2952                 :            : }
    2953                 :            : GEN
    2954                 :         40 : ffinit(GEN p, long n, long v)
    2955                 :            : {
    2956                 :         40 :   pari_sp av = avma;
    2957                 :         40 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    2958                 :            : }
    2959                 :            : 
    2960                 :            : GEN
    2961                 :        860 : ffnbirred(GEN p, long n)
    2962                 :            : {
    2963                 :        860 :   pari_sp av = avma;
    2964                 :            :   long j;
    2965                 :        860 :   GEN s = gen_0, dk, pd;
    2966                 :        860 :   dk = divisorsu(n);
    2967         [ +  + ]:       2855 :   for (j = 1; j < lg(dk); ++j)
    2968                 :            :   {
    2969                 :       1995 :     long d = dk[j];
    2970                 :       1995 :     long m = moebiusu(d);
    2971         [ +  + ]:       1995 :     if (!m) continue;
    2972                 :       1835 :     pd = powiu(p, n/d);
    2973         [ +  + ]:       1835 :     s = m>0 ? addii(s, pd): subii(s,pd);
    2974                 :            :   }
    2975                 :        860 :   return gerepileuptoint(av, divis(s, n));
    2976                 :            : }
    2977                 :            : 
    2978                 :            : GEN
    2979                 :        130 : ffsumnbirred(GEN p, long n)
    2980                 :            : {
    2981                 :        130 :   pari_sp av = avma;
    2982                 :            :   long i,j;
    2983                 :        130 :   GEN v,q, t = gen_0;
    2984                 :        130 :   v = cgetg(n+1,t_VECSMALL); v[1] = 1;
    2985                 :        130 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    2986         [ +  + ]:        610 :   for(i=2; i<=n; i++)
    2987                 :            :   {
    2988                 :        480 :     v[i] = moebiusu(i);
    2989                 :        480 :     gel(q,i) = mulii(gel(q,i-1), p);
    2990                 :            :   }
    2991         [ +  + ]:        740 :   for(i=1; i<=n; i++)
    2992                 :            :   {
    2993                 :        610 :     GEN s = gen_0;
    2994                 :        610 :     GEN dk = divisorsu(i);
    2995         [ +  + ]:       1940 :     for (j = 1; j < lg(dk); ++j)
    2996                 :            :     {
    2997                 :       1330 :       long d = dk[j], m = v[d];
    2998         [ +  + ]:       1330 :       if (!m) continue;
    2999         [ +  + ]:       1190 :       s = m>0 ? addii(s, gel(q, i/d)): subii(s, gel(q, i/d));
    3000                 :            :     }
    3001                 :        610 :     t = addii(t, divis(s, i));
    3002                 :            :   }
    3003                 :        130 :   return gerepileuptoint(av, t);
    3004                 :            : }
    3005                 :            : 
    3006                 :            : GEN
    3007                 :        100 : ffnbirred0(GEN p, long n, long flag)
    3008                 :            : {
    3009         [ -  + ]:        100 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    3010         [ -  + ]:        100 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    3011      [ +  +  - ]:        100 :   switch(flag)
    3012                 :            :   {
    3013                 :            :     case 0:
    3014                 :         50 :       return ffnbirred(p, n);
    3015                 :            :     case 1:
    3016                 :         50 :       return ffsumnbirred(p, n);
    3017                 :            :     default:
    3018                 :          0 :       pari_err_FLAG("ffnbirred");
    3019                 :            :   }
    3020                 :        100 :   return NULL; /* NOT REACHED */
    3021                 :            : }

Generated by: LCOV version 1.9