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 - RgX.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29454-e7b6e21fb9) Lines: 1620 1773 91.4 %
Date: 2024-07-18 09:03:47 Functions: 200 218 91.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; 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             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_pol
      19             : 
      20             : /*******************************************************************/
      21             : /*                                                                 */
      22             : /*                         GENERIC                                 */
      23             : /*                                                                 */
      24             : /*******************************************************************/
      25             : 
      26             : /* Return optimal parameter l for the evaluation of n/m polynomials of degree d
      27             :    Fractional values can be used if the evaluations are done with different
      28             :    accuracies, and thus have different weights.
      29             :  */
      30             : long
      31    19639733 : brent_kung_optpow(long d, long n, long m)
      32             : {
      33             :   long p, r;
      34    19639733 :   long pold=1, rold=n*(d-1);
      35    92456669 :   for(p=2; p<=d; p++)
      36             :   {
      37    72816936 :     r = m*(p-1) + n*((d-1)/p);
      38    72816936 :     if (r<rold) { pold=p; rold=r; }
      39             :   }
      40    19639733 :   return pold;
      41             : }
      42             : 
      43             : static GEN
      44    14446027 : gen_RgXQ_eval_powers(GEN P, GEN V, long a, long n, void *E, const struct bb_algebra *ff,
      45             :                                            GEN cmul(void *E, GEN P, long a, GEN x))
      46             : {
      47    14446027 :   pari_sp av = avma;
      48             :   long i;
      49    14446027 :   GEN z = cmul(E,P,a,ff->one(E));
      50    14421002 :   if (!z) z = gen_0;
      51    75475151 :   for (i=1; i<=n; i++)
      52             :   {
      53    61054786 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      54    61078171 :     if (t) {
      55    46149697 :       z = ff->add(E, z, t);
      56    46115066 :       if (gc_needed(av,2)) z = gerepileupto(av, z);
      57             :     }
      58             :   }
      59    14420365 :   return ff->red(E,z);
      60             : }
      61             : 
      62             : /* Brent & Kung
      63             :  * (Fast algorithms for manipulating formal power series, JACM 25:581-595, 1978)
      64             :  *
      65             :  * V as output by FpXQ_powers(x,l,T,p). For optimal performance, l is as given
      66             :  * by brent_kung_optpow */
      67             : GEN
      68    12396296 : gen_bkeval_powers(GEN P, long d, GEN V, void *E, const struct bb_algebra *ff,
      69             :                                      GEN cmul(void *E, GEN P, long a, GEN x))
      70             : {
      71    12396296 :   pari_sp av = avma;
      72    12396296 :   long l = lg(V)-1;
      73             :   GEN z, u;
      74             : 
      75    12396296 :   if (d < 0) return ff->zero(E);
      76    11790183 :   if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      77     1147539 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      78     1147539 :   if (DEBUGLEVEL>=8)
      79             :   {
      80           0 :     long cnt = 1 + (d - l) / (l-1);
      81           0 :     err_printf("RgX_RgXQV_eval(%ld/%ld): %ld RgXQ_mul\n", d, l-1, cnt);
      82             :   }
      83     1147539 :   d -= l;
      84     1147539 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      85     2654543 :   while (d >= l-1)
      86             :   {
      87     1505691 :     d -= l-1;
      88     1505691 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      89     1505756 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      90     1505660 :     if (gc_needed(av,2))
      91         105 :       z = gerepileupto(av, z);
      92             :   }
      93     1148852 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      94     1148852 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      95     1148846 :   return gerepileupto(av, ff->red(E,z));
      96             : }
      97             : 
      98             : GEN
      99      858752 : gen_bkeval(GEN Q, long d, GEN x, int use_sqr, void *E, const struct bb_algebra *ff,
     100             :                                       GEN cmul(void *E, GEN P, long a, GEN x))
     101             : {
     102      858752 :   pari_sp av = avma;
     103             :   GEN z, V;
     104             :   long rtd;
     105      858752 :   if (d < 0) return ff->zero(E);
     106      857807 :   rtd = (long) sqrt((double)d);
     107      857807 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     108      857831 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     109      857824 :   return gerepileupto(av, z);
     110             : }
     111             : 
     112             : static GEN
     113     2027287 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     114             : static GEN
     115    19230726 : _gen_add(void *E, GEN x, GEN y) { (void)E; return gadd(x, y); }
     116             : static GEN
     117           0 : _gen_sub(void *E, GEN x, GEN y) { (void)E; return gsub(x, y); }
     118             : static GEN
     119     1842326 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     120             : static GEN
     121      599488 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     122             : static GEN
     123     2068347 : _gen_one(void *E) { (void)E; return gen_1; }
     124             : static GEN
     125       24184 : _gen_zero(void *E) { (void)E; return gen_0; }
     126             : 
     127             : static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
     128             :               _gen_mul, _gen_sqr,_gen_one,_gen_zero };
     129             : 
     130             : static GEN
     131      512818 : _gen_cmul(void *E, GEN P, long a, GEN x)
     132      512818 : {(void)E; return gmul(gel(P,a+2), x);}
     133             : 
     134             : GEN
     135      166768 : RgX_RgV_eval(GEN Q, GEN x)
     136             : {
     137      166768 :   return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
     138             : }
     139             : 
     140             : GEN
     141           0 : RgX_Rg_eval_bk(GEN Q, GEN x)
     142             : {
     143           0 :   return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
     144             : }
     145             : 
     146             : GEN
     147        2947 : RgXV_RgV_eval(GEN Q, GEN x)
     148             : {
     149        2947 :   long i, l = lg(Q), vQ = gvar(Q);
     150        2947 :   GEN v = cgetg(l, t_VEC);
     151      248311 :   for (i = 1; i < l; i++)
     152             :   {
     153      245364 :     GEN Qi = gel(Q, i);
     154      245364 :     gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
     155             :   }
     156        2947 :   return v;
     157             : }
     158             : 
     159             : GEN
     160      312514 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
     161             : {
     162      312514 :   pari_sp av = avma, btop;
     163      312514 :   long i, d = degpol(P), o;
     164             :   GEN s;
     165      312514 :   if (signe(P)==0) return pol_0(varn(P));
     166      312514 :   s = gel(P, d+2);
     167      312514 :   if (d == 0) return gcopy(s);
     168      309238 :   o = RgX_deflate_order(P);
     169      309238 :   if (o > 1) A = gpowgs(A, o);
     170      309238 :   btop = avma;
     171     1303740 :   for (i = d-o; i >= 0; i-=o)
     172             :   {
     173      994504 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
     174      994502 :     if (gc_needed(btop,1))
     175             :     {
     176          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
     177          13 :       s = gerepileupto(btop, s);
     178             :     }
     179             :   }
     180      309236 :   return gerepileupto(av, s);
     181             : }
     182             : 
     183             : GEN
     184        1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
     185             : {
     186        1652 :   pari_sp av = avma;
     187        1652 :   long i, d = degpol(P), v = varn(A);
     188             :   GEN s;
     189        1652 :   if (signe(P)==0) return pol_0(v);
     190        1652 :   if (d == 0) return scalarpol(gel(P, d+2), v);
     191        1232 :   s = scalarpol_shallow(gel(P, d+2), v);
     192        4963 :   for (i = d-1; i >= 0; i--)
     193             :   {
     194        3731 :     GEN c = gel(P,i+2), b = gel(B,d+1-i);
     195        3731 :     s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
     196        3731 :     if (gc_needed(av,1))
     197             :     {
     198           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
     199           0 :       s = gerepileupto(av, s);
     200             :     }
     201             :   }
     202        1232 :   return gerepileupto(av, s);
     203             : }
     204             : 
     205             : const struct bb_algebra *
     206      281024 : get_Rg_algebra(void)
     207             : {
     208      281024 :   return &Rg_algebra;
     209             : }
     210             : 
     211             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     212             : 
     213             : static GEN
     214       30352 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     215             : {
     216             :   (void) E;
     217       30352 :   return RgX_divrem(x, y, r);
     218             : }
     219             : 
     220             : GEN
     221        6923 : RgX_digits(GEN x, GEN T)
     222             : {
     223        6923 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     224        6923 :   return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     225             : }
     226             : 
     227             : /*******************************************************************/
     228             : /*                                                                 */
     229             : /*                         RgX                                     */
     230             : /*                                                                 */
     231             : /*******************************************************************/
     232             : 
     233             : long
     234    24774028 : RgX_equal(GEN x, GEN y)
     235             : {
     236    24774028 :   long i = lg(x);
     237             : 
     238    24774028 :   if (i != lg(y)) return 0;
     239   106501685 :   for (i--; i > 1; i--)
     240    82032599 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     241    24469086 :   return 1;
     242             : }
     243             : 
     244             : /* Returns 1 in the base ring over which x is defined */
     245             : /* HACK: this also works for t_SER */
     246             : GEN
     247   156663862 : Rg_get_1(GEN x)
     248             : {
     249             :   GEN p, T;
     250   156663862 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     251   156663862 :   if (RgX_type_is_composite(tx))
     252    11623330 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     253   156663862 :   switch(tx)
     254             :   {
     255      795550 :     case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
     256         973 :     case t_PADIC: return cvtop(gen_1, p, lx);
     257        5481 :     case t_FFELT: return FF_1(T);
     258   155861858 :     default: return gen_1;
     259             :   }
     260             : }
     261             : /* Returns 0 in the base ring over which x is defined */
     262             : /* HACK: this also works for t_SER */
     263             : GEN
     264     6974331 : Rg_get_0(GEN x)
     265             : {
     266             :   GEN p, T;
     267     6974331 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     268     6974331 :   if (RgX_type_is_composite(tx))
     269       51170 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     270     6974331 :   switch(tx)
     271             :   {
     272         532 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     273          42 :     case t_PADIC: return zeropadic(p, lx);
     274         231 :     case t_FFELT: return FF_zero(T);
     275     6973526 :     default: return gen_0;
     276             :   }
     277             : }
     278             : 
     279             : GEN
     280        5642 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     281             : {
     282        5642 :   long i, n = degpol(P);
     283             :   GEN z, dz, dP;
     284        5642 :   if (n < 0) return gen_0;
     285        5642 :   P = Q_remove_denom(P, &dP);
     286        5642 :   z = gel(P,2); if (n == 0) return icopy(z);
     287        3353 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     288        3353 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     289        6930 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     290        3353 :   dz = mul_denom(dP, dV);
     291        3353 :   return dz? RgX_Rg_div(z, dz): z;
     292             : }
     293             : 
     294             : /* Return P(h * x), not memory clean */
     295             : GEN
     296       19672 : RgX_unscale(GEN P, GEN h)
     297             : {
     298       19672 :   long i, l = lg(P);
     299       19672 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     300       19672 :   Q[1] = P[1];
     301       19672 :   if (l == 2) return Q;
     302       19644 :   gel(Q,2) = gcopy(gel(P,2));
     303       42181 :   for (i=3; i<l; i++)
     304             :   {
     305       22537 :     hi = gmul(hi,h);
     306       22537 :     gel(Q,i) = gmul(gel(P,i), hi);
     307             :   }
     308       19644 :   return Q;
     309             : }
     310             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     311             : GEN
     312     1451102 : ZX_z_unscale(GEN P, long h)
     313             : {
     314     1451102 :   long i, l = lg(P);
     315     1451102 :   GEN Q = cgetg(l, t_POL);
     316     1451102 :   Q[1] = P[1];
     317     1451102 :   if (l == 2) return Q;
     318     1362598 :   gel(Q,2) = gel(P,2);
     319     1362598 :   if (l == 3) return Q;
     320     1336299 :   if (h == -1)
     321      230304 :     for (i = 3; i < l; i++)
     322             :     {
     323      191145 :       gel(Q,i) = negi(gel(P,i));
     324      191144 :       if (++i == l) break;
     325      141509 :       gel(Q,i) = gel(P,i);
     326             :     }
     327             :   else
     328             :   {
     329             :     GEN hi;
     330     1247504 :     gel(Q,3) = mulis(gel(P,3), h);
     331     1247503 :     hi = sqrs(h);
     332     5361626 :     for (i = 4; i < l; i++)
     333             :     {
     334     4114121 :       gel(Q,i) = mulii(gel(P,i), hi);
     335     4114119 :       if (i != l-1) hi = mulis(hi,h);
     336             :     }
     337             :   }
     338     1336299 :   return Q;
     339             : }
     340             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     341             : GEN
     342      999771 : ZX_unscale(GEN P, GEN h)
     343             : {
     344             :   long i, l;
     345             :   GEN Q, hi;
     346      999771 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     347         804 :   l = lg(P); Q = cgetg(l, t_POL);
     348         804 :   Q[1] = P[1];
     349         804 :   if (l == 2) return Q;
     350         804 :   gel(Q,2) = gel(P,2);
     351         804 :   if (l == 3) return Q;
     352         804 :   hi = h;
     353         804 :   gel(Q,3) = mulii(gel(P,3), hi);
     354        2531 :   for (i = 4; i < l; i++)
     355             :   {
     356        1727 :     hi = mulii(hi,h);
     357        1727 :     gel(Q,i) = mulii(gel(P,i), hi);
     358             :   }
     359         804 :   return Q;
     360             : }
     361             : /* P a ZX. Return P(x << n), not memory clean */
     362             : GEN
     363      216832 : ZX_unscale2n(GEN P, long n)
     364             : {
     365      216832 :   long i, ni = n, l = lg(P);
     366      216832 :   GEN Q = cgetg(l, t_POL);
     367      216834 :   Q[1] = P[1];
     368      216834 :   if (l == 2) return Q;
     369      216834 :   gel(Q,2) = gel(P,2);
     370      216834 :   if (l == 3) return Q;
     371      216834 :   gel(Q,3) = shifti(gel(P,3), ni);
     372      801251 :   for (i=4; i<l; i++)
     373             :   {
     374      584433 :     ni += n;
     375      584433 :     gel(Q,i) = shifti(gel(P,i), ni);
     376             :   }
     377      216818 :   return Q;
     378             : }
     379             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     380             : GEN
     381       12073 : ZX_unscale_div(GEN P, GEN h)
     382             : {
     383       12073 :   long i, l = lg(P);
     384       12073 :   GEN hi, Q = cgetg(l, t_POL);
     385       12073 :   Q[1] = P[1];
     386       12073 :   if (l == 2) return Q;
     387       12073 :   gel(Q,2) = diviiexact(gel(P,2), h);
     388       12073 :   if (l == 3) return Q;
     389       12073 :   gel(Q,3) = gel(P,3);
     390       12073 :   if (l == 4) return Q;
     391       12073 :   hi = h;
     392       12073 :   gel(Q,4) = mulii(gel(P,4), hi);
     393       61349 :   for (i=5; i<l; i++)
     394             :   {
     395       49276 :     hi = mulii(hi,h);
     396       49276 :     gel(Q,i) = mulii(gel(P,i), hi);
     397             :   }
     398       12073 :   return Q;
     399             : }
     400             : /* P(h*X) / h^k, assuming the result is a ZX */
     401             : GEN
     402         966 : ZX_unscale_divpow(GEN P, GEN h, long k)
     403             : {
     404         966 :   long i, j, l = lg(P);
     405         966 :   GEN H, Q = cgetg(l, t_POL);
     406         966 :   Q[1] = P[1]; if (l == 2) return Q;
     407         966 :   H = gpowers(h, maxss(k, l - 3 - k));
     408        3864 :   for (i = 2, j = k+1; j > 1 && i < l; i++)
     409        2898 :     gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
     410         966 :   if (i == l) return Q;
     411         966 :   gel(Q, i) = gel(P, i); i++;
     412        3542 :   for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
     413         966 :   return Q;
     414             : }
     415             : 
     416             : GEN
     417        5761 : RgXV_unscale(GEN x, GEN h)
     418             : {
     419        5761 :   if (isint1(h)) return RgX_copy(x);
     420       16465 :   pari_APPLY_same(RgX_unscale(gel(x,i), h));
     421             : }
     422             : 
     423             : /* Return h^degpol(P) P(x / h), not memory clean */
     424             : GEN
     425     4320414 : RgX_rescale(GEN P, GEN h)
     426             : {
     427     4320414 :   long i, l = lg(P);
     428     4320414 :   GEN Q = cgetg(l,t_POL), hi = h;
     429     4320411 :   gel(Q,l-1) = gel(P,l-1);
     430    11360608 :   for (i=l-2; i>=2; i--)
     431             :   {
     432    11357989 :     gel(Q,i) = gmul(gel(P,i), hi);
     433    11357857 :     if (i == 2) break;
     434     7039899 :     hi = gmul(hi,h);
     435             :   }
     436     4320577 :   Q[1] = P[1]; return Q;
     437             : }
     438             : 
     439             : GEN
     440        2401 : RgXV_rescale(GEN x, GEN h)
     441             : {
     442        2401 :   if (isint1(h)) return RgX_copy(x);
     443       16086 :   pari_APPLY_same(RgX_rescale(gel(x,i), h));
     444             : }
     445             : 
     446             : /* A(X^d) --> A(X) */
     447             : GEN
     448      976191 : RgX_deflate(GEN x0, long d)
     449             : {
     450             :   GEN z, y, x;
     451      976191 :   long i,id, dy, dx = degpol(x0);
     452      976190 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     453      371285 :   dy = dx/d;
     454      371285 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     455      371280 :   z = y + 2;
     456      371280 :   x = x0+ 2;
     457     1495595 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     458      371280 :   return y;
     459             : }
     460             : 
     461             : GEN
     462       18732 : RgX_homogenize(GEN P, long v)
     463             : {
     464             :   long i, l, d;
     465       18732 :   GEN Q = cgetg_copy(P, &l);
     466       18732 :   Q[1] = P[1]; d = l-3;
     467      165998 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
     468       18732 :   return Q;
     469             : }
     470             : 
     471             : /* F a t_RFRAC */
     472             : long
     473          98 : rfrac_deflate_order(GEN F)
     474             : {
     475          98 :   GEN N = gel(F,1), D = gel(F,2);
     476          98 :   long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
     477          98 :   if (m == 1) return 1;
     478          35 :   if (typ(N) == t_POL && varn(N) == varn(D))
     479          21 :     m = cgcd(m, RgX_deflate_order(N));
     480          35 :   return m;
     481             : }
     482             : /* F a t_RFRAC */
     483             : GEN
     484          98 : rfrac_deflate_max(GEN F, long *m)
     485             : {
     486          98 :   *m = rfrac_deflate_order(F);
     487          98 :   return rfrac_deflate(F, *m);
     488             : }
     489             : /* F a t_RFRAC */
     490             : GEN
     491          98 : rfrac_deflate(GEN F, long m)
     492             : {
     493          98 :   GEN N = gel(F,1), D = gel(F,2);
     494          98 :   if (m == 1) return F;
     495          35 :   if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
     496          35 :   D = RgX_deflate(D, m); return mkrfrac(N, D);
     497             : }
     498             : 
     499             : /* return x0(X^d) */
     500             : GEN
     501      797909 : RgX_inflate(GEN x0, long d)
     502             : {
     503      797909 :   long i, id, dy, dx = degpol(x0);
     504      797909 :   GEN x = x0 + 2, z, y;
     505      797909 :   if (dx <= 0) return leafcopy(x0);
     506      790307 :   dy = dx*d;
     507      790307 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     508      790307 :   z = y + 2;
     509    27845658 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     510    11405125 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     511      790307 :   return y;
     512             : }
     513             : 
     514             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     515             : static GEN
     516     5945935 : RgX_translate_basecase(GEN P, GEN c)
     517             : {
     518     5945935 :   pari_sp av = avma;
     519             :   GEN Q, R;
     520             :   long i, k, n;
     521             : 
     522     5945935 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     523     5943985 :   Q = leafcopy(P);
     524     5943984 :   R = Q+2; n = degpol(P);
     525     5943982 :   if (isint1(c))
     526             :   {
     527       26094 :     for (i=1; i<=n; i++)
     528             :     {
     529      134643 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
     530       22259 :       if (gc_needed(av,2))
     531             :       {
     532           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
     533           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     534             :       }
     535             :     }
     536             :   }
     537     5940147 :   else if (isintm1(c))
     538             :   {
     539      899506 :     for (i=1; i<=n; i++)
     540             :     {
     541     2841608 :       for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
     542      682111 :       if (gc_needed(av,2))
     543             :       {
     544           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
     545           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     546             :       }
     547             :     }
     548             :   }
     549             :   else
     550             :   {
     551    19553810 :     for (i=1; i<=n; i++)
     552             :     {
     553    45875254 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
     554    13831061 :       if (gc_needed(av,2))
     555             :       {
     556           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
     557           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     558             :       }
     559             :     }
     560             :   }
     561     5943661 :   return gerepilecopy(av, Q);
     562             : }
     563             : GEN
     564     5946359 : RgX_translate(GEN P, GEN c)
     565             : {
     566     5946359 :   pari_sp av = avma;
     567     5946359 :   long n = degpol(P);
     568     5946356 :   if (n < 40)
     569     5945936 :     return RgX_translate_basecase(P, c);
     570             :   else
     571             :   {
     572         420 :     long d = n >> 1;
     573         420 :     GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
     574         420 :     GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
     575         420 :     GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
     576         420 :     return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
     577             :   }
     578             : }
     579             : 
     580             : /* P(ax + b) */
     581             : GEN
     582           0 : RgX_affine(GEN P, GEN a, GEN b)
     583             : {
     584           0 :   if (signe(b)) P = RgX_translate(P, b);
     585           0 :   return RgX_unscale(P, a);
     586             : }
     587             : 
     588             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     589             : GEN
     590       23630 : RgXQX_translate(GEN P, GEN c, GEN T)
     591             : {
     592       23630 :   pari_sp av = avma;
     593             :   GEN Q, R;
     594             :   long i, k, n;
     595             : 
     596       23630 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     597       23322 :   Q = leafcopy(P);
     598       23322 :   R = Q+2; n = degpol(P);
     599       82416 :   for (i=1; i<=n; i++)
     600             :   {
     601      261092 :     for (k=n-i; k<n; k++)
     602             :     {
     603      201998 :       pari_sp av2 = avma;
     604      201998 :       gel(R,k) = gerepileupto(av2,
     605      201998 :                    RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
     606             :     }
     607       59094 :     if (gc_needed(av,2))
     608             :     {
     609           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     610           0 :       Q = gerepilecopy(av, Q); R = Q+2;
     611             :     }
     612             :   }
     613       23322 :   return gerepilecopy(av, Q);
     614             : }
     615             : 
     616             : /********************************************************************/
     617             : /**                                                                **/
     618             : /**                          CONVERSIONS                           **/
     619             : /**                       (not memory clean)                       **/
     620             : /**                                                                **/
     621             : /********************************************************************/
     622             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     623             :  * but everything else is */
     624             : static GEN
     625      166691 : QXQ_to_mod(GEN x, GEN T)
     626             : {
     627             :   long d;
     628      166691 :   switch(typ(x))
     629             :   {
     630       62332 :     case t_INT:  return icopy(x);
     631        2079 :     case t_FRAC: return gcopy(x);
     632      102280 :     case t_POL:
     633      102280 :       d = degpol(x);
     634      102280 :       if (d < 0) return gen_0;
     635      100012 :       if (d == 0) return gcopy(gel(x,2));
     636       97873 :       return mkpolmod(RgX_copy(x), T);
     637           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     638             :              return NULL;/* LCOV_EXCL_LINE */
     639             :   }
     640             : }
     641             : /* pure shallow version */
     642             : GEN
     643      709909 : QXQ_to_mod_shallow(GEN x, GEN T)
     644             : {
     645             :   long d;
     646      709909 :   switch(typ(x))
     647             :   {
     648      454715 :     case t_INT:
     649      454715 :     case t_FRAC: return x;
     650      255194 :     case t_POL:
     651      255194 :       d = degpol(x);
     652      255194 :       if (d < 0) return gen_0;
     653      209791 :       if (d == 0) return gel(x,2);
     654      194796 :       return mkpolmod(x, T);
     655           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     656             :              return NULL;/* LCOV_EXCL_LINE */
     657             :   }
     658             : }
     659             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
     660             :  * Not memory clean because T not copied, but everything else is */
     661             : static GEN
     662       37198 : QXQX_to_mod(GEN z, GEN T)
     663             : {
     664       37198 :   long i,l = lg(z);
     665       37198 :   GEN x = cgetg(l,t_POL);
     666      185864 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     667       37198 :   x[1] = z[1]; return normalizepol_lg(x,l);
     668             : }
     669             : /* pure shallow version */
     670             : GEN
     671      161781 : QXQX_to_mod_shallow(GEN z, GEN T)
     672             : {
     673      161781 :   long i,l = lg(z);
     674      161781 :   GEN x = cgetg(l,t_POL);
     675      790910 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
     676      161781 :   x[1] = z[1]; return normalizepol_lg(x,l);
     677             : }
     678             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     679             : GEN
     680       12502 : QXQXV_to_mod(GEN V, GEN T)
     681             : {
     682       12502 :   long i, l = lg(V);
     683       12502 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     684       49700 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     685       12502 :   return z;
     686             : }
     687             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     688             : GEN
     689       17840 : QXQV_to_mod(GEN V, GEN T)
     690             : {
     691       17840 :   long i, l = lg(V);
     692       17840 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     693       35865 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
     694       17840 :   return z;
     695             : }
     696             : 
     697             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     698             : GEN
     699       13412 : QXQC_to_mod_shallow(GEN V, GEN T)
     700             : {
     701       13412 :   long i, l = lg(V);
     702       13412 :   GEN z = cgetg(l, t_COL);
     703       94192 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
     704       13412 :   return z;
     705             : }
     706             : 
     707             : GEN
     708        6251 : QXQM_to_mod_shallow(GEN V, GEN T)
     709             : {
     710        6251 :   long i, l = lg(V);
     711        6251 :   GEN z = cgetg(l, t_MAT);
     712       19663 :   for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
     713        6251 :   return z;
     714             : }
     715             : 
     716             : GEN
     717     6178687 : RgX_renormalize_lg(GEN x, long lx)
     718             : {
     719             :   long i;
     720     9199547 :   for (i = lx-1; i>1; i--)
     721     8749883 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     722     6178687 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     723     6178687 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     724             : }
     725             : 
     726             : GEN
     727     1494752 : RgV_to_RgX(GEN x, long v)
     728             : {
     729     1494752 :   long i, k = lg(x);
     730             :   GEN p;
     731             : 
     732     4161595 :   while (--k && gequal0(gel(x,k)));
     733     1494751 :   if (!k) return pol_0(v);
     734     1487907 :   i = k+2; p = cgetg(i,t_POL);
     735     1487905 :   p[1] = evalsigne(1) | evalvarn(v);
     736    10246547 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     737     1487905 :   return p;
     738             : }
     739             : GEN
     740      189522 : RgV_to_RgX_reverse(GEN x, long v)
     741             : {
     742      189522 :   long j, k, l = lg(x);
     743             :   GEN p;
     744             : 
     745      191069 :   for (k = 1; k < l; k++)
     746      191069 :     if (!gequal0(gel(x,k))) break;
     747      189522 :   if (k == l) return pol_0(v);
     748      189522 :   k -= 1;
     749      189522 :   l -= k;
     750      189522 :   x += k;
     751      189522 :   p = cgetg(l+1,t_POL);
     752      189522 :   p[1] = evalsigne(1) | evalvarn(v);
     753      997162 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     754      189522 :   return p;
     755             : }
     756             : 
     757             : /* return the (N-dimensional) vector of coeffs of p */
     758             : GEN
     759    13147309 : RgX_to_RgC(GEN x, long N)
     760             : {
     761             :   long i, l;
     762             :   GEN z;
     763    13147309 :   l = lg(x)-1; x++;
     764    13147309 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     765    13147309 :   z = cgetg(N+1,t_COL);
     766    79956807 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     767    24051676 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     768    13147400 :   return z;
     769             : }
     770             : GEN
     771     1345228 : Rg_to_RgC(GEN x, long N)
     772             : {
     773     1345228 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     774             : }
     775             : 
     776             : /* vector of polynomials (in v) whose coefs are given by the columns of x */
     777             : GEN
     778      301265 : RgM_to_RgXV(GEN x, long v)
     779     1286747 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
     780             : GEN
     781        7097 : RgM_to_RgXV_reverse(GEN x, long v)
     782       28388 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
     783             : 
     784             : /* matrix whose entries are given by the coeffs of the polynomials in
     785             :  * vector v (considered as degree n-1 polynomials) */
     786             : GEN
     787      333594 : RgV_to_RgM(GEN x, long n)
     788     1676328 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
     789             : 
     790             : GEN
     791       77567 : RgXV_to_RgM(GEN x, long n)
     792      392972 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
     793             : 
     794             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     795             : GEN
     796       23466 : RgM_to_RgXX(GEN x, long v,long w)
     797             : {
     798       23466 :   long j, lx = lg(x);
     799       23466 :   GEN y = cgetg(lx+1, t_POL);
     800       23466 :   y[1] = evalsigne(1) | evalvarn(v);
     801       23466 :   y++;
     802      130418 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     803       23466 :   return normalizepol_lg(--y, lx+1);
     804             : }
     805             : 
     806             : /* matrix whose entries are given by the coeffs of the polynomial v in
     807             :  * two variables (considered as degree n-1 polynomials) */
     808             : GEN
     809         322 : RgXX_to_RgM(GEN v, long n)
     810             : {
     811         322 :   long j, N = lg(v)-1;
     812         322 :   GEN y = cgetg(N, t_MAT);
     813        1043 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     814         322 :   return y;
     815             : }
     816             : 
     817             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     818             : GEN
     819       31308 : RgXY_swapspec(GEN x, long n, long w, long nx)
     820             : {
     821       31308 :   long j, ly = n+3;
     822       31308 :   GEN y = cgetg(ly, t_POL);
     823       31308 :   y[1] = evalsigne(1);
     824      384972 :   for (j=2; j<ly; j++)
     825             :   {
     826             :     long k;
     827      353664 :     GEN a = cgetg(nx+2,t_POL);
     828      353664 :     a[1] = evalsigne(1) | evalvarn(w);
     829     1809625 :     for (k=0; k<nx; k++)
     830             :     {
     831     1455961 :       GEN xk = gel(x,k);
     832     1455961 :       if (typ(xk)==t_POL)
     833     1367445 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     834             :       else
     835       88516 :         gel(a,k+2) = j==2 ? xk: gen_0;
     836             :     }
     837      353664 :     gel(y,j) = normalizepol_lg(a, nx+2);
     838             :   }
     839       31308 :   return normalizepol_lg(y,ly);
     840             : }
     841             : 
     842             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     843             : GEN
     844         952 : RgXY_swap(GEN x, long n, long w)
     845             : {
     846         952 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     847         952 :   setvarn(z, varn(x)); return z;
     848             : }
     849             : 
     850             : long
     851         288 : RgXY_degreex(GEN b)
     852             : {
     853         288 :   long deg = 0, i;
     854         288 :   if (!signe(b)) return -1;
     855        1396 :   for (i = 2; i < lg(b); ++i)
     856             :   {
     857        1108 :     GEN bi = gel(b, i);
     858        1108 :     if (typ(bi) == t_POL)
     859        1030 :       deg = maxss(deg, degpol(bi));
     860             :   }
     861         288 :   return deg;
     862             : }
     863             : 
     864             : GEN
     865       38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
     866             : 
     867             : /* return (x % X^n). Shallow */
     868             : GEN
     869     7927068 : RgXn_red_shallow(GEN a, long n)
     870             : {
     871     7927068 :   long i, L = n+2, l = lg(a);
     872             :   GEN  b;
     873     7927068 :   if (L >= l) return a; /* deg(x) < n */
     874     5735935 :   b = cgetg(L, t_POL); b[1] = a[1];
     875    36582389 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     876     5735937 :   return normalizepol_lg(b,L);
     877             : }
     878             : 
     879             : GEN
     880         483 : RgXnV_red_shallow(GEN x, long n)
     881        2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
     882             : 
     883             : /* return (x * X^n). Shallow */
     884             : GEN
     885   168706865 : RgX_shift_shallow(GEN a, long n)
     886             : {
     887   168706865 :   long i, l = lg(a);
     888             :   GEN  b;
     889   168706865 :   if (l == 2 || !n) return a;
     890   110120357 :   l += n;
     891   110120357 :   if (n < 0)
     892             :   {
     893    55434364 :     if (l <= 2) return pol_0(varn(a));
     894    53963109 :     b = cgetg(l, t_POL); b[1] = a[1];
     895    53964050 :     a -= n;
     896   165082105 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     897             :   } else {
     898    54685993 :     b = cgetg(l, t_POL); b[1] = a[1];
     899    54689993 :     a -= n; n += 2;
     900   118713375 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     901   210435458 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     902             :   }
     903   108654043 :   return b;
     904             : }
     905             : /* return (x * X^n). */
     906             : GEN
     907     1567501 : RgX_shift(GEN a, long n)
     908             : {
     909     1567501 :   long i, l = lg(a);
     910             :   GEN  b;
     911     1567501 :   if (l == 2 || !n) return RgX_copy(a);
     912     1567137 :   l += n;
     913     1567137 :   if (n < 0)
     914             :   {
     915         847 :     if (l <= 2) return pol_0(varn(a));
     916         777 :     b = cgetg(l, t_POL); b[1] = a[1];
     917         777 :     a -= n;
     918        2947 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     919             :   } else {
     920     1566290 :     b = cgetg(l, t_POL); b[1] = a[1];
     921     1566290 :     a -= n; n += 2;
     922     3915481 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     923     4318405 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     924             :   }
     925     1567067 :   return b;
     926             : }
     927             : 
     928             : GEN
     929      317037 : RgX_rotate_shallow(GEN P, long k, long p)
     930             : {
     931      317037 :   long i, l = lgpol(P);
     932             :   GEN r;
     933      317037 :   if (signe(P)==0)
     934        1365 :     return pol_0(varn(P));
     935      315672 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     936     2100644 :   for(i=0; i<p; i++)
     937             :   {
     938     1784972 :     long s = 2+(i+k)%p;
     939     1784972 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
     940             :   }
     941      315672 :   return RgX_renormalize(r);
     942             : }
     943             : 
     944             : GEN
     945     2984903 : RgX_mulXn(GEN x, long d)
     946             : {
     947             :   pari_sp av;
     948             :   GEN z;
     949             :   long v;
     950     2984903 :   if (d >= 0) return RgX_shift(x, d);
     951     1472286 :   d = -d;
     952     1472286 :   v = RgX_val(x);
     953     1472286 :   if (v >= d) return RgX_shift(x, -d);
     954     1472279 :   av = avma;
     955     1472279 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
     956     1472279 :   return gerepileupto(av, z);
     957             : }
     958             : 
     959             : long
     960         588 : RgXV_maxdegree(GEN x)
     961             : {
     962         588 :   long d = -1, i, l = lg(x);
     963        4494 :   for (i = 1; i < l; i++)
     964        3906 :     d = maxss(d, degpol(gel(x,i)));
     965         588 :   return d;
     966             : }
     967             : 
     968             : long
     969     3647031 : RgX_val(GEN x)
     970             : {
     971     3647031 :   long i, lx = lg(x);
     972     3647031 :   if (lx == 2) return LONG_MAX;
     973     4576960 :   for (i = 2; i < lx; i++)
     974     4576904 :     if (!isexactzero(gel(x,i))) break;
     975     3646870 :   if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
     976     3646814 :   return i - 2;
     977             : }
     978             : long
     979    82046386 : RgX_valrem(GEN x, GEN *Z)
     980             : {
     981    82046386 :   long v, i, lx = lg(x);
     982    82046386 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
     983   127784021 :   for (i = 2; i < lx; i++)
     984   127783401 :     if (!isexactzero(gel(x,i))) break;
     985             :   /* possible with nonrational zeros */
     986    82046429 :   if (i == lx)
     987             :   {
     988          14 :     *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
     989          14 :     return LONG_MAX;
     990             :   }
     991    82046415 :   v = i - 2;
     992    82046415 :   *Z = RgX_shift_shallow(x, -v);
     993    82056600 :   return v;
     994             : }
     995             : long
     996      851424 : RgX_valrem_inexact(GEN x, GEN *Z)
     997             : {
     998             :   long v;
     999      851424 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
    1000      868404 :   for (v = 0;; v++)
    1001      868404 :     if (!gequal0(gel(x,2+v))) break;
    1002      851417 :   if (Z) *Z = RgX_shift_shallow(x, -v);
    1003      851417 :   return v;
    1004             : }
    1005             : 
    1006             : GEN
    1007       63875 : RgXQC_red(GEN x, GEN T)
    1008      394877 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
    1009             : 
    1010             : GEN
    1011        1204 : RgXQV_red(GEN x, GEN T)
    1012       27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
    1013             : 
    1014             : GEN
    1015       12257 : RgXQM_red(GEN x, GEN T)
    1016       76132 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
    1017             : 
    1018             : GEN
    1019         322 : RgXQM_mul(GEN P, GEN Q, GEN T)
    1020             : {
    1021         322 :   return RgXQM_red(RgM_mul(P, Q), T);
    1022             : }
    1023             : 
    1024             : GEN
    1025      457326 : RgXQX_red(GEN P, GEN T)
    1026             : {
    1027      457326 :   long i, l = lg(P);
    1028      457326 :   GEN Q = cgetg(l, t_POL);
    1029      457326 :   Q[1] = P[1];
    1030     2388018 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
    1031      457325 :   return normalizepol_lg(Q, l);
    1032             : }
    1033             : 
    1034             : GEN
    1035      819943 : RgX_deriv(GEN x)
    1036             : {
    1037      819943 :   long i,lx = lg(x)-1;
    1038             :   GEN y;
    1039             : 
    1040      819943 :   if (lx<3) return pol_0(varn(x));
    1041      816856 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
    1042     3618572 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
    1043      816851 :   y[1] = x[1]; return normalizepol_lg(y,i);
    1044             : }
    1045             : 
    1046             : GEN
    1047     1379583 : RgX_recipspec_shallow(GEN x, long l, long n)
    1048             : {
    1049             :   long i;
    1050     1379583 :   GEN z = cgetg(n+2,t_POL);
    1051     1379584 :   z[1] = 0; z += 2;
    1052    36094564 :   for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
    1053     1602779 :   for(   ; i<n; i++) gel(z, n-i-1) = gen_0;
    1054     1379584 :   return normalizepol_lg(z-2,n+2);
    1055             : }
    1056             : 
    1057             : GEN
    1058      639126 : RgXn_recip_shallow(GEN P, long n)
    1059             : {
    1060      639126 :   GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
    1061      639128 :   setvarn(Q, varn(P));
    1062      639128 :   return Q;
    1063             : }
    1064             : 
    1065             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1066             : GEN
    1067       31150 : RgX_recip(GEN x)
    1068             : {
    1069             :   long lx, i, j;
    1070       31150 :   GEN y = cgetg_copy(x, &lx);
    1071      272104 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
    1072       31150 :   return normalizepol_lg(y,lx);
    1073             : }
    1074             : /* shallow version */
    1075             : GEN
    1076       57568 : RgX_recip_shallow(GEN x)
    1077             : {
    1078             :   long lx, i, j;
    1079       57568 :   GEN y = cgetg_copy(x, &lx);
    1080      345247 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1081       57568 :   return normalizepol_lg(y,lx);
    1082             : }
    1083             : 
    1084             : GEN
    1085     3650583 : RgX_recip_i(GEN x)
    1086             : {
    1087             :   long lx, i, j;
    1088     3650583 :   GEN y = cgetg_copy(x, &lx);
    1089    19470815 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1090     3650585 :   return y;
    1091             : }
    1092             : /*******************************************************************/
    1093             : /*                                                                 */
    1094             : /*                      ADDITION / SUBTRACTION                     */
    1095             : /*                                                                 */
    1096             : /*******************************************************************/
    1097             : /* same variable */
    1098             : GEN
    1099    77269036 : RgX_add(GEN x, GEN y)
    1100             : {
    1101    77269036 :   long i, lx = lg(x), ly = lg(y);
    1102             :   GEN z;
    1103    77269036 :   if (ly <= lx) {
    1104    67041694 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1105   265604514 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1106   118914669 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1107    67017653 :     z = normalizepol_lg(z, lx);
    1108             :   } else {
    1109    10227342 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1110    43632369 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1111    31467557 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1112    10228634 :     z = normalizepol_lg(z, ly);
    1113             :   }
    1114    77260483 :   return z;
    1115             : }
    1116             : GEN
    1117    63770791 : RgX_sub(GEN x, GEN y)
    1118             : {
    1119    63770791 :   long i, lx = lg(x), ly = lg(y);
    1120             :   GEN z;
    1121    63770791 :   if (ly <= lx) {
    1122    30129335 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1123   157208694 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1124    56187681 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1125    30105791 :     z = normalizepol_lg(z, lx);
    1126             :   } else {
    1127    33641456 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1128   125123788 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1129    74479204 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1130    33608784 :     z = normalizepol_lg(z, ly);
    1131             :   }
    1132    63765985 :   return z;
    1133             : }
    1134             : GEN
    1135     5211025 : RgX_neg(GEN x)
    1136             : {
    1137     5211025 :   long i, lx = lg(x);
    1138     5211025 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    1139    40419180 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    1140     5210927 :   return y;
    1141             : }
    1142             : 
    1143             : GEN
    1144    18864059 : RgX_Rg_add(GEN y, GEN x)
    1145             : {
    1146             :   GEN z;
    1147    18864059 :   long lz = lg(y), i;
    1148    18864059 :   if (lz == 2) return scalarpol(x,varn(y));
    1149    15824372 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1150    15824331 :   gel(z,2) = gadd(gel(y,2),x);
    1151    61643132 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1152             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1153             :    * an inexact 0 */
    1154    15824379 :   return normalizepol_lg(z,lz);
    1155             : }
    1156             : GEN
    1157       64857 : RgX_Rg_add_shallow(GEN y, GEN x)
    1158             : {
    1159             :   GEN z;
    1160       64857 :   long lz = lg(y), i;
    1161       64857 :   if (lz == 2) return scalarpol(x,varn(y));
    1162       64857 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1163       64857 :   gel(z,2) = gadd(gel(y,2),x);
    1164      129847 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1165       64858 :   return normalizepol_lg(z,lz);
    1166             : }
    1167             : GEN
    1168      182649 : RgX_Rg_sub(GEN y, GEN x)
    1169             : {
    1170             :   GEN z;
    1171      182649 :   long lz = lg(y), i;
    1172      182649 :   if (lz == 2)
    1173             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1174        2548 :     long v = varn(y);
    1175        2548 :     if (isrationalzero(x)) return pol_0(v);
    1176          14 :     z = cgetg(3,t_POL);
    1177          14 :     z[1] = gequal0(x)? evalvarn(v)
    1178          14 :                    : evalvarn(v) | evalsigne(1);
    1179          14 :     gel(z,2) = gneg(x); return z;
    1180             :   }
    1181      180101 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1182      180101 :   gel(z,2) = gsub(gel(y,2),x);
    1183      468957 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1184      180103 :   return normalizepol_lg(z,lz);
    1185             : }
    1186             : GEN
    1187     4841475 : Rg_RgX_sub(GEN x, GEN y)
    1188             : {
    1189             :   GEN z;
    1190     4841475 :   long lz = lg(y), i;
    1191     4841475 :   if (lz == 2) return scalarpol(x,varn(y));
    1192     4840397 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1193     4840385 :   gel(z,2) = gsub(x, gel(y,2));
    1194     7258217 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1195     4840359 :   return normalizepol_lg(z,lz);
    1196             : }
    1197             : /*******************************************************************/
    1198             : /*                                                                 */
    1199             : /*                  KARATSUBA MULTIPLICATION                       */
    1200             : /*                                                                 */
    1201             : /*******************************************************************/
    1202             : #if 0
    1203             : /* to debug Karatsuba-like routines */
    1204             : GEN
    1205             : zx_debug_spec(GEN x, long nx)
    1206             : {
    1207             :   GEN z = cgetg(nx+2,t_POL);
    1208             :   long i;
    1209             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1210             :   z[1] = evalsigne(1); return z;
    1211             : }
    1212             : 
    1213             : GEN
    1214             : RgX_debug_spec(GEN x, long nx)
    1215             : {
    1216             :   GEN z = cgetg(nx+2,t_POL);
    1217             :   long i;
    1218             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1219             :   z[1] = evalsigne(1); return z;
    1220             : }
    1221             : #endif
    1222             : 
    1223             : /* generic multiplication */
    1224             : GEN
    1225     8675134 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1226             : {
    1227             :   GEN z, t;
    1228             :   long i;
    1229     8675134 :   if (nx == ny) {
    1230     1332260 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1231     4255070 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1232     1332251 :     return normalizepol_lg(z, nx+2);
    1233             :   }
    1234     7342874 :   if (ny < nx) {
    1235     7156412 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1236    26343820 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1237    17115067 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1238     7155942 :     return normalizepol_lg(z, nx+2);
    1239             :   } else {
    1240      186462 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1241     3673864 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1242      451157 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1243      186473 :     return normalizepol_lg(z, ny+2);
    1244             :   }
    1245             : }
    1246             : GEN
    1247      223331 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1248             : {
    1249             :   GEN z, t;
    1250             :   long i;
    1251      223331 :   if (nx == ny) {
    1252       12810 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1253     2183524 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1254       12810 :     return normalizepol_lg(z, nx+2);
    1255             :   }
    1256      210521 :   if (ny < nx) {
    1257      208743 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1258     3725832 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1259     2381590 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1260      208743 :     return normalizepol_lg(z, nx+2);
    1261             :   } else {
    1262        1778 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1263      330904 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1264       12222 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1265        1778 :     return normalizepol_lg(z, ny+2);
    1266             :   }
    1267             : }
    1268             : 
    1269             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1270             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1271             :  * t_VECSMALL, to hold this, which can be left on stack: gerepile
    1272             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1273             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1274             : static GEN
    1275    56452219 : RgXspec_kill0(GEN x, long lx)
    1276             : {
    1277    56452219 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
    1278             :   long i;
    1279   191740885 :   for (i=0; i <lx; i++)
    1280             :   {
    1281   135288671 :     GEN c = gel(x,i);
    1282   135288671 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1283             :   }
    1284    56452214 :   return z;
    1285             : }
    1286             : 
    1287             : INLINE GEN
    1288    94119896 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1289             : {
    1290    94119896 :   pari_sp av = avma;
    1291    94119896 :   GEN s = NULL;
    1292             :   long i;
    1293             : 
    1294   356926450 :   for (i=a; i<b; i++)
    1295   262815472 :     if (gel(y,i) && gel(x,-i))
    1296             :     {
    1297   184782160 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1298   184778093 :       s = s? gadd(s, t): t;
    1299             :     }
    1300    94110978 :   return s? gerepileupto(av, s): gen_0;
    1301             : }
    1302             : 
    1303             : /* assume nx >= ny > 0, return x * y * t^v */
    1304             : static GEN
    1305    21593956 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1306             : {
    1307             :   long i, lz, nz;
    1308             :   GEN z;
    1309             : 
    1310    21593956 :   x = RgXspec_kill0(x,nx);
    1311    21593902 :   y = RgXspec_kill0(y,ny);
    1312    21593921 :   lz = nx + ny + 1; nz = lz-2;
    1313    21593921 :   lz += v;
    1314    21593921 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1315    33327547 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1316    57235856 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1317    41189832 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1318    35641958 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1319    21593397 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1320             : }
    1321             : 
    1322             : /* return (x * X^d) + y. Assume d > 0 */
    1323             : GEN
    1324     9294633 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1325             : {
    1326             :   GEN x, y, xd, yd, zd;
    1327             :   long a, lz, nx, ny;
    1328             : 
    1329     9294633 :   if (!signe(x0)) return y0;
    1330     9188666 :   ny = lgpol(y0);
    1331     9188664 :   nx = lgpol(x0);
    1332     9188755 :   zd = (GEN)avma;
    1333     9188755 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1334     9188755 :   if (a <= 0)
    1335             :   {
    1336     1471692 :     lz = nx+d+2;
    1337     1471692 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1338     3382330 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1339     1471692 :     x = zd + a;
    1340     1489393 :     while (zd > x) gel(--zd,0) = gen_0;
    1341             :   }
    1342             :   else
    1343             :   {
    1344     7717063 :     xd = new_chunk(d); yd = y+d;
    1345     7717069 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1346     7716991 :     lz = (a>nx)? ny+2: lg(x)+d;
    1347    38782112 :     x += 2; while (xd > x) *--zd = *--xd;
    1348             :   }
    1349    22490245 :   while (yd > y) *--zd = *--yd;
    1350     9188683 :   *--zd = x0[1];
    1351     9188683 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1352             : }
    1353             : GEN
    1354      515890 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1355             : {
    1356             :   GEN x, y, xd, yd, zd;
    1357             :   long a, lz, nx, ny;
    1358             : 
    1359      515890 :   if (!signe(x0)) return RgX_copy(y0);
    1360      515120 :   nx = lgpol(x0);
    1361      515120 :   ny = lgpol(y0);
    1362      515120 :   zd = (GEN)avma;
    1363      515120 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1364      515120 :   if (a <= 0)
    1365             :   {
    1366      291789 :     lz = nx+d+2;
    1367      291789 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1368     4294362 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1369      291789 :     x = zd + a;
    1370      757955 :     while (zd > x) gel(--zd,0) = gen_0;
    1371             :   }
    1372             :   else
    1373             :   {
    1374      223331 :     xd = new_chunk(d); yd = y+d;
    1375      223331 :     x = RgX_addspec(x,yd, nx,a);
    1376      223331 :     lz = (a>nx)? ny+2: lg(x)+d;
    1377     8423551 :     x += 2; while (xd > x) *--zd = *--xd;
    1378             :   }
    1379     2622449 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1380      515120 :   *--zd = x0[1];
    1381      515120 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1382             : }
    1383             : 
    1384             : /* return x * y mod t^n */
    1385             : static GEN
    1386     6627805 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1387             : {
    1388     6627805 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1389             :   GEN z;
    1390     6627805 :   if (lx < 0) return pol_0(varn(x));
    1391     6627805 :   if (ly < 0) return pol_0(varn(x));
    1392     6627805 :   z = cgetg(lz, t_POL) + 2;
    1393     6627805 :   x+=2; if (lx > n) lx = n;
    1394     6627805 :   y+=2; if (ly > n) ly = n;
    1395     6627805 :   z[-1] = x[-1];
    1396     6627805 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1397     6627805 :   x = RgXspec_kill0(x, lx);
    1398     6627805 :   y = RgXspec_kill0(y, ly);
    1399             :   /* x:y:z [i] = term of degree i */
    1400    26217431 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1401    11825932 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1402     6673807 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1403     6627805 :   return normalizepol_lg(z - 2, lz);
    1404             : }
    1405             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1406             : static GEN
    1407     9413646 : RgXn_mul2(GEN f, GEN g, long n)
    1408             : {
    1409     9413646 :   pari_sp av = avma;
    1410             :   GEN fe,fo, ge,go, l,h,m;
    1411             :   long n0, n1;
    1412     9413646 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1413     6663225 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1414       35420 :   n0 = n>>1; n1 = n-n0;
    1415       35420 :   RgX_even_odd(f, &fe, &fo);
    1416       35420 :   RgX_even_odd(g, &ge, &go);
    1417       35420 :   l = RgXn_mul2(fe,ge,n1);
    1418       35420 :   h = RgXn_mul2(fo,go,n0);
    1419       35420 :   m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1420             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1421             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1422       35420 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1423             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1424       35420 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1425       35420 :   m = RgX_inflate(m,2);
    1426             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1427       35420 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1428       35420 :   h = RgX_inflate(h,2);
    1429       35420 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1430       35420 :   return gerepileupto(av, h);
    1431             : }
    1432             : /* (f*g) \/ x^n */
    1433             : static GEN
    1434     1589258 : RgX_mulhigh_i2(GEN f, GEN g, long n)
    1435             : {
    1436     1589258 :   long d = degpol(f)+degpol(g) + 1 - n;
    1437             :   GEN h;
    1438     1589258 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1439       29654 :   h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
    1440             :                                   RgX_recip_i(g), d));
    1441       29654 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1442             : }
    1443             : 
    1444             : /* (f*g) \/ x^n */
    1445             : static GEN
    1446           0 : RgX_sqrhigh_i2(GEN f, long n)
    1447             : {
    1448           0 :   long d = 2*degpol(f)+ 1 - n;
    1449             :   GEN h;
    1450           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1451           0 :   h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
    1452           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1453             : }
    1454             : 
    1455             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1456             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1457             :  * Only c, c0, c1, c2 are genuine GEN.
    1458             :  */
    1459             : GEN
    1460    22322984 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1461             : {
    1462             :   GEN a0, c, c0;
    1463    22322984 :   long n0, n0a, i, v = 0;
    1464             :   pari_sp av;
    1465             : 
    1466    28706525 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1467    28207214 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1468    22322959 :   if (na < nb) swapspec(a,b, na,nb);
    1469    22322959 :   if (!nb) return pol_0(0);
    1470             : 
    1471    22074367 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1472      480416 :   RgX_shift_inplace_init(v);
    1473      480431 :   i = (na>>1); n0 = na-i; na = i;
    1474      480431 :   av = avma; a0 = a+n0; n0a = n0;
    1475     1335947 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1476             : 
    1477      480431 :   if (nb > n0)
    1478             :   {
    1479             :     GEN b0,c1,c2;
    1480             :     long n0b;
    1481             : 
    1482      479032 :     nb -= n0; b0 = b+n0; n0b = n0;
    1483     1447011 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1484      479032 :     c = RgX_mulspec(a,b,n0a,n0b);
    1485      479032 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1486             : 
    1487      479032 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1488      479032 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1489             : 
    1490      479032 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1491      479032 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1492      479032 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1493             :   }
    1494             :   else
    1495             :   {
    1496        1399 :     c = RgX_mulspec(a,b,n0a,nb);
    1497        1399 :     c0 = RgX_mulspec(a0,b,na,nb);
    1498             :   }
    1499      480431 :   c0 = RgX_addmulXn(c0,c,n0);
    1500      480431 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1501             : }
    1502             : 
    1503             : INLINE GEN
    1504       50452 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1505             : {
    1506       50452 :   pari_sp av = avma;
    1507       50452 :   GEN s = NULL;
    1508       50452 :   long j, l = (i+1)>>1;
    1509      139622 :   for (j=a; j<l; j++)
    1510             :   {
    1511       89170 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1512       89170 :     if (xj && xx)
    1513             :     {
    1514       81526 :       GEN t = gmul(xj, xx);
    1515       81526 :       s = s? gadd(s, t): t;
    1516             :     }
    1517             :   }
    1518       50452 :   if (s) s = gshift(s,1);
    1519       50452 :   if ((i&1) == 0)
    1520             :   {
    1521       29673 :     GEN t = gel(x, i>>1);
    1522       29673 :     if (t) {
    1523       27279 :       t = gsqr(t);
    1524       27279 :       s = s? gadd(s, t): t;
    1525             :     }
    1526             :   }
    1527       50452 :   return s? gerepileupto(av,s): gen_0;
    1528             : }
    1529             : static GEN
    1530        8894 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1531             : {
    1532             :   long i, lz, nz;
    1533             :   GEN z;
    1534             : 
    1535        8894 :   if (!nx) return pol_0(0);
    1536        8894 :   x = RgXspec_kill0(x,nx);
    1537        8894 :   lz = (nx << 1) + 1, nz = lz-2;
    1538        8894 :   lz += v;
    1539        8894 :   z = cgetg(lz,t_POL) + 2;
    1540       13066 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1541       38567 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1542       29673 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1543        8894 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1544             : }
    1545             : /* return x^2 mod t^n */
    1546             : static GEN
    1547           0 : RgXn_sqr_basecase(GEN x, long n)
    1548             : {
    1549           0 :   long i, lz = n+2, lx = lgpol(x);
    1550             :   GEN z;
    1551           0 :   if (lx < 0) return pol_0(varn(x));
    1552           0 :   z = cgetg(lz, t_POL);
    1553           0 :   z[1] = x[1];
    1554           0 :   x+=2; if (lx > n) lx = n;
    1555           0 :   x = RgXspec_kill0(x,lx);
    1556           0 :   z+=2;/* x:z [i] = term of degree i */
    1557           0 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1558           0 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1559           0 :   z -= 2; return normalizepol_lg(z, lz);
    1560             : }
    1561             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1562             : static GEN
    1563         315 : RgXn_sqr2(GEN f, long n)
    1564             : {
    1565         315 :   pari_sp av = avma;
    1566             :   GEN fe,fo, l,h,m;
    1567             :   long n0, n1;
    1568         315 :   if (2*degpol(f) < n) return RgX_sqr_i(f);
    1569           0 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1570           0 :   n0 = n>>1; n1 = n-n0;
    1571           0 :   RgX_even_odd(f, &fe, &fo);
    1572           0 :   l = RgXn_sqr(fe,n1);
    1573           0 :   h = RgXn_sqr(fo,n0);
    1574           0 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1575             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1576             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1577           0 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1578             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1579           0 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1580           0 :   m = RgX_inflate(m,2);
    1581             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1582           0 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1583           0 :   h = RgX_inflate(h,2);
    1584           0 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1585           0 :   return gerepileupto(av, h);
    1586             : }
    1587             : GEN
    1588        8933 : RgX_sqrspec(GEN a, long na)
    1589             : {
    1590             :   GEN a0, c, c0, c1;
    1591        8933 :   long n0, n0a, i, v = 0;
    1592             :   pari_sp av;
    1593             : 
    1594       11019 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1595        8933 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1596          39 :   RgX_shift_inplace_init(v);
    1597          39 :   i = (na>>1); n0 = na-i; na = i;
    1598          39 :   av = avma; a0 = a+n0; n0a = n0;
    1599          39 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1600             : 
    1601          39 :   c = RgX_sqrspec(a,n0a);
    1602          39 :   c0 = RgX_sqrspec(a0,na);
    1603          39 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1604          39 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1605          39 :   c0 = RgX_addmulXn(c0,c,n0);
    1606          39 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1607             : }
    1608             : 
    1609             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1610             : GEN
    1611     1716966 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1612             : {
    1613     1716966 :   GEN z = RgX_mul(A, B);
    1614     1716956 :   if (a < b)
    1615       10414 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1616     1706542 :   else if (a > b)
    1617     1104512 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1618             :   else
    1619      602030 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1620     1716962 :   return z;
    1621             : }
    1622             : 
    1623             : GEN
    1624    20883108 : RgX_mul_i(GEN x, GEN y)
    1625             : {
    1626    20883108 :   GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
    1627    20882720 :   setvarn(z, varn(x)); return z;
    1628             : }
    1629             : 
    1630             : GEN
    1631        8855 : RgX_sqr_i(GEN x)
    1632             : {
    1633        8855 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1634        8855 :   setvarn(z,varn(x)); return z;
    1635             : }
    1636             : 
    1637             : /*******************************************************************/
    1638             : /*                                                                 */
    1639             : /*                               DIVISION                          */
    1640             : /*                                                                 */
    1641             : /*******************************************************************/
    1642             : GEN
    1643     3950300 : RgX_Rg_divexact(GEN x, GEN y) {
    1644     3950300 :   long i, lx = lg(x);
    1645             :   GEN z;
    1646     3950300 :   if (lx == 2) return gcopy(x);
    1647     3900525 :   switch(typ(y))
    1648             :   {
    1649     3870206 :     case t_INT:
    1650     3870206 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1651     3695933 :       break;
    1652        5215 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1653             :   }
    1654     3721037 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1655    23823958 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1656     3719107 :   return z;
    1657             : }
    1658             : GEN
    1659    29228754 : RgX_Rg_div(GEN x, GEN y) {
    1660    29228754 :   long i, lx = lg(x);
    1661             :   GEN z;
    1662    29228754 :   if (lx == 2) return gcopy(x);
    1663    28968696 :   switch(typ(y))
    1664             :   {
    1665    19775868 :     case t_INT:
    1666    19775868 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1667     3735406 :       break;
    1668       13125 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1669             :   }
    1670    12915109 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1671    47841032 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1672    12914732 :   return normalizepol_lg(z, lx);
    1673             : }
    1674             : GEN
    1675       38906 : RgX_normalize(GEN x)
    1676             : {
    1677       38906 :   GEN z, d = NULL;
    1678       38906 :   long i, n = lg(x)-1;
    1679       38906 :   for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
    1680       38906 :   if (i == 1) return pol_0(varn(x));
    1681       38906 :   if (i == n && isint1(d)) return x;
    1682       17584 :   n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
    1683       31570 :   for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
    1684       17584 :   gel(z,n) = Rg_get_1(d); return z;
    1685             : }
    1686             : GEN
    1687       10423 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
    1688             : GEN
    1689      248902 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1690             : {
    1691      248902 :   long l = lg(a), i;
    1692             :   GEN a0, z0, z;
    1693             : 
    1694      248902 :   if (l <= 3)
    1695             :   {
    1696           0 :     if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
    1697           0 :     return pol_0(varn(a));
    1698             :   }
    1699      248902 :   z = cgetg(l-1, t_POL);
    1700      248905 :   z[1] = a[1];
    1701      248905 :   a0 = a + l-1;
    1702      248905 :   z0 = z + l-2; *z0 = *a0--;
    1703     2948326 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1704             :   {
    1705     2699426 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1706     2699421 :     gel(z0,0) = t;
    1707             :   }
    1708      248900 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1709      248900 :   return z;
    1710             : }
    1711             : /* Polynomial division x / y:
    1712             :  *   if pr = ONLY_REM return remainder, otherwise return quotient
    1713             :  *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
    1714             :  *   if pr != NULL set *pr to remainder, as the last object on stack */
    1715             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1716             : static GEN
    1717    23330712 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
    1718             : {
    1719             :   pari_sp avy, av, av1;
    1720             :   long dx,dy,dz,i,j,sx,lr;
    1721             :   GEN z,p1,p2,rem,y_lead,mod,p;
    1722             :   GEN (*f)(GEN,GEN);
    1723             : 
    1724    23330712 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1725             : 
    1726    23330712 :   dy = degpol(y);
    1727    23330693 :   y_lead = gel(y,dy+2);
    1728    23330693 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1729             :   {
    1730           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1731           0 :     for (dy--; dy>=0; dy--)
    1732             :     {
    1733           0 :       y_lead = gel(y,dy+2);
    1734           0 :       if (!gequal0(y_lead)) break;
    1735             :     }
    1736             :   }
    1737    23330682 :   if (!dy) /* y is constant */
    1738             :   {
    1739        6801 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1740        6794 :     z = RgX_Rg_div(x, y_lead);
    1741        6794 :     if (pr == ONLY_DIVIDES) return z;
    1742        2580 :     if (pr) *pr = pol_0(varn(x));
    1743        2580 :     return z;
    1744             :   }
    1745    23323881 :   dx = degpol(x);
    1746    23323857 :   if (dx < dy)
    1747             :   {
    1748     3847577 :     if (pr == ONLY_REM) return RgX_copy(x);
    1749      374687 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1750      374666 :     z = pol_0(varn(x));
    1751      374666 :     if (pr) *pr = RgX_copy(x);
    1752      374666 :     return z;
    1753             :   }
    1754             : 
    1755             :   /* x,y in R[X], y non constant */
    1756    19476280 :   av = avma;
    1757    19476280 :   p = NULL;
    1758    19476280 :   if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1759             :   {
    1760      238784 :     z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
    1761      238784 :     if (!z) return gc_NULL(av);
    1762      238784 :     z = FpX_to_mod(z, p);
    1763      238784 :     if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
    1764      122528 :       return gerepileupto(av, z);
    1765      116256 :     *pr = FpX_to_mod(*pr, p);
    1766      116256 :     return gc_all(av, 2, &z, pr);
    1767             :   }
    1768    19237616 :   switch(typ(y_lead))
    1769             :   {
    1770           0 :     case t_REAL:
    1771           0 :       y_lead = ginv(y_lead);
    1772           0 :       f = gmul; mod = NULL;
    1773           0 :       break;
    1774        2151 :     case t_INTMOD:
    1775        2151 :     case t_POLMOD: y_lead = ginv(y_lead);
    1776        2151 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1777        2151 :       break;
    1778    19235465 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1779    19235461 :       f = gdiv; mod = NULL;
    1780             :   }
    1781             : 
    1782    19237612 :   if (y_lead == NULL)
    1783    17409620 :     p2 = gel(x,dx+2);
    1784             :   else {
    1785             :     for(;;) {
    1786     1827992 :       p2 = f(gel(x,dx+2),y_lead);
    1787     1827990 :       p2 = simplify_shallow(p2);
    1788     1827990 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1789             :     }
    1790     1827990 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1791             :     {
    1792           0 :       if (pr == ONLY_DIVIDES) {
    1793           0 :         set_avma(av);
    1794           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1795             :       }
    1796           0 :       if (pr == ONLY_REM)
    1797             :       {
    1798           0 :         if (dx < 0)
    1799           0 :           return gerepilecopy(av, scalarpol(p2, varn(x)));
    1800             :         else
    1801             :         {
    1802             :           GEN t;
    1803           0 :           set_avma(av);
    1804           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1805           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1806           0 :           return t;
    1807             :         }
    1808             :       }
    1809           0 :       if (pr) /* cf ONLY_REM above */
    1810             :       {
    1811           0 :         if (dx < 0)
    1812             :         {
    1813           0 :           p2 = gclone(p2);
    1814           0 :           set_avma(av);
    1815           0 :           z = pol_0(varn(x));
    1816           0 :           x = scalarpol(p2, varn(x));
    1817           0 :           gunclone(p2);
    1818             :         }
    1819             :         else
    1820             :         {
    1821             :           GEN t;
    1822           0 :           set_avma(av);
    1823           0 :           z = pol_0(varn(x));
    1824           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1825           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1826           0 :           x = t;
    1827             :         }
    1828           0 :         *pr = x;
    1829             :       }
    1830             :       else
    1831             :       {
    1832           0 :         set_avma(av);
    1833           0 :         z = pol_0(varn(x));
    1834             :       }
    1835           0 :       return z;
    1836             :     }
    1837             :   }
    1838             :   /* dx >= dy */
    1839    19237610 :   avy = avma;
    1840    19237610 :   dz = dx-dy;
    1841    19237610 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1842    19237599 :   x += 2;
    1843    19237599 :   z += 2;
    1844    19237599 :   y += 2;
    1845    19237599 :   gel(z,dz) = gcopy(p2);
    1846             : 
    1847    52567817 :   for (i=dx-1; i>=dy; i--)
    1848             :   {
    1849    33330421 :     av1=avma; p1=gel(x,i);
    1850  1136168286 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1851    33325361 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1852             : 
    1853    33325362 :     if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
    1854             :     else
    1855    23357388 :       p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
    1856    33329785 :     gel(z,i-dy) = p1;
    1857             :   }
    1858    19237396 :   if (!pr) return gerepileupto(av,z-2);
    1859             : 
    1860    12447913 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1861    12809592 :   for (sx=0; ; i--)
    1862             :   {
    1863    12809592 :     p1 = gel(x,i);
    1864             :     /* we always enter this loop at least once */
    1865    32142115 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1866    12808392 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1867    12808392 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
    1868     3584963 :     if (!isexactzero(p1)) break;
    1869     3434380 :     if (!i) break;
    1870      361594 :     set_avma(av1);
    1871             :   }
    1872    12447217 :   if (pr == ONLY_DIVIDES)
    1873             :   {
    1874        7917 :     if (sx) return gc_NULL(av);
    1875        7896 :     set_avma((pari_sp)rem); return gerepileupto(av,z-2);
    1876             :   }
    1877    12439300 :   lr=i+3; rem -= lr;
    1878    12439300 :   if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
    1879    12395681 :   else p1 = gerepileupto((pari_sp)rem,p1);
    1880    12439995 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    1881    12439995 :   rem[1] = z[-1];
    1882    12439995 :   rem += 2;
    1883    12439995 :   gel(rem,i) = p1;
    1884    17044798 :   for (i--; i>=0; i--)
    1885             :   {
    1886     4604807 :     av1=avma; p1 = gel(x,i);
    1887    13104456 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1888     4604358 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1889     4604569 :     gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
    1890             :   }
    1891    12439991 :   rem -= 2;
    1892    12439991 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1893    12440181 :   if (pr == ONLY_REM) return gerepileupto(av,rem);
    1894     6747108 :   z -= 2;
    1895             :   {
    1896     6747108 :     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
    1897     6747108 :     gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
    1898             :   }
    1899             : }
    1900             : 
    1901             : GEN
    1902    14164791 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1903             : {
    1904    14164791 :   if (pr == ONLY_REM) return RgX_rem(x, y);
    1905    14164791 :   return RgX_divrem_i(x, y, pr);
    1906             : }
    1907             : 
    1908             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1909             : GEN
    1910      136595 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1911             : {
    1912             :   long vx, dx, dy, dz, i, j, sx, lr;
    1913             :   pari_sp av0, av, tetpil;
    1914             :   GEN z,p1,rem,lead;
    1915             : 
    1916      136595 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1917      136595 :   vx = varn(x);
    1918      136595 :   dx = degpol(x);
    1919      136595 :   dy = degpol(y);
    1920      136595 :   if (dx < dy)
    1921             :   {
    1922       28288 :     if (pr)
    1923             :     {
    1924       28260 :       av0 = avma; x = RgXQX_red(x, T);
    1925       28260 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
    1926       28260 :       if (pr == ONLY_REM) return x;
    1927           0 :       *pr = x;
    1928             :     }
    1929          28 :     return pol_0(vx);
    1930             :   }
    1931      108307 :   lead = leading_coeff(y);
    1932      108307 :   if (!dy) /* y is constant */
    1933             :   {
    1934         602 :     if (pr && pr != ONLY_DIVIDES)
    1935             :     {
    1936           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1937           0 :       *pr = pol_0(vx);
    1938             :     }
    1939         602 :     if (gequal1(lead)) return RgX_copy(x);
    1940           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1941           0 :     return gerepile(av0,tetpil,RgXQX_red(x,T));
    1942             :   }
    1943      107705 :   av0 = avma; dz = dx-dy;
    1944      107705 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1945      107705 :   set_avma(av0);
    1946      107705 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1947      107705 :   x += 2; y += 2; z += 2;
    1948             : 
    1949      107705 :   p1 = gel(x,dx); av = avma;
    1950      107705 :   gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    1951      557728 :   for (i=dx-1; i>=dy; i--)
    1952             :   {
    1953      450023 :     av=avma; p1=gel(x,i);
    1954     2359397 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1955      449989 :     if (lead) p1 = gmul(grem(p1, T), lead);
    1956      449993 :     tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
    1957             :   }
    1958      107705 :   if (!pr) { guncloneNULL(lead); return z-2; }
    1959             : 
    1960      106431 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    1961      179483 :   for (sx=0; ; i--)
    1962             :   {
    1963      179483 :     p1 = gel(x,i);
    1964      619006 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1965      179481 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    1966      100135 :     if (!i) break;
    1967       73052 :     set_avma(av);
    1968             :   }
    1969      106431 :   if (pr == ONLY_DIVIDES)
    1970             :   {
    1971       27041 :     guncloneNULL(lead);
    1972       27041 :     if (sx) return gc_NULL(av0);
    1973       24535 :     return gc_const((pari_sp)rem, z-2);
    1974             :   }
    1975       79390 :   lr=i+3; rem -= lr;
    1976       79390 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    1977       79390 :   rem[1] = z[-1];
    1978       79390 :   p1 = gerepile((pari_sp)rem,tetpil,p1);
    1979       79390 :   rem += 2; gel(rem,i) = p1;
    1980      170189 :   for (i--; i>=0; i--)
    1981             :   {
    1982       90799 :     av=avma; p1 = gel(x,i);
    1983      307651 :     for (j=0; j<=i && j<=dz; j++)
    1984      216852 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1985       90799 :     tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
    1986             :   }
    1987       79390 :   rem -= 2;
    1988       79390 :   guncloneNULL(lead);
    1989       79390 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1990       79390 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
    1991         161 :   *pr = rem; return z-2;
    1992             : }
    1993             : 
    1994             : /*******************************************************************/
    1995             : /*                                                                 */
    1996             : /*                        PSEUDO-DIVISION                          */
    1997             : /*                                                                 */
    1998             : /*******************************************************************/
    1999             : INLINE GEN
    2000     1207242 : rem(GEN c, GEN T)
    2001             : {
    2002     1207242 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    2003     1207242 :   return c;
    2004             : }
    2005             : 
    2006             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    2007             : int
    2008       17055 : ZXQX_dvd(GEN x, GEN y, GEN T)
    2009             : {
    2010             :   long dx, dy, i, T_ismonic;
    2011       17055 :   pari_sp av = avma, av2;
    2012             :   GEN y_lead;
    2013             : 
    2014       17055 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    2015       17055 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2016       17055 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    2017             :   /* if monic, no point in using pseudo-division */
    2018       17055 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    2019       14542 :   T_ismonic = gequal1(leading_coeff(T));
    2020       14542 :   dx = degpol(x);
    2021       14542 :   if (dx < dy) return !signe(x);
    2022       14542 :   (void)new_chunk(2);
    2023       14542 :   x = RgX_recip_i(x)+2;
    2024       14542 :   y = RgX_recip_i(y)+2;
    2025             :   /* pay attention to sparse divisors */
    2026       29585 :   for (i = 1; i <= dy; i++)
    2027       15043 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    2028       14542 :   av2 = avma;
    2029             :   for (;;)
    2030       72589 :   {
    2031       87131 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    2032       87131 :     x0 = gneg(x0);
    2033       87131 :     m = gcdii(cx, y0);
    2034       87131 :     if (!equali1(m))
    2035             :     {
    2036       84636 :       x0 = gdiv(x0, m);
    2037       84636 :       y0 = diviiexact(y0, m);
    2038       84636 :       if (equali1(y0)) y0 = NULL;
    2039             :     }
    2040      178035 :     for (i=1; i<=dy; i++)
    2041             :     {
    2042       90904 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2043       90904 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    2044       90904 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2045       90904 :       gel(x,i) = c;
    2046             :     }
    2047      686642 :     for (   ; i<=dx; i++)
    2048             :     {
    2049      599511 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2050      599511 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2051      599511 :       gel(x,i) = c;
    2052             :     }
    2053      101467 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    2054       87131 :     if (dx < dy) break;
    2055       72589 :     if (gc_needed(av2,1))
    2056             :     {
    2057          28 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    2058          28 :       gerepilecoeffs(av2,x,dx+1);
    2059             :     }
    2060             :   }
    2061       14542 :   return gc_bool(av, dx < 0);
    2062             : }
    2063             : 
    2064             : /* T either NULL or a t_POL. */
    2065             : GEN
    2066      103327 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    2067             : {
    2068      103327 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    2069      103327 :   pari_sp av = avma, av2;
    2070             :   GEN y_lead;
    2071             : 
    2072      103327 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    2073      103327 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2074             :   /* if monic, no point in using pseudo-division */
    2075      103327 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    2076       97509 :   dx = degpol(x);
    2077       97509 :   if (dx < dy) return RgX_copy(x);
    2078       97509 :   (void)new_chunk(2);
    2079       97509 :   x = RgX_recip_i(x)+2;
    2080       97509 :   y = RgX_recip_i(y)+2;
    2081             :   /* pay attention to sparse divisors */
    2082      309017 :   for (i = 1; i <= dy; i++)
    2083      211508 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2084       97509 :   dz = dx-dy; p = dz+1;
    2085       97509 :   av2 = avma;
    2086             :   for (;;)
    2087             :   {
    2088      291243 :     gel(x,0) = gneg(gel(x,0)); p--;
    2089      935346 :     for (i=1; i<=dy; i++)
    2090             :     {
    2091      644096 :       GEN c = gmul(y_lead, gel(x,i));
    2092      644074 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    2093      644092 :       gel(x,i) = rem(c, T);
    2094             :     }
    2095      730244 :     for (   ; i<=dx; i++)
    2096             :     {
    2097      439000 :       GEN c = gmul(y_lead, gel(x,i));
    2098      438994 :       gel(x,i) = rem(c, T);
    2099             :     }
    2100      319949 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    2101      291244 :     if (dx < dy) break;
    2102      193734 :     if (gc_needed(av2,1))
    2103             :     {
    2104           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    2105           0 :       gerepilecoeffs(av2,x,dx+1);
    2106             :     }
    2107             :   }
    2108       97510 :   if (dx < 0) return pol_0(vx);
    2109       96327 :   lx = dx+3; x -= 2;
    2110       96327 :   x[0] = evaltyp(t_POL) | _evallg(lx);
    2111       96327 :   x[1] = evalsigne(1) | evalvarn(vx);
    2112       96327 :   x = RgX_recip_i(x);
    2113       96327 :   if (p)
    2114             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    2115       13348 :     GEN t = y_lead;
    2116       13348 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    2117           0 :       t = RgXQ_powu(t, p, T);
    2118             :     else
    2119       13348 :       t = gpowgs(t, p);
    2120       35718 :     for (i=2; i<lx; i++)
    2121             :     {
    2122       22370 :       GEN c = gmul(gel(x,i), t);
    2123       22370 :       gel(x,i) = rem(c,T);
    2124             :     }
    2125       13348 :     if (!T) return gerepileupto(av, x);
    2126             :   }
    2127       82979 :   return gerepilecopy(av, x);
    2128             : }
    2129             : 
    2130             : GEN
    2131      103327 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2132             : 
    2133             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2134             : GEN
    2135       12014 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2136             : {
    2137       12014 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2138       12014 :   pari_sp av = avma, av2;
    2139             :   GEN z, r, ypow, y_lead;
    2140             : 
    2141       12014 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2142       12014 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2143       12014 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2144        7589 :   dx = degpol(x);
    2145        7589 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2146        7589 :   if (dx == dy)
    2147             :   {
    2148          98 :     GEN x_lead = gel(x,lg(x)-1);
    2149          98 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2150          98 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2151          98 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2152          98 :     *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
    2153             :   }
    2154        7491 :   (void)new_chunk(2);
    2155        7491 :   x = RgX_recip_i(x)+2;
    2156        7491 :   y = RgX_recip_i(y)+2;
    2157             :   /* pay attention to sparse divisors */
    2158       39049 :   for (i = 1; i <= dy; i++)
    2159       31558 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2160        7491 :   dz = dx-dy; p = dz+1;
    2161        7491 :   lz = dz+3;
    2162        7491 :   z = cgetg(lz, t_POL);
    2163        7491 :   z[1] = evalsigne(1) | evalvarn(vx);
    2164       28745 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2165        7491 :   ypow = new_chunk(dz+1);
    2166        7491 :   gel(ypow,0) = gen_1;
    2167        7491 :   gel(ypow,1) = y_lead;
    2168       13763 :   for (i=2; i<=dz; i++)
    2169             :   {
    2170        6272 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2171        6272 :     gel(ypow,i) = rem(c,T);
    2172             :   }
    2173        7491 :   av2 = avma;
    2174        7491 :   for (iz=2;;)
    2175             :   {
    2176       15682 :     p--;
    2177       15682 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2178       70034 :     for (i=1; i<=dy; i++)
    2179             :     {
    2180       54352 :       GEN c = gmul(y_lead, gel(x,i));
    2181       54352 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2182       54352 :       gel(x,i) = rem(c, T);
    2183             :     }
    2184       41184 :     for (   ; i<=dx; i++)
    2185             :     {
    2186       25502 :       GEN c = gmul(y_lead, gel(x,i));
    2187       25502 :       gel(x,i) = rem(c,T);
    2188             :     }
    2189       15682 :     x++; dx--;
    2190       21254 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2191       15682 :     if (dx < dy) break;
    2192        8191 :     if (gc_needed(av2,1))
    2193             :     {
    2194           0 :       GEN X = x-2;
    2195           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2196           0 :       X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
    2197           0 :       gerepileall(av2,2, &X, &z); x = X+2;
    2198             :     }
    2199             :   }
    2200       14036 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2201        7491 :   if (dx < 0)
    2202         182 :     x = pol_0(vx);
    2203             :   else
    2204             :   {
    2205        7309 :     lx = dx+3; x -= 2;
    2206        7309 :     x[0] = evaltyp(t_POL) | _evallg(lx);
    2207        7309 :     x[1] = evalsigne(1) | evalvarn(vx);
    2208        7309 :     x = RgX_recip_i(x);
    2209             :   }
    2210        7491 :   z = RgX_recip_i(z);
    2211        7491 :   r = x;
    2212        7491 :   if (p)
    2213             :   {
    2214        3339 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2215        3339 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2216             :   }
    2217        7491 :   *ptr = r; return gc_all(av, 2, &z, ptr);
    2218             : }
    2219             : GEN
    2220       11776 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2221       11776 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2222             : 
    2223             : GEN
    2224           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2225           0 : { return RgXQX_red(RgX_mul(x,y), T); }
    2226             : GEN
    2227   500496680 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
    2228             : GEN
    2229      141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
    2230             : GEN
    2231       26145 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
    2232             : GEN
    2233          35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
    2234             : GEN
    2235         133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
    2236             : 
    2237             : GEN
    2238           0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
    2239             : 
    2240             : GEN
    2241           0 : RgXQX_powers(GEN P, long n, GEN T)
    2242             : {
    2243           0 :   GEN v = cgetg(n+2, t_VEC);
    2244             :   long i;
    2245           0 :   gel(v, 1) = pol_1(varn(T));
    2246           0 :   if (n==0) return v;
    2247           0 :   gel(v, 2) = gcopy(P);
    2248           0 :   for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
    2249           0 :   return v;
    2250             : }
    2251             : 
    2252             : static GEN
    2253      605085 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
    2254             : static GEN
    2255           0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
    2256             : static GEN
    2257       66945 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2258             : static GEN
    2259       86459 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
    2260             : static GEN
    2261      272626 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2262             : static GEN
    2263     1006576 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2264             : static GEN
    2265      984696 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2266             : static GEN
    2267         840 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2268             : static GEN
    2269      705040 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
    2270             : 
    2271             : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
    2272             :               _mul, _sqr, _one, _zero };
    2273             : 
    2274             : GEN
    2275           0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
    2276             : {
    2277           0 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
    2278             : }
    2279             : 
    2280             : GEN
    2281      401355 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2282             : {
    2283      401355 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2284      401355 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
    2285             : }
    2286             : 
    2287             : /* mod X^n */
    2288             : struct modXn {
    2289             :   long v; /* varn(X) */
    2290             :   long n;
    2291             : } ;
    2292             : static GEN
    2293       11893 : _sqrXn(void *data, GEN x) {
    2294       11893 :   struct modXn *S = (struct modXn*)data;
    2295       11893 :   return RgXn_sqr(x, S->n);
    2296             : }
    2297             : static GEN
    2298        4528 : _mulXn(void *data, GEN x, GEN y) {
    2299        4528 :   struct modXn *S = (struct modXn*)data;
    2300        4528 :   return RgXn_mul(x,y, S->n);
    2301             : }
    2302             : static GEN
    2303        1939 : _oneXn(void *data) {
    2304        1939 :   struct modXn *S = (struct modXn*)data;
    2305        1939 :   return pol_1(S->v);
    2306             : }
    2307             : static GEN
    2308           0 : _zeroXn(void *data) {
    2309           0 :   struct modXn *S = (struct modXn*)data;
    2310           0 :   return pol_0(S->v);
    2311             : }
    2312             : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
    2313             :                                           _oneXn, _zeroXn };
    2314             : 
    2315             : GEN
    2316         483 : RgXn_powers(GEN x, long m, long n)
    2317             : {
    2318         483 :   long d = degpol(x);
    2319         483 :   int use_sqr = (d<<1) >= n;
    2320             :   struct modXn S;
    2321         483 :   S.v = varn(x); S.n = n;
    2322         483 :   return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
    2323             : }
    2324             : 
    2325             : GEN
    2326        2286 : RgXn_powu_i(GEN x, ulong m, long n)
    2327             : {
    2328             :   struct modXn S;
    2329             :   long v;
    2330        2286 :   if (n == 0) return x;
    2331        2286 :   v = RgX_valrem(x, &x);
    2332        2286 :   if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
    2333        2265 :   S.v = varn(x); S.n = n;
    2334        2265 :   x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2335        2265 :   if (v) x = RgX_shift_shallow(x, m * v);
    2336        2265 :   return x;
    2337             : }
    2338             : GEN
    2339           0 : RgXn_powu(GEN x, ulong m, long n)
    2340             : {
    2341             :   pari_sp av;
    2342           0 :   if (n == 0) return gcopy(x);
    2343           0 :   av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
    2344             : }
    2345             : 
    2346             : GEN
    2347         966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
    2348             : {
    2349             :   struct modXn S;
    2350         966 :   S.v = varn(gel(x,2)); S.n = n;
    2351         966 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
    2352             : }
    2353             : 
    2354             : GEN
    2355           0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
    2356             : {
    2357           0 :   int use_sqr = 2*degpol(x) >= n;
    2358             :   struct modXn S;
    2359           0 :   S.v = varn(x); S.n = n;
    2360           0 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2361             : }
    2362             : 
    2363             : /* Q(x) mod t^n, x in R[t], n >= 1 */
    2364             : GEN
    2365        5145 : RgXn_eval(GEN Q, GEN x, long n)
    2366             : {
    2367        5145 :   long d = degpol(x);
    2368             :   int use_sqr;
    2369             :   struct modXn S;
    2370        5145 :   if (d == 1 && isrationalzero(gel(x,2)))
    2371             :   {
    2372        5138 :     GEN y = RgX_unscale(Q, gel(x,3));
    2373        5138 :     setvarn(y, varn(x)); return y;
    2374             :   }
    2375           7 :   S.v = varn(x);
    2376           7 :   S.n = n;
    2377           7 :   use_sqr = (d<<1) >= n;
    2378           7 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2379             : }
    2380             : 
    2381             : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
    2382             : static GEN
    2383     2613833 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2384             : {
    2385     2613833 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2386     2613833 :   return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
    2387             : }
    2388             : 
    2389             : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
    2390             : static GEN
    2391          14 : RgXn_sqrhigh(GEN f, long n2, long n)
    2392             : {
    2393          14 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2394          14 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2395             : }
    2396             : 
    2397             : static GEN
    2398     2188688 : RgXn_div_gen(GEN g, GEN f, long e)
    2399             : {
    2400             :   pari_sp av;
    2401             :   ulong mask;
    2402             :   GEN W, a;
    2403     2188688 :   long v = varn(f), n = 1;
    2404             : 
    2405     2188688 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2406     2188681 :   a = ginv(gel(f,2));
    2407     2188677 :   if (e == 1 && !g) return scalarpol(a, v);
    2408     2183805 :   else if (e == 2 && !g)
    2409             :   {
    2410             :     GEN b;
    2411      836926 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2412      246246 :     b = gneg(b);
    2413      246241 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2414      246243 :     return deg1pol(b, a, v);
    2415             :   }
    2416     1346879 :   av = avma;
    2417     1346879 :   W = scalarpol_shallow(a,v);
    2418     1346879 :   mask = quadratic_prec_mask(e);
    2419     3954818 :   while (mask > 1)
    2420             :   {
    2421             :     GEN u, fr;
    2422     2607939 :     long n2 = n;
    2423     2607939 :     n<<=1; if (mask & 1) n--;
    2424     2607939 :     mask >>= 1;
    2425     2607939 :     fr = RgXn_red_shallow(f, n);
    2426     2607939 :     if (mask>1 || !g)
    2427             :     {
    2428     1324176 :       u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2429     1324176 :       W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2430             :     }
    2431             :     else
    2432             :     {
    2433     1283763 :       GEN y = RgXn_mul(g, W, n), yt =  RgXn_red_shallow(y, n-n2);
    2434     1283763 :       u = RgXn_mul(yt, RgXn_mulhigh(fr,  W, n2, n), n-n2);
    2435     1283763 :       W = RgX_sub(y, RgX_shift_shallow(u, n2));
    2436             :     }
    2437     2607939 :     if (gc_needed(av,2))
    2438             :     {
    2439           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2440           0 :       W = gerepileupto(av, W);
    2441             :     }
    2442             :   }
    2443     1346879 :   return W;
    2444             : }
    2445             : 
    2446             : static GEN
    2447         119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
    2448             : {
    2449             :   GEN r;
    2450         119 :   if (lgefint(p) == 3)
    2451             :   {
    2452         119 :     ulong pp = uel(p, 2);
    2453         119 :     if (pp == 2)
    2454          14 :       r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
    2455             :     else
    2456         105 :       r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
    2457             :   }
    2458             :   else
    2459           0 :     r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
    2460         119 :   return FpX_to_mod(r, p);
    2461             : }
    2462             : 
    2463             : static GEN
    2464           7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
    2465             : {
    2466           7 :   GEN r, T = RgX_to_FpX(pol, p);
    2467           7 :   if (signe(T) == 0) pari_err_OP("/", x, y);
    2468           7 :   r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
    2469           7 :   return FpXQX_to_mod(r, T, p);
    2470             : }
    2471             : 
    2472             : static GEN
    2473          91 : RgXn_inv_FpX(GEN x, long e, GEN p)
    2474             : {
    2475             :   GEN r;
    2476          91 :   if (lgefint(p) == 3)
    2477             :   {
    2478          91 :     ulong pp = uel(p, 2);
    2479          91 :     if (pp == 2)
    2480          28 :       r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
    2481             :     else
    2482          63 :       r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
    2483             :   }
    2484             :   else
    2485           0 :     r = FpXn_inv(RgX_to_FpX(x, p), e, p);
    2486          91 :   return FpX_to_mod(r, p);
    2487             : }
    2488             : 
    2489             : static GEN
    2490           0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
    2491             : {
    2492           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2493           0 :   if (signe(T) == 0) pari_err_OP("/", gen_1, x);
    2494           0 :   r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
    2495           0 :   return FpXQX_to_mod(r, T, p);
    2496             : }
    2497             : 
    2498             : #define code(t1,t2) ((t1 << 6) | t2)
    2499             : 
    2500             : static GEN
    2501      905009 : RgXn_inv_fast(GEN x, long e)
    2502             : {
    2503             :   GEN p, pol;
    2504             :   long pa;
    2505      905009 :   long t = RgX_type(x,&p,&pol,&pa);
    2506      905016 :   switch(t)
    2507             :   {
    2508          91 :     case t_INTMOD: return RgXn_inv_FpX(x, e, p);
    2509           0 :     case code(t_POLMOD, t_INTMOD):
    2510           0 :                    return RgXn_inv_FpXQX(x, e, pol, p);
    2511      904925 :     default:       return NULL;
    2512             :   }
    2513             : }
    2514             : 
    2515             : static GEN
    2516     1283889 : RgXn_div_fast(GEN x, GEN y, long e)
    2517             : {
    2518             :   GEN p, pol;
    2519             :   long pa;
    2520     1283889 :   long t = RgX_type2(x,y,&p,&pol,&pa);
    2521     1283889 :   switch(t)
    2522             :   {
    2523         119 :     case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
    2524           7 :     case code(t_POLMOD, t_INTMOD):
    2525           7 :                    return RgXn_div_FpXQX(x, y, e, pol, p);
    2526     1283763 :     default:       return NULL;
    2527             :   }
    2528             : }
    2529             : #undef code
    2530             : 
    2531             : GEN
    2532     1283889 : RgXn_div_i(GEN g, GEN f, long e)
    2533             : {
    2534     1283889 :   GEN h = RgXn_div_fast(g, f, e);
    2535     1283889 :   if (h) return h;
    2536     1283763 :   return RgXn_div_gen(g, f, e);
    2537             : }
    2538             : 
    2539             : GEN
    2540      584436 : RgXn_div(GEN g, GEN f, long e)
    2541             : {
    2542      584436 :   pari_sp av = avma;
    2543      584436 :   return gerepileupto(av, RgXn_div_i(g, f, e));
    2544             : }
    2545             : 
    2546             : GEN
    2547      905010 : RgXn_inv_i(GEN f, long e)
    2548             : {
    2549      905010 :   GEN h = RgXn_inv_fast(f, e);
    2550      905016 :   if (h) return h;
    2551      904925 :   return RgXn_div_gen(NULL, f, e);
    2552             : }
    2553             : 
    2554             : GEN
    2555      807450 : RgXn_inv(GEN f, long e)
    2556             : {
    2557      807450 :   pari_sp av = avma;
    2558      807450 :   return gerepileupto(av, RgXn_inv_i(f, e));
    2559             : }
    2560             : 
    2561             : /* Compute intformal(x^n*S)/x^(n+1) */
    2562             : static GEN
    2563       58756 : RgX_integXn(GEN x, long n)
    2564             : {
    2565       58756 :   long i, lx = lg(x);
    2566             :   GEN y;
    2567       58756 :   if (lx == 2) return RgX_copy(x);
    2568       52722 :   y = cgetg(lx, t_POL); y[1] = x[1];
    2569      110573 :   for (i=2; i<lx; i++)
    2570       57853 :     gel(y,i) = gdivgs(gel(x,i), n+i-1);
    2571       52720 :   return RgX_renormalize_lg(y, lx);;
    2572             : }
    2573             : 
    2574             : GEN
    2575       52868 : RgXn_expint(GEN h, long e)
    2576             : {
    2577       52868 :   pari_sp av = avma, av2;
    2578       52868 :   long v = varn(h), n;
    2579       52868 :   GEN f = pol_1(v), g;
    2580             :   ulong mask;
    2581             : 
    2582       52869 :   if (!signe(h)) return f;
    2583       52869 :   g = pol_1(v);
    2584       52869 :   n = 1; mask = quadratic_prec_mask(e);
    2585       52869 :   av2 = avma;
    2586       58756 :   for (;mask>1;)
    2587             :   {
    2588             :     GEN u, w;
    2589       58756 :     long n2 = n;
    2590       58756 :     n<<=1; if (mask & 1) n--;
    2591       58756 :     mask >>= 1;
    2592       58756 :     u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
    2593       58757 :     u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
    2594       58756 :     w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
    2595       58755 :     f = RgX_add(f, RgX_shift_shallow(w, n2));
    2596       58753 :     if (mask<=1) break;
    2597        5887 :     u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
    2598        5887 :     g = RgX_sub(g, RgX_shift_shallow(u, n2));
    2599        5887 :     if (gc_needed(av2,2))
    2600             :     {
    2601           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
    2602           0 :       gerepileall(av2, 2, &f, &g);
    2603             :     }
    2604             :   }
    2605       52866 :   return gerepileupto(av, f);
    2606             : }
    2607             : 
    2608             : GEN
    2609           0 : RgXn_exp(GEN h, long e)
    2610             : {
    2611           0 :   long d = degpol(h);
    2612           0 :   if (d < 0) return pol_1(varn(h));
    2613           0 :   if (!d || !gequal0(gel(h,2)))
    2614           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2615           0 :   return RgXn_expint(RgX_deriv(h), e);
    2616             : }
    2617             : 
    2618             : GEN
    2619         154 : RgXn_reverse(GEN f, long e)
    2620             : {
    2621         154 :   pari_sp av = avma, av2;
    2622             :   ulong mask;
    2623             :   GEN fi, a, df, W, an;
    2624         154 :   long v = varn(f), n=1;
    2625         154 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2626           0 :     pari_err_INV("serreverse",f);
    2627         154 :   fi = ginv(gel(f,3));
    2628         154 :   a = deg1pol_shallow(fi,gen_0,v);
    2629         154 :   if (e <= 2) return gerepilecopy(av, a);
    2630         133 :   W = scalarpol(fi,v);
    2631         133 :   df = RgX_deriv(f);
    2632         133 :   mask = quadratic_prec_mask(e);
    2633         133 :   av2 = avma;
    2634         616 :   for (;mask>1;)
    2635             :   {
    2636             :     GEN u, fa, fr;
    2637         483 :     long n2 = n, rt;
    2638         483 :     n<<=1; if (mask & 1) n--;
    2639         483 :     mask >>= 1;
    2640         483 :     fr = RgXn_red_shallow(f, n);
    2641         483 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2642         483 :     an = RgXn_powers(a, rt, n);
    2643         483 :     if (n>1)
    2644             :     {
    2645         483 :       long n4 = (n2+1)>>1;
    2646         483 :       GEN dfr = RgXn_red_shallow(df, n2);
    2647         483 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2648         483 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2649         483 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2650             :     }
    2651         483 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2652         483 :     fa = RgX_shift(fa, -n2);
    2653         483 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2654         483 :     if (gc_needed(av2,2))
    2655             :     {
    2656           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2657           0 :       gerepileall(av2, 2, &a, &W);
    2658             :     }
    2659             :   }
    2660         133 :   return gerepileupto(av, a);
    2661             : }
    2662             : 
    2663             : GEN
    2664           7 : RgXn_sqrt(GEN h, long e)
    2665             : {
    2666           7 :   pari_sp av = avma, av2;
    2667           7 :   long v = varn(h), n = 1;
    2668           7 :   GEN f = scalarpol(gen_1, v), df = f;
    2669           7 :   ulong mask = quadratic_prec_mask(e);
    2670           7 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2671           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2672           7 :   av2 = avma;
    2673             :   while(1)
    2674           7 :   {
    2675          14 :     long n2 = n, m;
    2676             :     GEN g;
    2677          14 :     n<<=1; if (mask & 1) n--;
    2678          14 :     mask >>= 1;
    2679          14 :     m = n-n2;
    2680          14 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2681          14 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2682          14 :     if (mask==1) return gerepileupto(av, f);
    2683           7 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2684           7 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2685           7 :     if (gc_needed(av2,2))
    2686             :     {
    2687           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2688           0 :       gerepileall(av2, 2, &f, &df);
    2689             :     }
    2690             :   }
    2691             : }
    2692             : 
    2693             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2694             : GEN
    2695      116395 : RgXQ_powu(GEN x, ulong n, GEN T)
    2696             : {
    2697      116395 :   pari_sp av = avma;
    2698             : 
    2699      116395 :   if (!n) return pol_1(varn(x));
    2700       68211 :   if (n == 1) return RgX_copy(x);
    2701       19670 :   x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
    2702       19670 :   return gerepilecopy(av, x);
    2703             : }
    2704             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2705             : GEN
    2706      100438 : RgXQ_pow(GEN x, GEN n, GEN T)
    2707             : {
    2708             :   pari_sp av;
    2709      100438 :   long s = signe(n);
    2710             : 
    2711      100438 :   if (!s) return pol_1(varn(x));
    2712      100438 :   if (is_pm1(n) == 1)
    2713       86459 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2714       13979 :   av = avma;
    2715       13979 :   if (s < 0) x = RgXQ_inv(x, T);
    2716       13979 :   x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
    2717       13979 :   return gerepilecopy(av, x);
    2718             : }
    2719             : static GEN
    2720      213749 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
    2721             : static GEN
    2722      113877 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
    2723             : 
    2724             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2725             : GEN
    2726       10941 : ZXQ_powers(GEN x, long l, GEN T)
    2727             : {
    2728       10941 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2729       10941 :   return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
    2730             : }
    2731             : 
    2732             : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
    2733             : GEN
    2734      188556 : ZXQ_powu(GEN x, ulong n, GEN T)
    2735             : {
    2736      188556 :   pari_sp av = avma;
    2737             : 
    2738      188556 :   if (!n) return pol_1(varn(x));
    2739      188556 :   if (n == 1) return ZX_copy(x);
    2740      125044 :   x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
    2741      125048 :   return gerepilecopy(av, x);
    2742             : }
    2743             : 
    2744             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2745             : GEN
    2746        8498 : RgXQ_powers(GEN x, long l, GEN T)
    2747             : {
    2748        8498 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2749        8498 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2750             : }
    2751             : 
    2752             : GEN
    2753        7013 : RgXQV_factorback(GEN L, GEN e, GEN T)
    2754             : {
    2755        7013 :   return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
    2756             : }
    2757             : 
    2758             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2759             : GEN
    2760        8379 : QXQ_powers(GEN a, long n, GEN T)
    2761             : {
    2762             :   GEN den, v;
    2763        8379 :   if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
    2764        8358 :   v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
    2765             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2766        8358 :   if (den)
    2767             :   { /* restore denominators */
    2768        5257 :     GEN d = den;
    2769             :     long i;
    2770        5257 :     gel(v,2) = a;
    2771       26607 :     for (i=3; i<=n+1; i++) {
    2772       21350 :       d = mulii(d,den);
    2773       21350 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2774             :     }
    2775             :   }
    2776        8358 :   return v;
    2777             : }
    2778             : 
    2779             : static GEN
    2780        2527 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2781             : {
    2782        2527 :   long l, i, m = 0;
    2783             :   GEN dz, z;
    2784        2527 :   GEN V = cgetg_copy(v, &l);
    2785        8876 :   for (i = imin; i < l; i++)
    2786             :   {
    2787        6349 :     GEN c = gel(v, i);
    2788        6349 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2789             :   }
    2790        2527 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2791        2758 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2792        8876 :   for (i = imin; i < l; i++)
    2793             :   {
    2794        6349 :     GEN c = gel(v,i);
    2795        6349 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2796        6349 :     gel(V,i) = c;
    2797             :   }
    2798        2527 :   return V;
    2799             : }
    2800             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2801             : GEN
    2802        2296 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2803        2296 : { return do_QXQ_eval(v, 1, a, T); }
    2804             : 
    2805             : GEN
    2806         231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
    2807         231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2808             : 
    2809             : GEN
    2810        2891 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2811             : {
    2812        2891 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2813             : }
    2814             : 
    2815             : GEN
    2816        5144 : RgXQ_norm(GEN x, GEN T)
    2817             : {
    2818             :   pari_sp av;
    2819        5144 :   long dx = degpol(x);
    2820             :   GEN L, y;
    2821        5144 :   if (degpol(T)==0) return gpowgs(x,0);
    2822        5137 :   av = avma; y = resultant(T, x);
    2823        5137 :   L = leading_coeff(T);
    2824        5137 :   if (gequal1(L) || !signe(x)) return y;
    2825           0 :   return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
    2826             : }
    2827             : 
    2828             : GEN
    2829         448 : RgXQ_trace(GEN x, GEN T)
    2830             : {
    2831         448 :   pari_sp av = avma;
    2832             :   GEN dT, z;
    2833             :   long n;
    2834         448 :   if (degpol(T)==0) return gmulgs(x,0);
    2835         441 :   dT = RgX_deriv(T); n = degpol(dT);
    2836         441 :   z = RgXQ_mul(x, dT, T);
    2837         441 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2838         392 :   return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
    2839             : }
    2840             : 
    2841             : GEN
    2842     2786958 : RgX_blocks(GEN P, long n, long m)
    2843             : {
    2844     2786958 :   GEN z = cgetg(m+1,t_VEC);
    2845     2786958 :   long i,j, k=2, l = lg(P);
    2846     8554452 :   for(i=1; i<=m; i++)
    2847             :   {
    2848     5767494 :     GEN zi = cgetg(n+2,t_POL);
    2849     5767494 :     zi[1] = P[1];
    2850     5767494 :     gel(z,i) = zi;
    2851    16985676 :     for(j=2; j<n+2; j++)
    2852    11218182 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2853     5767494 :     zi = RgX_renormalize_lg(zi, n+2);
    2854             :   }
    2855     2786958 :   return z;
    2856             : }
    2857             : 
    2858             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2859             : void
    2860    49480010 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2861             : {
    2862    49480010 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2863             :   GEN p0, p1;
    2864             : 
    2865    49481028 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2866             : 
    2867    49387975 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2868    49387975 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2869    49392858 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2870   145966497 :   for (i=0; i<n1; i++)
    2871             :   {
    2872    96573418 :     p0[2+i] = p[2+(i<<1)];
    2873    96573418 :     p1[2+i] = p[3+(i<<1)];
    2874             :   }
    2875    49393079 :   if (n1 != n0)
    2876    17864494 :     p0[2+i] = p[2+(i<<1)];
    2877    49393079 :   *pe = normalizepol(p0);
    2878    49395897 :   *po = normalizepol(p1);
    2879             : }
    2880             : 
    2881             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2882             : GEN
    2883       43547 : RgX_splitting(GEN p, long k)
    2884             : {
    2885       43547 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2886             :   GEN r;
    2887             : 
    2888       43547 :   m = n/k;
    2889       43547 :   r = cgetg(k+1,t_VEC);
    2890      233807 :   for(i=1; i<=k; i++)
    2891             :   {
    2892      190260 :     gel(r,i) = cgetg(m+3, t_POL);
    2893      190260 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2894             :   }
    2895      571487 :   for (j=1, i=0, l=2; i<=n; i++)
    2896             :   {
    2897      527940 :     gmael(r,j,l) = gel(p,2+i);
    2898      527940 :     if (j==k) { j=1; l++; } else j++;
    2899             :   }
    2900      233807 :   for(i=1; i<=k; i++)
    2901      190260 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2902       43547 :   return r;
    2903             : }
    2904             : 
    2905             : /*******************************************************************/
    2906             : /*                                                                 */
    2907             : /*                        Kronecker form                           */
    2908             : /*                                                                 */
    2909             : /*******************************************************************/
    2910             : 
    2911             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2912             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2913             :  * normalized coefficients */
    2914             : GEN
    2915      179217 : Kronecker_to_mod(GEN z, GEN T)
    2916             : {
    2917      179217 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2918      179217 :   GEN x, t = cgetg(N,t_POL);
    2919      179217 :   t[1] = T[1];
    2920      179217 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2921      179217 :   x[1] = z[1];
    2922      179217 :   T = RgX_copy(T);
    2923     1321757 :   for (i=2; i<lx+2; i++, z+= N-2)
    2924             :   {
    2925     5487746 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2926     1142541 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2927             :   }
    2928      179216 :   N = (l-2) % (N-2) + 2;
    2929      459603 :   for (j=2; j<N; j++) t[j] = z[j];
    2930      179216 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2931      179217 :   return normalizepol_lg(x, i+1);
    2932             : }
    2933             : 
    2934             : /*******************************************************************/
    2935             : /*                                                                 */
    2936             : /*                        Domain detection                         */
    2937             : /*                                                                 */
    2938             : /*******************************************************************/
    2939             : 
    2940             : static GEN
    2941      152790 : zero_FpX_mod(GEN p, long v)
    2942             : {
    2943      152790 :   GEN r = cgetg(3,t_POL);
    2944      152790 :   r[1] = evalvarn(v);
    2945      152790 :   gel(r,2) = mkintmod(gen_0, icopy(p));
    2946      152790 :   return r;
    2947             : }
    2948             : 
    2949             : static GEN
    2950      818173 : RgX_mul_FpX(GEN x, GEN y, GEN p)
    2951             : {
    2952      818173 :   pari_sp av = avma;
    2953             :   GEN r;
    2954      818173 :   if (lgefint(p) == 3)
    2955             :   {
    2956      766817 :     ulong pp = uel(p, 2);
    2957      766817 :     r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
    2958             :                                   RgX_to_Flx(y, pp), pp));
    2959             :   }
    2960             :   else
    2961       51356 :     r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2962      818173 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2963      703120 :   return gerepileupto(av, FpX_to_mod(r, p));
    2964             : }
    2965             : 
    2966             : static GEN
    2967         532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
    2968             : {
    2969         532 :   GEN r = cgetg(3,t_POL);
    2970         532 :   r[1] = evalvarn(v);
    2971         532 :   gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
    2972         532 :   return r;
    2973             : }
    2974             : 
    2975             : static GEN
    2976        2044 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2977             : {
    2978        2044 :   pari_sp av = avma;
    2979             :   long dT;
    2980             :   GEN kx, ky, r;
    2981        2044 :   GEN T = RgX_to_FpX(pol, p);
    2982        2044 :   if (signe(T)==0) pari_err_OP("*", x, y);
    2983        2037 :   dT = degpol(T);
    2984        2037 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    2985        2037 :   ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
    2986        2037 :   r = FpX_mul(kx, ky, p);
    2987        2037 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    2988        1512 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    2989             : }
    2990             : 
    2991             : static GEN
    2992      418146 : RgX_liftred(GEN x, GEN T)
    2993      418146 : { return RgXQX_red(liftpol_shallow(x), T); }
    2994             : 
    2995             : static GEN
    2996      164098 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
    2997             : {
    2998      164098 :   pari_sp av = avma;
    2999      164098 :   long dT = degpol(T);
    3000      164098 :   GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
    3001             :                  RgXX_to_Kronecker(RgX_liftred(y, T), dT));
    3002      164098 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    3003             : }
    3004             : 
    3005             : static GEN
    3006        1421 : RgX_sqr_FpX(GEN x, GEN p)
    3007             : {
    3008        1421 :   pari_sp av = avma;
    3009             :   GEN r;
    3010        1421 :   if (lgefint(p) == 3)
    3011             :   {
    3012        1217 :     ulong pp = uel(p, 2);
    3013        1217 :     r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
    3014             :   }
    3015             :   else
    3016         204 :     r = FpX_sqr(RgX_to_FpX(x, p), p);
    3017        1421 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3018        1281 :   return gerepileupto(av, FpX_to_mod(r, p));
    3019             : }
    3020             : 
    3021             : static GEN
    3022         196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
    3023             : {
    3024         196 :   pari_sp av = avma;
    3025             :   long dT;
    3026         196 :   GEN kx, r, T = RgX_to_FpX(pol, p);
    3027         196 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3028         189 :   dT = degpol(T);
    3029         189 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    3030         189 :   r = FpX_sqr(kx, p);
    3031         189 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3032         189 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    3033             : }
    3034             : 
    3035             : static GEN
    3036       13418 : RgX_sqr_QXQX(GEN x, GEN T)
    3037             : {
    3038       13418 :   pari_sp av = avma;
    3039       13418 :   long dT = degpol(T);
    3040       13418 :   GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
    3041       13418 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    3042             : }
    3043             : 
    3044             : static GEN
    3045       68670 : RgX_rem_FpX(GEN x, GEN y, GEN p)
    3046             : {
    3047       68670 :   pari_sp av = avma;
    3048             :   GEN r;
    3049       68670 :   if (lgefint(p) == 3)
    3050             :   {
    3051       51880 :     ulong pp = uel(p, 2);
    3052       51880 :     r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
    3053             :                                   RgX_to_Flx(y, pp), pp));
    3054             :   }
    3055             :   else
    3056       16790 :     r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3057       68670 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3058       31073 :   return gerepileupto(av, FpX_to_mod(r, p));
    3059             : }
    3060             : 
    3061             : static GEN
    3062       38266 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
    3063             : {
    3064       38266 :   pari_sp av = avma;
    3065             :   GEN r;
    3066       38266 :   r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3067       38266 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3068             : }
    3069             : static GEN
    3070          70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3071             : {
    3072          70 :   pari_sp av = avma;
    3073             :   GEN r;
    3074          70 :   GEN T = RgX_to_FpX(pol, p);
    3075          70 :   if (signe(T) == 0) pari_err_OP("%", x, y);
    3076          63 :   if (lgefint(p) == 3)
    3077             :   {
    3078          55 :     ulong pp = uel(p, 2);
    3079          55 :     GEN Tp = ZX_to_Flx(T, pp);
    3080          55 :     r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
    3081             :                               RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3082             :   }
    3083             :   else
    3084           8 :     r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3085          63 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3086          56 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3087             : }
    3088             : 
    3089             : #define code(t1,t2) ((t1 << 6) | t2)
    3090             : static GEN
    3091    74933099 : RgX_mul_fast(GEN x, GEN y)
    3092             : {
    3093             :   GEN p, pol;
    3094             :   long pa;
    3095    74933099 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3096    74933371 :   switch(t)
    3097             :   {
    3098    38791206 :     case t_INT:    return ZX_mul(x,y);
    3099     3303045 :     case t_FRAC:   return QX_mul(x,y);
    3100      104735 :     case t_FFELT:  return FFX_mul(x, y, pol);
    3101      818173 :     case t_INTMOD: return RgX_mul_FpX(x, y, p);
    3102      164098 :     case code(t_POLMOD, t_INT):
    3103             :     case code(t_POLMOD, t_FRAC):
    3104      164098 :                    return RgX_mul_QXQX(x, y, pol);
    3105        2018 :     case code(t_POLMOD, t_INTMOD):
    3106        2018 :                    return RgX_mul_FpXQX(x, y, pol, p);
    3107    31750096 :     default:       return NULL;
    3108             :   }
    3109             : }
    3110             : static GEN
    3111     1218160 : RgX_sqr_fast(GEN x)
    3112             : {
    3113             :   GEN p, pol;
    3114             :   long pa;
    3115     1218160 :   long t = RgX_type(x,&p,&pol,&pa);
    3116     1218160 :   switch(t)
    3117             :   {
    3118     1161005 :     case t_INT:    return ZX_sqr(x);
    3119       30731 :     case t_FRAC:   return QX_sqr(x);
    3120        2534 :     case t_FFELT:  return FFX_sqr(x, pol);
    3121        1421 :     case t_INTMOD: return RgX_sqr_FpX(x, p);
    3122       13418 :     case code(t_POLMOD, t_INT):
    3123             :     case code(t_POLMOD, t_FRAC):
    3124       13418 :                    return RgX_sqr_QXQX(x, pol);
    3125         196 :     case code(t_POLMOD, t_INTMOD):
    3126         196 :                    return RgX_sqr_FpXQX(x, pol, p);
    3127        8855 :     default:       return NULL;
    3128             :   }
    3129             : }
    3130             : 
    3131             : static GEN
    3132    14205237 : RgX_rem_fast(GEN x, GEN y)
    3133             : {
    3134             :   GEN p, pol;
    3135             :   long pa;
    3136    14205237 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3137    14205384 :   switch(t)
    3138             :   {
    3139     3256722 :     case t_INT:    return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
    3140     1690954 :     case t_FRAC:   return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
    3141          84 :     case t_FFELT:  return FFX_rem(x, y, pol);
    3142       68670 :     case t_INTMOD: return RgX_rem_FpX(x, y, p);
    3143       38266 :     case code(t_POLMOD, t_INT):
    3144             :     case code(t_POLMOD, t_FRAC):
    3145       38266 :                    return RgX_rem_QXQX(x, y, pol);
    3146          60 :     case code(t_POLMOD, t_INTMOD):
    3147          60 :                    return RgX_rem_FpXQX(x, y, pol, p);
    3148     9150629 :     default:       return NULL;
    3149             :   }
    3150             : }
    3151             : 
    3152             : #undef code
    3153             : 
    3154             : GEN
    3155    61196127 : RgX_mul(GEN x, GEN y)
    3156             : {
    3157    61196127 :   GEN z = RgX_mul_fast(x,y);
    3158    61196287 :   if (!z) z = RgX_mul_i(x,y);
    3159    61195914 :   return z;
    3160             : }
    3161             : 
    3162             : GEN
    3163     1206008 : RgX_sqr(GEN x)
    3164             : {
    3165     1206008 :   GEN z = RgX_sqr_fast(x);
    3166     1206001 :   if (!z) z = RgX_sqr_i(x);
    3167     1206001 :   return z;
    3168             : }
    3169             : 
    3170             : GEN
    3171    14205250 : RgX_rem(GEN x, GEN y)
    3172             : {
    3173    14205250 :   GEN z = RgX_rem_fast(x, y);
    3174    14205358 :   if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
    3175    14205321 :   return z;
    3176             : }
    3177             : 
    3178             : static GEN
    3179           0 : _RgX_mul(void* E, GEN x, GEN y)
    3180           0 : { (void) E; return RgX_mul(x, y); }
    3181             : 
    3182             : GEN
    3183           0 : RgXV_prod(GEN V)
    3184           0 : { return gen_product(V, NULL, &_RgX_mul); }
    3185             : 
    3186             : GEN
    3187    11064369 : RgXn_mul(GEN f, GEN g, long n)
    3188             : {
    3189    11064369 :   pari_sp av = avma;
    3190    11064369 :   GEN h = RgX_mul_fast(f,g);
    3191    11064374 :   if (!h) return RgXn_mul2(f,g,n);
    3192     1786642 :   if (degpol(h) < n) return h;
    3193      382080 :   return gerepilecopy(av, RgXn_red_shallow(h, n));
    3194             : }
    3195             : 
    3196             : GEN
    3197       12152 : RgXn_sqr(GEN f, long n)
    3198             : {
    3199       12152 :   pari_sp av = avma;
    3200       12152 :   GEN g = RgX_sqr_fast(f);
    3201       12152 :   if (!g) return RgXn_sqr2(f,n);
    3202       11837 :   if (degpol(g) < n) return g;
    3203       10640 :   return gerepilecopy(av, RgXn_red_shallow(g, n));
    3204             : }
    3205             : 
    3206             : /* (f*g) \/ x^n */
    3207             : GEN
    3208     2672602 : RgX_mulhigh_i(GEN f, GEN g, long n)
    3209             : {
    3210     2672602 :   GEN h = RgX_mul_fast(f,g);
    3211     2672602 :   return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
    3212             : }
    3213             : 
    3214             : /* (f*g) \/ x^n */
    3215             : GEN
    3216           0 : RgX_sqrhigh_i(GEN f, long n)
    3217             : {
    3218           0 :   GEN h = RgX_sqr_fast(f);
    3219           0 :   return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
    3220             : }

Generated by: LCOV version 1.16