Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - polarit3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20783-cec4728) Lines: 1048 1232 85.1 %
Date: 2017-06-28 05:59:20 Functions: 114 128 89.1 %
Legend: Lines: hit not hit

          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        3311 : char_update_prime(struct charact *S, GEN p)
      34             : {
      35        3311 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      36        3311 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      37        3304 : }
      38             : static void
      39        3479 : char_update_int(struct charact *S, GEN n)
      40             : {
      41        3479 :   if (S->isprime)
      42             :   {
      43        3479 :     if (dvdii(n, S->q)) return;
      44           7 :     pari_err_MODULUS("characteristic", S->q, n);
      45             :   }
      46        3472 :   S->q = gcdii(S->q, n);
      47             : }
      48             : static void
      49      578228 : charact(struct charact *S, GEN x)
      50             : {
      51      578228 :   const long tx = typ(x);
      52             :   long i, l;
      53      578228 :   switch(tx)
      54             :   {
      55        3213 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      56        3248 :     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        9177 :       l = lg(x);
      61        9177 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      62        9163 :       break;
      63             :     case t_LIST:
      64           7 :       x = list_data(x);
      65           7 :       if (x) charact(S, x);
      66           7 :       break;
      67             :   }
      68      578200 : }
      69             : static void
      70       30296 : charact_res(struct charact *S, GEN x)
      71             : {
      72       30296 :   const long tx = typ(x);
      73             :   long i, l;
      74       30296 :   switch(tx)
      75             :   {
      76         266 :     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          63 :     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        9317 :       l = lg(x);
      83        9317 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      84        9317 :       break;
      85             :     case t_LIST:
      86           0 :       x = list_data(x);
      87           0 :       if (x) charact_res(S, x);
      88           0 :       break;
      89             :   }
      90       30296 : }
      91             : GEN
      92        6531 : characteristic(GEN x)
      93             : {
      94             :   struct charact S;
      95        6531 :   S.q = gen_0; S.isprime = 0;
      96        6531 :   charact(&S, x); return S.q;
      97             : }
      98             : GEN
      99        2324 : residual_characteristic(GEN x)
     100             : {
     101             :   struct charact S;
     102        2324 :   S.q = gen_0; S.isprime = 0;
     103        2324 :   charact_res(&S, x); return S.q;
     104             : }
     105             : 
     106             : int
     107    22926214 : Rg_is_Fp(GEN x, GEN *pp)
     108             : {
     109             :   GEN mod;
     110    22926214 :   switch(typ(x))
     111             :   {
     112             :   case t_INTMOD:
     113     3933559 :     mod = gel(x,1);
     114     3933559 :     if (!*pp) *pp = mod;
     115     3882480 :     else if (mod != *pp && !equalii(mod, *pp))
     116             :     {
     117           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_Fp");
     118           0 :       return 0;
     119             :     }
     120     3933559 :     return 1;
     121             :   case t_INT:
     122    16480268 :     return 1;
     123     2512387 :   default: return 0;
     124             :   }
     125             : }
     126             : 
     127             : int
     128     2561340 : RgX_is_FpX(GEN x, GEN *pp)
     129             : {
     130     2561340 :   long i, lx = lg(x);
     131    10139105 :   for (i=2; i<lx; i++)
     132     9326929 :     if (!Rg_is_Fp(gel(x, i), pp))
     133     1749164 :       return 0;
     134      812176 :   return 1;
     135             : }
     136             : 
     137             : int
     138     2113088 : RgV_is_FpV(GEN x, GEN *pp)
     139             : {
     140     2113088 :   long i, lx = lg(x);
     141    14937306 :   for (i=1; i<lx; i++)
     142    13587434 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     143     1349872 :   return 1;
     144             : }
     145             : 
     146             : int
     147      911390 : RgM_is_FpM(GEN x, GEN *pp)
     148             : {
     149      911390 :   long i, lx = lg(x);
     150     2258861 :   for (i=1; i<lx; i++)
     151     2110666 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     152      148195 :   return 1;
     153             : }
     154             : 
     155             : int
     156       81186 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     157             : {
     158             :   GEN pol, mod, p;
     159       81186 :   switch(typ(x))
     160             :   {
     161             :   case t_INTMOD:
     162       11781 :     return Rg_is_Fp(x, pp);
     163             :   case t_INT:
     164       20608 :     return 1;
     165             :   case t_POL:
     166        8624 :     return RgX_is_FpX(x, pp);
     167             :   case t_FFELT:
     168       10038 :     mod = FF_1(x); p = FF_p_i(x);
     169       10038 :     if (!*pp) *pp = p;
     170       10038 :     if (!*pT) *pT = mod;
     171       10038 :     if ((p != *pp && !equalii(p, *pp)) || (mod != *pT && !gequal(mod, *pT)))
     172             :     {
     173           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     174           0 :       return 0;
     175             :     }
     176       10038 :     return 1;
     177             :   case t_POLMOD:
     178        5544 :     mod = gel(x,1); pol = gel(x, 2);
     179        5544 :     if (!RgX_is_FpX(mod, pp)) return 0;
     180        5544 :     if (typ(pol)==t_POL)
     181             :     {
     182        5474 :       if (!RgX_is_FpX(pol, pp)) return 0;
     183             :     }
     184          70 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     185        5523 :     if (!*pT) *pT = mod;
     186        4816 :     else if (mod != *pT && !gequal(mod, *pT))
     187             :     {
     188           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     189           0 :       return 0;
     190             :     }
     191        5523 :     return 1;
     192             : 
     193       24591 :   default: return 0;
     194             :   }
     195             : }
     196             : 
     197             : int
     198       31528 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     199             : {
     200       31528 :   long i, lx = lg(x);
     201       86814 :   for (i = 2; i < lx; i++)
     202       80640 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     203        6174 :   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    30676573 : Rg_to_Fp(GEN x, GEN p)
     216             : {
     217    30676573 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     218      397737 :   switch(typ(x))
     219             :   {
     220       24296 :     case t_INT: return modii(x, p);
     221             :     case t_FRAC: {
     222          60 :       pari_sp av = avma;
     223          60 :       GEN z = modii(gel(x,1), p);
     224          60 :       if (z == gen_0) return gen_0;
     225          60 :       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      373381 :       GEN q = gel(x,1), a = gel(x,2);
     230      373381 :       if (equalii(q, p)) return icopy(a);
     231          14 :       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             :       return NULL; /* LCOV_EXCL_LINE */
     236             :   }
     237             : }
     238             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     239             : GEN
     240       45641 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     241             : {
     242       45641 :   long ta, tx = typ(x), v = get_FpX_var(T);
     243             :   GEN a, b;
     244       45641 :   if (is_const_t(tx))
     245             :   {
     246       41154 :     if (tx == t_FFELT)
     247             :     {
     248       21350 :       GEN z = FF_to_FpXQ(x);
     249       21350 :       setvarn(z, v);
     250       21350 :       return z;
     251             :     }
     252       19804 :     return scalar_ZX(Rg_to_Fp(x, p), v);
     253             :   }
     254        4487 :   switch(tx)
     255             :   {
     256             :     case t_POLMOD:
     257        3724 :       b = gel(x,1);
     258        3724 :       a = gel(x,2); ta = typ(a);
     259        3724 :       if (is_const_t(ta)) return scalar_ZX(Rg_to_Fp(a, p), v);
     260        3682 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     261        3682 :       a = RgX_to_FpX(a, p); if (ZX_equal(b,get_FpX_mod(T))) return a;
     262           0 :       if (signe(FpX_rem(b,T,p))==0) return FpX_rem(a, T, p);
     263           0 :       break;
     264             :     case t_POL:
     265         763 :       if (varn(x) != v) break;
     266         763 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     267             :     case t_RFRAC:
     268           0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     269           0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     270           0 :       return FpXQ_div(a,b, T,p);
     271             :   }
     272           0 :   pari_err_TYPE("Rg_to_FpXQ",x);
     273             :   return NULL; /* LCOV_EXCL_LINE */
     274             : }
     275             : GEN
     276      140147 : RgX_to_FpX(GEN x, GEN p)
     277             : {
     278             :   long i, l;
     279      140147 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     280      140147 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     281      140147 :   return FpX_renormalize(z, l);
     282             : }
     283             : 
     284             : GEN
     285        1099 : RgV_to_FpV(GEN x, GEN p)
     286             : {
     287        1099 :   long i, l = lg(x);
     288        1099 :   GEN z = cgetg(l, t_VEC);
     289        1099 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     290        1099 :   return z;
     291             : }
     292             : 
     293             : GEN
     294      915291 : RgC_to_FpC(GEN x, GEN p)
     295             : {
     296      915291 :   long i, l = lg(x);
     297      915291 :   GEN z = cgetg(l, t_COL);
     298      915291 :   for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     299      915291 :   return z;
     300             : }
     301             : 
     302             : GEN
     303       98864 : RgM_to_FpM(GEN x, GEN p)
     304             : {
     305       98864 :   long i, l = lg(x);
     306       98864 :   GEN z = cgetg(l, t_MAT);
     307       98864 :   for (i = 1; i < l; i++) gel(z,i) = RgC_to_FpC(gel(x,i), p);
     308       98864 :   return z;
     309             : }
     310             : GEN
     311       17116 : RgV_to_Flv(GEN x, ulong p)
     312             : {
     313       17116 :   long l = lg(x), i;
     314       17116 :   GEN a = cgetg(l, t_VECSMALL);
     315       17116 :   for (i=1; i<l; i++) a[i] = Rg_to_Fl(gel(x,i), p);
     316       17116 :   return a;
     317             : }
     318             : GEN
     319        2034 : RgM_to_Flm(GEN x, ulong p)
     320             : {
     321             :   long l, i;
     322        2034 :   GEN a = cgetg_copy(x, &l);
     323        2034 :   for (i=1; i<l; i++) gel(a,i) = RgV_to_Flv(gel(x,i), p);
     324        2034 :   return a;
     325             : }
     326             : 
     327             : GEN
     328         336 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     329             : {
     330         336 :   long i, l = lg(x);
     331         336 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     332         336 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     333         336 :   return FpXQX_renormalize(z, l);
     334             : }
     335             : GEN
     336         588 : RgX_to_FqX(GEN x, GEN T, GEN p)
     337             : {
     338         588 :   long i, l = lg(x);
     339         588 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     340         588 :   if (T)
     341       10563 :     for (i = 2; i < l; i++)
     342       10017 :       gel(z,i) = simplify_shallow(Rg_to_FpXQ(gel(x,i), T, p));
     343             :   else
     344        1554 :     for (i = 2; i < l; i++)
     345        1512 :       gel(z,i) = Rg_to_Fp(gel(x,i), p);
     346         588 :   return FpXQX_renormalize(z, l);
     347             : }
     348             : 
     349             : /* lg(V) > 1 */
     350             : GEN
     351      849765 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     352             : {
     353      849765 :   pari_sp av = avma;
     354      849765 :   long i, l = lg(V);
     355      849765 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     356     4181499 :   for(i=2; i<l; i++)
     357             :   {
     358     3331734 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     359     3331734 :     if ((i & 7) == 0) z = gerepileupto(av, z);
     360             :   }
     361      849765 :   return gerepileupto(av, FpX_red(z,p));
     362             : }
     363             : 
     364             : GEN
     365        2254 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     366             : {
     367        2254 :   long i, lz = lg(y);
     368             :   GEN z;
     369        2254 :   if (!T) return FpX_Fp_add(y, x, p);
     370        2254 :   if (lz == 2) return scalarpol(x, varn(y));
     371        2254 :   z = cgetg(lz,t_POL); z[1] = y[1];
     372        2254 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     373        2254 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     374             :   else
     375        1246 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     376        2254 :   return z;
     377             : }
     378             : 
     379             : GEN
     380       63361 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     381             : {
     382             :   long i, lP;
     383       63361 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     384       63361 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     385       63361 :   gel(res,lP-1) = gen_1; return res;
     386             : }
     387             : 
     388             : GEN
     389        2947 : FpXQX_normalize(GEN z, GEN T, GEN p)
     390             : {
     391             :   GEN lc;
     392        2947 :   if (lg(z) == 2) return z;
     393        2947 :   lc = leading_coeff(z);
     394        2947 :   if (typ(lc) == t_POL)
     395             :   {
     396        1267 :     if (lg(lc) > 3) /* non-constant */
     397        1236 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     398             :     /* constant */
     399          31 :     lc = gel(lc,2);
     400          31 :     z = shallowcopy(z);
     401          31 :     gel(z, lg(z)-1) = lc;
     402             :   }
     403             :   /* lc a t_INT */
     404        1711 :   if (equali1(lc)) return z;
     405          35 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     406             : }
     407             : 
     408             : GEN
     409      123403 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     410             : {
     411             :   pari_sp av;
     412             :   GEN p1, r;
     413      123403 :   long j, i=lg(x)-1;
     414      123403 :   if (i<=2)
     415       24136 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     416       99267 :   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      291207 :   for (i--; i>=2; i=j-1)
     420             :   {
     421      192255 :     for (j=i; !signe(gel(x,j)); j--)
     422         315 :       if (j==2)
     423             :       {
     424         182 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     425         182 :         return gerepileupto(av, Fq_mul(p1,y, T, p));
     426             :       }
     427      191940 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     428      191940 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     429             :   }
     430       99085 :   return gerepileupto(av, p1);
     431             : }
     432             : 
     433             : GEN
     434       30380 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     435             : {
     436       30380 :   long i, lb = lg(Q);
     437             :   GEN z;
     438       30380 :   if (!T) return FpXY_evalx(Q, x, p);
     439       20720 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     440      115969 :   for (i=2; i<lb; i++)
     441             :   {
     442       95249 :     GEN q = gel(Q,i);
     443       95249 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     444             :   }
     445       20720 :   return FpXQX_renormalize(z, lb);
     446             : }
     447             : 
     448             : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     449             : GEN
     450       12733 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     451             : {
     452       12733 :   pari_sp av = avma;
     453       12733 :   if (!T) return FpXY_eval(Q, y, x, p);
     454         336 :   return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     455             : }
     456             : 
     457             : /* a X^d */
     458             : GEN
     459      559293 : monomial(GEN a, long d, long v)
     460             : {
     461             :   long i, n;
     462             :   GEN P;
     463      559293 :   if (d < 0) {
     464           0 :     if (isrationalzero(a)) return pol_0(v);
     465           0 :     retmkrfrac(a, pol_xn(-d, v));
     466             :   }
     467      559293 :   if (gequal0(a))
     468             :   {
     469         364 :     if (isexactzero(a)) return scalarpol_shallow(a,v);
     470           0 :     n = d+2; P = cgetg(n+1, t_POL);
     471           0 :     P[1] = evalsigne(0) | evalvarn(v);
     472             :   }
     473             :   else
     474             :   {
     475      558929 :     n = d+2; P = cgetg(n+1, t_POL);
     476      558929 :     P[1] = evalsigne(1) | evalvarn(v);
     477             :   }
     478      558929 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     479      558929 :   gel(P,i) = a; return P;
     480             : }
     481             : GEN
     482     7580396 : monomialcopy(GEN a, long d, long v)
     483             : {
     484             :   long i, n;
     485             :   GEN P;
     486     7580396 :   if (d < 0) {
     487           7 :     if (isrationalzero(a)) return pol_0(v);
     488           7 :     retmkrfrac(gcopy(a), pol_xn(-d, v));
     489             :   }
     490     7580389 :   if (gequal0(a))
     491             :   {
     492           7 :     if (isexactzero(a)) return scalarpol(a,v);
     493           0 :     n = d+2; P = cgetg(n+1, t_POL);
     494           0 :     P[1] = evalsigne(0) | evalvarn(v);
     495             :   }
     496             :   else
     497             :   {
     498     7580382 :     n = d+2; P = cgetg(n+1, t_POL);
     499     7580382 :     P[1] = evalsigne(1) | evalvarn(v);
     500             :   }
     501     7580382 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     502     7580382 :   gel(P,i) = gcopy(a); return P;
     503             : }
     504             : GEN
     505       19866 : pol_x_powers(long N, long v)
     506             : {
     507       19866 :   GEN L = cgetg(N+1,t_VEC);
     508             :   long i;
     509       19866 :   for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
     510       19866 :   return L;
     511             : }
     512             : 
     513             : GEN
     514           0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     515             : {
     516           0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     517             : }
     518             : 
     519             : GEN
     520           0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     521             : {
     522           0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     523             : }
     524             : 
     525             : /*******************************************************************/
     526             : /*                                                                 */
     527             : /*                             Fq                                  */
     528             : /*                                                                 */
     529             : /*******************************************************************/
     530             : 
     531             : GEN
     532     6606086 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     533             : {
     534             :   (void)T;
     535     6606086 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     536             :   {
     537     2382232 :     case 0: return Fp_add(x,y,p);
     538      204162 :     case 1: return FpX_Fp_add(x,y,p);
     539      338045 :     case 2: return FpX_Fp_add(y,x,p);
     540     3681647 :     case 3: return FpX_add(x,y,p);
     541             :   }
     542           0 :   return NULL;
     543             : }
     544             : 
     545             : GEN
     546     4423614 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     547             : {
     548             :   (void)T;
     549     4423614 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     550             :   {
     551      155468 :     case 0: return Fp_sub(x,y,p);
     552        1188 :     case 1: return FpX_Fp_sub(x,y,p);
     553        8042 :     case 2: return Fp_FpX_sub(x,y,p);
     554     4258916 :     case 3: return FpX_sub(x,y,p);
     555             :   }
     556           0 :   return NULL;
     557             : }
     558             : 
     559             : GEN
     560      457795 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     561             : {
     562             :   (void)T;
     563      457795 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     564             : }
     565             : 
     566             : GEN
     567       11464 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     568             : {
     569             :   (void)T;
     570       11464 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     571             : }
     572             : 
     573             : /* If T==NULL do not reduce*/
     574             : GEN
     575    42207684 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     576             : {
     577    42207684 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     578             :   {
     579     2380747 :     case 0: return Fp_mul(x,y,p);
     580      187132 :     case 1: return FpX_Fp_mul(x,y,p);
     581      316439 :     case 2: return FpX_Fp_mul(y,x,p);
     582    39323366 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     583     2710088 :             else return FpX_mul(x,y,p);
     584             :   }
     585           0 :   return NULL;
     586             : }
     587             : 
     588             : /* If T==NULL do not reduce*/
     589             : GEN
     590      729973 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     591             : {
     592             :   (void) T;
     593      729973 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     594             : }
     595             : 
     596             : /* y t_INT */
     597             : GEN
     598       51652 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     599             : {
     600             :   (void)T;
     601      103304 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     602       51652 :                           : Fp_mul(x,y,p);
     603             : }
     604             : /* If T==NULL do not reduce*/
     605             : GEN
     606      261379 : Fq_sqr(GEN x, GEN T, GEN p)
     607             : {
     608      261379 :   if (typ(x) == t_POL)
     609             :   {
     610       11206 :     if (T) return FpXQ_sqr(x,T,p);
     611           0 :     else return FpX_sqr(x,p);
     612             :   }
     613             :   else
     614      250173 :     return Fp_sqr(x,p);
     615             : }
     616             : 
     617             : GEN
     618           0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     619             : {
     620           0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     621           0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     622             : }
     623             : 
     624             : GEN
     625           0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     626             : {
     627           0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     628           0 :   return FpXQ_invsafe(x,pol,p);
     629             : }
     630             : 
     631             : GEN
     632       20247 : Fq_inv(GEN x, GEN pol, GEN p)
     633             : {
     634       20247 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     635       15508 :   return FpXQ_inv(x,pol,p);
     636             : }
     637             : 
     638             : GEN
     639      479192 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     640             : {
     641      479192 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     642             :   {
     643      451976 :     case 0: return Fp_div(x,y,p);
     644       22806 :     case 1: return FpX_Fp_mul(x,Fp_inv(y,p),p);
     645         196 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     646        4214 :     case 3: return FpXQ_div(x,y,pol,p);
     647             :   }
     648           0 :   return NULL;
     649             : }
     650             : 
     651             : GEN
     652       24640 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     653             : {
     654       24640 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     655       12110 :   return FpXQ_pow(x,n,pol,p);
     656             : }
     657             : 
     658             : GEN
     659       12985 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     660             : {
     661       12985 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     662         553 :   return FpXQ_powu(x,n,pol,p);
     663             : }
     664             : 
     665             : GEN
     666      709049 : Fq_sqrt(GEN x, GEN T, GEN p)
     667             : {
     668      709049 :   if (typ(x) == t_INT)
     669             :   {
     670      698544 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     671         287 :     x = scalarpol_shallow(x, get_FpX_var(T));
     672             :   }
     673       10792 :   return FpXQ_sqrt(x,T,p);
     674             : }
     675             : GEN
     676       60600 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     677             : {
     678       60600 :   if (typ(x) == t_INT)
     679             :   {
     680             :     long d;
     681       60390 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     682         610 :     d = get_FpX_degree(T);
     683         610 :     if (ugcd(umodiu(n,d),d) == 1)
     684             :     {
     685         414 :       if (!zeta)
     686           7 :         return Fp_sqrtn(x,n,p,NULL);
     687             :       else
     688             :       {
     689             :         /* gcd(n,p-1)=gcd(n,p^d-1) <=> same number of solutions if Fp and F_{p^d} */
     690         407 :         if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     691         386 :           return Fp_sqrtn(x,n,p,zeta);
     692             :       }
     693             :     }
     694         217 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     695             :   }
     696         427 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     697             : }
     698             : 
     699             : struct _Fq_field
     700             : {
     701             :   GEN T, p;
     702             : };
     703             : 
     704             : static GEN
     705       47985 : _Fq_red(void *E, GEN x)
     706       47985 : { struct _Fq_field *s = (struct _Fq_field *)E;
     707       47985 :   return Fq_red(x, s->T, s->p);
     708             : }
     709             : 
     710             : static GEN
     711      601289 : _Fq_add(void *E, GEN x, GEN y)
     712             : {
     713             :   (void) E;
     714      601289 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     715             :   {
     716          14 :     case 0: return addii(x,y);
     717           0 :     case 1: return ZX_Z_add(x,y);
     718         216 :     case 2: return ZX_Z_add(y,x);
     719      601059 :     default: return ZX_add(x,y);
     720             :   }
     721             : }
     722             : 
     723             : static GEN
     724        1784 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     725             : 
     726             : static GEN
     727      622735 : _Fq_mul(void *E, GEN x, GEN y)
     728             : {
     729             :   (void) E;
     730      622735 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     731             :   {
     732          63 :     case 0: return mulii(x,y);
     733         315 :     case 1: return ZX_Z_mul(x,y);
     734          56 :     case 2: return ZX_Z_mul(y,x);
     735      622301 :     default: return ZX_mul(x,y);
     736             :   }
     737             : }
     738             : 
     739             : static GEN
     740        1441 : _Fq_inv(void *E, GEN x)
     741        1441 : { struct _Fq_field *s = (struct _Fq_field *)E;
     742        1441 :   return Fq_inv(x,s->T,s->p);
     743             : }
     744             : 
     745             : static int
     746       25714 : _Fq_equal0(GEN x) { return signe(x)==0; }
     747             : 
     748             : static GEN
     749        1891 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     750             : 
     751             : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     752             :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     753             : 
     754         219 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     755             : {
     756         219 :   GEN z = new_chunk(sizeof(struct _Fq_field));
     757         219 :   struct _Fq_field *e = (struct _Fq_field *) z;
     758         219 :   e->T = T; e->p  = p; *E = (void*)e;
     759         219 :   return &Fq_field;
     760             : }
     761             : 
     762             : /*******************************************************************/
     763             : /*                                                                 */
     764             : /*                             Fq[X]                               */
     765             : /*                                                                 */
     766             : /*******************************************************************/
     767             : /* P(X + c) */
     768             : GEN
     769           0 : FpX_translate(GEN P, GEN c, GEN p)
     770             : {
     771           0 :   pari_sp av = avma;
     772             :   GEN Q, *R;
     773             :   long i, k, n;
     774             : 
     775           0 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     776           0 :   Q = leafcopy(P);
     777           0 :   R = (GEN*)(Q+2); n = degpol(P);
     778           0 :   for (i=1; i<=n; i++)
     779             :   {
     780           0 :     for (k=n-i; k<n; k++)
     781           0 :       R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
     782             : 
     783           0 :     if (gc_needed(av,2))
     784             :     {
     785           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
     786           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     787             :     }
     788             :   }
     789           0 :   return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
     790             : }
     791             : /* P(X + c), c an Fq */
     792             : GEN
     793       39466 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
     794             : {
     795       39466 :   pari_sp av = avma;
     796             :   GEN Q, *R;
     797             :   long i, k, n;
     798             : 
     799             :   /* signe works for t_(INT|POL) */
     800       39466 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     801       39466 :   Q = leafcopy(P);
     802       39466 :   R = (GEN*)(Q+2); n = degpol(P);
     803      172935 :   for (i=1; i<=n; i++)
     804             :   {
     805      486717 :     for (k=n-i; k<n; k++)
     806      353248 :       R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
     807             : 
     808      133469 :     if (gc_needed(av,2))
     809             :     {
     810           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
     811           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     812             :     }
     813             :   }
     814       39466 :   return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
     815             : }
     816             : 
     817             : GEN
     818         630 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     819             : {
     820         630 :   pari_sp ltop = avma;
     821             :   long k;
     822             :   GEN W;
     823         630 :   if (lgefint(p) == 3)
     824             :   {
     825         591 :     ulong pp = p[2];
     826         591 :     GEN Tl = ZX_to_Flx(T, pp);
     827         591 :     GEN Vl = FqV_to_FlxV(V, T, p);
     828         591 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     829         591 :     return gerepileupto(ltop, FlxX_to_ZXX(Tl));
     830             :   }
     831          39 :   W = cgetg(lg(V),t_VEC);
     832         255 :   for(k=1; k < lg(V); k++)
     833         216 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     834          39 :   return gerepileupto(ltop, FpXQXV_prod(W, T, p));
     835             : }
     836             : 
     837             : GEN
     838      123277 : FqV_red(GEN z, GEN T, GEN p)
     839             : {
     840      123277 :   long i, l = lg(z);
     841      123277 :   GEN res = cgetg(l, typ(z));
     842      123277 :   for(i=1;i<l;i++) gel(res,i) = Fq_red(gel(z,i),T,p);
     843      123277 :   return res;
     844             : }
     845             : 
     846             : GEN
     847           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     848             : {
     849           0 :   long i, lx = lg(x);
     850             :   GEN z;
     851           0 :   if (!T) return FpC_add(x, y, p);
     852           0 :   z = cgetg(lx, t_COL);
     853           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_add(gel(x, i), gel(y, i), T, p);
     854           0 :   return z;
     855             : }
     856             : 
     857             : GEN
     858           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     859             : {
     860           0 :   long i, lx = lg(x);
     861             :   GEN z;
     862           0 :   if (!T) return FpC_sub(x, y, p);
     863           0 :   z = cgetg(lx, t_COL);
     864           0 :   for (i = 1; i < lx; i++) gel(z, i) = Fq_sub(gel(x, i), gel(y, i), T, p);
     865           0 :   return z;
     866             : }
     867             : 
     868             : GEN
     869           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
     870             : {
     871           0 :   long i, l = lg(x);
     872             :   GEN z;
     873           0 :   if (!T) return FpC_Fp_mul(x, y, p);
     874           0 :   z = cgetg(l, t_COL);
     875           0 :   for (i=1;i<l;i++) gel(z,i) = Fq_mul(gel(x,i),y,T,p);
     876           0 :   return z;
     877             : }
     878             : 
     879             : GEN
     880         591 : FqV_to_FlxV(GEN v, GEN T, GEN pp)
     881             : {
     882         591 :   long j, N = lg(v);
     883         591 :   long vT = evalvarn(get_FpX_var(T));
     884         591 :   ulong p = pp[2];
     885         591 :   GEN y = cgetg(N, t_VEC);
     886        2874 :   for (j=1; j<N; j++)
     887        4566 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     888        2283 :                                     : ZX_to_Flx(gel(v,j), p));
     889         591 :   return y;
     890             : }
     891             : 
     892             : GEN
     893       40623 : FqC_to_FlxC(GEN v, GEN T, GEN pp)
     894             : {
     895       40623 :   long j, N = lg(v);
     896       40623 :   long vT = evalvarn(get_FpX_var(T));
     897       40623 :   ulong p = pp[2];
     898       40623 :   GEN y = cgetg(N, t_COL);
     899     1079542 :   for (j=1; j<N; j++)
     900     2329407 :     gel(y,j) = (typ(gel(v,j))==t_INT?  Z_to_Flx(gel(v,j), p, vT)
     901     1290488 :                                     : ZX_to_Flx(gel(v,j), p));
     902       40623 :   return y;
     903             : }
     904             : 
     905             : GEN
     906        7733 : FqM_to_FlxM(GEN x, GEN T, GEN pp)
     907             : {
     908        7733 :   long j, n = lg(x);
     909        7733 :   GEN y = cgetg(n,t_MAT);
     910        7733 :   if (n == 1) return y;
     911       48356 :   for (j=1; j<n; j++)
     912       40623 :     gel(y,j) = FqC_to_FlxC(gel(x,j), T, pp);
     913        7733 :   return y;
     914             : }
     915             : 
     916             : /*******************************************************************/
     917             : /*                                                                 */
     918             : /*                          MODULAR GCD                            */
     919             : /*                                                                 */
     920             : /*******************************************************************/
     921             : /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
     922             : static GEN
     923    10773271 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
     924             : {
     925    10773271 :   ulong d, amod = umodiu(a, p);
     926    10773271 :   pari_sp av = avma;
     927             :   GEN ax;
     928             : 
     929    10773271 :   if (b == amod) return NULL;
     930     2236854 :   d = (b > amod)? b - amod: p - (amod - b); /* (b - a) mod p */
     931     2236854 :   (void)new_chunk(lgefint(pq)<<1); /* HACK */
     932     2236854 :   ax = mului(Fl_mul(d,qinv,p), q); /* d mod p, 0 mod q */
     933     2236854 :   avma = av; return addii(a, ax); /* in ]-q, pq[ assuming a in -]-q,q[ */
     934             : }
     935             : GEN
     936       16912 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
     937             : GEN
     938     3116759 : ZX_init_CRT(GEN Hp, ulong p, long v)
     939             : {
     940     3116759 :   long i, l = lg(Hp), lim = (long)(p>>1);
     941     3116759 :   GEN H = cgetg(l, t_POL);
     942     3116759 :   H[1] = evalsigne(1) | evalvarn(v);
     943    11015453 :   for (i=2; i<l; i++)
     944     7898694 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
     945     3116759 :   return H;
     946             : }
     947             : 
     948             : GEN
     949      217425 : ZM_init_CRT(GEN Hp, ulong p)
     950             : {
     951      217425 :   long i,j, m, l = lg(Hp), lim = (long)(p>>1);
     952      217425 :   GEN c, cp, H = cgetg(l, t_MAT);
     953      217425 :   if (l==1) return H;
     954      174018 :   m = lgcols(Hp);
     955     1073638 :   for (j=1; j<l; j++)
     956             :   {
     957      899620 :     cp = gel(Hp,j);
     958      899620 :     c = cgetg(m, t_COL);
     959      899620 :     gel(H,j) = c;
     960      899620 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
     961             :   }
     962      174018 :   return H;
     963             : }
     964             : 
     965             : int
     966       58087 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
     967             : {
     968       58087 :   GEN h, q = *ptq, qp = muliu(q,p), lim = shifti(qp,-1);
     969       58087 :   ulong qinv = Fl_inv(umodiu(q,p), p);
     970       58087 :   int stable = 1;
     971       58087 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp);
     972       58087 :   if (h)
     973             :   {
     974       14468 :     if (cmpii(h,lim) > 0) h = subii(h,qp);
     975       14468 :     *H = h; stable = 0;
     976             :   }
     977       58087 :   *ptq = qp; return stable;
     978             : }
     979             : 
     980             : static int
     981        5903 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
     982             : {
     983        5903 :   GEN H = *ptH, h, lim = shifti(qp,-1);
     984        5903 :   ulong qinv = Fl_inv(umodiu(q,p), p);
     985        5903 :   long i, l = lg(H), lp = lg(Hp);
     986        5903 :   int stable = 1;
     987             : 
     988        5903 :   if (l < lp)
     989             :   { /* degree increases */
     990           0 :     GEN x = cgetg(lp, t_POL);
     991           0 :     for (i=1; i<l; i++)  x[i] = H[i];
     992           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
     993           0 :     *ptH = H = x;
     994           0 :     stable = 0;
     995        5903 :   } else if (l > lp)
     996             :   { /* degree decreases */
     997           0 :     GEN x = cgetg(l, t_VECSMALL);
     998           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
     999           0 :     for (   ; i<l; i++) x[i] = 0;
    1000           0 :     Hp = x; lp = l;
    1001             :   }
    1002      210093 :   for (i=2; i<lp; i++)
    1003             :   {
    1004      204190 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp);
    1005      204190 :     if (h)
    1006             :     {
    1007       54985 :       if (cmpii(h,lim) > 0) h = subii(h,qp);
    1008       54985 :       gel(H,i) = h; stable = 0;
    1009             :     }
    1010             :   }
    1011        5903 :   return stable;
    1012             : }
    1013             : 
    1014             : int
    1015        4172 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1016             : {
    1017        4172 :   GEN q = *ptq, qp = muliu(q,p);
    1018        4172 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1019        4172 :   *ptq = qp; return stable;
    1020             : }
    1021             : 
    1022             : int
    1023      144813 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1024             : {
    1025      144813 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), lim = shifti(qp,-1);
    1026      144813 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1027      144813 :   long i,j, l = lg(H), m = lgcols(H);
    1028      144813 :   int stable = 1;
    1029     1059216 :   for (j=1; j<l; j++)
    1030    11282405 :     for (i=1; i<m; i++)
    1031             :     {
    1032    10368002 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
    1033    10368002 :       if (h)
    1034             :       {
    1035     2062213 :         if (cmpii(h,lim) > 0) h = subii(h,qp);
    1036     2062213 :         gcoeff(H,i,j) = h; stable = 0;
    1037             :       }
    1038             :     }
    1039      144813 :   *ptq = qp; return stable;
    1040             : }
    1041             : 
    1042             : GEN
    1043        1295 : ZVM_init_CRT(GEN Hp, ulong p)
    1044             : {
    1045             :   long i, j, k;
    1046             :   GEN c, cp, d, dp, H;
    1047        1295 :   long m, l = lg(Hp), lim = (long)(p>>1), n;
    1048        1295 :   H = cgetg(l, t_MAT);
    1049        1295 :   if (l==1) return H;
    1050        1295 :   m = lgcols(Hp);
    1051        1295 :   n = lg(gmael(Hp,1,1));
    1052        5572 :   for (j=1; j<l; j++)
    1053             :   {
    1054        4277 :     cp = gel(Hp,j);
    1055        4277 :     c = cgetg(m, t_COL);
    1056        4277 :     gel(H,j) = c;
    1057       65240 :     for (i=1; i<m; i++)
    1058             :     {
    1059       60963 :       dp = gel(cp, i);
    1060       60963 :       d = cgetg(n, t_VEC);
    1061       60963 :       gel(c, i) = d;
    1062      183631 :       for (k=1; k<n; k++)
    1063      122668 :         gel(d,k) = stoi(Fl_center(dp[k], p, lim));
    1064             :     }
    1065             :   }
    1066        1295 :   return H;
    1067             : }
    1068             : 
    1069             : int
    1070         572 : ZVM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1071             : {
    1072         572 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), lim = shifti(qp,-1);
    1073         572 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1074         572 :   long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
    1075         572 :   int stable = 1;
    1076        3773 :   for (j=1; j<l; j++)
    1077       74697 :     for (i=1; i<m; i++)
    1078      214488 :       for (k=1; k<n; k++)
    1079             :       {
    1080      142992 :         h = Fl_chinese_coprime(gmael3(H,j,i,k), mael3(Hp,j,i,k),q,p,qinv,qp);
    1081      142992 :         if (h)
    1082             :         {
    1083      105188 :           if (cmpii(h,lim) > 0) h = subii(h,qp);
    1084      105188 :           gmael3(H,j,i,k) = h; stable = 0;
    1085             :         }
    1086             :       }
    1087         572 :   *ptq = qp; return stable;
    1088             : }
    1089             : 
    1090             : /* record the degrees of Euclidean remainders (make them as large as
    1091             :  * possible : smaller values correspond to a degenerate sequence) */
    1092             : static void
    1093        1561 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1094             : {
    1095             :   long da,db,dc, ind;
    1096        1561 :   pari_sp av = avma;
    1097             : 
    1098        1561 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1099        1561 :   da = degpol(a);
    1100        1561 :   db = degpol(b);
    1101        1561 :   if (db > da)
    1102           0 :   { swapspec(a,b, da,db); }
    1103        1561 :   else if (!da) return;
    1104        1561 :   ind = 0;
    1105        9814 :   while (db)
    1106             :   {
    1107        6692 :     GEN c = Flx_rem(a,b, p);
    1108        6692 :     a = b; b = c; dc = degpol(c);
    1109        6692 :     if (dc < 0) break;
    1110             : 
    1111        6692 :     ind++;
    1112        6692 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1113        6692 :     if (gc_needed(av,2))
    1114             :     {
    1115           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1116           0 :       gerepileall(av, 2, &a,&b);
    1117             :     }
    1118        6692 :     db = dc; /* = degpol(b) */
    1119             :   }
    1120        1561 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1121        1561 :   avma = av; return;
    1122             : }
    1123             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1124             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1125             :  * resultant(a,b). Modular version of Collins's subresultant */
    1126             : static ulong
    1127        6996 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1128             : {
    1129             :   long da,db,dc, ind;
    1130        6996 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1131        6996 :   int s = 1;
    1132        6996 :   pari_sp av = avma;
    1133             : 
    1134        6996 :   *C0 = 1; *C1 = 0;
    1135        6996 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1136        6996 :   da = degpol(a);
    1137        6996 :   db = degpol(b);
    1138        6996 :   if (db > da)
    1139             :   {
    1140           0 :     swapspec(a,b, da,db);
    1141           0 :     if (both_odd(da,db)) s = -s;
    1142             :   }
    1143        6996 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1144        6996 :   ind = 0;
    1145       40727 :   while (db)
    1146             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1147             :      * da = deg a, db = deg b */
    1148       27099 :     GEN c = Flx_rem(a,b, p);
    1149       27099 :     long delta = da - db;
    1150             : 
    1151       27099 :     if (both_odd(da,db)) s = -s;
    1152       27099 :     lb = Fl_mul(b[db+2], cb, p);
    1153       27099 :     a = b; b = c; dc = degpol(c);
    1154       27099 :     ind++;
    1155       27099 :     if (dc != dglist[ind]) { avma = av; return 0; } /* degenerates */
    1156       26735 :     if (g == h)
    1157             :     { /* frequent */
    1158       24621 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1159       24621 :       ca = cb;
    1160       24621 :       cb = cc;
    1161             :     }
    1162             :     else
    1163             :     {
    1164        2114 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1165        2114 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1166        2114 :       ca = cb;
    1167        2114 :       cb = Fl_div(cc, ghdelta, p);
    1168             :     }
    1169       26735 :     da = db; /* = degpol(a) */
    1170       26735 :     db = dc; /* = degpol(b) */
    1171             : 
    1172       26735 :     g = lb;
    1173       26735 :     if (delta == 1)
    1174       17482 :       h = g; /* frequent */
    1175             :     else
    1176        9253 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1177             : 
    1178       26735 :     if (gc_needed(av,2))
    1179             :     {
    1180           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1181           0 :       gerepileall(av, 2, &a,&b);
    1182             :     }
    1183             :   }
    1184        6632 :   if (da > 1) return 0; /* Failure */
    1185             :   /* last non-constant polynomial has degree 1 */
    1186        6632 :   *C0 = Fl_mul(ca, a[2], p);
    1187        6632 :   *C1 = Fl_mul(ca, a[3], p);
    1188        6632 :   res = Fl_mul(cb, b[2], p);
    1189        6632 :   if (s == -1) res = p - res;
    1190        6632 :   avma = av; return res;
    1191             : }
    1192             : 
    1193             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1194             :  * Return 0 in case of degree drop. */
    1195             : static GEN
    1196        8557 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1197             : {
    1198             :   GEN z;
    1199        8557 :   long i, lb = lg(Q);
    1200        8557 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1201        8557 :   long vs=mael(Q,2,1);
    1202        8557 :   if (!leadz) return zero_Flx(vs);
    1203             : 
    1204        8557 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1205        8557 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1206        8557 :   z[i] = leadz; return z;
    1207             : }
    1208             : 
    1209             : GEN
    1210       17836 : FpXY_Fq_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1211             : {
    1212       17836 :   pari_sp av = avma;
    1213       17836 :   long i, lb = lg(Q);
    1214             :   GEN z;
    1215       17836 :   if (!T) return FpXY_evaly(Q, y, p, vx);
    1216        1148 :   if (lb == 2) return pol_0(vx);
    1217        1148 :   z = gel(Q, lb-1);
    1218        1148 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1219             : 
    1220        1148 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1221       28084 :   for (i=lb-2; i>=2; i--)
    1222             :   {
    1223       26936 :     GEN c = gel(Q,i);
    1224       26936 :     z = FqX_Fq_mul(z, y, T, p);
    1225       26936 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1226             :   }
    1227        1148 :   return gerepileupto(av, z);
    1228             : }
    1229             : 
    1230             : static GEN
    1231        5586 : ZX_norml1(GEN x)
    1232             : {
    1233        5586 :   long i, l = lg(x);
    1234             :   GEN s;
    1235             : 
    1236        5586 :   if (l == 2) return gen_0;
    1237        5586 :   s = gel(x, l-1); /* != 0 */
    1238       30394 :   for (i = l-2; i > 1; i--) {
    1239       24808 :     GEN xi = gel(x,i);
    1240       24808 :     if (!signe(x)) continue;
    1241       24808 :     s = addii_sign(s,1, xi,1);
    1242             :   }
    1243        5586 :   return s;
    1244             : }
    1245             : 
    1246             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1247             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1248             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1249             :  * Return e such that Res(A, B) < 2^e */
    1250             : ulong
    1251       56740 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1252             : {
    1253       56740 :   pari_sp av = avma;
    1254       56740 :   GEN a = gen_0, b = gen_0;
    1255       56740 :   long i , lA = lg(A), lB = lg(B);
    1256             :   double loga, logb;
    1257       56740 :   for (i=2; i<lA; i++) a = addii(a, sqri(gel(A,i)));
    1258      418961 :   for (i=2; i<lB; i++)
    1259             :   {
    1260      362221 :     GEN t = gel(B,i);
    1261      362221 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1262      362221 :     b = addii(b, sqri(t));
    1263             :   }
    1264       56740 :   loga = dbllog2(a);
    1265       56740 :   logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
    1266       56740 :   i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1267       56740 :   avma = av; return (i <= 0)? 1: 1 + (ulong)i;
    1268             : }
    1269             : 
    1270             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1271             : static ulong
    1272      189109 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
    1273             : {
    1274      189109 :   GEN ev = FlxY_evalx(b, n, p);
    1275      189109 :   long drop = lg(b) - lg(ev);
    1276      189109 :   ulong r = Flx_resultant(a, ev, p);
    1277      189109 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
    1278      189109 :   return r;
    1279             : }
    1280             : static GEN
    1281           4 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1282             : {
    1283           4 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1284           4 :   long drop = db-degpol(ev);
    1285           4 :   GEN r = FpX_resultant(a, ev, p);
    1286           4 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1287           4 :   return r;
    1288             : }
    1289             : 
    1290             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1291             : /* Return a Fly */
    1292             : static GEN
    1293        4462 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, ulong dres, long sx)
    1294             : {
    1295        4462 :   ulong i, n, la = Flx_lead(a);
    1296        4462 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1297        4462 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1298             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1299             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1300       97014 :   for (i=0,n = 1; i < dres; n++)
    1301             :   {
    1302       92552 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1303       92552 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1304             :   }
    1305        4462 :   if (i == dres)
    1306             :   {
    1307        4005 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
    1308             :   }
    1309        4462 :   return Flv_polint(x,y, p, sx);
    1310             : }
    1311             : 
    1312             : static GEN
    1313        9143 : FlxX_pseudorem(GEN x, GEN y, ulong p)
    1314             : {
    1315        9143 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1316        9143 :   pari_sp av = avma, av2;
    1317             : 
    1318        9143 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1319        9143 :   (void)new_chunk(2);
    1320        9147 :   dx=degpol(x); x = RgX_recip_shallow(x)+2;
    1321        9154 :   dy=degpol(y); y = RgX_recip_shallow(y)+2; dz=dx-dy; dp = dz+1;
    1322        9160 :   av2 = avma;
    1323             :   for (;;)
    1324             :   {
    1325       70767 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1326      272108 :     for (i=1; i<=dy; i++)
    1327      400406 :       gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
    1328      200203 :                               Flx_mul(gel(x,0), gel(y,i), p), p );
    1329     1251751 :     for (   ; i<=dx; i++)
    1330     1181007 :       gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
    1331       75003 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1332       70744 :     if (dx < dy) break;
    1333       61601 :     if (gc_needed(av2,1))
    1334             :     {
    1335           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1336           0 :       gerepilecoeffs(av2,x,dx+1);
    1337             :     }
    1338       61607 :   }
    1339        9143 :   if (dx < 0) return zero_Flx(0);
    1340        9143 :   lx = dx+3; x -= 2;
    1341        9143 :   x[0]=evaltyp(t_POL) | evallg(lx);
    1342        9148 :   x[1]=evalsigne(1) | evalvarn(vx);
    1343        9148 :   x = RgX_recip_shallow(x);
    1344        9150 :   if (dp)
    1345             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1346        2152 :     GEN t = Flx_powu(gel(y,0), dp, p);
    1347        8605 :     for (i=2; i<lx; i++)
    1348        6455 :       gel(x,i) = Flx_mul(gel(x,i), t, p);
    1349             :   }
    1350        9148 :   return gerepilecopy(av, x);
    1351             : }
    1352             : 
    1353             : /* return a Flx */
    1354             : GEN
    1355        2998 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
    1356             : {
    1357        2998 :   pari_sp av = avma, av2;
    1358             :   long degq,dx,dy,du,dv,dr,signh;
    1359             :   GEN z,g,h,r,p1;
    1360             : 
    1361        2998 :   dx=degpol(u); dy=degpol(v); signh=1;
    1362        3002 :   if (dx < dy)
    1363             :   {
    1364          28 :     swap(u,v); lswap(dx,dy);
    1365          28 :     if (both_odd(dx, dy)) signh = -signh;
    1366             :   }
    1367        3002 :   if (dy < 0) return zero_Flx(sx);
    1368        3002 :   if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
    1369             : 
    1370        3002 :   g = h = pol1_Flx(sx); av2 = avma;
    1371             :   for(;;)
    1372             :   {
    1373        9147 :     r = FlxX_pseudorem(u,v,p); dr = lg(r);
    1374        9145 :     if (dr == 2) { avma = av; return zero_Flx(sx); }
    1375        9145 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1376        9151 :     u = v; p1 = g; g = leading_coeff(u);
    1377        9157 :     switch(degq)
    1378             :     {
    1379           0 :       case 0: break;
    1380             :       case 1:
    1381        6745 :         p1 = Flx_mul(h,p1, p); h = g; break;
    1382             :       default:
    1383        2412 :         p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
    1384        2409 :         h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
    1385             :     }
    1386        9148 :     if (both_odd(du,dv)) signh = -signh;
    1387        9148 :     v = FlxY_Flx_div(r, p1, p);
    1388        9151 :     if (dr==3) break;
    1389        6145 :     if (gc_needed(av2,1))
    1390             :     {
    1391           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, dr = %ld",dr);
    1392           0 :       gerepileall(av2,4, &u, &v, &g, &h);
    1393             :     }
    1394        6143 :   }
    1395        3006 :   z = gel(v,2);
    1396        3006 :   if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
    1397        3006 :   if (signh < 0) z = Flx_neg(z,p);
    1398        3006 :   return gerepileupto(av, z);
    1399             : }
    1400             : 
    1401             : /* Warning:
    1402             :  * This function switches between valid and invalid variable ordering*/
    1403             : 
    1404             : static GEN
    1405        3106 : FlxY_to_FlyX(GEN b, long sv)
    1406             : {
    1407        3106 :   long i, n=-1;
    1408        3106 :   long sw = b[1]&VARNBITS;
    1409        3106 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1410        3107 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1411             : }
    1412             : 
    1413             : /* Return a Fly*/
    1414             : GEN
    1415        3101 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
    1416             : {
    1417        3101 :   pari_sp ltop=avma;
    1418        3101 :   long dres = degpol(a)*degpol(b);
    1419        3102 :   long sx=a[1], sy=b[1]&VARNBITS;
    1420             :   GEN z;
    1421        3102 :   b = FlxY_to_FlyX(b,sx);
    1422        3105 :   if ((ulong)dres >= pp)
    1423        3001 :     z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
    1424             :   else
    1425         104 :     z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
    1426        3111 :   return gerepileupto(ltop,z);
    1427             : }
    1428             : 
    1429             : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
    1430             :  * in variable v. This is an incorrect PARI object if initially varn(b) << v.
    1431             :  * We could return a vector of coeffs, but it is convenient to have degpol()
    1432             :  * and friends available. Even in that case, it will behave nicely with all
    1433             :  * functions treating a polynomial as a vector of coeffs (eg poleval).
    1434             :  * FOR INTERNAL USE! */
    1435             : GEN
    1436        5369 : swap_vars(GEN b0, long v)
    1437             : {
    1438        5369 :   long i, n = RgX_degree(b0, v);
    1439             :   GEN b, x;
    1440        5369 :   if (n < 0) return pol_0(v);
    1441        5369 :   b = cgetg(n+3, t_POL); x = b + 2;
    1442        5369 :   b[1] = evalsigne(1) | evalvarn(v);
    1443        5369 :   for (i=0; i<=n; i++) gel(x,i) = polcoeff_i(b0, i, v);
    1444        5369 :   return b;
    1445             : }
    1446             : 
    1447             : /* assume varn(b) << varn(a) */
    1448             : /* return a FpY*/
    1449             : GEN
    1450        3073 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    1451             : {
    1452        3073 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    1453             :   GEN la,x,y;
    1454             : 
    1455        3073 :   if (lgefint(p) == 3)
    1456             :   {
    1457        3072 :     ulong pp = uel(p,2);
    1458        3072 :     b = ZXX_to_FlxX(b, pp, vX);
    1459        3073 :     a = ZX_to_Flx(a, pp);
    1460        3074 :     x = Flx_FlxY_resultant(a, b, pp);
    1461        3082 :     return Flx_to_ZX(x);
    1462             :   }
    1463           1 :   db = RgXY_degreex(b);
    1464           1 :   dres = degpol(a)*degpol(b);
    1465           1 :   la = leading_coeff(a);
    1466           1 :   x = cgetg(dres+2, t_VEC);
    1467           1 :   y = cgetg(dres+2, t_VEC);
    1468             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1469             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1470           3 :   for (i=0,n = 1; i < dres; n++)
    1471             :   {
    1472           2 :     gel(x,++i) = utoipos(n);
    1473           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1474           2 :     gel(x,++i) = subiu(p,n);
    1475           2 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1476             :   }
    1477           1 :   if (i == dres)
    1478             :   {
    1479           0 :     gel(x,++i) = gen_0;
    1480           0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    1481             :   }
    1482           1 :   return FpV_polint(x,y, p, vY);
    1483             : }
    1484             : 
    1485             : GEN
    1486         441 : FpX_direct_compositum(GEN a, GEN b, GEN p)
    1487             : {
    1488         441 :   GEN x = deg1pol_shallow(gen_1, pol_x(varn(a)), fetch_var_higher()); /* x+y */
    1489         441 :   x = FpX_FpXY_resultant(a, poleval(b,x),p);
    1490         441 :   (void)delete_var(); return x;
    1491             : }
    1492             : 
    1493             : /* 0, 1, -1, 2, -2, ... */
    1494             : #define next_lambda(a) (a>0 ? -a : 1-a)
    1495             : GEN
    1496           0 : FpX_compositum(GEN a, GEN b, GEN p)
    1497             : {
    1498           0 :   long k, v = fetch_var_higher();
    1499           0 :   for (k = 1;; k = next_lambda(k))
    1500             :   {
    1501           0 :     GEN x = deg1pol_shallow(gen_1, gmulsg(k, pol_x(v)), 0); /* x + k y */
    1502           0 :     GEN C = FpX_FpXY_resultant(a, poleval(b,x),p);
    1503           0 :     if (FpX_is_squarefree(C, p)) { (void)delete_var(); return C; }
    1504           0 :   }
    1505             : }
    1506             : 
    1507             : /* Assume A in Z[Y], B in Q[Y][X], and Res_Y(A, B) in Z[X].
    1508             :  * If lambda = NULL, return Res_Y(A,B).
    1509             :  * Otherwise, find a small lambda (start from *lambda, use the sequence above)
    1510             :  * such that R(X) = Res_Y(A(Y), B(X + lambda Y)) is squarefree, reset *lambda
    1511             :  * to the chosen value and return R
    1512             :  *
    1513             :  * If LERS is non-NULL, set it to the Last non-constant polynomial in the
    1514             :  * Euclidean Remainder Sequence */
    1515             : GEN
    1516        4172 : ZX_ZXY_resultant_all(GEN A, GEN B0, long *plambda, GEN *LERS)
    1517             : {
    1518        4172 :   int checksqfree = plambda? 1: 0, stable;
    1519        4172 :   long lambda = plambda? *plambda: 0, cnt = 0;
    1520             :   ulong bound, dp;
    1521        4172 :   pari_sp av = avma, av2 = 0;
    1522        4172 :   long i,n, lb, degA = degpol(A), dres = degA*degpol(B0);
    1523        4172 :   long v = fetch_var_higher();
    1524        4172 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    1525        4172 :   long sX = evalvarn(vX);
    1526             :   GEN x, y, dglist, dB, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    1527             :   forprime_t S;
    1528             : 
    1529        4172 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    1530        4172 :   if (LERS)
    1531             :   {
    1532        1491 :     if (!checksqfree)
    1533           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    1534        1491 :     C0 = cgetg(dres+2, t_VECSMALL);
    1535        1491 :     C1 = cgetg(dres+2, t_VECSMALL);
    1536        1491 :     dglist = cgetg(dres+1, t_VECSMALL);
    1537             :   }
    1538        4172 :   x = cgetg(dres+2, t_VECSMALL);
    1539        4172 :   y = cgetg(dres+2, t_VECSMALL);
    1540        4172 :   B0 = Q_remove_denom(B0, &dB);
    1541        4172 :   if (!dB) B0 = leafcopy(B0);
    1542        4172 :   A = leafcopy(A);
    1543        4172 :   B = B0;
    1544        4172 :   setvarn(A,v);
    1545             :   /* make sure p large enough */
    1546             : INIT:
    1547             :   /* always except the first time */
    1548        5334 :   if (av2) { avma = av2; lambda = next_lambda(lambda); }
    1549        5334 :   if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
    1550        5334 :   B = swap_vars(B, vY); setvarn(B,v);
    1551             :   /* B0(lambda v + x, v) */
    1552        5334 :   if (DEBUGLEVEL>4 && checksqfree) err_printf("Trying lambda = %ld\n", lambda);
    1553        5334 :   av2 = avma;
    1554             : 
    1555        5334 :   if (degA <= 3)
    1556             :   { /* sub-resultant faster for small degrees */
    1557        3724 :     if (LERS)
    1558             :     { /* implies checksqfree */
    1559        1974 :       H = resultant_all(A,B,&q);
    1560        1974 :       if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    1561        1344 :       H0 = gel(q,2);
    1562        1344 :       if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    1563        1344 :       H1 = gel(q,3);
    1564        1344 :       if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    1565             :     }
    1566             :     else
    1567        1750 :       H = resultant(A,B);
    1568        3094 :     if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
    1569        2961 :     goto END;
    1570             :   }
    1571             : 
    1572        1610 :   H = H0 = H1 = NULL;
    1573        1610 :   lb = lg(B);
    1574        1610 :   bound = ZX_ZXY_ResBound(A, B, dB);
    1575        1610 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    1576        1610 :   dp = 1;
    1577        1610 :   init_modular_big(&S);
    1578             :   while (1)
    1579             :   {
    1580        4755 :     ulong p = u_forprime_next(&S);
    1581        4755 :     if (dB) { dp = umodiu(dB, p); if (!dp) continue; }
    1582             : 
    1583        4755 :     a = ZX_to_Flx(A, p);
    1584        4755 :     b = ZXX_to_FlxX(B, p, varn(A));
    1585        4755 :     if (LERS)
    1586             :     {
    1587             :       GEN Hi;
    1588         397 :       if (degpol(a) < degA || lg(b) < lb) continue; /* p | lc(A)lc(B) */
    1589         397 :       if (checksqfree)
    1590             :       { /* find degree list for generic Euclidean Remainder Sequence */
    1591         196 :         long goal = minss(degpol(a), degpol(b)); /* longest possible */
    1592         196 :         for (n=1; n <= goal; n++) dglist[n] = 0;
    1593         196 :         setlg(dglist, 1);
    1594        1645 :         for (n=0; n <= dres; n++)
    1595             :         {
    1596        1561 :           ev = FlxY_evalx_drop(b, n, p);
    1597        1561 :           Flx_resultant_set_dglist(a, ev, dglist, p);
    1598        1561 :           if (lg(dglist)-1 == goal) break;
    1599             :         }
    1600             :         /* last pol in ERS has degree > 1 ? */
    1601         196 :         goal = lg(dglist)-1;
    1602         196 :         if (degpol(B) == 1) { if (!goal) goto INIT; }
    1603             :         else
    1604             :         {
    1605         189 :           if (goal <= 1) goto INIT;
    1606         182 :           if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    1607             :         }
    1608         189 :         if (DEBUGLEVEL>4)
    1609           0 :           err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    1610             :       }
    1611             : 
    1612        7386 :       for (i=0,n = 0; i <= dres; n++)
    1613             :       {
    1614        6996 :         ev = FlxY_evalx_drop(b, n, p);
    1615        6996 :         x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    1616        6996 :         if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    1617             :       }
    1618         390 :       Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    1619         390 :       Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    1620             :     }
    1621             :     else
    1622             :     {
    1623        4358 :       long dropa = degA - degpol(a), dropb = lb - lg(b);
    1624        4358 :       Hp = Flx_FlxY_resultant_polint(a, b, p, (ulong)dres, sX);
    1625        4358 :       if (dropa && dropb)
    1626           0 :         Hp = zero_Flx(sX);
    1627             :       else {
    1628        4358 :         if (dropa)
    1629             :         { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1630           0 :           GEN c = gel(b,lb-1); /* lc(B) */
    1631           0 :           if (!odd(lb)) c = Flx_neg(c, p); /* deg B = lb - 3 */
    1632           0 :           if (!Flx_equal1(c)) {
    1633           0 :             c = Flx_powu(c, dropa, p);
    1634           0 :             if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
    1635             :           }
    1636             :         }
    1637        4358 :         else if (dropb)
    1638             :         { /* multiply by lc(A)^(deg B - deg b) */
    1639           0 :           ulong c = a[degA+2]; /* lc(A) */
    1640           0 :           c = Fl_powu(c, dropb, p);
    1641           0 :           if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
    1642             :         }
    1643             :       }
    1644             :     }
    1645        4748 :     if (!H && degpol(Hp) != dres) continue;
    1646        4748 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1647        4748 :     if (checksqfree) {
    1648        1050 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    1649         658 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    1650         658 :       checksqfree = 0;
    1651             :     }
    1652             : 
    1653        4356 :     if (!H)
    1654             :     { /* initialize */
    1655        1211 :       q = utoipos(p); stable = 0;
    1656        1211 :       H = ZX_init_CRT(Hp, p,vX);
    1657        1211 :       if (LERS) {
    1658         189 :         H0= ZX_init_CRT(H0p, p,vX);
    1659         189 :         H1= ZX_init_CRT(H1p, p,vX);
    1660             :       }
    1661             :     }
    1662             :     else
    1663             :     {
    1664        3145 :       if (LERS) {
    1665         201 :         GEN qp = muliu(q,p);
    1666         402 :         stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    1667         201 :                 & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    1668         201 :                 & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    1669         201 :         q = qp;
    1670             :       }
    1671             :       else
    1672        2944 :         stable = ZX_incremental_CRT(&H, Hp, &q, p);
    1673             :     }
    1674             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    1675             :      * Probabilistic anyway for H0, H1 */
    1676        4356 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    1677           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    1678        4356 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    1679        3145 :     if (gc_needed(av,2))
    1680             :     {
    1681           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    1682           0 :       gerepileall(av2, LERS? 4: 2, &H, &q, &H0, &H1);
    1683             :     }
    1684        3145 :   }
    1685             : END:
    1686        4172 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    1687        4172 :   setvarn(H, vX); (void)delete_var();
    1688        4172 :   if (plambda) *plambda = lambda;
    1689        4172 :   if (LERS)
    1690             :   {
    1691        1491 :     *LERS = mkvec2(H0,H1);
    1692        1491 :     gerepileall(av, 2, &H, LERS);
    1693        1491 :     return H;
    1694             :   }
    1695        2681 :   return gerepilecopy(av, H);
    1696             : }
    1697             : 
    1698             : GEN
    1699        1988 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    1700             : {
    1701        1988 :   return ZX_ZXY_resultant_all(A, B, lambda, NULL);
    1702             : }
    1703             : 
    1704             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    1705             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    1706             :  * squarefree */
    1707             : GEN
    1708        1862 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    1709             : {
    1710        1862 :   pari_sp av = avma;
    1711             :   GEN R, a;
    1712             :   long dA;
    1713             :   int delvar;
    1714             : 
    1715        1862 :   if (v < 0) v = 0;
    1716        1862 :   switch (typ(A))
    1717             :   {
    1718        1862 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    1719           0 :       A = constant_coeff(A);
    1720             :     default:
    1721           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    1722           0 :       return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    1723             :   }
    1724        1862 :   delvar = 0;
    1725        1862 :   if (varn(T) == 0)
    1726             :   {
    1727        1806 :     long v0 = fetch_var(); delvar = 1;
    1728        1806 :     T = leafcopy(T); setvarn(T,v0);
    1729        1806 :     A = leafcopy(A); setvarn(A,v0);
    1730             :   }
    1731        1862 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    1732        1862 :   if (delvar) (void)delete_var();
    1733        1862 :   setvarn(R, v); a = leading_coeff(T);
    1734        1862 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    1735        1862 :   return gerepileupto(av, R);
    1736             : }
    1737             : 
    1738             : 
    1739             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    1740             : GEN
    1741       11486 : ZXQ_charpoly(GEN A, GEN T, long v)
    1742             : {
    1743       11486 :   return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    1744             : }
    1745             : 
    1746             : GEN
    1747         819 : QXQ_charpoly(GEN A, GEN T, long v)
    1748             : {
    1749         819 :   pari_sp av = avma;
    1750         819 :   GEN den, B = Q_remove_denom(A, &den);
    1751         819 :   GEN P = ZXQ_charpoly(B, T, v);
    1752         819 :   return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
    1753             : }
    1754             : 
    1755             : static GEN
    1756      120361 : trivial_case(GEN A, GEN B)
    1757             : {
    1758             :   long d;
    1759      120361 :   if (typ(A) == t_INT) return powiu(A, degpol(B));
    1760      117504 :   d = degpol(A);
    1761      117504 :   if (d == 0) return trivial_case(gel(A,2),B);
    1762      116299 :   if (d < 0) return gen_0;
    1763      116284 :   return NULL;
    1764             : }
    1765             : 
    1766             : static long
    1767       57567 : get_nbprimes(ulong bound, ulong *pt_start)
    1768             : {
    1769             : #ifdef LONG_IS_64BIT
    1770       49218 :   ulong pstart = 4611686018427388039UL;
    1771             : #else
    1772        8349 :   ulong pstart = 1073741827UL;
    1773             : #endif
    1774       57567 :   *pt_start = pstart;
    1775       57567 :   return (bound/expu(pstart))+1;
    1776             : }
    1777             : 
    1778             : static ulong
    1779      331261 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    1780             : {
    1781      331261 :   pari_sp av = avma;
    1782             :   ulong H;
    1783             :   long dropa, dropb;
    1784      331261 :   ulong dp = dB ? umodiu(dB, p): 1;
    1785      331288 :   if (!b) b = Flx_deriv(a, p);
    1786      331300 :   dropa = degA - degpol(a);
    1787      331301 :   dropb = degB - degpol(b);
    1788      331304 :   if (dropa && dropb) /* p | lc(A), p | lc(B) */
    1789           0 :   { avma = av; return 0; }
    1790      331304 :   H = Flx_resultant(a, b, p);
    1791      331282 :   if (dropa)
    1792             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1793           0 :     ulong c = b[degB+2]; /* lc(B) */
    1794           0 :     if (odd(degB)) c = p - c;
    1795           0 :     c = Fl_powu(c, dropa, p);
    1796           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1797             :   }
    1798      331282 :   else if (dropb)
    1799             :   { /* multiply by lc(A)^(deg B - deg b) */
    1800           0 :     ulong c = a[degA+2]; /* lc(A) */
    1801           0 :     c = Fl_powu(c, dropb, p);
    1802           0 :     if (c != 1) H = Fl_mul(H, c, p);
    1803             :   }
    1804      331277 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    1805      331277 :   avma = av; return H;
    1806             : }
    1807             : 
    1808             : /* If B=NULL, assume B=A' */
    1809             : static GEN
    1810      217927 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    1811             : {
    1812      217927 :   pari_sp av = avma;
    1813      217927 :   long degA, degB, i, n = lg(P)-1;
    1814             :   GEN H, T;
    1815             : 
    1816      217927 :   degA = degpol(A);
    1817      217938 :   degB = B ? degpol(B): degA - 1;
    1818      217974 :   if (n == 1)
    1819             :   {
    1820      178238 :     ulong Hp, p = uel(P,1);
    1821             :     GEN a, b;
    1822      178238 :     a = ZX_to_Flx(A, p), b = B ? ZX_to_Flx(B, p): NULL;
    1823      178230 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1824      178241 :     avma = av;
    1825      178241 :     *mod = utoi(p); return utoi(Hp);
    1826             :   }
    1827       39736 :   T = ZV_producttree(P);
    1828       39736 :   A = ZX_nv_mod_tree(A, P, T);
    1829       39739 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    1830       39739 :   H = cgetg(n+1, t_VECSMALL);
    1831      192783 :   for(i=1; i <= n; i++)
    1832             :   {
    1833      153044 :     ulong p = P[i];
    1834      153044 :     GEN a = gel(A, i), b = B ? gel(B, i): NULL;
    1835      153044 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    1836             :   }
    1837       39739 :   H = ZV_chinese_tree(H, P, T, mod);
    1838       39737 :   gerepileall(av, 2, &H, mod);
    1839       39737 :   return H;
    1840             : }
    1841             : 
    1842             : GEN
    1843      178494 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    1844             : {
    1845      178494 :   GEN V = cgetg(3, t_VEC);
    1846      178531 :   if (isintzero(B)) B = NULL;
    1847      178510 :   if (isintzero(dB)) dB = NULL;
    1848      178526 :   gel(V,1) = ZX_resultant_slice(A,B,dB,P,&gel(V,2));
    1849      178461 :   return V;
    1850             : }
    1851             : 
    1852             : static GEN
    1853      218059 : primelist_disc(ulong *p, long n, GEN dB)
    1854             : {
    1855      218059 :   GEN P = cgetg(n+1, t_VECSMALL);
    1856             :   long i;
    1857      549430 :   for (i=1; i <= n; i++, *p = unextprime(*p+1))
    1858             :   {
    1859      331371 :     if (dB && umodiu(dB, *p)==0) { i--; continue; }
    1860      331371 :     P[i] = *p;
    1861             :   }
    1862      218059 :   return P;
    1863             : }
    1864             : 
    1865             : /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
    1866             : /* if B=NULL, take B = A' */
    1867             : GEN
    1868       60439 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    1869             : {
    1870             :   ulong p;
    1871       60439 :   pari_sp av = avma;
    1872             :   long n, m;
    1873             :   GEN  H, P, mod;
    1874       60439 :   int is_disc = !B;
    1875       60439 :   if (is_disc) B = ZX_deriv(A);
    1876             : 
    1877       60439 :   if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
    1878       57567 :   if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    1879       57567 :   n = get_nbprimes(bound+1, &p);/* +1 to account for sign */
    1880       57567 :   if (is_disc)
    1881       33442 :     B = NULL;
    1882       57567 :   m = minss(degpol(A)+(B ? degpol(B): 0), n);
    1883       57567 :   if (m == 1)
    1884             :   {
    1885       33242 :     GEN P = primelist_disc(&p, n, dB);
    1886       33242 :     H = ZX_resultant_slice(A, B, dB, P, &mod);
    1887             :   }
    1888             :   else
    1889             :   {
    1890       24325 :     long i, s = n/m, r = n - m*s, di = 0;
    1891       24325 :     GEN worker = strtoclosure("_ZX_resultant_worker", 3, A, B?B:gen_0, dB?dB:gen_0);
    1892             :     struct pari_mt pt;
    1893             :     long pending;
    1894       24325 :     if (DEBUGLEVEL > 4)
    1895           0 :       err_printf("ZX_resultant: bound 2^%ld, nb primes: %ld\n",bound, n);
    1896       24325 :     H = cgetg(m+1+!!r, t_VEC); P = cgetg(m+1+!!r, t_VEC);
    1897       24325 :     mt_queue_start_lim(&pt, worker, m);
    1898      219669 :     for (i=1; i<=m || pending; i++)
    1899             :     {
    1900             :       GEN done;
    1901      195344 :       mt_queue_submit(&pt, i, i<=m ? mkvec(primelist_disc(&p, s, dB)): NULL);
    1902      195344 :       done = mt_queue_get(&pt, NULL, &pending);
    1903      195344 :       if (done)
    1904             :       {
    1905      178591 :         di++;
    1906      178591 :         gel(H, di) = gel(done,1);
    1907      178591 :         gel(P, di) = gel(done,2);
    1908      178591 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
    1909             :       }
    1910             :     }
    1911       24325 :     mt_queue_end(&pt);
    1912       24325 :     if (r)
    1913             :     {
    1914        6226 :       GEN Pr = primelist_disc(&p, r, dB);
    1915        6226 :       gel(H, m+1) = ZX_resultant_slice(A, B, dB, Pr, &gel(P, m+1));
    1916             :     }
    1917       24325 :     H = ZV_chinese(H, P, &mod);
    1918       24325 :     if (DEBUGLEVEL>5) err_printf("done\n");
    1919             :   }
    1920       57567 :   H = Fp_center(H, mod, shifti(mod,-1));
    1921       57567 :   return gerepileuptoint(av, H);
    1922             : }
    1923             : 
    1924             : /* A0 and B0 in Q[X] */
    1925             : GEN
    1926       10530 : QX_resultant(GEN A0, GEN B0)
    1927             : {
    1928             :   GEN s, a, b, A, B;
    1929       10530 :   pari_sp av = avma;
    1930             : 
    1931       10530 :   A = Q_primitive_part(A0, &a);
    1932       10530 :   B = Q_primitive_part(B0, &b);
    1933       10530 :   s = ZX_resultant(A, B);
    1934       10530 :   if (!signe(s)) { avma = av; return gen_0; }
    1935       10530 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    1936       10530 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    1937       10530 :   return gerepileupto(av, s);
    1938             : }
    1939             : 
    1940             : GEN
    1941       26352 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    1942             : 
    1943             : GEN
    1944           0 : QXQ_intnorm(GEN A, GEN B)
    1945             : {
    1946             :   GEN c, n, R, lB;
    1947           0 :   long dA = degpol(A), dB = degpol(B);
    1948           0 :   pari_sp av = avma;
    1949           0 :   if (dA < 0) return gen_0;
    1950           0 :   A = Q_primitive_part(A, &c);
    1951           0 :   if (!c || typ(c) == t_INT) {
    1952           0 :     n = c;
    1953           0 :     R = ZX_resultant(B, A);
    1954             :   } else {
    1955           0 :     n = gel(c,1);
    1956           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    1957             :   }
    1958           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    1959           0 :   lB = leading_coeff(B);
    1960           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    1961           0 :   return gerepileuptoint(av, R);
    1962             : }
    1963             : 
    1964             : GEN
    1965           0 : QXQ_norm(GEN A, GEN B)
    1966             : {
    1967             :   GEN c, R, lB;
    1968           0 :   long dA = degpol(A), dB = degpol(B);
    1969           0 :   pari_sp av = avma;
    1970           0 :   if (dA < 0) return gen_0;
    1971           0 :   A = Q_primitive_part(A, &c);
    1972           0 :   R = ZX_resultant(B, A);
    1973           0 :   if (c) R = gmul(R, gpowgs(c, dB));
    1974           0 :   lB = leading_coeff(B);
    1975           0 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    1976           0 :   return gerepileupto(av, R);
    1977             : }
    1978             : 
    1979             : /* assume x has integral coefficients */
    1980             : GEN
    1981       34730 : ZX_disc_all(GEN x, ulong bound)
    1982             : {
    1983       34730 :   pari_sp av = avma;
    1984             :   GEN l, R;
    1985       34730 :   long s, d = degpol(x);
    1986       34730 :   if (d <= 1) return d ? gen_1: gen_0;
    1987       33442 :   s = (d & 2) ? -1: 1;
    1988       33442 :   l = leading_coeff(x);
    1989       33442 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    1990       33442 :   if (is_pm1(l))
    1991       30635 :   { if (signe(l) < 0) s = -s; }
    1992             :   else
    1993        2807 :     R = diviiexact(R,l);
    1994       33442 :   if (s == -1) togglesign_safe(&R);
    1995       33442 :   return gerepileuptoint(av,R);
    1996             : }
    1997       32286 : GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    1998             : 
    1999             : GEN
    2000           0 : QX_disc(GEN x)
    2001             : {
    2002           0 :   pari_sp av = avma;
    2003           0 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2004           0 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2005           0 :   return gerepileupto(av, d);
    2006             : }
    2007             : 
    2008             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2009             : GEN
    2010        9691 : QXQ_inv(GEN A, GEN B)
    2011             : {
    2012             :   GEN D, cU, q, U, V;
    2013             :   ulong p;
    2014        9691 :   pari_sp av2, av = avma;
    2015             :   forprime_t S;
    2016             :   pari_timer ti;
    2017        9691 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2018             :   /* A a QX, B a ZX */
    2019        9691 :   if (degpol(A) < 15) return RgXQ_inv(A,B);
    2020           7 :   A = Q_primitive_part(A, &D);
    2021             :   /* A, B in Z[X] */
    2022           7 :   init_modular_small(&S);
    2023           7 :   if (DEBUGLEVEL>5) timer_start(&ti);
    2024           7 :   av2 = avma; U = NULL;
    2025         578 :   while ((p = u_forprime_next(&S)))
    2026             :   {
    2027             :     GEN a, b, qp, Up, Vp;
    2028             :     int stable;
    2029             : 
    2030         571 :     a = ZX_to_Flx(A, p);
    2031         571 :     b = ZX_to_Flx(B, p);
    2032             :     /* if p | Res(A/G, B/G), discard */
    2033         578 :     if (!Flx_extresultant(b,a,p, &Vp,&Up)) continue;
    2034             : 
    2035         571 :     if (!U)
    2036             :     { /* First time */
    2037           7 :       U = ZX_init_CRT(Up,p,varn(A));
    2038           7 :       V = ZX_init_CRT(Vp,p,varn(A));
    2039           7 :       q = utoipos(p); continue;
    2040             :     }
    2041         564 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: mod %ld (bound 2^%ld)", p,expi(q));
    2042         564 :     qp = muliu(q,p);
    2043        1128 :     stable = ZX_incremental_CRT_raw(&U, Up, q,qp, p)
    2044         564 :            & ZX_incremental_CRT_raw(&V, Vp, q,qp, p);
    2045         564 :     if (stable)
    2046             :     { /* all stable: check divisibility */
    2047           7 :       GEN res = ZX_add(ZX_mul(A,U), ZX_mul(B,V));
    2048           7 :       if (degpol(res) == 0) {
    2049           7 :         res = gel(res,2);
    2050           7 :         D = D? gmul(D, res): res;
    2051          14 :         break;
    2052             :       } /* DONE */
    2053           0 :       if (DEBUGLEVEL) err_printf("QXQ_inv: char 0 check failed");
    2054             :     }
    2055         557 :     q = qp;
    2056         557 :     if (gc_needed(av,1))
    2057             :     {
    2058           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"QXQ_inv");
    2059           0 :       gerepileall(av2, 3, &q,&U,&V);
    2060             :     }
    2061             :   }
    2062           7 :   if (!p) pari_err_OVERFLOW("QXQ_inv [ran out of primes]");
    2063           7 :   cU = ZX_content(U);
    2064           7 :   if (!is_pm1(cU)) { U = Q_div_to_int(U, cU); D = gdiv(D, cU); }
    2065           7 :   return gerepileupto(av, RgX_Rg_div(U, D));
    2066             : }
    2067             : 
    2068             : /************************************************************************
    2069             :  *                                                                      *
    2070             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    2071             :  *                                                                      *
    2072             :  ************************************************************************/
    2073             : 
    2074             : /* irreducible (unitary) polynomial of degree n over Fp */
    2075             : GEN
    2076           0 : ffinit_rand(GEN p,long n)
    2077             : {
    2078             :   for(;;) {
    2079           0 :     pari_sp av = avma;
    2080           0 :     GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
    2081           0 :     if (FpX_is_irred(pol, p)) return pol;
    2082           0 :     avma = av;
    2083           0 :   }
    2084             : }
    2085             : 
    2086             : /* return an extension of degree 2^l of F_2, assume l > 0
    2087             :  * Not stack clean. */
    2088             : static GEN
    2089         587 : f2init(long l)
    2090             : {
    2091             :   GEN Q, T, S;
    2092             :   long i, v;
    2093             : 
    2094         587 :   if (l == 1) return polcyclo(3, 0);
    2095         552 :   v = fetch_var_higher();
    2096         550 :   S = mkpoln(4, gen_1,gen_1,gen_0,gen_0); /* y(y^2 + y) */
    2097         550 :   Q = mkpoln(3, gen_1,gen_1, S); /* x^2 + x + y(y^2+y) */
    2098         551 :   setvarn(Q, v);
    2099             : 
    2100             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    2101         551 :   T = mkpoln(5, gen_1,gen_0,gen_0,gen_1,gen_1);
    2102         551 :   setvarn(T, v);
    2103             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    2104             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    2105             :    * ==> x^2 + x + (b^2+b)b */
    2106         551 :   for (i=2; i<l; i++) T = FpX_FpXY_resultant(T, Q, gen_2); /* minpoly(b) / F2*/
    2107         554 :   (void)delete_var(); setvarn(T,0); return T;
    2108             : }
    2109             : 
    2110             : /* return an extension of degree p^l of F_p, assume l > 0
    2111             :  * Not stack clean. */
    2112             : GEN
    2113           0 : ffinit_Artin_Shreier(GEN ip, long l)
    2114             : {
    2115           0 :   long i, v, p = itos(ip);
    2116           0 :   GEN T, Q, xp = pol_xn(p,0); /* x^p */
    2117           0 :   T = ZX_sub(xp, deg1pol_shallow(gen_1,gen_1,0)); /* x^p - x - 1 */
    2118           0 :   if (l == 1) return T;
    2119             : 
    2120           0 :   v = fetch_var_higher();
    2121           0 :   setvarn(xp, v);
    2122           0 :   Q = ZX_sub(pol_xn(2*p-1,0), pol_xn(p,0));
    2123           0 :   Q = gsub(xp, deg1pol_shallow(gen_1, Q, v)); /* x^p - x - (y^(2p-1)-y^p) */
    2124           0 :   for (i = 2; i <= l; ++i) T = FpX_FpXY_resultant(T, Q, ip);
    2125           0 :   (void)delete_var(); setvarn(T,0); return T;
    2126             : }
    2127             : 
    2128             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    2129             : static long
    2130       12092 : fpinit_check(GEN p, long n, long l)
    2131             : {
    2132             :   ulong q;
    2133       12092 :   if (!uisprime(n)) return 0;
    2134        5873 :   q = umodiu(p,n); if (!q) return 0;
    2135        5306 :   return cgcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    2136             : }
    2137             : 
    2138             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    2139             :  * Return an irreducible polynomial of degree l over F_p.
    2140             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    2141             :  * finite fields", ACM, 1986 (5) 350--355.
    2142             :  * Not stack clean */
    2143             : static GEN
    2144        3010 : fpinit(GEN p, long l)
    2145             : {
    2146        3010 :   ulong n = 1+l;
    2147        3010 :   while (!fpinit_check(p,n,l)) n += l;
    2148        3010 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    2149        3010 :   return FpX_red(polsubcyclo(n,l,0),p);
    2150             : }
    2151             : 
    2152             : static GEN
    2153        3051 : ffinit_fact(GEN p, long n)
    2154             : {
    2155        3051 :   GEN P, F = gel(factoru_pow(n),3);
    2156             :   long i;
    2157        3052 :   if (!odd(n) && absequaliu(p, 2))
    2158         587 :     P = f2init(vals(n)); /* if n is even, F[1] = 2^vals(n)*/
    2159             :   else
    2160        2465 :     P = fpinit(p, F[1]);
    2161        3496 :   for (i = 2; i < lg(F); ++i)
    2162         441 :     P = FpX_direct_compositum(fpinit(p, F[i]), P, p);
    2163        3055 :   return P;
    2164             : }
    2165             : 
    2166             : static GEN
    2167         103 : ffinit_nofact(GEN p, long n)
    2168             : {
    2169         103 :   GEN P, Q = NULL;
    2170         103 :   if (lgefint(p)==3)
    2171             :   {
    2172           0 :     ulong pp = p[2], q;
    2173           0 :     long v = u_lvalrem(n,pp,&q);
    2174           0 :     if (v>0)
    2175             :     {
    2176           0 :       Q = (pp == 2)? f2init(v): fpinit(p,n/q);
    2177           0 :       n = q;
    2178             :     }
    2179             :   }
    2180             :   /* n coprime to p */
    2181         103 :   if (n==1) P = Q;
    2182             :   else
    2183             :   {
    2184         103 :     P = fpinit(p, n);
    2185         103 :     if (Q) P = FpX_direct_compositum(P, Q, p);
    2186             :   }
    2187         103 :   return P;
    2188             : }
    2189             : 
    2190             : static GEN
    2191        3945 : init_Fq_i(GEN p, long n, long v)
    2192             : {
    2193             :   GEN P;
    2194        3945 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    2195        3946 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    2196        3946 :   if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
    2197        3946 :   if (v < 0) v = 0;
    2198        3946 :   if (n == 1) return pol_x(v);
    2199        3757 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    2200        3155 :   if (lgefint(p)-2 <= expu(n))
    2201        3050 :     P = ffinit_fact(p,n);
    2202             :   else
    2203         103 :     P = ffinit_nofact(p,n);
    2204        3158 :   setvarn(P, v); return P;
    2205             : }
    2206             : GEN
    2207        3820 : init_Fq(GEN p, long n, long v)
    2208             : {
    2209        3820 :   pari_sp av = avma;
    2210        3820 :   return gerepileupto(av, init_Fq_i(p, n, v));
    2211             : }
    2212             : GEN
    2213         126 : ffinit(GEN p, long n, long v)
    2214             : {
    2215         126 :   pari_sp av = avma;
    2216         126 :   return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    2217             : }
    2218             : 
    2219             : GEN
    2220        3178 : ffnbirred(GEN p, long n)
    2221             : {
    2222        3178 :   pari_sp av = avma;
    2223             :   long j;
    2224        3178 :   GEN s = gen_0, dk, pd;
    2225        3178 :   dk = divisorsu(n);
    2226       10535 :   for (j = 1; j < lg(dk); ++j)
    2227             :   {
    2228        7357 :     long d = dk[j];
    2229        7357 :     long m = moebiusu(d);
    2230        7357 :     if (!m) continue;
    2231        6797 :     pd = powiu(p, n/d);
    2232        6797 :     s = m>0 ? addii(s, pd): subii(s,pd);
    2233             :   }
    2234        3178 :   return gerepileuptoint(av, divis(s, n));
    2235             : }
    2236             : 
    2237             : GEN
    2238         434 : ffsumnbirred(GEN p, long n)
    2239             : {
    2240         434 :   pari_sp av = avma;
    2241             :   long i,j;
    2242         434 :   GEN v,q, t = gen_0;
    2243         434 :   v = cgetg(n+1,t_VECSMALL); v[1] = 1;
    2244         434 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    2245        1547 :   for(i=2; i<=n; i++)
    2246             :   {
    2247        1113 :     v[i] = moebiusu(i);
    2248        1113 :     gel(q,i) = mulii(gel(q,i-1), p);
    2249             :   }
    2250        1981 :   for(i=1; i<=n; i++)
    2251             :   {
    2252        1547 :     GEN s = gen_0;
    2253        1547 :     GEN dk = divisorsu(i);
    2254        4725 :     for (j = 1; j < lg(dk); ++j)
    2255             :     {
    2256        3178 :       long d = dk[j], m = v[d];
    2257        3178 :       if (!m) continue;
    2258        2884 :       s = m>0 ? addii(s, gel(q, i/d)): subii(s, gel(q, i/d));
    2259             :     }
    2260        1547 :     t = addii(t, divis(s, i));
    2261             :   }
    2262         434 :   return gerepileuptoint(av, t);
    2263             : }
    2264             : 
    2265             : GEN
    2266         140 : ffnbirred0(GEN p, long n, long flag)
    2267             : {
    2268         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    2269         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    2270         140 :   switch(flag)
    2271             :   {
    2272             :     case 0:
    2273          70 :       return ffnbirred(p, n);
    2274             :     case 1:
    2275          70 :       return ffsumnbirred(p, n);
    2276             :     default:
    2277           0 :       pari_err_FLAG("ffnbirred");
    2278             :   }
    2279             :   return NULL; /* LCOV_EXCL_LINE */
    2280             : }

Generated by: LCOV version 1.11