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.14.0 lcov report (development 27775-aca467eab2) Lines: 1607 1768 90.9 %
Date: 2022-07-03 07:33:15 Functions: 198 217 91.2 %
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    19308667 : brent_kung_optpow(long d, long n, long m)
      32             : {
      33             :   long p, r;
      34    19308667 :   long pold=1, rold=n*(d-1);
      35    92514316 :   for(p=2; p<=d; p++)
      36             :   {
      37    73205649 :     r = m*(p-1) + n*((d-1)/p);
      38    73205649 :     if (r<rold) { pold=p; rold=r; }
      39             :   }
      40    19308667 :   return pold;
      41             : }
      42             : 
      43             : static GEN
      44    13054220 : 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    13054220 :   pari_sp av = avma;
      48             :   long i;
      49    13054220 :   GEN z = cmul(E,P,a,ff->one(E));
      50    13011562 :   if (!z) z = gen_0;
      51    67230285 :   for (i=1; i<=n; i++)
      52             :   {
      53    54211545 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      54    54214842 :     if (t) {
      55    41478654 :       z = ff->add(E, z, t);
      56    41465889 :       if (gc_needed(av,2)) z = gerepileupto(av, z);
      57             :     }
      58             :   }
      59    13018740 :   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    11614797 : 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    11614797 :   pari_sp av = avma;
      72    11614797 :   long l = lg(V)-1;
      73             :   GEN z, u;
      74             : 
      75    11614797 :   if (d < 0) return ff->zero(E);
      76    11051323 :   if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      77     1002437 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      78     1002437 :   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     1002437 :   d -= l;
      84     1002437 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      85     2001371 :   while (d >= l-1)
      86             :   {
      87      997619 :     d -= l-1;
      88      997619 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      89      997489 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      90      997569 :     if (gc_needed(av,2))
      91          70 :       z = gerepileupto(av, z);
      92             :   }
      93     1003752 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      94     1003781 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      95     1003771 :   return gerepileupto(av, ff->red(E,z));
      96             : }
      97             : 
      98             : GEN
      99      725542 : 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      725542 :   pari_sp av = avma;
     103             :   GEN z, V;
     104             :   long rtd;
     105      725542 :   if (d < 0) return ff->zero(E);
     106      724611 :   rtd = (long) sqrt((double)d);
     107      724611 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     108      724634 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     109      724641 :   return gerepileupto(av, z);
     110             : }
     111             : 
     112             : static GEN
     113     1187631 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     114             : static GEN
     115    15652276 : _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     1034136 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     120             : static GEN
     121      338360 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     122             : static GEN
     123     1211381 : _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      225316 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
     161             : {
     162      225316 :   pari_sp av = avma;
     163      225316 :   long i, d = degpol(P);
     164             :   GEN s;
     165      225316 :   if (signe(P)==0) return pol_0(varn(P));
     166      225316 :   s = gel(P, d+2);
     167      225316 :   if (d == 0) return gcopy(s);
     168     1076313 :   for (i = d-1; i >= 0; i--)
     169             :   {
     170      854343 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
     171      854343 :     if (gc_needed(av,1))
     172             :     {
     173           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
     174           0 :       s = gerepileupto(av, s);
     175             :     }
     176             :   }
     177      221970 :   return gerepileupto(av, s);
     178             : }
     179             : 
     180             : GEN
     181        1568 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
     182             : {
     183        1568 :   pari_sp av = avma;
     184        1568 :   long i, d = degpol(P), v = varn(A);
     185             :   GEN s;
     186        1568 :   if (signe(P)==0) return pol_0(v);
     187        1568 :   if (d == 0) return scalarpol(gel(P, d+2), v);
     188        1190 :   s = scalarpol_shallow(gel(P, d+2), v);
     189        4851 :   for (i = d-1; i >= 0; i--)
     190             :   {
     191        3661 :     GEN c = gel(P,i+2), b = gel(B,d+1-i);
     192        3661 :     s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
     193        3661 :     if (gc_needed(av,1))
     194             :     {
     195           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
     196           0 :       s = gerepileupto(av, s);
     197             :     }
     198             :   }
     199        1190 :   return gerepileupto(av, s);
     200             : }
     201             : 
     202             : const struct bb_algebra *
     203      153409 : get_Rg_algebra(void)
     204             : {
     205      153409 :   return &Rg_algebra;
     206             : }
     207             : 
     208             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     209             : 
     210             : static GEN
     211       18004 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     212             : {
     213             :   (void) E;
     214       18004 :   return RgX_divrem(x, y, r);
     215             : }
     216             : 
     217             : GEN
     218        4487 : RgX_digits(GEN x, GEN T)
     219             : {
     220        4487 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     221        4487 :   return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     222             : }
     223             : 
     224             : /*******************************************************************/
     225             : /*                                                                 */
     226             : /*                         RgX                                     */
     227             : /*                                                                 */
     228             : /*******************************************************************/
     229             : 
     230             : long
     231    22790495 : RgX_equal(GEN x, GEN y)
     232             : {
     233    22790495 :   long i = lg(x);
     234             : 
     235    22790495 :   if (i != lg(y)) return 0;
     236    96423093 :   for (i--; i > 1; i--)
     237    73936560 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     238    22486533 :   return 1;
     239             : }
     240             : 
     241             : /* Returns 1 in the base ring over which x is defined */
     242             : /* HACK: this also works for t_SER */
     243             : GEN
     244     2139127 : Rg_get_1(GEN x)
     245             : {
     246             :   GEN p, T;
     247     2139127 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     248     2139129 :   if (RgX_type_is_composite(tx))
     249       37726 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     250     2139129 :   switch(tx)
     251             :   {
     252         357 :     case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
     253           7 :     case t_PADIC: return cvtop(gen_1, p, lx);
     254         728 :     case t_FFELT: return FF_1(T);
     255     2138037 :     default: return gen_1;
     256             :   }
     257             : }
     258             : /* Returns 0 in the base ring over which x is defined */
     259             : /* HACK: this also works for t_SER */
     260             : GEN
     261     6973365 : Rg_get_0(GEN x)
     262             : {
     263             :   GEN p, T;
     264     6973365 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     265     6973365 :   if (RgX_type_is_composite(tx))
     266       51478 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     267     6973365 :   switch(tx)
     268             :   {
     269         434 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     270          42 :     case t_PADIC: return zeropadic(p, lx);
     271         231 :     case t_FFELT: return FF_zero(T);
     272     6972658 :     default: return gen_0;
     273             :   }
     274             : }
     275             : 
     276             : GEN
     277        4977 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     278             : {
     279        4977 :   long i, n = degpol(P);
     280             :   GEN z, dz, dP;
     281        4977 :   if (n < 0) return gen_0;
     282        4977 :   P = Q_remove_denom(P, &dP);
     283        4977 :   z = gel(P,2); if (n == 0) return icopy(z);
     284        2961 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     285        2961 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     286        5859 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     287        2961 :   dz = mul_denom(dP, dV);
     288        2961 :   return dz? RgX_Rg_div(z, dz): z;
     289             : }
     290             : 
     291             : /* Return P(h * x), not memory clean */
     292             : GEN
     293       17807 : RgX_unscale(GEN P, GEN h)
     294             : {
     295       17807 :   long i, l = lg(P);
     296       17807 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     297       17807 :   Q[1] = P[1];
     298       17807 :   if (l == 2) return Q;
     299       17779 :   gel(Q,2) = gcopy(gel(P,2));
     300       39483 :   for (i=3; i<l; i++)
     301             :   {
     302       21704 :     hi = gmul(hi,h);
     303       21704 :     gel(Q,i) = gmul(gel(P,i), hi);
     304             :   }
     305       17779 :   return Q;
     306             : }
     307             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     308             : GEN
     309     1452580 : ZX_z_unscale(GEN P, long h)
     310             : {
     311     1452580 :   long i, l = lg(P);
     312     1452580 :   GEN Q = cgetg(l, t_POL);
     313     1452580 :   Q[1] = P[1];
     314     1452580 :   if (l == 2) return Q;
     315     1364497 :   gel(Q,2) = gel(P,2);
     316     1364497 :   if (l == 3) return Q;
     317     1338079 :   if (h == -1)
     318      532122 :     for (i = 3; i < l; i++)
     319             :     {
     320      488417 :       gel(Q,i) = negi(gel(P,i));
     321      488417 :       if (++i == l) break;
     322      434907 :       gel(Q,i) = gel(P,i);
     323             :     }
     324             :   else
     325             :   {
     326             :     GEN hi;
     327     1240864 :     gel(Q,3) = mulis(gel(P,3), h);
     328     1240864 :     hi = sqrs(h);
     329     5342181 :     for (i = 4; i < l; i++)
     330             :     {
     331     4101316 :       gel(Q,i) = mulii(gel(P,i), hi);
     332     4101320 :       if (i != l-1) hi = mulis(hi,h);
     333             :     }
     334             :   }
     335     1338080 :   return Q;
     336             : }
     337             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     338             : GEN
     339      993698 : ZX_unscale(GEN P, GEN h)
     340             : {
     341             :   long i, l;
     342             :   GEN Q, hi;
     343      993698 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     344         734 :   l = lg(P); Q = cgetg(l, t_POL);
     345         734 :   Q[1] = P[1];
     346         734 :   if (l == 2) return Q;
     347         734 :   gel(Q,2) = gel(P,2);
     348         734 :   if (l == 3) return Q;
     349         734 :   hi = h;
     350         734 :   gel(Q,3) = mulii(gel(P,3), hi);
     351        2251 :   for (i = 4; i < l; i++)
     352             :   {
     353        1517 :     hi = mulii(hi,h);
     354        1517 :     gel(Q,i) = mulii(gel(P,i), hi);
     355             :   }
     356         734 :   return Q;
     357             : }
     358             : /* P a ZX. Return P(x << n), not memory clean */
     359             : GEN
     360      215980 : ZX_unscale2n(GEN P, long n)
     361             : {
     362      215980 :   long i, ni = n, l = lg(P);
     363      215980 :   GEN Q = cgetg(l, t_POL);
     364      215981 :   Q[1] = P[1];
     365      215981 :   if (l == 2) return Q;
     366      215981 :   gel(Q,2) = gel(P,2);
     367      215981 :   if (l == 3) return Q;
     368      215981 :   gel(Q,3) = shifti(gel(P,3), ni);
     369      796497 :   for (i=4; i<l; i++)
     370             :   {
     371      580552 :     ni += n;
     372      580552 :     gel(Q,i) = shifti(gel(P,i), ni);
     373             :   }
     374      215945 :   return Q;
     375             : }
     376             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     377             : GEN
     378        1281 : ZX_unscale_div(GEN P, GEN h)
     379             : {
     380        1281 :   long i, l = lg(P);
     381        1281 :   GEN hi, Q = cgetg(l, t_POL);
     382        1281 :   Q[1] = P[1];
     383        1281 :   if (l == 2) return Q;
     384        1281 :   gel(Q,2) = diviiexact(gel(P,2), h);
     385        1281 :   if (l == 3) return Q;
     386        1281 :   gel(Q,3) = gel(P,3);
     387        1281 :   if (l == 4) return Q;
     388        1281 :   hi = h;
     389        1281 :   gel(Q,4) = mulii(gel(P,4), hi);
     390        3787 :   for (i=5; i<l; i++)
     391             :   {
     392        2506 :     hi = mulii(hi,h);
     393        2506 :     gel(Q,i) = mulii(gel(P,i), hi);
     394             :   }
     395        1281 :   return Q;
     396             : }
     397             : /* P(h*X) / h^k, assuming the result is a ZX */
     398             : GEN
     399         966 : ZX_unscale_divpow(GEN P, GEN h, long k)
     400             : {
     401         966 :   long i, j, l = lg(P);
     402         966 :   GEN H, Q = cgetg(l, t_POL);
     403         966 :   Q[1] = P[1]; if (l == 2) return Q;
     404         966 :   H = gpowers(h, maxss(k, l - 3 - k));
     405        3864 :   for (i = 2, j = k+1; j > 1 && i < l; i++)
     406        2898 :     gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
     407         966 :   if (i == l) return Q;
     408         966 :   gel(Q, i) = gel(P, i); i++;
     409        3542 :   for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
     410         966 :   return Q;
     411             : }
     412             : 
     413             : GEN
     414        5474 : RgXV_unscale(GEN x, GEN h)
     415             : {
     416        5474 :   if (!h || isint1(h)) return x;
     417       15931 :   pari_APPLY_same(RgX_unscale(gel(x,i), h));
     418             : }
     419             : 
     420             : /* Return h^degpol(P) P(x / h), not memory clean */
     421             : GEN
     422     2418473 : RgX_rescale(GEN P, GEN h)
     423             : {
     424     2418473 :   long i, l = lg(P);
     425     2418473 :   GEN Q = cgetg(l,t_POL), hi = h;
     426     2418445 :   gel(Q,l-1) = gel(P,l-1);
     427     7581013 :   for (i=l-2; i>=2; i--)
     428             :   {
     429     7578327 :     gel(Q,i) = gmul(gel(P,i), hi);
     430     7577904 :     if (i == 2) break;
     431     5161928 :     hi = gmul(hi,h);
     432             :   }
     433     2418662 :   Q[1] = P[1]; return Q;
     434             : }
     435             : 
     436             : /* A(X^d) --> A(X) */
     437             : GEN
     438      717228 : RgX_deflate(GEN x0, long d)
     439             : {
     440             :   GEN z, y, x;
     441      717228 :   long i,id, dy, dx = degpol(x0);
     442      717227 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     443      330820 :   dy = dx/d;
     444      330820 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     445      330816 :   z = y + 2;
     446      330816 :   x = x0+ 2;
     447     1365614 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     448      330816 :   return y;
     449             : }
     450             : 
     451             : GEN
     452        3584 : RgX_homogenize(GEN P, long v)
     453             : {
     454             :   long i, l, d;
     455        3584 :   GEN Q = cgetg_copy(P, &l);
     456        3584 :   Q[1] = P[1]; d = l-3;
     457       31794 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
     458        3584 :   return Q;
     459             : }
     460             : 
     461             : /* F a t_RFRAC */
     462             : long
     463          77 : rfrac_deflate_order(GEN F)
     464             : {
     465          77 :   GEN N = gel(F,1), D = gel(F,2);
     466          77 :   long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
     467          77 :   if (m == 1) return 1;
     468          35 :   if (typ(N) == t_POL && varn(N) == varn(D))
     469          21 :     m = cgcd(m, RgX_deflate_order(N));
     470          35 :   return m;
     471             : }
     472             : /* F a t_RFRAC */
     473             : GEN
     474          77 : rfrac_deflate_max(GEN F, long *m)
     475             : {
     476          77 :   *m = rfrac_deflate_order(F);
     477          77 :   return rfrac_deflate(F, *m);
     478             : }
     479             : /* F a t_RFRAC */
     480             : GEN
     481          77 : rfrac_deflate(GEN F, long m)
     482             : {
     483          77 :   GEN N = gel(F,1), D = gel(F,2);
     484          77 :   if (m == 1) return F;
     485          35 :   if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
     486          35 :   D = RgX_deflate(D, m); return mkrfrac(N, D);
     487             : }
     488             : 
     489             : /* return x0(X^d) */
     490             : GEN
     491      663144 : RgX_inflate(GEN x0, long d)
     492             : {
     493      663144 :   long i, id, dy, dx = degpol(x0);
     494      663144 :   GEN x = x0 + 2, z, y;
     495      663144 :   if (dx <= 0) return leafcopy(x0);
     496      655542 :   dy = dx*d;
     497      655542 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     498      655542 :   z = y + 2;
     499    23458258 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     500    10789148 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     501      655542 :   return y;
     502             : }
     503             : 
     504             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     505             : static GEN
     506     4348490 : RgX_translate_basecase(GEN P, GEN c)
     507             : {
     508     4348490 :   pari_sp av = avma;
     509             :   GEN Q, R;
     510             :   long i, k, n;
     511             : 
     512     4348490 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     513     4346472 :   Q = leafcopy(P);
     514     4346478 :   R = Q+2; n = degpol(P);
     515     4346478 :   if (isint1(c))
     516             :   {
     517       26807 :     for (i=1; i<=n; i++)
     518             :     {
     519      138971 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
     520       22909 :       if (gc_needed(av,2))
     521             :       {
     522           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
     523           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     524             :       }
     525             :     }
     526             :   }
     527     4342570 :   else if (isintm1(c))
     528             :   {
     529      913884 :     for (i=1; i<=n; i++)
     530             :     {
     531     2909242 :       for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
     532      693821 :       if (gc_needed(av,2))
     533             :       {
     534           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
     535           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     536             :       }
     537             :     }
     538             :   }
     539             :   else
     540             :   {
     541    14687658 :     for (i=1; i<=n; i++)
     542             :     {
     543    37111475 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
     544    10565152 :       if (gc_needed(av,2))
     545             :       {
     546           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
     547           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     548             :       }
     549             :     }
     550             :   }
     551     4345779 :   return gerepilecopy(av, Q);
     552             : }
     553             : GEN
     554     4348914 : RgX_translate(GEN P, GEN c)
     555             : {
     556     4348914 :   pari_sp av = avma;
     557     4348914 :   long n = degpol(P);
     558     4348907 :   if (n < 40)
     559     4348490 :     return RgX_translate_basecase(P, c);
     560             :   else
     561             :   {
     562         417 :     long d = n >> 1;
     563         417 :     GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
     564         418 :     GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
     565         418 :     GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
     566         418 :     return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
     567             :   }
     568             : }
     569             : 
     570             : /* P(ax + b) */
     571             : GEN
     572           0 : RgX_affine(GEN P, GEN a, GEN b)
     573             : {
     574           0 :   if (signe(b)) P = RgX_translate(P, b);
     575           0 :   return RgX_unscale(P, a);
     576             : }
     577             : 
     578             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     579             : GEN
     580       22407 : RgXQX_translate(GEN P, GEN c, GEN T)
     581             : {
     582       22407 :   pari_sp av = avma;
     583             :   GEN Q, R;
     584             :   long i, k, n;
     585             : 
     586       22407 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     587       22099 :   Q = leafcopy(P);
     588       22099 :   R = Q+2; n = degpol(P);
     589       78617 :   for (i=1; i<=n; i++)
     590             :   {
     591      252560 :     for (k=n-i; k<n; k++)
     592             :     {
     593      196042 :       pari_sp av2 = avma;
     594      196042 :       gel(R,k) = gerepileupto(av2,
     595      196042 :                    RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
     596             :     }
     597       56518 :     if (gc_needed(av,2))
     598             :     {
     599           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     600           0 :       Q = gerepilecopy(av, Q); R = Q+2;
     601             :     }
     602             :   }
     603       22099 :   return gerepilecopy(av, Q);
     604             : }
     605             : 
     606             : /********************************************************************/
     607             : /**                                                                **/
     608             : /**                          CONVERSIONS                           **/
     609             : /**                       (not memory clean)                       **/
     610             : /**                                                                **/
     611             : /********************************************************************/
     612             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     613             :  * but everything else is */
     614             : static GEN
     615      160076 : QXQ_to_mod(GEN x, GEN T)
     616             : {
     617             :   long d;
     618      160076 :   switch(typ(x))
     619             :   {
     620       59668 :     case t_INT:  return icopy(x);
     621        1057 :     case t_FRAC: return gcopy(x);
     622       99351 :     case t_POL:
     623       99351 :       d = degpol(x);
     624       99351 :       if (d < 0) return gen_0;
     625       97083 :       if (d == 0) return gcopy(gel(x,2));
     626       94934 :       return mkpolmod(RgX_copy(x), T);
     627           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     628             :              return NULL;/* LCOV_EXCL_LINE */
     629             :   }
     630             : }
     631             : /* pure shallow version */
     632             : GEN
     633      686199 : QXQ_to_mod_shallow(GEN x, GEN T)
     634             : {
     635             :   long d;
     636      686199 :   switch(typ(x))
     637             :   {
     638      440083 :     case t_INT:
     639      440083 :     case t_FRAC: return x;
     640      246116 :     case t_POL:
     641      246116 :       d = degpol(x);
     642      246116 :       if (d < 0) return gen_0;
     643      202261 :       if (d == 0) return gel(x,2);
     644      188254 :       return mkpolmod(x, T);
     645           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     646             :              return NULL;/* LCOV_EXCL_LINE */
     647             :   }
     648             : }
     649             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
     650             :  * Not memory clean because T not copied, but everything else is */
     651             : static GEN
     652       37100 : QXQX_to_mod(GEN z, GEN T)
     653             : {
     654       37100 :   long i,l = lg(z);
     655       37100 :   GEN x = cgetg(l,t_POL);
     656      185402 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     657       37100 :   x[1] = z[1]; return normalizepol_lg(x,l);
     658             : }
     659             : /* pure shallow version */
     660             : GEN
     661      152651 : QXQX_to_mod_shallow(GEN z, GEN T)
     662             : {
     663      152651 :   long i,l = lg(z);
     664      152651 :   GEN x = cgetg(l,t_POL);
     665      758238 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
     666      152651 :   x[1] = z[1]; return normalizepol_lg(x,l);
     667             : }
     668             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     669             : GEN
     670       12439 : QXQXV_to_mod(GEN V, GEN T)
     671             : {
     672       12439 :   long i, l = lg(V);
     673       12439 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     674       49539 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     675       12439 :   return z;
     676             : }
     677             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     678             : GEN
     679       11654 : QXQV_to_mod(GEN V, GEN T)
     680             : {
     681       11654 :   long i, l = lg(V);
     682       11654 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     683       23428 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
     684       11654 :   return z;
     685             : }
     686             : 
     687             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     688             : GEN
     689       13328 : QXQC_to_mod_shallow(GEN V, GEN T)
     690             : {
     691       13328 :   long i, l = lg(V);
     692       13328 :   GEN z = cgetg(l, t_COL);
     693       93940 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
     694       13328 :   return z;
     695             : }
     696             : 
     697             : GEN
     698        6237 : QXQM_to_mod_shallow(GEN V, GEN T)
     699             : {
     700        6237 :   long i, l = lg(V);
     701        6237 :   GEN z = cgetg(l, t_MAT);
     702       19565 :   for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
     703        6237 :   return z;
     704             : }
     705             : 
     706             : GEN
     707     5594739 : RgX_renormalize_lg(GEN x, long lx)
     708             : {
     709             :   long i;
     710     8478885 :   for (i = lx-1; i>1; i--)
     711     8042939 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     712     5594739 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     713     5594741 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     714             : }
     715             : 
     716             : GEN
     717     1463868 : RgV_to_RgX(GEN x, long v)
     718             : {
     719     1463868 :   long i, k = lg(x);
     720             :   GEN p;
     721             : 
     722     4129165 :   while (--k && gequal0(gel(x,k)));
     723     1463869 :   if (!k) return pol_0(v);
     724     1457056 :   i = k+2; p = cgetg(i,t_POL);
     725     1457053 :   p[1] = evalsigne(1) | evalvarn(v);
     726    10056784 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     727     1457053 :   return p;
     728             : }
     729             : GEN
     730      180589 : RgV_to_RgX_reverse(GEN x, long v)
     731             : {
     732      180589 :   long j, k, l = lg(x);
     733             :   GEN p;
     734             : 
     735      181926 :   for (k = 1; k < l; k++)
     736      181926 :     if (!gequal0(gel(x,k))) break;
     737      180589 :   if (k == l) return pol_0(v);
     738      180589 :   k -= 1;
     739      180589 :   l -= k;
     740      180589 :   x += k;
     741      180589 :   p = cgetg(l+1,t_POL);
     742      180589 :   p[1] = evalsigne(1) | evalvarn(v);
     743      996628 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     744      180589 :   return p;
     745             : }
     746             : 
     747             : /* return the (N-dimensional) vector of coeffs of p */
     748             : GEN
     749    11807900 : RgX_to_RgC(GEN x, long N)
     750             : {
     751             :   long i, l;
     752             :   GEN z;
     753    11807900 :   l = lg(x)-1; x++;
     754    11807900 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     755    11807900 :   z = cgetg(N+1,t_COL);
     756    74454599 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     757    22461938 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     758    11808061 :   return z;
     759             : }
     760             : GEN
     761     1244111 : Rg_to_RgC(GEN x, long N)
     762             : {
     763     1244111 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     764             : }
     765             : 
     766             : /* vector of polynomials (in v) whose coefs are given by the columns of x */
     767             : GEN
     768      299720 : RgM_to_RgXV(GEN x, long v)
     769     1278491 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
     770             : GEN
     771        6951 : RgM_to_RgXV_reverse(GEN x, long v)
     772       27804 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
     773             : 
     774             : /* matrix whose entries are given by the coeffs of the polynomials in
     775             :  * vector v (considered as degree n-1 polynomials) */
     776             : GEN
     777      287237 : RgV_to_RgM(GEN x, long n)
     778     1529135 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
     779             : 
     780             : GEN
     781       43519 : RgXV_to_RgM(GEN x, long n)
     782      254269 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
     783             : 
     784             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     785             : GEN
     786       17965 : RgM_to_RgXX(GEN x, long v,long w)
     787             : {
     788       17965 :   long j, lx = lg(x);
     789       17965 :   GEN y = cgetg(lx+1, t_POL);
     790       17965 :   y[1] = evalsigne(1) | evalvarn(v);
     791       17965 :   y++;
     792      102985 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     793       17965 :   return normalizepol_lg(--y, lx+1);
     794             : }
     795             : 
     796             : /* matrix whose entries are given by the coeffs of the polynomial v in
     797             :  * two variables (considered as degree n-1 polynomials) */
     798             : GEN
     799         308 : RgXX_to_RgM(GEN v, long n)
     800             : {
     801         308 :   long j, N = lg(v)-1;
     802         308 :   GEN y = cgetg(N, t_MAT);
     803        1001 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     804         308 :   return y;
     805             : }
     806             : 
     807             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     808             : GEN
     809       21704 : RgXY_swapspec(GEN x, long n, long w, long nx)
     810             : {
     811       21704 :   long j, ly = n+3;
     812       21704 :   GEN y = cgetg(ly, t_POL);
     813       21704 :   y[1] = evalsigne(1);
     814      264591 :   for (j=2; j<ly; j++)
     815             :   {
     816             :     long k;
     817      242887 :     GEN a = cgetg(nx+2,t_POL);
     818      242887 :     a[1] = evalsigne(1) | evalvarn(w);
     819     1192666 :     for (k=0; k<nx; k++)
     820             :     {
     821      949779 :       GEN xk = gel(x,k);
     822      949779 :       if (typ(xk)==t_POL)
     823      872227 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     824             :       else
     825       77552 :         gel(a,k+2) = j==2 ? xk: gen_0;
     826             :     }
     827      242887 :     gel(y,j) = normalizepol_lg(a, nx+2);
     828             :   }
     829       21704 :   return normalizepol_lg(y,ly);
     830             : }
     831             : 
     832             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     833             : GEN
     834         952 : RgXY_swap(GEN x, long n, long w)
     835             : {
     836         952 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     837         952 :   setvarn(z, varn(x)); return z;
     838             : }
     839             : 
     840             : long
     841         274 : RgXY_degreex(GEN b)
     842             : {
     843         274 :   long deg = 0, i;
     844         274 :   if (!signe(b)) return -1;
     845        1354 :   for (i = 2; i < lg(b); ++i)
     846             :   {
     847        1080 :     GEN bi = gel(b, i);
     848        1080 :     if (typ(bi) == t_POL)
     849        1016 :       deg = maxss(deg, degpol(bi));
     850             :   }
     851         274 :   return deg;
     852             : }
     853             : 
     854             : GEN
     855       16632 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
     856             : 
     857             : /* return (x % X^n). Shallow */
     858             : GEN
     859     6489978 : RgXn_red_shallow(GEN a, long n)
     860             : {
     861     6489978 :   long i, L = n+2, l = lg(a);
     862             :   GEN  b;
     863     6489978 :   if (L >= l) return a; /* deg(x) < n */
     864     4426495 :   b = cgetg(L, t_POL); b[1] = a[1];
     865    27000750 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     866     4426495 :   return normalizepol_lg(b,L);
     867             : }
     868             : 
     869             : GEN
     870         483 : RgXnV_red_shallow(GEN x, long n)
     871        2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
     872             : 
     873             : /* return (x * X^n). Shallow */
     874             : GEN
     875   125680998 : RgX_shift_shallow(GEN a, long n)
     876             : {
     877   125680998 :   long i, l = lg(a);
     878             :   GEN  b;
     879   125680998 :   if (l == 2 || !n) return a;
     880    85794816 :   l += n;
     881    85794816 :   if (n < 0)
     882             :   {
     883    51329608 :     if (l <= 2) return pol_0(varn(a));
     884    50013023 :     b = cgetg(l, t_POL); b[1] = a[1];
     885    50013214 :     a -= n;
     886   145021038 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     887             :   } else {
     888    34465208 :     b = cgetg(l, t_POL); b[1] = a[1];
     889    34466688 :     a -= n; n += 2;
     890    73021391 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     891   380469000 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     892             :   }
     893    84479902 :   return b;
     894             : }
     895             : /* return (x * X^n). */
     896             : GEN
     897     1573311 : RgX_shift(GEN a, long n)
     898             : {
     899     1573311 :   long i, l = lg(a);
     900             :   GEN  b;
     901     1573311 :   if (l == 2 || !n) return RgX_copy(a);
     902     1572947 :   l += n;
     903     1572947 :   if (n < 0)
     904             :   {
     905         847 :     if (l <= 2) return pol_0(varn(a));
     906         777 :     b = cgetg(l, t_POL); b[1] = a[1];
     907         777 :     a -= n;
     908        2947 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     909             :   } else {
     910     1572100 :     b = cgetg(l, t_POL); b[1] = a[1];
     911     1572100 :     a -= n; n += 2;
     912     3936747 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     913     4379970 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     914             :   }
     915     1572877 :   return b;
     916             : }
     917             : 
     918             : GEN
     919      316113 : RgX_rotate_shallow(GEN P, long k, long p)
     920             : {
     921      316113 :   long i, l = lgpol(P);
     922             :   GEN r;
     923      316113 :   if (signe(P)==0)
     924        1365 :     return pol_0(varn(P));
     925      314748 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     926     2090452 :   for(i=0; i<p; i++)
     927             :   {
     928     1775704 :     long s = 2+(i+k)%p;
     929     1775704 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
     930             :   }
     931      314748 :   return RgX_renormalize(r);
     932             : }
     933             : 
     934             : GEN
     935     2979315 : RgX_mulXn(GEN x, long d)
     936             : {
     937             :   pari_sp av;
     938             :   GEN z;
     939             :   long v;
     940     2979315 :   if (d >= 0) return RgX_shift(x, d);
     941     1467405 :   d = -d;
     942     1467405 :   v = RgX_val(x);
     943     1467405 :   if (v >= d) return RgX_shift(x, -d);
     944     1467398 :   av = avma;
     945     1467398 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
     946     1467398 :   return gerepileupto(av, z);
     947             : }
     948             : 
     949             : long
     950         588 : RgXV_maxdegree(GEN x)
     951             : {
     952         588 :   long d = -1, i, l = lg(x);
     953        4494 :   for (i = 1; i < l; i++)
     954        3906 :     d = maxss(d, degpol(gel(x,i)));
     955         588 :   return d;
     956             : }
     957             : 
     958             : long
     959     3552599 : RgX_val(GEN x)
     960             : {
     961     3552599 :   long i, lx = lg(x);
     962     3552599 :   if (lx == 2) return LONG_MAX;
     963     4476907 :   for (i = 2; i < lx; i++)
     964     4476851 :     if (!isexactzero(gel(x,i))) break;
     965     3552445 :   if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
     966     3552389 :   return i - 2;
     967             : }
     968             : long
     969    62130869 : RgX_valrem(GEN x, GEN *Z)
     970             : {
     971    62130869 :   long v, i, lx = lg(x);
     972    62130869 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
     973   104078530 :   for (i = 2; i < lx; i++)
     974   104085872 :     if (!isexactzero(gel(x,i))) break;
     975             :   /* possible with nonrational zeros */
     976    62147968 :   if (i == lx)
     977             :   {
     978          14 :     *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
     979          14 :     return LONG_MAX;
     980             :   }
     981    62147954 :   v = i - 2;
     982    62147954 :   *Z = RgX_shift_shallow(x, -v);
     983    62162536 :   return v;
     984             : }
     985             : long
     986      516887 : RgX_valrem_inexact(GEN x, GEN *Z)
     987             : {
     988             :   long v;
     989      516887 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
     990      516880 :   for (v = 0;; v++)
     991      527921 :     if (!gequal0(gel(x,2+v))) break;
     992      516880 :   if (Z) *Z = RgX_shift_shallow(x, -v);
     993      516880 :   return v;
     994             : }
     995             : 
     996             : GEN
     997       63707 : RgXQC_red(GEN x, GEN T)
     998      393981 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
     999             : 
    1000             : GEN
    1001        1204 : RgXQV_red(GEN x, GEN T)
    1002       27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
    1003             : 
    1004             : GEN
    1005       12229 : RgXQM_red(GEN x, GEN T)
    1006       75936 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
    1007             : 
    1008             : GEN
    1009         322 : RgXQM_mul(GEN P, GEN Q, GEN T)
    1010             : {
    1011         322 :   return RgXQM_red(RgM_mul(P, Q), T);
    1012             : }
    1013             : 
    1014             : GEN
    1015      444212 : RgXQX_red(GEN P, GEN T)
    1016             : {
    1017      444212 :   long i, l = lg(P);
    1018      444212 :   GEN Q = cgetg(l, t_POL);
    1019      444212 :   Q[1] = P[1];
    1020     2341432 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
    1021      444210 :   return normalizepol_lg(Q, l);
    1022             : }
    1023             : 
    1024             : GEN
    1025      565610 : RgX_deriv(GEN x)
    1026             : {
    1027      565610 :   long i,lx = lg(x)-1;
    1028             :   GEN y;
    1029             : 
    1030      565610 :   if (lx<3) return pol_0(varn(x));
    1031      563839 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
    1032     2772732 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
    1033      563822 :   y[1] = x[1]; return normalizepol_lg(y,i);
    1034             : }
    1035             : 
    1036             : GEN
    1037     1118514 : RgX_recipspec_shallow(GEN x, long l, long n)
    1038             : {
    1039             :   long i;
    1040     1118514 :   GEN z = cgetg(n+2,t_POL);
    1041     1118518 :   z[1] = 0; z += 2;
    1042    19755625 :   for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
    1043     1336681 :   for(   ; i<n; i++) gel(z, n-i-1) = gen_0;
    1044     1118518 :   return normalizepol_lg(z-2,n+2);
    1045             : }
    1046             : 
    1047             : GEN
    1048      629967 : RgXn_recip_shallow(GEN P, long n)
    1049             : {
    1050      629967 :   GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
    1051      629972 :   setvarn(Q, varn(P));
    1052      629972 :   return Q;
    1053             : }
    1054             : 
    1055             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1056             : GEN
    1057       15449 : RgX_recip(GEN x)
    1058             : {
    1059             :   long lx, i, j;
    1060       15449 :   GEN y = cgetg_copy(x, &lx);
    1061      133301 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
    1062       15449 :   return normalizepol_lg(y,lx);
    1063             : }
    1064             : /* shallow version */
    1065             : GEN
    1066       57050 : RgX_recip_shallow(GEN x)
    1067             : {
    1068             :   long lx, i, j;
    1069       57050 :   GEN y = cgetg_copy(x, &lx);
    1070      342139 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1071       57050 :   return normalizepol_lg(y,lx);
    1072             : }
    1073             : 
    1074             : GEN
    1075     2271945 : RgX_recip_i(GEN x)
    1076             : {
    1077             :   long lx, i, j;
    1078     2271945 :   GEN y = cgetg_copy(x, &lx);
    1079    13153501 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1080     2271947 :   return y;
    1081             : }
    1082             : /*******************************************************************/
    1083             : /*                                                                 */
    1084             : /*                      ADDITION / SUBTRACTION                     */
    1085             : /*                                                                 */
    1086             : /*******************************************************************/
    1087             : /* same variable */
    1088             : GEN
    1089    73464137 : RgX_add(GEN x, GEN y)
    1090             : {
    1091    73464137 :   long i, lx = lg(x), ly = lg(y);
    1092             :   GEN z;
    1093    73464137 :   if (ly <= lx) {
    1094    64574667 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1095   261210644 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1096   126024058 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1097    64517425 :     z = normalizepol_lg(z, lx);
    1098             :   } else {
    1099     8889470 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1100    40945523 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1101    30138969 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1102     8890669 :     z = normalizepol_lg(z, ly);
    1103             :   }
    1104    73461611 :   return z;
    1105             : }
    1106             : GEN
    1107    43036350 : RgX_sub(GEN x, GEN y)
    1108             : {
    1109    43036350 :   long i, lx = lg(x), ly = lg(y);
    1110             :   GEN z;
    1111    43036350 :   if (ly <= lx) {
    1112    25343215 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1113   143195232 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1114    47020019 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1115    25307630 :     z = normalizepol_lg(z, lx);
    1116             :   } else {
    1117    17693135 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1118    72142204 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1119    41969629 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1120    17646189 :     z = normalizepol_lg(z, ly);
    1121             :   }
    1122    43033718 :   return z;
    1123             : }
    1124             : GEN
    1125     4565376 : RgX_neg(GEN x)
    1126             : {
    1127     4565376 :   long i, lx = lg(x);
    1128     4565376 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    1129    38190018 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    1130     4565147 :   return y;
    1131             : }
    1132             : 
    1133             : GEN
    1134    17249718 : RgX_Rg_add(GEN y, GEN x)
    1135             : {
    1136             :   GEN z;
    1137    17249718 :   long lz = lg(y), i;
    1138    17249718 :   if (lz == 2) return scalarpol(x,varn(y));
    1139    14226004 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1140    14225958 :   gel(z,2) = gadd(gel(y,2),x);
    1141    56198297 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1142             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1143             :    * an inexact 0 */
    1144    14226102 :   return normalizepol_lg(z,lz);
    1145             : }
    1146             : GEN
    1147       65380 : RgX_Rg_add_shallow(GEN y, GEN x)
    1148             : {
    1149             :   GEN z;
    1150       65380 :   long lz = lg(y), i;
    1151       65380 :   if (lz == 2) return scalarpol(x,varn(y));
    1152       65380 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1153       65380 :   gel(z,2) = gadd(gel(y,2),x);
    1154      130869 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1155       65378 :   return normalizepol_lg(z,lz);
    1156             : }
    1157             : GEN
    1158      184205 : RgX_Rg_sub(GEN y, GEN x)
    1159             : {
    1160             :   GEN z;
    1161      184205 :   long lz = lg(y), i;
    1162      184205 :   if (lz == 2)
    1163             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1164        2548 :     long v = varn(y);
    1165        2548 :     if (isrationalzero(x)) return pol_0(v);
    1166          14 :     z = cgetg(3,t_POL);
    1167          14 :     z[1] = gequal0(x)? evalvarn(v)
    1168          14 :                    : evalvarn(v) | evalsigne(1);
    1169          14 :     gel(z,2) = gneg(x); return z;
    1170             :   }
    1171      181657 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1172      181657 :   gel(z,2) = gsub(gel(y,2),x);
    1173      474689 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1174      181659 :   return normalizepol_lg(z,lz);
    1175             : }
    1176             : GEN
    1177     2948776 : Rg_RgX_sub(GEN x, GEN y)
    1178             : {
    1179             :   GEN z;
    1180     2948776 :   long lz = lg(y), i;
    1181     2948776 :   if (lz == 2) return scalarpol(x,varn(y));
    1182     2947698 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1183     2947683 :   gel(z,2) = gsub(x, gel(y,2));
    1184     5419582 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1185     2947719 :   return normalizepol_lg(z,lz);
    1186             : }
    1187             : /*******************************************************************/
    1188             : /*                                                                 */
    1189             : /*                  KARATSUBA MULTIPLICATION                       */
    1190             : /*                                                                 */
    1191             : /*******************************************************************/
    1192             : #if 0
    1193             : /* to debug Karatsuba-like routines */
    1194             : GEN
    1195             : zx_debug_spec(GEN x, long nx)
    1196             : {
    1197             :   GEN z = cgetg(nx+2,t_POL);
    1198             :   long i;
    1199             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1200             :   z[1] = evalsigne(1); return z;
    1201             : }
    1202             : 
    1203             : GEN
    1204             : RgX_debug_spec(GEN x, long nx)
    1205             : {
    1206             :   GEN z = cgetg(nx+2,t_POL);
    1207             :   long i;
    1208             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1209             :   z[1] = evalsigne(1); return z;
    1210             : }
    1211             : #endif
    1212             : 
    1213             : /* generic multiplication */
    1214             : GEN
    1215     7389271 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1216             : {
    1217             :   GEN z, t;
    1218             :   long i;
    1219     7389271 :   if (nx == ny) {
    1220     1302541 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1221     4294052 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1222     1302534 :     return normalizepol_lg(z, nx+2);
    1223             :   }
    1224     6086730 :   if (ny < nx) {
    1225     5900195 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1226    23605143 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1227    14653351 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1228     5899568 :     return normalizepol_lg(z, nx+2);
    1229             :   } else {
    1230      186535 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1231     3674312 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1232      452510 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1233      186555 :     return normalizepol_lg(z, ny+2);
    1234             :   }
    1235             : }
    1236             : GEN
    1237      232948 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1238             : {
    1239             :   GEN z, t;
    1240             :   long i;
    1241      232948 :   if (nx == ny) {
    1242       12810 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1243     2183524 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1244       12810 :     return normalizepol_lg(z, nx+2);
    1245             :   }
    1246      220138 :   if (ny < nx) {
    1247      218360 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1248     3787494 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1249     2513584 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1250      218360 :     return normalizepol_lg(z, nx+2);
    1251             :   } else {
    1252        1778 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1253      330904 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1254       12222 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1255        1778 :     return normalizepol_lg(z, ny+2);
    1256             :   }
    1257             : }
    1258             : 
    1259             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1260             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1261             :  * t_VECSMALL, to hold this, which can be left on stack: gerepile
    1262             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1263             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1264             : static GEN
    1265    47186107 : RgXspec_kill0(GEN x, long lx)
    1266             : {
    1267    47186107 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
    1268             :   long i;
    1269   167658520 :   for (i=0; i <lx; i++)
    1270             :   {
    1271   120472408 :     GEN c = gel(x,i);
    1272   120472408 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1273             :   }
    1274    47186112 :   return z;
    1275             : }
    1276             : 
    1277             : INLINE GEN
    1278    84481397 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1279             : {
    1280    84481397 :   pari_sp av = avma;
    1281    84481397 :   GEN s = NULL;
    1282             :   long i;
    1283             : 
    1284   336885968 :   for (i=a; i<b; i++)
    1285   252417138 :     if (gel(y,i) && gel(x,-i))
    1286             :     {
    1287   174389233 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1288   174381819 :       s = s? gadd(s, t): t;
    1289             :     }
    1290    84468830 :   return s? gerepileupto(av, s): gen_0;
    1291             : }
    1292             : 
    1293             : /* assume nx >= ny > 0, return x * y * t^v */
    1294             : static GEN
    1295    17319722 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1296             : {
    1297             :   long i, lz, nz;
    1298             :   GEN z;
    1299             : 
    1300    17319722 :   x = RgXspec_kill0(x,nx);
    1301    17319651 :   y = RgXspec_kill0(y,ny);
    1302    17319710 :   lz = nx + ny + 1; nz = lz-2;
    1303    17319710 :   lz += v;
    1304    17319710 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1305    29051532 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1306    48713539 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1307    32835528 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1308    31393858 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1309    17319134 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1310             : }
    1311             : 
    1312             : /* return (x * X^d) + y. Assume d > 0 */
    1313             : GEN
    1314     7355270 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1315             : {
    1316             :   GEN x, y, xd, yd, zd;
    1317             :   long a, lz, nx, ny;
    1318             : 
    1319     7355270 :   if (!signe(x0)) return y0;
    1320     7253511 :   ny = lgpol(y0);
    1321     7253513 :   nx = lgpol(x0);
    1322     7253641 :   zd = (GEN)avma;
    1323     7253641 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1324     7253641 :   if (a <= 0)
    1325             :   {
    1326      841744 :     lz = nx+d+2;
    1327      841744 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1328     2121831 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1329      841750 :     x = zd + a;
    1330      857309 :     while (zd > x) gel(--zd,0) = gen_0;
    1331             :   }
    1332             :   else
    1333             :   {
    1334     6411897 :     xd = new_chunk(d); yd = y+d;
    1335     6411911 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1336     6411832 :     lz = (a>nx)? ny+2: lg(x)+d;
    1337    34736087 :     x += 2; while (xd > x) *--zd = *--xd;
    1338             :   }
    1339    18798459 :   while (yd > y) *--zd = *--yd;
    1340     7253582 :   *--zd = x0[1];
    1341     7253582 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1342             : }
    1343             : GEN
    1344      525507 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1345             : {
    1346             :   GEN x, y, xd, yd, zd;
    1347             :   long a, lz, nx, ny;
    1348             : 
    1349      525507 :   if (!signe(x0)) return RgX_copy(y0);
    1350      524737 :   nx = lgpol(x0);
    1351      524737 :   ny = lgpol(y0);
    1352      524737 :   zd = (GEN)avma;
    1353      524737 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1354      524737 :   if (a <= 0)
    1355             :   {
    1356      291789 :     lz = nx+d+2;
    1357      291789 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1358     4294362 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1359      291789 :     x = zd + a;
    1360      757955 :     while (zd > x) gel(--zd,0) = gen_0;
    1361             :   }
    1362             :   else
    1363             :   {
    1364      232948 :     xd = new_chunk(d); yd = y+d;
    1365      232948 :     x = RgX_addspec(x,yd, nx,a);
    1366      232948 :     lz = (a>nx)? ny+2: lg(x)+d;
    1367     8607590 :     x += 2; while (xd > x) *--zd = *--xd;
    1368             :   }
    1369     2693813 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1370      524737 :   *--zd = x0[1];
    1371      524737 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1372             : }
    1373             : 
    1374             : /* return x * y mod t^n */
    1375             : static GEN
    1376     6269741 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1377             : {
    1378     6269741 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1379             :   GEN z;
    1380     6269741 :   if (lx < 0) return pol_0(varn(x));
    1381     6269741 :   if (ly < 0) return pol_0(varn(x));
    1382     6269741 :   z = cgetg(lz, t_POL) + 2;
    1383     6269741 :   x+=2; if (lx > n) lx = n;
    1384     6269741 :   y+=2; if (ly > n) ly = n;
    1385     6269741 :   z[-1] = x[-1];
    1386     6269741 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1387     6269741 :   x = RgXspec_kill0(x, lx);
    1388     6269741 :   y = RgXspec_kill0(y, ly);
    1389             :   /* x:y:z [i] = term of degree i */
    1390    24956444 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1391    11035618 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1392     6314553 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1393     6269741 :   return normalizepol_lg(z - 2, lz);
    1394             : }
    1395             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1396             : static GEN
    1397     8856583 : RgXn_mul2(GEN f, GEN g, long n)
    1398             : {
    1399     8856583 :   pari_sp av = avma;
    1400             :   GEN fe,fo, ge,go, l,h,m;
    1401             :   long n0, n1;
    1402     8856583 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1403     6305161 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1404       35420 :   n0 = n>>1; n1 = n-n0;
    1405       35420 :   RgX_even_odd(f, &fe, &fo);
    1406       35420 :   RgX_even_odd(g, &ge, &go);
    1407       35420 :   l = RgXn_mul2(fe,ge,n1);
    1408       35420 :   h = RgXn_mul2(fo,go,n0);
    1409       35420 :   m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1410             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1411             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1412       35420 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1413             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1414       35420 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1415       35420 :   m = RgX_inflate(m,2);
    1416             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1417       35420 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1418       35420 :   h = RgX_inflate(h,2);
    1419       35420 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1420       35420 :   return gerepileupto(av, h);
    1421             : }
    1422             : /* (f*g) \/ x^n */
    1423             : static GEN
    1424     1506881 : RgX_mulhigh_i2(GEN f, GEN g, long n)
    1425             : {
    1426     1506881 :   long d = degpol(f)+degpol(g) + 1 - n;
    1427             :   GEN h;
    1428     1506881 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1429       21310 :   h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
    1430             :                                   RgX_recip_i(g), d));
    1431       21310 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1432             : }
    1433             : 
    1434             : /* (f*g) \/ x^n */
    1435             : static GEN
    1436           0 : RgX_sqrhigh_i2(GEN f, long n)
    1437             : {
    1438           0 :   long d = 2*degpol(f)+ 1 - n;
    1439             :   GEN h;
    1440           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1441           0 :   h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
    1442           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1443             : }
    1444             : 
    1445             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1446             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1447             :  * Only c, c0, c1, c2 are genuine GEN.
    1448             :  */
    1449             : GEN
    1450    18057127 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1451             : {
    1452             :   GEN a0, c, c0;
    1453    18057127 :   long n0, n0a, i, v = 0;
    1454             :   pari_sp av;
    1455             : 
    1456    24441162 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1457    23939779 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1458    18057110 :   if (na < nb) swapspec(a,b, na,nb);
    1459    18057110 :   if (!nb) return pol_0(0);
    1460             : 
    1461    17809771 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1462      490013 :   RgX_shift_inplace_init(v);
    1463      490048 :   i = (na>>1); n0 = na-i; na = i;
    1464      490048 :   av = avma; a0 = a+n0; n0a = n0;
    1465     1345548 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1466             : 
    1467      490048 :   if (nb > n0)
    1468             :   {
    1469             :     GEN b0,c1,c2;
    1470             :     long n0b;
    1471             : 
    1472      488677 :     nb -= n0; b0 = b+n0; n0b = n0;
    1473     1456640 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1474      488677 :     c = RgX_mulspec(a,b,n0a,n0b);
    1475      488677 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1476             : 
    1477      488677 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1478      488677 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1479             : 
    1480      488677 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1481      488677 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1482      488677 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1483             :   }
    1484             :   else
    1485             :   {
    1486        1371 :     c = RgX_mulspec(a,b,n0a,nb);
    1487        1371 :     c0 = RgX_mulspec(a0,b,na,nb);
    1488             :   }
    1489      490048 :   c0 = RgX_addmulXn(c0,c,n0);
    1490      490048 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1491             : }
    1492             : 
    1493             : INLINE GEN
    1494       48160 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1495             : {
    1496       48160 :   pari_sp av = avma;
    1497       48160 :   GEN s = NULL;
    1498       48160 :   long j, l = (i+1)>>1;
    1499      136314 :   for (j=a; j<l; j++)
    1500             :   {
    1501       88154 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1502       88154 :     if (xj && xx)
    1503             :     {
    1504       80888 :       GEN t = gmul(xj, xx);
    1505       80888 :       s = s? gadd(s, t): t;
    1506             :     }
    1507             :   }
    1508       48160 :   if (s) s = gshift(s,1);
    1509       48160 :   if ((i&1) == 0)
    1510             :   {
    1511       27816 :     GEN t = gel(x, i>>1);
    1512       27816 :     if (t) {
    1513       25520 :       t = gsqr(t);
    1514       25520 :       s = s? gadd(s, t): t;
    1515             :     }
    1516             :   }
    1517       48160 :   return s? gerepileupto(av,s): gen_0;
    1518             : }
    1519             : static GEN
    1520        7472 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1521             : {
    1522             :   long i, lz, nz;
    1523             :   GEN z;
    1524             : 
    1525        7472 :   if (!nx) return pol_0(0);
    1526        7472 :   x = RgXspec_kill0(x,nx);
    1527        7472 :   lz = (nx << 1) + 1, nz = lz-2;
    1528        7472 :   lz += v;
    1529        7472 :   z = cgetg(lz,t_POL) + 2;
    1530       11140 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1531       35288 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1532       27816 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1533        7472 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1534             : }
    1535             : /* return x^2 mod t^n */
    1536             : static GEN
    1537           0 : RgXn_sqr_basecase(GEN x, long n)
    1538             : {
    1539           0 :   long i, lz = n+2, lx = lgpol(x);
    1540             :   GEN z;
    1541           0 :   if (lx < 0) return pol_0(varn(x));
    1542           0 :   z = cgetg(lz, t_POL);
    1543           0 :   z[1] = x[1];
    1544           0 :   x+=2; if (lx > n) lx = n;
    1545           0 :   x = RgXspec_kill0(x,lx);
    1546           0 :   z+=2;/* x:z [i] = term of degree i */
    1547           0 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1548           0 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1549           0 :   z -= 2; return normalizepol_lg(z, lz);
    1550             : }
    1551             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1552             : static GEN
    1553         315 : RgXn_sqr2(GEN f, long n)
    1554             : {
    1555         315 :   pari_sp av = avma;
    1556             :   GEN fe,fo, l,h,m;
    1557             :   long n0, n1;
    1558         315 :   if (2*degpol(f) < n) return RgX_sqr_i(f);
    1559           0 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1560           0 :   n0 = n>>1; n1 = n-n0;
    1561           0 :   RgX_even_odd(f, &fe, &fo);
    1562           0 :   l = RgXn_sqr(fe,n1);
    1563           0 :   h = RgXn_sqr(fo,n0);
    1564           0 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1565             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1566             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1567           0 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1568             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1569           0 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1570           0 :   m = RgX_inflate(m,2);
    1571             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1572           0 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1573           0 :   h = RgX_inflate(h,2);
    1574           0 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1575           0 :   return gerepileupto(av, h);
    1576             : }
    1577             : GEN
    1578        7511 : RgX_sqrspec(GEN a, long na)
    1579             : {
    1580             :   GEN a0, c, c0, c1;
    1581        7511 :   long n0, n0a, i, v = 0;
    1582             :   pari_sp av;
    1583             : 
    1584        9345 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1585        7511 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1586          39 :   RgX_shift_inplace_init(v);
    1587          39 :   i = (na>>1); n0 = na-i; na = i;
    1588          39 :   av = avma; a0 = a+n0; n0a = n0;
    1589          39 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1590             : 
    1591          39 :   c = RgX_sqrspec(a,n0a);
    1592          39 :   c0 = RgX_sqrspec(a0,na);
    1593          39 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1594          39 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1595          39 :   c0 = RgX_addmulXn(c0,c,n0);
    1596          39 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1597             : }
    1598             : 
    1599             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1600             : GEN
    1601     1684320 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1602             : {
    1603     1684320 :   GEN z = RgX_mul(A, B);
    1604     1684309 :   if (a < b)
    1605       10423 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1606     1673886 :   else if (a > b)
    1607     1056373 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1608             :   else
    1609      617513 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1610     1684308 :   return z;
    1611             : }
    1612             : 
    1613             : GEN
    1614    16588393 : RgX_mul_i(GEN x, GEN y)
    1615             : {
    1616    16588393 :   GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
    1617    16587801 :   setvarn(z, varn(x)); return z;
    1618             : }
    1619             : 
    1620             : GEN
    1621        7433 : RgX_sqr_i(GEN x)
    1622             : {
    1623        7433 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1624        7433 :   setvarn(z,varn(x)); return z;
    1625             : }
    1626             : 
    1627             : /*******************************************************************/
    1628             : /*                                                                 */
    1629             : /*                               DIVISION                          */
    1630             : /*                                                                 */
    1631             : /*******************************************************************/
    1632             : GEN
    1633     3929976 : RgX_Rg_divexact(GEN x, GEN y) {
    1634     3929976 :   long i, lx = lg(x);
    1635             :   GEN z;
    1636     3929976 :   if (lx == 2) return gcopy(x);
    1637     3880874 :   switch(typ(y))
    1638             :   {
    1639     3851785 :     case t_INT:
    1640     3851785 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1641     3682376 :       break;
    1642        4123 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1643             :   }
    1644     3707342 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1645    23638422 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1646     3693651 :   return z;
    1647             : }
    1648             : GEN
    1649    28191462 : RgX_Rg_div(GEN x, GEN y) {
    1650    28191462 :   long i, lx = lg(x);
    1651             :   GEN z;
    1652    28191462 :   if (lx == 2) return gcopy(x);
    1653    27932071 :   switch(typ(y))
    1654             :   {
    1655    18787038 :     case t_INT:
    1656    18787038 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1657     3729688 :       break;
    1658        4837 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1659             :   }
    1660    12869884 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1661    48660550 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1662    12869464 :   return normalizepol_lg(z, lx);
    1663             : }
    1664             : GEN
    1665       13237 : RgX_normalize(GEN x)
    1666             : {
    1667       13237 :   GEN z, d = NULL;
    1668       13237 :   long i, n = lg(x)-1;
    1669       13237 :   for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
    1670       13237 :   if (i == 1) return pol_0(varn(x));
    1671       13237 :   if (i == n && isint1(d)) return x;
    1672        5586 :   n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
    1673       13034 :   for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
    1674        5586 :   gel(z,n) = Rg_get_1(d); return z;
    1675             : }
    1676             : GEN
    1677       10318 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
    1678             : GEN
    1679      229569 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1680             : {
    1681      229569 :   long l = lg(a), i;
    1682             :   GEN a0, z0, z;
    1683             : 
    1684      229569 :   if (l <= 3)
    1685             :   {
    1686           0 :     if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
    1687           0 :     return pol_0(0);
    1688             :   }
    1689      229569 :   z = cgetg(l-1, t_POL);
    1690      229575 :   z[1] = a[1];
    1691      229575 :   a0 = a + l-1;
    1692      229575 :   z0 = z + l-2; *z0 = *a0--;
    1693     2872953 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1694             :   {
    1695     2643384 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1696     2643378 :     gel(z0,0) = t;
    1697             :   }
    1698      229569 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1699      229569 :   return z;
    1700             : }
    1701             : /* Polynomial division x / y:
    1702             :  *   if pr = ONLY_REM return remainder, otherwise return quotient
    1703             :  *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
    1704             :  *   if pr != NULL set *pr to remainder, as the last object on stack */
    1705             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1706             : static GEN
    1707    18288231 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
    1708             : {
    1709             :   pari_sp avy, av, av1;
    1710             :   long dx,dy,dz,i,j,sx,lr;
    1711             :   GEN z,p1,p2,rem,y_lead,mod,p;
    1712             :   GEN (*f)(GEN,GEN);
    1713             : 
    1714    18288231 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1715             : 
    1716    18288231 :   dy = degpol(y);
    1717    18288184 :   y_lead = gel(y,dy+2);
    1718    18288184 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1719             :   {
    1720           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1721           0 :     for (dy--; dy>=0; dy--)
    1722             :     {
    1723           0 :       y_lead = gel(y,dy+2);
    1724           0 :       if (!gequal0(y_lead)) break;
    1725             :     }
    1726             :   }
    1727    18288138 :   if (!dy) /* y is constant */
    1728             :   {
    1729        3114 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1730        3107 :     z = RgX_Rg_div(x, y_lead);
    1731        3107 :     if (pr == ONLY_DIVIDES) return z;
    1732        1224 :     if (pr) *pr = pol_0(varn(x));
    1733        1224 :     return z;
    1734             :   }
    1735    18285024 :   dx = degpol(x);
    1736    18284963 :   if (dx < dy)
    1737             :   {
    1738     1915352 :     if (pr == ONLY_REM) return RgX_copy(x);
    1739      374673 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1740      374652 :     z = pol_0(varn(x));
    1741      374652 :     if (pr) *pr = RgX_copy(x);
    1742      374652 :     return z;
    1743             :   }
    1744             : 
    1745             :   /* x,y in R[X], y non constant */
    1746    16369611 :   av = avma;
    1747    16369611 :   p = NULL;
    1748    16369611 :   if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1749             :   {
    1750      226359 :     z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
    1751      226359 :     if (!z) return gc_NULL(av);
    1752      226359 :     z = FpX_to_mod(z, p);
    1753      226359 :     if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
    1754      121821 :       return gerepileupto(av, z);
    1755      104538 :     *pr = FpX_to_mod(*pr, p);
    1756      104538 :     return gc_all(av, 2, &z, pr);
    1757             :   }
    1758    16143412 :   switch(typ(y_lead))
    1759             :   {
    1760           0 :     case t_REAL:
    1761           0 :       y_lead = ginv(y_lead);
    1762           0 :       f = gmul; mod = NULL;
    1763           0 :       break;
    1764        2143 :     case t_INTMOD:
    1765        2143 :     case t_POLMOD: y_lead = ginv(y_lead);
    1766        2143 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1767        2143 :       break;
    1768    16141269 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1769    16141243 :       f = gdiv; mod = NULL;
    1770             :   }
    1771             : 
    1772    16143386 :   if (y_lead == NULL)
    1773    14320960 :     p2 = gel(x,dx+2);
    1774             :   else {
    1775             :     for(;;) {
    1776     1822426 :       p2 = f(gel(x,dx+2),y_lead);
    1777     1822432 :       p2 = simplify_shallow(p2);
    1778     1822432 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1779             :     }
    1780     1822432 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1781             :     {
    1782           0 :       if (pr == ONLY_DIVIDES) {
    1783           0 :         set_avma(av);
    1784           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1785             :       }
    1786           0 :       if (pr == ONLY_REM)
    1787             :       {
    1788           0 :         if (dx < 0)
    1789           0 :           return gerepilecopy(av, scalarpol(p2, varn(x)));
    1790             :         else
    1791             :         {
    1792             :           GEN t;
    1793           0 :           set_avma(av);
    1794           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1795           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1796           0 :           return t;
    1797             :         }
    1798             :       }
    1799           0 :       if (pr) /* cf ONLY_REM above */
    1800             :       {
    1801           0 :         if (dx < 0)
    1802             :         {
    1803           0 :           p2 = gclone(p2);
    1804           0 :           set_avma(av);
    1805           0 :           z = pol_0(varn(x));
    1806           0 :           x = scalarpol(p2, varn(x));
    1807           0 :           gunclone(p2);
    1808             :         }
    1809             :         else
    1810             :         {
    1811             :           GEN t;
    1812           0 :           set_avma(av);
    1813           0 :           z = pol_0(varn(x));
    1814           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1815           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1816           0 :           x = t;
    1817             :         }
    1818           0 :         *pr = x;
    1819             :       }
    1820             :       else
    1821             :       {
    1822           0 :         set_avma(av);
    1823           0 :         z = pol_0(varn(x));
    1824             :       }
    1825           0 :       return z;
    1826             :     }
    1827             :   }
    1828             :   /* dx >= dy */
    1829    16143392 :   avy = avma;
    1830    16143392 :   dz = dx-dy;
    1831    16143392 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1832    16143329 :   x += 2;
    1833    16143329 :   z += 2;
    1834    16143329 :   y += 2;
    1835    16143329 :   gel(z,dz) = gcopy(p2);
    1836             : 
    1837    44670505 :   for (i=dx-1; i>=dy; i--)
    1838             :   {
    1839    28527300 :     av1=avma; p1=gel(x,i);
    1840  1093145729 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1841    28521621 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1842             : 
    1843    28521621 :     if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
    1844             :     else
    1845    18956818 :       p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
    1846    28526528 :     gel(z,i-dy) = p1;
    1847             :   }
    1848    16143205 :   if (!pr) return gerepileupto(av,z-2);
    1849             : 
    1850     9428033 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1851     9428021 :   for (sx=0; ; i--)
    1852             :   {
    1853      360171 :     p1 = gel(x,i);
    1854             :     /* we always enter this loop at least once */
    1855    26209855 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1856     9786762 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1857     9786762 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
    1858     3500727 :     if (!isexactzero(p1)) break;
    1859     3423735 :     if (!i) break;
    1860      360165 :     set_avma(av1);
    1861             :   }
    1862     9427096 :   if (pr == ONLY_DIVIDES)
    1863             :   {
    1864        1890 :     if (sx) return gc_NULL(av);
    1865        1869 :     set_avma((pari_sp)rem); return gerepileupto(av,z-2);
    1866             :   }
    1867     9425206 :   lr=i+3; rem -= lr;
    1868     9425206 :   if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
    1869     9381608 :   else p1 = gerepileupto((pari_sp)rem,p1);
    1870     9426119 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1871     9426080 :   rem[1] = z[-1];
    1872     9426080 :   rem += 2;
    1873     9426080 :   gel(rem,i) = p1;
    1874    14179040 :   for (i--; i>=0; i--)
    1875             :   {
    1876     4752947 :     av1=avma; p1 = gel(x,i);
    1877    13443728 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1878     4752364 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1879     4752650 :     gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
    1880             :   }
    1881     9426093 :   rem -= 2;
    1882     9426093 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1883     9426337 :   if (pr == ONLY_REM) return gerepileupto(av,rem);
    1884     5478030 :   z -= 2;
    1885             :   {
    1886     5478030 :     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
    1887     5478030 :     gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
    1888             :   }
    1889             : }
    1890             : 
    1891             : GEN
    1892    12799286 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1893             : {
    1894    12799286 :   if (pr == ONLY_REM) return RgX_rem(x, y);
    1895    12799286 :   return RgX_divrem_i(x, y, pr);
    1896             : }
    1897             : 
    1898             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1899             : GEN
    1900      133569 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1901             : {
    1902             :   long vx, dx, dy, dz, i, j, sx, lr;
    1903             :   pari_sp av0, av, tetpil;
    1904             :   GEN z,p1,rem,lead;
    1905             : 
    1906      133569 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1907      133569 :   vx = varn(x);
    1908      133569 :   dx = degpol(x);
    1909      133569 :   dy = degpol(y);
    1910      133569 :   if (dx < dy)
    1911             :   {
    1912       26389 :     if (pr)
    1913             :     {
    1914       26389 :       av0 = avma; x = RgXQX_red(x, T);
    1915       26389 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
    1916       26389 :       if (pr == ONLY_REM) return x;
    1917           0 :       *pr = x;
    1918             :     }
    1919           0 :     return pol_0(vx);
    1920             :   }
    1921      107180 :   lead = leading_coeff(y);
    1922      107180 :   if (!dy) /* y is constant */
    1923             :   {
    1924         546 :     if (pr && pr != ONLY_DIVIDES)
    1925             :     {
    1926           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1927           0 :       *pr = pol_0(vx);
    1928             :     }
    1929         546 :     if (gequal1(lead)) return RgX_copy(x);
    1930           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1931           0 :     return gerepile(av0,tetpil,RgXQX_red(x,T));
    1932             :   }
    1933      106634 :   av0 = avma; dz = dx-dy;
    1934      106634 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1935      106634 :   set_avma(av0);
    1936      106634 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1937      106634 :   x += 2; y += 2; z += 2;
    1938             : 
    1939      106634 :   p1 = gel(x,dx); av = avma;
    1940      106634 :   gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    1941      548235 :   for (i=dx-1; i>=dy; i--)
    1942             :   {
    1943      441601 :     av=avma; p1=gel(x,i);
    1944     2331178 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1945      441573 :     if (lead) p1 = gmul(grem(p1, T), lead);
    1946      441564 :     tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
    1947             :   }
    1948      106634 :   if (!pr) { guncloneNULL(lead); return z-2; }
    1949             : 
    1950      105388 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    1951      105388 :   for (sx=0; ; i--)
    1952             :   {
    1953       72467 :     p1 = gel(x,i);
    1954      613390 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1955      177851 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    1956       99480 :     if (!i) break;
    1957       72467 :     set_avma(av);
    1958             :   }
    1959      105388 :   if (pr == ONLY_DIVIDES)
    1960             :   {
    1961       27006 :     guncloneNULL(lead);
    1962       27006 :     if (sx) return gc_NULL(av0);
    1963       24500 :     return gc_const((pari_sp)rem, z-2);
    1964             :   }
    1965       78382 :   lr=i+3; rem -= lr;
    1966       78382 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1967       78382 :   rem[1] = z[-1];
    1968       78382 :   p1 = gerepile((pari_sp)rem,tetpil,p1);
    1969       78382 :   rem += 2; gel(rem,i) = p1;
    1970      168034 :   for (i--; i>=0; i--)
    1971             :   {
    1972       89652 :     av=avma; p1 = gel(x,i);
    1973      304586 :     for (j=0; j<=i && j<=dz; j++)
    1974      214934 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1975       89652 :     tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
    1976             :   }
    1977       78382 :   rem -= 2;
    1978       78382 :   guncloneNULL(lead);
    1979       78382 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1980       78382 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
    1981         161 :   *pr = rem; return z-2;
    1982             : }
    1983             : 
    1984             : /*******************************************************************/
    1985             : /*                                                                 */
    1986             : /*                        PSEUDO-DIVISION                          */
    1987             : /*                                                                 */
    1988             : /*******************************************************************/
    1989             : INLINE GEN
    1990     1305765 : rem(GEN c, GEN T)
    1991             : {
    1992     1305765 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    1993     1305766 :   return c;
    1994             : }
    1995             : 
    1996             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    1997             : int
    1998       12414 : ZXQX_dvd(GEN x, GEN y, GEN T)
    1999             : {
    2000             :   long dx, dy, dz, i, p, T_ismonic;
    2001       12414 :   pari_sp av = avma, av2;
    2002             :   GEN y_lead;
    2003             : 
    2004       12414 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    2005       12414 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2006       12414 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    2007             :   /* if monic, no point in using pseudo-division */
    2008       12414 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    2009        9936 :   T_ismonic = gequal1(leading_coeff(T));
    2010        9936 :   dx = degpol(x);
    2011        9936 :   if (dx < dy) return !signe(x);
    2012        9936 :   (void)new_chunk(2);
    2013        9936 :   x = RgX_recip_i(x)+2;
    2014        9936 :   y = RgX_recip_i(y)+2;
    2015             :   /* pay attention to sparse divisors */
    2016       20424 :   for (i = 1; i <= dy; i++)
    2017       10488 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    2018        9936 :   dz = dx-dy; p = dz+1;
    2019        9936 :   av2 = avma;
    2020             :   for (;;)
    2021       59543 :   {
    2022       69479 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    2023       69479 :     x0 = gneg(x0); p--;
    2024       69479 :     m = gcdii(cx, y0);
    2025       69478 :     if (!equali1(m))
    2026             :     {
    2027       67605 :       x0 = gdiv(x0, m);
    2028       67605 :       y0 = diviiexact(y0, m);
    2029       67604 :       if (equali1(y0)) y0 = NULL;
    2030             :     }
    2031      143173 :     for (i=1; i<=dy; i++)
    2032             :     {
    2033       73694 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2034       73694 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    2035       73696 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2036       73696 :       gel(x,i) = c;
    2037             :     }
    2038      625510 :     for (   ; i<=dx; i++)
    2039             :     {
    2040      556031 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2041      556031 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2042      556031 :       gel(x,i) = c;
    2043             :     }
    2044       79570 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    2045       69479 :     if (dx < dy) break;
    2046       59543 :     if (gc_needed(av2,1))
    2047             :     {
    2048          28 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    2049          28 :       gerepilecoeffs(av2,x,dx+1);
    2050             :     }
    2051             :   }
    2052        9936 :   return gc_bool(av, dx < 0);
    2053             : }
    2054             : 
    2055             : /* T either NULL or a t_POL. */
    2056             : GEN
    2057      133918 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    2058             : {
    2059      133918 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    2060      133918 :   pari_sp av = avma, av2;
    2061             :   GEN y_lead;
    2062             : 
    2063      133918 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    2064      133918 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2065             :   /* if monic, no point in using pseudo-division */
    2066      133919 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    2067      128235 :   dx = degpol(x);
    2068      128235 :   if (dx < dy) return RgX_copy(x);
    2069      128235 :   (void)new_chunk(2);
    2070      128235 :   x = RgX_recip_i(x)+2;
    2071      128235 :   y = RgX_recip_i(y)+2;
    2072             :   /* pay attention to sparse divisors */
    2073      365745 :   for (i = 1; i <= dy; i++)
    2074      237510 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2075      128235 :   dz = dx-dy; p = dz+1;
    2076      128235 :   av2 = avma;
    2077             :   for (;;)
    2078             :   {
    2079      356159 :     gel(x,0) = gneg(gel(x,0)); p--;
    2080     1055697 :     for (i=1; i<=dy; i++)
    2081             :     {
    2082      699532 :       GEN c = gmul(y_lead, gel(x,i));
    2083      699483 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    2084      699523 :       gel(x,i) = rem(c, T);
    2085             :     }
    2086      845844 :     for (   ; i<=dx; i++)
    2087             :     {
    2088      489686 :       GEN c = gmul(y_lead, gel(x,i));
    2089      489668 :       gel(x,i) = rem(c, T);
    2090             :     }
    2091      382979 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    2092      356160 :     if (dx < dy) break;
    2093      227924 :     if (gc_needed(av2,1))
    2094             :     {
    2095           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    2096           0 :       gerepilecoeffs(av2,x,dx+1);
    2097             :     }
    2098             :   }
    2099      128236 :   if (dx < 0) return pol_0(vx);
    2100      126717 :   lx = dx+3; x -= 2;
    2101      126717 :   x[0] = evaltyp(t_POL) | evallg(lx);
    2102      126717 :   x[1] = evalsigne(1) | evalvarn(vx);
    2103      126717 :   x = RgX_recip_i(x);
    2104      126717 :   if (p)
    2105             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    2106       12011 :     GEN t = y_lead;
    2107       12011 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    2108           0 :       t = RgXQ_powu(t, p, T);
    2109             :     else
    2110       12011 :       t = gpowgs(t, p);
    2111       32282 :     for (i=2; i<lx; i++)
    2112             :     {
    2113       20272 :       GEN c = gmul(gel(x,i), t);
    2114       20270 :       gel(x,i) = rem(c,T);
    2115             :     }
    2116       12010 :     if (!T) return gerepileupto(av, x);
    2117             :   }
    2118      114706 :   return gerepilecopy(av, x);
    2119             : }
    2120             : 
    2121             : GEN
    2122      133918 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2123             : 
    2124             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2125             : GEN
    2126       11705 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2127             : {
    2128       11705 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2129       11705 :   pari_sp av = avma, av2;
    2130             :   GEN z, r, ypow, y_lead;
    2131             : 
    2132       11705 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2133       11705 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2134       11705 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2135        7407 :   dx = degpol(x);
    2136        7407 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2137        7407 :   if (dx == dy)
    2138             :   {
    2139          84 :     GEN x_lead = gel(x,lg(x)-1);
    2140          84 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2141          84 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2142          84 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2143          84 :     *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
    2144             :   }
    2145        7323 :   (void)new_chunk(2);
    2146        7323 :   x = RgX_recip_i(x)+2;
    2147        7323 :   y = RgX_recip_i(y)+2;
    2148             :   /* pay attention to sparse divisors */
    2149       37438 :   for (i = 1; i <= dy; i++)
    2150       30115 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2151        7323 :   dz = dx-dy; p = dz+1;
    2152        7323 :   lz = dz+3;
    2153        7323 :   z = cgetg(lz, t_POL);
    2154        7323 :   z[1] = evalsigne(1) | evalvarn(vx);
    2155       27969 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2156        7323 :   ypow = new_chunk(dz+1);
    2157        7323 :   gel(ypow,0) = gen_1;
    2158        7323 :   gel(ypow,1) = y_lead;
    2159       13323 :   for (i=2; i<=dz; i++)
    2160             :   {
    2161        6000 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2162        6000 :     gel(ypow,i) = rem(c,T);
    2163             :   }
    2164        7323 :   av2 = avma;
    2165        7323 :   for (iz=2;;)
    2166             :   {
    2167        7911 :     p--;
    2168       15234 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2169       66064 :     for (i=1; i<=dy; i++)
    2170             :     {
    2171       50830 :       GEN c = gmul(y_lead, gel(x,i));
    2172       50830 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2173       50830 :       gel(x,i) = rem(c, T);
    2174             :     }
    2175       39526 :     for (   ; i<=dx; i++)
    2176             :     {
    2177       24292 :       GEN c = gmul(y_lead, gel(x,i));
    2178       24292 :       gel(x,i) = rem(c,T);
    2179             :     }
    2180       15234 :     x++; dx--;
    2181       20646 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2182       15234 :     if (dx < dy) break;
    2183        7911 :     if (gc_needed(av2,1))
    2184             :     {
    2185           0 :       GEN X = x-2;
    2186           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2187           0 :       X[0] = evaltyp(t_POL)|evallg(dx+3); X[1] = z[1]; /* hack */
    2188           0 :       gerepileall(av2,2, &X, &z); x = X+2;
    2189             :     }
    2190             :   }
    2191       13203 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2192        7323 :   if (dx < 0)
    2193          98 :     x = pol_0(vx);
    2194             :   else
    2195             :   {
    2196        7225 :     lx = dx+3; x -= 2;
    2197        7225 :     x[0] = evaltyp(t_POL) | evallg(lx);
    2198        7225 :     x[1] = evalsigne(1) | evalvarn(vx);
    2199        7225 :     x = RgX_recip_i(x);
    2200             :   }
    2201        7323 :   z = RgX_recip_i(z);
    2202        7323 :   r = x;
    2203        7323 :   if (p)
    2204             :   {
    2205        3269 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2206        3269 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2207             :   }
    2208        7323 :   *ptr = r; return gc_all(av, 2, &z, ptr);
    2209             : }
    2210             : GEN
    2211       11467 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2212       11467 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2213             : 
    2214             : GEN
    2215           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2216           0 : { return RgXQX_red(RgX_mul(x,y), T); }
    2217             : GEN
    2218   486254225 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
    2219             : GEN
    2220      141344 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
    2221             : GEN
    2222       19285 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
    2223             : GEN
    2224          35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
    2225             : GEN
    2226         133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
    2227             : 
    2228             : GEN
    2229           0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
    2230             : 
    2231             : GEN
    2232           0 : RgXQX_powers(GEN P, long n, GEN T)
    2233             : {
    2234           0 :   GEN v = cgetg(n+2, t_VEC);
    2235             :   long i;
    2236           0 :   gel(v, 1) = pol_1(varn(T));
    2237           0 :   if (n==0) return v;
    2238           0 :   gel(v, 2) = gcopy(P);
    2239           0 :   for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
    2240           0 :   return v;
    2241             : }
    2242             : 
    2243             : static GEN
    2244      598557 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
    2245             : static GEN
    2246           0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
    2247             : static GEN
    2248       71702 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2249             : static GEN
    2250       86982 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
    2251             : static GEN
    2252      268928 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2253             : static GEN
    2254      997491 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2255             : static GEN
    2256      972308 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2257             : static GEN
    2258         826 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2259             : static GEN
    2260      698180 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
    2261             : 
    2262             : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
    2263             :               _mul, _sqr, _one, _zero };
    2264             : 
    2265             : GEN
    2266           0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
    2267             : {
    2268           0 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
    2269             : }
    2270             : 
    2271             : GEN
    2272      398784 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2273             : {
    2274      398784 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2275      398783 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
    2276             : }
    2277             : 
    2278             : /* mod X^n */
    2279             : struct modXn {
    2280             :   long v; /* varn(X) */
    2281             :   long n;
    2282             : } ;
    2283             : static GEN
    2284       11893 : _sqrXn(void *data, GEN x) {
    2285       11893 :   struct modXn *S = (struct modXn*)data;
    2286       11893 :   return RgXn_sqr(x, S->n);
    2287             : }
    2288             : static GEN
    2289        4528 : _mulXn(void *data, GEN x, GEN y) {
    2290        4528 :   struct modXn *S = (struct modXn*)data;
    2291        4528 :   return RgXn_mul(x,y, S->n);
    2292             : }
    2293             : static GEN
    2294        1939 : _oneXn(void *data) {
    2295        1939 :   struct modXn *S = (struct modXn*)data;
    2296        1939 :   return pol_1(S->v);
    2297             : }
    2298             : static GEN
    2299           0 : _zeroXn(void *data) {
    2300           0 :   struct modXn *S = (struct modXn*)data;
    2301           0 :   return pol_0(S->v);
    2302             : }
    2303             : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
    2304             :                                           _oneXn, _zeroXn };
    2305             : 
    2306             : GEN
    2307         483 : RgXn_powers(GEN x, long m, long n)
    2308             : {
    2309         483 :   long d = degpol(x);
    2310         483 :   int use_sqr = (d<<1) >= n;
    2311             :   struct modXn S;
    2312         483 :   S.v = varn(x); S.n = n;
    2313         483 :   return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
    2314             : }
    2315             : 
    2316             : GEN
    2317        2286 : RgXn_powu_i(GEN x, ulong m, long n)
    2318             : {
    2319             :   struct modXn S;
    2320             :   long v;
    2321        2286 :   if (n == 0) return x;
    2322        2286 :   v = RgX_valrem(x, &x);
    2323        2286 :   if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
    2324        2265 :   S.v = varn(x); S.n = n;
    2325        2265 :   x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2326        2265 :   if (v) x = RgX_shift_shallow(x, m * v);
    2327        2265 :   return x;
    2328             : }
    2329             : GEN
    2330           0 : RgXn_powu(GEN x, ulong m, long n)
    2331             : {
    2332             :   pari_sp av;
    2333           0 :   if (n == 0) return gcopy(x);
    2334           0 :   av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
    2335             : }
    2336             : 
    2337             : GEN
    2338         966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
    2339             : {
    2340             :   struct modXn S;
    2341         966 :   S.v = varn(gel(x,2)); S.n = n;
    2342         966 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
    2343             : }
    2344             : 
    2345             : GEN
    2346           0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
    2347             : {
    2348           0 :   int use_sqr = 2*degpol(x) >= n;
    2349             :   struct modXn S;
    2350           0 :   S.v = varn(x); S.n = n;
    2351           0 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2352             : }
    2353             : 
    2354             : /* Q(x) mod t^n, x in R[t], n >= 1 */
    2355             : GEN
    2356        4067 : RgXn_eval(GEN Q, GEN x, long n)
    2357             : {
    2358        4067 :   long d = degpol(x);
    2359             :   int use_sqr;
    2360             :   struct modXn S;
    2361        4067 :   if (d == 1 && isrationalzero(gel(x,2)))
    2362             :   {
    2363        4060 :     GEN y = RgX_unscale(Q, gel(x,3));
    2364        4060 :     setvarn(y, varn(x)); return y;
    2365             :   }
    2366           7 :   S.v = varn(x);
    2367           7 :   S.n = n;
    2368           7 :   use_sqr = (d<<1) >= n;
    2369           7 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2370             : }
    2371             : 
    2372             : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
    2373             : static GEN
    2374     2351682 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2375             : {
    2376     2351682 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2377     2351682 :   return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
    2378             : }
    2379             : 
    2380             : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
    2381             : static GEN
    2382          14 : RgXn_sqrhigh(GEN f, long n2, long n)
    2383             : {
    2384          14 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2385          14 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2386             : }
    2387             : 
    2388             : static GEN
    2389     1365800 : RgXn_div_gen(GEN g, GEN f, long e)
    2390             : {
    2391             :   pari_sp av;
    2392             :   ulong mask;
    2393             :   GEN W, a;
    2394     1365800 :   long v = varn(f), n = 1;
    2395             : 
    2396     1365800 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2397     1365793 :   a = ginv(gel(f,2));
    2398     1365772 :   if (e == 1 && !g) return scalarpol(a, v);
    2399     1361992 :   else if (e == 2 && !g)
    2400             :   {
    2401             :     GEN b;
    2402      152693 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2403      129168 :     b = gneg(b);
    2404      129172 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2405      129173 :     return deg1pol(b, a, v);
    2406             :   }
    2407     1209299 :   av = avma;
    2408     1209299 :   W = scalarpol_shallow(a,v);
    2409     1209300 :   mask = quadratic_prec_mask(e);
    2410     3555116 :   while (mask > 1)
    2411             :   {
    2412             :     GEN u, fr;
    2413     2345816 :     long n2 = n;
    2414     2345816 :     n<<=1; if (mask & 1) n--;
    2415     2345816 :     mask >>= 1;
    2416     2345816 :     fr = RgXn_red_shallow(f, n);
    2417     2345816 :     if (mask>1 || !g)
    2418             :     {
    2419     1189209 :       u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2420     1189209 :       W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2421             :     }
    2422             :     else
    2423             :     {
    2424     1156607 :       GEN y = RgXn_mul(g, W, n), yt =  RgXn_red_shallow(y, n-n2);
    2425     1156607 :       u = RgXn_mul(yt, RgXn_mulhigh(fr,  W, n2, n), n-n2);
    2426     1156607 :       W = RgX_sub(y, RgX_shift_shallow(u, n2));
    2427             :     }
    2428     2345816 :     if (gc_needed(av,2))
    2429             :     {
    2430           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2431           0 :       W = gerepileupto(av, W);
    2432             :     }
    2433             :   }
    2434     1209300 :   return W;
    2435             : }
    2436             : 
    2437             : static GEN
    2438         119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
    2439             : {
    2440             :   GEN r;
    2441         119 :   if (lgefint(p) == 3)
    2442             :   {
    2443         119 :     ulong pp = uel(p, 2);
    2444         119 :     if (pp == 2)
    2445          14 :       r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
    2446             :     else
    2447         105 :       r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
    2448             :   }
    2449             :   else
    2450           0 :     r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
    2451         119 :   return FpX_to_mod(r, p);
    2452             : }
    2453             : 
    2454             : static GEN
    2455           0 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
    2456             : {
    2457           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2458           0 :   if (signe(T) == 0) pari_err_OP("/", x, y);
    2459           0 :   r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
    2460           0 :   return FpXQX_to_mod(r, T, p);
    2461             : }
    2462             : 
    2463             : static GEN
    2464          56 : RgXn_inv_FpX(GEN x, long e, GEN p)
    2465             : {
    2466             :   GEN r;
    2467          56 :   if (lgefint(p) == 3)
    2468             :   {
    2469          56 :     ulong pp = uel(p, 2);
    2470          56 :     if (pp == 2)
    2471          14 :       r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
    2472             :     else
    2473          42 :       r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
    2474             :   }
    2475             :   else
    2476           0 :     r = FpXn_inv(RgX_to_FpX(x, p), e, p);
    2477          56 :   return FpX_to_mod(r, p);
    2478             : }
    2479             : 
    2480             : static GEN
    2481           0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
    2482             : {
    2483           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2484           0 :   if (signe(T) == 0) pari_err_OP("/", gen_1, x);
    2485           0 :   r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
    2486           0 :   return FpXQX_to_mod(r, T, p);
    2487             : }
    2488             : 
    2489             : #define code(t1,t2) ((t1 << 6) | t2)
    2490             : 
    2491             : static GEN
    2492      209245 : RgXn_inv_fast(GEN x, long e)
    2493             : {
    2494             :   GEN p, pol;
    2495             :   long pa;
    2496      209245 :   long t = RgX_type(x,&p,&pol,&pa);
    2497      209251 :   switch(t)
    2498             :   {
    2499          56 :     case t_INTMOD: return RgXn_inv_FpX(x, e, p);
    2500           1 :     case code(t_POLMOD, t_INTMOD):
    2501           1 :                    return RgXn_inv_FpXQX(x, e, pol, p);
    2502      209194 :     default:       return NULL;
    2503             :   }
    2504             : }
    2505             : 
    2506             : static GEN
    2507     1156726 : RgXn_div_fast(GEN x, GEN y, long e)
    2508             : {
    2509             :   GEN p, pol;
    2510             :   long pa;
    2511     1156726 :   long t = RgX_type2(x,y,&p,&pol,&pa);
    2512     1156726 :   switch(t)
    2513             :   {
    2514         119 :     case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
    2515           0 :     case code(t_POLMOD, t_INTMOD):
    2516           0 :                    return RgXn_div_FpXQX(x, y, e, pol, p);
    2517     1156607 :     default:       return NULL;
    2518             :   }
    2519             : }
    2520             : #undef code
    2521             : 
    2522             : GEN
    2523     1156726 : RgXn_div_i(GEN g, GEN f, long e)
    2524             : {
    2525     1156726 :   GEN h = RgXn_div_fast(g, f, e);
    2526     1156726 :   if (h) return h;
    2527     1156607 :   return RgXn_div_gen(g, f, e);
    2528             : }
    2529             : 
    2530             : GEN
    2531      478259 : RgXn_div(GEN g, GEN f, long e)
    2532             : {
    2533      478259 :   pari_sp av = avma;
    2534      478259 :   return gerepileupto(av, RgXn_div_i(g, f, e));
    2535             : }
    2536             : 
    2537             : GEN
    2538      209244 : RgXn_inv_i(GEN f, long e)
    2539             : {
    2540      209244 :   GEN h = RgXn_inv_fast(f, e);
    2541      209249 :   if (h) return h;
    2542      209193 :   return RgXn_div_gen(NULL, f, e);
    2543             : }
    2544             : 
    2545             : GEN
    2546      119496 : RgXn_inv(GEN f, long e)
    2547             : {
    2548      119496 :   pari_sp av = avma;
    2549      119496 :   return gerepileupto(av, RgXn_inv_i(f, e));
    2550             : }
    2551             : 
    2552             : /* Compute intformal(x^n*S)/x^(n+1) */
    2553             : static GEN
    2554       58715 : RgX_integXn(GEN x, long n)
    2555             : {
    2556       58715 :   long i, lx = lg(x);
    2557             :   GEN y;
    2558       58715 :   if (lx == 2) return RgX_copy(x);
    2559       52682 :   y = cgetg(lx, t_POL); y[1] = x[1];
    2560      110479 :   for (i=2; i<lx; i++)
    2561       57798 :     gel(y,i) = gdivgs(gel(x,i), n+i-1);
    2562       52681 :   return RgX_renormalize_lg(y, lx);;
    2563             : }
    2564             : 
    2565             : GEN
    2566       52855 : RgXn_expint(GEN h, long e)
    2567             : {
    2568       52855 :   pari_sp av = avma, av2;
    2569       52855 :   long v = varn(h), n;
    2570       52855 :   GEN f = pol_1(v), g;
    2571             :   ulong mask;
    2572             : 
    2573       52855 :   if (!signe(h)) return f;
    2574       52855 :   g = pol_1(v);
    2575       52855 :   n = 1; mask = quadratic_prec_mask(e);
    2576       52855 :   av2 = avma;
    2577       58714 :   for (;mask>1;)
    2578             :   {
    2579             :     GEN u, w;
    2580       58714 :     long n2 = n;
    2581       58714 :     n<<=1; if (mask & 1) n--;
    2582       58714 :     mask >>= 1;
    2583       58714 :     u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
    2584       58713 :     u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
    2585       58715 :     w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
    2586       58711 :     f = RgX_add(f, RgX_shift_shallow(w, n2));
    2587       58709 :     if (mask<=1) break;
    2588        5859 :     u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
    2589        5859 :     g = RgX_sub(g, RgX_shift_shallow(u, n2));
    2590        5859 :     if (gc_needed(av2,2))
    2591             :     {
    2592           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
    2593           0 :       gerepileall(av2, 2, &f, &g);
    2594             :     }
    2595             :   }
    2596       52850 :   return gerepileupto(av, f);
    2597             : }
    2598             : 
    2599             : GEN
    2600           0 : RgXn_exp(GEN h, long e)
    2601             : {
    2602           0 :   long d = degpol(h);
    2603           0 :   if (d < 0) return pol_1(varn(h));
    2604           0 :   if (!d || !gequal0(gel(h,2)))
    2605           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2606           0 :   return RgXn_expint(RgX_deriv(h), e);
    2607             : }
    2608             : 
    2609             : GEN
    2610         154 : RgXn_reverse(GEN f, long e)
    2611             : {
    2612         154 :   pari_sp av = avma, av2;
    2613             :   ulong mask;
    2614             :   GEN fi, a, df, W, an;
    2615         154 :   long v = varn(f), n=1;
    2616         154 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2617           0 :     pari_err_INV("serreverse",f);
    2618         154 :   fi = ginv(gel(f,3));
    2619         154 :   a = deg1pol_shallow(fi,gen_0,v);
    2620         154 :   if (e <= 2) return gerepilecopy(av, a);
    2621         133 :   W = scalarpol(fi,v);
    2622         133 :   df = RgX_deriv(f);
    2623         133 :   mask = quadratic_prec_mask(e);
    2624         133 :   av2 = avma;
    2625         616 :   for (;mask>1;)
    2626             :   {
    2627             :     GEN u, fa, fr;
    2628         483 :     long n2 = n, rt;
    2629         483 :     n<<=1; if (mask & 1) n--;
    2630         483 :     mask >>= 1;
    2631         483 :     fr = RgXn_red_shallow(f, n);
    2632         483 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2633         483 :     an = RgXn_powers(a, rt, n);
    2634         483 :     if (n>1)
    2635             :     {
    2636         483 :       long n4 = (n2+1)>>1;
    2637         483 :       GEN dfr = RgXn_red_shallow(df, n2);
    2638         483 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2639         483 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2640         483 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2641             :     }
    2642         483 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2643         483 :     fa = RgX_shift(fa, -n2);
    2644         483 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2645         483 :     if (gc_needed(av2,2))
    2646             :     {
    2647           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2648           0 :       gerepileall(av2, 2, &a, &W);
    2649             :     }
    2650             :   }
    2651         133 :   return gerepileupto(av, a);
    2652             : }
    2653             : 
    2654             : GEN
    2655           7 : RgXn_sqrt(GEN h, long e)
    2656             : {
    2657           7 :   pari_sp av = avma, av2;
    2658           7 :   long v = varn(h), n = 1;
    2659           7 :   GEN f = scalarpol(gen_1, v), df = f;
    2660           7 :   ulong mask = quadratic_prec_mask(e);
    2661           7 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2662           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2663           7 :   av2 = avma;
    2664             :   while(1)
    2665           7 :   {
    2666          14 :     long n2 = n, m;
    2667             :     GEN g;
    2668          14 :     n<<=1; if (mask & 1) n--;
    2669          14 :     mask >>= 1;
    2670          14 :     m = n-n2;
    2671          14 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2672          14 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2673          14 :     if (mask==1) return gerepileupto(av, f);
    2674           7 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2675           7 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2676           7 :     if (gc_needed(av2,2))
    2677             :     {
    2678           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2679           0 :       gerepileall(av2, 2, &f, &df);
    2680             :     }
    2681             :   }
    2682             : }
    2683             : 
    2684             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2685             : GEN
    2686      114025 : RgXQ_powu(GEN x, ulong n, GEN T)
    2687             : {
    2688      114025 :   pari_sp av = avma;
    2689             : 
    2690      114025 :   if (!n) return pol_1(varn(x));
    2691       66326 :   if (n == 1) return RgX_copy(x);
    2692       18312 :   x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
    2693       18312 :   return gerepilecopy(av, x);
    2694             : }
    2695             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2696             : GEN
    2697      109018 : RgXQ_pow(GEN x, GEN n, GEN T)
    2698             : {
    2699             :   pari_sp av;
    2700      109018 :   long s = signe(n);
    2701             : 
    2702      109018 :   if (!s) return pol_1(varn(x));
    2703      109018 :   if (is_pm1(n) == 1)
    2704       86982 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2705       22036 :   av = avma;
    2706       22036 :   if (s < 0) x = RgXQ_inv(x, T);
    2707       22036 :   x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
    2708       22036 :   return gerepilecopy(av, x);
    2709             : }
    2710             : static GEN
    2711      188601 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
    2712             : static GEN
    2713      102843 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
    2714             : 
    2715             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2716             : GEN
    2717       10514 : ZXQ_powers(GEN x, long l, GEN T)
    2718             : {
    2719       10514 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2720       10514 :   return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
    2721             : }
    2722             : 
    2723             : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
    2724             : GEN
    2725      163346 : ZXQ_powu(GEN x, ulong n, GEN T)
    2726             : {
    2727      163346 :   pari_sp av = avma;
    2728             : 
    2729      163346 :   if (!n) return pol_1(varn(x));
    2730      163346 :   if (n == 1) return ZX_copy(x);
    2731      108804 :   x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
    2732      108803 :   return gerepilecopy(av, x);
    2733             : }
    2734             : 
    2735             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2736             : GEN
    2737        4200 : RgXQ_powers(GEN x, long l, GEN T)
    2738             : {
    2739        4200 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2740        4200 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2741             : }
    2742             : 
    2743             : GEN
    2744        6895 : RgXQV_factorback(GEN L, GEN e, GEN T)
    2745             : {
    2746        6895 :   return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
    2747             : }
    2748             : 
    2749             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2750             : GEN
    2751        7945 : QXQ_powers(GEN a, long n, GEN T)
    2752             : {
    2753             :   GEN den, v;
    2754        7945 :   if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
    2755        7931 :   v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
    2756             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2757        7931 :   if (den)
    2758             :   { /* restore denominators */
    2759        5033 :     GEN d = den;
    2760             :     long i;
    2761        5033 :     gel(v,2) = a;
    2762       25627 :     for (i=3; i<=n+1; i++) {
    2763       20594 :       d = mulii(d,den);
    2764       20594 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2765             :     }
    2766             :   }
    2767        7931 :   return v;
    2768             : }
    2769             : 
    2770             : static GEN
    2771        2247 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2772             : {
    2773        2247 :   long l, i, m = 0;
    2774             :   GEN dz, z;
    2775        2247 :   GEN V = cgetg_copy(v, &l);
    2776        7924 :   for (i = imin; i < l; i++)
    2777             :   {
    2778        5677 :     GEN c = gel(v, i);
    2779        5677 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2780             :   }
    2781        2247 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2782        2478 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2783        7924 :   for (i = imin; i < l; i++)
    2784             :   {
    2785        5677 :     GEN c = gel(v,i);
    2786        5677 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2787        5677 :     gel(V,i) = c;
    2788             :   }
    2789        2247 :   return V;
    2790             : }
    2791             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2792             : GEN
    2793        2016 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2794        2016 : { return do_QXQ_eval(v, 1, a, T); }
    2795             : 
    2796             : GEN
    2797         231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
    2798         231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2799             : 
    2800             : GEN
    2801        2891 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2802             : {
    2803        2891 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2804             : }
    2805             : 
    2806             : GEN
    2807        5096 : RgXQ_norm(GEN x, GEN T)
    2808             : {
    2809             :   pari_sp av;
    2810        5096 :   long dx = degpol(x);
    2811             :   GEN L, y;
    2812        5096 :   if (degpol(T)==0) return gpowgs(x,0);
    2813        5089 :   av = avma; y = resultant(T, x);
    2814        5089 :   L = leading_coeff(T);
    2815        5089 :   if (gequal1(L) || !signe(x)) return y;
    2816           0 :   return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
    2817             : }
    2818             : 
    2819             : GEN
    2820         413 : RgXQ_trace(GEN x, GEN T)
    2821             : {
    2822         413 :   pari_sp av = avma;
    2823             :   GEN dT, z;
    2824             :   long n;
    2825         413 :   if (degpol(T)==0) return gmulgs(x,0);
    2826         406 :   dT = RgX_deriv(T); n = degpol(dT);
    2827         406 :   z = RgXQ_mul(x, dT, T);
    2828         406 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2829         357 :   return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
    2830             : }
    2831             : 
    2832             : GEN
    2833     2498732 : RgX_blocks(GEN P, long n, long m)
    2834             : {
    2835     2498732 :   GEN z = cgetg(m+1,t_VEC);
    2836     2498732 :   long i,j, k=2, l = lg(P);
    2837     7689242 :   for(i=1; i<=m; i++)
    2838             :   {
    2839     5190510 :     GEN zi = cgetg(n+2,t_POL);
    2840     5190510 :     zi[1] = P[1];
    2841     5190510 :     gel(z,i) = zi;
    2842    15456560 :     for(j=2; j<n+2; j++)
    2843    10266050 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2844     5190510 :     zi = RgX_renormalize_lg(zi, n+2);
    2845             :   }
    2846     2498732 :   return z;
    2847             : }
    2848             : 
    2849             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2850             : void
    2851    29360619 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2852             : {
    2853    29360619 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2854             :   GEN p0, p1;
    2855             : 
    2856    29360105 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2857             : 
    2858    29267172 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2859    29267172 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2860    29268948 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2861    95257524 :   for (i=0; i<n1; i++)
    2862             :   {
    2863    65985997 :     p0[2+i] = p[2+(i<<1)];
    2864    65985997 :     p1[2+i] = p[3+(i<<1)];
    2865             :   }
    2866    29271527 :   if (n1 != n0)
    2867    13812427 :     p0[2+i] = p[2+(i<<1)];
    2868    29271527 :   *pe = normalizepol(p0);
    2869    29264027 :   *po = normalizepol(p1);
    2870             : }
    2871             : 
    2872             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2873             : GEN
    2874       43442 : RgX_splitting(GEN p, long k)
    2875             : {
    2876       43442 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2877             :   GEN r;
    2878             : 
    2879       43442 :   m = n/k;
    2880       43442 :   r = cgetg(k+1,t_VEC);
    2881      233226 :   for(i=1; i<=k; i++)
    2882             :   {
    2883      189784 :     gel(r,i) = cgetg(m+3, t_POL);
    2884      189784 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2885             :   }
    2886      570101 :   for (j=1, i=0, l=2; i<=n; i++)
    2887             :   {
    2888      526659 :     gmael(r,j,l) = gel(p,2+i);
    2889      526659 :     if (j==k) { j=1; l++; } else j++;
    2890             :   }
    2891      233226 :   for(i=1; i<=k; i++)
    2892      189784 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2893       43442 :   return r;
    2894             : }
    2895             : 
    2896             : /*******************************************************************/
    2897             : /*                                                                 */
    2898             : /*                        Kronecker form                           */
    2899             : /*                                                                 */
    2900             : /*******************************************************************/
    2901             : 
    2902             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2903             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2904             :  * normalized coefficients */
    2905             : GEN
    2906      176711 : Kronecker_to_mod(GEN z, GEN T)
    2907             : {
    2908      176711 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2909      176711 :   GEN x, t = cgetg(N,t_POL);
    2910      176711 :   t[1] = T[1];
    2911      176711 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2912      176711 :   x[1] = z[1];
    2913      176711 :   T = RgX_copy(T);
    2914     1311938 :   for (i=2; i<lx+2; i++, z+= N-2)
    2915             :   {
    2916     5426036 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2917     1135227 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2918             :   }
    2919      176711 :   N = (l-2) % (N-2) + 2;
    2920      455603 :   for (j=2; j<N; j++) t[j] = z[j];
    2921      176711 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2922      176711 :   return normalizepol_lg(x, i+1);
    2923             : }
    2924             : 
    2925             : /*******************************************************************/
    2926             : /*                                                                 */
    2927             : /*                        Domain detection                         */
    2928             : /*                                                                 */
    2929             : /*******************************************************************/
    2930             : 
    2931             : static GEN
    2932      151593 : zero_FpX_mod(GEN p, long v)
    2933             : {
    2934      151593 :   GEN r = cgetg(3,t_POL);
    2935      151593 :   r[1] = evalvarn(v);
    2936      151593 :   gel(r,2) = mkintmod(gen_0, icopy(p));
    2937      151593 :   return r;
    2938             : }
    2939             : 
    2940             : static GEN
    2941      500611 : RgX_mul_FpX(GEN x, GEN y, GEN p)
    2942             : {
    2943      500611 :   pari_sp av = avma;
    2944             :   GEN r;
    2945      500611 :   if (lgefint(p) == 3)
    2946             :   {
    2947      470479 :     ulong pp = uel(p, 2);
    2948      470479 :     r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
    2949             :                                   RgX_to_Flx(y, pp), pp));
    2950             :   }
    2951             :   else
    2952       30132 :     r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2953      500611 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2954      386727 :   return gerepileupto(av, FpX_to_mod(r, p));
    2955             : }
    2956             : 
    2957             : static GEN
    2958         532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
    2959             : {
    2960         532 :   GEN r = cgetg(3,t_POL);
    2961         532 :   r[1] = evalvarn(v);
    2962         532 :   gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
    2963         532 :   return r;
    2964             : }
    2965             : 
    2966             : static GEN
    2967        2037 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2968             : {
    2969        2037 :   pari_sp av = avma;
    2970             :   long dT;
    2971             :   GEN kx, ky, r;
    2972        2037 :   GEN T = RgX_to_FpX(pol, p);
    2973        2037 :   if (signe(T)==0) pari_err_OP("*", x, y);
    2974        2030 :   dT = degpol(T);
    2975        2030 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    2976        2030 :   ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
    2977        2030 :   r = FpX_mul(kx, ky, p);
    2978        2030 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    2979        1505 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    2980             : }
    2981             : 
    2982             : static GEN
    2983      410318 : RgX_liftred(GEN x, GEN T)
    2984      410318 : { return RgXQX_red(liftpol_shallow(x), T); }
    2985             : 
    2986             : static GEN
    2987      161530 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
    2988             : {
    2989      161530 :   pari_sp av = avma;
    2990      161530 :   long dT = degpol(T);
    2991      161530 :   GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
    2992             :                  RgXX_to_Kronecker(RgX_liftred(y, T), dT));
    2993      161530 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    2994             : }
    2995             : 
    2996             : static GEN
    2997        1358 : RgX_sqr_FpX(GEN x, GEN p)
    2998             : {
    2999        1358 :   pari_sp av = avma;
    3000             :   GEN r;
    3001        1358 :   if (lgefint(p) == 3)
    3002             :   {
    3003        1154 :     ulong pp = uel(p, 2);
    3004        1154 :     r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
    3005             :   }
    3006             :   else
    3007         204 :     r = FpX_sqr(RgX_to_FpX(x, p), p);
    3008        1358 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3009        1239 :   return gerepileupto(av, FpX_to_mod(r, p));
    3010             : }
    3011             : 
    3012             : static GEN
    3013         196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
    3014             : {
    3015         196 :   pari_sp av = avma;
    3016             :   long dT;
    3017         196 :   GEN kx, r, T = RgX_to_FpX(pol, p);
    3018         196 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3019         189 :   dT = degpol(T);
    3020         189 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    3021         189 :   r = FpX_sqr(kx, p);
    3022         189 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3023         189 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    3024             : }
    3025             : 
    3026             : static GEN
    3027       13487 : RgX_sqr_QXQX(GEN x, GEN T)
    3028             : {
    3029       13487 :   pari_sp av = avma;
    3030       13487 :   long dT = degpol(T);
    3031       13487 :   GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
    3032       13487 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    3033             : }
    3034             : 
    3035             : static GEN
    3036       67466 : RgX_rem_FpX(GEN x, GEN y, GEN p)
    3037             : {
    3038       67466 :   pari_sp av = avma;
    3039             :   GEN r;
    3040       67466 :   if (lgefint(p) == 3)
    3041             :   {
    3042       50676 :     ulong pp = uel(p, 2);
    3043       50676 :     r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
    3044             :                                   RgX_to_Flx(y, pp), pp));
    3045             :   }
    3046             :   else
    3047       16790 :     r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3048       67466 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3049       29876 :   return gerepileupto(av, FpX_to_mod(r, p));
    3050             : }
    3051             : 
    3052             : static GEN
    3053       36885 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
    3054             : {
    3055       36885 :   pari_sp av = avma;
    3056             :   GEN r;
    3057       36885 :   r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3058       36885 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3059             : }
    3060             : static GEN
    3061          70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3062             : {
    3063          70 :   pari_sp av = avma;
    3064             :   GEN r;
    3065          70 :   GEN T = RgX_to_FpX(pol, p);
    3066          70 :   if (signe(T) == 0) pari_err_OP("%", x, y);
    3067          63 :   if (lgefint(p) == 3)
    3068             :   {
    3069          55 :     ulong pp = uel(p, 2);
    3070          55 :     GEN Tp = ZX_to_Flx(T, pp);
    3071          55 :     r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
    3072             :                               RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3073             :   }
    3074             :   else
    3075           8 :     r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3076          63 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3077          56 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3078             : }
    3079             : 
    3080             : #define code(t1,t2) ((t1 << 6) | t2)
    3081             : static GEN
    3082    68253444 : RgX_mul_fast(GEN x, GEN y)
    3083             : {
    3084             :   GEN p, pol;
    3085             :   long pa;
    3086    68253444 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3087    68253639 :   switch(t)
    3088             :   {
    3089    37889342 :     case t_INT:    return ZX_mul(x,y);
    3090     2771146 :     case t_FRAC:   return QX_mul(x,y);
    3091      104665 :     case t_FFELT:  return FFX_mul(x, y, pol);
    3092      500611 :     case t_INTMOD: return RgX_mul_FpX(x, y, p);
    3093      161530 :     case code(t_POLMOD, t_INT):
    3094             :     case code(t_POLMOD, t_FRAC):
    3095      161530 :                    return RgX_mul_QXQX(x, y, pol);
    3096        2029 :     case code(t_POLMOD, t_INTMOD):
    3097        2029 :                    return RgX_mul_FpXQX(x, y, pol, p);
    3098    26824316 :     default:       return NULL;
    3099             :   }
    3100             : }
    3101             : static GEN
    3102      792560 : RgX_sqr_fast(GEN x)
    3103             : {
    3104             :   GEN p, pol;
    3105             :   long pa;
    3106      792560 :   long t = RgX_type(x,&p,&pol,&pa);
    3107      792560 :   switch(t)
    3108             :   {
    3109      751122 :     case t_INT:    return ZX_sqr(x);
    3110       16668 :     case t_FRAC:   return QX_sqr(x);
    3111        2296 :     case t_FFELT:  return FFX_sqr(x, pol);
    3112        1358 :     case t_INTMOD: return RgX_sqr_FpX(x, p);
    3113       13487 :     case code(t_POLMOD, t_INT):
    3114             :     case code(t_POLMOD, t_FRAC):
    3115       13487 :                    return RgX_sqr_QXQX(x, pol);
    3116         196 :     case code(t_POLMOD, t_INTMOD):
    3117         196 :                    return RgX_sqr_FpXQX(x, pol, p);
    3118        7433 :     default:       return NULL;
    3119             :   }
    3120             : }
    3121             : 
    3122             : static GEN
    3123    10452618 : RgX_rem_fast(GEN x, GEN y)
    3124             : {
    3125             :   GEN p, pol;
    3126             :   long pa;
    3127    10452618 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3128    10452850 :   switch(t)
    3129             :   {
    3130     3215137 :     case t_INT:    return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
    3131     1659505 :     case t_FRAC:   return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
    3132          84 :     case t_FFELT:  return FFX_rem(x, y, pol);
    3133       67466 :     case t_INTMOD: return RgX_rem_FpX(x, y, p);
    3134       36885 :     case code(t_POLMOD, t_INT):
    3135             :     case code(t_POLMOD, t_FRAC):
    3136       36885 :                    return RgX_rem_QXQX(x, y, pol);
    3137          82 :     case code(t_POLMOD, t_INTMOD):
    3138          82 :                    return RgX_rem_FpXQX(x, y, pol, p);
    3139     5473691 :     default:       return NULL;
    3140             :   }
    3141             : }
    3142             : 
    3143             : #undef code
    3144             : 
    3145             : GEN
    3146    55617623 : RgX_mul(GEN x, GEN y)
    3147             : {
    3148    55617623 :   GEN z = RgX_mul_fast(x,y);
    3149    55617678 :   if (!z) z = RgX_mul_i(x,y);
    3150    55617028 :   return z;
    3151             : }
    3152             : 
    3153             : GEN
    3154      780408 : RgX_sqr(GEN x)
    3155             : {
    3156      780408 :   GEN z = RgX_sqr_fast(x);
    3157      780402 :   if (!z) z = RgX_sqr_i(x);
    3158      780402 :   return z;
    3159             : }
    3160             : 
    3161             : GEN
    3162    10452635 : RgX_rem(GEN x, GEN y)
    3163             : {
    3164    10452635 :   GEN z = RgX_rem_fast(x, y);
    3165    10452790 :   if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
    3166    10452669 :   return z;
    3167             : }
    3168             : 
    3169             : static GEN
    3170           0 : _RgX_mul(void* E, GEN x, GEN y)
    3171           0 : { (void) E; return RgX_mul(x, y); }
    3172             : 
    3173             : GEN
    3174           0 : RgXV_prod(GEN V)
    3175           0 : { return gen_product(V, NULL, &_RgX_mul); }
    3176             : 
    3177             : GEN
    3178    10225417 : RgXn_mul(GEN f, GEN g, long n)
    3179             : {
    3180    10225417 :   pari_sp av = avma;
    3181    10225417 :   GEN h = RgX_mul_fast(f,g);
    3182    10225419 :   if (!h) return RgXn_mul2(f,g,n);
    3183     1496406 :   if (degpol(h) < n) return h;
    3184      288819 :   return gerepilecopy(av, RgXn_red_shallow(h, n));
    3185             : }
    3186             : 
    3187             : GEN
    3188       12152 : RgXn_sqr(GEN f, long n)
    3189             : {
    3190       12152 :   pari_sp av = avma;
    3191       12152 :   GEN g = RgX_sqr_fast(f);
    3192       12152 :   if (!g) return RgXn_sqr2(f,n);
    3193       11837 :   if (degpol(g) < n) return g;
    3194       10640 :   return gerepilecopy(av, RgXn_red_shallow(g, n));
    3195             : }
    3196             : 
    3197             : /* (f*g) \/ x^n */
    3198             : GEN
    3199     2410409 : RgX_mulhigh_i(GEN f, GEN g, long n)
    3200             : {
    3201     2410409 :   GEN h = RgX_mul_fast(f,g);
    3202     2410409 :   return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
    3203             : }
    3204             : 
    3205             : /* (f*g) \/ x^n */
    3206             : GEN
    3207           0 : RgX_sqrhigh_i(GEN f, long n)
    3208             : {
    3209           0 :   GEN h = RgX_sqr_fast(f);
    3210           0 :   return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
    3211             : }

Generated by: LCOV version 1.13