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.18.1 lcov report (development 30166-834405808f) Lines: 1628 1783 91.3 %
Date: 2025-04-21 09:18:53 Functions: 204 222 91.9 %
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    19439423 : brent_kung_optpow(long d, long n, long m)
      32             : {
      33             :   long p, r;
      34    19439423 :   long pold=1, rold=n*(d-1);
      35    91147987 :   for(p=2; p<=d; p++)
      36             :   {
      37    71708564 :     r = m*(p-1) + n*((d-1)/p);
      38    71708564 :     if (r<rold) { pold=p; rold=r; }
      39             :   }
      40    19439423 :   return pold;
      41             : }
      42             : 
      43             : static GEN
      44    14562316 : 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    14562316 :   pari_sp av = avma;
      48             :   long i;
      49    14562316 :   GEN z = cmul(E,P,a,ff->one(E));
      50    14538322 :   if (!z) z = gen_0;
      51    76414094 :   for (i=1; i<=n; i++)
      52             :   {
      53    61876200 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      54    61900434 :     if (t) {
      55    46952361 :       z = ff->add(E, z, t);
      56    46917966 :       if (gc_needed(av,2)) z = gc_upto(av, z);
      57             :     }
      58             :   }
      59    14537894 :   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    12468657 : 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    12468657 :   pari_sp av = avma;
      72    12468657 :   long l = lg(V)-1;
      73             :   GEN z, u;
      74             : 
      75    12468657 :   if (d < 0) return ff->zero(E);
      76    11860882 :   if (d < l) return gc_upto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      77     1157664 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      78     1157664 :   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     1157664 :   d -= l;
      84     1157664 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      85     2700672 :   while (d >= l-1)
      86             :   {
      87     1541857 :     d -= l-1;
      88     1541857 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      89     1541836 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      90     1541819 :     if (gc_needed(av,2))
      91          91 :       z = gc_upto(av, z);
      92             :   }
      93     1158815 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      94     1158819 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      95     1158826 :   return gc_upto(av, ff->red(E,z));
      96             : }
      97             : 
      98             : GEN
      99      873988 : 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      873988 :   pari_sp av = avma;
     103             :   GEN z, V;
     104             :   long rtd;
     105      873988 :   if (d < 0) return ff->zero(E);
     106      872693 :   rtd = (long) sqrt((double)d);
     107      872693 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     108      872710 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     109      872710 :   return gc_upto(av, z);
     110             : }
     111             : 
     112             : static GEN
     113     2027824 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     114             : static GEN
     115    19247437 : _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     1839969 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     120             : static GEN
     121      591667 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     122             : static GEN
     123     2068835 : _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      668656 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
     161             : {
     162      668656 :   pari_sp av = avma, btop;
     163      668656 :   long i, d = degpol(P), o;
     164             :   GEN s;
     165      668653 :   if (signe(P)==0) return pol_0(varn(P));
     166      668653 :   s = gel(P, d+2);
     167      668653 :   if (d == 0) return gcopy(s);
     168      665377 :   o = RgX_deflate_order(P);
     169      665393 :   if (o > 1) A = gpowgs(A, o);
     170      665405 :   btop = avma;
     171     2313533 :   for (i = d-o; i >= 0; i-=o)
     172             :   {
     173     1648139 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
     174     1648128 :     if (gc_needed(btop,1))
     175             :     {
     176          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
     177          13 :       s = gc_upto(btop, s);
     178             :     }
     179             :   }
     180      665394 :   return gc_upto(av, s);
     181             : }
     182             : 
     183             : GEN
     184        1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
     185             : {
     186        1652 :   pari_sp av = avma;
     187        1652 :   long i, d = degpol(P), v = varn(A);
     188             :   GEN s;
     189        1652 :   if (signe(P)==0) return pol_0(v);
     190        1652 :   if (d == 0) return scalarpol(gel(P, d+2), v);
     191        1232 :   s = scalarpol_shallow(gel(P, d+2), v);
     192        4963 :   for (i = d-1; i >= 0; i--)
     193             :   {
     194        3731 :     GEN c = gel(P,i+2), b = gel(B,d+1-i);
     195        3731 :     s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
     196        3731 :     if (gc_needed(av,1))
     197             :     {
     198           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
     199           0 :       s = gc_upto(av, s);
     200             :     }
     201             :   }
     202        1232 :   return gc_upto(av, s);
     203             : }
     204             : 
     205             : const struct bb_algebra *
     206      281056 : get_Rg_algebra(void)
     207             : {
     208      281056 :   return &Rg_algebra;
     209             : }
     210             : 
     211             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     212             : 
     213             : static GEN
     214       11333 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     215             : {
     216             :   (void) E;
     217       11333 :   return RgX_divrem(x, y, r);
     218             : }
     219             : 
     220             : GEN
     221        3157 : RgX_digits(GEN x, GEN T)
     222             : {
     223        3157 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     224        3157 :   if (signe(x)==0) return(cgetg(1, t_VEC));
     225        3157 :   return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     226             : }
     227             : 
     228             : /*******************************************************************/
     229             : /*                                                                 */
     230             : /*                         RgX                                     */
     231             : /*                                                                 */
     232             : /*******************************************************************/
     233             : 
     234             : long
     235    24838206 : RgX_equal(GEN x, GEN y)
     236             : {
     237    24838206 :   long i = lg(x);
     238             : 
     239    24838206 :   if (i != lg(y)) return 0;
     240   106915591 :   for (i--; i > 1; i--)
     241    82405340 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     242    24510251 :   return 1;
     243             : }
     244             : 
     245             : /* Returns 1 in the base ring over which x is defined */
     246             : /* HACK: this also works for t_SER */
     247             : GEN
     248   156946165 : Rg_get_1(GEN x)
     249             : {
     250             :   GEN p, T;
     251   156946165 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     252   156946166 :   if (RgX_type_is_composite(tx))
     253    11657560 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     254   156946166 :   switch(tx)
     255             :   {
     256      427406 :     case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
     257         945 :     case t_PADIC: return cvtop(gen_1, p, lx);
     258        5950 :     case t_FFELT: return FF_1(T);
     259   156511865 :     default: return gen_1;
     260             :   }
     261             : }
     262             : /* Returns 0 in the base ring over which x is defined */
     263             : /* HACK: this also works for t_SER */
     264             : GEN
     265     6942912 : Rg_get_0(GEN x)
     266             : {
     267             :   GEN p, T;
     268     6942912 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     269     6942912 :   if (RgX_type_is_composite(tx))
     270       61334 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     271     6942912 :   switch(tx)
     272             :   {
     273         497 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     274          42 :     case t_PADIC: return zeropadic(p, lx);
     275         210 :     case t_FFELT: return FF_zero(T);
     276     6942163 :     default: return gen_0;
     277             :   }
     278             : }
     279             : 
     280             : GEN
     281        6874 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     282             : {
     283        6874 :   long i, n = degpol(P);
     284             :   GEN z, dz, dP;
     285        6874 :   if (n < 0) return gen_0;
     286        6874 :   P = Q_remove_denom(P, &dP);
     287        6874 :   z = gel(P,2); if (n == 0) return icopy(z);
     288        3878 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     289        3878 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     290        7483 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     291        3878 :   dz = mul_denom(dP, dV);
     292        3878 :   return dz? RgX_Rg_div(z, dz): z;
     293             : }
     294             : 
     295             : /* Return P(h * x), not memory clean */
     296             : GEN
     297       51365 : RgX_unscale(GEN P, GEN h)
     298             : {
     299       51365 :   long i, l = lg(P);
     300       51365 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     301       51365 :   Q[1] = P[1];
     302       51365 :   if (l == 2) return Q;
     303       40123 :   gel(Q,2) = gcopy(gel(P,2));
     304       94488 :   for (i=3; i<l; i++)
     305             :   {
     306       54365 :     hi = gmul(hi,h);
     307       54365 :     gel(Q,i) = gmul(gel(P,i), hi);
     308             :   }
     309       40123 :   return Q;
     310             : }
     311             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     312             : GEN
     313     1467269 : ZX_z_unscale(GEN P, long h)
     314             : {
     315     1467269 :   long i, l = lg(P);
     316     1467269 :   GEN Q = cgetg(l, t_POL);
     317     1467269 :   Q[1] = P[1];
     318     1467269 :   if (l == 2) return Q;
     319     1378723 :   gel(Q,2) = gel(P,2);
     320     1378723 :   if (l == 3) return Q;
     321     1352417 :   if (h == -1)
     322      236970 :     for (i = 3; i < l; i++)
     323             :     {
     324      196387 :       gel(Q,i) = negi(gel(P,i));
     325      196387 :       if (++i == l) break;
     326      145384 :       gel(Q,i) = gel(P,i);
     327             :     }
     328             :   else
     329             :   {
     330             :     GEN hi;
     331     1260831 :     gel(Q,3) = mulis(gel(P,3), h);
     332     1260830 :     hi = sqrs(h);
     333     5415785 :     for (i = 4; i < l; i++)
     334             :     {
     335     4154953 :       gel(Q,i) = mulii(gel(P,i), hi);
     336     4154956 :       if (i != l-1) hi = mulis(hi,h);
     337             :     }
     338             :   }
     339     1352418 :   return Q;
     340             : }
     341             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     342             : GEN
     343     1013197 : ZX_unscale(GEN P, GEN h)
     344             : {
     345             :   long i, l;
     346             :   GEN Q, hi;
     347     1013197 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     348         881 :   l = lg(P); Q = cgetg(l, t_POL);
     349         881 :   Q[1] = P[1];
     350         881 :   if (l == 2) return Q;
     351         881 :   gel(Q,2) = gel(P,2);
     352         881 :   if (l == 3) return Q;
     353         881 :   hi = h;
     354         881 :   gel(Q,3) = mulii(gel(P,3), hi);
     355        2741 :   for (i = 4; i < l; i++)
     356             :   {
     357        1860 :     hi = mulii(hi,h);
     358        1860 :     gel(Q,i) = mulii(gel(P,i), hi);
     359             :   }
     360         881 :   return Q;
     361             : }
     362             : /* P a ZX. Return P(x << n), not memory clean */
     363             : GEN
     364      924326 : ZX_unscale2n(GEN P, long n)
     365             : {
     366      924326 :   long i, ni = n, l = lg(P);
     367      924326 :   GEN Q = cgetg(l, t_POL);
     368      924329 :   Q[1] = P[1];
     369      924329 :   if (l == 2) return Q;
     370      924329 :   gel(Q,2) = gel(P,2);
     371      924329 :   if (l == 3) return Q;
     372      924329 :   gel(Q,3) = shifti(gel(P,3), ni);
     373     3166624 :   for (i=4; i<l; i++)
     374             :   {
     375     2242339 :     ni += n;
     376     2242339 :     gel(Q,i) = shifti(gel(P,i), ni);
     377             :   }
     378      924285 :   return Q;
     379             : }
     380             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     381             : GEN
     382       12997 : ZX_unscale_div(GEN P, GEN h)
     383             : {
     384       12997 :   long i, l = lg(P);
     385       12997 :   GEN hi, Q = cgetg(l, t_POL);
     386       12997 :   Q[1] = P[1];
     387       12997 :   if (l == 2) return Q;
     388       12997 :   gel(Q,2) = diviiexact(gel(P,2), h);
     389       12997 :   if (l == 3) return Q;
     390       12997 :   gel(Q,3) = gel(P,3);
     391       12997 :   if (l == 4) return Q;
     392       12997 :   hi = h;
     393       12997 :   gel(Q,4) = mulii(gel(P,4), hi);
     394       64121 :   for (i=5; i<l; i++)
     395             :   {
     396       51124 :     hi = mulii(hi,h);
     397       51124 :     gel(Q,i) = mulii(gel(P,i), hi);
     398             :   }
     399       12997 :   return Q;
     400             : }
     401             : /* P(h*X) / h^k, assuming the result is a ZX */
     402             : GEN
     403        1393 : ZX_unscale_divpow(GEN P, GEN h, long k)
     404             : {
     405        1393 :   long i, j, l = lg(P);
     406        1393 :   GEN H, Q = cgetg(l, t_POL);
     407        1393 :   Q[1] = P[1]; if (l == 2) return Q;
     408        1393 :   H = gpowers(h, maxss(k, l - 3 - k));
     409        5572 :   for (i = 2, j = k+1; j > 1 && i < l; i++)
     410        4179 :     gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
     411        1393 :   if (i == l) return Q;
     412        1393 :   gel(Q, i) = gel(P, i); i++;
     413        5082 :   for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
     414        1393 :   return Q;
     415             : }
     416             : 
     417             : GEN
     418        6223 : RgXV_unscale(GEN x, GEN h)
     419             : {
     420        6223 :   if (isint1(h)) return gcopy(x);
     421       18182 :   pari_APPLY_same(RgX_unscale(gel(x,i), h));
     422             : }
     423             : 
     424             : /* Return h^degpol(P) P(x / h), not memory clean */
     425             : GEN
     426     4342224 : RgX_rescale(GEN P, GEN h)
     427             : {
     428     4342224 :   long i, l = lg(P);
     429     4342224 :   GEN Q = cgetg(l,t_POL), hi = h;
     430     4342221 :   gel(Q,l-1) = gel(P,l-1);
     431    11427189 :   for (i=l-2; i>=2; i--)
     432             :   {
     433    11424568 :     gel(Q,i) = gmul(gel(P,i), hi);
     434    11424411 :     if (i == 2) break;
     435     7084648 :     hi = gmul(hi,h);
     436             :   }
     437     4342384 :   Q[1] = P[1]; return Q;
     438             : }
     439             : 
     440             : GEN
     441        2401 : RgXV_rescale(GEN x, GEN h)
     442             : {
     443        2401 :   if (isint1(h)) return RgX_copy(x);
     444       16086 :   pari_APPLY_same(RgX_rescale(gel(x,i), h));
     445             : }
     446             : 
     447             : /* A(X^d) --> A(X) */
     448             : GEN
     449     1175180 : RgX_deflate(GEN x0, long d)
     450             : {
     451             :   GEN z, y, x;
     452     1175180 :   long i,id, dy, dx = degpol(x0);
     453     1175177 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     454      459089 :   dy = dx/d;
     455      459089 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     456      459089 :   z = y + 2;
     457      459089 :   x = x0+ 2;
     458     1762688 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     459      459089 :   return y;
     460             : }
     461             : 
     462             : GEN
     463        1260 : RgX_homogenize_deg(GEN P, long d, long v)
     464             : {
     465             :   long i, l;
     466        1260 :   GEN Q = cgetg_copy(P, &l);
     467        1260 :   Q[1] = P[1];
     468        4018 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
     469        1260 :   return Q;
     470             : }
     471             : 
     472             : GEN
     473       18732 : RgX_homogenize(GEN P, long v)
     474             : {
     475             :   long i, l, d;
     476       18732 :   GEN Q = cgetg_copy(P, &l);
     477       18732 :   Q[1] = P[1]; d = l-3;
     478      165998 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
     479       18732 :   return Q;
     480             : }
     481             : 
     482             : /* F a t_RFRAC */
     483             : long
     484         126 : rfrac_deflate_order(GEN F)
     485             : {
     486         126 :   GEN N = gel(F,1), D = gel(F,2);
     487         126 :   long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
     488         126 :   if (m == 1) return 1;
     489          42 :   if (typ(N) == t_POL && varn(N) == varn(D))
     490          28 :     m = cgcd(m, RgX_deflate_order(N));
     491          42 :   return m;
     492             : }
     493             : /* F a t_RFRAC */
     494             : GEN
     495         126 : rfrac_deflate_max(GEN F, long *m)
     496             : {
     497         126 :   *m = rfrac_deflate_order(F);
     498         126 :   return rfrac_deflate(F, *m);
     499             : }
     500             : /* F a t_RFRAC */
     501             : GEN
     502         126 : rfrac_deflate(GEN F, long m)
     503             : {
     504         126 :   GEN N = gel(F,1), D = gel(F,2);
     505         126 :   if (m == 1) return F;
     506          42 :   if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
     507          42 :   D = RgX_deflate(D, m); return mkrfrac(N, D);
     508             : }
     509             : 
     510             : /* return x0(X^d) */
     511             : GEN
     512      881773 : RgX_inflate(GEN x0, long d)
     513             : {
     514      881773 :   long i, id, dy, dx = degpol(x0);
     515      881773 :   GEN x = x0 + 2, z, y;
     516      881773 :   if (dx <= 0) return leafcopy(x0);
     517      804483 :   dy = dx*d;
     518      804483 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     519      804483 :   z = y + 2;
     520    28079645 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     521    11488640 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     522      804483 :   return y;
     523             : }
     524             : 
     525             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     526             : static GEN
     527     5701668 : RgX_translate_basecase(GEN P, GEN c)
     528             : {
     529     5701668 :   pari_sp av = avma;
     530             :   GEN Q, R;
     531             :   long i, k, n;
     532             : 
     533     5701668 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     534     5699740 :   Q = leafcopy(P);
     535     5699735 :   R = Q+2; n = degpol(P);
     536     5699734 :   if (isint1(c))
     537             :   {
     538        7595 :     for (i=1; i<=n; i++)
     539             :     {
     540       20265 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
     541        5390 :       if (gc_needed(av,2))
     542             :       {
     543           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
     544           0 :         Q = gc_GEN(av, Q); R = Q+2;
     545             :       }
     546             :     }
     547             :   }
     548     5697526 :   else if (isintm1(c))
     549             :   {
     550       15974 :     for (i=1; i<=n; i++)
     551             :     {
     552       49630 :       for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
     553       12299 :       if (gc_needed(av,2))
     554             :       {
     555           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
     556           0 :         Q = gc_GEN(av, Q); R = Q+2;
     557             :       }
     558             :     }
     559             :   }
     560             :   else
     561             :   {
     562    19424660 :     for (i=1; i<=n; i++)
     563             :     {
     564    45450449 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
     565    13730810 :       if (gc_needed(av,2))
     566             :       {
     567           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
     568           0 :         Q = gc_GEN(av, Q); R = Q+2;
     569             :       }
     570             :     }
     571             :   }
     572     5699466 :   return gc_GEN(av, Q);
     573             : }
     574             : 
     575             : static GEN
     576      102627 : zero_FpX_mod(GEN p, long v)
     577             : {
     578      102627 :   GEN r = cgetg(3,t_POL);
     579      102627 :   r[1] = evalvarn(v);
     580      102627 :   gel(r,2) = mkintmod(gen_0, icopy(p));
     581      102627 :   return r;
     582             : }
     583             : 
     584             : static GEN
     585           0 : RgX_translate_FpX(GEN P, GEN c, GEN p)
     586             : {
     587           0 :   pari_sp av = avma;
     588             :   GEN r;
     589             : #if 0
     590             :   if (lgefint(p) == 3)
     591             :   {
     592             :     ulong pp = uel(p, 2);
     593             :     r = Flx_to_ZX_inplace(Flx_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
     594             :   }
     595             :   else
     596             : #endif
     597           0 :     r = FpX_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
     598           0 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
     599           0 :   return gc_upto(av, FpX_to_mod(r, p));
     600             : }
     601             : 
     602             : static GEN
     603     6063361 : RgX_translate_fast(GEN P, GEN c)
     604             : {
     605             :   GEN p, pol;
     606             :   long pa;
     607     6063361 :   long t = RgX_Rg_type(P, c, &p,&pol,&pa);
     608     6063364 :   switch(t)
     609             :   {
     610      361884 :     case t_INT:    return ZX_translate(P, c);
     611           1 :     case t_INTMOD: return RgX_translate_FpX(P, c, p);
     612     5701479 :     default:       return NULL;
     613             :   }
     614             : }
     615             : 
     616             : static GEN
     617     5701857 : RgX_translate_i(GEN P, GEN c)
     618             : {
     619     5701857 :   pari_sp av = avma;
     620             :   long n;
     621     5701857 :   n = degpol(P);
     622     5701857 :   if (n < 40)
     623     5701668 :     return RgX_translate_basecase(P, c);
     624             :   else
     625             :   {
     626         189 :     long d = n >> 1;
     627         189 :     GEN Q = RgX_translate_i(RgX_shift_shallow(P, -d), c);
     628         189 :     GEN R = RgX_translate_i(RgXn_red_shallow(P, d), c);
     629         189 :     GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
     630         189 :     return gc_upto(av, RgX_add(RgX_mul(Q, S), R));
     631             :   }
     632             : }
     633             : 
     634             : GEN
     635     6063359 : RgX_translate(GEN P, GEN c)
     636             : {
     637     6063359 :   GEN R = RgX_translate_fast(P, c);
     638     6063363 :   return R ? R: RgX_translate_i(P,c);
     639             : }
     640             : /* P(ax + b) */
     641             : GEN
     642       30212 : RgX_affine(GEN P, GEN a, GEN b)
     643             : {
     644       30212 :   if (!gequal0(b)) P = RgX_translate(P, b);
     645       30212 :   return RgX_unscale(P, a);
     646             : }
     647             : 
     648             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     649             : GEN
     650       32464 : RgXQX_translate(GEN P, GEN c, GEN T)
     651             : {
     652       32464 :   pari_sp av = avma;
     653             :   GEN Q, R;
     654             :   long i, k, n;
     655             : 
     656       32464 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     657       32121 :   Q = leafcopy(P);
     658       32121 :   R = Q+2; n = degpol(P);
     659      103762 :   for (i=1; i<=n; i++)
     660             :   {
     661      301670 :     for (k=n-i; k<n; k++)
     662             :     {
     663      230029 :       pari_sp av2 = avma;
     664      230029 :       gel(R,k) = gc_upto(av2,
     665      230029 :                    RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
     666             :     }
     667       71641 :     if (gc_needed(av,2))
     668             :     {
     669           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     670           0 :       Q = gc_GEN(av, Q); R = Q+2;
     671             :     }
     672             :   }
     673       32121 :   return gc_GEN(av, Q);
     674             : }
     675             : 
     676             : /********************************************************************/
     677             : /**                                                                **/
     678             : /**                          CONVERSIONS                           **/
     679             : /**                       (not memory clean)                       **/
     680             : /**                                                                **/
     681             : /********************************************************************/
     682             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     683             :  * but everything else is */
     684             : static GEN
     685      166661 : QXQ_to_mod(GEN x, GEN T)
     686             : {
     687             :   long d;
     688      166661 :   switch(typ(x))
     689             :   {
     690       62394 :     case t_INT:  return icopy(x);
     691        2079 :     case t_FRAC: return gcopy(x);
     692      102188 :     case t_POL:
     693      102188 :       d = degpol(x);
     694      102188 :       if (d < 0) return gen_0;
     695       99927 :       if (d == 0) return gcopy(gel(x,2));
     696       97816 :       return mkpolmod(RgX_copy(x), T);
     697           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     698             :              return NULL;/* LCOV_EXCL_LINE */
     699             :   }
     700             : }
     701             : /* pure shallow version */
     702             : GEN
     703      839309 : QXQ_to_mod_shallow(GEN x, GEN T)
     704             : {
     705             :   long d;
     706      839309 :   switch(typ(x))
     707             :   {
     708      544308 :     case t_INT:
     709      544308 :     case t_FRAC: return x;
     710      295001 :     case t_POL:
     711      295001 :       d = degpol(x);
     712      295001 :       if (d < 0) return gen_0;
     713      248750 :       if (d == 0) return gel(x,2);
     714      232150 :       return mkpolmod(x, T);
     715           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     716             :              return NULL;/* LCOV_EXCL_LINE */
     717             :   }
     718             : }
     719             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
     720             :  * Not memory clean because T not copied, but everything else is */
     721             : static GEN
     722       37183 : QXQX_to_mod(GEN z, GEN T)
     723             : {
     724       37183 :   long i,l = lg(z);
     725       37183 :   GEN x = cgetg(l,t_POL);
     726      185594 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     727       37182 :   x[1] = z[1]; return normalizepol_lg(x,l);
     728             : }
     729             : /* pure shallow version */
     730             : GEN
     731      203487 : QXQX_to_mod_shallow(GEN z, GEN T)
     732             : {
     733      203487 :   long i,l = lg(z);
     734      203487 :   GEN x = cgetg(l,t_POL);
     735      959608 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
     736      203487 :   x[1] = z[1]; return normalizepol_lg(x,l);
     737             : }
     738             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     739             : GEN
     740       12516 : QXQXV_to_mod(GEN V, GEN T)
     741             : {
     742       12516 :   long i, l = lg(V);
     743       12516 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     744       49699 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     745       12516 :   return z;
     746             : }
     747             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     748             : GEN
     749       18358 : QXQV_to_mod(GEN V, GEN T)
     750             : {
     751       18358 :   long i, l = lg(V);
     752       18358 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     753       36607 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
     754       18358 :   return z;
     755             : }
     756             : 
     757             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     758             : GEN
     759       14644 : QXQC_to_mod_shallow(GEN V, GEN T)
     760             : {
     761       14644 :   long i, l = lg(V);
     762       14644 :   GEN z = cgetg(l, t_COL);
     763       97832 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
     764       14644 :   return z;
     765             : }
     766             : 
     767             : GEN
     768        6622 : QXQM_to_mod_shallow(GEN V, GEN T)
     769             : {
     770        6622 :   long i, l = lg(V);
     771        6622 :   GEN z = cgetg(l, t_MAT);
     772       21266 :   for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
     773        6622 :   return z;
     774             : }
     775             : 
     776             : GEN
     777     7885676 : RgX_renormalize_lg(GEN x, long lx)
     778             : {
     779             :   long i;
     780    10927277 :   for (i = lx-1; i>1; i--)
     781    10474337 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     782     7885676 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     783     7885676 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     784             : }
     785             : 
     786             : GEN
     787     1492318 : RgV_to_RgX(GEN x, long v)
     788             : {
     789     1492318 :   long i, k = lg(x);
     790             :   GEN p;
     791             : 
     792     4152610 :   while (--k && gequal0(gel(x,k)));
     793     1492318 :   if (!k) return pol_0(v);
     794     1485068 :   i = k+2; p = cgetg(i,t_POL);
     795     1485061 :   p[1] = evalsigne(1) | evalvarn(v);
     796    10228897 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     797     1485061 :   return p;
     798             : }
     799             : GEN
     800      188283 : RgV_to_RgX_reverse(GEN x, long v)
     801             : {
     802      188283 :   long j, k, l = lg(x);
     803             :   GEN p;
     804             : 
     805      189851 :   for (k = 1; k < l; k++)
     806      189851 :     if (!gequal0(gel(x,k))) break;
     807      188283 :   if (k == l) return pol_0(v);
     808      188283 :   k -= 1;
     809      188283 :   l -= k;
     810      188283 :   x += k;
     811      188283 :   p = cgetg(l+1,t_POL);
     812      188283 :   p[1] = evalsigne(1) | evalvarn(v);
     813      985038 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     814      188283 :   return p;
     815             : }
     816             : 
     817             : /* return the (N-dimensional) vector of coeffs of p */
     818             : GEN
     819    13407655 : RgX_to_RgC(GEN x, long N)
     820             : {
     821             :   long i, l;
     822             :   GEN z;
     823    13407655 :   l = lg(x)-1; x++;
     824    13407655 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     825    13407655 :   z = cgetg(N+1,t_COL);
     826    80851169 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     827    24469331 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     828    13407746 :   return z;
     829             : }
     830             : GEN
     831     1355017 : Rg_to_RgC(GEN x, long N)
     832             : {
     833     1355017 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     834             : }
     835             : 
     836             : /* vector of polynomials (in v) whose coefs are given by the columns of x */
     837             : GEN
     838      301959 : RgM_to_RgXV(GEN x, long v)
     839     1284546 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
     840             : GEN
     841        7202 : RgM_to_RgXV_reverse(GEN x, long v)
     842       28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
     843             : 
     844             : /* matrix whose entries are given by the coeffs of the polynomials in
     845             :  * vector v (considered as degree n-1 polynomials) */
     846             : GEN
     847      335569 : RgV_to_RgM(GEN x, long n)
     848     1686055 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
     849             : 
     850             : GEN
     851       78528 : RgXV_to_RgM(GEN x, long n)
     852      399387 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
     853             : 
     854             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     855             : GEN
     856       23469 : RgM_to_RgXX(GEN x, long v,long w)
     857             : {
     858       23469 :   long j, lx = lg(x);
     859       23469 :   GEN y = cgetg(lx+1, t_POL);
     860       23469 :   y[1] = evalsigne(1) | evalvarn(v);
     861       23469 :   y++;
     862      130426 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     863       23469 :   return normalizepol_lg(--y, lx+1);
     864             : }
     865             : 
     866             : /* matrix whose entries are given by the coeffs of the polynomial v in
     867             :  * two variables (considered as degree n-1 polynomials) */
     868             : GEN
     869         322 : RgXX_to_RgM(GEN v, long n)
     870             : {
     871         322 :   long j, N = lg(v)-1;
     872         322 :   GEN y = cgetg(N, t_MAT);
     873        1043 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     874         322 :   return y;
     875             : }
     876             : 
     877             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     878             : GEN
     879       32769 : RgXY_swapspec(GEN x, long n, long w, long nx)
     880             : {
     881       32769 :   long j, ly = n+3;
     882       32769 :   GEN y = cgetg(ly, t_POL);
     883       32769 :   y[1] = evalsigne(1);
     884      401026 :   for (j=2; j<ly; j++)
     885             :   {
     886             :     long k;
     887      368257 :     GEN a = cgetg(nx+2,t_POL);
     888      368257 :     a[1] = evalsigne(1) | evalvarn(w);
     889     2028285 :     for (k=0; k<nx; k++)
     890             :     {
     891     1660028 :       GEN xk = gel(x,k);
     892     1660028 :       if (typ(xk)==t_POL && varn(xk)==w)
     893     1565653 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     894             :       else
     895       94375 :         gel(a,k+2) = j==2 ? xk: gen_0;
     896             :     }
     897      368257 :     gel(y,j) = normalizepol_lg(a, nx+2);
     898             :   }
     899       32769 :   return normalizepol_lg(y,ly);
     900             : }
     901             : 
     902             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     903             : GEN
     904        2051 : RgXY_swap(GEN x, long n, long w)
     905             : {
     906        2051 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     907        2051 :   setvarn(z, varn(x)); return z;
     908             : }
     909             : 
     910             : long
     911        1387 : RgXY_degreex(GEN b)
     912             : {
     913        1387 :   long deg = 0, i;
     914        1387 :   if (!signe(b)) return -1;
     915       15592 :   for (i = 2; i < lg(b); ++i)
     916             :   {
     917       14205 :     GEN bi = gel(b, i);
     918       14205 :     if (typ(bi) == t_POL)
     919       13742 :       deg = maxss(deg, degpol(bi));
     920             :   }
     921        1387 :   return deg;
     922             : }
     923             : 
     924             : GEN
     925       38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
     926             : 
     927             : /* return (x % X^n). Shallow */
     928             : GEN
     929     7943806 : RgXn_red_shallow(GEN a, long n)
     930             : {
     931     7943806 :   long i, L = n+2, l = lg(a);
     932             :   GEN  b;
     933     7943806 :   if (L >= l) return a; /* deg(x) < n */
     934     5751553 :   b = cgetg(L, t_POL); b[1] = a[1];
     935    37323243 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     936     5751553 :   return normalizepol_lg(b,L);
     937             : }
     938             : 
     939             : GEN
     940         483 : RgXnV_red_shallow(GEN x, long n)
     941        2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
     942             : 
     943             : /* return (x * X^n). Shallow */
     944             : GEN
     945   175587709 : RgX_shift_shallow(GEN a, long n)
     946             : {
     947   175587709 :   long i, l = lg(a);
     948             :   GEN  b;
     949   175587709 :   if (l == 2 || !n) return a;
     950   109163726 :   l += n;
     951   109163726 :   if (n < 0)
     952             :   {
     953    54203261 :     if (l <= 2) return pol_0(varn(a));
     954    52729164 :     b = cgetg(l, t_POL); b[1] = a[1];
     955    52729902 :     a -= n;
     956   163992919 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     957             :   } else {
     958    54960465 :     b = cgetg(l, t_POL); b[1] = a[1];
     959    54964444 :     a -= n; n += 2;
     960   119463251 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     961   211986137 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     962             :   }
     963   107694346 :   return b;
     964             : }
     965             : /* return (x * X^n). */
     966             : GEN
     967     1578029 : RgX_shift(GEN a, long n)
     968             : {
     969     1578029 :   long i, l = lg(a);
     970             :   GEN  b;
     971     1578029 :   if (l == 2 || !n) return RgX_copy(a);
     972     1577077 :   l += n;
     973     1577077 :   if (n < 0)
     974             :   {
     975        1442 :     if (l <= 2) return pol_0(varn(a));
     976        1372 :     b = cgetg(l, t_POL); b[1] = a[1];
     977        1372 :     a -= n;
     978        9534 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     979             :   } else {
     980     1575635 :     b = cgetg(l, t_POL); b[1] = a[1];
     981     1575635 :     a -= n; n += 2;
     982     3979076 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     983     4346090 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     984             :   }
     985     1577007 :   return b;
     986             : }
     987             : 
     988             : GEN
     989      317037 : RgX_rotate_shallow(GEN P, long k, long p)
     990             : {
     991      317037 :   long i, l = lgpol(P);
     992             :   GEN r;
     993      317037 :   if (signe(P)==0)
     994        1365 :     return pol_0(varn(P));
     995      315672 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     996     2100644 :   for(i=0; i<p; i++)
     997             :   {
     998     1784972 :     long s = 2+(i+k)%p;
     999     1784972 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
    1000             :   }
    1001      315672 :   return RgX_renormalize(r);
    1002             : }
    1003             : 
    1004             : GEN
    1005     2997279 : RgX_mulXn(GEN x, long d)
    1006             : {
    1007             :   pari_sp av;
    1008             :   GEN z;
    1009             :   long v;
    1010     2997279 :   if (d >= 0) return RgX_shift(x, d);
    1011     1472517 :   d = -d;
    1012     1472517 :   v = RgX_val(x);
    1013     1472517 :   if (v >= d) return RgX_shift(x, -d);
    1014     1472503 :   av = avma;
    1015     1472503 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
    1016     1472503 :   return gc_upto(av, z);
    1017             : }
    1018             : 
    1019             : long
    1020         588 : RgXV_maxdegree(GEN x)
    1021             : {
    1022         588 :   long d = -1, i, l = lg(x);
    1023        4494 :   for (i = 1; i < l; i++)
    1024        3906 :     d = maxss(d, degpol(gel(x,i)));
    1025         588 :   return d;
    1026             : }
    1027             : 
    1028             : long
    1029     3667765 : RgX_val(GEN x)
    1030             : {
    1031     3667765 :   long i, lx = lg(x);
    1032     3667765 :   if (lx == 2) return LONG_MAX;
    1033     4619541 :   for (i = 2; i < lx; i++)
    1034     4619485 :     if (!isexactzero(gel(x,i))) break;
    1035     3656838 :   if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
    1036     3656782 :   return i - 2;
    1037             : }
    1038             : long
    1039    80776371 : RgX_valrem(GEN x, GEN *Z)
    1040             : {
    1041    80776371 :   long v, i, lx = lg(x);
    1042    80776371 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
    1043   125187910 :   for (i = 2; i < lx; i++)
    1044   125189411 :     if (!isexactzero(gel(x,i))) break;
    1045             :   /* possible with nonrational zeros */
    1046    80776362 :   if (i == lx)
    1047             :   {
    1048          14 :     *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
    1049          14 :     return LONG_MAX;
    1050             :   }
    1051    80776348 :   v = i - 2;
    1052    80776348 :   *Z = RgX_shift_shallow(x, -v);
    1053    80788765 :   return v;
    1054             : }
    1055             : long
    1056      851993 : RgX_valrem_inexact(GEN x, GEN *Z)
    1057             : {
    1058             :   long v;
    1059      851993 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
    1060      868980 :   for (v = 0;; v++)
    1061      868980 :     if (!gequal0(gel(x,2+v))) break;
    1062      851986 :   if (Z) *Z = RgX_shift_shallow(x, -v);
    1063      851986 :   return v;
    1064             : }
    1065             : 
    1066             : GEN
    1067       67683 : RgXQC_red(GEN x, GEN T)
    1068      412979 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
    1069             : 
    1070             : GEN
    1071        1211 : RgXQV_red(GEN x, GEN T)
    1072       27832 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
    1073             : 
    1074             : GEN
    1075       12999 : RgXQM_red(GEN x, GEN T)
    1076       80682 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
    1077             : 
    1078             : GEN
    1079         322 : RgXQM_mul(GEN P, GEN Q, GEN T)
    1080             : {
    1081         322 :   return RgXQM_red(RgM_mul(P, Q), T);
    1082             : }
    1083             : 
    1084             : GEN
    1085      508195 : RgXQX_red(GEN P, GEN T)
    1086             : {
    1087      508195 :   long i, l = lg(P);
    1088      508195 :   GEN Q = cgetg(l, t_POL);
    1089      508195 :   Q[1] = P[1];
    1090     2626816 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
    1091      508194 :   return normalizepol_lg(Q, l);
    1092             : }
    1093             : 
    1094             : GEN
    1095      825254 : RgX_deriv(GEN x)
    1096             : {
    1097      825254 :   long i,lx = lg(x)-1;
    1098             :   GEN y;
    1099             : 
    1100      825254 :   if (lx<3) return pol_0(varn(x));
    1101      822153 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
    1102     3639405 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
    1103      822138 :   y[1] = x[1]; return normalizepol_lg(y,i);
    1104             : }
    1105             : 
    1106             : GEN
    1107     2583784 : RgX_recipspec_shallow(GEN x, long l, long n)
    1108             : {
    1109             :   long i;
    1110     2583784 :   GEN z = cgetg(n+2,t_POL);
    1111     2583791 :   z[1] = 0; z += 2;
    1112   144314959 :   for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
    1113     2833526 :   for(   ; i<n; i++) gel(z, n-i-1) = gen_0;
    1114     2583791 :   return normalizepol_lg(z-2,n+2);
    1115             : }
    1116             : 
    1117             : GEN
    1118      642858 : RgXn_recip_shallow(GEN P, long n)
    1119             : {
    1120      642858 :   GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
    1121      642866 :   setvarn(Q, varn(P));
    1122      642866 :   return Q;
    1123             : }
    1124             : 
    1125             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1126             : GEN
    1127       31367 : RgX_recip(GEN x)
    1128             : {
    1129             :   long lx, i, j;
    1130       31367 :   GEN y = cgetg_copy(x, &lx);
    1131      274442 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
    1132       31367 :   return normalizepol_lg(y,lx);
    1133             : }
    1134             : /* shallow version */
    1135             : GEN
    1136       59388 : RgX_recip_shallow(GEN x)
    1137             : {
    1138             :   long lx, i, j;
    1139       59388 :   GEN y = cgetg_copy(x, &lx);
    1140      356118 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1141       59388 :   return normalizepol_lg(y,lx);
    1142             : }
    1143             : 
    1144             : GEN
    1145     5112553 : RgX_recip_i(GEN x)
    1146             : {
    1147             :   long lx, i, j;
    1148     5112553 :   GEN y = cgetg_copy(x, &lx);
    1149    25943121 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1150     5112547 :   return y;
    1151             : }
    1152             : /*******************************************************************/
    1153             : /*                                                                 */
    1154             : /*                      ADDITION / SUBTRACTION                     */
    1155             : /*                                                                 */
    1156             : /*******************************************************************/
    1157             : /* same variable */
    1158             : GEN
    1159   144726598 : RgX_add(GEN x, GEN y)
    1160             : {
    1161   144726598 :   long i, lx = lg(x), ly = lg(y);
    1162             :   GEN z;
    1163   144726598 :   if (ly <= lx) {
    1164   130476014 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1165   508051903 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1166   190221220 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1167   130451959 :     z = normalizepol_lg(z, lx);
    1168             :   } else {
    1169    14250584 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1170    54482507 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1171    41633837 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1172    14251809 :     z = normalizepol_lg(z, ly);
    1173             :   }
    1174   144718126 :   return z;
    1175             : }
    1176             : GEN
    1177    72773695 : RgX_sub(GEN x, GEN y)
    1178             : {
    1179    72773695 :   long i, lx = lg(x), ly = lg(y);
    1180             :   GEN z;
    1181    72773695 :   if (ly <= lx) {
    1182    38939243 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1183   188189509 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1184    66413264 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1185    38915789 :     z = normalizepol_lg(z, lx);
    1186             :   } else {
    1187    33834452 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1188   125718707 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1189    74933384 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1190    33800281 :     z = normalizepol_lg(z, ly);
    1191             :   }
    1192    72769681 :   return z;
    1193             : }
    1194             : GEN
    1195     8820865 : RgX_neg(GEN x)
    1196    53145032 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
    1197             : 
    1198             : GEN
    1199    44377518 : RgX_Rg_add(GEN y, GEN x)
    1200             : {
    1201             :   GEN z;
    1202    44377518 :   long lz = lg(y), i;
    1203    44377518 :   if (lz == 2) return scalarpol(x,varn(y));
    1204    22633651 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1205    22633618 :   gel(z,2) = gadd(gel(y,2),x);
    1206    78670959 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1207             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1208             :    * an inexact 0 */
    1209    22633562 :   return normalizepol_lg(z,lz);
    1210             : }
    1211             : GEN
    1212       65073 : RgX_Rg_add_shallow(GEN y, GEN x)
    1213             : {
    1214             :   GEN z;
    1215       65073 :   long lz = lg(y), i;
    1216       65073 :   if (lz == 2) return scalarpol(x,varn(y));
    1217       65073 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1218       65073 :   gel(z,2) = gadd(gel(y,2),x);
    1219      130290 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1220       65073 :   return normalizepol_lg(z,lz);
    1221             : }
    1222             : GEN
    1223     1222883 : RgX_Rg_sub(GEN y, GEN x)
    1224             : {
    1225             :   GEN z;
    1226     1222883 :   long lz = lg(y), i;
    1227     1222883 :   if (lz == 2)
    1228             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1229      922867 :     long v = varn(y);
    1230      922867 :     if (isrationalzero(x)) return pol_0(v);
    1231        2235 :     z = cgetg(3,t_POL);
    1232        2235 :     z[1] = gequal0(x)? evalvarn(v)
    1233        2235 :                    : evalvarn(v) | evalsigne(1);
    1234        2235 :     gel(z,2) = gneg(x); return z;
    1235             :   }
    1236      300016 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1237      300016 :   gel(z,2) = gsub(gel(y,2),x);
    1238      686822 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1239      300016 :   return normalizepol_lg(z,lz);
    1240             : }
    1241             : GEN
    1242     4951641 : Rg_RgX_sub(GEN x, GEN y)
    1243             : {
    1244             :   GEN z;
    1245     4951641 :   long lz = lg(y), i;
    1246     4951641 :   if (lz == 2) return scalarpol(x,varn(y));
    1247     4881247 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1248     4881230 :   gel(z,2) = gsub(x, gel(y,2));
    1249     7333194 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1250     4881202 :   return normalizepol_lg(z,lz);
    1251             : }
    1252             : /*******************************************************************/
    1253             : /*                                                                 */
    1254             : /*                  KARATSUBA MULTIPLICATION                       */
    1255             : /*                                                                 */
    1256             : /*******************************************************************/
    1257             : #if 0
    1258             : /* to debug Karatsuba-like routines */
    1259             : GEN
    1260             : zx_debug_spec(GEN x, long nx)
    1261             : {
    1262             :   GEN z = cgetg(nx+2,t_POL);
    1263             :   long i;
    1264             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1265             :   z[1] = evalsigne(1); return z;
    1266             : }
    1267             : 
    1268             : GEN
    1269             : RgX_debug_spec(GEN x, long nx)
    1270             : {
    1271             :   GEN z = cgetg(nx+2,t_POL);
    1272             :   long i;
    1273             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1274             :   z[1] = evalsigne(1); return z;
    1275             : }
    1276             : #endif
    1277             : 
    1278             : /* generic multiplication */
    1279             : GEN
    1280     8907728 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1281             : {
    1282             :   GEN z, t;
    1283             :   long i;
    1284     8907728 :   if (nx == ny) {
    1285     1536128 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1286     4953190 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1287     1536118 :     return normalizepol_lg(z, nx+2);
    1288             :   }
    1289     7371600 :   if (ny < nx) {
    1290     7186372 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1291    26452576 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1292    17172257 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1293     7185836 :     return normalizepol_lg(z, nx+2);
    1294             :   } else {
    1295      185228 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1296     3672041 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1297      447238 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1298      185239 :     return normalizepol_lg(z, ny+2);
    1299             :   }
    1300             : }
    1301             : GEN
    1302      222319 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1303             : {
    1304             :   GEN z, t;
    1305             :   long i;
    1306      222319 :   if (nx == ny) {
    1307       12824 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1308     2185778 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1309       12824 :     return normalizepol_lg(z, nx+2);
    1310             :   }
    1311      209495 :   if (ny < nx) {
    1312      207717 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1313     3724789 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1314     2371725 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1315      207717 :     return normalizepol_lg(z, nx+2);
    1316             :   } else {
    1317        1778 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1318      330904 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1319       12222 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1320        1778 :     return normalizepol_lg(z, ny+2);
    1321             :   }
    1322             : }
    1323             : 
    1324             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1325             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1326             :  * t_VECSMALL, to hold this, which can be left on stack: gc_GEN_unsafe
    1327             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1328             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1329             : static GEN
    1330    81694057 : RgXspec_kill0(GEN x, long lx)
    1331             : {
    1332    81694057 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gc_GEN_unsafe-wise */
    1333             :   long i;
    1334   246568526 :   for (i=0; i <lx; i++)
    1335             :   {
    1336   164874465 :     GEN c = gel(x,i);
    1337   164874465 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1338             :   }
    1339    81694061 :   return z;
    1340             : }
    1341             : 
    1342             : INLINE GEN
    1343   111053993 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1344             : {
    1345   111053993 :   pari_sp av = avma;
    1346   111053993 :   GEN s = NULL;
    1347             :   long i;
    1348             : 
    1349   391653101 :   for (i=a; i<b; i++)
    1350   280608433 :     if (gel(y,i) && gel(x,-i))
    1351             :     {
    1352   202283801 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1353   202279292 :       s = s? gadd(s, t): t;
    1354             :     }
    1355   111044668 :   return s? gc_upto(av, s): gen_0;
    1356             : }
    1357             : 
    1358             : /* assume nx >= ny > 0, return x * y * t^v */
    1359             : static GEN
    1360    34200163 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1361             : {
    1362             :   long i, lz, nz;
    1363             :   GEN z;
    1364             : 
    1365    34200163 :   x = RgXspec_kill0(x,nx);
    1366    34200118 :   y = RgXspec_kill0(y,ny);
    1367    34200132 :   lz = nx + ny + 1; nz = lz-2;
    1368    34200132 :   lz += v;
    1369    34200132 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1370    72777285 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1371    83192635 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1372    56628541 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1373    48992536 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1374    34199566 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1375             : }
    1376             : 
    1377             : /* return (x * X^d) + y. Assume d > 0 */
    1378             : GEN
    1379    10066115 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1380             : {
    1381             :   GEN x, y, xd, yd, zd;
    1382             :   long a, lz, nx, ny;
    1383             : 
    1384    10066115 :   if (!signe(x0)) return y0;
    1385     9431284 :   ny = lgpol(y0);
    1386     9431283 :   nx = lgpol(x0);
    1387     9431366 :   zd = (GEN)avma;
    1388     9431366 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1389     9431366 :   if (a <= 0)
    1390             :   {
    1391     1479446 :     lz = nx+d+2;
    1392     1479446 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1393     3397428 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1394     1479446 :     x = zd + a;
    1395     1497143 :     while (zd > x) gel(--zd,0) = gen_0;
    1396             :   }
    1397             :   else
    1398             :   {
    1399     7951920 :     xd = new_chunk(d); yd = y+d;
    1400     7951925 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1401     7951841 :     lz = (a>nx)? ny+2: lg(x)+d;
    1402    39625065 :     x += 2; while (xd > x) *--zd = *--xd;
    1403             :   }
    1404    22974012 :   while (yd > y) *--zd = *--yd;
    1405     9431287 :   *--zd = x0[1];
    1406     9431287 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1407             : }
    1408             : GEN
    1409      514789 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1410             : {
    1411             :   GEN x, y, xd, yd, zd;
    1412             :   long a, lz, nx, ny;
    1413             : 
    1414      514789 :   if (!signe(x0)) return RgX_copy(y0);
    1415      514005 :   nx = lgpol(x0);
    1416      514005 :   ny = lgpol(y0);
    1417      514005 :   zd = (GEN)avma;
    1418      514005 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1419      514005 :   if (a <= 0)
    1420             :   {
    1421      291686 :     lz = nx+d+2;
    1422      291686 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1423     4293251 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1424      291686 :     x = zd + a;
    1425      757796 :     while (zd > x) gel(--zd,0) = gen_0;
    1426             :   }
    1427             :   else
    1428             :   {
    1429      222319 :     xd = new_chunk(d); yd = y+d;
    1430      222319 :     x = RgX_addspec(x,yd, nx,a);
    1431      222319 :     lz = (a>nx)? ny+2: lg(x)+d;
    1432     8415923 :     x += 2; while (xd > x) *--zd = *--xd;
    1433             :   }
    1434     2615106 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1435      514005 :   *--zd = x0[1];
    1436      514005 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1437             : }
    1438             : 
    1439             : /* return x * y mod t^n */
    1440             : static GEN
    1441     6627959 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1442             : {
    1443     6627959 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1444             :   GEN z;
    1445     6627959 :   if (lx < 0) return pol_0(varn(x));
    1446     6627959 :   if (ly < 0) return pol_0(varn(x));
    1447     6627959 :   z = cgetg(lz, t_POL) + 2;
    1448     6627959 :   x+=2; if (lx > n) lx = n;
    1449     6627959 :   y+=2; if (ly > n) ly = n;
    1450     6627959 :   z[-1] = x[-1];
    1451     6627959 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1452     6627959 :   x = RgXspec_kill0(x, lx);
    1453     6627959 :   y = RgXspec_kill0(y, ly);
    1454             :   /* x:y:z [i] = term of degree i */
    1455    26223885 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1456    11826184 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1457     6674031 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1458     6627959 :   return normalizepol_lg(z - 2, lz);
    1459             : }
    1460             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1461             : static GEN
    1462     9414318 : RgXn_mul2(GEN f, GEN g, long n)
    1463             : {
    1464     9414318 :   pari_sp av = avma;
    1465             :   GEN fe,fo, ge,go, l,h,m;
    1466             :   long n0, n1;
    1467     9414318 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1468     6663449 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1469       35490 :   n0 = n>>1; n1 = n-n0;
    1470       35490 :   RgX_even_odd(f, &fe, &fo);
    1471       35490 :   RgX_even_odd(g, &ge, &go);
    1472       35490 :   l = RgXn_mul2(fe,ge,n1);
    1473       35490 :   h = RgXn_mul2(fo,go,n0);
    1474       35490 :   m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1475             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1476             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1477       35490 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1478             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1479       35490 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1480       35490 :   m = RgX_inflate(m,2);
    1481             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1482       35490 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1483       35490 :   h = RgX_inflate(h,2);
    1484       35490 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1485       35490 :   return gc_upto(av, h);
    1486             : }
    1487             : /* (f*g) \/ x^n */
    1488             : static GEN
    1489     1589447 : RgX_mulhigh_i2(GEN f, GEN g, long n)
    1490             : {
    1491     1589447 :   long d = degpol(f)+degpol(g) + 1 - n;
    1492             :   GEN h;
    1493     1589447 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1494       29654 :   h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
    1495             :                                   RgX_recip_i(g), d));
    1496       29654 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1497             : }
    1498             : 
    1499             : /* (f*g) \/ x^n */
    1500             : static GEN
    1501           0 : RgX_sqrhigh_i2(GEN f, long n)
    1502             : {
    1503           0 :   long d = 2*degpol(f)+ 1 - n;
    1504             :   GEN h;
    1505           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1506           0 :   h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
    1507           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1508             : }
    1509             : 
    1510             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1511             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1512             :  * Only c, c0, c1, c2 are genuine GEN.
    1513             :  */
    1514             : GEN
    1515    40266375 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1516             : {
    1517             :   GEN a0, c, c0;
    1518    40266375 :   long n0, n0a, i, v = 0;
    1519             :   pari_sp av;
    1520             : 
    1521    65264683 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1522    61957057 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1523    40266345 :   if (na < nb) swapspec(a,b, na,nb);
    1524    40266345 :   if (!nb) return pol_0(0);
    1525             : 
    1526    34679404 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1527      479255 :   RgX_shift_inplace_init(v);
    1528      479260 :   i = (na>>1); n0 = na-i; na = i;
    1529      479260 :   av = avma; a0 = a+n0; n0a = n0;
    1530     1334119 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1531             : 
    1532      479260 :   if (nb > n0)
    1533             :   {
    1534             :     GEN b0,c1,c2;
    1535             :     long n0b;
    1536             : 
    1537      477901 :     nb -= n0; b0 = b+n0; n0b = n0;
    1538     1445031 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1539      477901 :     c = RgX_mulspec(a,b,n0a,n0b);
    1540      477901 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1541             : 
    1542      477901 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1543      477901 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1544             : 
    1545      477901 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1546      477901 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1547      477901 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1548             :   }
    1549             :   else
    1550             :   {
    1551        1359 :     c = RgX_mulspec(a,b,n0a,nb);
    1552        1359 :     c0 = RgX_mulspec(a0,b,na,nb);
    1553             :   }
    1554      479260 :   c0 = RgX_addmulXn(c0,c,n0);
    1555      479260 :   return RgX_shift_inplace(gc_upto(av,c0), v);
    1556             : }
    1557             : 
    1558             : INLINE GEN
    1559      100082 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1560             : {
    1561      100082 :   pari_sp av = avma;
    1562      100082 :   GEN s = NULL;
    1563      100082 :   long j, l = (i+1)>>1;
    1564      203924 :   for (j=a; j<l; j++)
    1565             :   {
    1566      103842 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1567      103842 :     if (xj && xx)
    1568             :     {
    1569       94707 :       GEN t = gmul(xj, xx);
    1570       94707 :       s = s? gadd(s, t): t;
    1571             :     }
    1572             :   }
    1573      100082 :   if (s) s = gshift(s,1);
    1574      100082 :   if ((i&1) == 0)
    1575             :   {
    1576       69041 :     GEN t = gel(x, i>>1);
    1577       69041 :     if (t) {
    1578       65982 :       t = gsqr(t);
    1579       65982 :       s = s? gadd(s, t): t;
    1580             :     }
    1581             :   }
    1582      100082 :   return s? gc_upto(av,s): gen_0;
    1583             : }
    1584             : static GEN
    1585       38000 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1586             : {
    1587             :   long i, lz, nz;
    1588             :   GEN z;
    1589             : 
    1590       38000 :   if (!nx) return pol_0(0);
    1591       38000 :   x = RgXspec_kill0(x,nx);
    1592       38000 :   lz = (nx << 1) + 1, nz = lz-2;
    1593       38000 :   lz += v;
    1594       38000 :   z = cgetg(lz,t_POL) + 2;
    1595       94602 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1596      107041 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1597       69041 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1598       38000 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1599             : }
    1600             : /* return x^2 mod t^n */
    1601             : static GEN
    1602           0 : RgXn_sqr_basecase(GEN x, long n)
    1603             : {
    1604           0 :   long i, lz = n+2, lx = lgpol(x);
    1605             :   GEN z;
    1606           0 :   if (lx < 0) return pol_0(varn(x));
    1607           0 :   z = cgetg(lz, t_POL);
    1608           0 :   z[1] = x[1];
    1609           0 :   x+=2; if (lx > n) lx = n;
    1610           0 :   x = RgXspec_kill0(x,lx);
    1611           0 :   z+=2;/* x:z [i] = term of degree i */
    1612           0 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1613           0 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1614           0 :   z -= 2; return normalizepol_lg(z, lz);
    1615             : }
    1616             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1617             : static GEN
    1618         315 : RgXn_sqr2(GEN f, long n)
    1619             : {
    1620         315 :   pari_sp av = avma;
    1621             :   GEN fe,fo, l,h,m;
    1622             :   long n0, n1;
    1623         315 :   if (2*degpol(f) < n) return RgX_sqr_i(f);
    1624           0 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1625           0 :   n0 = n>>1; n1 = n-n0;
    1626           0 :   RgX_even_odd(f, &fe, &fo);
    1627           0 :   l = RgXn_sqr(fe,n1);
    1628           0 :   h = RgXn_sqr(fo,n0);
    1629           0 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1630             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1631             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1632           0 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1633             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1634           0 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1635           0 :   m = RgX_inflate(m,2);
    1636             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1637           0 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1638           0 :   h = RgX_inflate(h,2);
    1639           0 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1640           0 :   return gc_upto(av, h);
    1641             : }
    1642             : GEN
    1643       38039 : RgX_sqrspec(GEN a, long na)
    1644             : {
    1645             :   GEN a0, c, c0, c1;
    1646       38039 :   long n0, n0a, i, v = 0;
    1647             :   pari_sp av;
    1648             : 
    1649       66340 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1650       38039 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1651          39 :   RgX_shift_inplace_init(v);
    1652          39 :   i = (na>>1); n0 = na-i; na = i;
    1653          39 :   av = avma; a0 = a+n0; n0a = n0;
    1654          39 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1655             : 
    1656          39 :   c = RgX_sqrspec(a,n0a);
    1657          39 :   c0 = RgX_sqrspec(a0,na);
    1658          39 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1659          39 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1660          39 :   c0 = RgX_addmulXn(c0,c,n0);
    1661          39 :   return RgX_shift_inplace(gc_upto(av,c0), v);
    1662             : }
    1663             : 
    1664             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1665             : GEN
    1666     1719019 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1667             : {
    1668     1719019 :   GEN z = RgX_mul(A, B);
    1669     1719005 :   if (a < b)
    1670       10372 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1671     1708633 :   else if (a > b)
    1672     1104943 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1673             :   else
    1674      603690 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1675     1719012 :   return z;
    1676             : }
    1677             : 
    1678             : GEN
    1679    38829987 : RgX_mul_i(GEN x, GEN y)
    1680             : {
    1681    38829987 :   GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
    1682    38829569 :   setvarn(z, varn(x)); return z;
    1683             : }
    1684             : 
    1685             : GEN
    1686       37961 : RgX_sqr_i(GEN x)
    1687             : {
    1688       37961 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1689       37961 :   setvarn(z,varn(x)); return z;
    1690             : }
    1691             : 
    1692             : /*******************************************************************/
    1693             : /*                                                                 */
    1694             : /*                               DIVISION                          */
    1695             : /*                                                                 */
    1696             : /*******************************************************************/
    1697             : GEN
    1698     7176935 : RgX_Rg_divexact(GEN x, GEN y) {
    1699     7176935 :   long i, lx = lg(x);
    1700             :   GEN z;
    1701     7176935 :   if (lx == 2) return gcopy(x);
    1702     7125976 :   switch(typ(y))
    1703             :   {
    1704     6964739 :     case t_INT:
    1705     6964739 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1706     6589397 :       break;
    1707        5243 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1708             :   }
    1709     6745391 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1710    34434704 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1711     6743724 :   return z;
    1712             : }
    1713             : GEN
    1714    30320927 : RgX_Rg_div(GEN x, GEN y) {
    1715    30320927 :   long i, lx = lg(x);
    1716             :   GEN z;
    1717    30320927 :   if (lx == 2) return gcopy(x);
    1718    30048050 :   switch(typ(y))
    1719             :   {
    1720    20815503 :     case t_INT:
    1721    20815503 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1722     4192047 :       break;
    1723        6055 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1724             :   }
    1725    13418539 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1726    50019608 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1727    13418117 :   return normalizepol_lg(z, lx);
    1728             : }
    1729             : GEN
    1730       39130 : RgX_normalize(GEN x)
    1731             : {
    1732       39130 :   GEN z, d = NULL;
    1733       39130 :   long i, n = lg(x)-1;
    1734       39130 :   for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
    1735       39130 :   if (i == 1) return pol_0(varn(x));
    1736       39130 :   if (i == n && isint1(d)) return x;
    1737       17612 :   n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
    1738       31612 :   for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
    1739       17612 :   gel(z,n) = Rg_get_1(d); return z;
    1740             : }
    1741             : GEN
    1742       10458 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
    1743             : GEN
    1744      251003 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1745             : {
    1746      251003 :   long l = lg(a), i;
    1747             :   GEN a0, z0, z;
    1748             : 
    1749      251003 :   if (l <= 3)
    1750             :   {
    1751           0 :     if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
    1752           0 :     return pol_0(varn(a));
    1753             :   }
    1754      251003 :   z = cgetg(l-1, t_POL);
    1755      251007 :   z[1] = a[1];
    1756      251007 :   a0 = a + l-1;
    1757      251007 :   z0 = z + l-2; *z0 = *a0--;
    1758     2954660 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1759             :   {
    1760     2703663 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1761     2703653 :     gel(z0,0) = t;
    1762             :   }
    1763      250997 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1764      250997 :   return z;
    1765             : }
    1766             : /* Polynomial division x / y:
    1767             :  *   if pr = ONLY_REM return remainder, otherwise return quotient
    1768             :  *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
    1769             :  *   if pr != NULL set *pr to remainder, as the last object on stack */
    1770             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1771             : static GEN
    1772    23512200 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
    1773             : {
    1774             :   pari_sp avy, av, av1;
    1775             :   long dx,dy,dz,i,j,sx,lr;
    1776             :   GEN z,p1,p2,rem,y_lead,mod,p;
    1777             :   GEN (*f)(GEN,GEN);
    1778             : 
    1779    23512200 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1780             : 
    1781    23512200 :   dy = degpol(y);
    1782    23512188 :   y_lead = gel(y,dy+2);
    1783    23512188 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1784             :   {
    1785           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1786           0 :     for (dy--; dy>=0; dy--)
    1787             :     {
    1788           0 :       y_lead = gel(y,dy+2);
    1789           0 :       if (!gequal0(y_lead)) break;
    1790             :     }
    1791             :   }
    1792    23512182 :   if (!dy) /* y is constant */
    1793             :   {
    1794        6844 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1795        6844 :     z = RgX_Rg_div(x, y_lead);
    1796        6844 :     if (pr == ONLY_DIVIDES) return z;
    1797        2630 :     if (pr) *pr = pol_0(varn(x));
    1798        2630 :     return z;
    1799             :   }
    1800    23505338 :   dx = degpol(x);
    1801    23505307 :   if (dx < dy)
    1802             :   {
    1803     3847534 :     if (pr == ONLY_REM) return RgX_copy(x);
    1804      348542 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1805      348521 :     z = pol_0(varn(x));
    1806      348521 :     if (pr) *pr = RgX_copy(x);
    1807      348521 :     return z;
    1808             :   }
    1809             : 
    1810             :   /* x,y in R[X], y non constant */
    1811    19657773 :   av = avma;
    1812    19657773 :   p = NULL;
    1813    19657773 :   if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1814             :   {
    1815      130277 :     z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
    1816      130277 :     if (!z) return gc_NULL(av);
    1817      130277 :     z = FpX_to_mod(z, p);
    1818      130277 :     if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
    1819       71134 :       return gc_upto(av, z);
    1820       59143 :     *pr = FpX_to_mod(*pr, p);
    1821       59143 :     return gc_all(av, 2, &z, pr);
    1822             :   }
    1823    19527630 :   switch(typ(y_lead))
    1824             :   {
    1825           0 :     case t_REAL:
    1826           0 :       y_lead = ginv(y_lead);
    1827           0 :       f = gmul; mod = NULL;
    1828           0 :       break;
    1829        1878 :     case t_INTMOD:
    1830        1878 :     case t_POLMOD: y_lead = ginv(y_lead);
    1831        1878 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1832        1878 :       break;
    1833    19525752 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1834    19525740 :       f = gdiv; mod = NULL;
    1835             :   }
    1836             : 
    1837    19527618 :   if (y_lead == NULL)
    1838    17515766 :     p2 = gel(x,dx+2);
    1839             :   else {
    1840             :     for(;;) {
    1841     2011852 :       p2 = f(gel(x,dx+2),y_lead);
    1842     2011853 :       p2 = simplify_shallow(p2);
    1843     2011853 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1844             :     }
    1845     2011853 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1846             :     {
    1847           0 :       if (pr == ONLY_DIVIDES) {
    1848           0 :         set_avma(av);
    1849           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1850             :       }
    1851           0 :       if (pr == ONLY_REM)
    1852             :       {
    1853           0 :         if (dx < 0)
    1854           0 :           return gc_GEN(av, scalarpol(p2, varn(x)));
    1855             :         else
    1856             :         {
    1857             :           GEN t;
    1858           0 :           set_avma(av);
    1859           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1860           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1861           0 :           return t;
    1862             :         }
    1863             :       }
    1864           0 :       if (pr) /* cf ONLY_REM above */
    1865             :       {
    1866           0 :         if (dx < 0)
    1867             :         {
    1868           0 :           p2 = gclone(p2);
    1869           0 :           set_avma(av);
    1870           0 :           z = pol_0(varn(x));
    1871           0 :           x = scalarpol(p2, varn(x));
    1872           0 :           gunclone(p2);
    1873             :         }
    1874             :         else
    1875             :         {
    1876             :           GEN t;
    1877           0 :           set_avma(av);
    1878           0 :           z = pol_0(varn(x));
    1879           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1880           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1881           0 :           x = t;
    1882             :         }
    1883           0 :         *pr = x;
    1884             :       }
    1885             :       else
    1886             :       {
    1887           0 :         set_avma(av);
    1888           0 :         z = pol_0(varn(x));
    1889             :       }
    1890           0 :       return z;
    1891             :     }
    1892             :   }
    1893             :   /* dx >= dy */
    1894    19527619 :   avy = avma;
    1895    19527619 :   dz = dx-dy;
    1896    19527619 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1897    19527577 :   x += 2;
    1898    19527577 :   z += 2;
    1899    19527577 :   y += 2;
    1900    19527577 :   gel(z,dz) = gcopy(p2);
    1901             : 
    1902    54085011 :   for (i=dx-1; i>=dy; i--)
    1903             :   {
    1904    34557588 :     av1=avma; p1=gel(x,i);
    1905  1140249587 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1906    34552289 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1907             : 
    1908    34552289 :     if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
    1909             :     else
    1910    24563311 :       p1 = avma==av1? gcopy(p1): gc_upto(av1,p1);
    1911    34557024 :     gel(z,i-dy) = p1;
    1912             :   }
    1913    19527423 :   if (!pr) return gc_upto(av,z-2);
    1914             : 
    1915    12550230 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1916    12907253 :   for (sx=0; ; i--)
    1917             :   {
    1918    12907253 :     p1 = gel(x,i);
    1919             :     /* we always enter this loop at least once */
    1920    32218666 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1921    12906054 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1922    12906054 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
    1923     3576633 :     if (!isexactzero(p1)) break;
    1924     3426506 :     if (!i) break;
    1925      356989 :     set_avma(av1);
    1926             :   }
    1927    12549435 :   if (pr == ONLY_DIVIDES)
    1928             :   {
    1929        7931 :     if (sx) return gc_NULL(av);
    1930        7910 :     set_avma((pari_sp)rem); return gc_upto(av,z-2);
    1931             :   }
    1932    12541504 :   lr=i+3; rem -= lr;
    1933    12541504 :   if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
    1934    12499215 :   else p1 = gc_upto((pari_sp)rem,p1);
    1935    12542272 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    1936    12542272 :   rem[1] = z[-1];
    1937    12542272 :   rem += 2;
    1938    12542272 :   gel(rem,i) = p1;
    1939    17282589 :   for (i--; i>=0; i--)
    1940             :   {
    1941     4740349 :     av1=avma; p1 = gel(x,i);
    1942    13636227 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1943     4739842 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1944     4740074 :     gel(rem,i) = avma==av1? gcopy(p1):gc_upto(av1,p1);
    1945             :   }
    1946    12542240 :   rem -= 2;
    1947    12542240 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1948    12542463 :   if (pr == ONLY_REM) return gc_upto(av,rem);
    1949     6760690 :   z -= 2; *pr = rem; return gc_all_unsafe(av,avy,2,&z,pr);
    1950             : }
    1951             : 
    1952             : GEN
    1953    14231495 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1954             : {
    1955    14231495 :   if (pr == ONLY_REM) return RgX_rem(x, y);
    1956    14231495 :   return RgX_divrem_i(x, y, pr);
    1957             : }
    1958             : 
    1959             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1960             : GEN
    1961      155548 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1962             : {
    1963             :   long vx, dx, dy, dz, i, j, sx, lr;
    1964             :   pari_sp av0, av, tetpil;
    1965             :   GEN z,p1,rem,lead;
    1966             : 
    1967      155548 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1968      155548 :   vx = varn(x);
    1969      155548 :   dx = degpol(x);
    1970      155548 :   dy = degpol(y);
    1971      155548 :   if (dx < dy)
    1972             :   {
    1973       41112 :     if (pr)
    1974             :     {
    1975       41084 :       av0 = avma; x = RgXQX_red(x, T);
    1976       41084 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
    1977       41077 :       if (pr == ONLY_REM) return x;
    1978           0 :       *pr = x;
    1979             :     }
    1980          28 :     return pol_0(vx);
    1981             :   }
    1982      114436 :   lead = leading_coeff(y);
    1983      114436 :   if (!dy) /* y is constant */
    1984             :   {
    1985         602 :     if (pr && pr != ONLY_DIVIDES)
    1986             :     {
    1987           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1988           0 :       *pr = pol_0(vx);
    1989             :     }
    1990         602 :     if (gequal1(lead)) return RgX_copy(x);
    1991           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1992           0 :     return gc_GEN_unsafe(av0,tetpil,RgXQX_red(x,T));
    1993             :   }
    1994      113834 :   av0 = avma; dz = dx-dy;
    1995      113834 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1996      113834 :   set_avma(av0);
    1997      113834 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1998      113834 :   x += 2; y += 2; z += 2;
    1999             : 
    2000      113834 :   p1 = gel(x,dx); av = avma;
    2001      113834 :   gel(z,dz) = lead? gc_upto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    2002      579252 :   for (i=dx-1; i>=dy; i--)
    2003             :   {
    2004      465418 :     av=avma; p1=gel(x,i);
    2005     2412488 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    2006      465391 :     if (lead) p1 = gmul(grem(p1, T), lead);
    2007      465397 :     tetpil=avma; gel(z,i-dy) = gc_GEN_unsafe(av,tetpil, grem(p1, T));
    2008             :   }
    2009      113834 :   if (!pr) { guncloneNULL(lead); return z-2; }
    2010             : 
    2011      112560 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    2012      189213 :   for (sx=0; ; i--)
    2013             :   {
    2014      189213 :     p1 = gel(x,i);
    2015      654707 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    2016      189216 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    2017      103737 :     if (!i) break;
    2018       76654 :     set_avma(av);
    2019             :   }
    2020      112560 :   if (pr == ONLY_DIVIDES)
    2021             :   {
    2022       27066 :     guncloneNULL(lead);
    2023       27066 :     if (sx) return gc_NULL(av0);
    2024       24507 :     return gc_const((pari_sp)rem, z-2);
    2025             :   }
    2026       85494 :   lr=i+3; rem -= lr;
    2027       85494 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    2028       85494 :   rem[1] = z[-1];
    2029       85494 :   p1 = gc_GEN_unsafe((pari_sp)rem,tetpil,p1);
    2030       85494 :   rem += 2; gel(rem,i) = p1;
    2031      186658 :   for (i--; i>=0; i--)
    2032             :   {
    2033      101164 :     av=avma; p1 = gel(x,i);
    2034      337033 :     for (j=0; j<=i && j<=dz; j++)
    2035      235869 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    2036      101164 :     tetpil=avma; gel(rem,i) = gc_GEN_unsafe(av,tetpil, grem(p1, T));
    2037             :   }
    2038       85494 :   rem -= 2;
    2039       85494 :   guncloneNULL(lead);
    2040       85494 :   if (!sx) (void)normalizepol_lg(rem, lr);
    2041       85494 :   if (pr == ONLY_REM) return gc_upto(av0,rem);
    2042         161 :   *pr = rem; return z-2;
    2043             : }
    2044             : 
    2045             : /*******************************************************************/
    2046             : /*                                                                 */
    2047             : /*                        PSEUDO-DIVISION                          */
    2048             : /*                                                                 */
    2049             : /*******************************************************************/
    2050             : INLINE GEN
    2051     4042753 : rem(GEN c, GEN T)
    2052             : {
    2053     4042753 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    2054     4042753 :   return c;
    2055             : }
    2056             : 
    2057             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    2058             : int
    2059       17335 : ZXQX_dvd(GEN x, GEN y, GEN T)
    2060             : {
    2061             :   long dx, dy, i, T_ismonic;
    2062       17335 :   pari_sp av = avma, av2;
    2063             :   GEN y_lead;
    2064             : 
    2065       17335 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    2066       17335 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2067       17335 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    2068             :   /* if monic, no point in using pseudo-division */
    2069       17335 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    2070       14794 :   T_ismonic = gequal1(leading_coeff(T));
    2071       14794 :   dx = degpol(x);
    2072       14794 :   if (dx < dy) return !signe(x);
    2073       14794 :   (void)new_chunk(2);
    2074       14794 :   x = RgX_recip_i(x)+2;
    2075       14794 :   y = RgX_recip_i(y)+2;
    2076             :   /* pay attention to sparse divisors */
    2077       30089 :   for (i = 1; i <= dy; i++)
    2078       15295 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    2079       14794 :   av2 = avma;
    2080             :   for (;;)
    2081       73009 :   {
    2082       87803 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    2083       87803 :     x0 = gneg(x0);
    2084       87803 :     m = gcdii(cx, y0);
    2085       87803 :     if (!equali1(m))
    2086             :     {
    2087       85301 :       x0 = gdiv(x0, m);
    2088       85301 :       y0 = diviiexact(y0, m);
    2089       85301 :       if (equali1(y0)) y0 = NULL;
    2090             :     }
    2091      179379 :     for (i=1; i<=dy; i++)
    2092             :     {
    2093       91576 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2094       91576 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    2095       91576 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2096       91576 :       gel(x,i) = c;
    2097             :     }
    2098      688168 :     for (   ; i<=dx; i++)
    2099             :     {
    2100      600365 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2101      600365 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2102      600365 :       gel(x,i) = c;
    2103             :     }
    2104      102349 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    2105       87803 :     if (dx < dy) break;
    2106       73009 :     if (gc_needed(av2,1))
    2107             :     {
    2108          28 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    2109          28 :       gc_slice(av2,x,dx+1);
    2110             :     }
    2111             :   }
    2112       14794 :   return gc_bool(av, dx < 0);
    2113             : }
    2114             : 
    2115             : /* T either NULL or a t_POL. */
    2116             : GEN
    2117      627554 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    2118             : {
    2119      627554 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    2120      627554 :   pari_sp av = avma, av2;
    2121             :   GEN y_lead;
    2122             : 
    2123      627554 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    2124      627554 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2125             :   /* if monic, no point in using pseudo-division */
    2126      627554 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    2127      579241 :   dx = degpol(x);
    2128      579240 :   if (dx < dy) return RgX_copy(x);
    2129      579233 :   (void)new_chunk(2);
    2130      579233 :   x = RgX_recip_i(x)+2;
    2131      579235 :   y = RgX_recip_i(y)+2;
    2132             :   /* pay attention to sparse divisors */
    2133     1951521 :   for (i = 1; i <= dy; i++)
    2134     1372286 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2135      579235 :   dz = dx-dy; p = dz+1;
    2136      579235 :   av2 = avma;
    2137             :   for (;;)
    2138             :   {
    2139     1246540 :     gel(x,0) = gneg(gel(x,0)); p--;
    2140     4192131 :     for (i=1; i<=dy; i++)
    2141             :     {
    2142     2945609 :       GEN c = gmul(y_lead, gel(x,i));
    2143     2945525 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    2144     2945579 :       gel(x,i) = rem(c, T);
    2145             :     }
    2146     2182455 :     for (   ; i<=dx; i++)
    2147             :     {
    2148      935921 :       GEN c = gmul(y_lead, gel(x,i));
    2149      935908 :       gel(x,i) = rem(c, T);
    2150             :     }
    2151     1289342 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    2152     1246533 :     if (dx < dy) break;
    2153      667305 :     if (gc_needed(av2,1))
    2154             :     {
    2155           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    2156           0 :       gc_slice(av2,x,dx+1);
    2157             :     }
    2158             :   }
    2159      579228 :   if (dx < 0) return pol_0(vx);
    2160      579067 :   lx = dx+3; x -= 2;
    2161      579067 :   x[0] = evaltyp(t_POL) | _evallg(lx);
    2162      579067 :   x[1] = evalsigne(1) | evalvarn(vx);
    2163      579067 :   x = RgX_recip_i(x);
    2164      579070 :   if (p)
    2165             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    2166       28676 :     GEN t = y_lead;
    2167       28676 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    2168           0 :       t = RgXQ_powu(t, p, T);
    2169             :     else
    2170       28676 :       t = gpowgs(t, p);
    2171       88381 :     for (i=2; i<lx; i++)
    2172             :     {
    2173       59705 :       GEN c = gmul(gel(x,i), t);
    2174       59705 :       gel(x,i) = rem(c,T);
    2175             :     }
    2176       28676 :     if (!T) return gc_upto(av, x);
    2177             :   }
    2178      550394 :   return gc_GEN(av, x);
    2179             : }
    2180             : 
    2181             : GEN
    2182      627554 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2183             : 
    2184             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2185             : GEN
    2186       12056 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2187             : {
    2188       12056 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2189       12056 :   pari_sp av = avma, av2;
    2190             :   GEN z, r, ypow, y_lead;
    2191             : 
    2192       12056 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2193       12056 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2194       12056 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2195        7561 :   dx = degpol(x);
    2196        7561 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2197        7561 :   if (dx == dy)
    2198             :   {
    2199          98 :     GEN x_lead = gel(x,lg(x)-1);
    2200          98 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2201          98 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2202          98 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2203          98 :     *ptr = gc_upto(av, r); return scalarpol(x_lead, vx);
    2204             :   }
    2205        7463 :   (void)new_chunk(2);
    2206        7463 :   x = RgX_recip_i(x)+2;
    2207        7463 :   y = RgX_recip_i(y)+2;
    2208             :   /* pay attention to sparse divisors */
    2209       39000 :   for (i = 1; i <= dy; i++)
    2210       31537 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2211        7463 :   dz = dx-dy; p = dz+1;
    2212        7463 :   lz = dz+3;
    2213        7463 :   z = cgetg(lz, t_POL);
    2214        7463 :   z[1] = evalsigne(1) | evalvarn(vx);
    2215       28640 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2216        7463 :   ypow = new_chunk(dz+1);
    2217        7463 :   gel(ypow,0) = gen_1;
    2218        7463 :   gel(ypow,1) = y_lead;
    2219       13714 :   for (i=2; i<=dz; i++)
    2220             :   {
    2221        6251 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2222        6251 :     gel(ypow,i) = rem(c,T);
    2223             :   }
    2224        7463 :   av2 = avma;
    2225        7463 :   for (iz=2;;)
    2226             :   {
    2227       15605 :     p--;
    2228       15605 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2229       69915 :     for (i=1; i<=dy; i++)
    2230             :     {
    2231       54310 :       GEN c = gmul(y_lead, gel(x,i));
    2232       54310 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2233       54310 :       gel(x,i) = rem(c, T);
    2234             :     }
    2235       41058 :     for (   ; i<=dx; i++)
    2236             :     {
    2237       25453 :       GEN c = gmul(y_lead, gel(x,i));
    2238       25453 :       gel(x,i) = rem(c,T);
    2239             :     }
    2240       15605 :     x++; dx--;
    2241       21177 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2242       15605 :     if (dx < dy) break;
    2243        8142 :     if (gc_needed(av2,1))
    2244             :     {
    2245           0 :       GEN X = x-2;
    2246           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2247           0 :       X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
    2248           0 :       (void)gc_all(av2,2, &X, &z); x = X+2;
    2249             :     }
    2250             :   }
    2251       14008 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2252        7463 :   if (dx < 0)
    2253         182 :     x = pol_0(vx);
    2254             :   else
    2255             :   {
    2256        7281 :     lx = dx+3; x -= 2;
    2257        7281 :     x[0] = evaltyp(t_POL) | _evallg(lx);
    2258        7281 :     x[1] = evalsigne(1) | evalvarn(vx);
    2259        7281 :     x = RgX_recip_i(x);
    2260             :   }
    2261        7463 :   z = RgX_recip_i(z);
    2262        7463 :   r = x;
    2263        7463 :   if (p)
    2264             :   {
    2265        3339 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2266        3339 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2267             :   }
    2268        7463 :   *ptr = r; return gc_all(av, 2, &z, ptr);
    2269             : }
    2270             : GEN
    2271       11818 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2272       11818 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2273             : 
    2274             : GEN
    2275           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2276           0 : { return RgXQX_red(RgX_mul(x,y), T); }
    2277             : GEN
    2278   748054547 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
    2279             : GEN
    2280      141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
    2281             : GEN
    2282       26320 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
    2283             : GEN
    2284          35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
    2285             : GEN
    2286         133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
    2287             : 
    2288             : GEN
    2289           0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
    2290             : 
    2291             : GEN
    2292           0 : RgXQX_powers(GEN P, long n, GEN T)
    2293             : {
    2294           0 :   GEN v = cgetg(n+2, t_VEC);
    2295             :   long i;
    2296           0 :   gel(v, 1) = pol_1(varn(T));
    2297           0 :   if (n==0) return v;
    2298           0 :   gel(v, 2) = gcopy(P);
    2299           0 :   for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
    2300           0 :   return v;
    2301             : }
    2302             : 
    2303             : static GEN
    2304      623028 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
    2305             : static GEN
    2306           0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
    2307             : static GEN
    2308       67312 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2309             : static GEN
    2310       88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
    2311             : static GEN
    2312      276613 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2313             : static GEN
    2314     1039953 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2315             : static GEN
    2316     1018589 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2317             : static GEN
    2318        1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2319             : static GEN
    2320      723917 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
    2321             : 
    2322             : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
    2323             :               _mul, _sqr, _one, _zero };
    2324             : 
    2325             : GEN
    2326           0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
    2327             : {
    2328           0 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
    2329             : }
    2330             : 
    2331             : GEN
    2332      417191 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2333             : {
    2334      417191 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2335      417189 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
    2336             : }
    2337             : 
    2338             : /* mod X^n */
    2339             : struct modXn {
    2340             :   long v; /* varn(X) */
    2341             :   long n;
    2342             : } ;
    2343             : static GEN
    2344       11893 : _sqrXn(void *data, GEN x) {
    2345       11893 :   struct modXn *S = (struct modXn*)data;
    2346       11893 :   return RgXn_sqr(x, S->n);
    2347             : }
    2348             : static GEN
    2349        4528 : _mulXn(void *data, GEN x, GEN y) {
    2350        4528 :   struct modXn *S = (struct modXn*)data;
    2351        4528 :   return RgXn_mul(x,y, S->n);
    2352             : }
    2353             : static GEN
    2354        1939 : _oneXn(void *data) {
    2355        1939 :   struct modXn *S = (struct modXn*)data;
    2356        1939 :   return pol_1(S->v);
    2357             : }
    2358             : static GEN
    2359           0 : _zeroXn(void *data) {
    2360           0 :   struct modXn *S = (struct modXn*)data;
    2361           0 :   return pol_0(S->v);
    2362             : }
    2363             : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
    2364             :                                           _oneXn, _zeroXn };
    2365             : 
    2366             : GEN
    2367         483 : RgXn_powers(GEN x, long m, long n)
    2368             : {
    2369         483 :   long d = degpol(x);
    2370         483 :   int use_sqr = (d<<1) >= n;
    2371             :   struct modXn S;
    2372         483 :   S.v = varn(x); S.n = n;
    2373         483 :   return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
    2374             : }
    2375             : 
    2376             : GEN
    2377        2286 : RgXn_powu_i(GEN x, ulong m, long n)
    2378             : {
    2379             :   struct modXn S;
    2380             :   long v;
    2381        2286 :   if (n == 0) return x;
    2382        2286 :   v = RgX_valrem(x, &x);
    2383        2286 :   if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
    2384        2265 :   S.v = varn(x); S.n = n;
    2385        2265 :   x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2386        2265 :   if (v) x = RgX_shift_shallow(x, m * v);
    2387        2265 :   return x;
    2388             : }
    2389             : GEN
    2390           0 : RgXn_powu(GEN x, ulong m, long n)
    2391             : {
    2392             :   pari_sp av;
    2393           0 :   if (n == 0) return gcopy(x);
    2394           0 :   av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
    2395             : }
    2396             : 
    2397             : GEN
    2398         966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
    2399             : {
    2400             :   struct modXn S;
    2401         966 :   S.v = varn(gel(x,2)); S.n = n;
    2402         966 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
    2403             : }
    2404             : 
    2405             : GEN
    2406           0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
    2407             : {
    2408           0 :   int use_sqr = 2*degpol(x) >= n;
    2409             :   struct modXn S;
    2410           0 :   S.v = varn(x); S.n = n;
    2411           0 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2412             : }
    2413             : 
    2414             : /* Q(x) mod t^n, x in R[t], n >= 1 */
    2415             : GEN
    2416        5159 : RgXn_eval(GEN Q, GEN x, long n)
    2417             : {
    2418        5159 :   long d = degpol(x);
    2419             :   int use_sqr;
    2420             :   struct modXn S;
    2421        5159 :   if (d == 1 && isrationalzero(gel(x,2)))
    2422             :   {
    2423        5152 :     GEN y = RgX_unscale(Q, gel(x,3));
    2424        5152 :     setvarn(y, varn(x)); return y;
    2425             :   }
    2426           7 :   S.v = varn(x);
    2427           7 :   S.n = n;
    2428           7 :   use_sqr = (d<<1) >= n;
    2429           7 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2430             : }
    2431             : 
    2432             : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
    2433             : static GEN
    2434     2614673 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2435             : {
    2436     2614673 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2437     2614673 :   return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
    2438             : }
    2439             : 
    2440             : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
    2441             : static GEN
    2442          14 : RgXn_sqrhigh(GEN f, long n2, long n)
    2443             : {
    2444          14 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2445          14 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2446             : }
    2447             : 
    2448             : static GEN
    2449     2188945 : RgXn_div_gen(GEN g, GEN f, long e)
    2450             : {
    2451             :   pari_sp av;
    2452             :   ulong mask;
    2453             :   GEN W, a;
    2454     2188945 :   long v = varn(f), n = 1;
    2455             : 
    2456     2188945 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2457     2188938 :   a = ginv(gel(f,2));
    2458     2188938 :   if (e == 1 && !g) return scalarpol(a, v);
    2459     2184052 :   else if (e == 2 && !g)
    2460             :   {
    2461             :     GEN b;
    2462      836935 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2463      246253 :     b = gneg(b);
    2464      246251 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2465      246251 :     return deg1pol(b, a, v);
    2466             :   }
    2467     1347117 :   av = avma;
    2468     1347117 :   W = scalarpol_shallow(a,v);
    2469     1347117 :   mask = quadratic_prec_mask(e);
    2470     3955833 :   while (mask > 1)
    2471             :   {
    2472             :     GEN u, fr;
    2473     2608716 :     long n2 = n;
    2474     2608716 :     n<<=1; if (mask & 1) n--;
    2475     2608716 :     mask >>= 1;
    2476     2608716 :     fr = RgXn_red_shallow(f, n);
    2477     2608716 :     if (mask>1 || !g)
    2478             :     {
    2479     1324736 :       u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2480     1324736 :       W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2481             :     }
    2482             :     else
    2483             :     {
    2484     1283980 :       GEN y = RgXn_mul(g, W, n), yt =  RgXn_red_shallow(y, n-n2);
    2485     1283980 :       u = RgXn_mul(yt, RgXn_mulhigh(fr,  W, n2, n), n-n2);
    2486     1283980 :       W = RgX_sub(y, RgX_shift_shallow(u, n2));
    2487             :     }
    2488     2608716 :     if (gc_needed(av,2))
    2489             :     {
    2490           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2491           0 :       W = gc_upto(av, W);
    2492             :     }
    2493             :   }
    2494     1347117 :   return W;
    2495             : }
    2496             : 
    2497             : static GEN
    2498         119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
    2499             : {
    2500             :   GEN r;
    2501         119 :   if (lgefint(p) == 3)
    2502             :   {
    2503         119 :     ulong pp = uel(p, 2);
    2504         119 :     if (pp == 2)
    2505          14 :       r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
    2506             :     else
    2507         105 :       r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
    2508             :   }
    2509             :   else
    2510           0 :     r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
    2511         119 :   return FpX_to_mod(r, p);
    2512             : }
    2513             : 
    2514             : static GEN
    2515           7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
    2516             : {
    2517           7 :   GEN r, T = RgX_to_FpX(pol, p);
    2518           7 :   if (signe(T) == 0) pari_err_OP("/", x, y);
    2519           7 :   r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
    2520           7 :   return FpXQX_to_mod(r, T, p);
    2521             : }
    2522             : 
    2523             : static GEN
    2524          91 : RgXn_inv_FpX(GEN x, long e, GEN p)
    2525             : {
    2526             :   GEN r;
    2527          91 :   if (lgefint(p) == 3)
    2528             :   {
    2529          91 :     ulong pp = uel(p, 2);
    2530          91 :     if (pp == 2)
    2531          28 :       r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
    2532             :     else
    2533          63 :       r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
    2534             :   }
    2535             :   else
    2536           0 :     r = FpXn_inv(RgX_to_FpX(x, p), e, p);
    2537          91 :   return FpX_to_mod(r, p);
    2538             : }
    2539             : 
    2540             : static GEN
    2541           0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
    2542             : {
    2543           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2544           0 :   if (signe(T) == 0) pari_err_OP("/", gen_1, x);
    2545           0 :   r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
    2546           0 :   return FpXQX_to_mod(r, T, p);
    2547             : }
    2548             : 
    2549             : static GEN
    2550      905057 : RgXn_inv_fast(GEN x, long e)
    2551             : {
    2552             :   GEN p, pol;
    2553             :   long pa;
    2554      905057 :   long t = RgX_type(x,&p,&pol,&pa);
    2555      905056 :   switch(t)
    2556             :   {
    2557          91 :     case t_INTMOD: return RgXn_inv_FpX(x, e, p);
    2558           0 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    2559           0 :                    return RgXn_inv_FpXQX(x, e, pol, p);
    2560      904965 :     default:       return NULL;
    2561             :   }
    2562             : }
    2563             : 
    2564             : static GEN
    2565     1284106 : RgXn_div_fast(GEN x, GEN y, long e)
    2566             : {
    2567             :   GEN p, pol;
    2568             :   long pa;
    2569     1284106 :   long t = RgX_type2(x,y,&p,&pol,&pa);
    2570     1284106 :   switch(t)
    2571             :   {
    2572         119 :     case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
    2573           7 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    2574           7 :                    return RgXn_div_FpXQX(x, y, e, pol, p);
    2575     1283980 :     default:       return NULL;
    2576             :   }
    2577             : }
    2578             : 
    2579             : GEN
    2580     1284106 : RgXn_div_i(GEN g, GEN f, long e)
    2581             : {
    2582     1284106 :   GEN h = RgXn_div_fast(g, f, e);
    2583     1284106 :   if (h) return h;
    2584     1283980 :   return RgXn_div_gen(g, f, e);
    2585             : }
    2586             : 
    2587             : GEN
    2588      584443 : RgXn_div(GEN g, GEN f, long e)
    2589             : {
    2590      584443 :   pari_sp av = avma;
    2591      584443 :   return gc_upto(av, RgXn_div_i(g, f, e));
    2592             : }
    2593             : 
    2594             : GEN
    2595      905057 : RgXn_inv_i(GEN f, long e)
    2596             : {
    2597      905057 :   GEN h = RgXn_inv_fast(f, e);
    2598      905056 :   if (h) return h;
    2599      904965 :   return RgXn_div_gen(NULL, f, e);
    2600             : }
    2601             : 
    2602             : GEN
    2603      807474 : RgXn_inv(GEN f, long e)
    2604             : {
    2605      807474 :   pari_sp av = avma;
    2606      807474 :   return gc_upto(av, RgXn_inv_i(f, e));
    2607             : }
    2608             : 
    2609             : /* intformal(x^n*S) / x^(n+1) */
    2610             : static GEN
    2611       56377 : RgX_integXn(GEN x, long n)
    2612      114357 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
    2613             : 
    2614             : GEN
    2615       52919 : RgXn_expint(GEN h, long e)
    2616             : {
    2617       52919 :   pari_sp av = avma, av2;
    2618       52919 :   long v = varn(h), n;
    2619       52919 :   GEN f = pol_1(v), g;
    2620             :   ulong mask;
    2621             : 
    2622       52919 :   if (!signe(h)) return f;
    2623       50427 :   g = pol_1(v);
    2624       50427 :   n = 1; mask = quadratic_prec_mask(e);
    2625       50427 :   av2 = avma;
    2626       56377 :   for (;mask>1;)
    2627             :   {
    2628             :     GEN u, w;
    2629       56377 :     long n2 = n;
    2630       56377 :     n<<=1; if (mask & 1) n--;
    2631       56377 :     mask >>= 1;
    2632       56377 :     u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
    2633       56377 :     u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
    2634       56377 :     w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
    2635       56377 :     f = RgX_add(f, RgX_shift_shallow(w, n2));
    2636       56378 :     if (mask<=1) break;
    2637        5950 :     u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
    2638        5950 :     g = RgX_sub(g, RgX_shift_shallow(u, n2));
    2639        5950 :     if (gc_needed(av2,2))
    2640             :     {
    2641           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
    2642           0 :       (void)gc_all(av2, 2, &f, &g);
    2643             :     }
    2644             :   }
    2645       50428 :   return gc_upto(av, f);
    2646             : }
    2647             : 
    2648             : GEN
    2649           0 : RgXn_exp(GEN h, long e)
    2650             : {
    2651           0 :   long d = degpol(h);
    2652           0 :   if (d < 0) return pol_1(varn(h));
    2653           0 :   if (!d || !gequal0(gel(h,2)))
    2654           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2655           0 :   return RgXn_expint(RgX_deriv(h), e);
    2656             : }
    2657             : 
    2658             : GEN
    2659         154 : RgXn_reverse(GEN f, long e)
    2660             : {
    2661         154 :   pari_sp av = avma, av2;
    2662             :   ulong mask;
    2663             :   GEN fi, a, df, W, an;
    2664         154 :   long v = varn(f), n=1;
    2665         154 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2666           0 :     pari_err_INV("serreverse",f);
    2667         154 :   fi = ginv(gel(f,3));
    2668         154 :   a = deg1pol_shallow(fi,gen_0,v);
    2669         154 :   if (e <= 2) return gc_GEN(av, a);
    2670         133 :   W = scalarpol(fi,v);
    2671         133 :   df = RgX_deriv(f);
    2672         133 :   mask = quadratic_prec_mask(e);
    2673         133 :   av2 = avma;
    2674         616 :   for (;mask>1;)
    2675             :   {
    2676             :     GEN u, fa, fr;
    2677         483 :     long n2 = n, rt;
    2678         483 :     n<<=1; if (mask & 1) n--;
    2679         483 :     mask >>= 1;
    2680         483 :     fr = RgXn_red_shallow(f, n);
    2681         483 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2682         483 :     an = RgXn_powers(a, rt, n);
    2683         483 :     if (n>1)
    2684             :     {
    2685         483 :       long n4 = (n2+1)>>1;
    2686         483 :       GEN dfr = RgXn_red_shallow(df, n2);
    2687         483 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2688         483 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2689         483 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2690             :     }
    2691         483 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2692         483 :     fa = RgX_shift(fa, -n2);
    2693         483 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2694         483 :     if (gc_needed(av2,2))
    2695             :     {
    2696           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2697           0 :       (void)gc_all(av2, 2, &a, &W);
    2698             :     }
    2699             :   }
    2700         133 :   return gc_upto(av, a);
    2701             : }
    2702             : 
    2703             : GEN
    2704           7 : RgXn_sqrt(GEN h, long e)
    2705             : {
    2706           7 :   pari_sp av = avma, av2;
    2707           7 :   long v = varn(h), n = 1;
    2708           7 :   GEN f = scalarpol(gen_1, v), df = f;
    2709           7 :   ulong mask = quadratic_prec_mask(e);
    2710           7 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2711           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2712           7 :   av2 = avma;
    2713             :   while(1)
    2714           7 :   {
    2715          14 :     long n2 = n, m;
    2716             :     GEN g;
    2717          14 :     n<<=1; if (mask & 1) n--;
    2718          14 :     mask >>= 1;
    2719          14 :     m = n-n2;
    2720          14 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2721          14 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2722          14 :     if (mask==1) return gc_upto(av, f);
    2723           7 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2724           7 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2725           7 :     if (gc_needed(av2,2))
    2726             :     {
    2727           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2728           0 :       (void)gc_all(av2, 2, &f, &df);
    2729             :     }
    2730             :   }
    2731             : }
    2732             : 
    2733             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2734             : GEN
    2735      116695 : RgXQ_powu(GEN x, ulong n, GEN T)
    2736             : {
    2737      116695 :   pari_sp av = avma;
    2738             : 
    2739      116695 :   if (!n) return pol_1(varn(x));
    2740       68452 :   if (n == 1) return RgX_copy(x);
    2741       19824 :   x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
    2742       19824 :   return gc_GEN(av, x);
    2743             : }
    2744             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2745             : GEN
    2746      102160 : RgXQ_pow(GEN x, GEN n, GEN T)
    2747             : {
    2748             :   pari_sp av;
    2749      102160 :   long s = signe(n);
    2750             : 
    2751      102160 :   if (!s) return pol_1(varn(x));
    2752      102160 :   if (is_pm1(n) == 1)
    2753       88307 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2754       13853 :   av = avma;
    2755       13853 :   if (s < 0) x = RgXQ_inv(x, T);
    2756       13853 :   x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
    2757       13853 :   return gc_GEN(av, x);
    2758             : }
    2759             : static GEN
    2760      200577 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
    2761             : static GEN
    2762      111154 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
    2763             : 
    2764             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2765             : GEN
    2766       11746 : ZXQ_powers(GEN x, long l, GEN T)
    2767             : {
    2768       11746 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2769       11746 :   return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
    2770             : }
    2771             : 
    2772             : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
    2773             : GEN
    2774      168513 : ZXQ_powu(GEN x, ulong n, GEN T)
    2775             : {
    2776      168513 :   pari_sp av = avma;
    2777             : 
    2778      168513 :   if (!n) return pol_1(varn(x));
    2779      168513 :   if (n == 1) return ZX_copy(x);
    2780      113877 :   x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
    2781      113878 :   return gc_GEN(av, x);
    2782             : }
    2783             : 
    2784             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2785             : GEN
    2786        8554 : RgXQ_powers(GEN x, long l, GEN T)
    2787             : {
    2788        8554 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2789        8554 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2790             : }
    2791             : 
    2792             : GEN
    2793        7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
    2794             : {
    2795        7118 :   return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
    2796             : }
    2797             : 
    2798             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2799             : GEN
    2800        9191 : QXQ_powers(GEN a, long n, GEN T)
    2801             : {
    2802             :   GEN den, v;
    2803        9191 :   if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
    2804        9163 :   v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
    2805             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2806        9163 :   if (den)
    2807             :   { /* restore denominators */
    2808        5620 :     GEN d = den;
    2809             :     long i;
    2810        5620 :     gel(v,2) = a;
    2811       27280 :     for (i=3; i<=n+1; i++) {
    2812       21660 :       d = mulii(d,den);
    2813       21660 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2814             :     }
    2815             :   }
    2816        9163 :   return v;
    2817             : }
    2818             : 
    2819             : static GEN
    2820        3241 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2821             : {
    2822        3241 :   long l, i, m = 0;
    2823             :   GEN dz, z;
    2824        3241 :   GEN V = cgetg_copy(v, &l);
    2825       10794 :   for (i = imin; i < l; i++)
    2826             :   {
    2827        7553 :     GEN c = gel(v, i);
    2828        7553 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2829             :   }
    2830        3241 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2831        3479 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2832       10794 :   for (i = imin; i < l; i++)
    2833             :   {
    2834        7553 :     GEN c = gel(v,i);
    2835        7553 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2836        7553 :     gel(V,i) = c;
    2837             :   }
    2838        3241 :   return V;
    2839             : }
    2840             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2841             : GEN
    2842        3003 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2843        3003 : { return do_QXQ_eval(v, 1, a, T); }
    2844             : 
    2845             : GEN
    2846         238 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
    2847         238 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2848             : 
    2849             : GEN
    2850        2940 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2851             : {
    2852        2940 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2853             : }
    2854             : 
    2855             : GEN
    2856        5144 : RgXQ_norm(GEN x, GEN T)
    2857             : {
    2858             :   pari_sp av;
    2859        5144 :   long dx = degpol(x);
    2860             :   GEN L, y;
    2861        5144 :   if (degpol(T)==0) return gpowgs(x,0);
    2862        5137 :   av = avma; y = resultant(T, x);
    2863        5137 :   L = leading_coeff(T);
    2864        5137 :   if (gequal1(L) || !signe(x)) return y;
    2865           0 :   return gc_upto(av, gdiv(y, gpowgs(L, dx)));
    2866             : }
    2867             : 
    2868             : GEN
    2869         476 : RgXQ_trace(GEN x, GEN T)
    2870             : {
    2871         476 :   pari_sp av = avma;
    2872             :   GEN dT, z;
    2873             :   long n;
    2874         476 :   if (degpol(T)==0) return gmulgs(x,0);
    2875         469 :   dT = RgX_deriv(T); n = degpol(dT);
    2876         469 :   z = RgXQ_mul(x, dT, T);
    2877         469 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2878         420 :   return gc_upto(av, gdiv(gel(z,2+n), gel(T,3+n)));
    2879             : }
    2880             : 
    2881             : GEN
    2882     2788358 : RgX_blocks(GEN P, long n, long m)
    2883             : {
    2884     2788358 :   GEN z = cgetg(m+1,t_VEC);
    2885     2788358 :   long i,j, k=2, l = lg(P);
    2886     8558652 :   for(i=1; i<=m; i++)
    2887             :   {
    2888     5770294 :     GEN zi = cgetg(n+2,t_POL);
    2889     5770294 :     zi[1] = P[1];
    2890     5770294 :     gel(z,i) = zi;
    2891    17014348 :     for(j=2; j<n+2; j++)
    2892    11244054 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2893     5770294 :     zi = RgX_renormalize_lg(zi, n+2);
    2894             :   }
    2895     2788358 :   return z;
    2896             : }
    2897             : 
    2898             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2899             : void
    2900    49737995 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2901             : {
    2902    49737995 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2903             :   GEN p0, p1;
    2904             : 
    2905    49738388 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2906             : 
    2907    49645111 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2908    49645111 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2909    49648573 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2910   146734466 :   for (i=0; i<n1; i++)
    2911             :   {
    2912    97084642 :     p0[2+i] = p[2+(i<<1)];
    2913    97084642 :     p1[2+i] = p[3+(i<<1)];
    2914             :   }
    2915    49649824 :   if (n1 != n0)
    2916    18009339 :     p0[2+i] = p[2+(i<<1)];
    2917    49649824 :   *pe = normalizepol(p0);
    2918    49654513 :   *po = normalizepol(p1);
    2919             : }
    2920             : 
    2921             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2922             : GEN
    2923       43631 : RgX_splitting(GEN p, long k)
    2924             : {
    2925       43631 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2926             :   GEN r;
    2927             : 
    2928       43631 :   m = n/k;
    2929       43631 :   r = cgetg(k+1,t_VEC);
    2930      234059 :   for(i=1; i<=k; i++)
    2931             :   {
    2932      190428 :     gel(r,i) = cgetg(m+3, t_POL);
    2933      190428 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2934             :   }
    2935      571991 :   for (j=1, i=0, l=2; i<=n; i++)
    2936             :   {
    2937      528360 :     gmael(r,j,l) = gel(p,2+i);
    2938      528360 :     if (j==k) { j=1; l++; } else j++;
    2939             :   }
    2940      234059 :   for(i=1; i<=k; i++)
    2941      190428 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2942       43631 :   return r;
    2943             : }
    2944             : 
    2945             : /*******************************************************************/
    2946             : /*                                                                 */
    2947             : /*                        Kronecker form                           */
    2948             : /*                                                                 */
    2949             : /*******************************************************************/
    2950             : 
    2951             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2952             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2953             :  * normalized coefficients */
    2954             : GEN
    2955      186168 : Kronecker_to_mod(GEN z, GEN T)
    2956             : {
    2957      186168 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2958      186168 :   GEN x, t = cgetg(N,t_POL);
    2959      186168 :   t[1] = T[1];
    2960      186168 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2961      186168 :   x[1] = z[1];
    2962      186168 :   T = RgX_copy(T);
    2963     1399389 :   for (i=2; i<lx+2; i++, z+= N-2)
    2964             :   {
    2965     5761158 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2966     1213221 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2967             :   }
    2968      186168 :   N = (l-2) % (N-2) + 2;
    2969      472176 :   for (j=2; j<N; j++) t[j] = z[j];
    2970      186168 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2971      186168 :   return normalizepol_lg(x, i+1);
    2972             : }
    2973             : 
    2974             : /*******************************************************************/
    2975             : /*                                                                 */
    2976             : /*                        Domain detection                         */
    2977             : /*                                                                 */
    2978             : /*******************************************************************/
    2979             : 
    2980             : static GEN
    2981      609812 : RgX_mul_FpX(GEN x, GEN y, GEN p)
    2982             : {
    2983      609812 :   pari_sp av = avma;
    2984             :   GEN r;
    2985      609812 :   if (lgefint(p) == 3)
    2986             :   {
    2987      558456 :     ulong pp = uel(p, 2);
    2988      558456 :     r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
    2989             :                                   RgX_to_Flx(y, pp), pp));
    2990             :   }
    2991             :   else
    2992       51356 :     r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2993      609812 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2994      544614 :   return gc_upto(av, FpX_to_mod(r, p));
    2995             : }
    2996             : 
    2997             : static GEN
    2998         616 : zero_FpXQX_mod(GEN pol, GEN p, long v)
    2999             : {
    3000         616 :   GEN r = cgetg(3,t_POL);
    3001         616 :   r[1] = evalvarn(v);
    3002         616 :   gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
    3003         616 :   return r;
    3004             : }
    3005             : 
    3006             : static GEN
    3007        2079 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3008             : {
    3009        2079 :   pari_sp av = avma;
    3010             :   long dT;
    3011             :   GEN kx, ky, r;
    3012        2079 :   GEN T = RgX_to_FpX(pol, p);
    3013        2079 :   if (signe(T)==0) pari_err_OP("*", x, y);
    3014        2072 :   dT = degpol(T);
    3015        2072 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    3016        2072 :   ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
    3017        2072 :   r = FpX_mul(kx, ky, p);
    3018        2072 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3019        1463 :   return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    3020             : }
    3021             : 
    3022             : static GEN
    3023      455141 : RgX_liftred(GEN x, GEN T)
    3024      455141 : { return RgXQX_red(liftpol_shallow(x), T); }
    3025             : 
    3026             : static GEN
    3027      171091 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
    3028             : {
    3029      171091 :   pari_sp av = avma;
    3030      171091 :   long dT = degpol(T);
    3031      171091 :   GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
    3032             :                  RgXX_to_Kronecker(RgX_liftred(y, T), dT));
    3033      171091 :   return gc_upto(av, Kronecker_to_mod(r, T));
    3034             : }
    3035             : 
    3036             : static GEN
    3037        1407 : RgX_sqr_FpX(GEN x, GEN p)
    3038             : {
    3039        1407 :   pari_sp av = avma;
    3040             :   GEN r;
    3041        1407 :   if (lgefint(p) == 3)
    3042             :   {
    3043        1203 :     ulong pp = uel(p, 2);
    3044        1203 :     r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
    3045             :   }
    3046             :   else
    3047         204 :     r = FpX_sqr(RgX_to_FpX(x, p), p);
    3048        1407 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3049        1274 :   return gc_upto(av, FpX_to_mod(r, p));
    3050             : }
    3051             : 
    3052             : static GEN
    3053         196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
    3054             : {
    3055         196 :   pari_sp av = avma;
    3056             :   long dT;
    3057         196 :   GEN kx, r, T = RgX_to_FpX(pol, p);
    3058         196 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3059         189 :   dT = degpol(T);
    3060         189 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    3061         189 :   r = FpX_sqr(kx, p);
    3062         189 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3063         189 :   return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    3064             : }
    3065             : 
    3066             : static GEN
    3067       13425 : RgX_sqr_QXQX(GEN x, GEN T)
    3068             : {
    3069       13425 :   pari_sp av = avma;
    3070       13425 :   long dT = degpol(T);
    3071       13425 :   GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
    3072       13425 :   return gc_upto(av, Kronecker_to_mod(r, T));
    3073             : }
    3074             : 
    3075             : static GEN
    3076       68320 : RgX_rem_FpX(GEN x, GEN y, GEN p)
    3077             : {
    3078       68320 :   pari_sp av = avma;
    3079             :   GEN r;
    3080       68320 :   if (lgefint(p) == 3)
    3081             :   {
    3082       51684 :     ulong pp = uel(p, 2);
    3083       51684 :     r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
    3084             :                                   RgX_to_Flx(y, pp), pp));
    3085             :   }
    3086             :   else
    3087       16636 :     r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3088       68320 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3089       31024 :   return gc_upto(av, FpX_to_mod(r, p));
    3090             : }
    3091             : 
    3092             : static GEN
    3093       49767 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
    3094             : {
    3095       49767 :   pari_sp av = avma;
    3096             :   GEN r;
    3097       49767 :   r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3098       49767 :   return gc_GEN(av, QXQX_to_mod_shallow(r, T));
    3099             : }
    3100             : static GEN
    3101          70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3102             : {
    3103          70 :   pari_sp av = avma;
    3104             :   GEN r;
    3105          70 :   GEN T = RgX_to_FpX(pol, p);
    3106          70 :   if (signe(T) == 0) pari_err_OP("%", x, y);
    3107          63 :   if (lgefint(p) == 3)
    3108             :   {
    3109          55 :     ulong pp = uel(p, 2);
    3110          55 :     GEN Tp = ZX_to_Flx(T, pp);
    3111          55 :     r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
    3112             :                               RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3113             :   }
    3114             :   else
    3115           8 :     r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3116          63 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3117          56 :   return gc_upto(av, FpXQX_to_mod(r, T, p));
    3118             : }
    3119             : 
    3120             : static GEN
    3121   121131203 : RgX_mul_fast(GEN x, GEN y)
    3122             : {
    3123             :   GEN p, pol;
    3124             :   long pa;
    3125   121131203 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3126   121131480 :   switch(t)
    3127             :   {
    3128    67189526 :     case t_INT:    return ZX_mul(x,y);
    3129     3356236 :     case t_FRAC:   return QX_mul(x,y);
    3130      105127 :     case t_FFELT:  return FFX_mul(x, y, pol);
    3131      609812 :     case t_INTMOD: return RgX_mul_FpX(x, y, p);
    3132      171091 :     case RgX_type_code(t_POLMOD, t_INT):
    3133             :     case RgX_type_code(t_POLMOD, t_FRAC):
    3134      171091 :                    return RgX_mul_QXQX(x, y, pol);
    3135        2060 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3136        2060 :                    return RgX_mul_FpXQX(x, y, pol, p);
    3137    49697628 :     default:       return NULL;
    3138             :   }
    3139             : }
    3140             : static GEN
    3141     1384542 : RgX_sqr_fast(GEN x)
    3142             : {
    3143             :   GEN p, pol;
    3144             :   long pa;
    3145     1384542 :   long t = RgX_type(x,&p,&pol,&pa);
    3146     1384545 :   switch(t)
    3147             :   {
    3148     1245800 :     case t_INT:    return ZX_sqr(x);
    3149       83261 :     case t_FRAC:   return QX_sqr(x);
    3150        2499 :     case t_FFELT:  return FFX_sqr(x, pol);
    3151        1407 :     case t_INTMOD: return RgX_sqr_FpX(x, p);
    3152       13425 :     case RgX_type_code(t_POLMOD, t_INT):
    3153             :     case RgX_type_code(t_POLMOD, t_FRAC):
    3154       13425 :                    return RgX_sqr_QXQX(x, pol);
    3155         192 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3156         192 :                    return RgX_sqr_FpXQX(x, pol, p);
    3157       37961 :     default:       return NULL;
    3158             :   }
    3159             : }
    3160             : 
    3161             : static GEN
    3162    14689577 : RgX_rem_fast(GEN x, GEN y)
    3163             : {
    3164             :   GEN p, pol;
    3165             :   long pa;
    3166    14689577 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3167    14689685 :   switch(t)
    3168             :   {
    3169     3471688 :     case t_INT:    return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
    3170     1835477 :     case t_FRAC:   return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
    3171          84 :     case t_FFELT:  return FFX_rem(x, y, pol);
    3172       68320 :     case t_INTMOD: return RgX_rem_FpX(x, y, p);
    3173       49767 :     case RgX_type_code(t_POLMOD, t_INT):
    3174             :     case RgX_type_code(t_POLMOD, t_FRAC):
    3175       49767 :                    return RgX_rem_QXQX(x, y, pol);
    3176          60 :     case RgX_type_code(t_POLMOD, t_INTMOD):
    3177          60 :                    return RgX_rem_FpXQX(x, y, pol, p);
    3178     9264289 :     default:       return NULL;
    3179             :   }
    3180             : }
    3181             : 
    3182             : GEN
    3183   107398495 : RgX_mul(GEN x, GEN y)
    3184             : {
    3185   107398495 :   GEN z = RgX_mul_fast(x,y);
    3186   107398645 :   if (!z) z = RgX_mul_i(x,y);
    3187   107398209 :   return z;
    3188             : }
    3189             : 
    3190             : GEN
    3191     1372390 : RgX_sqr(GEN x)
    3192             : {
    3193     1372390 :   GEN z = RgX_sqr_fast(x);
    3194     1372391 :   if (!z) z = RgX_sqr_i(x);
    3195     1372391 :   return z;
    3196             : }
    3197             : 
    3198             : GEN
    3199    14689588 : RgX_rem(GEN x, GEN y)
    3200             : {
    3201    14689588 :   GEN z = RgX_rem_fast(x, y);
    3202    14689672 :   if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
    3203    14689621 :   return z;
    3204             : }
    3205             : 
    3206             : static GEN
    3207           0 : _RgX_mul(void* E, GEN x, GEN y)
    3208           0 : { (void) E; return RgX_mul(x, y); }
    3209             : 
    3210             : GEN
    3211           0 : RgXV_prod(GEN V)
    3212           0 : { return gen_product(V, NULL, &_RgX_mul); }
    3213             : 
    3214             : GEN
    3215    11061657 : RgXn_mul(GEN f, GEN g, long n)
    3216             : {
    3217    11061657 :   pari_sp av = avma;
    3218    11061657 :   GEN h = RgX_mul_fast(f,g);
    3219    11061654 :   if (!h) return RgXn_mul2(f,g,n);
    3220     1783460 :   if (degpol(h) < n) return h;
    3221      382975 :   return gc_GEN(av, RgXn_red_shallow(h, n));
    3222             : }
    3223             : 
    3224             : GEN
    3225       12152 : RgXn_sqr(GEN f, long n)
    3226             : {
    3227       12152 :   pari_sp av = avma;
    3228       12152 :   GEN g = RgX_sqr_fast(f);
    3229       12152 :   if (!g) return RgXn_sqr2(f,n);
    3230       11837 :   if (degpol(g) < n) return g;
    3231       10640 :   return gc_GEN(av, RgXn_red_shallow(g, n));
    3232             : }
    3233             : 
    3234             : /* (f*g) \/ x^n */
    3235             : GEN
    3236     2671064 : RgX_mulhigh_i(GEN f, GEN g, long n)
    3237             : {
    3238     2671064 :   GEN h = RgX_mul_fast(f,g);
    3239     2671064 :   return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
    3240             : }
    3241             : 
    3242             : /* (f*g) \/ x^n */
    3243             : GEN
    3244           0 : RgX_sqrhigh_i(GEN f, long n)
    3245             : {
    3246           0 :   GEN h = RgX_sqr_fast(f);
    3247           0 :   return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
    3248             : }

Generated by: LCOV version 1.16