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 to exceed 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.18.1 lcov report (development 30401-87099e65a6) Lines: 1903 2156 88.3 %
Date: 2025-07-18 09:23:12 Functions: 205 225 91.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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /***********************************************************************/
      16             : /**                                                                   **/
      17             : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      18             : /**                         (third part)                              **/
      19             : /**                                                                   **/
      20             : /***********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_pol
      25             : 
      26             : /************************************************************************
      27             :  **                                                                    **
      28             :  **                      Ring membership                               **
      29             :  **                                                                    **
      30             :  ************************************************************************/
      31             : struct charact {
      32             :   GEN q;
      33             :   int isprime;
      34             : };
      35             : static void
      36        1225 : char_update_prime(struct charact *S, GEN p)
      37             : {
      38        1225 :   if (!S->isprime) { S->isprime = 1; S->q = p; }
      39        1225 :   if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
      40        1218 : }
      41             : static void
      42        6594 : char_update_int(struct charact *S, GEN n)
      43             : {
      44        6594 :   if (S->isprime)
      45             :   {
      46           7 :     if (dvdii(n, S->q)) return;
      47           7 :     pari_err_MODULUS("characteristic", S->q, n);
      48             :   }
      49        6587 :   S->q = gcdii(S->q, n);
      50             : }
      51             : static void
      52      176946 : charact(struct charact *S, GEN x)
      53             : {
      54      176946 :   const long tx = typ(x);
      55             :   long i, l;
      56      176946 :   switch(tx)
      57             :   {
      58        5145 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      59        1134 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      60       25711 :     case t_COMPLEX: case t_QUAD:
      61             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      62             :     case t_VEC: case t_COL: case t_MAT:
      63       25711 :       l = lg(x);
      64      175021 :       for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
      65       25697 :       break;
      66           7 :     case t_LIST:
      67           7 :       x = list_data(x);
      68           7 :       if (x) charact(S, x);
      69           7 :       break;
      70             :   }
      71      176918 : }
      72             : static void
      73        4739 : charact_res(struct charact *S, GEN x)
      74             : {
      75        4739 :   const long tx = typ(x);
      76             :   long i, l;
      77        4739 :   switch(tx)
      78             :   {
      79        1449 :     case t_INTMOD:char_update_int(S, gel(x,1)); break;
      80           0 :     case t_FFELT: char_update_prime(S, gel(x,4)); break;
      81          91 :     case t_PADIC: char_update_prime(S, padic_p(x)); break;
      82        1722 :     case t_COMPLEX: case t_QUAD:
      83             :     case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
      84             :     case t_VEC: case t_COL: case t_MAT:
      85        1722 :       l = lg(x);
      86        6132 :       for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
      87        1722 :       break;
      88           0 :     case t_LIST:
      89           0 :       x = list_data(x);
      90           0 :       if (x) charact_res(S, x);
      91           0 :       break;
      92             :   }
      93        4739 : }
      94             : GEN
      95       27622 : characteristic(GEN x)
      96             : {
      97             :   struct charact S;
      98       27622 :   S.q = gen_0; S.isprime = 0;
      99       27622 :   charact(&S, x); return S.q;
     100             : }
     101             : GEN
     102         329 : residual_characteristic(GEN x)
     103             : {
     104             :   struct charact S;
     105         329 :   S.q = gen_0; S.isprime = 0;
     106         329 :   charact_res(&S, x); return S.q;
     107             : }
     108             : 
     109             : int
     110    70276169 : Rg_is_Fp(GEN x, GEN *pp)
     111             : {
     112             :   GEN mod;
     113    70276169 :   switch(typ(x))
     114             :   {
     115     2482676 :   case t_INTMOD:
     116     2482676 :     mod = gel(x,1);
     117     2482676 :     if (!*pp) *pp = mod;
     118     2341976 :     else if (mod != *pp && !equalii(mod, *pp))
     119             :     {
     120           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_Fp");
     121           0 :       return 0;
     122             :     }
     123     2482676 :     return 1;
     124    56415657 :   case t_INT:
     125    56415657 :     return 1;
     126    11377836 :   default: return 0;
     127             :   }
     128             : }
     129             : 
     130             : int
     131    28008859 : RgX_is_FpX(GEN x, GEN *pp)
     132             : {
     133    28008859 :   long i, lx = lg(x);
     134    86881055 :   for (i=2; i<lx; i++)
     135    70250028 :     if (!Rg_is_Fp(gel(x, i), pp))
     136    11377827 :       return 0;
     137    16631027 :   return 1;
     138             : }
     139             : 
     140             : int
     141           0 : RgV_is_FpV(GEN x, GEN *pp)
     142             : {
     143           0 :   long i, lx = lg(x);
     144           0 :   for (i=1; i<lx; i++)
     145           0 :     if (!Rg_is_Fp(gel(x,i), pp)) return 0;
     146           0 :   return 1;
     147             : }
     148             : 
     149             : int
     150           0 : RgM_is_FpM(GEN x, GEN *pp)
     151             : {
     152           0 :   long i, lx = lg(x);
     153           0 :   for (i=1; i<lx; i++)
     154           0 :     if (!RgV_is_FpV(gel(x, i), pp)) return 0;
     155           0 :   return 1;
     156             : }
     157             : 
     158             : int
     159       60802 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
     160             : {
     161             :   GEN pol, mod, p;
     162       60802 :   switch(typ(x))
     163             :   {
     164       26131 :   case t_INTMOD:
     165       26131 :     return Rg_is_Fp(x, pp);
     166        8561 :   case t_INT:
     167        8561 :     return 1;
     168          21 :   case t_POL:
     169          21 :     return RgX_is_FpX(x, pp);
     170       21350 :   case t_FFELT:
     171       21350 :     mod = x; p = FF_p_i(x);
     172       21350 :     if (!*pp) *pp = p;
     173       21350 :     if (!*pT) *pT = mod;
     174       19824 :     else if (typ(*pT)!=t_FFELT || !FF_samefield(*pT,mod))
     175             :     {
     176          42 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     177          42 :       return 0;
     178             :     }
     179       21308 :     return 1;
     180        4585 :   case t_POLMOD:
     181        4585 :     mod = gel(x,1); pol = gel(x, 2);
     182        4585 :     if (!RgX_is_FpX(mod, pp)) return 0;
     183        4585 :     if (typ(pol)==t_POL)
     184             :     {
     185        4578 :       if (!RgX_is_FpX(pol, pp)) return 0;
     186             :     }
     187           7 :     else if (!Rg_is_Fp(pol, pp)) return 0;
     188        4585 :     if (!*pT) *pT = mod;
     189        4431 :     else if (mod != *pT && !gequal(mod, *pT))
     190             :     {
     191           0 :       if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
     192           0 :       return 0;
     193             :     }
     194        4585 :     return 1;
     195             : 
     196         154 :   default: return 0;
     197             :   }
     198             : }
     199             : 
     200             : int
     201     2021494 : RgX_is_ZXX(GEN x, long *v)
     202             : {
     203             :   long i;
     204     8485807 :   for (i = lg(x)-1; i > 1; i--)
     205             :   {
     206     6464522 :     GEN xi = gel(x,i);
     207     6464522 :     long t = typ(xi), vi;
     208     6464522 :     if (t==t_INT) continue;
     209     1013317 :     if (t!=t_POL || !RgX_is_ZX(xi)) return 0;
     210     1013143 :     vi = varn(xi);
     211     1013143 :     if (*v < 0) *v = vi;
     212        2401 :     else if (vi!=*v) return 0;
     213             :   }
     214     2021285 :   return 1;
     215             : }
     216             : 
     217             : int
     218        3381 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
     219             : {
     220        3381 :   long i, lx = lg(x);
     221       63427 :   for (i = 2; i < lx; i++)
     222       60144 :     if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
     223        3283 :   return 1;
     224             : }
     225             : 
     226             : /************************************************************************
     227             :  **                                                                    **
     228             :  **                      Ring conversion                               **
     229             :  **                                                                    **
     230             :  ************************************************************************/
     231             : 
     232             : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
     233             :  * If x is an INTMOD, assume modulus is a multiple of p. */
     234             : GEN
     235    52334410 : Rg_to_Fp(GEN x, GEN p)
     236             : {
     237    52334410 :   if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
     238     4555516 :   switch(typ(x))
     239             :   {
     240      289107 :     case t_INT: return modii(x, p);
     241       18790 :     case t_FRAC: {
     242       18790 :       pari_sp av = avma;
     243       18790 :       GEN z = modii(gel(x,1), p);
     244       18790 :       if (z == gen_0) return gen_0;
     245       18785 :       return gc_INT(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
     246             :     }
     247          70 :     case t_PADIC: return padic_to_Fp(x, p);
     248     4247562 :     case t_INTMOD: {
     249     4247562 :       GEN q = gel(x,1), a = gel(x,2);
     250     4247562 :       if (equalii(q, p)) return icopy(a);
     251          14 :       if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
     252           0 :       return remii(a, p);
     253             :     }
     254           0 :     default: pari_err_TYPE("Rg_to_Fp",x);
     255             :       return NULL; /* LCOV_EXCL_LINE */
     256             :   }
     257             : }
     258             : /* If x is a POLMOD, assume modulus is a multiple of T. */
     259             : GEN
     260     1291986 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
     261             : {
     262     1291986 :   long ta, tx = typ(x), v = get_FpX_var(T);
     263             :   GEN a, b;
     264     1291986 :   if (is_const_t(tx))
     265             :   {
     266       59182 :     if (tx == t_FFELT)
     267             :     {
     268       17355 :       GEN z = FF_to_FpXQ(x);
     269       17355 :       setvarn(z, v);
     270       17355 :       return z;
     271             :     }
     272       41827 :     return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
     273             :   }
     274     1232804 :   switch(tx)
     275             :   {
     276     1230697 :     case t_POLMOD:
     277     1230697 :       b = gel(x,1);
     278     1230697 :       a = gel(x,2); ta = typ(a);
     279     1230697 :       if (is_const_t(ta))
     280        3885 :         return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
     281     1226812 :       b = RgX_to_FpX(b, p); if (varn(b) != v) break;
     282     1226812 :       a = RgX_to_FpX(a, p);
     283     1226812 :       if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
     284     1226812 :         return FpX_rem(a, T, p);
     285           0 :       break;
     286        2107 :     case t_POL:
     287        2107 :       if (varn(x) != v) break;
     288        2100 :       return FpX_rem(RgX_to_FpX(x,p), T, p);
     289           0 :     case t_RFRAC:
     290           0 :       a = Rg_to_FpXQ(gel(x,1), T,p);
     291           0 :       b = Rg_to_FpXQ(gel(x,2), T,p);
     292           0 :       return FpXQ_div(a,b, T,p);
     293             :   }
     294           7 :   pari_err_TYPE("Rg_to_FpXQ",x);
     295             :   return NULL; /* LCOV_EXCL_LINE */
     296             : }
     297             : GEN
     298     3335720 : RgX_to_FpX(GEN x, GEN p)
     299             : {
     300             :   long i, l;
     301     3335720 :   GEN z = cgetg_copy(x, &l); z[1] = x[1];
     302    14765324 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     303     3335720 :   return FpX_renormalize(z, l);
     304             : }
     305             : 
     306             : GEN
     307         140 : RgV_to_FpV(GEN x, GEN p)
     308         483 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
     309             : 
     310             : GEN
     311     1751090 : RgC_to_FpC(GEN x, GEN p)
     312    28499485 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
     313             : 
     314             : GEN
     315      222349 : RgM_to_FpM(GEN x, GEN p)
     316     1973397 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
     317             : 
     318             : GEN
     319      369001 : RgV_to_Flv(GEN x, ulong p)
     320     1676894 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
     321             : 
     322             : GEN
     323      118296 : RgM_to_Flm(GEN x, ulong p)
     324      422998 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
     325             : 
     326             : GEN
     327        5105 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
     328             : {
     329        5105 :   long i, l = lg(x);
     330        5105 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     331       43394 :   for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
     332        5105 :   return FpXQX_renormalize(z, l);
     333             : }
     334             : GEN
     335       49186 : RgX_to_FqX(GEN x, GEN T, GEN p)
     336             : {
     337       49186 :   long i, l = lg(x);
     338       49186 :   GEN z = cgetg(l, t_POL); z[1] = x[1];
     339       49186 :   if (T)
     340       10990 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
     341             :   else
     342      791394 :     for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     343       49186 :   return FpXQX_renormalize(z, l);
     344             : }
     345             : 
     346             : GEN
     347      219128 : RgC_to_FqC(GEN x, GEN T, GEN p)
     348             : {
     349      219128 :   long i, l = lg(x);
     350      219128 :   GEN z = cgetg(l, t_COL);
     351      219128 :   if (T)
     352     1430310 :     for (i = 1; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
     353             :   else
     354           0 :     for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
     355      219128 :   return z;
     356             : }
     357             : 
     358             : GEN
     359       52325 : RgM_to_FqM(GEN x, GEN T, GEN p)
     360      271418 : { pari_APPLY_same(RgC_to_FqC(gel(x, i), T, p)) }
     361             : 
     362             : /* lg(V) > 1 */
     363             : GEN
     364      851487 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
     365             : {
     366      851487 :   pari_sp av = avma;
     367      851487 :   long i, l = lg(V);
     368      851487 :   GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
     369     4201029 :   for(i=2; i<l; i++)
     370             :   {
     371     3349542 :     z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
     372     3349542 :     if ((i & 7) == 0) z = gc_upto(av, z);
     373             :   }
     374      851487 :   return gc_upto(av, FpX_red(z,p));
     375             : }
     376             : 
     377             : GEN
     378       55832 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
     379             : {
     380       55832 :   long i, lz = lg(y);
     381             :   GEN z;
     382       55832 :   if (!T) return FpX_Fp_add(y, x, p);
     383        8666 :   if (lz == 2) return scalarpol(x, varn(y));
     384        7952 :   z = cgetg(lz,t_POL); z[1] = y[1];
     385        7952 :   gel(z,2) = Fq_add(gel(y,2),x, T, p);
     386        7952 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     387             :   else
     388       33145 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     389        7952 :   return z;
     390             : }
     391             : 
     392             : GEN
     393        1059 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
     394             : {
     395        1059 :   long i, lz = lg(y);
     396             :   GEN z;
     397        1059 :   if (!T) return FpX_Fp_sub(y, x, p);
     398        1059 :   if (lz == 2) return scalarpol(x, varn(y));
     399        1059 :   z = cgetg(lz,t_POL); z[1] = y[1];
     400        1059 :   gel(z,2) = Fq_sub(gel(y,2), x, T, p);
     401        1059 :   if (lz == 3) z = FpXX_renormalize(z,lz);
     402             :   else
     403       10278 :     for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
     404        1059 :   return z;
     405             : }
     406             : 
     407             : GEN
     408      149052 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
     409             : {
     410             :   long i, lP;
     411      149052 :   GEN res = cgetg_copy(P, &lP); res[1] = P[1];
     412      918799 :   for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
     413      149052 :   gel(res,lP-1) = gen_1; return res;
     414             : }
     415             : 
     416             : GEN
     417       38189 : FpXQX_normalize(GEN z, GEN T, GEN p)
     418             : {
     419             :   GEN lc;
     420       38189 :   if (lg(z) == 2) return z;
     421       38175 :   lc = leading_coeff(z);
     422       38175 :   if (typ(lc) == t_POL)
     423             :   {
     424        2167 :     if (lg(lc) > 3) /* nonconstant */
     425        1902 :       return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
     426             :     /* constant */
     427         265 :     lc = gel(lc,2);
     428         265 :     z = shallowcopy(z);
     429         265 :     gel(z, lg(z)-1) = lc;
     430             :   }
     431             :   /* lc a t_INT */
     432       36273 :   if (equali1(lc)) return z;
     433          87 :   return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
     434             : }
     435             : 
     436             : GEN
     437      390934 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
     438             : {
     439             :   pari_sp av;
     440             :   GEN p1, r;
     441      390934 :   long j, i=lg(x)-1;
     442      390934 :   if (i<=2)
     443       45107 :     return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
     444      345827 :   av=avma; p1=gel(x,i);
     445             :   /* specific attention to sparse polynomials (see poleval)*/
     446             :   /*You've guessed it! It's a copy-paste(tm)*/
     447     1150438 :   for (i--; i>=2; i=j-1)
     448             :   {
     449      805296 :     for (j=i; !signe(gel(x,j)); j--)
     450         686 :       if (j==2)
     451             :       {
     452         490 :         if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
     453         490 :         return gc_upto(av, Fq_mul(p1,y, T, p));
     454             :       }
     455      804610 :     r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
     456      804610 :     p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
     457             :   }
     458      345338 :   return gc_upto(av, p1);
     459             : }
     460             : 
     461             : GEN
     462       97321 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
     463             : {
     464       97321 :   long i, lb = lg(Q);
     465             :   GEN z;
     466       97321 :   if (!T) return FpXY_evalx(Q, x, p);
     467       86961 :   z = cgetg(lb, t_POL); z[1] = Q[1];
     468      462945 :   for (i=2; i<lb; i++)
     469             :   {
     470      375984 :     GEN q = gel(Q,i);
     471      375984 :     gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
     472             :   }
     473       86961 :   return FpXQX_renormalize(z, lb);
     474             : }
     475             : 
     476             : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
     477             : GEN
     478       14623 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
     479             : {
     480       14623 :   pari_sp av = avma;
     481       14623 :   if (!T) return FpXY_eval(Q, y, x, p);
     482         966 :   return gc_upto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
     483             : }
     484             : 
     485             : /* a X^d */
     486             : GEN
     487    13640946 : monomial(GEN a, long d, long v)
     488             : {
     489             :   long i, n;
     490             :   GEN P;
     491    13640946 :   if (d < 0) {
     492          14 :     if (isrationalzero(a)) return pol_0(v);
     493          14 :     retmkrfrac(a, pol_xn(-d, v));
     494             :   }
     495    13640932 :   if (gequal0(a))
     496             :   {
     497       93989 :     if (isexactzero(a)) return scalarpol_shallow(a,v);
     498           0 :     n = d+2; P = cgetg(n+1, t_POL);
     499           0 :     P[1] = evalsigne(0) | evalvarn(v);
     500             :   }
     501             :   else
     502             :   {
     503    13546941 :     n = d+2; P = cgetg(n+1, t_POL);
     504    13546943 :     P[1] = evalsigne(1) | evalvarn(v);
     505             :   }
     506    32924128 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     507    13546943 :   gel(P,i) = a; return P;
     508             : }
     509             : GEN
     510      157969 : monomialcopy(GEN a, long d, long v)
     511             : {
     512             :   long i, n;
     513             :   GEN P;
     514      157969 :   if (d < 0) {
     515          14 :     if (isrationalzero(a)) return pol_0(v);
     516          14 :     retmkrfrac(gcopy(a), pol_xn(-d, v));
     517             :   }
     518      157955 :   if (gequal0(a))
     519             :   {
     520          14 :     if (isexactzero(a)) return scalarpol(a,v);
     521           7 :     n = d+2; P = cgetg(n+1, t_POL);
     522           7 :     P[1] = evalsigne(0) | evalvarn(v);
     523             :   }
     524             :   else
     525             :   {
     526      157941 :     n = d+2; P = cgetg(n+1, t_POL);
     527      157941 :     P[1] = evalsigne(1) | evalvarn(v);
     528             :   }
     529      314678 :   for (i = 2; i < n; i++) gel(P,i) = gen_0;
     530      157948 :   gel(P,i) = gcopy(a); return P;
     531             : }
     532             : GEN
     533      325931 : pol_x_powers(long N, long v)
     534             : {
     535      325931 :   GEN L = cgetg(N+1,t_VEC);
     536             :   long i;
     537      789016 :   for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
     538      325939 :   return L;
     539             : }
     540             : 
     541             : GEN
     542           0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
     543             : {
     544           0 :   return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
     545             : }
     546             : 
     547             : GEN
     548           0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
     549             : {
     550           0 :   return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
     551             : }
     552             : 
     553             : /*******************************************************************/
     554             : /*                                                                 */
     555             : /*                             Fq                                  */
     556             : /*                                                                 */
     557             : /*******************************************************************/
     558             : 
     559             : GEN
     560    11592895 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
     561             : {
     562             :   (void)T;
     563    11592895 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     564             :   {
     565     1143687 :     case 0: return Fp_add(x,y,p);
     566      748136 :     case 1: return FpX_Fp_add(x,y,p);
     567       91991 :     case 2: return FpX_Fp_add(y,x,p);
     568     9609081 :     case 3: return FpX_add(x,y,p);
     569             :   }
     570             :   return NULL;/*LCOV_EXCL_LINE*/
     571             : }
     572             : 
     573             : GEN
     574     8563061 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
     575             : {
     576             :   (void)T;
     577     8563061 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     578             :   {
     579      256368 :     case 0: return Fp_sub(x,y,p);
     580       24540 :     case 1: return FpX_Fp_sub(x,y,p);
     581       23896 :     case 2: return Fp_FpX_sub(x,y,p);
     582     8258257 :     case 3: return FpX_sub(x,y,p);
     583             :   }
     584             :   return NULL;/*LCOV_EXCL_LINE*/
     585             : }
     586             : 
     587             : GEN
     588     1080478 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
     589             : {
     590             :   (void)T;
     591     1080478 :   return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
     592             : }
     593             : 
     594             : GEN
     595       81354 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
     596             : {
     597             :   (void)T;
     598       81354 :   return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
     599             : }
     600             : 
     601             : /* If T==NULL do not reduce*/
     602             : GEN
     603     8608662 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
     604             : {
     605     8608662 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     606             :   {
     607     1037917 :     case 0: return Fp_mul(x,y,p);
     608      128565 :     case 1: return FpX_Fp_mul(x,y,p);
     609      394685 :     case 2: return FpX_Fp_mul(y,x,p);
     610     7047495 :     case 3: if (T) return FpXQ_mul(x,y,T,p);
     611     4478770 :             else return FpX_mul(x,y,p);
     612             :   }
     613             :   return NULL;/*LCOV_EXCL_LINE*/
     614             : }
     615             : 
     616             : /* If T==NULL do not reduce*/
     617             : GEN
     618      490543 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
     619             : {
     620             :   (void) T;
     621      490543 :   return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
     622             : }
     623             : 
     624             : /* y t_INT */
     625             : GEN
     626       44041 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
     627             : {
     628             :   (void)T;
     629        6822 :   return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
     630       50863 :                           : Fp_mul(x,y,p);
     631             : }
     632             : /* If T==NULL do not reduce*/
     633             : GEN
     634      611418 : Fq_sqr(GEN x, GEN T, GEN p)
     635             : {
     636      611418 :   if (typ(x) == t_POL)
     637             :   {
     638       70585 :     if (T) return FpXQ_sqr(x,T,p);
     639           0 :     else return FpX_sqr(x,p);
     640             :   }
     641             :   else
     642      540833 :     return Fp_sqr(x,p);
     643             : }
     644             : 
     645             : GEN
     646           0 : Fq_neg_inv(GEN x, GEN T, GEN p)
     647             : {
     648           0 :   if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
     649           0 :   return FpXQ_inv(FpX_neg(x,p),T,p);
     650             : }
     651             : 
     652             : GEN
     653           0 : Fq_invsafe(GEN x, GEN pol, GEN p)
     654             : {
     655           0 :   if (typ(x) == t_INT) return Fp_invsafe(x,p);
     656           0 :   return FpXQ_invsafe(x,pol,p);
     657             : }
     658             : 
     659             : GEN
     660       89312 : Fq_inv(GEN x, GEN pol, GEN p)
     661             : {
     662       89312 :   if (typ(x) == t_INT) return Fp_inv(x,p);
     663       81546 :   return FpXQ_inv(x,pol,p);
     664             : }
     665             : 
     666             : GEN
     667      343791 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
     668             : {
     669      343791 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     670             :   {
     671      318402 :     case 0: return Fp_div(x,y,p);
     672       16702 :     case 1: return FpX_Fp_div(x,y,p);
     673         406 :     case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
     674        8281 :     case 3: return FpXQ_div(x,y,pol,p);
     675             :   }
     676             :   return NULL;/*LCOV_EXCL_LINE*/
     677             : }
     678             : 
     679             : GEN
     680     1088154 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
     681             : {
     682     1088154 :   if (typ(x) == t_INT) return Fp_pow(x,n,p);
     683      137458 :   return FpXQ_pow(x,n,pol,p);
     684             : }
     685             : 
     686             : GEN
     687       15050 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
     688             : {
     689       15050 :   if (typ(x) == t_INT) return Fp_powu(x,n,p);
     690        1267 :   return FpXQ_powu(x,n,pol,p);
     691             : }
     692             : 
     693             : GEN
     694     1892926 : Fq_sqrt(GEN x, GEN T, GEN p)
     695             : {
     696     1892926 :   if (typ(x) == t_INT)
     697             :   {
     698     1825064 :     if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
     699        9621 :     x = scalarpol_shallow(x, get_FpX_var(T));
     700             :   }
     701       77483 :   return FpXQ_sqrt(x,T,p);
     702             : }
     703             : GEN
     704      170786 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
     705             : {
     706      170786 :   if (typ(x) == t_INT)
     707             :   {
     708             :     long d;
     709      170415 :     if (!T) return Fp_sqrtn(x,n,p,zeta);
     710         126 :     d = get_FpX_degree(T);
     711         126 :     if (ugcdiu(n,d) == 1)
     712             :     {
     713          56 :       if (!zeta) return Fp_sqrtn(x,n,p,NULL);
     714             :       /* gcd(n,p-1)=gcd(n,q-1): same number of solutions in Fp and F_q */
     715          21 :       if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
     716          14 :         return Fp_sqrtn(x,n,p,zeta);
     717             :     }
     718          77 :     x = scalarpol(x, get_FpX_var(T)); /* left on stack */
     719             :   }
     720         448 :   return FpXQ_sqrtn(x,n,T,p,zeta);
     721             : }
     722             : 
     723             : struct _Fq_field
     724             : {
     725             :   GEN T, p;
     726             : };
     727             : 
     728             : static GEN
     729      303247 : _Fq_red(void *E, GEN x)
     730      303247 : { struct _Fq_field *s = (struct _Fq_field *)E;
     731      303247 :   return Fq_red(x, s->T, s->p);
     732             : }
     733             : 
     734             : static GEN
     735      362523 : _Fq_add(void *E, GEN x, GEN y)
     736             : {
     737             :   (void) E;
     738      362523 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     739             :   {
     740          14 :     case 0: return addii(x,y);
     741           0 :     case 1: return ZX_Z_add(x,y);
     742       15918 :     case 2: return ZX_Z_add(y,x);
     743      346591 :     default: return ZX_add(x,y);
     744             :   }
     745             : }
     746             : 
     747             : static GEN
     748      315028 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
     749             : 
     750             : static GEN
     751      243341 : _Fq_mul(void *E, GEN x, GEN y)
     752             : {
     753             :   (void) E;
     754      243341 :   switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
     755             :   {
     756         679 :     case 0: return mulii(x,y);
     757        1085 :     case 1: return ZX_Z_mul(x,y);
     758        1043 :     case 2: return ZX_Z_mul(y,x);
     759      240534 :     default: return ZX_mul(x,y);
     760             :   }
     761             : }
     762             : 
     763             : static GEN
     764        9331 : _Fq_inv(void *E, GEN x)
     765        9331 : { struct _Fq_field *s = (struct _Fq_field *)E;
     766        9331 :   return Fq_inv(x,s->T,s->p);
     767             : }
     768             : 
     769             : static int
     770       38388 : _Fq_equal0(GEN x) { return signe(x)==0; }
     771             : 
     772             : static GEN
     773        4151 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
     774             : 
     775             : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
     776             :                                        _Fq_inv,_Fq_equal0,_Fq_s};
     777             : 
     778        4725 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
     779             : {
     780        4725 :   if (!T)
     781           0 :     return get_Fp_field(E, p);
     782             :   else
     783             :   {
     784        4725 :     GEN z = new_chunk(sizeof(struct _Fq_field));
     785        4725 :     struct _Fq_field *e = (struct _Fq_field *) z;
     786        4725 :     e->T = T; e->p  = p; *E = (void*)e;
     787        4725 :     return &Fq_field;
     788             :   }
     789             : }
     790             : 
     791             : /*******************************************************************/
     792             : /*                                                                 */
     793             : /*                             Fq[X]                               */
     794             : /*                                                                 */
     795             : /*******************************************************************/
     796             : /* FpV_red(vecbinomial(n), p) */
     797             : static GEN
     798           0 : vecbinomial_Fp(long n, GEN p)
     799             : {
     800           0 :   GEN C = vecbinomial(n) + 1;
     801           0 :   long k, d = (n + 1) >> 1;
     802           0 :   for (k = d; k >= 1; k--)
     803             :   {
     804           0 :     GEN a = gel(C,k);
     805           0 :     if (cmpii(a, p) < 0) break;
     806           0 :     gel(C,k) = gel(C, n-k) = remii(a, p);
     807             :   }
     808           0 :   return C - 1;
     809             : }
     810             : /* (X+u)^n */
     811             : static GEN
     812         434 : Fp_XpN_powu(GEN u, ulong n, GEN p, long v)
     813             : {
     814             :   pari_sp av;
     815             :   ulong k;
     816             :   GEN B, C, V;
     817         434 :   if (!n) return pol_1(v);
     818         434 :   if (is_pm1(u))
     819         434 :     return Xpm1_powu(n, signe(u), v);
     820           0 :   av = avma;
     821           0 :   V = Fp_powers(u, n, p);
     822           0 :   B = vecbinomial_Fp(n, p);
     823           0 :   C = cgetg(n+3, t_POL);
     824           0 :   C[1] = evalsigne(1)| evalvarn(v);
     825           0 :   for (k=1; k <= n+1; k++)
     826           0 :     gel(C,k+1) = Fp_mul(gel(V,n+2-k), gel(B,k), p);
     827           0 :   return gc_upto(av, C);
     828             : }
     829             : 
     830             : static GEN
     831         700 : FpX_Fp_translate_basecase(GEN P, GEN c, GEN p)
     832             : {
     833         700 :   pari_sp av = avma;
     834             :   GEN Q;
     835             :   long i, k, n;
     836             : 
     837         700 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     838         560 :   Q = leafcopy(P); n = degpol(P);
     839        1316 :   for (i=1; i<=n; i++)
     840             :   {
     841        2016 :     for (k=n-i; k<n; k++)
     842        1260 :       gel(Q,2+k) = Fp_add(gel(Q,2+k), Fp_mul(c, gel(Q,2+k+1), p), p);
     843             : 
     844         756 :     if (gc_needed(av,2))
     845             :     {
     846           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FpX_Fp_translate, i = %ld/%ld", i,n);
     847           0 :       Q = gc_GEN(av, Q);
     848             :     }
     849             :   }
     850         560 :   return gc_GEN(av, FpX_renormalize(Q, lg(Q)));
     851             : }
     852             : 
     853             : GEN
     854        1134 : FpX_Fp_translate(GEN P, GEN c, GEN p)
     855             : {
     856        1134 :   pari_sp av = avma;
     857        1134 :   long n = degpol(P);
     858        1134 :   if (n<=3 || 25*(n-3) < expi(p))
     859         700 :     return FpX_Fp_translate_basecase(P, c, p);
     860             :   else
     861             :   {
     862         434 :     long d = n >> 1;
     863         434 :     GEN Q = FpX_Fp_translate(RgX_shift_shallow(P, -d), c, p);
     864         434 :     GEN R = FpX_Fp_translate(RgXn_red_shallow(P, d), c, p);
     865         434 :     GEN S = Fp_XpN_powu(c, d, p, varn(P));
     866         434 :     return gc_upto(av, FpX_add(FpX_mul(Q, S, p), R, p));
     867             :   }
     868             : }
     869             : 
     870             : /* (X+u)^n mod (T,p) */
     871             : static GEN
     872           0 : FpXQX_XpN_powu(GEN u, ulong n, GEN T, GEN p, long v)
     873             : {
     874             :   pari_sp av;
     875             :   ulong k;
     876             :   GEN B, C, V;
     877           0 :   if (!n) return pol_1(v);
     878           0 :   av = avma;
     879           0 :   V = FpXQ_powers(u, n, T, p);
     880           0 :   B = vecbinomial_Fp(n, p);
     881           0 :   C = cgetg(n+3, t_POL);
     882           0 :   C[1] = evalsigne(1)| evalvarn(v);
     883           0 :   for (k=1; k <= n+1; k++)
     884           0 :     gel(C,k+1) = Fq_mul(gel(V,n+2-k), gel(B,k), T, p);
     885           0 :   return gc_upto(av, C);
     886             : }
     887             : 
     888             : static GEN
     889       33887 : FpXQX_FpXQ_translate_basecase(GEN P, GEN c, GEN T, GEN p)
     890             : {
     891       33887 :   pari_sp av = avma;
     892             :   GEN Q;
     893             :   long i, k, n;
     894             : 
     895             :   /* signe works for t_(INT|POL) */
     896       33887 :   if (!signe(P) || !signe(c)) return RgX_copy(P);
     897       33887 :   Q = leafcopy(P); n = degpol(P);
     898      150080 :   for (i=1; i<=n; i++)
     899             :   {
     900      433594 :     for (k=n-i; k<n; k++)
     901      317401 :       gel(Q,2+k) = Fq_add(gel(Q,2+k), Fq_mul(c, gel(Q,2+k+1), T, p), T, p);
     902             : 
     903      116193 :     if (gc_needed(av,2))
     904             :     {
     905           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FqX_Fq_translate, i = %ld/%ld", i,n);
     906           0 :       Q = gc_GEN(av, Q);
     907             :     }
     908             :   }
     909       33887 :   return gc_GEN(av, FpXQX_renormalize(Q, lg(Q)));
     910             : }
     911             : 
     912             : GEN
     913       33887 : FpXQX_FpXQ_translate(GEN P, GEN c, GEN T, GEN p)
     914             : {
     915       33887 :   pari_sp av = avma;
     916       33887 :   long n = degpol(P);
     917       33887 :   if (n < 220)
     918       33887 :     return FpXQX_FpXQ_translate_basecase(P, c, T, p);
     919             :   else
     920             :   {
     921           0 :     long d = n >> 1;
     922           0 :     GEN Q = FpXQX_FpXQ_translate(RgX_shift_shallow(P, -d), c, T, p);
     923           0 :     GEN R = FpXQX_FpXQ_translate(RgXn_red_shallow(P, d), c, T, p);
     924           0 :     GEN S = FpXQX_XpN_powu(c, d, T, p, varn(P));
     925           0 :     return gc_upto(av, FpXX_add(FpXQX_mul(Q, S, T, p), R, p));
     926             :   }
     927             : }
     928             : 
     929             : GEN
     930       33880 : FqX_Fq_translate(GEN x, GEN y, GEN T, GEN p)
     931             : {
     932       33880 :   if (!T) return FpX_Fp_translate(x,y,p);
     933       33880 :   if (typ(y)==t_INT)
     934             :   {
     935           0 :     pari_sp av = avma;
     936           0 :     y = scalarpol_shallow(y,varn(T));
     937           0 :     return gc_upto(av, FpXQX_FpXQ_translate(x,y,T,p));
     938             :   }
     939       33880 :   return FpXQX_FpXQ_translate(x,y,T,p);
     940             : }
     941             : 
     942             : GEN
     943       40453 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
     944             : {
     945       40453 :   pari_sp ltop = avma;
     946             :   long k;
     947             :   GEN W;
     948       40453 :   if (lgefint(p) == 3)
     949             :   {
     950       31742 :     ulong pp = p[2];
     951       31742 :     GEN Tl = ZX_to_Flx(T, pp);
     952       31746 :     GEN Vl = FqC_to_FlxqC(V, Tl, pp);
     953       31745 :     Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
     954       31746 :     return gc_upto(ltop, FlxX_to_ZXX(Tl));
     955             :   }
     956        8711 :   W = cgetg(lg(V),t_VEC);
     957       78117 :   for(k=1; k < lg(V); k++)
     958       69406 :     gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
     959        8711 :   return gc_upto(ltop, FpXQXV_prod(W, T, p));
     960             : }
     961             : 
     962             : GEN
     963      109458 : FqV_red(GEN x, GEN T, GEN p)
     964      778133 : { pari_APPLY_type(t_VEC, Fq_red(gel(x,i), T, p)) }
     965             : 
     966             : GEN
     967       24058 : FqC_red(GEN x, GEN T, GEN p)
     968      163384 : { pari_APPLY_type(t_COL, Fq_red(gel(x,i), T, p)) }
     969             : 
     970             : GEN
     971        1701 : FqM_red(GEN x, GEN T, GEN p)
     972        5411 : { pari_APPLY_same(FqC_red(gel(x,i), T, p)) }
     973             : 
     974             : GEN
     975           0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
     976             : {
     977           0 :   if (!T) return FpC_add(x, y, p);
     978           0 :   pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
     979             : }
     980             : 
     981             : GEN
     982           0 : FqM_add(GEN x, GEN y, GEN T, GEN p)
     983             : {
     984           0 :   if (!T) return FpM_add(x, y, p);
     985           0 :   pari_APPLY_same(FqC_add(gel(x,i), gel(y,i), T, p))
     986             : }
     987             : 
     988             : GEN
     989           0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
     990             : {
     991           0 :   if (!T) return FpC_sub(x, y, p);
     992           0 :   pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
     993             : }
     994             : 
     995             : GEN
     996           0 : FqM_sub(GEN x, GEN y, GEN T, GEN p)
     997             : {
     998           0 :   if (!T) return FpM_sub(x, y, p);
     999           0 :   pari_APPLY_same(FqC_sub(gel(x,i), gel(y,i), T, p))
    1000             : }
    1001             : 
    1002             : GEN
    1003           0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
    1004             : {
    1005           0 :   if (!T) return FpC_Fp_mul(x, y, p);
    1006           0 :   pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
    1007             : }
    1008             : 
    1009             : GEN
    1010         105 : FqC_FqV_mul(GEN x, GEN y, GEN T, GEN p)
    1011             : {
    1012         105 :   long i,j, lx=lg(x), ly=lg(y);
    1013             :   GEN z;
    1014         105 :   if (ly==1) return cgetg(1,t_MAT);
    1015         105 :   z = cgetg(ly,t_MAT);
    1016         819 :   for (j=1; j < ly; j++)
    1017             :   {
    1018         714 :     GEN zj = cgetg(lx,t_COL);
    1019        4200 :     for (i=1; i<lx; i++) gel(zj,i) = Fq_mul(gel(x,i),gel(y,j), T, p);
    1020         714 :     gel(z, j) = zj;
    1021             :   }
    1022         105 :   return z;
    1023             : }
    1024             : 
    1025             : GEN
    1026        5467 : FpXC_center(GEN x, GEN p, GEN pov2)
    1027       41476 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
    1028             : 
    1029             : GEN
    1030      109023 : FqC_to_FlxqC(GEN x, GEN T, ulong p)
    1031      109023 : { long sv = get_Flx_var(T);
    1032     4834761 :   pari_APPLY_type(t_COL,typ(gel(x,i))==t_INT ?
    1033             :                   Z_to_Flx(gel(x,i), p, sv): ZX_to_Flx(gel(x,i), p)) }
    1034             : 
    1035             : GEN
    1036        8708 : FqM_to_FlxqM(GEN x, GEN T, ulong p)
    1037       85985 : { pari_APPLY_same(FqC_to_FlxqC(gel(x,i), T, p)) }
    1038             : 
    1039             : GEN
    1040        1800 : FpXM_center(GEN x, GEN p, GEN pov2)
    1041        7267 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
    1042             : 
    1043             : /*******************************************************************/
    1044             : /*                                                                 */
    1045             : /*                          GENERIC CRT                            */
    1046             : /*                                                                 */
    1047             : /*******************************************************************/
    1048             : static GEN
    1049     9414376 : primelist(forprime_t *S, long n, GEN dB)
    1050             : {
    1051     9414376 :   GEN P = cgetg(n+1, t_VECSMALL);
    1052     9414352 :   long i = 1;
    1053             :   ulong p;
    1054    22385065 :   while (i <= n && (p = u_forprime_next(S)))
    1055    12970713 :     if (!dB || umodiu(dB, p)) P[i++] = p;
    1056     9414350 :   return P;
    1057             : }
    1058             : 
    1059             : void
    1060     8826280 : gen_inccrt_i(const char *str, GEN worker, GEN dB, long n, long mmin,
    1061             :              forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
    1062             :              GEN center(GEN, GEN, GEN))
    1063             : {
    1064     8826280 :   long m = mmin? minss(mmin, n): usqrt(n);
    1065             :   GEN  H, P, mod;
    1066             :   pari_timer ti;
    1067     8826285 :   if (DEBUGLEVEL > 4)
    1068             :   {
    1069           0 :     timer_start(&ti);
    1070           0 :     err_printf("%s: nb primes: %ld\n",str, n);
    1071             :   }
    1072     8826272 :   if (m == 1)
    1073             :   {
    1074     8510282 :     GEN P = primelist(S, n, dB);
    1075     8510241 :     GEN done = closure_callgen1(worker, P);
    1076     8510229 :     H = gel(done,1);
    1077     8510229 :     mod = gel(done,2);
    1078     8510229 :     if (!*pH && center) H = center(H, mod, shifti(mod,-1));
    1079     8510138 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
    1080             :   }
    1081             :   else
    1082             :   {
    1083      315990 :     long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
    1084             :     struct pari_mt pt;
    1085      315990 :     long pending = 0;
    1086      315990 :     H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
    1087      315990 :     mt_queue_start_lim(&pt, worker, m);
    1088     1286080 :     for (i=1; i<=m || pending; i++)
    1089             :     {
    1090             :       GEN done;
    1091      970090 :       GEN pr = i <= m ? mkvec(primelist(S, i<=r ? s: s-1, dB)): NULL;
    1092      970092 :       mt_queue_submit(&pt, i, pr);
    1093      970090 :       done = mt_queue_get(&pt, NULL, &pending);
    1094      970090 :       if (done)
    1095             :       {
    1096      904098 :         di++;
    1097      904098 :         gel(H, di) = gel(done,1);
    1098      904098 :         gel(P, di) = gel(done,2);
    1099      904098 :         if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
    1100             :       }
    1101             :     }
    1102      315990 :     mt_queue_end(&pt);
    1103      315990 :     if (DEBUGLEVEL>5) err_printf("\n");
    1104      315990 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
    1105      315990 :     H = crt(H, P, &mod);
    1106      315990 :     if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
    1107             :   }
    1108     8826128 :   if (*pH) H = crt(mkvec2(*pH, H), mkvec2(*pmod, mod), &mod);
    1109     8826128 :   *pH = H; *pmod = mod;
    1110     8826128 : }
    1111             : void
    1112     3072995 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
    1113             :            forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
    1114             :            GEN center(GEN, GEN, GEN))
    1115             : {
    1116     3072995 :   pari_sp av = avma;
    1117     3072995 :   gen_inccrt_i(str, worker, dB, n, mmin, S, pH, pmod, crt, center);
    1118     3072897 :   (void)gc_all(av, 2, pH, pmod);
    1119     3073070 : }
    1120             : 
    1121             : GEN
    1122     2285212 : gen_crt(const char *str, GEN worker, forprime_t *S, GEN dB, ulong bound, long mmin, GEN *pmod,
    1123             :         GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
    1124             : {
    1125     2285212 :   GEN mod = gen_1, H = NULL;
    1126             :   ulong e;
    1127             : 
    1128     2285212 :   bound++;
    1129     4570488 :   while (bound > (e = expi(mod)))
    1130             :   {
    1131     2285145 :     long n = (bound - e) / expu(S->p) + 1;
    1132     2285183 :     gen_inccrt(str, worker, dB, n, mmin, S, &H, &mod, crt, center);
    1133             :   }
    1134     2285247 :   if (pmod) *pmod = mod;
    1135     2285247 :   return H;
    1136             : }
    1137             : 
    1138             : /*******************************************************************/
    1139             : /*                                                                 */
    1140             : /*                          MODULAR GCD                            */
    1141             : /*                                                                 */
    1142             : /*******************************************************************/
    1143             : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
    1144             : static GEN
    1145     5162722 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
    1146             : {
    1147     5162722 :   ulong d, amod = umodiu(a, p);
    1148     5162732 :   pari_sp av = avma;
    1149             :   GEN ax;
    1150             : 
    1151     5162732 :   if (b == amod) return NULL;
    1152     2126722 :   d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
    1153     2127116 :   if (d >= 1 + (p>>1))
    1154     1037849 :     ax = subii(a, mului(p-d, q));
    1155             :   else
    1156             :   {
    1157     1089267 :     ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
    1158     1088854 :     if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
    1159             :   }
    1160     2126348 :   return gc_INT(av, ax);
    1161             : }
    1162             : GEN
    1163         406 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
    1164             : GEN
    1165       32487 : ZX_init_CRT(GEN Hp, ulong p, long v)
    1166             : {
    1167       32487 :   long i, l = lg(Hp), lim = (long)(p>>1);
    1168       32487 :   GEN H = cgetg(l, t_POL);
    1169       32487 :   H[1] = evalsigne(1) | evalvarn(v);
    1170      801135 :   for (i=2; i<l; i++)
    1171      768648 :     gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
    1172       32487 :   return ZX_renormalize(H,l);
    1173             : }
    1174             : 
    1175             : GEN
    1176        5978 : ZM_init_CRT(GEN Hp, ulong p)
    1177             : {
    1178        5978 :   long i,j, m, l = lg(Hp), lim = (long)(p>>1);
    1179        5978 :   GEN c, cp, H = cgetg(l, t_MAT);
    1180        5978 :   if (l==1) return H;
    1181        5978 :   m = lgcols(Hp);
    1182       20223 :   for (j=1; j<l; j++)
    1183             :   {
    1184       14245 :     cp = gel(Hp,j);
    1185       14245 :     c = cgetg(m, t_COL);
    1186       14245 :     gel(H,j) = c;
    1187      169316 :     for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
    1188             :   }
    1189        5978 :   return H;
    1190             : }
    1191             : 
    1192             : int
    1193        7742 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
    1194             : {
    1195        7742 :   GEN h, q = *ptq, qp = muliu(q,p);
    1196        7742 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1197        7742 :   int stable = 1;
    1198        7742 :   h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
    1199        7742 :   if (h) { *H = h; stable = 0; }
    1200        7742 :   *ptq = qp; return stable;
    1201             : }
    1202             : 
    1203             : static int
    1204      148373 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
    1205             : {
    1206      148373 :   GEN H = *ptH, h, qp2 = shifti(qp,-1);
    1207      148374 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1208      148378 :   long i, l = lg(H), lp = lg(Hp);
    1209      148378 :   int stable = 1;
    1210             : 
    1211      148378 :   if (l < lp)
    1212             :   { /* degree increases */
    1213           0 :     GEN x = cgetg(lp, t_POL);
    1214           0 :     for (i=1; i<l; i++)  x[i] = H[i];
    1215           0 :     for (   ; i<lp; i++) gel(x,i) = gen_0;
    1216           0 :     *ptH = H = x;
    1217           0 :     stable = 0;
    1218      148378 :   } else if (l > lp)
    1219             :   { /* degree decreases */
    1220           0 :     GEN x = cgetg(l, t_VECSMALL);
    1221           0 :     for (i=1; i<lp; i++)  x[i] = Hp[i];
    1222           0 :     for (   ; i<l; i++) x[i] = 0;
    1223           0 :     Hp = x; lp = l;
    1224             :   }
    1225     4939667 :   for (i=2; i<lp; i++)
    1226             :   {
    1227     4791354 :     h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
    1228     4791289 :     if (h) { gel(H,i) = h; stable = 0; }
    1229             :   }
    1230      148313 :   (void)ZX_renormalize(H,lp);
    1231      148379 :   return stable;
    1232             : }
    1233             : 
    1234             : int
    1235           0 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
    1236             : {
    1237           0 :   GEN q = *ptq, qp = muliu(q,p);
    1238           0 :   int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
    1239           0 :   *ptq = qp; return stable;
    1240             : }
    1241             : 
    1242             : int
    1243        7787 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1244             : {
    1245        7787 :   GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1246        7787 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1247        7787 :   long i,j, l = lg(H), m = lgcols(H);
    1248        7787 :   int stable = 1;
    1249       27451 :   for (j=1; j<l; j++)
    1250      205516 :     for (i=1; i<m; i++)
    1251             :     {
    1252      185852 :       h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
    1253      185852 :       if (h) { gcoeff(H,i,j) = h; stable = 0; }
    1254             :     }
    1255        7787 :   *ptq = qp; return stable;
    1256             : }
    1257             : 
    1258             : GEN
    1259         679 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
    1260             : {
    1261             :   long i, j, k;
    1262             :   GEN H;
    1263         679 :   long m, l = lg(Hp), lim = (long)(p>>1), n;
    1264         679 :   H = cgetg(l, t_MAT);
    1265         679 :   if (l==1) return H;
    1266         679 :   m = lgcols(Hp);
    1267         679 :   n = deg + 3;
    1268        2268 :   for (j=1; j<l; j++)
    1269             :   {
    1270        1589 :     GEN cp = gel(Hp,j);
    1271        1589 :     GEN c = cgetg(m, t_COL);
    1272        1589 :     gel(H,j) = c;
    1273       24465 :     for (i=1; i<m; i++)
    1274             :     {
    1275       22876 :       GEN dp = gel(cp, i);
    1276       22876 :       long l = lg(dp);
    1277       22876 :       GEN d = cgetg(n, t_POL);
    1278       22876 :       gel(c, i) = d;
    1279       22876 :       d[1] = dp[1] | evalsigne(1);
    1280       46459 :       for (k=2; k<l; k++)
    1281       23583 :         gel(d,k) = stoi(Fl_center(dp[k], p, lim));
    1282       45493 :       for (   ; k<n; k++)
    1283       22617 :         gel(d,k) = gen_0;
    1284             :     }
    1285             :   }
    1286         679 :   return H;
    1287             : }
    1288             : 
    1289             : int
    1290         653 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
    1291             : {
    1292         653 :   GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
    1293         653 :   ulong qinv = Fl_inv(umodiu(q,p), p);
    1294         653 :   long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
    1295         653 :   int stable = 1;
    1296        2225 :   for (j=1; j<l; j++)
    1297       90418 :     for (i=1; i<m; i++)
    1298             :     {
    1299       88846 :       GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
    1300       88846 :       long lh = lg(hp);
    1301      246641 :       for (k=2; k<lh; k++)
    1302             :       {
    1303      157795 :         v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
    1304      157795 :         if (v) { gel(h,k) = v; stable = 0; }
    1305             :       }
    1306      108763 :       for (; k<n; k++)
    1307             :       {
    1308       19917 :         v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
    1309       19917 :         if (v) { gel(h,k) = v; stable = 0; }
    1310             :       }
    1311             :     }
    1312         653 :   *ptq = qp; return stable;
    1313             : }
    1314             : 
    1315             : /* record the degrees of Euclidean remainders (make them as large as
    1316             :  * possible : smaller values correspond to a degenerate sequence) */
    1317             : static void
    1318       23740 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
    1319             : {
    1320             :   long da,db,dc, ind;
    1321       23740 :   pari_sp av = avma;
    1322             : 
    1323       23740 :   if (lgpol(a)==0 || lgpol(b)==0) return;
    1324       22473 :   da = degpol(a);
    1325       22473 :   db = degpol(b);
    1326       22473 :   if (db > da)
    1327           0 :   { swapspec(a,b, da,db); }
    1328       22473 :   else if (!da) return;
    1329       22473 :   ind = 0;
    1330      145768 :   while (db)
    1331             :   {
    1332      123296 :     GEN c = Flx_rem(a,b, p);
    1333      123295 :     a = b; b = c; dc = degpol(c);
    1334      123295 :     if (dc < 0) break;
    1335             : 
    1336      123295 :     ind++;
    1337      123295 :     if (dc > dglist[ind]) dglist[ind] = dc;
    1338      123295 :     if (gc_needed(av,2))
    1339             :     {
    1340           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1341           0 :       (void)gc_all(av, 2, &a,&b);
    1342             :     }
    1343      123295 :     db = dc; /* = degpol(b) */
    1344             :   }
    1345       22472 :   if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
    1346       22472 :   set_avma(av);
    1347             : }
    1348             : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
    1349             :  * "generic" degree sequence as given by dglist, set *Ci and return
    1350             :  * resultant(a,b). Modular version of Collins's subresultant */
    1351             : static ulong
    1352     2089843 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
    1353             : {
    1354             :   long da,db,dc, ind;
    1355     2089843 :   ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
    1356     2089843 :   int s = 1;
    1357     2089843 :   pari_sp av = avma;
    1358             : 
    1359     2089843 :   *C0 = 1; *C1 = 0;
    1360     2089843 :   if (lgpol(a)==0 || lgpol(b)==0) return 0;
    1361     2080368 :   da = degpol(a);
    1362     2080368 :   db = degpol(b);
    1363     2080373 :   if (db > da)
    1364             :   {
    1365           0 :     swapspec(a,b, da,db);
    1366           0 :     if (both_odd(da,db)) s = -s;
    1367             :   }
    1368     2080373 :   else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
    1369     2080373 :   ind = 0;
    1370    19819316 :   while (db)
    1371             :   { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
    1372             :      * da = deg a, db = deg b */
    1373    17743680 :     GEN c = Flx_rem(a,b, p);
    1374    17651247 :     long delta = da - db;
    1375             : 
    1376    17651247 :     if (both_odd(da,db)) s = -s;
    1377    17649594 :     lb = Fl_mul(b[db+2], cb, p);
    1378    17678457 :     a = b; b = c; dc = degpol(c);
    1379    17686063 :     ind++;
    1380    17686063 :     if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
    1381    17681088 :     if (g == h)
    1382             :     { /* frequent */
    1383    17621118 :       ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
    1384    17681287 :       ca = cb;
    1385    17681287 :       cb = cc;
    1386             :     }
    1387             :     else
    1388             :     {
    1389       59970 :       ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
    1390       59970 :       ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
    1391       59970 :       ca = cb;
    1392       59970 :       cb = Fl_div(cc, ghdelta, p);
    1393             :     }
    1394    17739835 :     da = db; /* = degpol(a) */
    1395    17739835 :     db = dc; /* = degpol(b) */
    1396             : 
    1397    17739835 :     g = lb;
    1398    17739835 :     if (delta == 1)
    1399    17639160 :       h = g; /* frequent */
    1400             :     else
    1401      100675 :       h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
    1402             : 
    1403    17740351 :     if (gc_needed(av,2))
    1404             :     {
    1405           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
    1406           0 :       (void)gc_all(av, 2, &a,&b);
    1407             :     }
    1408             :   }
    1409     2075636 :   if (da > 1) return 0; /* Failure */
    1410             :   /* last nonconstant polynomial has degree 1 */
    1411     2075636 :   *C0 = Fl_mul(ca, a[2], p);
    1412     2075595 :   *C1 = Fl_mul(ca, a[3], p);
    1413     2075615 :   res = Fl_mul(cb, b[2], p);
    1414     2075615 :   if (s == -1) res = p - res;
    1415     2075615 :   return gc_ulong(av,res);
    1416             : }
    1417             : 
    1418             : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
    1419             :  * Return 0 in case of degree drop. */
    1420             : static GEN
    1421     2113703 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
    1422             : {
    1423             :   GEN z;
    1424     2113703 :   long i, lb = lg(Q);
    1425     2113703 :   ulong leadz = Flx_eval(leading_coeff(Q), x, p);
    1426     2113498 :   long vs=mael(Q,2,1);
    1427     2113498 :   if (!leadz) return zero_Flx(vs);
    1428             : 
    1429     2102838 :   z = cgetg(lb, t_VECSMALL); z[1] = vs;
    1430    20093387 :   for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
    1431     2102520 :   z[i] = leadz; return z;
    1432             : }
    1433             : 
    1434             : GEN
    1435        2072 : FpXY_FpXQ_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
    1436             : {
    1437        2072 :   pari_sp av = avma;
    1438        2072 :   long i, lb = lg(Q);
    1439             :   GEN z;
    1440        2072 :   if (lb == 2) return pol_0(vx);
    1441        2072 :   z = gel(Q, lb-1);
    1442        2072 :   if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
    1443             : 
    1444        2072 :   if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
    1445       48636 :   for (i=lb-2; i>=2; i--)
    1446             :   {
    1447       46564 :     GEN c = gel(Q,i);
    1448       46564 :     z = FqX_Fq_mul(z, y, T, p);
    1449       46564 :     z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
    1450             :   }
    1451        2072 :   return gc_upto(av, z);
    1452             : }
    1453             : 
    1454             : static GEN
    1455     1302814 : ZX_norml1(GEN x)
    1456             : {
    1457     1302814 :   long i, l = lg(x);
    1458             :   GEN s;
    1459             : 
    1460     1302814 :   if (l == 2) return gen_0;
    1461     1210232 :   s = gel(x, l-1); /* != 0 */
    1462     2719631 :   for (i = l-2; i > 1; i--) {
    1463     1509423 :     GEN xi = gel(x,i);
    1464     1509423 :     if (!signe(xi)) continue;
    1465     1224359 :     s = addii_sign(s,1, xi,1);
    1466             :   }
    1467     1210208 :   return s;
    1468             : }
    1469             : /* x >= 0, y != 0, return x + |y| */
    1470             : static GEN
    1471       25555 : addii_abs(GEN x, GEN y)
    1472             : {
    1473       25555 :   if (!signe(x)) return absi_shallow(y);
    1474       16043 :   return addii_sign(x,1, y,1);
    1475             : }
    1476             : 
    1477             : /* x a ZX, return sum_{i >= k} |x[i]| binomial(i, k) */
    1478             : static GEN
    1479       31644 : ZX_norml1_1(GEN x, long k)
    1480             : {
    1481       31644 :   long i, d = degpol(x);
    1482             :   GEN s, C; /* = binomial(i, k) */
    1483             : 
    1484       31644 :   if (!d || k > d) return gen_0;
    1485       31644 :   s = absi_shallow(gel(x, k+2)); /* may be 0 */
    1486       31644 :   C = gen_1;
    1487       68052 :   for (i = k+1; i <= d; i++) {
    1488       36403 :     GEN xi = gel(x,i+2);
    1489       36403 :     if (k) C = diviuexact(muliu(C, i), i-k);
    1490       36408 :     if (signe(xi)) s = addii_abs(s, mulii(C, xi));
    1491             :   }
    1492       31649 :   return s;
    1493             : }
    1494             : /* x has non-negative real coefficients */
    1495             : static GEN
    1496        3283 : RgX_norml1_1(GEN x, long k)
    1497             : {
    1498        3283 :   long i, d = degpol(x);
    1499             :   GEN s, C; /* = binomial(i, k) */
    1500             : 
    1501        3283 :   if (!d || k > d) return gen_0;
    1502        3283 :   s = gel(x, k+2); /* may be 0 */
    1503        3283 :   C = gen_1;
    1504        9198 :   for (i = k+1; i <= d; i++) {
    1505        5915 :     GEN xi = gel(x,i+2);
    1506        5915 :     if (k) C = diviuexact(muliu(C, i), i-k);
    1507        5915 :     if (!gequal0(xi)) s = gadd(s, gmul(C, xi));
    1508             :   }
    1509        3283 :   return s;
    1510             : }
    1511             : 
    1512             : /* N_2(A)^2 */
    1513             : static GEN
    1514        9019 : sqrN2(GEN A, long prec)
    1515             : {
    1516        9019 :   pari_sp av = avma;
    1517        9019 :   long i, l = lg(A);
    1518        9019 :   GEN a = gen_0;
    1519       43727 :   for (i = 2; i < l; i++)
    1520             :   {
    1521       34708 :     a = gadd(a, gabs(gnorm(gel(A,i)), prec));
    1522       34708 :     if (gc_needed(av,1))
    1523             :     {
    1524           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
    1525           0 :       a = gc_upto(av, a);
    1526             :     }
    1527             :   }
    1528        9019 :   return a;
    1529             : }
    1530             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1531             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1532             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1533             :  * Return e such that Res(A, B) < 2^e */
    1534             : static GEN
    1535        8165 : RgX_RgXY_ResBound(GEN A, GEN B, long prec)
    1536             : {
    1537        8165 :   pari_sp av = avma;
    1538        8165 :   GEN b = gen_0, bnd;
    1539        8165 :   long i, lB = lg(B);
    1540       31747 :   for (i=2; i<lB; i++)
    1541             :   {
    1542       23582 :     GEN t = gel(B,i);
    1543       23582 :     if (typ(t) == t_POL) t = gnorml1(t, prec);
    1544       23582 :     b = gadd(b, gabs(gsqr(t), prec));
    1545       23582 :     if (gc_needed(av,1))
    1546             :     {
    1547           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
    1548           0 :       b = gc_upto(av, b);
    1549             :     }
    1550             :   }
    1551        8165 :   bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
    1552             :                    gpowgs(b, degpol(A))), prec);
    1553        8165 :   return gc_upto(av, bnd);
    1554             : }
    1555             : /* A,B in C[X] return RgX_RgXY_ResBound(A, B(x+y)) */
    1556             : static GEN
    1557         854 : RgX_RgXY_ResBound_1(GEN A, GEN B, long prec)
    1558             : {
    1559         854 :   pari_sp av = avma, av2;
    1560         854 :   GEN b = gen_0, bnd;
    1561         854 :   long i, lB = lg(B);
    1562         854 :   B = shallowcopy(B);
    1563        4137 :   for (i=2; i<lB; i++) gel(B,i) = gabs(gel(B,i), prec);
    1564         854 :   av2 = avma;
    1565        4137 :   for (i=2; i<lB; i++)
    1566             :   {
    1567        3283 :     b = gadd(b, gsqr(RgX_norml1_1(B, i-2)));
    1568        3283 :     if (gc_needed(av2,1))
    1569             :     {
    1570           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
    1571           0 :       b = gc_upto(av2, b);
    1572             :     }
    1573             :   }
    1574         854 :   bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
    1575             :                    gpowgs(b, degpol(A))), prec);
    1576         854 :   return gc_upto(av, bnd);
    1577             : }
    1578             : 
    1579             : /* log2 N_2(A)^2 */
    1580             : static double
    1581      176992 : log2N2(GEN A)
    1582             : {
    1583      176992 :   pari_sp av = avma;
    1584      176992 :   long i, l = lg(A);
    1585      176992 :   GEN a = gen_0;
    1586     1337038 :   for (i=2; i < l; i++)
    1587             :   {
    1588     1160049 :     a = addii(a, sqri(gel(A,i)));
    1589     1160046 :     if (gc_needed(av,1))
    1590             :     {
    1591           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1592           0 :       a = gc_upto(av, a);
    1593             :     }
    1594             :   }
    1595      176989 :   return gc_double(av, dbllog2(a));
    1596             : }
    1597             : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
    1598             :  *   bound = N_2(A)^degpol B N_2(B)^degpol(A),  where
    1599             :  *     N_2(A) = sqrt(sum (N_1(Ai))^2)
    1600             :  * Return e such that Res(A, B) < 2^e */
    1601             : 
    1602             : static double
    1603     2188065 : resbound(GEN B)
    1604             : {
    1605     2188065 :   pari_sp av = avma;
    1606     2188065 :   long i, lB = lg(B);
    1607     2188065 :   GEN b = gen_0;
    1608     9743495 :   for (i=2; i<lB; i++)
    1609             :   {
    1610     7555580 :     GEN t = gel(B,i);
    1611     7555580 :     if (typ(t) == t_POL) t = ZX_norml1(t);
    1612     7555571 :     b = addii(b, sqri(t));
    1613     7555447 :     if (gc_needed(av,1))
    1614             :     {
    1615           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1616           0 :       b = gc_upto(av, b);
    1617             :     }
    1618             :   }
    1619     2187915 :   return gc_double(av, dbllog2(b));
    1620             : }
    1621             : 
    1622             : ulong
    1623      166910 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
    1624             : {
    1625      166910 :   pari_sp av = avma;
    1626      166910 :   double logb = resbound(B);
    1627             :   long i;
    1628      166908 :   if (dB) logb -= 2 * dbllog2(dB);
    1629      166908 :   i = (long)((degpol(B) * log2N2(A) + degpol(A) * logb) / 2);
    1630      166909 :   return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
    1631             : }
    1632             : static ulong
    1633     1010585 : ZXX_ResBound(GEN A, GEN B)
    1634             : {
    1635     1010585 :   pari_sp av = avma;
    1636     1010585 :   double loga  = resbound(A), logb = resbound(B);
    1637     1010576 :   long i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
    1638     1010584 :   return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
    1639             : }
    1640             : 
    1641             : /* A,B ZX. Return ZX_ZXY_ResBound(A(x), B(x+y)) */
    1642             : static ulong
    1643       10084 : ZX_ZXY_ResBound_1(GEN A, GEN B)
    1644             : {
    1645       10084 :   pari_sp av = avma;
    1646       10084 :   GEN b = gen_0;
    1647       10084 :   long i, lB = lg(B);
    1648       41732 :   for (i=2; i<lB; i++)
    1649             :   {
    1650       31644 :     b = addii(b, sqri(ZX_norml1_1(B, i-2)));
    1651       31648 :     if (gc_needed(av,1))
    1652             :     {
    1653           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
    1654           0 :       b = gc_upto(av, b);
    1655             :     }
    1656             :   }
    1657       10088 :   i = (long)((degpol(B) * log2N2(A) + degpol(A) * dbllog2(b)) / 2);
    1658       10085 :   return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
    1659             : }
    1660             : /* special case B = A' */
    1661             : static ulong
    1662     1134901 : ZX_discbound(GEN A)
    1663             : {
    1664     1134901 :   pari_sp av = avma;
    1665     1134901 :   GEN a = gen_0, b = gen_0;
    1666     1134901 :   long i , lA = lg(A), dA = degpol(A);
    1667             :   double loga, logb;
    1668     6771564 :   for (i = 2; i < lA; i++)
    1669             :   {
    1670     5636952 :     GEN c = sqri(gel(A,i));
    1671     5636651 :     a = addii(a, c);
    1672     5636626 :     if (i > 2) b = addii(b, mulii(c, sqru(i-2)));
    1673     5636611 :     if (gc_needed(av,1))
    1674             :     {
    1675           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZX_discbound i = %ld",i);
    1676           0 :       (void)gc_all(av, 2, &a, &b);
    1677             :     }
    1678             :   }
    1679     1134612 :   loga = dbllog2(a);
    1680     1134742 :   logb = dbllog2(b); set_avma(av);
    1681     1134787 :   i = (long)(((dA-1) * loga + dA * logb) / 2);
    1682     1134787 :   return (i <= 0)? 1: 1 + (ulong)i;
    1683             : }
    1684             : 
    1685             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1686             : static ulong
    1687     5539714 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong pi, ulong la)
    1688             : {
    1689     5539714 :   GEN ev = FlxY_evalx_pre(b, n, p, pi);
    1690     5540437 :   long drop = lg(b) - lg(ev);
    1691     5540437 :   ulong r = Flx_resultant_pre(a, ev, p, pi);
    1692     5539457 :   if (drop && la != 1) r = Fl_mul(r, Fl_powu_pre(la, drop, p, pi), p);
    1693     5539499 :   return r;
    1694             : }
    1695             : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
    1696             : static ulong
    1697     4347846 : FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong pi, GEN la, GEN lb)
    1698             : {
    1699     4347846 :   GEN av = FlxY_evalx_pre(a, n, p, pi);
    1700     4347841 :   GEN bv = FlxY_evalx_pre(b, n, p, pi);
    1701     4347845 :   long lav = lgpol(av), lbv = lgpol(bv), dropa, dropb;
    1702             :   ulong r;
    1703     4347836 :   if (lav==0 && lbv==0) return 1UL;
    1704     4347829 :   dropa = lgpol(a) - lav;
    1705     4347825 :   dropb = lgpol(b) - lbv;
    1706     4347825 :   r = Flx_resultant_pre(av, bv, p, pi);
    1707     4347828 :   if (dropa)
    1708             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    1709          14 :     ulong c = Flx_eval_pre(lb, n,p ,pi); /* lc(B) */
    1710          14 :     if (odd(degpol(b))) c = p - c;
    1711          14 :     c = Fl_powu(c, dropa, p);
    1712          14 :     if (c != 1UL) r = Fl_mul(r, c, p);
    1713             :   }
    1714     4347814 :   else if (dropb)
    1715             :   { /* multiply by lc(A)^(deg B - deg b) */
    1716           7 :     ulong c = Flx_eval_pre(la, n,p ,pi); /* lc(B) */
    1717           7 :     c = Fl_powu(c, dropb, p);
    1718           7 :     if (c != 1UL) r = Fl_mul(r, c, p);
    1719             :   }
    1720     4347828 :   return r;
    1721             : }
    1722             : 
    1723             : static GEN
    1724         284 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
    1725             : {
    1726         284 :   GEN ev = FpXY_evaly(b, n, p, vX);
    1727         284 :   long drop = db-degpol(ev);
    1728         284 :   GEN r = FpX_resultant(a, ev, p);
    1729         284 :   if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
    1730         284 :   return r;
    1731             : }
    1732             : 
    1733             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1734             : /* Return a Fly */
    1735             : static GEN
    1736      368515 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, ulong pi, long dres, long sx)
    1737             : {
    1738             :   long i;
    1739      368515 :   ulong n, la = Flx_lead(a);
    1740      368515 :   GEN  x = cgetg(dres+2, t_VECSMALL);
    1741      368513 :   GEN  y = cgetg(dres+2, t_VECSMALL);
    1742             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1743             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1744     2957624 :   for (i=0,n = 1; i < dres; n++)
    1745             :   {
    1746     2589108 :     x[++i] = n;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
    1747     2588987 :     x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
    1748             :   }
    1749      368516 :   if (i == dres)
    1750             :   {
    1751      363010 :     x[++i] = 0;   y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
    1752             :   }
    1753      368517 :   return Flv_polint(x,y, p, sx);
    1754             : }
    1755             : 
    1756             : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
    1757             : /* Return a Fly */
    1758             : static GEN
    1759     1118028 : FlxX_resultant_polint(GEN a, GEN b, ulong p, ulong pi, long dres, long sx)
    1760             : {
    1761             :   long i;
    1762             :   ulong n;
    1763     1118028 :   GEN la = leading_coeff(a), lb = leading_coeff(b);
    1764     1118027 :   GEN x = cgetg(dres+2, t_VECSMALL);
    1765     1118028 :   GEN y = cgetg(dres+2, t_VECSMALL);
    1766             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1767             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1768     2998275 :   for (i=0,n = 1; i < dres; n++)
    1769             :   {
    1770     1880244 :     x[++i] = n;   y[i] = FlxY_eval_resultant(a,b, x[i], p,pi, la,lb);
    1771     1880242 :     x[++i] = p-n; y[i] = FlxY_eval_resultant(a,b, x[i], p,pi, la,lb);
    1772             :   }
    1773     1118031 :   if (i == dres)
    1774             :   {
    1775      587380 :     x[++i] = 0;   y[i] = FlxY_eval_resultant(a,b, x[i], p,pi, la,lb);
    1776             :   }
    1777     1118033 :   return Flv_polint(x,y, p, sx);
    1778             : }
    1779             : 
    1780             : static GEN
    1781        7498 : FlxX_pseudorem(GEN x, GEN y, ulong p, ulong pi)
    1782             : {
    1783        7498 :   long vx = varn(x), dx, dy, dz, i, lx, dp;
    1784        7498 :   pari_sp av = avma, av2;
    1785             : 
    1786        7498 :   if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
    1787        7498 :   (void)new_chunk(2);
    1788        7499 :   dx=degpol(x); x = RgX_recip_i(x)+2;
    1789        7500 :   dy=degpol(y); y = RgX_recip_i(y)+2; dz=dx-dy; dp = dz+1;
    1790        7500 :   av2 = avma;
    1791             :   for (;;)
    1792             :   {
    1793       61446 :     gel(x,0) = Flx_neg(gel(x,0), p); dp--;
    1794      230122 :     for (i=1; i<=dy; i++)
    1795      168444 :       gel(x,i) = Flx_add( Flx_mul_pre(gel(y,0), gel(x,i), p, pi),
    1796      168656 :                           Flx_mul_pre(gel(x,0), gel(y,i), p, pi), p );
    1797     1112232 :     for (   ; i<=dx; i++)
    1798     1055286 :       gel(x,i) = Flx_mul_pre(gel(y,0), gel(x,i), p, pi);
    1799       65312 :     do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
    1800       56946 :     if (dx < dy) break;
    1801       49446 :     if (gc_needed(av2,1))
    1802             :     {
    1803           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
    1804           0 :       gc_slice(av2,x,dx+1);
    1805             :     }
    1806             :   }
    1807        7500 :   if (dx < 0) return zero_Flx(0);
    1808        7500 :   lx = dx+3; x -= 2;
    1809        7500 :   x[0]=evaltyp(t_POL) | _evallg(lx);
    1810        7500 :   x[1]=evalsigne(1) | evalvarn(vx);
    1811        7500 :   x = RgX_recip_i(x);
    1812        7501 :   if (dp)
    1813             :   { /* multiply by y[0]^dp   [beware dummy vars from FpX_FpXY_resultant] */
    1814        1959 :     GEN t = Flx_powu_pre(gel(y,0), dp, p, pi);
    1815        7837 :     for (i=2; i<lx; i++) gel(x,i) = Flx_mul_pre(gel(x,i), t, p, pi);
    1816             :   }
    1817        7500 :   return gc_GEN(av, x);
    1818             : }
    1819             : 
    1820             : static GEN
    1821        1973 : Flx_Lazard(GEN x, GEN y, long n, ulong p, ulong pi)
    1822             : {
    1823             :   long a;
    1824             :   GEN c;
    1825        1973 :   if (n == 1) return x;
    1826        1973 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    1827        1973 :   c=x; n-=a;
    1828        1973 :   y = Flx_get_red_pre(y, p, pi);
    1829        8830 :   while (a>1)
    1830             :   {
    1831        6857 :     a>>=1; c = Flx_div_pre(Flx_sqr_pre(c, p, pi), y, p, pi);
    1832        6856 :     if (n>=a) { c = Flx_div_pre(Flx_mul_pre(c,x,p,pi),y, p,pi); n -= a; }
    1833             :   }
    1834        1973 :   return c;
    1835             : }
    1836             : 
    1837             : /* return a Flx */
    1838             : static GEN
    1839        2510 : FlxX_resultant_subres(GEN u, GEN v, ulong p, ulong pi, long sx)
    1840             : {
    1841        2510 :   pari_sp av = avma, av2;
    1842             :   long degq, dx, dy, du, dv, dr, signh;
    1843             :   GEN z, g, h, r, p1;
    1844             : 
    1845        2510 :   dx = degpol(u); dy = degpol(v); signh = 1;
    1846        2510 :   if (dx < dy)
    1847             :   {
    1848           7 :     swap(u,v); lswap(dx,dy);
    1849           7 :     if (both_odd(dx, dy)) signh = -signh;
    1850             :   }
    1851        2510 :   if (dy < 0) return zero_Flx(sx);
    1852        2510 :   if (dy==0) return gc_upto(av, Flx_powu_pre(gel(v,2),dx,p,pi));
    1853             : 
    1854        2510 :   g = h = pol1_Flx(sx); av2 = avma;
    1855             :   for(;;)
    1856             :   {
    1857        7501 :     r = FlxX_pseudorem(u,v,p,pi); dr = lg(r);
    1858        7501 :     if (dr == 2) { set_avma(av); return zero_Flx(sx); }
    1859        7501 :     du = degpol(u); dv = degpol(v); degq = du-dv;
    1860        7500 :     u = v; p1 = g; g = leading_coeff(u);
    1861        7500 :     switch(degq)
    1862             :     {
    1863           0 :       case 0: break;
    1864        5527 :       case 1:
    1865        5527 :         p1 = Flx_mul_pre(h,p1, p, pi); h = g; break;
    1866        1973 :       default:
    1867        1973 :         p1 = Flx_mul_pre(Flx_powu_pre(h,degq,p,pi), p1, p, pi);
    1868        1973 :         h = Flx_Lazard(g,h,degq, p, pi);
    1869             :     }
    1870        7499 :     if (both_odd(du,dv)) signh = -signh;
    1871        7499 :     v = FlxY_Flx_div_pre(r, p1, p, pi);
    1872        7499 :     if (dr==3) break;
    1873        4989 :     if (gc_needed(av2,1))
    1874             :     {
    1875           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
    1876           0 :       (void)gc_all(av2,4, &u, &v, &g, &h);
    1877             :     }
    1878             :   }
    1879        2510 :   z = gel(v,2);
    1880        2510 :   if (dv > 1) z = Flx_Lazard(z,h,dv,p,pi);
    1881        2510 :   if (signh < 0) z = Flx_neg(z,p);
    1882        2510 :   return gc_upto(av, z);
    1883             : }
    1884             : 
    1885             : /* Return a Flx*/
    1886             : GEN
    1887           0 : FlxX_resultant_pre(GEN a, GEN b, ulong p, ulong pi, long sx)
    1888             : {
    1889           0 :   pari_sp ltop = avma;
    1890           0 :   long da = degpol(a), db = degpol(b);
    1891             :   ulong dres;
    1892             :   GEN z;
    1893           0 :   if (da<0 || db<0) return pol0_Flx(sx);
    1894           0 :   dres = FlxY_degreex(a)*db+FlxY_degreex(b)*da;
    1895           0 :   z = dres >= p ? FlxX_resultant_subres(a, b, p, pi, sx)
    1896           0 :                 : FlxX_resultant_polint(a, b, p, pi, dres, sx);
    1897           0 :   return gc_upto(ltop, z);
    1898             : }
    1899             : 
    1900             : GEN
    1901           0 : FlxX_resultant(GEN a, GEN b, ulong p, long sx)
    1902           0 : { return FlxX_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), sx); }
    1903             : 
    1904             : /* Warning:
    1905             :  * This function switches between valid and invalid variable ordering*/
    1906             : static GEN
    1907        6120 : FlxY_to_FlyX(GEN b, long sv)
    1908             : {
    1909        6120 :   long i, n=-1;
    1910        6120 :   long sw = b[1]&VARNBITS;
    1911       20877 :   for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
    1912        6120 :   return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
    1913             : }
    1914             : 
    1915             : /* Return a Fly*/
    1916             : GEN
    1917        6120 : Flx_FlxY_resultant(GEN a, GEN b, ulong p)
    1918             : {
    1919        6120 :   pari_sp ltop=avma;
    1920        6120 :   long dres = degpol(a)*degpol(b);
    1921        6120 :   long sx=a[1], sy=b[1]&VARNBITS;
    1922        6120 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    1923             :   GEN z;
    1924        6120 :   b = FlxY_to_FlyX(b,sx);
    1925        6120 :   if ((ulong)dres >= p)
    1926        2510 :     z = FlxX_resultant_subres(Fly_to_FlxY(a, sy), b, p, pi, sy);
    1927             :   else
    1928        3610 :     z = Flx_FlxY_resultant_polint(a, b, p, pi, (ulong)dres, sy);
    1929        6120 :   return gc_upto(ltop,z);
    1930             : }
    1931             : 
    1932             : /* Return a t_POL in variable vc whose coeffs are the coeffs of b in
    1933             :  * variable v; vc must have higher priority than all variables occuring in b. */
    1934             : GEN
    1935      146885 : swap_vars(GEN b, long v, long vc)
    1936             : {
    1937      146885 :   long i, n = RgX_degree(b, v);
    1938             :   GEN c, x;
    1939      146884 :   if (n < 0) return pol_0(vc);
    1940      146884 :   c = cgetg(n+3, t_POL); x = c + 2;
    1941      146884 :   c[1] = evalsigne(1) | evalvarn(vc);
    1942      971237 :   for (i = 0; i <= n; i++) gel(x,i) = polcoef_i(b, i, v);
    1943      146884 :   return c;
    1944             : }
    1945             : 
    1946             : /* assume varn(b) << varn(a) */
    1947             : /* return a FpY*/
    1948             : GEN
    1949          15 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
    1950             : {
    1951          15 :   long i,n,dres, db, vY = varn(b), vX = varn(a);
    1952             :   GEN la,x,y;
    1953             : 
    1954          15 :   if (lgefint(p) == 3)
    1955             :   {
    1956           0 :     ulong pp = uel(p,2);
    1957           0 :     b = ZXX_to_FlxX(b, pp, vX);
    1958           0 :     a = ZX_to_Flx(a, pp);
    1959           0 :     x = Flx_FlxY_resultant(a, b, pp);
    1960           0 :     return Flx_to_ZX(x);
    1961             :   }
    1962          15 :   db = RgXY_degreex(b);
    1963          15 :   dres = degpol(a)*degpol(b);
    1964          15 :   la = leading_coeff(a);
    1965          15 :   x = cgetg(dres+2, t_VEC);
    1966          15 :   y = cgetg(dres+2, t_VEC);
    1967             :  /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
    1968             :   * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
    1969         157 :   for (i=0,n = 1; i < dres; n++)
    1970             :   {
    1971         142 :     gel(x,++i) = utoipos(n);
    1972         142 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1973         142 :     gel(x,++i) = subiu(p,n);
    1974         142 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
    1975             :   }
    1976          15 :   if (i == dres)
    1977             :   {
    1978           0 :     gel(x,++i) = gen_0;
    1979           0 :     gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
    1980             :   }
    1981          15 :   return FpV_polint(x,y, p, vY);
    1982             : }
    1983             : 
    1984             : GEN
    1985         191 : FpX_composedsum(GEN P, GEN Q, GEN p)
    1986             : {
    1987         191 :   pari_sp av = avma;
    1988         191 :   if (lgefint(p)==3)
    1989             :   {
    1990           0 :     ulong pp = p[2];
    1991           0 :     GEN z = Flx_composedsum(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
    1992           0 :     return gc_upto(av, Flx_to_ZX(z));
    1993             :   }
    1994             :   else
    1995             :   {
    1996         191 :     long n = 1+ degpol(P)*degpol(Q);
    1997         191 :     GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
    1998         191 :     GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
    1999         191 :     GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
    2000         191 :     GEN lead = Fp_mul(Fp_powu(leading_coeff(P),degpol(Q), p),
    2001         191 :         Fp_powu(leading_coeff(Q),degpol(P), p), p);
    2002         191 :     GEN R = FpX_fromNewton(L, p);
    2003         191 :     return gc_upto(av, FpX_Fp_mul(R, lead, p));
    2004             :   }
    2005             : }
    2006             : 
    2007             : GEN
    2008           0 : FpX_composedprod(GEN P, GEN Q, GEN p)
    2009             : {
    2010           0 :   pari_sp av = avma;
    2011           0 :   if (lgefint(p)==3)
    2012             :   {
    2013           0 :     ulong pp = p[2];
    2014           0 :     GEN z = Flx_composedprod(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
    2015           0 :     return gc_upto(av, Flx_to_ZX(z));
    2016             :   }
    2017             :   else
    2018             :   {
    2019           0 :     long n = 1+ degpol(P)*degpol(Q);
    2020           0 :     GEN L = FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
    2021           0 :     return gc_upto(av,FpX_fromNewton(L, p));
    2022             :   }
    2023             : }
    2024             : 
    2025             : static GEN
    2026         191 : _FpX_composedsum(void *E, GEN a, GEN b)
    2027         191 : { return FpX_composedsum(a,b, (GEN)E); }
    2028             : 
    2029             : GEN
    2030        1637 : FpXV_composedsum(GEN V, GEN p)
    2031             : {
    2032        1637 :   if (lgefint(p)==3)
    2033             :   {
    2034           0 :     ulong pp = p[2];
    2035           0 :     return Flx_to_ZX(FlxV_composedsum(ZXV_to_FlxV(V, pp), pp));
    2036             :   }
    2037        1637 :   return gen_product(V, (void *)p, &_FpX_composedsum);
    2038             : }
    2039             : 
    2040             : /* 0, 1, -1, 2, -2, ... */
    2041             : #define next_lambda(a) (a>0 ? -a : 1-a)
    2042             : 
    2043             : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
    2044             :  * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
    2045             :  * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
    2046             :  * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
    2047             :  * the Last nonconstant polynomial in the Euclidean Remainder Sequence */
    2048             : static GEN
    2049       22421 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
    2050             : {
    2051             :   ulong bound, dp;
    2052       22421 :   pari_sp av = avma, av2 = 0;
    2053       22421 :   long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
    2054             :   long stable, checksqfree, i,n, cnt, degB;
    2055       22421 :   long v, vX = varn(B0), vY = varn(A); /* vY < vX */
    2056             :   GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
    2057             :   forprime_t S;
    2058             : 
    2059       22421 :   if (degA == 1)
    2060             :   {
    2061        1260 :     GEN a1 = gel(A,3), a0 = gel(A,2);
    2062        1260 :     B = lambda? RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
    2063        1260 :     H = gsubst(B, vY, gdiv(gneg(a0),a1));
    2064        1260 :    if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
    2065        1260 :     *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
    2066        1260 :     return gc_all(av, 2, &H, LERS);
    2067             :   }
    2068             : 
    2069       21161 :   dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
    2070       21161 :   C0 = cgetg(dres+2, t_VECSMALL);
    2071       21161 :   C1 = cgetg(dres+2, t_VECSMALL);
    2072       21161 :   dglist = cgetg(dres+1, t_VECSMALL);
    2073       21161 :   x = cgetg(dres+2, t_VECSMALL);
    2074       21161 :   y = cgetg(dres+2, t_VECSMALL);
    2075       21161 :   B0 = leafcopy(B0);
    2076       21161 :   A = leafcopy(A);
    2077       21161 :   B = B0;
    2078       21161 :   v = fetch_var_higher(); setvarn(A,v);
    2079             :   /* make sure p large enough */
    2080       22320 : INIT:
    2081             :   /* always except the first time */
    2082       22320 :   if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
    2083       22320 :   if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
    2084       22320 :   B = swap_vars(B, vY, v);
    2085             :   /* B0(lambda v + x, v) */
    2086       22320 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    2087       22320 :   av2 = avma;
    2088             : 
    2089       22320 :   if (degA <= 3)
    2090             :   { /* sub-resultant faster for small degrees */
    2091       11319 :     H = RgX_resultant_all(A,B,&q);
    2092       11319 :     if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
    2093       10374 :     H0 = gel(q,2);
    2094       10374 :     if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
    2095       10374 :     H1 = gel(q,3);
    2096       10374 :     if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
    2097       10374 :     if (!ZX_is_squarefree(H)) goto INIT;
    2098       10332 :     goto END;
    2099             :   }
    2100             : 
    2101       11001 :   H = H0 = H1 = NULL;
    2102       11001 :   degB = degpol(B);
    2103       11001 :   bound = ZX_ZXY_ResBound(A, B, NULL);
    2104       11001 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    2105       11001 :   dp = 1;
    2106       11001 :   init_modular_big(&S);
    2107       11001 :   for(cnt = 0, checksqfree = 1;;)
    2108       49460 :   {
    2109       60461 :     ulong p = u_forprime_next(&S);
    2110             :     GEN Hi;
    2111       60461 :     a = ZX_to_Flx(A, p);
    2112       60461 :     b = ZXX_to_FlxX(B, p, varn(A));
    2113       60461 :     if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
    2114       60461 :     if (checksqfree)
    2115             :     { /* find degree list for generic Euclidean Remainder Sequence */
    2116       11001 :       long goal = minss(degpol(a), degpol(b)); /* longest possible */
    2117       74267 :       for (n=1; n <= goal; n++) dglist[n] = 0;
    2118       11001 :       setlg(dglist, 1);
    2119       24154 :       for (n=0; n <= dres; n++)
    2120             :       {
    2121       23741 :         ev = FlxY_evalx_drop(b, n, p);
    2122       23740 :         Flx_resultant_set_dglist(a, ev, dglist, p);
    2123       23740 :         if (lg(dglist)-1 == goal) break;
    2124             :       }
    2125             :       /* last pol in ERS has degree > 1 ? */
    2126       11000 :       goal = lg(dglist)-1;
    2127       11000 :       if (degpol(B) == 1) { if (!goal) goto INIT; }
    2128             :       else
    2129             :       {
    2130       10938 :         if (goal <= 1) goto INIT;
    2131       10854 :         if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
    2132             :       }
    2133       10916 :       if (DEBUGLEVEL>4)
    2134           0 :         err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
    2135             :     }
    2136             : 
    2137     2150352 :     for (i=0,n = 0; i <= dres; n++)
    2138             :     {
    2139     2089976 :       ev = FlxY_evalx_drop(b, n, p);
    2140     2089817 :       x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
    2141     2089975 :       if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
    2142             :     }
    2143       60376 :     Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
    2144       60377 :     Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
    2145       60377 :     if (!H && degpol(Hp) != dres) continue;
    2146       60377 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2147       60377 :     if (checksqfree) {
    2148       10917 :       if (!Flx_is_squarefree(Hp, p)) goto INIT;
    2149       10829 :       if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    2150       10829 :       checksqfree = 0;
    2151             :     }
    2152             : 
    2153       60289 :     if (!H)
    2154             :     { /* initialize */
    2155       10829 :       q = utoipos(p); stable = 0;
    2156       10829 :       H = ZX_init_CRT(Hp, p,vX);
    2157       10829 :       H0= ZX_init_CRT(H0p, p,vX);
    2158       10829 :       H1= ZX_init_CRT(H1p, p,vX);
    2159             :     }
    2160             :     else
    2161             :     {
    2162       49460 :       GEN qp = muliu(q,p);
    2163       49460 :       stable  = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
    2164       49460 :               & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
    2165       49459 :               & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
    2166       49460 :       q = qp;
    2167             :     }
    2168             :     /* could make it probabilistic for H ? [e.g if stable twice, etc]
    2169             :      * Probabilistic anyway for H0, H1 */
    2170       60289 :     if (DEBUGLEVEL>5 && (stable ||  ++cnt==100))
    2171           0 :     { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
    2172       60289 :     if (stable && (ulong)expi(q) >= bound) break; /* DONE */
    2173       49460 :     if (gc_needed(av,2))
    2174             :     {
    2175           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
    2176           0 :       (void)gc_all(av2, 4, &H, &q, &H0, &H1);
    2177             :     }
    2178             :   }
    2179       21161 : END:
    2180       21161 :   if (DEBUGLEVEL>5) err_printf(" done\n");
    2181       21161 :   setvarn(H, vX); (void)delete_var();
    2182       21161 :   *LERS = mkvec2(H0,H1);
    2183       21161 :   *plambda = lambda; return gc_all(av, 2, &H, LERS);
    2184             : }
    2185             : 
    2186             : GEN
    2187       60214 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
    2188             : {
    2189       60214 :   if (LERS)
    2190             :   {
    2191       22421 :     if (!plambda)
    2192           0 :       pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
    2193       22421 :     return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
    2194             :   }
    2195       37793 :   return ZX_ZXY_rnfequation(A, B, plambda);
    2196             : }
    2197             : 
    2198             : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
    2199             :  * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
    2200             :  * squarefree */
    2201             : GEN
    2202       22595 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
    2203             : {
    2204       22595 :   pari_sp av = avma;
    2205             :   GEN R, a;
    2206             :   long dA;
    2207             :   int delvar;
    2208             : 
    2209       22595 :   if (v < 0) v = 0;
    2210       22595 :   switch (typ(A))
    2211             :   {
    2212       22595 :     case t_POL: dA = degpol(A); if (dA > 0) break;
    2213           0 :       A = constant_coeff(A);
    2214           0 :     default:
    2215           0 :       if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
    2216           0 :       return gc_upto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
    2217             :   }
    2218       22595 :   delvar = 0;
    2219       22595 :   if (varncmp(varn(T), 0) <= 0)
    2220             :   {
    2221        3681 :     long v0 = fetch_var(); delvar = 1;
    2222        3681 :     T = leafcopy(T); setvarn(T,v0);
    2223        3681 :     A = leafcopy(A); setvarn(A,v0);
    2224             :   }
    2225       22595 :   R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
    2226       22595 :   if (delvar) (void)delete_var();
    2227       22595 :   setvarn(R, v); a = leading_coeff(T);
    2228       22595 :   if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
    2229       22595 :   return gc_upto(av, R);
    2230             : }
    2231             : 
    2232             : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
    2233             : GEN
    2234     1247424 : ZXQ_charpoly(GEN A, GEN T, long v)
    2235             : {
    2236     1247424 :   return (degpol(T) < 16) ? RgXQ_charpoly_i(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
    2237             : }
    2238             : 
    2239             : GEN
    2240        9772 : QXQ_charpoly(GEN A, GEN T, long v)
    2241             : {
    2242        9772 :   pari_sp av = avma;
    2243        9772 :   GEN den, B = Q_remove_denom(A, &den);
    2244        9772 :   GEN P = ZXQ_charpoly(B, T, v);
    2245        9772 :   return gc_GEN(av, den ? RgX_rescale(P, ginv(den)): P);
    2246             : }
    2247             : 
    2248             : static ulong
    2249     3864951 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
    2250             : {
    2251     3864951 :   pari_sp av = avma;
    2252     3864951 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2253             :   ulong H, dp;
    2254     3864851 :   if (dropa && dropb) return 0; /* p | lc(A), p | lc(B) */
    2255     3864851 :   H = Flx_resultant(a, b, p);
    2256     3864613 :   if (dropa)
    2257             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2258           0 :     ulong c = b[degB+2]; /* lc(B) */
    2259           0 :     if (odd(degB)) c = p - c;
    2260           0 :     c = Fl_powu(c, dropa, p);
    2261           0 :     if (c != 1) H = Fl_mul(H, c, p);
    2262             :   }
    2263     3864613 :   else if (dropb)
    2264             :   { /* multiply by lc(A)^(deg B - deg b) */
    2265           0 :     ulong c = a[degA+2]; /* lc(A) */
    2266           0 :     c = Fl_powu(c, dropb, p);
    2267           0 :     if (c != 1) H = Fl_mul(H, c, p);
    2268             :   }
    2269     3864615 :   dp = dB ? umodiu(dB, p): 1;
    2270     3864615 :   if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
    2271     3864616 :   return gc_ulong(av, H);
    2272             : }
    2273             : 
    2274             : /* If B=NULL, assume B=A' */
    2275             : static GEN
    2276     1495325 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
    2277             : {
    2278     1495325 :   pari_sp av = avma, av2;
    2279     1495325 :   long degA, degB, i, n = lg(P)-1;
    2280             :   GEN H, T;
    2281             : 
    2282     1495325 :   degA = degpol(A);
    2283     1495319 :   degB = B? degpol(B): degA - 1;
    2284     1495320 :   if (n == 1)
    2285             :   {
    2286      811323 :     ulong Hp, p = uel(P,1);
    2287      811323 :     GEN a = ZX_to_Flx(A, p), b = B? ZX_to_Flx(B, p): Flx_deriv(a, p);
    2288      811306 :     Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
    2289      811315 :     set_avma(av); *mod = utoipos(p); return utoi(Hp);
    2290             :   }
    2291      683997 :   T = ZV_producttree(P);
    2292      683997 :   A = ZX_nv_mod_tree(A, P, T);
    2293      683991 :   if (B) B = ZX_nv_mod_tree(B, P, T);
    2294      683991 :   H = cgetg(n+1, t_VECSMALL); av2 = avma;
    2295     3737332 :   for(i=1; i <= n; i++, set_avma(av2))
    2296             :   {
    2297     3053346 :     ulong p = P[i];
    2298     3053346 :     GEN a = gel(A,i), b = B? gel(B,i): Flx_deriv(a, p);
    2299     3053653 :     H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
    2300             :   }
    2301      683986 :   H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
    2302      683994 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2303             : }
    2304             : 
    2305             : GEN
    2306     1495329 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
    2307             : {
    2308     1495329 :   GEN V = cgetg(3, t_VEC);
    2309     1495325 :   if (typ(B) == t_INT) B = NULL;
    2310     1495325 :   if (!signe(dB)) dB = NULL;
    2311     1495325 :   gel(V,1) = ZX_resultant_slice(A, B, dB, P, &gel(V,2));
    2312     1495309 :   return V;
    2313             : }
    2314             : 
    2315             : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
    2316             :  * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
    2317             : GEN
    2318     1350821 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
    2319             : {
    2320     1350821 :   pari_sp av = avma;
    2321             :   forprime_t S;
    2322             :   GEN  H, worker;
    2323     1350821 :   if (!B && degpol(A)==2)
    2324             :   {
    2325      114096 :     GEN a = gel(A,4), b = gel(A,3), c = gel(A,2);
    2326      114096 :     H = mulii(a, subii(shifti(mulii(a, c), 2), sqri(b)));
    2327      114088 :     if (dB) H = diviiexact(H, sqri(dB));
    2328      114088 :     return gc_INT(av, H);
    2329             :   }
    2330     1236720 :   if (B)
    2331             :   {
    2332      153991 :     long a = degpol(A), b = degpol(B);
    2333      153991 :     if (a < 0 || b < 0) return gen_0;
    2334      153961 :     if (!a) return powiu(gel(A,2), b);
    2335      153961 :     if (!b) return powiu(gel(B,2), a);
    2336      153336 :     if (minss(a, b) <= 1)
    2337             :     {
    2338       76624 :       H = RgX_resultant_all(A, B, NULL);
    2339       76624 :       if (dB) H = diviiexact(H, powiu(dB, a));
    2340       76624 :       return gc_INT(av, H);
    2341             :     }
    2342       76712 :     if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
    2343             :   }
    2344     1159449 :   worker = snm_closure(is_entry("_ZX_resultant_worker"),
    2345             :                        mkvec3(A, B? B: gen_0, dB? dB: gen_0));
    2346     1159577 :   init_modular_big(&S);
    2347     1159527 :   H = gen_crt("ZX_resultant_all", worker, &S, dB, bound, 0, NULL,
    2348             :               ZV_chinese_center, Fp_center);
    2349     1159542 :   return gc_INT(av, H);
    2350             : }
    2351             : 
    2352             : /* A0 and B0 in Q[X] */
    2353             : GEN
    2354          56 : QX_resultant(GEN A0, GEN B0)
    2355             : {
    2356             :   GEN s, a, b, A, B;
    2357          56 :   pari_sp av = avma;
    2358             : 
    2359          56 :   A = Q_primitive_part(A0, &a);
    2360          56 :   B = Q_primitive_part(B0, &b);
    2361          56 :   s = ZX_resultant(A, B);
    2362          56 :   if (!signe(s)) { set_avma(av); return gen_0; }
    2363          56 :   if (a) s = gmul(s, gpowgs(a,degpol(B)));
    2364          56 :   if (b) s = gmul(s, gpowgs(b,degpol(A)));
    2365          56 :   return gc_upto(av, s);
    2366             : }
    2367             : 
    2368             : GEN
    2369       56077 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
    2370             : 
    2371             : GEN
    2372           0 : QXQ_intnorm(GEN A, GEN B)
    2373             : {
    2374             :   GEN c, n, R, lB;
    2375           0 :   long dA = degpol(A), dB = degpol(B);
    2376           0 :   pari_sp av = avma;
    2377           0 :   if (dA < 0) return gen_0;
    2378           0 :   A = Q_primitive_part(A, &c);
    2379           0 :   if (!c || typ(c) == t_INT) {
    2380           0 :     n = c;
    2381           0 :     R = ZX_resultant(B, A);
    2382             :   } else {
    2383           0 :     n = gel(c,1);
    2384           0 :     R = ZX_resultant_all(B, A, gel(c,2), 0);
    2385             :   }
    2386           0 :   if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
    2387           0 :   lB = leading_coeff(B);
    2388           0 :   if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
    2389           0 :   return gc_INT(av, R);
    2390             : }
    2391             : 
    2392             : GEN
    2393       19418 : QXQ_norm(GEN A, GEN B)
    2394             : {
    2395             :   GEN c, R, lB;
    2396       19418 :   long dA = degpol(A), dB = degpol(B);
    2397       19418 :   pari_sp av = avma;
    2398       19418 :   if (dA < 0) return gen_0;
    2399       19418 :   A = Q_primitive_part(A, &c);
    2400       19418 :   R = ZX_resultant(B, A);
    2401       19418 :   if (c) R = gmul(R, gpowgs(c, dB));
    2402       19418 :   lB = leading_coeff(B);
    2403       19418 :   if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
    2404       19418 :   return gc_upto(av, R);
    2405             : }
    2406             : 
    2407             : /* assume x has integral coefficients */
    2408             : GEN
    2409     1200211 : ZX_disc_all(GEN x, ulong bound)
    2410             : {
    2411     1200211 :   pari_sp av = avma;
    2412     1200211 :   long s, d = degpol(x);
    2413             :   GEN l, R;
    2414             : 
    2415     1200209 :   if (d <= 1) return d == 1? gen_1: gen_0;
    2416     1196916 :   s = (d & 2) ? -1: 1;
    2417     1196916 :   l = leading_coeff(x);
    2418     1196915 :   if (!bound) bound = ZX_discbound(x);
    2419     1196801 :   R = ZX_resultant_all(x, NULL, NULL, bound);
    2420     1196885 :   if (is_pm1(l))
    2421     1017846 :   { if (signe(l) < 0) s = -s; }
    2422             :   else
    2423      179037 :     R = diviiexact(R,l);
    2424     1196883 :   if (s == -1) togglesign_safe(&R);
    2425     1196878 :   return gc_INT(av,R);
    2426             : }
    2427             : 
    2428             : GEN
    2429     1138151 : ZX_disc(GEN x) { return ZX_disc_all(x,0); }
    2430             : 
    2431             : static GEN
    2432       11010 : ZXQX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, GEN T, ulong p)
    2433             : {
    2434       11010 :   pari_sp av = avma;
    2435       11010 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2436             :   GEN H, dp;
    2437       11010 :   if (dropa && dropb) return pol0_Flx(T[1]); /* p | lc(A), p | lc(B) */
    2438       11010 :   H = FlxqX_saferesultant(a, b, T, p);
    2439       11008 :   if (!H) return NULL;
    2440       11008 :   if (dropa)
    2441             :   { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2442           0 :     GEN c = gel(b,degB+2); /* lc(B) */
    2443           0 :     if (odd(degB)) c = Flx_neg(c, p);
    2444           0 :     c = Flxq_powu(c, dropa, T, p);
    2445           0 :     if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
    2446             :   }
    2447       11008 :   else if (dropb)
    2448             :   { /* multiply by lc(A)^(deg B - deg b) */
    2449           0 :     GEN c = gel(a,degA+2); /* lc(A) */
    2450           0 :     c = Flxq_powu(c, dropb, T, p);
    2451           0 :     if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
    2452             :   }
    2453       11008 :   dp = dB ? ZX_to_Flx(dB, p): pol1_Flx(T[1]);
    2454       11009 :   if (!Flx_equal1(dp))
    2455             :   {
    2456           0 :     GEN idp = Flxq_invsafe(dp, T, p);
    2457           0 :     if (!idp) return NULL;
    2458           0 :     H = Flxq_mul(H, Flxq_powu(idp, degA, T, p), T, p);
    2459             :   }
    2460       11009 :   return gc_leaf(av, H);
    2461             : }
    2462             : 
    2463             : /* If B=NULL, assume B=A' */
    2464             : static GEN
    2465        4911 : ZXQX_resultant_slice(GEN A, GEN B, GEN U, GEN dB, GEN P, GEN *mod)
    2466             : {
    2467        4911 :   pari_sp av = avma;
    2468        4911 :   long degA, degB, i, n = lg(P)-1;
    2469             :   GEN H, T;
    2470        4911 :   long v = varn(U), redo = 0;
    2471             : 
    2472        4911 :   degA = degpol(A);
    2473        4911 :   degB = B? degpol(B): degA - 1;
    2474        4911 :   if (n == 1)
    2475             :   {
    2476        3177 :     ulong p = uel(P,1);
    2477        3177 :     GEN a = ZXX_to_FlxX(A, p, v), b = B? ZXX_to_FlxX(B, p, v): FlxX_deriv(a, p);
    2478        3177 :     GEN u = ZX_to_Flx(U, p);
    2479        3177 :     GEN Hp = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
    2480        3177 :     if (!Hp) { set_avma(av); *mod = gen_1; return pol_0(v); }
    2481        3177 :     Hp = gc_upto(av, Flx_to_ZX(Hp)); *mod = utoipos(p); return Hp;
    2482             :   }
    2483        1734 :   T = ZV_producttree(P);
    2484        1734 :   A = ZXX_nv_mod_tree(A, P, T, v);
    2485        1734 :   if (B) B = ZXX_nv_mod_tree(B, P, T, v);
    2486        1734 :   U = ZX_nv_mod_tree(U, P, T);
    2487        1734 :   H = cgetg(n+1, t_VEC);
    2488        9566 :   for(i=1; i <= n; i++)
    2489             :   {
    2490        7832 :     ulong p = P[i];
    2491        7832 :     GEN a = gel(A,i), b = B? gel(B,i): FlxX_deriv(a, p), u = gel(U, i);
    2492        7833 :     GEN h = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
    2493        7832 :     if (!h)
    2494             :     {
    2495           0 :       gel(H,i) = pol_0(v);
    2496           0 :       P[i] = 1; redo = 1;
    2497             :     }
    2498             :     else
    2499        7832 :       gel(H,i) = h;
    2500             :   }
    2501        1734 :   if (redo) T = ZV_producttree(P);
    2502        1734 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2503        1734 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2504             : }
    2505             : 
    2506             : GEN
    2507        4911 : ZXQX_resultant_worker(GEN P, GEN A, GEN B, GEN T, GEN dB)
    2508             : {
    2509        4911 :   GEN V = cgetg(3, t_VEC);
    2510        4911 :   if (isintzero(B)) B = NULL;
    2511        4911 :   if (!signe(dB)) dB = NULL;
    2512        4911 :   gel(V,1) = ZXQX_resultant_slice(A, B, T, dB, P, &gel(V,2));
    2513        4911 :   return V;
    2514             : }
    2515             : 
    2516             : static ulong
    2517        4315 : ZXQX_resultant_bound_i(GEN nf, GEN A, GEN B, GEN (*f)(GEN,GEN,long))
    2518             : {
    2519        4315 :   pari_sp av = avma;
    2520        4315 :   GEN r, M = nf_L2_bound(nf, NULL, &r);
    2521        4315 :   long v = nf_get_varn(nf), i, l = lg(r);
    2522        4315 :   GEN a = cgetg(l, t_COL);
    2523       13334 :   for (i = 1; i < l; i++)
    2524        9019 :     gel(a, i) = f(gsubst(A, v, gel(r,i)), gsubst(B, v, gel(r,i)), DEFAULTPREC);
    2525        4315 :   return gc_ulong(av, (ulong) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
    2526             : }
    2527             : static ulong
    2528        4000 : ZXQX_resultant_bound(GEN nf, GEN A, GEN B)
    2529        4000 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound); }
    2530             : 
    2531             : static GEN
    2532          56 : _ZXQ_powu(GEN x, ulong u, GEN T)
    2533          56 : { return typ(x) == t_INT? powiu(x, u): ZXQ_powu(x, u, T); }
    2534             : 
    2535             : /* Compute Res(A, B/dB) in Z[X]/T, assuming A,B in Z[X,Y], dB in Z or NULL (= 1)
    2536             :  * If B=NULL, take B = A' and assume deg A > 1 */
    2537             : static GEN
    2538        3997 : ZXQX_resultant_all(GEN A, GEN B, GEN T, GEN dB, ulong bound)
    2539             : {
    2540        3997 :   pari_sp av = avma;
    2541             :   forprime_t S;
    2542             :   GEN  H, worker;
    2543        3997 :   if (B)
    2544             :   {
    2545          63 :     long a = degpol(A), b = degpol(B);
    2546          63 :     if (a < 0 || b < 0) return gen_0;
    2547          63 :     if (!a) return _ZXQ_powu(gel(A,2), b, T);
    2548          63 :     if (!b) return _ZXQ_powu(gel(B,2), a, T);
    2549             :   } else
    2550        3934 :     if (!bound) B = RgX_deriv(A);
    2551        3997 :   if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, B);
    2552        3997 :   worker = snm_closure(is_entry("_ZXQX_resultant_worker"),
    2553             :                        mkvec4(A, B? B: gen_0, T, dB? dB: gen_0));
    2554        3997 :   init_modular_big(&S);
    2555        3997 :   H = gen_crt("ZXQX_resultant_all", worker, &S, dB, bound, 0, NULL,
    2556             :               nxV_chinese_center, FpX_center);
    2557        3997 :   if (DEBUGLEVEL)
    2558           0 :     err_printf("ZXQX_resultant_all: a priori bound: %lu, a posteriori: %lu\n",
    2559             :                bound, expi(gsupnorm(H, DEFAULTPREC)));
    2560        3997 :   return gc_upto(av, H);
    2561             : }
    2562             : 
    2563             : GEN
    2564         119 : nfX_resultant(GEN nf, GEN x, GEN y)
    2565             : {
    2566         119 :   pari_sp av = avma;
    2567         119 :   GEN cx, cy, D, T = nf_get_pol(nf);
    2568         119 :   long dx = degpol(x), dy = degpol(y);
    2569         119 :   if (dx < 0 || dy < 0) return gen_0;
    2570         119 :   x = Q_primitive_part(x, &cx); if (cx) cx = gpowgs(cx, dy);
    2571         119 :   y = Q_primitive_part(y, &cy); if (cy) cy = gpowgs(cy, dx);
    2572         119 :   if (!dx)      D = _ZXQ_powu(gel(x,2), dy, T);
    2573         119 :   else if (!dy) D = _ZXQ_powu(gel(y,2), dx, T);
    2574             :   else
    2575             :   {
    2576          63 :     ulong bound = ZXQX_resultant_bound(nf, x, y);
    2577          63 :     D = ZXQX_resultant_all(x, y, T, NULL, bound);
    2578             :   }
    2579         119 :   cx = mul_content(cx, cy); if (cx) D = gmul(D, cx);
    2580         119 :   return gc_upto(av, D);
    2581             : }
    2582             : 
    2583             : static GEN
    2584         252 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    2585             : 
    2586             : static GEN
    2587        3934 : ZXQX_disc_all(GEN x, GEN T, ulong bound)
    2588             : {
    2589        3934 :   pari_sp av = avma;
    2590        3934 :   long s, d = degpol(x), v = varn(T);
    2591             :   GEN l, R;
    2592             : 
    2593        3934 :   if (d <= 1) return d == 1? pol_1(v): pol_0(v);
    2594        3934 :   s = (d & 2) ? -1: 1;
    2595        3934 :   l = leading_coeff(x);
    2596        3934 :   R = ZXQX_resultant_all(x, NULL, T, NULL, bound);
    2597        3934 :   if (!gequal1(l)) R = QXQ_div(R, to_ZX(l,v), T);
    2598        3934 :   if (s == -1) R = RgX_neg(R);
    2599        3934 :   return gc_upto(av, R);
    2600             : }
    2601             : 
    2602             : GEN
    2603           7 : QX_disc(GEN x)
    2604             : {
    2605           7 :   pari_sp av = avma;
    2606           7 :   GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
    2607           7 :   if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
    2608           7 :   return gc_upto(av, d);
    2609             : }
    2610             : 
    2611             : GEN
    2612        4165 : nfX_disc(GEN nf, GEN x)
    2613             : {
    2614        4165 :   pari_sp av = avma;
    2615        4165 :   GEN c, D, T = nf_get_pol(nf);
    2616             :   ulong bound;
    2617        4165 :   long d = degpol(x), v = varn(T);
    2618        4165 :   if (d <= 1) return d == 1? pol_1(v): pol_0(v);
    2619        3934 :   x = Q_primitive_part(x, &c);
    2620        3934 :   bound = ZXQX_resultant_bound(nf, x, RgX_deriv(x));
    2621        3934 :   D = ZXQX_disc_all(x, T, bound);
    2622        3934 :   if (c) D = gmul(D, gpowgs(c, 2*d - 2));
    2623        3934 :   return gc_upto(av, D);
    2624             : }
    2625             : 
    2626             : GEN
    2627      846032 : QXQ_mul(GEN x, GEN y, GEN T)
    2628             : {
    2629      846032 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2630      846033 :   GEN dy, ny = Q_primitive_part(y, &dy);
    2631      846033 :   GEN z = ZXQ_mul(nx, ny, T);
    2632      846034 :   if (dx || dy)
    2633             :   {
    2634      843234 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
    2635      843234 :     if (!gequal1(d)) z = ZX_Q_mul(z, d);
    2636             :   }
    2637      846034 :   return z;
    2638             : }
    2639             : 
    2640             : GEN
    2641      407041 : QXQ_sqr(GEN x, GEN T)
    2642             : {
    2643      407041 :   GEN dx, nx = Q_primitive_part(x, &dx);
    2644      407041 :   GEN z = ZXQ_sqr(nx, T);
    2645      407041 :   if (dx)
    2646      405305 :     z = ZX_Q_mul(z, gsqr(dx));
    2647      407041 :   return z;
    2648             : }
    2649             : 
    2650             : static GEN
    2651      212909 : QXQ_inv_slice(GEN A, GEN B, GEN P, GEN *mod)
    2652             : {
    2653      212909 :   pari_sp av = avma;
    2654      212909 :   long i, n = lg(P)-1, v = varn(A), redo = 0;
    2655             :   GEN H, T;
    2656      212909 :   if (n == 1)
    2657             :   {
    2658      165821 :     ulong p = uel(P,1);
    2659      165821 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    2660      165821 :     GEN U = Flxq_invsafe(a, b, p);
    2661      165821 :     if (!U)
    2662             :     {
    2663          24 :       set_avma(av);
    2664          24 :       *mod = gen_1; return pol_0(v);
    2665             :     }
    2666      165797 :     H = gc_GEN(av, Flx_to_ZX(U));
    2667      165797 :     *mod = utoipos(p); return H;
    2668             :   }
    2669       47088 :   T = ZV_producttree(P);
    2670       47088 :   A = ZX_nv_mod_tree(A, P, T);
    2671       47088 :   B = ZX_nv_mod_tree(B, P, T);
    2672       47087 :   H = cgetg(n+1, t_VEC);
    2673      238084 :   for(i=1; i <= n; i++)
    2674             :   {
    2675      190996 :     ulong p = P[i];
    2676      190996 :     GEN a = gel(A,i), b = gel(B,i);
    2677      190996 :     GEN U = Flxq_invsafe(a, b, p);
    2678      190998 :     if (!U)
    2679             :     {
    2680         601 :       gel(H,i) = pol_0(v);
    2681         601 :       P[i] = 1; redo = 1;
    2682             :     }
    2683             :     else
    2684      190397 :       gel(H,i) = U;
    2685             :   }
    2686       47088 :   if (redo) T = ZV_producttree(P);
    2687       47088 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2688       47088 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2689             : }
    2690             : 
    2691             : GEN
    2692      212909 : QXQ_inv_worker(GEN P, GEN A, GEN B)
    2693             : {
    2694      212909 :   GEN V = cgetg(3, t_VEC);
    2695      212909 :   gel(V,1) = QXQ_inv_slice(A, B, P, &gel(V,2));
    2696      212909 :   return V;
    2697             : }
    2698             : 
    2699             : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
    2700             : GEN
    2701      146108 : QXQ_inv(GEN A, GEN B)
    2702             : {
    2703             :   GEN D, Ap, Bp;
    2704             :   ulong pp;
    2705      146108 :   pari_sp av2, av = avma;
    2706             :   forprime_t S;
    2707      146108 :   GEN worker, U, H = NULL, mod = gen_1;
    2708             :   pari_timer ti;
    2709             :   long k, dA, dB;
    2710      146108 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2711             :   /* A a QX, B a ZX */
    2712      146108 :   A = Q_primitive_part(A, &D);
    2713      146107 :   dA = degpol(A); dB= degpol(B);
    2714             :   /* A, B in Z[X] */
    2715      146107 :   init_modular_small(&S);
    2716             :   do {
    2717      146108 :     pp = u_forprime_next(&S);
    2718      146107 :     Ap = ZX_to_Flx(A, pp);
    2719      146108 :     Bp = ZX_to_Flx(B, pp);
    2720      146108 :   } while (degpol(Ap) != dA || degpol(Bp) != dB);
    2721      146108 :   if (degpol(Flx_gcd(Ap, Bp, pp)) != 0 && degpol(ZX_gcd(A,B))!=0)
    2722          14 :     pari_err_INV("QXQ_inv",mkpolmod(A,B));
    2723      146094 :   worker = snm_closure(is_entry("_QXQ_inv_worker"), mkvec2(A, B));
    2724      146094 :   av2 = avma;
    2725      146094 :   for (k = 1; ;k *= 2)
    2726       42591 :   {
    2727             :     GEN res, b, N, den;
    2728      188685 :     gen_inccrt_i("QXQ_inv", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
    2729             :                  nxV_chinese_center, FpX_center);
    2730      188685 :     (void)gc_all(av2, 2, &H, &mod);
    2731      188685 :     b = sqrti(shifti(mod,-1));
    2732      188685 :     if (DEBUGLEVEL>5) timer_start(&ti);
    2733      188685 :     U = FpX_ratlift(H, mod, b, b, NULL);
    2734      188685 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: ratlift");
    2735      194431 :     if (!U) continue;
    2736      151840 :     N = Q_remove_denom(U, &den); if (!den) den = gen_1;
    2737      151840 :     res = Flx_rem(Flx_Fl_sub(Flx_mul(Ap, ZX_to_Flx(N,pp), pp),
    2738             :                   umodiu(den, pp), pp), Bp, pp);
    2739      151840 :     if (degpol(res) >= 0) continue;
    2740      146094 :     res = ZX_Z_sub(ZX_mul(A, N), den);
    2741      146094 :     res = ZX_is_monic(B) ? ZX_rem(res, B): RgX_pseudorem(res, B);
    2742      146094 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: final check");
    2743      146094 :     if (degpol(res)<0)
    2744             :     {
    2745      146094 :       if (D) U = RgX_Rg_div(U, D);
    2746      146094 :       return gc_GEN(av, U);
    2747             :     }
    2748             :   }
    2749             : }
    2750             : 
    2751             : static GEN
    2752      121107 : QXQ_div_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
    2753             : {
    2754      121107 :   pari_sp av = avma;
    2755      121107 :   long i, n = lg(P)-1, v = varn(A), redo = 0;
    2756             :   GEN H, T;
    2757      121107 :   if (n == 1)
    2758             :   {
    2759       44701 :     ulong p = uel(P,1);
    2760       44701 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), c = ZX_to_Flx(C, p);
    2761       44700 :     GEN bi = Flxq_invsafe(b, c, p), U;
    2762       44700 :     if (!bi)
    2763             :     {
    2764           0 :       set_avma(av);
    2765           0 :       *mod = gen_1; return pol_0(v);
    2766             :     }
    2767       44700 :     U = Flxq_mul(a, bi, c, p);
    2768       44700 :     H = gc_GEN(av, Flx_to_ZX(U));
    2769       44701 :     *mod = utoipos(p); return H;
    2770             :   }
    2771       76406 :   T = ZV_producttree(P);
    2772       76406 :   A = ZX_nv_mod_tree(A, P, T);
    2773       76406 :   B = ZX_nv_mod_tree(B, P, T);
    2774       76406 :   C = ZX_nv_mod_tree(C, P, T);
    2775       76406 :   H = cgetg(n+1, t_VEC);
    2776      337444 :   for(i=1; i <= n; i++)
    2777             :   {
    2778      261038 :     ulong p = P[i];
    2779      261038 :     GEN a = gel(A,i), b = gel(B,i), c = gel(C, i);
    2780      261038 :     GEN bi = Flxq_invsafe(b, c, p);
    2781      261040 :     if (!bi)
    2782             :     {
    2783           4 :       gel(H,i) = pol_0(v);
    2784           4 :       P[i] = 1; redo = 1;
    2785             :     }
    2786             :     else
    2787      261036 :       gel(H,i) = Flxq_mul(a, bi, c, p);
    2788             :   }
    2789       76406 :   if (redo) T = ZV_producttree(P);
    2790       76406 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2791       76406 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2792             : }
    2793             : 
    2794             : GEN
    2795      121107 : QXQ_div_worker(GEN P, GEN A, GEN B, GEN C)
    2796             : {
    2797      121107 :   GEN V = cgetg(3, t_VEC);
    2798      121107 :   gel(V,1) = QXQ_div_slice(A, B, C, P, &gel(V,2));
    2799      121107 :   return V;
    2800             : }
    2801             : 
    2802             : /* lift(Mod(A/B, C)). C a ZX, A, B a scalar or a QX */
    2803             : GEN
    2804       33207 : QXQ_div(GEN A, GEN B, GEN C)
    2805             : {
    2806             :   GEN DA, DB, Ap, Bp, Cp;
    2807             :   ulong pp;
    2808       33207 :   pari_sp av2, av = avma;
    2809             :   forprime_t S;
    2810       33207 :   GEN worker, U, H = NULL, mod = gen_1;
    2811             :   pari_timer ti;
    2812             :   long k, dA, dB, dC;
    2813       33207 :   if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
    2814             :   /* A a QX, B a ZX */
    2815       33207 :   A = Q_primitive_part(A, &DA);
    2816       33206 :   B = Q_primitive_part(B, &DB);
    2817       33206 :   dA = degpol(A); dB = degpol(B); dC = degpol(C);
    2818             :   /* A, B in Z[X] */
    2819       33206 :   init_modular_small(&S);
    2820             :   do {
    2821       33207 :     pp = u_forprime_next(&S);
    2822       33207 :     Ap = ZX_to_Flx(A, pp);
    2823       33206 :     Bp = ZX_to_Flx(B, pp);
    2824       33206 :     Cp = ZX_to_Flx(C, pp);
    2825       33206 :   } while (degpol(Ap) != dA || degpol(Bp) != dB || degpol(Cp) != dC);
    2826       33206 :   if (degpol(Flx_gcd(Bp, Cp, pp)) != 0 && degpol(ZX_gcd(B,C))!=0)
    2827           0 :     pari_err_INV("QXQ_div",mkpolmod(B,C));
    2828       33206 :   worker = snm_closure(is_entry("_QXQ_div_worker"), mkvec3(A, B, C));
    2829       33207 :   av2 = avma;
    2830       33207 :   for (k = 1; ;k *= 2)
    2831       46720 :   {
    2832             :     GEN res, b, N, den;
    2833       79927 :     gen_inccrt_i("QXQ_div", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
    2834             :                  nxV_chinese_center, FpX_center);
    2835       79927 :     (void)gc_all(av2, 2, &H, &mod);
    2836       79927 :     b = sqrti(shifti(mod,-1));
    2837       79927 :     if (DEBUGLEVEL>5) timer_start(&ti);
    2838       79927 :     U = FpX_ratlift(H, mod, b, b, NULL);
    2839       79927 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: ratlift");
    2840       90556 :     if (!U) continue;
    2841       43836 :     N = Q_remove_denom(U, &den); if (!den) den = gen_1;
    2842       43836 :     res = Flx_rem(Flx_sub(Flx_mul(Bp, ZX_to_Flx(N,pp), pp),
    2843             :                           Flx_Fl_mul(Ap, umodiu(den, pp), pp), pp), Cp, pp);
    2844       43836 :     if (degpol(res) >= 0) continue;
    2845       33207 :     res = ZX_sub(ZX_mul(B, N), ZX_Z_mul(A,den));
    2846       33207 :     res = ZX_is_monic(C) ? ZX_rem(res, C): RgX_pseudorem(res, C);
    2847       33207 :     if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: final check");
    2848       33207 :     if (degpol(res)<0)
    2849             :     {
    2850       33207 :       if (DA && DB) U = RgX_Rg_mul(U, gdiv(DA,DB));
    2851       28069 :       else if (DA) U = RgX_Rg_mul(U, DA);
    2852       15981 :       else if (DB) U = RgX_Rg_div(U, DB);
    2853       33207 :       return gc_GEN(av, U);
    2854             :     }
    2855             :   }
    2856             : }
    2857             : 
    2858             : /************************************************************************
    2859             :  *                                                                      *
    2860             :  *                           ZXQ_minpoly                                *
    2861             :  *                                                                      *
    2862             :  ************************************************************************/
    2863             : 
    2864             : static GEN
    2865        3523 : ZXQ_minpoly_slice(GEN A, GEN B, long d, GEN P, GEN *mod)
    2866             : {
    2867        3523 :   pari_sp av = avma;
    2868        3523 :   long i, n = lg(P)-1, v = evalvarn(varn(B));
    2869             :   GEN H, T;
    2870        3523 :   if (n == 1)
    2871             :   {
    2872         716 :     ulong p = uel(P,1);
    2873         716 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    2874         716 :     GEN Hp = Flxq_minpoly(a, b, p);
    2875         716 :     if (degpol(Hp) != d) { p = 1; Hp = pol0_Flx(v); }
    2876         716 :     H = gc_upto(av, Flx_to_ZX(Hp));
    2877         716 :     *mod = utoipos(p); return H;
    2878             :   }
    2879        2807 :   T = ZV_producttree(P);
    2880        2807 :   A = ZX_nv_mod_tree(A, P, T);
    2881        2807 :   B = ZX_nv_mod_tree(B, P, T);
    2882        2807 :   H = cgetg(n+1, t_VEC);
    2883       16838 :   for(i=1; i <= n; i++)
    2884             :   {
    2885       14031 :     ulong p = P[i];
    2886       14031 :     GEN a = gel(A,i), b = gel(B,i);
    2887       14031 :     GEN m = Flxq_minpoly(a, b, p);
    2888       14031 :     if (degpol(m) != d) { P[i] = 1; m = pol0_Flx(v); }
    2889       14031 :     gel(H, i) = m;
    2890             :   }
    2891        2807 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2892        2807 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2893             : }
    2894             : 
    2895             : GEN
    2896        3523 : ZXQ_minpoly_worker(GEN P, GEN A, GEN B, long d)
    2897             : {
    2898        3523 :   GEN V = cgetg(3, t_VEC);
    2899        3523 :   gel(V,1) = ZXQ_minpoly_slice(A, B, d, P, &gel(V,2));
    2900        3523 :   return V;
    2901             : }
    2902             : 
    2903             : GEN
    2904        1701 : ZXQ_minpoly(GEN A, GEN B, long d, ulong bound)
    2905             : {
    2906        1701 :   pari_sp av = avma;
    2907             :   GEN worker, H, dB;
    2908             :   forprime_t S;
    2909        1701 :   B = Q_remove_denom(B, &dB);
    2910        1701 :   worker = strtoclosure("_ZXQ_minpoly_worker", 3, A, B, stoi(d));
    2911        1701 :   init_modular_big(&S);
    2912        1701 :   H = gen_crt("ZXQ_minpoly", worker, &S, dB, bound, 0, NULL,
    2913             :                nxV_chinese_center, FpX_center_i);
    2914        1701 :   return gc_GEN(av, H);
    2915             : }
    2916             : 
    2917             : /************************************************************************
    2918             :  *                                                                      *
    2919             :  *                   ZX_ZXY_resultant                                   *
    2920             :  *                                                                      *
    2921             :  ************************************************************************/
    2922             : 
    2923             : static GEN
    2924      364909 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
    2925             :                        long degA, long degB, long dres, long sX)
    2926             : {
    2927      364909 :   pari_sp av = avma;
    2928      364909 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    2929      364905 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2930      364905 :   GEN Hp = Flx_FlxY_resultant_polint(a, b, p, pi, dres, sX);
    2931      364908 :   if (dropa && dropb)
    2932           0 :     Hp = zero_Flx(sX);
    2933             :   else {
    2934      364908 :     if (dropa)
    2935             :     { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    2936           0 :       GEN c = gel(b,degB+2); /* lc(B) */
    2937           0 :       if (odd(degB)) c = Flx_neg(c, p);
    2938           0 :       if (!Flx_equal1(c)) {
    2939           0 :         c = Flx_powu_pre(c, dropa, p, pi);
    2940           0 :         if (!Flx_equal1(c)) Hp = Flx_mul_pre(Hp, c, p, pi);
    2941             :       }
    2942             :     }
    2943      364908 :     else if (dropb)
    2944             :     { /* multiply by lc(A)^(deg B - deg b) */
    2945           0 :       ulong c = uel(a, degA+2); /* lc(A) */
    2946           0 :       c = Fl_powu(c, dropb, p);
    2947           0 :       if (c != 1) Hp = Flx_Fl_mul_pre(Hp, c, p, pi);
    2948             :     }
    2949             :   }
    2950      364908 :   if (dp != 1) Hp = Flx_Fl_mul_pre(Hp, Fl_powu_pre(Fl_inv(dp,p), degA, p, pi), p, pi);
    2951      364907 :   return gc_leaf(av, Hp);
    2952             : }
    2953             : 
    2954             : static GEN
    2955      124963 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
    2956             :                        GEN P, GEN *mod, long sX, long vY)
    2957             : {
    2958      124963 :   pari_sp av = avma;
    2959      124963 :   long i, n = lg(P)-1;
    2960             :   GEN H, T, D;
    2961      124963 :   if (n == 1)
    2962             :   {
    2963       40164 :     ulong p = uel(P,1);
    2964       40164 :     ulong dp = dB ? umodiu(dB, p): 1;
    2965       40164 :     GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
    2966       40165 :     GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2967       40165 :     H = gc_upto(av, Flx_to_ZX(Hp));
    2968       40165 :     *mod = utoipos(p); return H;
    2969             :   }
    2970       84799 :   T = ZV_producttree(P);
    2971       84799 :   A = ZX_nv_mod_tree(A, P, T);
    2972       84799 :   B = ZXX_nv_mod_tree(B, P, T, vY);
    2973       84799 :   D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
    2974       84799 :   H = cgetg(n+1, t_VEC);
    2975      364209 :   for(i=1; i <= n; i++)
    2976             :   {
    2977      279411 :     ulong p = P[i];
    2978      279411 :     GEN a = gel(A,i), b = gel(B,i);
    2979      279411 :     ulong dp = D ? uel(D, i): 1;
    2980      279411 :     gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    2981             :   }
    2982       84798 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    2983       84799 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    2984             : }
    2985             : 
    2986             : GEN
    2987      124963 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
    2988             : {
    2989      124963 :   GEN V = cgetg(3, t_VEC);
    2990      124963 :   if (isintzero(dB)) dB = NULL;
    2991      124963 :   gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
    2992      124964 :   return V;
    2993             : }
    2994             : 
    2995             : GEN
    2996       79197 : ZX_ZXY_resultant(GEN A, GEN B)
    2997             : {
    2998       79197 :   pari_sp av = avma;
    2999             :   forprime_t S;
    3000             :   ulong bound;
    3001       79197 :   long v = fetch_var_higher();
    3002       79197 :   long degA = degpol(A), degB, dres = degA * degpol(B);
    3003       79197 :   long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
    3004       79197 :   long sX = evalvarn(vX);
    3005             :   GEN worker, H, dB;
    3006       79197 :   B = Q_remove_denom(B, &dB);
    3007       79197 :   if (!dB) B = leafcopy(B);
    3008       79197 :   A = leafcopy(A); setvarn(A,v);
    3009       79197 :   B = swap_vars(B, vY, v); degB = degpol(B);
    3010       79197 :   bound = ZX_ZXY_ResBound(A, B, dB);
    3011       79195 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    3012      158390 :   worker = snm_closure(is_entry("_ZX_ZXY_resultant_worker"),
    3013       79195 :                        mkvec4(A, B, dB? dB: gen_0,
    3014             :                               mkvecsmall5(degA, degB, dres, sX, vY)));
    3015       79197 :   init_modular_big(&S);
    3016       79197 :   H = gen_crt("ZX_ZXY_resultant_all", worker, &S, dB, bound, 0, NULL,
    3017             :                nxV_chinese_center, FpX_center_i);
    3018       79197 :   setvarn(H, vX); (void)delete_var();
    3019       79197 :   return gc_GEN(av, H);
    3020             : }
    3021             : 
    3022             : static GEN
    3023     1118026 : ZXX_resultant_prime(GEN a, GEN b, ulong p, long degA, long degB, long dres, long sX)
    3024             : {
    3025     1118026 :   pari_sp av = avma;
    3026     1118026 :   long dropa = degA - degpol(a), dropb = degB - degpol(b);
    3027     1118025 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    3028     1118025 :   GEN Hp = FlxX_resultant_polint(a, b, p, pi, dres, sX);
    3029     1118029 :   if (dropa && dropb)
    3030           0 :     Hp = zero_Flx(sX);
    3031             :   else {
    3032     1118029 :     if (dropa)
    3033             :     { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
    3034           0 :       GEN c = gel(b,degB+2); /* lc(B) */
    3035           0 :       if (odd(degB)) c = Flx_neg(c, p);
    3036           0 :       if (!Flx_equal1(c)) {
    3037           0 :         c = Flx_powu_pre(c, dropa, p, pi);
    3038           0 :         if (!Flx_equal1(c)) Hp = Flx_mul_pre(Hp, c, p, pi);
    3039             :       }
    3040             :     }
    3041     1118029 :     else if (dropb)
    3042             :     { /* multiply by lc(A)^(deg B - deg b) */
    3043           0 :       ulong c = uel(a, degA+2); /* lc(A) */
    3044           0 :       c = Fl_powu(c, dropb, p);
    3045           0 :       if (c != 1) Hp = Flx_Fl_mul_pre(Hp, c, p, pi);
    3046             :     }
    3047             :   }
    3048     1118029 :   return gc_leaf(av, Hp);
    3049             : }
    3050             : 
    3051             : static GEN
    3052     1015861 : ZXX_resultant_slice(GEN A, GEN B, long degA, long degB, long dres,
    3053             :                        GEN P, GEN *mod, long sX, long vY)
    3054             : {
    3055     1015861 :   pari_sp av = avma;
    3056     1015861 :   long i, n = lg(P)-1;
    3057             :   GEN H, T;
    3058     1015861 :   if (n == 1)
    3059             :   {
    3060      931872 :     ulong p = uel(P,1);
    3061      931872 :     GEN a = ZXX_to_FlxX(A, p, vY), b = ZXX_to_FlxX(B, p, vY);
    3062      931867 :     GEN Hp = ZXX_resultant_prime(a, b, p, degA, degB, dres, sX);
    3063      931871 :     H = gc_upto(av, Flx_to_ZX(Hp));
    3064      931872 :     *mod = utoipos(p); return H;
    3065             :   }
    3066       83989 :   T = ZV_producttree(P);
    3067       83989 :   A = ZXX_nv_mod_tree(A, P, T, vY);
    3068       83989 :   B = ZXX_nv_mod_tree(B, P, T, vY);
    3069       83989 :   H = cgetg(n+1, t_VEC);
    3070      270148 :   for(i=1; i <= n; i++)
    3071             :   {
    3072      186159 :     ulong p = P[i];
    3073      186159 :     GEN a = gel(A,i), b = gel(B,i);
    3074      186159 :     gel(H,i) = ZXX_resultant_prime(a, b, p, degA, degB, dres, sX);
    3075             :   }
    3076       83989 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    3077       83989 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    3078             : }
    3079             : 
    3080             : GEN
    3081     1015862 : ZXX_resultant_worker(GEN P, GEN A, GEN B, GEN v)
    3082             : {
    3083     1015862 :   GEN V = cgetg(3, t_VEC);
    3084     1015861 :   gel(V,1) = ZXX_resultant_slice(A, B, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
    3085     1015861 :   return V;
    3086             : }
    3087             : 
    3088             : static int
    3089     2021231 : ZXX_is_sparse(GEN x, long r)
    3090             : {
    3091     2021231 :   long c = 0, i, l = lg(x), ly = 2, thr;
    3092     8483448 :   for (i = 2; i < l; i++)
    3093             :   {
    3094     6462217 :     GEN xi = gel(x,i);
    3095     6462217 :     if (typ(xi)==t_INT)
    3096     5450543 :       c += !!signe(xi);
    3097             :     else
    3098             :     {
    3099     1011674 :       long j, li = lg(xi);
    3100     1011674 :       ly = maxss(ly, li);
    3101     3072668 :       for (j = 2; j < li; j++)
    3102     2060994 :         c += !!signe(gel(xi,j));
    3103             :     }
    3104             :   }
    3105     2021231 :   thr = (l-2)*(ly-2)/r;
    3106     2021231 :   if (DEBUGLEVEL >= 5)
    3107           0 :     err_printf("ZXX_is_sparse: lx=%ld ly=%ld r=%ld c=%ld < thr=%ld \n",l-2, ly-2,r,c,thr);
    3108     2021232 :   return c < thr;
    3109             : }
    3110             : 
    3111             : 
    3112             : static GEN
    3113     1010589 : ZXX_resultant_interp(GEN A, GEN B, long vX)
    3114             : {
    3115     1010589 :   pari_sp av = avma;
    3116             :   forprime_t S;
    3117             :   ulong bound;
    3118     1010589 :   long degA = degpol(A), degB = degpol(B), dres = degA * RgXY_degreex(B) + RgXY_degreex(A)*degB;
    3119     1010585 :   long vY = varn(A), sX = evalvarn(vX);
    3120             :   GEN worker, H;
    3121     1010585 :   bound = ZXX_ResBound(A, B);
    3122     1010589 :   if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
    3123     1010589 :   worker = snm_closure(is_entry("_ZXX_resultant_worker"),
    3124             :                        mkvec3(A, B, mkvecsmall5(degA, degB, dres, sX, vY)));
    3125     1010592 :   init_modular_big(&S);
    3126     1010591 :   H = gen_crt("ZXX_resultant", worker, &S, NULL, bound, 0, NULL, nxV_chinese_center, FpX_center_i);
    3127     1010591 :   return gc_GEN(av, H);
    3128             : }
    3129             : 
    3130             : GEN
    3131     1010624 : ZXX_resultant(GEN A, GEN B, long vX)
    3132             : {
    3133     1010624 :   if (!ZXX_is_sparse(A,6) && !ZXX_is_sparse(B,6))
    3134     1010589 :     return ZXX_resultant_interp(A, B, vX);
    3135          35 :   return RgX_resultant_all(A, B, NULL);
    3136             : }
    3137             : 
    3138             : static long
    3139       40523 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
    3140             : {
    3141       40523 :   pari_sp av = avma;
    3142       40523 :   long degA = degpol(A), degB, dres = degA*degpol(B0);
    3143       40523 :   long v = fetch_var_higher();
    3144       40523 :   long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
    3145       40523 :   long sX = evalvarn(vX);
    3146             :   GEN dB, B, a, b, Hp;
    3147             :   forprime_t S;
    3148             : 
    3149       40523 :   B0 = Q_remove_denom(B0, &dB);
    3150       40523 :   if (!dB) B0 = leafcopy(B0);
    3151       40523 :   A = leafcopy(A);
    3152       40523 :   B = B0;
    3153       40523 :   setvarn(A,v);
    3154       45334 : INIT:
    3155       45334 :   if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
    3156       45334 :   B = swap_vars(B, vY, v);
    3157             :   /* B0(lambda v + x, v) */
    3158       45333 :   if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    3159             : 
    3160       45333 :   degB = degpol(B);
    3161       45333 :   init_modular_big(&S);
    3162             :   while (1)
    3163           0 :   {
    3164       45333 :     ulong p = u_forprime_next(&S);
    3165       45333 :     ulong dp = dB ? umodiu(dB, p): 1;
    3166       45333 :     if (!dp) continue;
    3167       45333 :     a = ZX_to_Flx(A, p);
    3168       45334 :     b = ZXX_to_FlxX(B, p, v);
    3169       45334 :     Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
    3170       45334 :     if (degpol(Hp) != dres) continue;
    3171       45334 :     if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
    3172       45334 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
    3173       40523 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    3174       40523 :     (void)delete_var(); return gc_long(av,lambda);
    3175             :   }
    3176             : }
    3177             : 
    3178             : GEN
    3179       60563 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
    3180             : {
    3181       60563 :   if (lambda)
    3182             :   {
    3183       40523 :     *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
    3184       40523 :     if (*lambda) B = RgX_Rg_translate(B, monomial(stoi(*lambda), 1, varn(A)));
    3185             :   }
    3186       60563 :   return ZX_ZXY_resultant(A,B);
    3187             : }
    3188             : 
    3189             : static GEN
    3190       10350 : ZX_composedsum_slice(GEN A, GEN B, GEN P, GEN *mod)
    3191             : {
    3192       10350 :   pari_sp av = avma;
    3193       10350 :   long i, n = lg(P)-1;
    3194             :   GEN H, T;
    3195       10350 :   if (n == 1)
    3196             :   {
    3197        9848 :     ulong p = uel(P,1);
    3198        9848 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    3199        9846 :     GEN Hp = Flx_composedsum(a, b, p);
    3200        9842 :     H = gc_upto(av, Flx_to_ZX(Hp));
    3201        9850 :     *mod = utoipos(p); return H;
    3202             :   }
    3203         502 :   T = ZV_producttree(P);
    3204         502 :   A = ZX_nv_mod_tree(A, P, T);
    3205         502 :   B = ZX_nv_mod_tree(B, P, T);
    3206         502 :   H = cgetg(n+1, t_VEC);
    3207        4526 :   for(i=1; i <= n; i++)
    3208             :   {
    3209        4024 :     ulong p = P[i];
    3210        4024 :     GEN a = gel(A,i), b = gel(B,i);
    3211        4024 :     gel(H,i) = Flx_composedsum(a, b, p);
    3212             :   }
    3213         502 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    3214         502 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    3215             : }
    3216             : 
    3217             : GEN
    3218       10349 : ZX_composedsum_worker(GEN P, GEN A, GEN B)
    3219             : {
    3220       10349 :   GEN V = cgetg(3, t_VEC);
    3221       10350 :   gel(V,1) = ZX_composedsum_slice(A, B, P, &gel(V,2));
    3222       10351 :   return V;
    3223             : }
    3224             : 
    3225             : static GEN
    3226       10085 : ZX_composedsum_i(GEN A, GEN B, GEN lead)
    3227             : {
    3228       10085 :   pari_sp av = avma;
    3229             :   forprime_t S;
    3230             :   ulong bound;
    3231             :   GEN H, worker, mod;
    3232       10085 :   if (degpol(A) < degpol(B)) swap(A, B);
    3233       10084 :   if (!lead) lead  = mulii(leading_coeff(A),leading_coeff(B));
    3234       10084 :   bound = ZX_ZXY_ResBound_1(A, B);
    3235       10085 :   worker = snm_closure(is_entry("_ZX_composedsum_worker"), mkvec2(A,B));
    3236       10086 :   init_modular_big(&S);
    3237       10084 :   H = gen_crt("ZX_composedsum", worker, &S, lead, bound, 0, &mod,
    3238             :               nxV_chinese_center, FpX_center);
    3239       10087 :   return gc_upto(av, H);
    3240             : }
    3241             : 
    3242             : static long
    3243        9696 : ZX_compositum_lambda(GEN A, GEN B, GEN lead, long lambda)
    3244             : {
    3245        9696 :   pari_sp av = avma;
    3246             :   forprime_t S;
    3247             :   ulong p;
    3248        9696 :   init_modular_big(&S);
    3249        9699 :   p = u_forprime_next(&S);
    3250             :   while (1)
    3251         112 :   {
    3252             :     GEN Hp, a;
    3253        9811 :     if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    3254        9811 :     if (lead && dvdiu(lead,p)) { p = u_forprime_next(&S); continue; }
    3255        9804 :     a = ZX_to_Flx(ZX_rescale(A, stoi(-lambda)), p);
    3256        9804 :     Hp = Flx_composedsum(a, ZX_to_Flx(B, p), p);
    3257        9802 :     if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); continue; }
    3258        9693 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    3259        9693 :     return gc_long(av, lambda);
    3260             :   }
    3261             : }
    3262             : 
    3263             : GEN
    3264        9699 : ZX_compositum(GEN A, GEN B, long *lambda)
    3265             : {
    3266        9699 :   GEN lead  = mulii(leading_coeff(A),leading_coeff(B));
    3267        9696 :   if (lambda)
    3268             :   {
    3269        9696 :     *lambda = ZX_compositum_lambda(A, B, lead, *lambda);
    3270        9693 :     A = ZX_rescale(A, stoi(-*lambda));
    3271             :   }
    3272        9700 :   return ZX_composedsum_i(A, B, lead);
    3273             : }
    3274             : 
    3275             : GEN
    3276         385 : ZX_composedsum(GEN A, GEN B)
    3277         385 : { return ZX_composedsum_i(A, B, NULL); }
    3278             : 
    3279             : static GEN
    3280         359 : ZXQX_composedsum_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
    3281             : {
    3282         359 :   pari_sp av = avma;
    3283         359 :   long i, n = lg(P)-1, dC = degpol(C), v = varn(C);
    3284             :   GEN H, T;
    3285         359 :   if (n == 1)
    3286             :   {
    3287         181 :     ulong p = uel(P,1);
    3288         181 :     GEN a = ZXX_to_FlxX(A, p, v), b = ZXX_to_FlxX(B, p, v);
    3289         181 :     GEN c = ZX_to_Flx(C, p);
    3290         181 :     GEN Hp = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
    3291         181 :     H = gc_upto(av, Flm_to_ZM(Hp));
    3292         181 :     *mod = utoipos(p); return H;
    3293             :   }
    3294         178 :   T = ZV_producttree(P);
    3295         178 :   A = ZXX_nv_mod_tree(A, P, T, v);
    3296         178 :   B = ZXX_nv_mod_tree(B, P, T, v);
    3297         178 :   C = ZX_nv_mod_tree(C, P, T);
    3298         178 :   H = cgetg(n+1, t_VEC);
    3299         660 :   for(i=1; i <= n; i++)
    3300             :   {
    3301         482 :     ulong p = P[i];
    3302         482 :     GEN a = gel(A,i), b = gel(B,i), c = gel(C,i);
    3303         482 :     gel(H,i) = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
    3304             :   }
    3305         178 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
    3306         178 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    3307             : }
    3308             : 
    3309             : GEN
    3310         359 : ZXQX_composedsum_worker(GEN P, GEN A, GEN B, GEN C)
    3311             : {
    3312         359 :   GEN V = cgetg(3, t_VEC);
    3313         359 :   gel(V,1) = ZXQX_composedsum_slice(A, B, C, P, &gel(V,2));
    3314         359 :   return V;
    3315             : }
    3316             : 
    3317             : static GEN
    3318         315 : ZXQX_composedsum(GEN A, GEN B, GEN T, ulong bound)
    3319             : {
    3320         315 :   pari_sp av = avma;
    3321             :   forprime_t S;
    3322             :   GEN H, worker, mod;
    3323         315 :   GEN lead = mulii(Q_content(leading_coeff(A)), Q_content(leading_coeff(B)));
    3324         315 :   worker = snm_closure(is_entry("_ZXQX_composedsum_worker")
    3325             :                       , mkvec3(A,B,T));
    3326         315 :   init_modular_big(&S);
    3327         315 :   H = gen_crt("ZXQX_composedsum", worker, &S, lead, bound, 0, &mod,
    3328             :               nmV_chinese_center, FpM_center);
    3329         315 :   if (DEBUGLEVEL > 4)
    3330           0 :     err_printf("nfcompositum: a priori bound: %lu, a posteriori: %lu\n",
    3331             :                bound, expi(gsupnorm(H, DEFAULTPREC)));
    3332         315 :   return gc_GEN(av, RgM_to_RgXX(H, varn(A), varn(T)));
    3333             : }
    3334             : 
    3335             : static long
    3336         315 : ZXQX_composedsum_bound(GEN nf, GEN A, GEN B)
    3337         315 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound_1); }
    3338             : 
    3339             : GEN
    3340         315 : nf_direct_compositum(GEN nf, GEN A, GEN B)
    3341             : {
    3342         315 :   ulong bnd = ZXQX_composedsum_bound(nf, A, B);
    3343         315 :   return ZXQX_composedsum(A, B, nf_get_pol(nf), bnd);
    3344             : }
    3345             : 
    3346             : /************************************************************************
    3347             :  *                                                                      *
    3348             :  *                   IRREDUCIBLE POLYNOMIAL / Fp                        *
    3349             :  *                                                                      *
    3350             :  ************************************************************************/
    3351             : 
    3352             : /* irreducible (unitary) polynomial of degree n over Fp */
    3353             : GEN
    3354           0 : ffinit_rand(GEN p,long n)
    3355             : {
    3356           0 :   for(;;) {
    3357           0 :     pari_sp av = avma;
    3358           0 :     GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
    3359           0 :     if (FpX_is_irred(pol, p)) return pol;
    3360           0 :     set_avma(av);
    3361             :   }
    3362             : }
    3363             : 
    3364             : /* return an extension of degree 2^l of F_2, assume l > 0
    3365             :  * Not stack clean. */
    3366             : static GEN
    3367         600 : ffinit_Artin_Schreier_2(long l)
    3368             : {
    3369             :   GEN Q, T, S;
    3370             :   long i, v;
    3371             : 
    3372         600 :   if (l == 1) return mkvecsmall4(0,1,1,1); /*x^2 + x + 1*/
    3373         551 :   v = fetch_var_higher();
    3374         551 :   S = mkvecsmall5(0, 0, 0, 1, 1); /* y(y^2 + y) */
    3375         551 :   Q = mkpoln(3, pol1_Flx(0), pol1_Flx(0), S); /* x^2 + x + y(y^2+y) */
    3376         551 :   setvarn(Q, v);
    3377             : 
    3378             :   /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
    3379         551 :   T = mkvecsmalln(6,evalvarn(v),1UL,1UL,0UL,0UL,1UL);
    3380             :   /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
    3381             :    * ==> x^2 + x + a(y) b irred. over K for any root b of Q
    3382             :    * ==> x^2 + x + (b^2+b)b */
    3383        3040 :   for (i=2; i<l; i++) T = Flx_FlxY_resultant(T, Q, 2); /* minpoly(b) / F2*/
    3384         551 :   (void)delete_var(); T[1] = 0; return T;
    3385             : }
    3386             : 
    3387             : /* return an extension of degree p^l of F_p, assume l > 0
    3388             :  * Not stack clean. */
    3389             : GEN
    3390         957 : ffinit_Artin_Schreier(ulong p, long l)
    3391             : {
    3392             :   long i, v;
    3393             :   GEN Q, R, S, T, xp;
    3394         957 :   if (p==2) return ffinit_Artin_Schreier_2(l);
    3395         357 :   xp = polxn_Flx(p,0); /* x^p */
    3396         357 :   T = Flx_sub(xp, mkvecsmall3(0,1,1),p); /* x^p - x - 1 */
    3397         357 :   if (l == 1) return T;
    3398             : 
    3399           7 :   v = evalvarn(fetch_var_higher());
    3400           7 :   xp[1] = v;
    3401           7 :   R = Flx_sub(polxn_Flx(2*p-1,0), polxn_Flx(p,0),p);
    3402           7 :   S = Flx_sub(xp, polx_Flx(0), p);
    3403           7 :   Q = FlxY_Flx_sub(Flx_to_FlxX(S, v), R, p); /* x^p - x - (y^(2p-1)-y^p) */
    3404          14 :   for (i = 2; i <= l; ++i) T = Flx_FlxY_resultant(T, Q, p);
    3405           7 :   (void)delete_var(); T[1] = 0; return T;
    3406             : }
    3407             : 
    3408             : static long
    3409      149711 : flinit_check(ulong p, long n, long l)
    3410             : {
    3411             :   ulong q;
    3412      149711 :   if (!uisprime(n)) return 0;
    3413      102496 :   q = p % n; if (!q) return 0;
    3414       99885 :   return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    3415             : }
    3416             : 
    3417             : static GEN
    3418       31944 : flinit(ulong p, long l)
    3419             : {
    3420       31944 :   ulong n = 1+l;
    3421       96782 :   while (!flinit_check(p,n,l)) n += l;
    3422       31944 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    3423       31944 :   return ZX_to_Flx(polsubcyclo(n,l,0), p);
    3424             : }
    3425             : 
    3426             : static GEN
    3427       28991 : ffinit_fact_Flx(ulong p, long n)
    3428             : {
    3429       28991 :   GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
    3430       28991 :   long i, l = lg(Fm);
    3431       28991 :   P = cgetg(l, t_VEC);
    3432       61892 :   for (i = 1; i < l; i++)
    3433       32901 :     gel(P,i) = p==uel(Fp,i) ? ffinit_Artin_Schreier(p, Fe[i])
    3434       32901 :                             : flinit(p, uel(Fm,i));
    3435       28991 :   return FlxV_composedsum(P, p);
    3436             : }
    3437             : 
    3438             : static GEN
    3439       52936 : init_Flxq_i(ulong p, long n, long sv)
    3440             : {
    3441             :   GEN P;
    3442       52936 :   if (!odd(p) && p != 2) pari_err_PRIME("ffinit", utoi(p));
    3443       52929 :   if (n == 1) return polx_Flx(sv);
    3444       52929 :   if (flinit_check(p, n+1, n))
    3445             :   {
    3446       23938 :     P = const_vecsmall(n+2,1);
    3447       23938 :     P[1] = sv; return P;
    3448             :   }
    3449       28991 :   P = ffinit_fact_Flx(p,n);
    3450       28991 :   P[1] = sv; return P;
    3451             : }
    3452             : 
    3453             : GEN
    3454           0 : init_Flxq(ulong p, long n, long v)
    3455             : {
    3456           0 :   pari_sp av = avma;
    3457           0 :   return gc_upto(av, init_Flxq_i(p, n, v));
    3458             : }
    3459             : 
    3460             : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
    3461             : static long
    3462        8207 : fpinit_check(GEN p, long n, long l)
    3463             : {
    3464             :   ulong q;
    3465        8207 :   if (!uisprime(n)) return 0;
    3466        4842 :   q = umodiu(p,n); if (!q) return 0;
    3467        4842 :   return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
    3468             : }
    3469             : 
    3470             : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
    3471             :  * Return an irreducible polynomial of degree l over F_p.
    3472             :  * Variant of Adleman and Lenstra "Finding irreducible polynomials over
    3473             :  * finite fields", ACM, 1986 (5) 350--355.
    3474             :  * Not stack clean */
    3475             : static GEN
    3476        1828 : fpinit(GEN p, long l)
    3477             : {
    3478        1828 :   ulong n = 1+l;
    3479        6168 :   while (!fpinit_check(p,n,l)) n += l;
    3480        1828 :   if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
    3481        1828 :   return FpX_red(polsubcyclo(n,l,0),p);
    3482             : }
    3483             : 
    3484             : static GEN
    3485        1637 : ffinit_fact(GEN p, long n)
    3486             : {
    3487        1637 :   GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
    3488        1637 :   long i, l = lg(Fm);
    3489        1637 :   P = cgetg(l, t_VEC);
    3490        3465 :   for (i = 1; i < l; ++i)
    3491        3656 :     gel(P,i) = absequaliu(p, Fp[i]) ?
    3492           0 :                  Flx_to_ZX(ffinit_Artin_Schreier(Fp[i], Fe[i]))
    3493        1828 :                : fpinit(p, Fm[i]);
    3494        1637 :   return FpXV_composedsum(P, p);
    3495             : }
    3496             : 
    3497             : static GEN
    3498       55241 : init_Fq_i(GEN p, long n, long v)
    3499             : {
    3500             :   GEN P;
    3501       55241 :   if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
    3502       55241 :   if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
    3503       55241 :   if (cmpiu(p, 2) < 0) pari_err_PRIME("ffinit",p);
    3504       55234 :   if (v < 0) v = 0;
    3505       55234 :   if (n == 1) return pol_x(v);
    3506       54982 :   if (lgefint(p) == 3)
    3507       52936 :     return Flx_to_ZX(init_Flxq_i(p[2], n, evalvarn(v)));
    3508        2046 :   if (!mpodd(p)) pari_err_PRIME("ffinit", p);
    3509        2039 :   if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
    3510        1637 :   P = ffinit_fact(p,n);
    3511        1637 :   setvarn(P, v); return P;
    3512             : }
    3513             : GEN
    3514       54674 : init_Fq(GEN p, long n, long v)
    3515             : {
    3516       54674 :   pari_sp av = avma;
    3517       54674 :   return gc_upto(av, init_Fq_i(p, n, v));
    3518             : }
    3519             : GEN
    3520         567 : ffinit(GEN p, long n, long v)
    3521             : {
    3522         567 :   pari_sp av = avma;
    3523         567 :   return gc_upto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
    3524             : }
    3525             : 
    3526             : GEN
    3527        3178 : ffnbirred(GEN p, long n)
    3528             : {
    3529        3178 :   pari_sp av = avma;
    3530        3178 :   GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
    3531        3178 :   long j, l = lg(D);
    3532        6797 :   for (j = 2; j < l; j++) /* skip d = 1 */
    3533             :   {
    3534        3619 :     long md = D[j]; /* mu(d) * d, d squarefree */
    3535        3619 :     GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
    3536        3619 :     s = md > 0? addii(s, pd): subii(s,pd);
    3537             :   }
    3538        3178 :   return gc_INT(av, diviuexact(s, n));
    3539             : }
    3540             : 
    3541             : GEN
    3542         616 : ffsumnbirred(GEN p, long n)
    3543             : {
    3544         616 :   pari_sp av = avma, av2;
    3545         616 :   GEN q, t = p, v = vecfactoru_i(1, n);
    3546             :   long i;
    3547         616 :   q = cgetg(n+1,t_VEC); gel(q,1) = p;
    3548        1764 :   for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
    3549         616 :   av2 = avma;
    3550        1764 :   for (i=2; i<=n; i++)
    3551             :   {
    3552        1148 :     GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
    3553        1148 :     long j, l = lg(D);
    3554        2534 :     for (j = 2; j < l; j++) /* skip 1 */
    3555             :     {
    3556        1386 :       long md = D[j];
    3557        1386 :       GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
    3558        1386 :       s = md > 0? addii(s, pd): subii(s, pd);
    3559             :     }
    3560        1148 :     t = gc_INT(av2, addii(t, diviuexact(s, i)));
    3561             :   }
    3562         616 :   return gc_INT(av, t);
    3563             : }
    3564             : 
    3565             : GEN
    3566         140 : ffnbirred0(GEN p, long n, long flag)
    3567             : {
    3568         140 :   if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
    3569         140 :   if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
    3570         140 :   switch(flag)
    3571             :   {
    3572          70 :     case 0: return ffnbirred(p, n);
    3573          70 :     case 1: return ffsumnbirred(p, n);
    3574             :   }
    3575           0 :   pari_err_FLAG("ffnbirred");
    3576             :   return NULL; /* LCOV_EXCL_LINE */
    3577             : }
    3578             : 
    3579             : static void
    3580        2261 : checkmap(GEN m, const char *s)
    3581             : {
    3582        2261 :   if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
    3583           0 :     pari_err_TYPE(s,m);
    3584        2261 : }
    3585             : 
    3586             : GEN
    3587         189 : ffembed(GEN a, GEN b)
    3588             : {
    3589         189 :   pari_sp av = avma;
    3590         189 :   GEN p, Ta, Tb, g, r = NULL;
    3591         189 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
    3592         189 :   if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
    3593         189 :   p = FF_p_i(a); g = FF_gen(a);
    3594         189 :   if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
    3595         189 :   Ta = FF_mod(a);
    3596         189 :   Tb = FF_mod(b);
    3597         189 :   if (degpol(Tb)%degpol(Ta)!=0)
    3598           7 :     pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
    3599         182 :   r = gel(FFX_roots(Ta, b), 1);
    3600         182 :   return gc_GEN(av, mkvec2(g,r));
    3601             : }
    3602             : 
    3603             : GEN
    3604          91 : ffextend(GEN a, GEN P, long v)
    3605             : {
    3606          91 :   pari_sp av = avma;
    3607             :   long n;
    3608             :   GEN p, T, R, g, m;
    3609          91 :   if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
    3610          91 :   T = a; p = FF_p_i(a);
    3611          91 :   if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
    3612          49 :   if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
    3613          49 :   if (v < 0) v = varn(P);
    3614          49 :   n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
    3615          49 :   m = ffembed(a, g);
    3616          49 :   R = FFX_roots(ffmap(m, P),g);
    3617          49 :   return gc_GEN(av, mkvec2(gel(R,1), m));
    3618             : }
    3619             : 
    3620             : GEN
    3621          42 : fffrobenius(GEN a, long n)
    3622             : {
    3623          42 :   if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
    3624          42 :   retmkvec2(FF_gen(a), FF_Frobenius(a, n));
    3625             : }
    3626             : 
    3627             : GEN
    3628         133 : ffinvmap(GEN m)
    3629             : {
    3630         133 :   pari_sp av = avma;
    3631             :   long i, l;
    3632         133 :   GEN T, F, a, g, r, f = NULL;
    3633         133 :   checkmap(m, "ffinvmap");
    3634         133 :   a = gel(m,1); r = gel(m,2);
    3635         133 :   if (typ(r) != t_FFELT)
    3636           7 :    pari_err_TYPE("ffinvmap", m);
    3637         126 :   g = FF_gen(a);
    3638         126 :   T = FF_mod(r);
    3639         126 :   F = gel(FFX_factor(T, a), 1);
    3640         126 :   l = lg(F);
    3641         490 :   for(i=1; i<l; i++)
    3642             :   {
    3643         490 :     GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
    3644         490 :     if (degpol(s)==0 && gequal(constant_coeff(s),g)) { f = gel(F, i); break; }
    3645             :   }
    3646         126 :   if (f==NULL) pari_err_TYPE("ffinvmap", m);
    3647         126 :   if (degpol(f)==1) f = FF_neg_i(gel(f,2));
    3648         126 :   return gc_GEN(av, mkvec2(FF_gen(r),f));
    3649             : }
    3650             : 
    3651             : static GEN
    3652        1260 : ffpartmapimage(const char *s, GEN r)
    3653             : {
    3654        1260 :    GEN a = NULL, p = NULL;
    3655        1260 :    if (typ(r)==t_POL && degpol(r) >= 1
    3656        1260 :       && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
    3657           0 :    pari_err_TYPE(s, r);
    3658             :    return NULL; /* LCOV_EXCL_LINE */
    3659             : }
    3660             : 
    3661             : static GEN
    3662        2709 : ffeltmap_i(GEN m, GEN x)
    3663             : {
    3664        2709 :    GEN r = gel(m,2);
    3665        2709 :    if (!FF_samefield(x, gel(m,1)))
    3666          84 :      pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
    3667        2625 :    if (typ(r)==t_FFELT)
    3668        1659 :      return FF_map(r, x);
    3669             :    else
    3670         966 :      return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
    3671             : }
    3672             : 
    3673             : static GEN
    3674        4459 : ffmap_i(GEN m, GEN x)
    3675             : {
    3676             :   GEN y;
    3677        4459 :   long i, lx, tx = typ(x);
    3678        4459 :   switch(tx)
    3679             :   {
    3680        2541 :     case t_FFELT:
    3681        2541 :       return ffeltmap_i(m, x);
    3682        1267 :     case t_POL: case t_RFRAC: case t_SER:
    3683             :     case t_VEC: case t_COL: case t_MAT:
    3684        1267 :       y = cgetg_copy(x, &lx);
    3685        1988 :       for (i = 1; i < lontyp[tx]; i++) y[i] = x[i];
    3686        4564 :       for (; i < lx; i++)
    3687             :       {
    3688        3339 :         GEN yi = ffmap_i(m, gel(x,i));
    3689        3297 :         if (!yi) return NULL;
    3690        3297 :         gel(y,i) = yi;
    3691             :       }
    3692        1225 :       return y;
    3693             :   }
    3694         651 :   return gcopy(x);
    3695             : }
    3696             : 
    3697             : GEN
    3698        1036 : ffmap(GEN m, GEN x)
    3699             : {
    3700        1036 :   pari_sp ltop = avma;
    3701             :   GEN y;
    3702        1036 :   checkmap(m, "ffmap");
    3703        1036 :   y = ffmap_i(m, x);
    3704        1036 :   if (y) return y;
    3705          42 :   retgc_const(ltop, cgetg(1, t_VEC));
    3706             : }
    3707             : 
    3708             : static GEN
    3709         252 : ffeltmaprel_i(GEN m, GEN x)
    3710             : {
    3711         252 :    GEN g = gel(m,1), r = gel(m,2);
    3712         252 :    if (!FF_samefield(x, g))
    3713           0 :      pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
    3714         252 :    if (typ(r)==t_FFELT)
    3715          84 :      retmkpolmod(FF_map(r, x), pol_x(FF_var(g)));
    3716             :    else
    3717         168 :      retmkpolmod(FFX_preimagerel(x, r, ffpartmapimage("ffmap", r)), gcopy(r));
    3718             : }
    3719             : 
    3720             : static GEN
    3721         252 : ffmaprel_i(GEN m, GEN x)
    3722             : {
    3723         252 :   switch(typ(x))
    3724             :   {
    3725         252 :     case t_FFELT:
    3726         252 :       return ffeltmaprel_i(m, x);
    3727           0 :     case t_POL: pari_APPLY_pol_normalized(ffmaprel_i(m, gel(x,i)));
    3728           0 :     case t_SER: pari_APPLY_ser_normalized(ffmaprel_i(m, gel(x,i)));
    3729           0 :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
    3730           0 :       pari_APPLY_same(ffmaprel_i(m, gel(x,i)));
    3731             :   }
    3732           0 :   return gcopy(x);
    3733             : }
    3734             : GEN
    3735         252 : ffmaprel(GEN m, GEN x) { checkmap(m, "ffmaprel"); return ffmaprel_i(m, x); }
    3736             : 
    3737             : static void
    3738          84 : err_compo(GEN m, GEN n)
    3739          84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
    3740             : 
    3741             : GEN
    3742         420 : ffcompomap(GEN m, GEN n)
    3743             : {
    3744         420 :   pari_sp av = avma;
    3745         420 :   GEN g = gel(n,1), r, m2, n2;
    3746         420 :   checkmap(m, "ffcompomap");
    3747         420 :   checkmap(n, "ffcompomap");
    3748         420 :   m2 = gel(m,2); n2 = gel(n,2);
    3749         420 :   switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
    3750             :   {
    3751          84 :     case 0:
    3752          84 :       if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
    3753          42 :       r = FF_map(gel(m,2), n2);
    3754          42 :       break;
    3755          84 :     case 2:
    3756          84 :       r = ffmap_i(m, n2);
    3757          42 :       if (lg(r) == 1) err_compo(m,n);
    3758          42 :       break;
    3759         168 :     case 1:
    3760         168 :       r = ffeltmap_i(m, n2);
    3761         126 :       if (!r)
    3762             :       {
    3763             :         GEN a, A, R, M;
    3764             :         long dm, dn;
    3765          42 :         a = ffpartmapimage("ffcompomap",m2);
    3766          42 :         A = FF_to_FpXQ_i(FF_neg(n2));
    3767          42 :         setvarn(A, 1);
    3768          42 :         R = deg1pol(gen_1, A, 0);
    3769          42 :         setvarn(R, 0);
    3770          42 :         M = gcopy(m2);
    3771          42 :         setvarn(M, 1);
    3772          42 :         r = polresultant0(R, M, 1, 0);
    3773          42 :         dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
    3774          42 :         if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
    3775          42 :         setvarn(r, varn(FF_mod(g)));
    3776             :       }
    3777         126 :       break;
    3778          84 :     case 3:
    3779             :     {
    3780             :       GEN M, R, T, p, a;
    3781          84 :       a = ffpartmapimage("ffcompomap",n2);
    3782          84 :       if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
    3783          42 :       p = FF_p_i(gel(n,1));
    3784          42 :       T = FF_mod(gel(n,1));
    3785          42 :       setvarn(T, 1);
    3786          42 :       R = RgX_to_FpXQX(n2,T,p);
    3787          42 :       setvarn(R, 0);
    3788          42 :       M = gcopy(m2);
    3789          42 :       setvarn(M, 1);
    3790          42 :       r = polresultant0(R, M, 1, 0);
    3791          42 :       setvarn(r, varn(n2));
    3792             :     }
    3793             :   }
    3794         252 :   return gc_GEN(av, mkvec2(g,r));
    3795             : }

Generated by: LCOV version 1.16