Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - RgX.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 1610 1770 91.0 %
Date: 2024-03-18 08:03:28 Functions: 199 218 91.3 %
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    19605090 : brent_kung_optpow(long d, long n, long m)
      32             : {
      33             :   long p, r;
      34    19605090 :   long pold=1, rold=n*(d-1);
      35    92264160 :   for(p=2; p<=d; p++)
      36             :   {
      37    72659070 :     r = m*(p-1) + n*((d-1)/p);
      38    72659070 :     if (r<rold) { pold=p; rold=r; }
      39             :   }
      40    19605090 :   return pold;
      41             : }
      42             : 
      43             : static GEN
      44    13709059 : 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    13709059 :   pari_sp av = avma;
      48             :   long i;
      49    13709059 :   GEN z = cmul(E,P,a,ff->one(E));
      50    13679447 :   if (!z) z = gen_0;
      51    69489093 :   for (i=1; i<=n; i++)
      52             :   {
      53    55808074 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      54    55828534 :     if (t) {
      55    43045417 :       z = ff->add(E, z, t);
      56    43014858 :       if (gc_needed(av,2)) z = gerepileupto(av, z);
      57             :     }
      58             :   }
      59    13681019 :   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    12266302 : 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    12266302 :   pari_sp av = avma;
      72    12266302 :   long l = lg(V)-1;
      73             :   GEN z, u;
      74             : 
      75    12266302 :   if (d < 0) return ff->zero(E);
      76    11660406 :   if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      77     1030209 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      78     1030209 :   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     1030209 :   d -= l;
      84     1030209 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      85     2047686 :   while (d >= l-1)
      86             :   {
      87     1016147 :     d -= l-1;
      88     1016147 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      89     1016129 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      90     1016125 :     if (gc_needed(av,2))
      91         105 :       z = gerepileupto(av, z);
      92             :   }
      93     1031539 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      94     1031555 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      95     1031553 :   return gerepileupto(av, ff->red(E,z));
      96             : }
      97             : 
      98             : GEN
      99      730866 : 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      730866 :   pari_sp av = avma;
     103             :   GEN z, V;
     104             :   long rtd;
     105      730866 :   if (d < 0) return ff->zero(E);
     106      729921 :   rtd = (long) sqrt((double)d);
     107      729921 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     108      729940 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     109      729931 :   return gerepileupto(av, z);
     110             : }
     111             : 
     112             : static GEN
     113     1193489 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     114             : static GEN
     115    15694499 : _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     1041758 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     120             : static GEN
     121      345072 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     122             : static GEN
     123     1217225 : _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      226034 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
     161             : {
     162      226034 :   pari_sp av = avma;
     163      226034 :   long i, d = degpol(P);
     164             :   GEN s;
     165      226034 :   if (signe(P)==0) return pol_0(varn(P));
     166      226034 :   s = gel(P, d+2);
     167      226034 :   if (d == 0) return gcopy(s);
     168     1079687 :   for (i = d-1; i >= 0; i--)
     169             :   {
     170      856936 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
     171      856936 :     if (gc_needed(av,1))
     172             :     {
     173          13 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
     174          13 :       s = gerepileupto(av, s);
     175             :     }
     176             :   }
     177      222751 :   return gerepileupto(av, s);
     178             : }
     179             : 
     180             : GEN
     181        1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
     182             : {
     183        1652 :   pari_sp av = avma;
     184        1652 :   long i, d = degpol(P), v = varn(A);
     185             :   GEN s;
     186        1652 :   if (signe(P)==0) return pol_0(v);
     187        1652 :   if (d == 0) return scalarpol(gel(P, d+2), v);
     188        1232 :   s = scalarpol_shallow(gel(P, d+2), v);
     189        4963 :   for (i = d-1; i >= 0; i--)
     190             :   {
     191        3731 :     GEN c = gel(P,i+2), b = gel(B,d+1-i);
     192        3731 :     s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
     193        3731 :     if (gc_needed(av,1))
     194             :     {
     195           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
     196           0 :       s = gerepileupto(av, s);
     197             :     }
     198             :   }
     199        1232 :   return gerepileupto(av, s);
     200             : }
     201             : 
     202             : const struct bb_algebra *
     203      154104 : get_Rg_algebra(void)
     204             : {
     205      154104 :   return &Rg_algebra;
     206             : }
     207             : 
     208             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     209             : 
     210             : static GEN
     211       30352 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     212             : {
     213             :   (void) E;
     214       30352 :   return RgX_divrem(x, y, r);
     215             : }
     216             : 
     217             : GEN
     218        6923 : RgX_digits(GEN x, GEN T)
     219             : {
     220        6923 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     221        6923 :   return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     222             : }
     223             : 
     224             : /*******************************************************************/
     225             : /*                                                                 */
     226             : /*                         RgX                                     */
     227             : /*                                                                 */
     228             : /*******************************************************************/
     229             : 
     230             : long
     231    24761439 : RgX_equal(GEN x, GEN y)
     232             : {
     233    24761439 :   long i = lg(x);
     234             : 
     235    24761439 :   if (i != lg(y)) return 0;
     236   106406540 :   for (i--; i > 1; i--)
     237    81949644 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     238    24456896 :   return 1;
     239             : }
     240             : 
     241             : /* Returns 1 in the base ring over which x is defined */
     242             : /* HACK: this also works for t_SER */
     243             : GEN
     244   156646093 : Rg_get_1(GEN x)
     245             : {
     246             :   GEN p, T;
     247   156646093 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     248   156646093 :   if (RgX_type_is_composite(tx))
     249    11621909 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     250   156646093 :   switch(tx)
     251             :   {
     252      795214 :     case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
     253         973 :     case t_PADIC: return cvtop(gen_1, p, lx);
     254        5481 :     case t_FFELT: return FF_1(T);
     255   155844425 :     default: return gen_1;
     256             :   }
     257             : }
     258             : /* Returns 0 in the base ring over which x is defined */
     259             : /* HACK: this also works for t_SER */
     260             : GEN
     261     6974331 : Rg_get_0(GEN x)
     262             : {
     263             :   GEN p, T;
     264     6974331 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     265     6974331 :   if (RgX_type_is_composite(tx))
     266       51170 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     267     6974331 :   switch(tx)
     268             :   {
     269         525 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     270          42 :     case t_PADIC: return zeropadic(p, lx);
     271         231 :     case t_FFELT: return FF_zero(T);
     272     6973533 :     default: return gen_0;
     273             :   }
     274             : }
     275             : 
     276             : GEN
     277        5502 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     278             : {
     279        5502 :   long i, n = degpol(P);
     280             :   GEN z, dz, dP;
     281        5502 :   if (n < 0) return gen_0;
     282        5502 :   P = Q_remove_denom(P, &dP);
     283        5502 :   z = gel(P,2); if (n == 0) return icopy(z);
     284        3248 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     285        3248 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     286        6720 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     287        3248 :   dz = mul_denom(dP, dV);
     288        3248 :   return dz? RgX_Rg_div(z, dz): z;
     289             : }
     290             : 
     291             : /* Return P(h * x), not memory clean */
     292             : GEN
     293       19156 : RgX_unscale(GEN P, GEN h)
     294             : {
     295       19156 :   long i, l = lg(P);
     296       19156 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     297       19156 :   Q[1] = P[1];
     298       19156 :   if (l == 2) return Q;
     299       19128 :   gel(Q,2) = gcopy(gel(P,2));
     300       40354 :   for (i=3; i<l; i++)
     301             :   {
     302       21226 :     hi = gmul(hi,h);
     303       21226 :     gel(Q,i) = gmul(gel(P,i), hi);
     304             :   }
     305       19128 :   return Q;
     306             : }
     307             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     308             : GEN
     309     1450722 : ZX_z_unscale(GEN P, long h)
     310             : {
     311     1450722 :   long i, l = lg(P);
     312     1450722 :   GEN Q = cgetg(l, t_POL);
     313     1450722 :   Q[1] = P[1];
     314     1450722 :   if (l == 2) return Q;
     315     1362288 :   gel(Q,2) = gel(P,2);
     316     1362288 :   if (l == 3) return Q;
     317     1335968 :   if (h == -1)
     318      229138 :     for (i = 3; i < l; i++)
     319             :     {
     320      190204 :       gel(Q,i) = negi(gel(P,i));
     321      190204 :       if (++i == l) break;
     322      140527 :       gel(Q,i) = gel(P,i);
     323             :     }
     324             :   else
     325             :   {
     326             :     GEN hi;
     327     1247357 :     gel(Q,3) = mulis(gel(P,3), h);
     328     1247357 :     hi = sqrs(h);
     329     5361087 :     for (i = 4; i < l; i++)
     330             :     {
     331     4113730 :       gel(Q,i) = mulii(gel(P,i), hi);
     332     4113730 :       if (i != l-1) hi = mulis(hi,h);
     333             :     }
     334             :   }
     335     1335968 :   return Q;
     336             : }
     337             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     338             : GEN
     339      999575 : ZX_unscale(GEN P, GEN h)
     340             : {
     341             :   long i, l;
     342             :   GEN Q, hi;
     343      999575 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     344         804 :   l = lg(P); Q = cgetg(l, t_POL);
     345         804 :   Q[1] = P[1];
     346         804 :   if (l == 2) return Q;
     347         804 :   gel(Q,2) = gel(P,2);
     348         804 :   if (l == 3) return Q;
     349         804 :   hi = h;
     350         804 :   gel(Q,3) = mulii(gel(P,3), hi);
     351        2531 :   for (i = 4; i < l; i++)
     352             :   {
     353        1727 :     hi = mulii(hi,h);
     354        1727 :     gel(Q,i) = mulii(gel(P,i), hi);
     355             :   }
     356         804 :   return Q;
     357             : }
     358             : /* P a ZX. Return P(x << n), not memory clean */
     359             : GEN
     360      216988 : ZX_unscale2n(GEN P, long n)
     361             : {
     362      216988 :   long i, ni = n, l = lg(P);
     363      216988 :   GEN Q = cgetg(l, t_POL);
     364      216988 :   Q[1] = P[1];
     365      216988 :   if (l == 2) return Q;
     366      216988 :   gel(Q,2) = gel(P,2);
     367      216988 :   if (l == 3) return Q;
     368      216988 :   gel(Q,3) = shifti(gel(P,3), ni);
     369      798898 :   for (i=4; i<l; i++)
     370             :   {
     371      581922 :     ni += n;
     372      581922 :     gel(Q,i) = shifti(gel(P,i), ni);
     373             :   }
     374      216976 :   return Q;
     375             : }
     376             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     377             : GEN
     378       12073 : ZX_unscale_div(GEN P, GEN h)
     379             : {
     380       12073 :   long i, l = lg(P);
     381       12073 :   GEN hi, Q = cgetg(l, t_POL);
     382       12073 :   Q[1] = P[1];
     383       12073 :   if (l == 2) return Q;
     384       12073 :   gel(Q,2) = diviiexact(gel(P,2), h);
     385       12073 :   if (l == 3) return Q;
     386       12073 :   gel(Q,3) = gel(P,3);
     387       12073 :   if (l == 4) return Q;
     388       12073 :   hi = h;
     389       12073 :   gel(Q,4) = mulii(gel(P,4), hi);
     390       61349 :   for (i=5; i<l; i++)
     391             :   {
     392       49276 :     hi = mulii(hi,h);
     393       49276 :     gel(Q,i) = mulii(gel(P,i), hi);
     394             :   }
     395       12073 :   return Q;
     396             : }
     397             : /* P(h*X) / h^k, assuming the result is a ZX */
     398             : GEN
     399         966 : ZX_unscale_divpow(GEN P, GEN h, long k)
     400             : {
     401         966 :   long i, j, l = lg(P);
     402         966 :   GEN H, Q = cgetg(l, t_POL);
     403         966 :   Q[1] = P[1]; if (l == 2) return Q;
     404         966 :   H = gpowers(h, maxss(k, l - 3 - k));
     405        3864 :   for (i = 2, j = k+1; j > 1 && i < l; i++)
     406        2898 :     gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
     407         966 :   if (i == l) return Q;
     408         966 :   gel(Q, i) = gel(P, i); i++;
     409        3542 :   for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
     410         966 :   return Q;
     411             : }
     412             : 
     413             : GEN
     414        5754 : RgXV_unscale(GEN x, GEN h)
     415             : {
     416        5754 :   if (isint1(h)) return RgX_copy(x);
     417       16510 :   pari_APPLY_same(RgX_unscale(gel(x,i), h));
     418             : }
     419             : 
     420             : /* Return h^degpol(P) P(x / h), not memory clean */
     421             : GEN
     422     4141604 : RgX_rescale(GEN P, GEN h)
     423             : {
     424     4141604 :   long i, l = lg(P);
     425     4141604 :   GEN Q = cgetg(l,t_POL), hi = h;
     426     4141599 :   gel(Q,l-1) = gel(P,l-1);
     427    10970592 :   for (i=l-2; i>=2; i--)
     428             :   {
     429    10967983 :     gel(Q,i) = gmul(gel(P,i), hi);
     430    10967814 :     if (i == 2) break;
     431     6828677 :     hi = gmul(hi,h);
     432             :   }
     433     4141746 :   Q[1] = P[1]; return Q;
     434             : }
     435             : 
     436             : GEN
     437        2401 : RgXV_rescale(GEN x, GEN h)
     438             : {
     439        2401 :   if (isint1(h)) return RgX_copy(x);
     440       16086 :   pari_APPLY_same(RgX_rescale(gel(x,i), h));
     441             : }
     442             : 
     443             : /* A(X^d) --> A(X) */
     444             : GEN
     445      959258 : RgX_deflate(GEN x0, long d)
     446             : {
     447             :   GEN z, y, x;
     448      959258 :   long i,id, dy, dx = degpol(x0);
     449      959256 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     450      371227 :   dy = dx/d;
     451      371227 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     452      371226 :   z = y + 2;
     453      371226 :   x = x0+ 2;
     454     1493827 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     455      371226 :   return y;
     456             : }
     457             : 
     458             : GEN
     459       18732 : RgX_homogenize(GEN P, long v)
     460             : {
     461             :   long i, l, d;
     462       18732 :   GEN Q = cgetg_copy(P, &l);
     463       18732 :   Q[1] = P[1]; d = l-3;
     464      165998 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
     465       18732 :   return Q;
     466             : }
     467             : 
     468             : /* F a t_RFRAC */
     469             : long
     470          98 : rfrac_deflate_order(GEN F)
     471             : {
     472          98 :   GEN N = gel(F,1), D = gel(F,2);
     473          98 :   long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
     474          98 :   if (m == 1) return 1;
     475          35 :   if (typ(N) == t_POL && varn(N) == varn(D))
     476          21 :     m = cgcd(m, RgX_deflate_order(N));
     477          35 :   return m;
     478             : }
     479             : /* F a t_RFRAC */
     480             : GEN
     481          98 : rfrac_deflate_max(GEN F, long *m)
     482             : {
     483          98 :   *m = rfrac_deflate_order(F);
     484          98 :   return rfrac_deflate(F, *m);
     485             : }
     486             : /* F a t_RFRAC */
     487             : GEN
     488          98 : rfrac_deflate(GEN F, long m)
     489             : {
     490          98 :   GEN N = gel(F,1), D = gel(F,2);
     491          98 :   if (m == 1) return F;
     492          35 :   if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
     493          35 :   D = RgX_deflate(D, m); return mkrfrac(N, D);
     494             : }
     495             : 
     496             : /* return x0(X^d) */
     497             : GEN
     498      796535 : RgX_inflate(GEN x0, long d)
     499             : {
     500      796535 :   long i, id, dy, dx = degpol(x0);
     501      796535 :   GEN x = x0 + 2, z, y;
     502      796535 :   if (dx <= 0) return leafcopy(x0);
     503      788933 :   dy = dx*d;
     504      788933 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     505      788933 :   z = y + 2;
     506    27826068 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     507    11394799 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     508      788933 :   return y;
     509             : }
     510             : 
     511             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     512             : static GEN
     513     5806729 : RgX_translate_basecase(GEN P, GEN c)
     514             : {
     515     5806729 :   pari_sp av = avma;
     516             :   GEN Q, R;
     517             :   long i, k, n;
     518             : 
     519     5806729 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     520     5804690 :   Q = leafcopy(P);
     521     5804697 :   R = Q+2; n = degpol(P);
     522     5804694 :   if (isint1(c))
     523             :   {
     524       26094 :     for (i=1; i<=n; i++)
     525             :     {
     526      134643 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
     527       22259 :       if (gc_needed(av,2))
     528             :       {
     529           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
     530           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     531             :       }
     532             :     }
     533             :   }
     534     5800860 :   else if (isintm1(c))
     535             :   {
     536      892245 :     for (i=1; i<=n; i++)
     537             :     {
     538     2828676 :       for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
     539      677249 :       if (gc_needed(av,2))
     540             :       {
     541           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
     542           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     543             :       }
     544             :     }
     545             :   }
     546             :   else
     547             :   {
     548    19111606 :     for (i=1; i<=n; i++)
     549             :     {
     550    44911728 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
     551    13525746 :       if (gc_needed(av,2))
     552             :       {
     553           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
     554           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     555             :       }
     556             :     }
     557             :   }
     558     5804348 :   return gerepilecopy(av, Q);
     559             : }
     560             : GEN
     561     5807156 : RgX_translate(GEN P, GEN c)
     562             : {
     563     5807156 :   pari_sp av = avma;
     564     5807156 :   long n = degpol(P);
     565     5807151 :   if (n < 40)
     566     5806731 :     return RgX_translate_basecase(P, c);
     567             :   else
     568             :   {
     569         420 :     long d = n >> 1;
     570         420 :     GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
     571         420 :     GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
     572         420 :     GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
     573         420 :     return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
     574             :   }
     575             : }
     576             : 
     577             : /* P(ax + b) */
     578             : GEN
     579           0 : RgX_affine(GEN P, GEN a, GEN b)
     580             : {
     581           0 :   if (signe(b)) P = RgX_translate(P, b);
     582           0 :   return RgX_unscale(P, a);
     583             : }
     584             : 
     585             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     586             : GEN
     587       23651 : RgXQX_translate(GEN P, GEN c, GEN T)
     588             : {
     589       23651 :   pari_sp av = avma;
     590             :   GEN Q, R;
     591             :   long i, k, n;
     592             : 
     593       23651 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     594       23322 :   Q = leafcopy(P);
     595       23322 :   R = Q+2; n = degpol(P);
     596       82367 :   for (i=1; i<=n; i++)
     597             :   {
     598      260749 :     for (k=n-i; k<n; k++)
     599             :     {
     600      201704 :       pari_sp av2 = avma;
     601      201704 :       gel(R,k) = gerepileupto(av2,
     602      201704 :                    RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
     603             :     }
     604       59045 :     if (gc_needed(av,2))
     605             :     {
     606           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     607           0 :       Q = gerepilecopy(av, Q); R = Q+2;
     608             :     }
     609             :   }
     610       23322 :   return gerepilecopy(av, Q);
     611             : }
     612             : 
     613             : /********************************************************************/
     614             : /**                                                                **/
     615             : /**                          CONVERSIONS                           **/
     616             : /**                       (not memory clean)                       **/
     617             : /**                                                                **/
     618             : /********************************************************************/
     619             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     620             :  * but everything else is */
     621             : static GEN
     622      166649 : QXQ_to_mod(GEN x, GEN T)
     623             : {
     624             :   long d;
     625      166649 :   switch(typ(x))
     626             :   {
     627       62304 :     case t_INT:  return icopy(x);
     628        2079 :     case t_FRAC: return gcopy(x);
     629      102266 :     case t_POL:
     630      102266 :       d = degpol(x);
     631      102266 :       if (d < 0) return gen_0;
     632       99998 :       if (d == 0) return gcopy(gel(x,2));
     633       97859 :       return mkpolmod(RgX_copy(x), T);
     634           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     635             :              return NULL;/* LCOV_EXCL_LINE */
     636             :   }
     637             : }
     638             : /* pure shallow version */
     639             : GEN
     640      709351 : QXQ_to_mod_shallow(GEN x, GEN T)
     641             : {
     642             :   long d;
     643      709351 :   switch(typ(x))
     644             :   {
     645      454645 :     case t_INT:
     646      454645 :     case t_FRAC: return x;
     647      254706 :     case t_POL:
     648      254706 :       d = degpol(x);
     649      254706 :       if (d < 0) return gen_0;
     650      209296 :       if (d == 0) return gel(x,2);
     651      194294 :       return mkpolmod(x, T);
     652           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     653             :              return NULL;/* LCOV_EXCL_LINE */
     654             :   }
     655             : }
     656             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
     657             :  * Not memory clean because T not copied, but everything else is */
     658             : static GEN
     659       37184 : QXQX_to_mod(GEN z, GEN T)
     660             : {
     661       37184 :   long i,l = lg(z);
     662       37184 :   GEN x = cgetg(l,t_POL);
     663      185808 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     664       37184 :   x[1] = z[1]; return normalizepol_lg(x,l);
     665             : }
     666             : /* pure shallow version */
     667             : GEN
     668      161699 : QXQX_to_mod_shallow(GEN z, GEN T)
     669             : {
     670      161699 :   long i,l = lg(z);
     671      161699 :   GEN x = cgetg(l,t_POL);
     672      790270 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
     673      161699 :   x[1] = z[1]; return normalizepol_lg(x,l);
     674             : }
     675             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     676             : GEN
     677       12495 : QXQXV_to_mod(GEN V, GEN T)
     678             : {
     679       12495 :   long i, l = lg(V);
     680       12495 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     681       49679 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     682       12495 :   return z;
     683             : }
     684             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     685             : GEN
     686       18177 : QXQV_to_mod(GEN V, GEN T)
     687             : {
     688       18177 :   long i, l = lg(V);
     689       18177 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     690       36202 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
     691       18177 :   return z;
     692             : }
     693             : 
     694             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     695             : GEN
     696       13412 : QXQC_to_mod_shallow(GEN V, GEN T)
     697             : {
     698       13412 :   long i, l = lg(V);
     699       13412 :   GEN z = cgetg(l, t_COL);
     700       94192 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
     701       13412 :   return z;
     702             : }
     703             : 
     704             : GEN
     705        6251 : QXQM_to_mod_shallow(GEN V, GEN T)
     706             : {
     707        6251 :   long i, l = lg(V);
     708        6251 :   GEN z = cgetg(l, t_MAT);
     709       19663 :   for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
     710        6251 :   return z;
     711             : }
     712             : 
     713             : GEN
     714     6151719 : RgX_renormalize_lg(GEN x, long lx)
     715             : {
     716             :   long i;
     717     9169552 :   for (i = lx-1; i>1; i--)
     718     8720564 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     719     6151719 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     720     6151718 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     721             : }
     722             : 
     723             : GEN
     724     1488756 : RgV_to_RgX(GEN x, long v)
     725             : {
     726     1488756 :   long i, k = lg(x);
     727             :   GEN p;
     728             : 
     729     4137857 :   while (--k && gequal0(gel(x,k)));
     730     1488755 :   if (!k) return pol_0(v);
     731     1481911 :   i = k+2; p = cgetg(i,t_POL);
     732     1481910 :   p[1] = evalsigne(1) | evalvarn(v);
     733    10207160 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     734     1481910 :   return p;
     735             : }
     736             : GEN
     737      189375 : RgV_to_RgX_reverse(GEN x, long v)
     738             : {
     739      189375 :   long j, k, l = lg(x);
     740             :   GEN p;
     741             : 
     742      190852 :   for (k = 1; k < l; k++)
     743      190852 :     if (!gequal0(gel(x,k))) break;
     744      189375 :   if (k == l) return pol_0(v);
     745      189375 :   k -= 1;
     746      189375 :   l -= k;
     747      189375 :   x += k;
     748      189375 :   p = cgetg(l+1,t_POL);
     749      189375 :   p[1] = evalsigne(1) | evalvarn(v);
     750      993263 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     751      189375 :   return p;
     752             : }
     753             : 
     754             : /* return the (N-dimensional) vector of coeffs of p */
     755             : GEN
     756    13110421 : RgX_to_RgC(GEN x, long N)
     757             : {
     758             :   long i, l;
     759             :   GEN z;
     760    13110421 :   l = lg(x)-1; x++;
     761    13110421 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     762    13110421 :   z = cgetg(N+1,t_COL);
     763    79800780 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     764    23986841 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     765    13110546 :   return z;
     766             : }
     767             : GEN
     768     1344064 : Rg_to_RgC(GEN x, long N)
     769             : {
     770     1344064 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     771             : }
     772             : 
     773             : /* vector of polynomials (in v) whose coefs are given by the columns of x */
     774             : GEN
     775      300636 : RgM_to_RgXV(GEN x, long v)
     776     1282353 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
     777             : GEN
     778        7097 : RgM_to_RgXV_reverse(GEN x, long v)
     779       28388 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
     780             : 
     781             : /* matrix whose entries are given by the coeffs of the polynomials in
     782             :  * vector v (considered as degree n-1 polynomials) */
     783             : GEN
     784      333260 : RgV_to_RgM(GEN x, long n)
     785     1674829 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
     786             : 
     787             : GEN
     788       77517 : RgXV_to_RgM(GEN x, long n)
     789      392271 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
     790             : 
     791             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     792             : GEN
     793       23116 : RgM_to_RgXX(GEN x, long v,long w)
     794             : {
     795       23116 :   long j, lx = lg(x);
     796       23116 :   GEN y = cgetg(lx+1, t_POL);
     797       23116 :   y[1] = evalsigne(1) | evalvarn(v);
     798       23116 :   y++;
     799      127478 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     800       23116 :   return normalizepol_lg(--y, lx+1);
     801             : }
     802             : 
     803             : /* matrix whose entries are given by the coeffs of the polynomial v in
     804             :  * two variables (considered as degree n-1 polynomials) */
     805             : GEN
     806         322 : RgXX_to_RgM(GEN v, long n)
     807             : {
     808         322 :   long j, N = lg(v)-1;
     809         322 :   GEN y = cgetg(N, t_MAT);
     810        1043 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     811         322 :   return y;
     812             : }
     813             : 
     814             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     815             : GEN
     816       31244 : RgXY_swapspec(GEN x, long n, long w, long nx)
     817             : {
     818       31244 :   long j, ly = n+3;
     819       31244 :   GEN y = cgetg(ly, t_POL);
     820       31244 :   y[1] = evalsigne(1);
     821      384620 :   for (j=2; j<ly; j++)
     822             :   {
     823             :     long k;
     824      353376 :     GEN a = cgetg(nx+2,t_POL);
     825      353376 :     a[1] = evalsigne(1) | evalvarn(w);
     826     1809209 :     for (k=0; k<nx; k++)
     827             :     {
     828     1455833 :       GEN xk = gel(x,k);
     829     1455833 :       if (typ(xk)==t_POL)
     830     1367317 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     831             :       else
     832       88516 :         gel(a,k+2) = j==2 ? xk: gen_0;
     833             :     }
     834      353376 :     gel(y,j) = normalizepol_lg(a, nx+2);
     835             :   }
     836       31244 :   return normalizepol_lg(y,ly);
     837             : }
     838             : 
     839             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     840             : GEN
     841         952 : RgXY_swap(GEN x, long n, long w)
     842             : {
     843         952 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     844         952 :   setvarn(z, varn(x)); return z;
     845             : }
     846             : 
     847             : long
     848         288 : RgXY_degreex(GEN b)
     849             : {
     850         288 :   long deg = 0, i;
     851         288 :   if (!signe(b)) return -1;
     852        1396 :   for (i = 2; i < lg(b); ++i)
     853             :   {
     854        1108 :     GEN bi = gel(b, i);
     855        1108 :     if (typ(bi) == t_POL)
     856        1030 :       deg = maxss(deg, degpol(bi));
     857             :   }
     858         288 :   return deg;
     859             : }
     860             : 
     861             : GEN
     862       32578 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
     863             : 
     864             : /* return (x % X^n). Shallow */
     865             : GEN
     866     7905450 : RgXn_red_shallow(GEN a, long n)
     867             : {
     868     7905450 :   long i, L = n+2, l = lg(a);
     869             :   GEN  b;
     870     7905450 :   if (L >= l) return a; /* deg(x) < n */
     871     5720040 :   b = cgetg(L, t_POL); b[1] = a[1];
     872    36560167 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     873     5720042 :   return normalizepol_lg(b,L);
     874             : }
     875             : 
     876             : GEN
     877         483 : RgXnV_red_shallow(GEN x, long n)
     878        2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
     879             : 
     880             : /* return (x * X^n). Shallow */
     881             : GEN
     882   164731989 : RgX_shift_shallow(GEN a, long n)
     883             : {
     884   164731989 :   long i, l = lg(a);
     885             :   GEN  b;
     886   164731989 :   if (l == 2 || !n) return a;
     887   107873449 :   l += n;
     888   107873449 :   if (n < 0)
     889             :   {
     890    55066577 :     if (l <= 2) return pol_0(varn(a));
     891    53601602 :     b = cgetg(l, t_POL); b[1] = a[1];
     892    53602342 :     a -= n;
     893   163698210 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     894             :   } else {
     895    52806872 :     b = cgetg(l, t_POL); b[1] = a[1];
     896    52811303 :     a -= n; n += 2;
     897   114967577 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     898   204266895 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     899             :   }
     900   106413645 :   return b;
     901             : }
     902             : /* return (x * X^n). */
     903             : GEN
     904     1567501 : RgX_shift(GEN a, long n)
     905             : {
     906     1567501 :   long i, l = lg(a);
     907             :   GEN  b;
     908     1567501 :   if (l == 2 || !n) return RgX_copy(a);
     909     1567137 :   l += n;
     910     1567137 :   if (n < 0)
     911             :   {
     912         847 :     if (l <= 2) return pol_0(varn(a));
     913         777 :     b = cgetg(l, t_POL); b[1] = a[1];
     914         777 :     a -= n;
     915        2947 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     916             :   } else {
     917     1566290 :     b = cgetg(l, t_POL); b[1] = a[1];
     918     1566290 :     a -= n; n += 2;
     919     3915502 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     920     4318398 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     921             :   }
     922     1567067 :   return b;
     923             : }
     924             : 
     925             : GEN
     926      317037 : RgX_rotate_shallow(GEN P, long k, long p)
     927             : {
     928      317037 :   long i, l = lgpol(P);
     929             :   GEN r;
     930      317037 :   if (signe(P)==0)
     931        1365 :     return pol_0(varn(P));
     932      315672 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     933     2100644 :   for(i=0; i<p; i++)
     934             :   {
     935     1784972 :     long s = 2+(i+k)%p;
     936     1784972 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
     937             :   }
     938      315672 :   return RgX_renormalize(r);
     939             : }
     940             : 
     941             : GEN
     942     2984861 : RgX_mulXn(GEN x, long d)
     943             : {
     944             :   pari_sp av;
     945             :   GEN z;
     946             :   long v;
     947     2984861 :   if (d >= 0) return RgX_shift(x, d);
     948     1472244 :   d = -d;
     949     1472244 :   v = RgX_val(x);
     950     1472244 :   if (v >= d) return RgX_shift(x, -d);
     951     1472237 :   av = avma;
     952     1472237 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
     953     1472237 :   return gerepileupto(av, z);
     954             : }
     955             : 
     956             : long
     957         588 : RgXV_maxdegree(GEN x)
     958             : {
     959         588 :   long d = -1, i, l = lg(x);
     960        4494 :   for (i = 1; i < l; i++)
     961        3906 :     d = maxss(d, degpol(gel(x,i)));
     962         588 :   return d;
     963             : }
     964             : 
     965             : long
     966     3633368 : RgX_val(GEN x)
     967             : {
     968     3633368 :   long i, lx = lg(x);
     969     3633368 :   if (lx == 2) return LONG_MAX;
     970     4562932 :   for (i = 2; i < lx; i++)
     971     4562876 :     if (!isexactzero(gel(x,i))) break;
     972     3633214 :   if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
     973     3633158 :   return i - 2;
     974             : }
     975             : long
     976    80182599 : RgX_valrem(GEN x, GEN *Z)
     977             : {
     978    80182599 :   long v, i, lx = lg(x);
     979    80182599 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
     980   125521769 :   for (i = 2; i < lx; i++)
     981   125522764 :     if (!isexactzero(gel(x,i))) break;
     982             :   /* possible with nonrational zeros */
     983    80182489 :   if (i == lx)
     984             :   {
     985          14 :     *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
     986          14 :     return LONG_MAX;
     987             :   }
     988    80182475 :   v = i - 2;
     989    80182475 :   *Z = RgX_shift_shallow(x, -v);
     990    80193394 :   return v;
     991             : }
     992             : long
     993      831791 : RgX_valrem_inexact(GEN x, GEN *Z)
     994             : {
     995             :   long v;
     996      831791 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
     997      847559 :   for (v = 0;; v++)
     998      847559 :     if (!gequal0(gel(x,2+v))) break;
     999      831784 :   if (Z) *Z = RgX_shift_shallow(x, -v);
    1000      831784 :   return v;
    1001             : }
    1002             : 
    1003             : GEN
    1004       63875 : RgXQC_red(GEN x, GEN T)
    1005      394877 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
    1006             : 
    1007             : GEN
    1008        1204 : RgXQV_red(GEN x, GEN T)
    1009       27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
    1010             : 
    1011             : GEN
    1012       12257 : RgXQM_red(GEN x, GEN T)
    1013       76132 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
    1014             : 
    1015             : GEN
    1016         322 : RgXQM_mul(GEN P, GEN Q, GEN T)
    1017             : {
    1018         322 :   return RgXQM_red(RgM_mul(P, Q), T);
    1019             : }
    1020             : 
    1021             : GEN
    1022      456154 : RgXQX_red(GEN P, GEN T)
    1023             : {
    1024      456154 :   long i, l = lg(P);
    1025      456154 :   GEN Q = cgetg(l, t_POL);
    1026      456154 :   Q[1] = P[1];
    1027     2383598 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
    1028      456154 :   return normalizepol_lg(Q, l);
    1029             : }
    1030             : 
    1031             : GEN
    1032      793097 : RgX_deriv(GEN x)
    1033             : {
    1034      793097 :   long i,lx = lg(x)-1;
    1035             :   GEN y;
    1036             : 
    1037      793097 :   if (lx<3) return pol_0(varn(x));
    1038      790402 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
    1039     3517094 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
    1040      790393 :   y[1] = x[1]; return normalizepol_lg(y,i);
    1041             : }
    1042             : 
    1043             : GEN
    1044     1377452 : RgX_recipspec_shallow(GEN x, long l, long n)
    1045             : {
    1046             :   long i;
    1047     1377452 :   GEN z = cgetg(n+2,t_POL);
    1048     1377464 :   z[1] = 0; z += 2;
    1049    36053375 :   for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
    1050     1600568 :   for(   ; i<n; i++) gel(z, n-i-1) = gen_0;
    1051     1377464 :   return normalizepol_lg(z-2,n+2);
    1052             : }
    1053             : 
    1054             : GEN
    1055      637320 : RgXn_recip_shallow(GEN P, long n)
    1056             : {
    1057      637320 :   GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
    1058      637331 :   setvarn(Q, varn(P));
    1059      637331 :   return Q;
    1060             : }
    1061             : 
    1062             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1063             : GEN
    1064       31150 : RgX_recip(GEN x)
    1065             : {
    1066             :   long lx, i, j;
    1067       31150 :   GEN y = cgetg_copy(x, &lx);
    1068      272104 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
    1069       31150 :   return normalizepol_lg(y,lx);
    1070             : }
    1071             : /* shallow version */
    1072             : GEN
    1073       57526 : RgX_recip_shallow(GEN x)
    1074             : {
    1075             :   long lx, i, j;
    1076       57526 :   GEN y = cgetg_copy(x, &lx);
    1077      344995 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1078       57526 :   return normalizepol_lg(y,lx);
    1079             : }
    1080             : 
    1081             : GEN
    1082     3629079 : RgX_recip_i(GEN x)
    1083             : {
    1084             :   long lx, i, j;
    1085     3629079 :   GEN y = cgetg_copy(x, &lx);
    1086    19216310 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1087     3629086 :   return y;
    1088             : }
    1089             : /*******************************************************************/
    1090             : /*                                                                 */
    1091             : /*                      ADDITION / SUBTRACTION                     */
    1092             : /*                                                                 */
    1093             : /*******************************************************************/
    1094             : /* same variable */
    1095             : GEN
    1096    76929232 : RgX_add(GEN x, GEN y)
    1097             : {
    1098    76929232 :   long i, lx = lg(x), ly = lg(y);
    1099             :   GEN z;
    1100    76929232 :   if (ly <= lx) {
    1101    66737495 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1102   264536980 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1103   118618385 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1104    66713636 :     z = normalizepol_lg(z, lx);
    1105             :   } else {
    1106    10191737 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1107    43486359 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1108    31404943 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1109    10193118 :     z = normalizepol_lg(z, ly);
    1110             :   }
    1111    76920783 :   return z;
    1112             : }
    1113             : GEN
    1114    61843517 : RgX_sub(GEN x, GEN y)
    1115             : {
    1116    61843517 :   long i, lx = lg(x), ly = lg(y);
    1117             :   GEN z;
    1118    61843517 :   if (ly <= lx) {
    1119    29633238 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1120   155366441 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1121    55217438 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1122    29609879 :     z = normalizepol_lg(z, lx);
    1123             :   } else {
    1124    32210279 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1125   120268473 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1126    71563503 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1127    32176359 :     z = normalizepol_lg(z, ly);
    1128             :   }
    1129    61838255 :   return z;
    1130             : }
    1131             : GEN
    1132     5172674 : RgX_neg(GEN x)
    1133             : {
    1134     5172674 :   long i, lx = lg(x);
    1135     5172674 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    1136    40227381 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    1137     5172594 :   return y;
    1138             : }
    1139             : 
    1140             : GEN
    1141    18860129 : RgX_Rg_add(GEN y, GEN x)
    1142             : {
    1143             :   GEN z;
    1144    18860129 :   long lz = lg(y), i;
    1145    18860129 :   if (lz == 2) return scalarpol(x,varn(y));
    1146    15822099 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1147    15822046 :   gel(z,2) = gadd(gel(y,2),x);
    1148    61628824 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1149             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1150             :    * an inexact 0 */
    1151    15822035 :   return normalizepol_lg(z,lz);
    1152             : }
    1153             : GEN
    1154       64697 : RgX_Rg_add_shallow(GEN y, GEN x)
    1155             : {
    1156             :   GEN z;
    1157       64697 :   long lz = lg(y), i;
    1158       64697 :   if (lz == 2) return scalarpol(x,varn(y));
    1159       64697 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1160       64697 :   gel(z,2) = gadd(gel(y,2),x);
    1161      129544 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1162       64697 :   return normalizepol_lg(z,lz);
    1163             : }
    1164             : GEN
    1165      182848 : RgX_Rg_sub(GEN y, GEN x)
    1166             : {
    1167             :   GEN z;
    1168      182848 :   long lz = lg(y), i;
    1169      182848 :   if (lz == 2)
    1170             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1171        2548 :     long v = varn(y);
    1172        2548 :     if (isrationalzero(x)) return pol_0(v);
    1173          14 :     z = cgetg(3,t_POL);
    1174          14 :     z[1] = gequal0(x)? evalvarn(v)
    1175          14 :                    : evalvarn(v) | evalsigne(1);
    1176          14 :     gel(z,2) = gneg(x); return z;
    1177             :   }
    1178      180300 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1179      180299 :   gel(z,2) = gsub(gel(y,2),x);
    1180      470638 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1181      180300 :   return normalizepol_lg(z,lz);
    1182             : }
    1183             : GEN
    1184     4684885 : Rg_RgX_sub(GEN x, GEN y)
    1185             : {
    1186             :   GEN z;
    1187     4684885 :   long lz = lg(y), i;
    1188     4684885 :   if (lz == 2) return scalarpol(x,varn(y));
    1189     4683807 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1190     4683788 :   gel(z,2) = gsub(x, gel(y,2));
    1191     7076637 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1192     4683791 :   return normalizepol_lg(z,lz);
    1193             : }
    1194             : /*******************************************************************/
    1195             : /*                                                                 */
    1196             : /*                  KARATSUBA MULTIPLICATION                       */
    1197             : /*                                                                 */
    1198             : /*******************************************************************/
    1199             : #if 0
    1200             : /* to debug Karatsuba-like routines */
    1201             : GEN
    1202             : zx_debug_spec(GEN x, long nx)
    1203             : {
    1204             :   GEN z = cgetg(nx+2,t_POL);
    1205             :   long i;
    1206             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1207             :   z[1] = evalsigne(1); return z;
    1208             : }
    1209             : 
    1210             : GEN
    1211             : RgX_debug_spec(GEN x, long nx)
    1212             : {
    1213             :   GEN z = cgetg(nx+2,t_POL);
    1214             :   long i;
    1215             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1216             :   z[1] = evalsigne(1); return z;
    1217             : }
    1218             : #endif
    1219             : 
    1220             : /* generic multiplication */
    1221             : GEN
    1222     8506518 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1223             : {
    1224             :   GEN z, t;
    1225             :   long i;
    1226     8506518 :   if (nx == ny) {
    1227     1326905 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1228     4240020 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1229     1326897 :     return normalizepol_lg(z, nx+2);
    1230             :   }
    1231     7179613 :   if (ny < nx) {
    1232     6993058 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1233    25833582 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1234    16787185 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1235     6992455 :     return normalizepol_lg(z, nx+2);
    1236             :   } else {
    1237      186555 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1238     3674649 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1239      451027 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1240      186560 :     return normalizepol_lg(z, ny+2);
    1241             :   }
    1242             : }
    1243             : GEN
    1244      222775 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1245             : {
    1246             :   GEN z, t;
    1247             :   long i;
    1248      222775 :   if (nx == ny) {
    1249       12810 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1250     2183524 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1251       12810 :     return normalizepol_lg(z, nx+2);
    1252             :   }
    1253      209965 :   if (ny < nx) {
    1254      208187 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1255     3723510 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1256     2377417 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1257      208187 :     return normalizepol_lg(z, nx+2);
    1258             :   } else {
    1259        1778 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1260      330904 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1261       12222 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1262        1778 :     return normalizepol_lg(z, ny+2);
    1263             :   }
    1264             : }
    1265             : 
    1266             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1267             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1268             :  * t_VECSMALL, to hold this, which can be left on stack: gerepile
    1269             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1270             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1271             : static GEN
    1272    55650894 : RgXspec_kill0(GEN x, long lx)
    1273             : {
    1274    55650894 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
    1275             :   long i;
    1276   189612059 :   for (i=0; i <lx; i++)
    1277             :   {
    1278   133961168 :     GEN c = gel(x,i);
    1279   133961168 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1280             :   }
    1281    55650891 :   return z;
    1282             : }
    1283             : 
    1284             : INLINE GEN
    1285    93180251 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1286             : {
    1287    93180251 :   pari_sp av = avma;
    1288    93180251 :   GEN s = NULL;
    1289             :   long i;
    1290             : 
    1291   354955903 :   for (i=a; i<b; i++)
    1292   261785428 :     if (gel(y,i) && gel(x,-i))
    1293             :     {
    1294   183752044 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1295   183746537 :       s = s? gadd(s, t): t;
    1296             :     }
    1297    93170475 :   return s? gerepileupto(av, s): gen_0;
    1298             : }
    1299             : 
    1300             : /* assume nx >= ny > 0, return x * y * t^v */
    1301             : static GEN
    1302    21219774 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1303             : {
    1304             :   long i, lz, nz;
    1305             :   GEN z;
    1306             : 
    1307    21219774 :   x = RgXspec_kill0(x,nx);
    1308    21219747 :   y = RgXspec_kill0(y,ny);
    1309    21219752 :   lz = nx + ny + 1; nz = lz-2;
    1310    21219752 :   lz += v;
    1311    21219752 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1312    32953129 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1313    56362011 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1314    40467546 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1315    35142434 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1316    21219247 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1317             : }
    1318             : 
    1319             : /* return (x * X^d) + y. Assume d > 0 */
    1320             : GEN
    1321     9067675 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1322             : {
    1323             :   GEN x, y, xd, yd, zd;
    1324             :   long a, lz, nx, ny;
    1325             : 
    1326     9067675 :   if (!signe(x0)) return y0;
    1327     8962080 :   ny = lgpol(y0);
    1328     8962091 :   nx = lgpol(x0);
    1329     8962187 :   zd = (GEN)avma;
    1330     8962187 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1331     8962187 :   if (a <= 0)
    1332             :   {
    1333     1412810 :     lz = nx+d+2;
    1334     1412810 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1335     3264491 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1336     1412810 :     x = zd + a;
    1337     1430402 :     while (zd > x) gel(--zd,0) = gen_0;
    1338             :   }
    1339             :   else
    1340             :   {
    1341     7549377 :     xd = new_chunk(d); yd = y+d;
    1342     7549379 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1343     7549287 :     lz = (a>nx)? ny+2: lg(x)+d;
    1344    38097265 :     x += 2; while (xd > x) *--zd = *--xd;
    1345             :   }
    1346    21985264 :   while (yd > y) *--zd = *--yd;
    1347     8962097 :   *--zd = x0[1];
    1348     8962097 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1349             : }
    1350             : GEN
    1351      515334 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1352             : {
    1353             :   GEN x, y, xd, yd, zd;
    1354             :   long a, lz, nx, ny;
    1355             : 
    1356      515334 :   if (!signe(x0)) return RgX_copy(y0);
    1357      514564 :   nx = lgpol(x0);
    1358      514564 :   ny = lgpol(y0);
    1359      514564 :   zd = (GEN)avma;
    1360      514564 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1361      514564 :   if (a <= 0)
    1362             :   {
    1363      291789 :     lz = nx+d+2;
    1364      291789 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1365     4294362 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1366      291789 :     x = zd + a;
    1367      757955 :     while (zd > x) gel(--zd,0) = gen_0;
    1368             :   }
    1369             :   else
    1370             :   {
    1371      222775 :     xd = new_chunk(d); yd = y+d;
    1372      222775 :     x = RgX_addspec(x,yd, nx,a);
    1373      222775 :     lz = (a>nx)? ny+2: lg(x)+d;
    1374     8417612 :     x += 2; while (xd > x) *--zd = *--xd;
    1375             :   }
    1376     2619571 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1377      514564 :   *--zd = x0[1];
    1378      514564 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1379             : }
    1380             : 
    1381             : /* return x * y mod t^n */
    1382             : static GEN
    1383     6601346 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1384             : {
    1385     6601346 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1386             :   GEN z;
    1387     6601346 :   if (lx < 0) return pol_0(varn(x));
    1388     6601346 :   if (ly < 0) return pol_0(varn(x));
    1389     6601346 :   z = cgetg(lz, t_POL) + 2;
    1390     6601346 :   x+=2; if (lx > n) lx = n;
    1391     6601346 :   y+=2; if (ly > n) ly = n;
    1392     6601346 :   z[-1] = x[-1];
    1393     6601346 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1394     6601346 :   x = RgXspec_kill0(x, lx);
    1395     6601346 :   y = RgXspec_kill0(y, ly);
    1396             :   /* x:y:z [i] = term of degree i */
    1397    26177251 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1398    11846647 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1399     6647372 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1400     6601346 :   return normalizepol_lg(z - 2, lz);
    1401             : }
    1402             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1403             : static GEN
    1404     9375638 : RgXn_mul2(GEN f, GEN g, long n)
    1405             : {
    1406     9375638 :   pari_sp av = avma;
    1407             :   GEN fe,fo, ge,go, l,h,m;
    1408             :   long n0, n1;
    1409     9375638 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1410     6636766 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1411       35420 :   n0 = n>>1; n1 = n-n0;
    1412       35420 :   RgX_even_odd(f, &fe, &fo);
    1413       35420 :   RgX_even_odd(g, &ge, &go);
    1414       35420 :   l = RgXn_mul2(fe,ge,n1);
    1415       35420 :   h = RgXn_mul2(fo,go,n0);
    1416       35420 :   m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1417             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1418             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1419       35420 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1420             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1421       35420 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1422       35420 :   m = RgX_inflate(m,2);
    1423             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1424       35420 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1425       35420 :   h = RgX_inflate(h,2);
    1426       35420 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1427       35420 :   return gerepileupto(av, h);
    1428             : }
    1429             : /* (f*g) \/ x^n */
    1430             : static GEN
    1431     1583632 : RgX_mulhigh_i2(GEN f, GEN g, long n)
    1432             : {
    1433     1583632 :   long d = degpol(f)+degpol(g) + 1 - n;
    1434             :   GEN h;
    1435     1583632 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1436       31376 :   h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
    1437             :                                   RgX_recip_i(g), d));
    1438       31376 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1439             : }
    1440             : 
    1441             : /* (f*g) \/ x^n */
    1442             : static GEN
    1443           0 : RgX_sqrhigh_i2(GEN f, long n)
    1444             : {
    1445           0 :   long d = 2*degpol(f)+ 1 - n;
    1446             :   GEN h;
    1447           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1448           0 :   h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
    1449           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1450             : }
    1451             : 
    1452             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1453             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1454             :  * Only c, c0, c1, c2 are genuine GEN.
    1455             :  */
    1456             : GEN
    1457    21947811 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1458             : {
    1459             :   GEN a0, c, c0;
    1460    21947811 :   long n0, n0a, i, v = 0;
    1461             :   pari_sp av;
    1462             : 
    1463    28331352 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1464    27831783 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1465    21947796 :   if (na < nb) swapspec(a,b, na,nb);
    1466    21947796 :   if (!nb) return pol_0(0);
    1467             : 
    1468    21699607 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1469      479858 :   RgX_shift_inplace_init(v);
    1470      479875 :   i = (na>>1); n0 = na-i; na = i;
    1471      479875 :   av = avma; a0 = a+n0; n0a = n0;
    1472     1335391 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1473             : 
    1474      479875 :   if (nb > n0)
    1475             :   {
    1476             :     GEN b0,c1,c2;
    1477             :     long n0b;
    1478             : 
    1479      478566 :     nb -= n0; b0 = b+n0; n0b = n0;
    1480     1446545 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1481      478566 :     c = RgX_mulspec(a,b,n0a,n0b);
    1482      478566 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1483             : 
    1484      478566 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1485      478566 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1486             : 
    1487      478566 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1488      478566 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1489      478566 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1490             :   }
    1491             :   else
    1492             :   {
    1493        1309 :     c = RgX_mulspec(a,b,n0a,nb);
    1494        1309 :     c0 = RgX_mulspec(a0,b,na,nb);
    1495             :   }
    1496      479875 :   c0 = RgX_addmulXn(c0,c,n0);
    1497      479875 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1498             : }
    1499             : 
    1500             : INLINE GEN
    1501       50354 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1502             : {
    1503       50354 :   pari_sp av = avma;
    1504       50354 :   GEN s = NULL;
    1505       50354 :   long j, l = (i+1)>>1;
    1506      139496 :   for (j=a; j<l; j++)
    1507             :   {
    1508       89142 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1509       89142 :     if (xj && xx)
    1510             :     {
    1511       81498 :       GEN t = gmul(xj, xx);
    1512       81498 :       s = s? gadd(s, t): t;
    1513             :     }
    1514             :   }
    1515       50354 :   if (s) s = gshift(s,1);
    1516       50354 :   if ((i&1) == 0)
    1517             :   {
    1518       29603 :     GEN t = gel(x, i>>1);
    1519       29603 :     if (t) {
    1520       27209 :       t = gsqr(t);
    1521       27209 :       s = s? gadd(s, t): t;
    1522             :     }
    1523             :   }
    1524       50354 :   return s? gerepileupto(av,s): gen_0;
    1525             : }
    1526             : static GEN
    1527        8852 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1528             : {
    1529             :   long i, lz, nz;
    1530             :   GEN z;
    1531             : 
    1532        8852 :   if (!nx) return pol_0(0);
    1533        8852 :   x = RgXspec_kill0(x,nx);
    1534        8852 :   lz = (nx << 1) + 1, nz = lz-2;
    1535        8852 :   lz += v;
    1536        8852 :   z = cgetg(lz,t_POL) + 2;
    1537       13024 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1538       38455 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1539       29603 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1540        8852 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1541             : }
    1542             : /* return x^2 mod t^n */
    1543             : static GEN
    1544           0 : RgXn_sqr_basecase(GEN x, long n)
    1545             : {
    1546           0 :   long i, lz = n+2, lx = lgpol(x);
    1547             :   GEN z;
    1548           0 :   if (lx < 0) return pol_0(varn(x));
    1549           0 :   z = cgetg(lz, t_POL);
    1550           0 :   z[1] = x[1];
    1551           0 :   x+=2; if (lx > n) lx = n;
    1552           0 :   x = RgXspec_kill0(x,lx);
    1553           0 :   z+=2;/* x:z [i] = term of degree i */
    1554           0 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1555           0 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1556           0 :   z -= 2; return normalizepol_lg(z, lz);
    1557             : }
    1558             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1559             : static GEN
    1560         315 : RgXn_sqr2(GEN f, long n)
    1561             : {
    1562         315 :   pari_sp av = avma;
    1563             :   GEN fe,fo, l,h,m;
    1564             :   long n0, n1;
    1565         315 :   if (2*degpol(f) < n) return RgX_sqr_i(f);
    1566           0 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1567           0 :   n0 = n>>1; n1 = n-n0;
    1568           0 :   RgX_even_odd(f, &fe, &fo);
    1569           0 :   l = RgXn_sqr(fe,n1);
    1570           0 :   h = RgXn_sqr(fo,n0);
    1571           0 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1572             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1573             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1574           0 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1575             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1576           0 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1577           0 :   m = RgX_inflate(m,2);
    1578             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1579           0 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1580           0 :   h = RgX_inflate(h,2);
    1581           0 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1582           0 :   return gerepileupto(av, h);
    1583             : }
    1584             : GEN
    1585        8891 : RgX_sqrspec(GEN a, long na)
    1586             : {
    1587             :   GEN a0, c, c0, c1;
    1588        8891 :   long n0, n0a, i, v = 0;
    1589             :   pari_sp av;
    1590             : 
    1591       10977 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1592        8891 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1593          39 :   RgX_shift_inplace_init(v);
    1594          39 :   i = (na>>1); n0 = na-i; na = i;
    1595          39 :   av = avma; a0 = a+n0; n0a = n0;
    1596          39 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1597             : 
    1598          39 :   c = RgX_sqrspec(a,n0a);
    1599          39 :   c0 = RgX_sqrspec(a0,na);
    1600          39 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1601          39 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1602          39 :   c0 = RgX_addmulXn(c0,c,n0);
    1603          39 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1604             : }
    1605             : 
    1606             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1607             : GEN
    1608     1692206 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1609             : {
    1610     1692206 :   GEN z = RgX_mul(A, B);
    1611     1692202 :   if (a < b)
    1612       10419 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1613     1681783 :   else if (a > b)
    1614     1100371 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1615             :   else
    1616      581412 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1617     1692201 :   return z;
    1618             : }
    1619             : 
    1620             : GEN
    1621    20509517 : RgX_mul_i(GEN x, GEN y)
    1622             : {
    1623    20509517 :   GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
    1624    20509174 :   setvarn(z, varn(x)); return z;
    1625             : }
    1626             : 
    1627             : GEN
    1628        8813 : RgX_sqr_i(GEN x)
    1629             : {
    1630        8813 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1631        8813 :   setvarn(z,varn(x)); return z;
    1632             : }
    1633             : 
    1634             : /*******************************************************************/
    1635             : /*                                                                 */
    1636             : /*                               DIVISION                          */
    1637             : /*                                                                 */
    1638             : /*******************************************************************/
    1639             : GEN
    1640     3948138 : RgX_Rg_divexact(GEN x, GEN y) {
    1641     3948138 :   long i, lx = lg(x);
    1642             :   GEN z;
    1643     3948138 :   if (lx == 2) return gcopy(x);
    1644     3898318 :   switch(typ(y))
    1645             :   {
    1646     3868596 :     case t_INT:
    1647     3868596 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1648     3694790 :       break;
    1649        4655 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1650             :   }
    1651     3719857 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1652    23802926 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1653     3719000 :   return z;
    1654             : }
    1655             : GEN
    1656    29210116 : RgX_Rg_div(GEN x, GEN y) {
    1657    29210116 :   long i, lx = lg(x);
    1658             :   GEN z;
    1659    29210116 :   if (lx == 2) return gcopy(x);
    1660    28950072 :   switch(typ(y))
    1661             :   {
    1662    19760268 :     case t_INT:
    1663    19760268 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1664     3724627 :       break;
    1665       12929 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1666             :   }
    1667    12901502 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1668    47811741 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1669    12901193 :   return normalizepol_lg(z, lx);
    1670             : }
    1671             : GEN
    1672       18998 : RgX_normalize(GEN x)
    1673             : {
    1674       18998 :   GEN z, d = NULL;
    1675       18998 :   long i, n = lg(x)-1;
    1676       18998 :   for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
    1677       18998 :   if (i == 1) return pol_0(varn(x));
    1678       18998 :   if (i == n && isint1(d)) return x;
    1679        7322 :   n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
    1680       16177 :   for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
    1681        7322 :   gel(z,n) = Rg_get_1(d); return z;
    1682             : }
    1683             : GEN
    1684       10423 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
    1685             : GEN
    1686      242797 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1687             : {
    1688      242797 :   long l = lg(a), i;
    1689             :   GEN a0, z0, z;
    1690             : 
    1691      242797 :   if (l <= 3)
    1692             :   {
    1693           0 :     if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
    1694           0 :     return pol_0(varn(a));
    1695             :   }
    1696      242797 :   z = cgetg(l-1, t_POL);
    1697      242800 :   z[1] = a[1];
    1698      242800 :   a0 = a + l-1;
    1699      242800 :   z0 = z + l-2; *z0 = *a0--;
    1700     2928187 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1701             :   {
    1702     2685388 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1703     2685387 :     gel(z0,0) = t;
    1704             :   }
    1705      242799 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1706      242799 :   return z;
    1707             : }
    1708             : /* Polynomial division x / y:
    1709             :  *   if pr = ONLY_REM return remainder, otherwise return quotient
    1710             :  *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
    1711             :  *   if pr != NULL set *pr to remainder, as the last object on stack */
    1712             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1713             : static GEN
    1714    22903764 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
    1715             : {
    1716             :   pari_sp avy, av, av1;
    1717             :   long dx,dy,dz,i,j,sx,lr;
    1718             :   GEN z,p1,p2,rem,y_lead,mod,p;
    1719             :   GEN (*f)(GEN,GEN);
    1720             : 
    1721    22903764 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1722             : 
    1723    22903764 :   dy = degpol(y);
    1724    22903728 :   y_lead = gel(y,dy+2);
    1725    22903728 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1726             :   {
    1727           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1728           0 :     for (dy--; dy>=0; dy--)
    1729             :     {
    1730           0 :       y_lead = gel(y,dy+2);
    1731           0 :       if (!gequal0(y_lead)) break;
    1732             :     }
    1733             :   }
    1734    22903729 :   if (!dy) /* y is constant */
    1735             :   {
    1736        5373 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1737        5366 :     z = RgX_Rg_div(x, y_lead);
    1738        5366 :     if (pr == ONLY_DIVIDES) return z;
    1739        2580 :     if (pr) *pr = pol_0(varn(x));
    1740        2580 :     return z;
    1741             :   }
    1742    22898356 :   dx = degpol(x);
    1743    22898296 :   if (dx < dy)
    1744             :   {
    1745     3703365 :     if (pr == ONLY_REM) return RgX_copy(x);
    1746      374645 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1747      374624 :     z = pol_0(varn(x));
    1748      374624 :     if (pr) *pr = RgX_copy(x);
    1749      374624 :     return z;
    1750             :   }
    1751             : 
    1752             :   /* x,y in R[X], y non constant */
    1753    19194931 :   av = avma;
    1754    19194931 :   p = NULL;
    1755    19194931 :   if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1756             :   {
    1757      238693 :     z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
    1758      238693 :     if (!z) return gc_NULL(av);
    1759      238693 :     z = FpX_to_mod(z, p);
    1760      238693 :     if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
    1761      122458 :       return gerepileupto(av, z);
    1762      116235 :     *pr = FpX_to_mod(*pr, p);
    1763      116235 :     return gc_all(av, 2, &z, pr);
    1764             :   }
    1765    18956382 :   switch(typ(y_lead))
    1766             :   {
    1767           0 :     case t_REAL:
    1768           0 :       y_lead = ginv(y_lead);
    1769           0 :       f = gmul; mod = NULL;
    1770           0 :       break;
    1771        2137 :     case t_INTMOD:
    1772        2137 :     case t_POLMOD: y_lead = ginv(y_lead);
    1773        2137 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1774        2137 :       break;
    1775    18954245 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1776    18954245 :       f = gdiv; mod = NULL;
    1777             :   }
    1778             : 
    1779    18956382 :   if (y_lead == NULL)
    1780    17131839 :     p2 = gel(x,dx+2);
    1781             :   else {
    1782             :     for(;;) {
    1783     1824543 :       p2 = f(gel(x,dx+2),y_lead);
    1784     1824542 :       p2 = simplify_shallow(p2);
    1785     1824542 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1786             :     }
    1787     1824542 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1788             :     {
    1789           0 :       if (pr == ONLY_DIVIDES) {
    1790           0 :         set_avma(av);
    1791           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1792             :       }
    1793           0 :       if (pr == ONLY_REM)
    1794             :       {
    1795           0 :         if (dx < 0)
    1796           0 :           return gerepilecopy(av, scalarpol(p2, varn(x)));
    1797             :         else
    1798             :         {
    1799             :           GEN t;
    1800           0 :           set_avma(av);
    1801           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1802           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1803           0 :           return t;
    1804             :         }
    1805             :       }
    1806           0 :       if (pr) /* cf ONLY_REM above */
    1807             :       {
    1808           0 :         if (dx < 0)
    1809             :         {
    1810           0 :           p2 = gclone(p2);
    1811           0 :           set_avma(av);
    1812           0 :           z = pol_0(varn(x));
    1813           0 :           x = scalarpol(p2, varn(x));
    1814           0 :           gunclone(p2);
    1815             :         }
    1816             :         else
    1817             :         {
    1818             :           GEN t;
    1819           0 :           set_avma(av);
    1820           0 :           z = pol_0(varn(x));
    1821           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1822           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1823           0 :           x = t;
    1824             :         }
    1825           0 :         *pr = x;
    1826             :       }
    1827             :       else
    1828             :       {
    1829           0 :         set_avma(av);
    1830           0 :         z = pol_0(varn(x));
    1831             :       }
    1832           0 :       return z;
    1833             :     }
    1834             :   }
    1835             :   /* dx >= dy */
    1836    18956381 :   avy = avma;
    1837    18956381 :   dz = dx-dy;
    1838    18956381 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1839    18956360 :   x += 2;
    1840    18956360 :   z += 2;
    1841    18956360 :   y += 2;
    1842    18956360 :   gel(z,dz) = gcopy(p2);
    1843             : 
    1844    51846191 :   for (i=dx-1; i>=dy; i--)
    1845             :   {
    1846    32890029 :     av1=avma; p1=gel(x,i);
    1847  1135114355 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1848    32885022 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1849             : 
    1850    32885023 :     if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
    1851             :     else
    1852    22917015 :       p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
    1853    32889358 :     gel(z,i-dy) = p1;
    1854             :   }
    1855    18956162 :   if (!pr) return gerepileupto(av,z-2);
    1856             : 
    1857    12167362 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1858    12529001 :   for (sx=0; ; i--)
    1859             :   {
    1860    12529001 :     p1 = gel(x,i);
    1861             :     /* we always enter this loop at least once */
    1862    31527830 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1863    12527739 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1864    12527739 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
    1865     3573938 :     if (!isexactzero(p1)) break;
    1866     3429144 :     if (!i) break;
    1867      361580 :     set_avma(av1);
    1868             :   }
    1869    12166605 :   if (pr == ONLY_DIVIDES)
    1870             :   {
    1871        2709 :     if (sx) return gc_NULL(av);
    1872        2688 :     set_avma((pari_sp)rem); return gerepileupto(av,z-2);
    1873             :   }
    1874    12163896 :   lr=i+3; rem -= lr;
    1875    12163896 :   if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
    1876    12120277 :   else p1 = gerepileupto((pari_sp)rem,p1);
    1877    12164631 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    1878    12164631 :   rem[1] = z[-1];
    1879    12164631 :   rem += 2;
    1880    12164631 :   gel(rem,i) = p1;
    1881    16703919 :   for (i--; i>=0; i--)
    1882             :   {
    1883     4539302 :     av1=avma; p1 = gel(x,i);
    1884    12907575 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1885     4538867 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1886     4539062 :     gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
    1887             :   }
    1888    12164617 :   rem -= 2;
    1889    12164617 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1890    12164816 :   if (pr == ONLY_REM) return gerepileupto(av,rem);
    1891     6638913 :   z -= 2;
    1892             :   {
    1893     6638913 :     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
    1894     6638913 :     gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
    1895             :   }
    1896             : }
    1897             : 
    1898             : GEN
    1899    14049154 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1900             : {
    1901    14049154 :   if (pr == ONLY_REM) return RgX_rem(x, y);
    1902    14049154 :   return RgX_divrem_i(x, y, pr);
    1903             : }
    1904             : 
    1905             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1906             : GEN
    1907      136541 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1908             : {
    1909             :   long vx, dx, dy, dz, i, j, sx, lr;
    1910             :   pari_sp av0, av, tetpil;
    1911             :   GEN z,p1,rem,lead;
    1912             : 
    1913      136541 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1914      136541 :   vx = varn(x);
    1915      136541 :   dx = degpol(x);
    1916      136541 :   dy = degpol(y);
    1917      136541 :   if (dx < dy)
    1918             :   {
    1919       28316 :     if (pr)
    1920             :     {
    1921       28288 :       av0 = avma; x = RgXQX_red(x, T);
    1922       28288 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
    1923       28288 :       if (pr == ONLY_REM) return x;
    1924           0 :       *pr = x;
    1925             :     }
    1926          28 :     return pol_0(vx);
    1927             :   }
    1928      108225 :   lead = leading_coeff(y);
    1929      108225 :   if (!dy) /* y is constant */
    1930             :   {
    1931         602 :     if (pr && pr != ONLY_DIVIDES)
    1932             :     {
    1933           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1934           0 :       *pr = pol_0(vx);
    1935             :     }
    1936         602 :     if (gequal1(lead)) return RgX_copy(x);
    1937           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1938           0 :     return gerepile(av0,tetpil,RgXQX_red(x,T));
    1939             :   }
    1940      107623 :   av0 = avma; dz = dx-dy;
    1941      107623 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1942      107623 :   set_avma(av0);
    1943      107623 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1944      107623 :   x += 2; y += 2; z += 2;
    1945             : 
    1946      107623 :   p1 = gel(x,dx); av = avma;
    1947      107623 :   gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    1948      556146 :   for (i=dx-1; i>=dy; i--)
    1949             :   {
    1950      448523 :     av=avma; p1=gel(x,i);
    1951     2352621 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1952      448503 :     if (lead) p1 = gmul(grem(p1, T), lead);
    1953      448506 :     tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
    1954             :   }
    1955      107623 :   if (!pr) { guncloneNULL(lead); return z-2; }
    1956             : 
    1957      106349 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    1958      179407 :   for (sx=0; ; i--)
    1959             :   {
    1960      179407 :     p1 = gel(x,i);
    1961      618427 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1962      179404 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    1963      100134 :     if (!i) break;
    1964       73058 :     set_avma(av);
    1965             :   }
    1966      106349 :   if (pr == ONLY_DIVIDES)
    1967             :   {
    1968       27034 :     guncloneNULL(lead);
    1969       27034 :     if (sx) return gc_NULL(av0);
    1970       24528 :     return gc_const((pari_sp)rem, z-2);
    1971             :   }
    1972       79315 :   lr=i+3; rem -= lr;
    1973       79315 :   rem[0] = evaltyp(t_POL) | _evallg(lr);
    1974       79315 :   rem[1] = z[-1];
    1975       79315 :   p1 = gerepile((pari_sp)rem,tetpil,p1);
    1976       79315 :   rem += 2; gel(rem,i) = p1;
    1977      169687 :   for (i--; i>=0; i--)
    1978             :   {
    1979       90372 :     av=avma; p1 = gel(x,i);
    1980      306377 :     for (j=0; j<=i && j<=dz; j++)
    1981      216005 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1982       90372 :     tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
    1983             :   }
    1984       79315 :   rem -= 2;
    1985       79315 :   guncloneNULL(lead);
    1986       79315 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1987       79315 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
    1988         161 :   *pr = rem; return z-2;
    1989             : }
    1990             : 
    1991             : /*******************************************************************/
    1992             : /*                                                                 */
    1993             : /*                        PSEUDO-DIVISION                          */
    1994             : /*                                                                 */
    1995             : /*******************************************************************/
    1996             : INLINE GEN
    1997     1318926 : rem(GEN c, GEN T)
    1998             : {
    1999     1318926 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    2000     1318928 :   return c;
    2001             : }
    2002             : 
    2003             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    2004             : int
    2005       17083 : ZXQX_dvd(GEN x, GEN y, GEN T)
    2006             : {
    2007             :   long dx, dy, i, T_ismonic;
    2008       17083 :   pari_sp av = avma, av2;
    2009             :   GEN y_lead;
    2010             : 
    2011       17083 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    2012       17083 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2013       17083 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    2014             :   /* if monic, no point in using pseudo-division */
    2015       17083 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    2016       14570 :   T_ismonic = gequal1(leading_coeff(T));
    2017       14570 :   dx = degpol(x);
    2018       14570 :   if (dx < dy) return !signe(x);
    2019       14570 :   (void)new_chunk(2);
    2020       14570 :   x = RgX_recip_i(x)+2;
    2021       14570 :   y = RgX_recip_i(y)+2;
    2022             :   /* pay attention to sparse divisors */
    2023       29641 :   for (i = 1; i <= dy; i++)
    2024       15071 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    2025       14570 :   av2 = avma;
    2026             :   for (;;)
    2027       72617 :   {
    2028       87187 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    2029       87187 :     x0 = gneg(x0);
    2030       87187 :     m = gcdii(cx, y0);
    2031       87187 :     if (!equali1(m))
    2032             :     {
    2033       84692 :       x0 = gdiv(x0, m);
    2034       84692 :       y0 = diviiexact(y0, m);
    2035       84692 :       if (equali1(y0)) y0 = NULL;
    2036             :     }
    2037      178147 :     for (i=1; i<=dy; i++)
    2038             :     {
    2039       90960 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2040       90960 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    2041       90960 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2042       90960 :       gel(x,i) = c;
    2043             :     }
    2044      686726 :     for (   ; i<=dx; i++)
    2045             :     {
    2046      599539 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2047      599539 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2048      599539 :       gel(x,i) = c;
    2049             :     }
    2050      101523 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    2051       87187 :     if (dx < dy) break;
    2052       72617 :     if (gc_needed(av2,1))
    2053             :     {
    2054          28 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    2055          28 :       gerepilecoeffs(av2,x,dx+1);
    2056             :     }
    2057             :   }
    2058       14570 :   return gc_bool(av, dx < 0);
    2059             : }
    2060             : 
    2061             : /* T either NULL or a t_POL. */
    2062             : GEN
    2063      138508 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    2064             : {
    2065      138508 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    2066      138508 :   pari_sp av = avma, av2;
    2067             :   GEN y_lead;
    2068             : 
    2069      138508 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    2070      138508 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2071             :   /* if monic, no point in using pseudo-division */
    2072      138508 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    2073      131567 :   dx = degpol(x);
    2074      131567 :   if (dx < dy) return RgX_copy(x);
    2075      131567 :   (void)new_chunk(2);
    2076      131567 :   x = RgX_recip_i(x)+2;
    2077      131567 :   y = RgX_recip_i(y)+2;
    2078             :   /* pay attention to sparse divisors */
    2079      376771 :   for (i = 1; i <= dy; i++)
    2080      245204 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2081      131567 :   dz = dx-dy; p = dz+1;
    2082      131567 :   av2 = avma;
    2083             :   for (;;)
    2084             :   {
    2085      362346 :     gel(x,0) = gneg(gel(x,0)); p--;
    2086     1076372 :     for (i=1; i<=dy; i++)
    2087             :     {
    2088      714033 :       GEN c = gmul(y_lead, gel(x,i));
    2089      713998 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    2090      714014 :       gel(x,i) = rem(c, T);
    2091             :     }
    2092      843959 :     for (   ; i<=dx; i++)
    2093             :     {
    2094      481614 :       GEN c = gmul(y_lead, gel(x,i));
    2095      481608 :       gel(x,i) = rem(c, T);
    2096             :     }
    2097      391317 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    2098      362345 :     if (dx < dy) break;
    2099      230779 :     if (gc_needed(av2,1))
    2100             :     {
    2101           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    2102           0 :       gerepilecoeffs(av2,x,dx+1);
    2103             :     }
    2104             :   }
    2105      131566 :   if (dx < 0) return pol_0(vx);
    2106      130383 :   lx = dx+3; x -= 2;
    2107      130383 :   x[0] = evaltyp(t_POL) | _evallg(lx);
    2108      130383 :   x[1] = evalsigne(1) | evalvarn(vx);
    2109      130383 :   x = RgX_recip_i(x);
    2110      130383 :   if (p)
    2111             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    2112       13435 :     GEN t = y_lead;
    2113       13435 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    2114           0 :       t = RgXQ_powu(t, p, T);
    2115             :     else
    2116       13435 :       t = gpowgs(t, p);
    2117       35924 :     for (i=2; i<lx; i++)
    2118             :     {
    2119       22489 :       GEN c = gmul(gel(x,i), t);
    2120       22489 :       gel(x,i) = rem(c,T);
    2121             :     }
    2122       13435 :     if (!T) return gerepileupto(av, x);
    2123             :   }
    2124      116948 :   return gerepilecopy(av, x);
    2125             : }
    2126             : 
    2127             : GEN
    2128      138508 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2129             : 
    2130             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2131             : GEN
    2132       11888 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2133             : {
    2134       11888 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2135       11888 :   pari_sp av = avma, av2;
    2136             :   GEN z, r, ypow, y_lead;
    2137             : 
    2138       11888 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2139       11888 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2140       11888 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2141        7526 :   dx = degpol(x);
    2142        7526 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2143        7526 :   if (dx == dy)
    2144             :   {
    2145          98 :     GEN x_lead = gel(x,lg(x)-1);
    2146          98 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2147          98 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2148          98 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2149          98 :     *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
    2150             :   }
    2151        7428 :   (void)new_chunk(2);
    2152        7428 :   x = RgX_recip_i(x)+2;
    2153        7428 :   y = RgX_recip_i(y)+2;
    2154             :   /* pay attention to sparse divisors */
    2155       38664 :   for (i = 1; i <= dy; i++)
    2156       31236 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2157        7428 :   dz = dx-dy; p = dz+1;
    2158        7428 :   lz = dz+3;
    2159        7428 :   z = cgetg(lz, t_POL);
    2160        7428 :   z[1] = evalsigne(1) | evalvarn(vx);
    2161       28521 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2162        7428 :   ypow = new_chunk(dz+1);
    2163        7428 :   gel(ypow,0) = gen_1;
    2164        7428 :   gel(ypow,1) = y_lead;
    2165       13665 :   for (i=2; i<=dz; i++)
    2166             :   {
    2167        6237 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2168        6237 :     gel(ypow,i) = rem(c,T);
    2169             :   }
    2170        7428 :   av2 = avma;
    2171        7428 :   for (iz=2;;)
    2172             :   {
    2173       15549 :     p--;
    2174       15549 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2175       69229 :     for (i=1; i<=dy; i++)
    2176             :     {
    2177       53680 :       GEN c = gmul(y_lead, gel(x,i));
    2178       53680 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2179       53680 :       gel(x,i) = rem(c, T);
    2180             :     }
    2181       40932 :     for (   ; i<=dx; i++)
    2182             :     {
    2183       25383 :       GEN c = gmul(y_lead, gel(x,i));
    2184       25383 :       gel(x,i) = rem(c,T);
    2185             :     }
    2186       15549 :     x++; dx--;
    2187       21093 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2188       15549 :     if (dx < dy) break;
    2189        8121 :     if (gc_needed(av2,1))
    2190             :     {
    2191           0 :       GEN X = x-2;
    2192           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2193           0 :       X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
    2194           0 :       gerepileall(av2,2, &X, &z); x = X+2;
    2195             :     }
    2196             :   }
    2197       13938 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2198        7428 :   if (dx < 0)
    2199         182 :     x = pol_0(vx);
    2200             :   else
    2201             :   {
    2202        7246 :     lx = dx+3; x -= 2;
    2203        7246 :     x[0] = evaltyp(t_POL) | _evallg(lx);
    2204        7246 :     x[1] = evalsigne(1) | evalvarn(vx);
    2205        7246 :     x = RgX_recip_i(x);
    2206             :   }
    2207        7428 :   z = RgX_recip_i(z);
    2208        7428 :   r = x;
    2209        7428 :   if (p)
    2210             :   {
    2211        3311 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2212        3311 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2213             :   }
    2214        7428 :   *ptr = r; return gc_all(av, 2, &z, ptr);
    2215             : }
    2216             : GEN
    2217       11650 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2218       11650 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2219             : 
    2220             : GEN
    2221           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2222           0 : { return RgXQX_red(RgX_mul(x,y), T); }
    2223             : GEN
    2224   497700009 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
    2225             : GEN
    2226      141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
    2227             : GEN
    2228       26145 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
    2229             : GEN
    2230          35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
    2231             : GEN
    2232         133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
    2233             : 
    2234             : GEN
    2235           0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
    2236             : 
    2237             : GEN
    2238           0 : RgXQX_powers(GEN P, long n, GEN T)
    2239             : {
    2240           0 :   GEN v = cgetg(n+2, t_VEC);
    2241             :   long i;
    2242           0 :   gel(v, 1) = pol_1(varn(T));
    2243           0 :   if (n==0) return v;
    2244           0 :   gel(v, 2) = gcopy(P);
    2245           0 :   for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
    2246           0 :   return v;
    2247             : }
    2248             : 
    2249             : static GEN
    2250      603801 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
    2251             : static GEN
    2252           0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
    2253             : static GEN
    2254       65142 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2255             : static GEN
    2256       86452 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
    2257             : static GEN
    2258      270600 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2259             : static GEN
    2260     1004461 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2261             : static GEN
    2262      978219 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2263             : static GEN
    2264         840 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2265             : static GEN
    2266      703299 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
    2267             : 
    2268             : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
    2269             :               _mul, _sqr, _one, _zero };
    2270             : 
    2271             : GEN
    2272           0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
    2273             : {
    2274           0 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
    2275             : }
    2276             : 
    2277             : GEN
    2278      400528 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2279             : {
    2280      400528 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2281      400527 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
    2282             : }
    2283             : 
    2284             : /* mod X^n */
    2285             : struct modXn {
    2286             :   long v; /* varn(X) */
    2287             :   long n;
    2288             : } ;
    2289             : static GEN
    2290       11893 : _sqrXn(void *data, GEN x) {
    2291       11893 :   struct modXn *S = (struct modXn*)data;
    2292       11893 :   return RgXn_sqr(x, S->n);
    2293             : }
    2294             : static GEN
    2295        4528 : _mulXn(void *data, GEN x, GEN y) {
    2296        4528 :   struct modXn *S = (struct modXn*)data;
    2297        4528 :   return RgXn_mul(x,y, S->n);
    2298             : }
    2299             : static GEN
    2300        1939 : _oneXn(void *data) {
    2301        1939 :   struct modXn *S = (struct modXn*)data;
    2302        1939 :   return pol_1(S->v);
    2303             : }
    2304             : static GEN
    2305           0 : _zeroXn(void *data) {
    2306           0 :   struct modXn *S = (struct modXn*)data;
    2307           0 :   return pol_0(S->v);
    2308             : }
    2309             : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
    2310             :                                           _oneXn, _zeroXn };
    2311             : 
    2312             : GEN
    2313         483 : RgXn_powers(GEN x, long m, long n)
    2314             : {
    2315         483 :   long d = degpol(x);
    2316         483 :   int use_sqr = (d<<1) >= n;
    2317             :   struct modXn S;
    2318         483 :   S.v = varn(x); S.n = n;
    2319         483 :   return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
    2320             : }
    2321             : 
    2322             : GEN
    2323        2286 : RgXn_powu_i(GEN x, ulong m, long n)
    2324             : {
    2325             :   struct modXn S;
    2326             :   long v;
    2327        2286 :   if (n == 0) return x;
    2328        2286 :   v = RgX_valrem(x, &x);
    2329        2286 :   if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
    2330        2265 :   S.v = varn(x); S.n = n;
    2331        2265 :   x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2332        2265 :   if (v) x = RgX_shift_shallow(x, m * v);
    2333        2265 :   return x;
    2334             : }
    2335             : GEN
    2336           0 : RgXn_powu(GEN x, ulong m, long n)
    2337             : {
    2338             :   pari_sp av;
    2339           0 :   if (n == 0) return gcopy(x);
    2340           0 :   av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
    2341             : }
    2342             : 
    2343             : GEN
    2344         966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
    2345             : {
    2346             :   struct modXn S;
    2347         966 :   S.v = varn(gel(x,2)); S.n = n;
    2348         966 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
    2349             : }
    2350             : 
    2351             : GEN
    2352           0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
    2353             : {
    2354           0 :   int use_sqr = 2*degpol(x) >= n;
    2355             :   struct modXn S;
    2356           0 :   S.v = varn(x); S.n = n;
    2357           0 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2358             : }
    2359             : 
    2360             : /* Q(x) mod t^n, x in R[t], n >= 1 */
    2361             : GEN
    2362        5092 : RgXn_eval(GEN Q, GEN x, long n)
    2363             : {
    2364        5092 :   long d = degpol(x);
    2365             :   int use_sqr;
    2366             :   struct modXn S;
    2367        5092 :   if (d == 1 && isrationalzero(gel(x,2)))
    2368             :   {
    2369        5085 :     GEN y = RgX_unscale(Q, gel(x,3));
    2370        5085 :     setvarn(y, varn(x)); return y;
    2371             :   }
    2372           7 :   S.v = varn(x);
    2373           7 :   S.n = n;
    2374           7 :   use_sqr = (d<<1) >= n;
    2375           7 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2376             : }
    2377             : 
    2378             : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
    2379             : static GEN
    2380     2602526 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2381             : {
    2382     2602526 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2383     2602526 :   return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
    2384             : }
    2385             : 
    2386             : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
    2387             : static GEN
    2388          14 : RgXn_sqrhigh(GEN f, long n2, long n)
    2389             : {
    2390          14 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2391          14 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2392             : }
    2393             : 
    2394             : static GEN
    2395     2180798 : RgXn_div_gen(GEN g, GEN f, long e)
    2396             : {
    2397             :   pari_sp av;
    2398             :   ulong mask;
    2399             :   GEN W, a;
    2400     2180798 :   long v = varn(f), n = 1;
    2401             : 
    2402     2180798 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2403     2180791 :   a = ginv(gel(f,2));
    2404     2180792 :   if (e == 1 && !g) return scalarpol(a, v);
    2405     2176036 :   else if (e == 2 && !g)
    2406             :   {
    2407             :     GEN b;
    2408      835472 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2409      244793 :     b = gneg(b);
    2410      244792 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2411      244793 :     return deg1pol(b, a, v);
    2412             :   }
    2413     1340564 :   av = avma;
    2414     1340564 :   W = scalarpol_shallow(a,v);
    2415     1340564 :   mask = quadratic_prec_mask(e);
    2416     3937224 :   while (mask > 1)
    2417             :   {
    2418             :     GEN u, fr;
    2419     2596660 :     long n2 = n;
    2420     2596660 :     n<<=1; if (mask & 1) n--;
    2421     2596660 :     mask >>= 1;
    2422     2596660 :     fr = RgXn_red_shallow(f, n);
    2423     2596660 :     if (mask>1 || !g)
    2424             :     {
    2425     1317299 :       u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2426     1317299 :       W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2427             :     }
    2428             :     else
    2429             :     {
    2430     1279361 :       GEN y = RgXn_mul(g, W, n), yt =  RgXn_red_shallow(y, n-n2);
    2431     1279361 :       u = RgXn_mul(yt, RgXn_mulhigh(fr,  W, n2, n), n-n2);
    2432     1279361 :       W = RgX_sub(y, RgX_shift_shallow(u, n2));
    2433             :     }
    2434     2596660 :     if (gc_needed(av,2))
    2435             :     {
    2436           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2437           0 :       W = gerepileupto(av, W);
    2438             :     }
    2439             :   }
    2440     1340564 :   return W;
    2441             : }
    2442             : 
    2443             : static GEN
    2444         119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
    2445             : {
    2446             :   GEN r;
    2447         119 :   if (lgefint(p) == 3)
    2448             :   {
    2449         119 :     ulong pp = uel(p, 2);
    2450         119 :     if (pp == 2)
    2451          14 :       r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
    2452             :     else
    2453         105 :       r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
    2454             :   }
    2455             :   else
    2456           0 :     r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
    2457         119 :   return FpX_to_mod(r, p);
    2458             : }
    2459             : 
    2460             : static GEN
    2461           0 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
    2462             : {
    2463           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2464           0 :   if (signe(T) == 0) pari_err_OP("/", x, y);
    2465           0 :   r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
    2466           0 :   return FpXQX_to_mod(r, T, p);
    2467             : }
    2468             : 
    2469             : static GEN
    2470          91 : RgXn_inv_FpX(GEN x, long e, GEN p)
    2471             : {
    2472             :   GEN r;
    2473          91 :   if (lgefint(p) == 3)
    2474             :   {
    2475          91 :     ulong pp = uel(p, 2);
    2476          91 :     if (pp == 2)
    2477          28 :       r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
    2478             :     else
    2479          63 :       r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
    2480             :   }
    2481             :   else
    2482           0 :     r = FpXn_inv(RgX_to_FpX(x, p), e, p);
    2483          91 :   return FpX_to_mod(r, p);
    2484             : }
    2485             : 
    2486             : static GEN
    2487           0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
    2488             : {
    2489           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2490           0 :   if (signe(T) == 0) pari_err_OP("/", gen_1, x);
    2491           0 :   r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
    2492           0 :   return FpXQX_to_mod(r, T, p);
    2493             : }
    2494             : 
    2495             : #define code(t1,t2) ((t1 << 6) | t2)
    2496             : 
    2497             : static GEN
    2498      901530 : RgXn_inv_fast(GEN x, long e)
    2499             : {
    2500             :   GEN p, pol;
    2501             :   long pa;
    2502      901530 :   long t = RgX_type(x,&p,&pol,&pa);
    2503      901528 :   switch(t)
    2504             :   {
    2505          91 :     case t_INTMOD: return RgXn_inv_FpX(x, e, p);
    2506           0 :     case code(t_POLMOD, t_INTMOD):
    2507           0 :                    return RgXn_inv_FpXQX(x, e, pol, p);
    2508      901437 :     default:       return NULL;
    2509             :   }
    2510             : }
    2511             : 
    2512             : static GEN
    2513     1279480 : RgXn_div_fast(GEN x, GEN y, long e)
    2514             : {
    2515             :   GEN p, pol;
    2516             :   long pa;
    2517     1279480 :   long t = RgX_type2(x,y,&p,&pol,&pa);
    2518     1279480 :   switch(t)
    2519             :   {
    2520         119 :     case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
    2521           0 :     case code(t_POLMOD, t_INTMOD):
    2522           0 :                    return RgXn_div_FpXQX(x, y, e, pol, p);
    2523     1279361 :     default:       return NULL;
    2524             :   }
    2525             : }
    2526             : #undef code
    2527             : 
    2528             : GEN
    2529     1279480 : RgXn_div_i(GEN g, GEN f, long e)
    2530             : {
    2531     1279480 :   GEN h = RgXn_div_fast(g, f, e);
    2532     1279480 :   if (h) return h;
    2533     1279361 :   return RgXn_div_gen(g, f, e);
    2534             : }
    2535             : 
    2536             : GEN
    2537      583134 : RgXn_div(GEN g, GEN f, long e)
    2538             : {
    2539      583134 :   pari_sp av = avma;
    2540      583134 :   return gerepileupto(av, RgXn_div_i(g, f, e));
    2541             : }
    2542             : 
    2543             : GEN
    2544      901530 : RgXn_inv_i(GEN f, long e)
    2545             : {
    2546      901530 :   GEN h = RgXn_inv_fast(f, e);
    2547      901527 :   if (h) return h;
    2548      901436 :   return RgXn_div_gen(NULL, f, e);
    2549             : }
    2550             : 
    2551             : GEN
    2552      807453 : RgXn_inv(GEN f, long e)
    2553             : {
    2554      807453 :   pari_sp av = avma;
    2555      807453 :   return gerepileupto(av, RgXn_inv_i(f, e));
    2556             : }
    2557             : 
    2558             : /* Compute intformal(x^n*S)/x^(n+1) */
    2559             : static GEN
    2560       58716 : RgX_integXn(GEN x, long n)
    2561             : {
    2562       58716 :   long i, lx = lg(x);
    2563             :   GEN y;
    2564       58716 :   if (lx == 2) return RgX_copy(x);
    2565       52682 :   y = cgetg(lx, t_POL); y[1] = x[1];
    2566      110481 :   for (i=2; i<lx; i++)
    2567       57799 :     gel(y,i) = gdivgs(gel(x,i), n+i-1);
    2568       52682 :   return RgX_renormalize_lg(y, lx);;
    2569             : }
    2570             : 
    2571             : GEN
    2572       52856 : RgXn_expint(GEN h, long e)
    2573             : {
    2574       52856 :   pari_sp av = avma, av2;
    2575       52856 :   long v = varn(h), n;
    2576       52856 :   GEN f = pol_1(v), g;
    2577             :   ulong mask;
    2578             : 
    2579       52856 :   if (!signe(h)) return f;
    2580       52856 :   g = pol_1(v);
    2581       52856 :   n = 1; mask = quadratic_prec_mask(e);
    2582       52856 :   av2 = avma;
    2583       58715 :   for (;mask>1;)
    2584             :   {
    2585             :     GEN u, w;
    2586       58714 :     long n2 = n;
    2587       58714 :     n<<=1; if (mask & 1) n--;
    2588       58714 :     mask >>= 1;
    2589       58714 :     u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
    2590       58716 :     u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
    2591       58716 :     w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
    2592       58715 :     f = RgX_add(f, RgX_shift_shallow(w, n2));
    2593       58714 :     if (mask<=1) break;
    2594        5859 :     u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
    2595        5859 :     g = RgX_sub(g, RgX_shift_shallow(u, n2));
    2596        5859 :     if (gc_needed(av2,2))
    2597             :     {
    2598           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
    2599           0 :       gerepileall(av2, 2, &f, &g);
    2600             :     }
    2601             :   }
    2602       52856 :   return gerepileupto(av, f);
    2603             : }
    2604             : 
    2605             : GEN
    2606           0 : RgXn_exp(GEN h, long e)
    2607             : {
    2608           0 :   long d = degpol(h);
    2609           0 :   if (d < 0) return pol_1(varn(h));
    2610           0 :   if (!d || !gequal0(gel(h,2)))
    2611           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2612           0 :   return RgXn_expint(RgX_deriv(h), e);
    2613             : }
    2614             : 
    2615             : GEN
    2616         154 : RgXn_reverse(GEN f, long e)
    2617             : {
    2618         154 :   pari_sp av = avma, av2;
    2619             :   ulong mask;
    2620             :   GEN fi, a, df, W, an;
    2621         154 :   long v = varn(f), n=1;
    2622         154 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2623           0 :     pari_err_INV("serreverse",f);
    2624         154 :   fi = ginv(gel(f,3));
    2625         154 :   a = deg1pol_shallow(fi,gen_0,v);
    2626         154 :   if (e <= 2) return gerepilecopy(av, a);
    2627         133 :   W = scalarpol(fi,v);
    2628         133 :   df = RgX_deriv(f);
    2629         133 :   mask = quadratic_prec_mask(e);
    2630         133 :   av2 = avma;
    2631         616 :   for (;mask>1;)
    2632             :   {
    2633             :     GEN u, fa, fr;
    2634         483 :     long n2 = n, rt;
    2635         483 :     n<<=1; if (mask & 1) n--;
    2636         483 :     mask >>= 1;
    2637         483 :     fr = RgXn_red_shallow(f, n);
    2638         483 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2639         483 :     an = RgXn_powers(a, rt, n);
    2640         483 :     if (n>1)
    2641             :     {
    2642         483 :       long n4 = (n2+1)>>1;
    2643         483 :       GEN dfr = RgXn_red_shallow(df, n2);
    2644         483 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2645         483 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2646         483 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2647             :     }
    2648         483 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2649         483 :     fa = RgX_shift(fa, -n2);
    2650         483 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2651         483 :     if (gc_needed(av2,2))
    2652             :     {
    2653           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2654           0 :       gerepileall(av2, 2, &a, &W);
    2655             :     }
    2656             :   }
    2657         133 :   return gerepileupto(av, a);
    2658             : }
    2659             : 
    2660             : GEN
    2661           7 : RgXn_sqrt(GEN h, long e)
    2662             : {
    2663           7 :   pari_sp av = avma, av2;
    2664           7 :   long v = varn(h), n = 1;
    2665           7 :   GEN f = scalarpol(gen_1, v), df = f;
    2666           7 :   ulong mask = quadratic_prec_mask(e);
    2667           7 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2668           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2669           7 :   av2 = avma;
    2670             :   while(1)
    2671           7 :   {
    2672          14 :     long n2 = n, m;
    2673             :     GEN g;
    2674          14 :     n<<=1; if (mask & 1) n--;
    2675          14 :     mask >>= 1;
    2676          14 :     m = n-n2;
    2677          14 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2678          14 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2679          14 :     if (mask==1) return gerepileupto(av, f);
    2680           7 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2681           7 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2682           7 :     if (gc_needed(av2,2))
    2683             :     {
    2684           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2685           0 :       gerepileall(av2, 2, &f, &df);
    2686             :     }
    2687             :   }
    2688             : }
    2689             : 
    2690             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2691             : GEN
    2692      115505 : RgXQ_powu(GEN x, ulong n, GEN T)
    2693             : {
    2694      115505 :   pari_sp av = avma;
    2695             : 
    2696      115505 :   if (!n) return pol_1(varn(x));
    2697       67766 :   if (n == 1) return RgX_copy(x);
    2698       19670 :   x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
    2699       19670 :   return gerepilecopy(av, x);
    2700             : }
    2701             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2702             : GEN
    2703      100165 : RgXQ_pow(GEN x, GEN n, GEN T)
    2704             : {
    2705             :   pari_sp av;
    2706      100165 :   long s = signe(n);
    2707             : 
    2708      100165 :   if (!s) return pol_1(varn(x));
    2709      100165 :   if (is_pm1(n) == 1)
    2710       86452 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2711       13713 :   av = avma;
    2712       13713 :   if (s < 0) x = RgXQ_inv(x, T);
    2713       13713 :   x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
    2714       13713 :   return gerepilecopy(av, x);
    2715             : }
    2716             : static GEN
    2717      213627 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
    2718             : static GEN
    2719      113806 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
    2720             : 
    2721             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2722             : GEN
    2723       10906 : ZXQ_powers(GEN x, long l, GEN T)
    2724             : {
    2725       10906 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2726       10906 :   return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
    2727             : }
    2728             : 
    2729             : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
    2730             : GEN
    2731      188475 : ZXQ_powu(GEN x, ulong n, GEN T)
    2732             : {
    2733      188475 :   pari_sp av = avma;
    2734             : 
    2735      188475 :   if (!n) return pol_1(varn(x));
    2736      188475 :   if (n == 1) return ZX_copy(x);
    2737      124983 :   x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
    2738      124988 :   return gerepilecopy(av, x);
    2739             : }
    2740             : 
    2741             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2742             : GEN
    2743        4214 : RgXQ_powers(GEN x, long l, GEN T)
    2744             : {
    2745        4214 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2746        4214 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2747             : }
    2748             : 
    2749             : GEN
    2750        7013 : RgXQV_factorback(GEN L, GEN e, GEN T)
    2751             : {
    2752        7013 :   return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
    2753             : }
    2754             : 
    2755             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2756             : GEN
    2757        8344 : QXQ_powers(GEN a, long n, GEN T)
    2758             : {
    2759             :   GEN den, v;
    2760        8344 :   if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
    2761        8323 :   v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
    2762             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2763        8323 :   if (den)
    2764             :   { /* restore denominators */
    2765        5222 :     GEN d = den;
    2766             :     long i;
    2767        5222 :     gel(v,2) = a;
    2768       26502 :     for (i=3; i<=n+1; i++) {
    2769       21280 :       d = mulii(d,den);
    2770       21280 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2771             :     }
    2772             :   }
    2773        8323 :   return v;
    2774             : }
    2775             : 
    2776             : static GEN
    2777        2492 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2778             : {
    2779        2492 :   long l, i, m = 0;
    2780             :   GEN dz, z;
    2781        2492 :   GEN V = cgetg_copy(v, &l);
    2782        8701 :   for (i = imin; i < l; i++)
    2783             :   {
    2784        6209 :     GEN c = gel(v, i);
    2785        6209 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2786             :   }
    2787        2492 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2788        2723 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2789        8701 :   for (i = imin; i < l; i++)
    2790             :   {
    2791        6209 :     GEN c = gel(v,i);
    2792        6209 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2793        6209 :     gel(V,i) = c;
    2794             :   }
    2795        2492 :   return V;
    2796             : }
    2797             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2798             : GEN
    2799        2261 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2800        2261 : { return do_QXQ_eval(v, 1, a, T); }
    2801             : 
    2802             : GEN
    2803         231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
    2804         231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2805             : 
    2806             : GEN
    2807        2891 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2808             : {
    2809        2891 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2810             : }
    2811             : 
    2812             : GEN
    2813        5179 : RgXQ_norm(GEN x, GEN T)
    2814             : {
    2815             :   pari_sp av;
    2816        5179 :   long dx = degpol(x);
    2817             :   GEN L, y;
    2818        5179 :   if (degpol(T)==0) return gpowgs(x,0);
    2819        5172 :   av = avma; y = resultant(T, x);
    2820        5172 :   L = leading_coeff(T);
    2821        5172 :   if (gequal1(L) || !signe(x)) return y;
    2822           0 :   return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
    2823             : }
    2824             : 
    2825             : GEN
    2826         455 : RgXQ_trace(GEN x, GEN T)
    2827             : {
    2828         455 :   pari_sp av = avma;
    2829             :   GEN dT, z;
    2830             :   long n;
    2831         455 :   if (degpol(T)==0) return gmulgs(x,0);
    2832         448 :   dT = RgX_deriv(T); n = degpol(dT);
    2833         448 :   z = RgXQ_mul(x, dT, T);
    2834         448 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2835         399 :   return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
    2836             : }
    2837             : 
    2838             : GEN
    2839     2775644 : RgX_blocks(GEN P, long n, long m)
    2840             : {
    2841     2775644 :   GEN z = cgetg(m+1,t_VEC);
    2842     2775644 :   long i,j, k=2, l = lg(P);
    2843     8520510 :   for(i=1; i<=m; i++)
    2844             :   {
    2845     5744866 :     GEN zi = cgetg(n+2,t_POL);
    2846     5744866 :     zi[1] = P[1];
    2847     5744866 :     gel(z,i) = zi;
    2848    16945748 :     for(j=2; j<n+2; j++)
    2849    11200882 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2850     5744866 :     zi = RgX_renormalize_lg(zi, n+2);
    2851             :   }
    2852     2775644 :   return z;
    2853             : }
    2854             : 
    2855             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2856             : void
    2857    47588428 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2858             : {
    2859    47588428 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2860             :   GEN p0, p1;
    2861             : 
    2862    47588596 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2863             : 
    2864    47495627 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2865    47495627 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2866    47497653 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2867   140935286 :   for (i=0; i<n1; i++)
    2868             :   {
    2869    93437237 :     p0[2+i] = p[2+(i<<1)];
    2870    93437237 :     p1[2+i] = p[3+(i<<1)];
    2871             :   }
    2872    47498049 :   if (n1 != n0)
    2873    17417224 :     p0[2+i] = p[2+(i<<1)];
    2874    47498049 :   *pe = normalizepol(p0);
    2875    47501944 :   *po = normalizepol(p1);
    2876             : }
    2877             : 
    2878             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2879             : GEN
    2880       43547 : RgX_splitting(GEN p, long k)
    2881             : {
    2882       43547 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2883             :   GEN r;
    2884             : 
    2885       43547 :   m = n/k;
    2886       43547 :   r = cgetg(k+1,t_VEC);
    2887      233807 :   for(i=1; i<=k; i++)
    2888             :   {
    2889      190260 :     gel(r,i) = cgetg(m+3, t_POL);
    2890      190260 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2891             :   }
    2892      571487 :   for (j=1, i=0, l=2; i<=n; i++)
    2893             :   {
    2894      527940 :     gmael(r,j,l) = gel(p,2+i);
    2895      527940 :     if (j==k) { j=1; l++; } else j++;
    2896             :   }
    2897      233807 :   for(i=1; i<=k; i++)
    2898      190260 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2899       43547 :   return r;
    2900             : }
    2901             : 
    2902             : /*******************************************************************/
    2903             : /*                                                                 */
    2904             : /*                        Kronecker form                           */
    2905             : /*                                                                 */
    2906             : /*******************************************************************/
    2907             : 
    2908             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2909             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2910             :  * normalized coefficients */
    2911             : GEN
    2912      178559 : Kronecker_to_mod(GEN z, GEN T)
    2913             : {
    2914      178559 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2915      178559 :   GEN x, t = cgetg(N,t_POL);
    2916      178559 :   t[1] = T[1];
    2917      178559 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2918      178559 :   x[1] = z[1];
    2919      178559 :   T = RgX_copy(T);
    2920     1318040 :   for (i=2; i<lx+2; i++, z+= N-2)
    2921             :   {
    2922     5462080 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2923     1139481 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2924             :   }
    2925      178559 :   N = (l-2) % (N-2) + 2;
    2926      459058 :   for (j=2; j<N; j++) t[j] = z[j];
    2927      178559 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2928      178559 :   return normalizepol_lg(x, i+1);
    2929             : }
    2930             : 
    2931             : /*******************************************************************/
    2932             : /*                                                                 */
    2933             : /*                        Domain detection                         */
    2934             : /*                                                                 */
    2935             : /*******************************************************************/
    2936             : 
    2937             : static GEN
    2938      152783 : zero_FpX_mod(GEN p, long v)
    2939             : {
    2940      152783 :   GEN r = cgetg(3,t_POL);
    2941      152783 :   r[1] = evalvarn(v);
    2942      152783 :   gel(r,2) = mkintmod(gen_0, icopy(p));
    2943      152783 :   return r;
    2944             : }
    2945             : 
    2946             : static GEN
    2947      818061 : RgX_mul_FpX(GEN x, GEN y, GEN p)
    2948             : {
    2949      818061 :   pari_sp av = avma;
    2950             :   GEN r;
    2951      818061 :   if (lgefint(p) == 3)
    2952             :   {
    2953      766705 :     ulong pp = uel(p, 2);
    2954      766705 :     r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
    2955             :                                   RgX_to_Flx(y, pp), pp));
    2956             :   }
    2957             :   else
    2958       51356 :     r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2959      818061 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2960      703001 :   return gerepileupto(av, FpX_to_mod(r, p));
    2961             : }
    2962             : 
    2963             : static GEN
    2964         532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
    2965             : {
    2966         532 :   GEN r = cgetg(3,t_POL);
    2967         532 :   r[1] = evalvarn(v);
    2968         532 :   gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
    2969         532 :   return r;
    2970             : }
    2971             : 
    2972             : static GEN
    2973        2044 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2974             : {
    2975        2044 :   pari_sp av = avma;
    2976             :   long dT;
    2977             :   GEN kx, ky, r;
    2978        2044 :   GEN T = RgX_to_FpX(pol, p);
    2979        2044 :   if (signe(T)==0) pari_err_OP("*", x, y);
    2980        2037 :   dT = degpol(T);
    2981        2037 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    2982        2037 :   ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
    2983        2037 :   r = FpX_mul(kx, ky, p);
    2984        2037 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    2985        1512 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    2986             : }
    2987             : 
    2988             : static GEN
    2989      416946 : RgX_liftred(GEN x, GEN T)
    2990      416946 : { return RgXQX_red(liftpol_shallow(x), T); }
    2991             : 
    2992             : static GEN
    2993      163468 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
    2994             : {
    2995      163468 :   pari_sp av = avma;
    2996      163468 :   long dT = degpol(T);
    2997      163468 :   GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
    2998             :                  RgXX_to_Kronecker(RgX_liftred(y, T), dT));
    2999      163468 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    3000             : }
    3001             : 
    3002             : static GEN
    3003        1365 : RgX_sqr_FpX(GEN x, GEN p)
    3004             : {
    3005        1365 :   pari_sp av = avma;
    3006             :   GEN r;
    3007        1365 :   if (lgefint(p) == 3)
    3008             :   {
    3009        1161 :     ulong pp = uel(p, 2);
    3010        1161 :     r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
    3011             :   }
    3012             :   else
    3013         204 :     r = FpX_sqr(RgX_to_FpX(x, p), p);
    3014        1365 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3015        1239 :   return gerepileupto(av, FpX_to_mod(r, p));
    3016             : }
    3017             : 
    3018             : static GEN
    3019         196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
    3020             : {
    3021         196 :   pari_sp av = avma;
    3022             :   long dT;
    3023         196 :   GEN kx, r, T = RgX_to_FpX(pol, p);
    3024         196 :   if (signe(T)==0) pari_err_OP("*",x,x);
    3025         189 :   dT = degpol(T);
    3026         189 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    3027         189 :   r = FpX_sqr(kx, p);
    3028         189 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3029         189 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    3030             : }
    3031             : 
    3032             : static GEN
    3033       13390 : RgX_sqr_QXQX(GEN x, GEN T)
    3034             : {
    3035       13390 :   pari_sp av = avma;
    3036       13390 :   long dT = degpol(T);
    3037       13390 :   GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
    3038       13390 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    3039             : }
    3040             : 
    3041             : static GEN
    3042       68670 : RgX_rem_FpX(GEN x, GEN y, GEN p)
    3043             : {
    3044       68670 :   pari_sp av = avma;
    3045             :   GEN r;
    3046       68670 :   if (lgefint(p) == 3)
    3047             :   {
    3048       51880 :     ulong pp = uel(p, 2);
    3049       51880 :     r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
    3050             :                                   RgX_to_Flx(y, pp), pp));
    3051             :   }
    3052             :   else
    3053       16790 :     r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    3054       68670 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    3055       31073 :   return gerepileupto(av, FpX_to_mod(r, p));
    3056             : }
    3057             : 
    3058             : static GEN
    3059       38310 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
    3060             : {
    3061       38310 :   pari_sp av = avma;
    3062             :   GEN r;
    3063       38310 :   r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
    3064       38310 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    3065             : }
    3066             : static GEN
    3067          70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    3068             : {
    3069          70 :   pari_sp av = avma;
    3070             :   GEN r;
    3071          70 :   GEN T = RgX_to_FpX(pol, p);
    3072          70 :   if (signe(T) == 0) pari_err_OP("%", x, y);
    3073          63 :   if (lgefint(p) == 3)
    3074             :   {
    3075          55 :     ulong pp = uel(p, 2);
    3076          55 :     GEN Tp = ZX_to_Flx(T, pp);
    3077          55 :     r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
    3078             :                               RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3079             :   }
    3080             :   else
    3081           8 :     r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3082          63 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3083          56 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3084             : }
    3085             : 
    3086             : #define code(t1,t2) ((t1 << 6) | t2)
    3087             : static GEN
    3088    74502695 : RgX_mul_fast(GEN x, GEN y)
    3089             : {
    3090             :   GEN p, pol;
    3091             :   long pa;
    3092    74502695 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3093    74502929 :   switch(t)
    3094             :   {
    3095    38806117 :     case t_INT:    return ZX_mul(x,y);
    3096     3277376 :     case t_FRAC:   return QX_mul(x,y);
    3097      104742 :     case t_FFELT:  return FFX_mul(x, y, pol);
    3098      818061 :     case t_INTMOD: return RgX_mul_FpX(x, y, p);
    3099      163468 :     case code(t_POLMOD, t_INT):
    3100             :     case code(t_POLMOD, t_FRAC):
    3101      163468 :                    return RgX_mul_QXQX(x, y, pol);
    3102        2011 :     case code(t_POLMOD, t_INTMOD):
    3103        2011 :                    return RgX_mul_FpXQX(x, y, pol, p);
    3104    31331154 :     default:       return NULL;
    3105             :   }
    3106             : }
    3107             : static GEN
    3108     1206604 : RgX_sqr_fast(GEN x)
    3109             : {
    3110             :   GEN p, pol;
    3111             :   long pa;
    3112     1206604 :   long t = RgX_type(x,&p,&pol,&pa);
    3113     1206604 :   switch(t)
    3114             :   {
    3115     1149596 :     case t_INT:    return ZX_sqr(x);
    3116       30710 :     case t_FRAC:   return QX_sqr(x);
    3117        2534 :     case t_FFELT:  return FFX_sqr(x, pol);
    3118        1365 :     case t_INTMOD: return RgX_sqr_FpX(x, p);
    3119       13390 :     case code(t_POLMOD, t_INT):
    3120             :     case code(t_POLMOD, t_FRAC):
    3121       13390 :                    return RgX_sqr_QXQX(x, pol);
    3122         196 :     case code(t_POLMOD, t_INTMOD):
    3123         196 :                    return RgX_sqr_FpXQX(x, pol, p);
    3124        8813 :     default:       return NULL;
    3125             :   }
    3126             : }
    3127             : 
    3128             : static GEN
    3129    13878566 : RgX_rem_fast(GEN x, GEN y)
    3130             : {
    3131             :   GEN p, pol;
    3132             :   long pa;
    3133    13878566 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3134    13878669 :   switch(t)
    3135             :   {
    3136     3243972 :     case t_INT:    return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
    3137     1688264 :     case t_FRAC:   return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
    3138          84 :     case t_FFELT:  return FFX_rem(x, y, pol);
    3139       68670 :     case t_INTMOD: return RgX_rem_FpX(x, y, p);
    3140       38310 :     case code(t_POLMOD, t_INT):
    3141             :     case code(t_POLMOD, t_FRAC):
    3142       38310 :                    return RgX_rem_QXQX(x, y, pol);
    3143          64 :     case code(t_POLMOD, t_INTMOD):
    3144          64 :                    return RgX_rem_FpXQX(x, y, pol, p);
    3145     8839305 :     default:       return NULL;
    3146             :   }
    3147             : }
    3148             : 
    3149             : #undef code
    3150             : 
    3151             : GEN
    3152    60828521 : RgX_mul(GEN x, GEN y)
    3153             : {
    3154    60828521 :   GEN z = RgX_mul_fast(x,y);
    3155    60828633 :   if (!z) z = RgX_mul_i(x,y);
    3156    60828275 :   return z;
    3157             : }
    3158             : 
    3159             : GEN
    3160     1194452 : RgX_sqr(GEN x)
    3161             : {
    3162     1194452 :   GEN z = RgX_sqr_fast(x);
    3163     1194445 :   if (!z) z = RgX_sqr_i(x);
    3164     1194445 :   return z;
    3165             : }
    3166             : 
    3167             : GEN
    3168    13878583 : RgX_rem(GEN x, GEN y)
    3169             : {
    3170    13878583 :   GEN z = RgX_rem_fast(x, y);
    3171    13878654 :   if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
    3172    13878569 :   return z;
    3173             : }
    3174             : 
    3175             : static GEN
    3176           0 : _RgX_mul(void* E, GEN x, GEN y)
    3177           0 : { (void) E; return RgX_mul(x, y); }
    3178             : 
    3179             : GEN
    3180           0 : RgXV_prod(GEN V)
    3181           0 : { return gen_product(V, NULL, &_RgX_mul); }
    3182             : 
    3183             : GEN
    3184    11012927 : RgXn_mul(GEN f, GEN g, long n)
    3185             : {
    3186    11012927 :   pari_sp av = avma;
    3187    11012927 :   GEN h = RgX_mul_fast(f,g);
    3188    11012927 :   if (!h) return RgXn_mul2(f,g,n);
    3189     1774925 :   if (degpol(h) < n) return h;
    3190      376108 :   return gerepilecopy(av, RgXn_red_shallow(h, n));
    3191             : }
    3192             : 
    3193             : GEN
    3194       12152 : RgXn_sqr(GEN f, long n)
    3195             : {
    3196       12152 :   pari_sp av = avma;
    3197       12152 :   GEN g = RgX_sqr_fast(f);
    3198       12152 :   if (!g) return RgXn_sqr2(f,n);
    3199       11837 :   if (degpol(g) < n) return g;
    3200       10640 :   return gerepilecopy(av, RgXn_red_shallow(g, n));
    3201             : }
    3202             : 
    3203             : /* (f*g) \/ x^n */
    3204             : GEN
    3205     2661254 : RgX_mulhigh_i(GEN f, GEN g, long n)
    3206             : {
    3207     2661254 :   GEN h = RgX_mul_fast(f,g);
    3208     2661256 :   return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
    3209             : }
    3210             : 
    3211             : /* (f*g) \/ x^n */
    3212             : GEN
    3213           0 : RgX_sqrhigh_i(GEN f, long n)
    3214             : {
    3215           0 :   GEN h = RgX_sqr_fast(f);
    3216           0 :   return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
    3217             : }

Generated by: LCOV version 1.14