Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - RgX.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 26712-590d837a1c) Lines: 1553 1734 89.6 %
Date: 2021-06-22 07:13:04 Functions: 187 205 91.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_pol
      19             : 
      20             : /*******************************************************************/
      21             : /*                                                                 */
      22             : /*                         GENERIC                                 */
      23             : /*                                                                 */
      24             : /*******************************************************************/
      25             : 
      26             : /* Return optimal parameter l for the evaluation of n/m polynomials of degree d
      27             :    Fractional values can be used if the evaluations are done with different
      28             :    accuracies, and thus have different weights.
      29             :  */
      30             : long
      31    19675304 : brent_kung_optpow(long d, long n, long m)
      32             : {
      33             :   long p, r;
      34    19675304 :   long pold=1, rold=n*(d-1);
      35    97994607 :   for(p=2; p<=d; p++)
      36             :   {
      37    78319303 :     r = m*(p-1) + n*((d-1)/p);
      38    78319303 :     if (r<rold) { pold=p; rold=r; }
      39             :   }
      40    19675304 :   return pold;
      41             : }
      42             : 
      43             : static GEN
      44    12652245 : 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    12652245 :   pari_sp av = avma;
      48             :   long i;
      49    12652245 :   GEN z = cmul(E,P,a,ff->one(E));
      50    12611958 :   if (!z) z = gen_0;
      51    63351905 :   for (i=1; i<=n; i++)
      52             :   {
      53    50733823 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      54    50738802 :     if (t) {
      55    39475189 :       z = ff->add(E, z, t);
      56    39460448 :       if (gc_needed(av,2)) z = gerepileupto(av, z);
      57             :     }
      58             :   }
      59    12618082 :   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    11290726 : 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    11290726 :   pari_sp av = avma;
      72    11290726 :   long l = lg(V)-1;
      73             :   GEN z, u;
      74             : 
      75    11290726 :   if (d < 0) return ff->zero(E);
      76    10737915 :   if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      77      983230 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      78      983230 :   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      983230 :   d -= l;
      84      983230 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      85     1912799 :   while (d >= l-1)
      86             :   {
      87      928332 :     d -= l-1;
      88      928332 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      89      928280 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      90      928323 :     if (gc_needed(av,2))
      91          95 :       z = gerepileupto(av, z);
      92             :   }
      93      984467 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      94      984500 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      95      984481 :   return gerepileupto(av, ff->red(E,z));
      96             : }
      97             : 
      98             : GEN
      99      710827 : 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      710827 :   pari_sp av = avma;
     103             :   GEN z, V;
     104             :   long rtd;
     105      710827 :   if (d < 0) return ff->zero(E);
     106      709896 :   rtd = (long) sqrt((double)d);
     107      709896 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     108      709919 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     109      709915 :   return gerepileupto(av, z);
     110             : }
     111             : 
     112             : static GEN
     113     1109376 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     114             : static GEN
     115    14057676 : _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      951542 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     120             : static GEN
     121      311609 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     122             : static GEN
     123     1132011 : _gen_one(void *E) { (void)E; return gen_1; }
     124             : static GEN
     125       24102 : _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      511166 : _gen_cmul(void *E, GEN P, long a, GEN x)
     132      511166 : {(void)E; return gmul(gel(P,a+2), x);}
     133             : 
     134             : GEN
     135      166434 : RgX_RgV_eval(GEN Q, GEN x)
     136             : {
     137      166434 :   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        2954 : RgXV_RgV_eval(GEN Q, GEN x)
     148             : {
     149        2954 :   long i, l = lg(Q), vQ = gvar(Q);
     150        2954 :   GEN v = cgetg(l, t_VEC);
     151      248395 :   for (i = 1; i < l; i++)
     152             :   {
     153      245441 :     GEN Qi = gel(Q, i);
     154      245441 :     gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
     155             :   }
     156        2954 :   return v;
     157             : }
     158             : 
     159             : GEN
     160       12376 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
     161             : {
     162       12376 :   pari_sp av = avma;
     163             :   long d, i;
     164             :   GEN s;
     165       12376 :   if (typ(P)!=t_POL)
     166        1547 :     return mkvec2(P, gen_1);
     167       10829 :   d = degpol(P);
     168       10829 :   s = gel(P, d+2);
     169       54635 :   for (i = d-1; i >= 0; i--)
     170             :   {
     171       43806 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
     172       43806 :     if (gc_needed(av,1))
     173             :     {
     174           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
     175           0 :       s = gerepileupto(av, s);
     176             :     }
     177             :   }
     178       10829 :   s = gerepileupto(av, s);
     179       10829 :   return mkvec2(s, gel(B,d+1));
     180             : }
     181             : 
     182             : GEN
     183        1792 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
     184             : {
     185        1792 :   pari_sp av = avma;
     186        1792 :   long i, d = degpol(P), v = varn(A);
     187             :   GEN s;
     188        1792 :   if (signe(P)==0) return mkvec2(pol_0(v), pol_1(v));
     189        1568 :   s = scalarpol_shallow(gel(P, d+2), v);
     190        5229 :   for (i = d-1; i >= 0; i--)
     191             :   {
     192        3661 :     GEN c = gel(P,i+2), b = gel(B,d+1-i);
     193        3661 :     s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
     194        3661 :     if (gc_needed(av,1))
     195             :     {
     196           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
     197           0 :       s = gerepileupto(av, s);
     198             :     }
     199             :   }
     200        1568 :   s = gerepileupto(av, s);
     201        1568 :   return mkvec2(s, gel(B,d+1));
     202             : }
     203             : 
     204             : const struct bb_algebra *
     205      143824 : get_Rg_algebra(void)
     206             : {
     207      143824 :   return &Rg_algebra;
     208             : }
     209             : 
     210             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     211             : 
     212             : static GEN
     213       17780 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     214             : {
     215             :   (void) E;
     216       17780 :   return RgX_divrem(x, y, r);
     217             : }
     218             : 
     219             : GEN
     220        4382 : RgX_digits(GEN x, GEN T)
     221             : {
     222        4382 :   pari_sp av = avma;
     223        4382 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     224        4382 :   GEN z = gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     225        4382 :   return gerepileupto(av, z);
     226             : }
     227             : 
     228             : /*******************************************************************/
     229             : /*                                                                 */
     230             : /*                         RgX                                     */
     231             : /*                                                                 */
     232             : /*******************************************************************/
     233             : 
     234             : long
     235    23071716 : RgX_equal(GEN x, GEN y)
     236             : {
     237    23071716 :   long i = lg(x);
     238             : 
     239    23071716 :   if (i != lg(y)) return 0;
     240    97496655 :   for (i--; i > 1; i--)
     241    74723277 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     242    22773378 :   return 1;
     243             : }
     244             : 
     245             : /* Returns 1 in the base ring over which x is defined */
     246             : /* HACK: this also works for t_SER */
     247             : GEN
     248     2155819 : Rg_get_1(GEN x)
     249             : {
     250             :   GEN p, T;
     251     2155819 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     252     2155821 :   if (RgX_type_is_composite(tx))
     253       36027 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     254     2155821 :   switch(tx)
     255             :   {
     256         357 :     case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
     257           7 :     case t_PADIC: return cvtop(gen_1, p, lx);
     258         728 :     case t_FFELT: return FF_1(T);
     259     2154729 :     default: return gen_1;
     260             :   }
     261             : }
     262             : /* Returns 0 in the base ring over which x is defined */
     263             : /* HACK: this also works for t_SER */
     264             : GEN
     265     6973645 : Rg_get_0(GEN x)
     266             : {
     267             :   GEN p, T;
     268     6973645 :   long i, lx, tx = Rg_type(x, &p, &T, &lx);
     269     6973645 :   if (RgX_type_is_composite(tx))
     270       51436 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     271     6973645 :   switch(tx)
     272             :   {
     273         441 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     274           0 :     case t_PADIC: return zeropadic(p, lx);
     275         266 :     case t_FFELT: return FF_zero(T);
     276     6972938 :     default: return gen_0;
     277             :   }
     278             : }
     279             : 
     280             : GEN
     281        5131 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     282             : {
     283        5131 :   long i, n = degpol(P);
     284             :   GEN z, dz, dP;
     285        5131 :   if (n < 0) return gen_0;
     286        5131 :   P = Q_remove_denom(P, &dP);
     287        5131 :   z = gel(P,2); if (n == 0) return icopy(z);
     288        3108 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     289        3108 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     290        6181 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     291        3108 :   dz = mul_denom(dP, dV);
     292        3108 :   return dz? RgX_Rg_div(z, dz): z;
     293             : }
     294             : 
     295             : /* Return P(h * x), not memory clean */
     296             : GEN
     297       19408 : RgX_unscale(GEN P, GEN h)
     298             : {
     299       19408 :   long i, l = lg(P);
     300       19408 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     301       19408 :   Q[1] = P[1];
     302       19408 :   if (l == 2) return Q;
     303       19380 :   gel(Q,2) = gcopy(gel(P,2));
     304       90949 :   for (i=3; i<l; i++)
     305             :   {
     306       71569 :     hi = gmul(hi,h);
     307       71569 :     gel(Q,i) = gmul(gel(P,i), hi);
     308             :   }
     309       19380 :   return Q;
     310             : }
     311             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     312             : GEN
     313      515333 : ZX_z_unscale(GEN P, long h)
     314             : {
     315      515333 :   long i, l = lg(P);
     316      515333 :   GEN Q = cgetg(l, t_POL);
     317      515333 :   Q[1] = P[1];
     318      515333 :   if (l == 2) return Q;
     319      515333 :   gel(Q,2) = gel(P,2);
     320      515333 :   if (l == 3) return Q;
     321      515333 :   if (h == -1)
     322      529439 :     for (i = 3; i < l; i++)
     323             :     {
     324      486719 :       gel(Q,i) = negi(gel(P,i));
     325      486717 :       if (++i == l) break;
     326      433896 :       gel(Q,i) = gel(P,i);
     327             :     }
     328             :   else
     329             :   {
     330             :     GEN hi;
     331      419790 :     gel(Q,3) = mulis(gel(P,3), h);
     332      419790 :     hi = sqrs(h);
     333     2129507 :     for (i = 4; i < l; i++)
     334             :     {
     335     1709720 :       gel(Q,i) = mulii(gel(P,i), hi);
     336     1709714 :       if (i != l-1) hi = mulis(hi,h);
     337             :     }
     338             :   }
     339      515328 :   return Q;
     340             : }
     341             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     342             : GEN
     343       59353 : ZX_unscale(GEN P, GEN h)
     344             : {
     345             :   long i, l;
     346             :   GEN Q, hi;
     347       59353 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     348         733 :   l = lg(P); Q = cgetg(l, t_POL);
     349         733 :   Q[1] = P[1];
     350         733 :   if (l == 2) return Q;
     351         733 :   gel(Q,2) = gel(P,2);
     352         733 :   if (l == 3) return Q;
     353         733 :   hi = h;
     354         733 :   gel(Q,3) = mulii(gel(P,3), hi);
     355        2248 :   for (i = 4; i < l; i++)
     356             :   {
     357        1515 :     hi = mulii(hi,h);
     358        1515 :     gel(Q,i) = mulii(gel(P,i), hi);
     359             :   }
     360         733 :   return Q;
     361             : }
     362             : /* P a ZX. Return P(x << n), not memory clean */
     363             : GEN
     364      212749 : ZX_unscale2n(GEN P, long n)
     365             : {
     366      212749 :   long i, ni = n, l = lg(P);
     367      212749 :   GEN Q = cgetg(l, t_POL);
     368      212747 :   Q[1] = P[1];
     369      212747 :   if (l == 2) return Q;
     370      212747 :   gel(Q,2) = gel(P,2);
     371      212747 :   if (l == 3) return Q;
     372      212747 :   gel(Q,3) = shifti(gel(P,3), ni);
     373      791811 :   for (i=4; i<l; i++)
     374             :   {
     375      579103 :     ni += n;
     376      579103 :     gel(Q,i) = shifti(gel(P,i), ni);
     377             :   }
     378      212708 :   return Q;
     379             : }
     380             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     381             : GEN
     382        1197 : ZX_unscale_div(GEN P, GEN h)
     383             : {
     384        1197 :   long i, l = lg(P);
     385        1197 :   GEN hi, Q = cgetg(l, t_POL);
     386        1197 :   Q[1] = P[1];
     387        1197 :   if (l == 2) return Q;
     388        1197 :   gel(Q,2) = diviiexact(gel(P,2), h);
     389        1197 :   if (l == 3) return Q;
     390        1197 :   gel(Q,3) = gel(P,3);
     391        1197 :   if (l == 4) return Q;
     392        1197 :   hi = h;
     393        1197 :   gel(Q,4) = mulii(gel(P,4), hi);
     394        5264 :   for (i=5; i<l; i++)
     395             :   {
     396        4067 :     hi = mulii(hi,h);
     397        4067 :     gel(Q,i) = mulii(gel(P,i), hi);
     398             :   }
     399        1197 :   return Q;
     400             : }
     401             : 
     402             : GEN
     403        4662 : RgXV_unscale(GEN v, GEN h)
     404             : {
     405             :   long i, l;
     406             :   GEN w;
     407        4662 :   if (!h || isint1(h)) return v;
     408        3834 :   w = cgetg_copy(v, &l);
     409       13358 :   for (i=1; i<l; i++) gel(w,i) = RgX_unscale(gel(v,i), h);
     410        3834 :   return w;
     411             : }
     412             : 
     413             : /* Return h^degpol(P) P(x / h), not memory clean */
     414             : GEN
     415     2970139 : RgX_rescale(GEN P, GEN h)
     416             : {
     417     2970139 :   long i, l = lg(P);
     418     2970139 :   GEN Q = cgetg(l,t_POL), hi = h;
     419     2970116 :   gel(Q,l-1) = gel(P,l-1);
     420     9511691 :   for (i=l-2; i>=2; i--)
     421             :   {
     422     9509005 :     gel(Q,i) = gmul(gel(P,i), hi);
     423     9508554 :     if (i == 2) break;
     424     6540936 :     hi = gmul(hi,h);
     425             :   }
     426     2970304 :   Q[1] = P[1]; return Q;
     427             : }
     428             : 
     429             : /* A(X^d) --> A(X) */
     430             : GEN
     431      781606 : RgX_deflate(GEN x0, long d)
     432             : {
     433             :   GEN z, y, x;
     434      781606 :   long i,id, dy, dx = degpol(x0);
     435      781605 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     436      356806 :   dy = dx/d;
     437      356806 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     438      356804 :   z = y + 2;
     439      356804 :   x = x0+ 2;
     440     1540813 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     441      356804 :   return y;
     442             : }
     443             : 
     444             : /* F a t_RFRAC */
     445             : long
     446          70 : rfrac_deflate_order(GEN F)
     447             : {
     448          70 :   GEN N = gel(F,1), D = gel(F,2);
     449          70 :   long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
     450          70 :   if (m == 1) return 1;
     451          28 :   if (typ(N) == t_POL && varn(N) == varn(D))
     452          14 :     m = cgcd(m, RgX_deflate_order(N));
     453          28 :   return m;
     454             : }
     455             : /* F a t_RFRAC */
     456             : GEN
     457          70 : rfrac_deflate_max(GEN F, long *m)
     458             : {
     459          70 :   *m = rfrac_deflate_order(F);
     460          70 :   return rfrac_deflate(F, *m);
     461             : }
     462             : /* F a t_RFRAC */
     463             : GEN
     464          70 : rfrac_deflate(GEN F, long m)
     465             : {
     466          70 :   GEN N = gel(F,1), D = gel(F,2);
     467          70 :   if (m == 1) return F;
     468          28 :   if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
     469          28 :   D = RgX_deflate(D, m); return mkrfrac(N, D);
     470             : }
     471             : 
     472             : /* return x0(X^d) */
     473             : GEN
     474      671532 : RgX_inflate(GEN x0, long d)
     475             : {
     476      671532 :   long i, id, dy, dx = degpol(x0);
     477      671532 :   GEN x = x0 + 2, z, y;
     478      671532 :   if (dx <= 0) return leafcopy(x0);
     479      663958 :   dy = dx*d;
     480      663958 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     481      663958 :   z = y + 2;
     482    22909630 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     483    10561390 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     484      663958 :   return y;
     485             : }
     486             : 
     487             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     488             : static GEN
     489     4029623 : RgX_translate_basecase(GEN P, GEN c)
     490             : {
     491     4029623 :   pari_sp av = avma;
     492             :   GEN Q, R;
     493             :   long i, k, n;
     494             : 
     495     4029623 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     496     4027664 :   Q = leafcopy(P);
     497     4027659 :   R = Q+2; n = degpol(P);
     498     4027656 :   if (isint1(c))
     499             :   {
     500       27017 :     for (i=1; i<=n; i++)
     501             :     {
     502      139596 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
     503       23077 :       if (gc_needed(av,2))
     504             :       {
     505           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
     506           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     507             :       }
     508             :     }
     509             :   }
     510     4023703 :   else if (isintm1(c))
     511             :   {
     512      900973 :     for (i=1; i<=n; i++)
     513             :     {
     514     2889548 :       for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
     515      685148 :       if (gc_needed(av,2))
     516             :       {
     517           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
     518           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     519             :       }
     520             :     }
     521             :   }
     522             :   else
     523             :   {
     524    14440142 :     for (i=1; i<=n; i++)
     525             :     {
     526    39168130 :       for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
     527    10632262 :       if (gc_needed(av,2))
     528             :       {
     529           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
     530           0 :         Q = gerepilecopy(av, Q); R = Q+2;
     531             :       }
     532             :     }
     533             :   }
     534     4026995 :   return gerepilecopy(av, Q);
     535             : }
     536             : GEN
     537     4030043 : RgX_translate(GEN P, GEN c)
     538             : {
     539     4030043 :   pari_sp av = avma;
     540     4030043 :   long n = degpol(P);
     541     4030037 :   if (n < 40)
     542     4029621 :     return RgX_translate_basecase(P, c);
     543             :   else
     544             :   {
     545         416 :     long d = n >> 1;
     546         416 :     GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
     547         418 :     GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
     548         418 :     GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
     549         418 :     return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
     550             :   }
     551             : }
     552             : 
     553             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     554             : GEN
     555       22421 : RgXQX_translate(GEN P, GEN c, GEN T)
     556             : {
     557       22421 :   pari_sp av = avma;
     558             :   GEN Q, R;
     559             :   long i, k, n;
     560             : 
     561       22421 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     562       22092 :   Q = leafcopy(P);
     563       22092 :   R = Q+2; n = degpol(P);
     564       78431 :   for (i=1; i<=n; i++)
     565             :   {
     566      251433 :     for (k=n-i; k<n; k++)
     567             :     {
     568      195094 :       pari_sp av2 = avma;
     569      195094 :       gel(R,k) = gerepileupto(av2,
     570      195094 :                    RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
     571             :     }
     572       56339 :     if (gc_needed(av,2))
     573             :     {
     574           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     575           0 :       Q = gerepilecopy(av, Q); R = Q+2;
     576             :     }
     577             :   }
     578       22092 :   return gerepilecopy(av, Q);
     579             : }
     580             : 
     581             : /********************************************************************/
     582             : /**                                                                **/
     583             : /**                          CONVERSIONS                           **/
     584             : /**                       (not memory clean)                       **/
     585             : /**                                                                **/
     586             : /********************************************************************/
     587             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     588             :  * but everything else is */
     589             : static GEN
     590      157766 : QXQ_to_mod(GEN x, GEN T)
     591             : {
     592             :   long d;
     593      157766 :   switch(typ(x))
     594             :   {
     595       58324 :     case t_INT:  return icopy(x);
     596        1001 :     case t_FRAC: return gcopy(x);
     597       98441 :     case t_POL:
     598       98441 :       d = degpol(x);
     599       98441 :       if (d < 0) return gen_0;
     600       96173 :       if (d == 0) return gcopy(gel(x,2));
     601       94052 :       return mkpolmod(RgX_copy(x), T);
     602           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     603             :              return NULL;/* LCOV_EXCL_LINE */
     604             :   }
     605             : }
     606             : /* pure shallow version */
     607             : GEN
     608      681880 : QXQ_to_mod_shallow(GEN x, GEN T)
     609             : {
     610             :   long d;
     611      681880 :   switch(typ(x))
     612             :   {
     613      439183 :     case t_INT:
     614      439183 :     case t_FRAC: return x;
     615      242697 :     case t_POL:
     616      242697 :       d = degpol(x);
     617      242697 :       if (d < 0) return gen_0;
     618      200074 :       if (d == 0) return gel(x,2);
     619      186397 :       return mkpolmod(x, T);
     620           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     621             :              return NULL;/* LCOV_EXCL_LINE */
     622             :   }
     623             : }
     624             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
     625             :  * Not memory clean because T not copied, but everything else is */
     626             : static GEN
     627       36267 : QXQX_to_mod(GEN z, GEN T)
     628             : {
     629       36267 :   long i,l = lg(z);
     630       36267 :   GEN x = cgetg(l,t_POL);
     631      180999 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     632       36267 :   x[1] = z[1]; return normalizepol_lg(x,l);
     633             : }
     634             : /* pure shallow version */
     635             : GEN
     636      151643 : QXQX_to_mod_shallow(GEN z, GEN T)
     637             : {
     638      151643 :   long i,l = lg(z);
     639      151643 :   GEN x = cgetg(l,t_POL);
     640      754465 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
     641      151643 :   x[1] = z[1]; return normalizepol_lg(x,l);
     642             : }
     643             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     644             : GEN
     645       12313 : QXQXV_to_mod(GEN V, GEN T)
     646             : {
     647       12313 :   long i, l = lg(V);
     648       12313 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     649       48580 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     650       12313 :   return z;
     651             : }
     652             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     653             : GEN
     654       10820 : QXQV_to_mod(GEN V, GEN T)
     655             : {
     656       10820 :   long i, l = lg(V);
     657       10820 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     658       23854 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
     659       10820 :   return z;
     660             : }
     661             : 
     662             : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
     663             : GEN
     664       13300 : QXQC_to_mod_shallow(GEN V, GEN T)
     665             : {
     666       13300 :   long i, l = lg(V);
     667       13300 :   GEN z = cgetg(l, t_COL);
     668       92358 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
     669       13300 :   return z;
     670             : }
     671             : 
     672             : GEN
     673        6223 : QXQM_to_mod_shallow(GEN V, GEN T)
     674             : {
     675        6223 :   long i, l = lg(V);
     676        6223 :   GEN z = cgetg(l, t_MAT);
     677       19523 :   for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
     678        6223 :   return z;
     679             : }
     680             : 
     681             : GEN
     682     4210024 : RgX_renormalize_lg(GEN x, long lx)
     683             : {
     684             :   long i;
     685     6817783 :   for (i = lx-1; i>1; i--)
     686     6439650 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     687     4210024 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     688     4210024 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     689             : }
     690             : 
     691             : GEN
     692     1367287 : RgV_to_RgX(GEN x, long v)
     693             : {
     694     1367287 :   long i, k = lg(x);
     695             :   GEN p;
     696             : 
     697     4022952 :   while (--k && gequal0(gel(x,k)));
     698     1367287 :   if (!k) return pol_0(v);
     699     1360550 :   i = k+2; p = cgetg(i,t_POL);
     700     1360553 :   p[1] = evalsigne(1) | evalvarn(v);
     701     9429968 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     702     1360553 :   return p;
     703             : }
     704             : GEN
     705      168283 : RgV_to_RgX_reverse(GEN x, long v)
     706             : {
     707      168283 :   long j, k, l = lg(x);
     708             :   GEN p;
     709             : 
     710      168822 :   for (k = 1; k < l; k++)
     711      168822 :     if (!gequal0(gel(x,k))) break;
     712      168283 :   if (k == l) return pol_0(v);
     713      168283 :   k -= 1;
     714      168283 :   l -= k;
     715      168283 :   x += k;
     716      168283 :   p = cgetg(l+1,t_POL);
     717      168283 :   p[1] = evalsigne(1) | evalvarn(v);
     718      948202 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     719      168283 :   return p;
     720             : }
     721             : 
     722             : /* return the (N-dimensional) vector of coeffs of p */
     723             : GEN
     724    11600676 : RgX_to_RgC(GEN x, long N)
     725             : {
     726             :   long i, l;
     727             :   GEN z;
     728    11600676 :   l = lg(x)-1; x++;
     729    11600676 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     730    11600676 :   z = cgetg(N+1,t_COL);
     731    73273985 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     732    22294997 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     733    11600812 :   return z;
     734             : }
     735             : GEN
     736     1233927 : Rg_to_RgC(GEN x, long N)
     737             : {
     738     1233927 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     739             : }
     740             : 
     741             : /* vector of polynomials (in v) whose coefs are given by the columns of x */
     742             : GEN
     743      298491 : RgM_to_RgXV(GEN x, long v)
     744     1280722 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
     745             : GEN
     746        2849 : RgM_to_RgXV_reverse(GEN x, long v)
     747       11396 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
     748             : 
     749             : /* matrix whose entries are given by the coeffs of the polynomials in
     750             :  * vector v (considered as degree n-1 polynomials) */
     751             : GEN
     752      283294 : RgV_to_RgM(GEN x, long n)
     753     1515016 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
     754             : 
     755             : GEN
     756       43513 : RgXV_to_RgM(GEN x, long n)
     757      253754 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
     758             : 
     759             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     760             : GEN
     761       17778 : RgM_to_RgXX(GEN x, long v,long w)
     762             : {
     763       17778 :   long j, lx = lg(x);
     764       17778 :   GEN y = cgetg(lx+1, t_POL);
     765       17778 :   y[1] = evalsigne(1) | evalvarn(v);
     766       17778 :   y++;
     767      102035 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     768       17778 :   return normalizepol_lg(--y, lx+1);
     769             : }
     770             : 
     771             : /* matrix whose entries are given by the coeffs of the polynomial v in
     772             :  * two variables (considered as degree n-1 polynomials) */
     773             : GEN
     774         308 : RgXX_to_RgM(GEN v, long n)
     775             : {
     776         308 :   long j, N = lg(v)-1;
     777         308 :   GEN y = cgetg(N, t_MAT);
     778        1001 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     779         308 :   return y;
     780             : }
     781             : 
     782             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     783             : GEN
     784       28670 : RgXY_swapspec(GEN x, long n, long w, long nx)
     785             : {
     786       28670 :   long j, ly = n+3;
     787       28670 :   GEN y = cgetg(ly, t_POL);
     788       28670 :   y[1] = evalsigne(1);
     789      315570 :   for (j=2; j<ly; j++)
     790             :   {
     791             :     long k;
     792      286900 :     GEN a = cgetg(nx+2,t_POL);
     793      286900 :     a[1] = evalsigne(1) | evalvarn(w);
     794     1395223 :     for (k=0; k<nx; k++)
     795             :     {
     796     1108323 :       GEN xk = gel(x,k);
     797     1108323 :       if (typ(xk)==t_POL)
     798     1030722 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     799             :       else
     800       77601 :         gel(a,k+2) = j==2 ? xk: gen_0;
     801             :     }
     802      286900 :     gel(y,j) = normalizepol_lg(a, nx+2);
     803             :   }
     804       28670 :   return normalizepol_lg(y,ly);
     805             : }
     806             : 
     807             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     808             : GEN
     809         952 : RgXY_swap(GEN x, long n, long w)
     810             : {
     811         952 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     812         952 :   setvarn(z, varn(x)); return z;
     813             : }
     814             : 
     815             : long
     816         190 : RgXY_degreex(GEN b)
     817             : {
     818         190 :   long deg = 0, i;
     819         190 :   if (!signe(b)) return -1;
     820        1032 :   for (i = 2; i < lg(b); ++i)
     821             :   {
     822         842 :     GEN bi = gel(b, i);
     823         842 :     if (typ(bi) == t_POL)
     824         799 :       deg = maxss(deg, degpol(bi));
     825             :   }
     826         190 :   return deg;
     827             : }
     828             : 
     829             : GEN
     830        1232 : RgXY_derivx(GEN x)
     831             : {
     832             :   long i,lx;
     833        1232 :   GEN y = cgetg_copy(x,&lx);
     834       16632 :   for (i=2; i<lx ; i++)
     835       15400 :     gel(y,i) = RgX_deriv(gel(x,i));
     836        1232 :   y[1] = x[1]; return normalizepol_lg(y,i);
     837             : }
     838             : 
     839             : /* return (x % X^n). Shallow */
     840             : GEN
     841     4808777 : RgXn_red_shallow(GEN a, long n)
     842             : {
     843     4808777 :   long i, L = n+2, l = lg(a);
     844             :   GEN  b;
     845     4808777 :   if (L >= l) return a; /* deg(x) < n */
     846     3127052 :   b = cgetg(L, t_POL); b[1] = a[1];
     847    23592582 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     848     3127052 :   return normalizepol_lg(b,L);
     849             : }
     850             : 
     851             : GEN
     852         483 : RgXnV_red_shallow(GEN x, long n)
     853        2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
     854             : 
     855             : /* return (x * X^n). Shallow */
     856             : GEN
     857   137385910 : RgX_shift_shallow(GEN a, long n)
     858             : {
     859   137385910 :   long i, l = lg(a);
     860             :   GEN  b;
     861   137385910 :   if (l == 2 || !n) return a;
     862    92989748 :   l += n;
     863    92989748 :   if (n < 0)
     864             :   {
     865    52421010 :     if (l <= 2) return pol_0(varn(a));
     866    51539103 :     b = cgetg(l, t_POL); b[1] = a[1];
     867    51539301 :     a -= n;
     868   153382116 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     869             :   } else {
     870    40568738 :     b = cgetg(l, t_POL); b[1] = a[1];
     871    40569423 :     a -= n; n += 2;
     872    84957683 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     873   407700260 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     874             :   }
     875    92108724 :   return b;
     876             : }
     877             : /* return (x * X^n). */
     878             : GEN
     879     1572450 : RgX_shift(GEN a, long n)
     880             : {
     881     1572450 :   long i, l = lg(a);
     882             :   GEN  b;
     883     1572450 :   if (l == 2 || !n) return RgX_copy(a);
     884     1572086 :   l += n;
     885     1572086 :   if (n < 0)
     886             :   {
     887         854 :     if (l <= 2) return pol_0(varn(a));
     888         784 :     b = cgetg(l, t_POL); b[1] = a[1];
     889         784 :     a -= n;
     890        2982 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     891             :   } else {
     892     1571232 :     b = cgetg(l, t_POL); b[1] = a[1];
     893     1571232 :     a -= n; n += 2;
     894     3935011 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     895     4375630 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     896             :   }
     897     1572016 :   return b;
     898             : }
     899             : 
     900             : GEN
     901      316113 : RgX_rotate_shallow(GEN P, long k, long p)
     902             : {
     903      316113 :   long i, l = lgpol(P);
     904             :   GEN r;
     905      316113 :   if (signe(P)==0)
     906        1365 :     return pol_0(varn(P));
     907      314748 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     908     2090452 :   for(i=0; i<p; i++)
     909             :   {
     910     1775704 :     long s = 2+(i+k)%p;
     911     1775704 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
     912             :   }
     913      314748 :   return RgX_renormalize(r);
     914             : }
     915             : 
     916             : GEN
     917     2977754 : RgX_mulXn(GEN x, long d)
     918             : {
     919             :   pari_sp av;
     920             :   GEN z;
     921             :   long v;
     922     2977754 :   if (d >= 0) return RgX_shift(x, d);
     923     1466215 :   d = -d;
     924     1466215 :   v = RgX_val(x);
     925     1466215 :   if (v >= d) return RgX_shift(x, -d);
     926     1466201 :   av = avma;
     927     1466201 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
     928     1466201 :   return gerepileupto(av, z);
     929             : }
     930             : 
     931             : long
     932         588 : RgXV_maxdegree(GEN x)
     933             : {
     934         588 :   long d = -1, i, l = lg(x);
     935        4494 :   for (i = 1; i < l; i++)
     936        3906 :     d = maxss(d, degpol(gel(x,i)));
     937         588 :   return d;
     938             : }
     939             : 
     940             : long
     941     2903806 : RgX_val(GEN x)
     942             : {
     943     2903806 :   long i, lx = lg(x);
     944     2903806 :   if (lx == 2) return LONG_MAX;
     945     3828149 :   for (i = 2; i < lx; i++)
     946     3828107 :     if (!isexactzero(gel(x,i))) break;
     947     2903666 :   if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
     948     2903624 :   return i - 2;
     949             : }
     950             : long
     951    68792121 : RgX_valrem(GEN x, GEN *Z)
     952             : {
     953    68792121 :   long v, i, lx = lg(x);
     954    68792121 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
     955   113131001 :   for (i = 2; i < lx; i++)
     956   113137276 :     if (!isexactzero(gel(x,i))) break;
     957             :   /* possible with nonrational zeros */
     958    68811890 :   if (i == lx) { *Z = pol_0(varn(x)); return LONG_MAX; }
     959    68811890 :   v = i - 2;
     960    68811890 :   *Z = RgX_shift_shallow(x, -v);
     961    68824897 :   return v;
     962             : }
     963             : long
     964       34258 : RgX_valrem_inexact(GEN x, GEN *Z)
     965             : {
     966             :   long v;
     967       34258 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
     968       34251 :   for (v = 0;; v++)
     969       35506 :     if (!gequal0(gel(x,2+v))) break;
     970       34251 :   if (Z) *Z = RgX_shift_shallow(x, -v);
     971       34251 :   return v;
     972             : }
     973             : 
     974             : GEN
     975       63651 : RgXQC_red(GEN x, GEN T)
     976      392315 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
     977             : 
     978             : GEN
     979        1204 : RgXQV_red(GEN x, GEN T)
     980       27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
     981             : 
     982             : GEN
     983       12201 : RgXQM_red(GEN x, GEN T)
     984       75852 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
     985             : 
     986             : GEN
     987         322 : RgXQM_mul(GEN P, GEN Q, GEN T)
     988             : {
     989         322 :   return RgXQM_red(RgM_mul(P, Q), T);
     990             : }
     991             : 
     992             : GEN
     993      443113 : RgXQX_red(GEN P, GEN T)
     994             : {
     995      443113 :   long i, l = lg(P);
     996      443113 :   GEN Q = cgetg(l, t_POL);
     997      443113 :   Q[1] = P[1];
     998     2382601 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
     999      443113 :   return normalizepol_lg(Q, l);
    1000             : }
    1001             : 
    1002             : GEN
    1003      631875 : RgX_deriv(GEN x)
    1004             : {
    1005      631875 :   long i,lx = lg(x)-1;
    1006             :   GEN y;
    1007             : 
    1008      631875 :   if (lx<3) return pol_0(varn(x));
    1009      630174 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
    1010     3135747 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
    1011      630163 :   y[1] = x[1]; return normalizepol_lg(y,i);
    1012             : }
    1013             : 
    1014             : GEN
    1015     1011646 : RgX_recipspec_shallow(GEN x, long l, long n)
    1016             : {
    1017             :   long i;
    1018     1011646 :   GEN z = cgetg(n+2,t_POL);
    1019     1011645 :   z[1] = 0; z += 2;
    1020    18034345 :   for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
    1021     1139740 :   for(   ; i<n; i++) gel(z, n-i-1) = gen_0;
    1022     1011645 :   return normalizepol_lg(z-2,n+2);
    1023             : }
    1024             : 
    1025             : GEN
    1026      532851 : RgXn_recip_shallow(GEN P, long n)
    1027             : {
    1028      532851 :   GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
    1029      532857 :   setvarn(Q, varn(P));
    1030      532857 :   return Q;
    1031             : }
    1032             : 
    1033             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1034             : GEN
    1035       13531 : RgX_recip(GEN x)
    1036             : {
    1037             :   long lx, i, j;
    1038       13531 :   GEN y = cgetg_copy(x, &lx);
    1039      119441 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
    1040       13531 :   return normalizepol_lg(y,lx);
    1041             : }
    1042             : /* shallow version */
    1043             : GEN
    1044          70 : RgX_recip_shallow(GEN x)
    1045             : {
    1046             :   long lx, i, j;
    1047          70 :   GEN y = cgetg_copy(x, &lx);
    1048         259 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1049          70 :   return normalizepol_lg(y,lx);
    1050             : }
    1051             : 
    1052             : GEN
    1053     2702188 : RgX_recip_i(GEN x)
    1054             : {
    1055             :   long lx, i, j;
    1056     2702188 :   GEN y = cgetg_copy(x, &lx);
    1057    15862237 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
    1058     2702193 :   return y;
    1059             : }
    1060             : /*******************************************************************/
    1061             : /*                                                                 */
    1062             : /*                      ADDITION / SUBTRACTION                     */
    1063             : /*                                                                 */
    1064             : /*******************************************************************/
    1065             : /* same variable */
    1066             : GEN
    1067    78158112 : RgX_add(GEN x, GEN y)
    1068             : {
    1069    78158112 :   long i, lx = lg(x), ly = lg(y);
    1070             :   GEN z;
    1071    78158112 :   if (ly <= lx) {
    1072    70196314 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1073   285147657 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1074   147610288 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1075    70143722 :     z = normalizepol_lg(z, lx);
    1076             :   } else {
    1077     7961798 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1078    40141469 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1079    27502703 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1080     7962922 :     z = normalizepol_lg(z, ly);
    1081             :   }
    1082    78155913 :   return z;
    1083             : }
    1084             : GEN
    1085    57961632 : RgX_sub(GEN x, GEN y)
    1086             : {
    1087    57961632 :   long i, lx = lg(x), ly = lg(y);
    1088             :   GEN z;
    1089    57961632 :   if (ly <= lx) {
    1090    36204711 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1091   197467707 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1092    79653167 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1093    36171849 :     z = normalizepol_lg(z, lx);
    1094             :   } else {
    1095    21756921 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1096    89329461 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1097    50567607 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1098    21713354 :     z = normalizepol_lg(z, ly);
    1099             :   }
    1100    57957247 :   return z;
    1101             : }
    1102             : GEN
    1103     5094764 : RgX_neg(GEN x)
    1104             : {
    1105     5094764 :   long i, lx = lg(x);
    1106     5094764 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    1107    41848637 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    1108     5094531 :   return y;
    1109             : }
    1110             : 
    1111             : GEN
    1112    21338011 : RgX_Rg_add(GEN y, GEN x)
    1113             : {
    1114             :   GEN z;
    1115    21338011 :   long lz = lg(y), i;
    1116    21338011 :   if (lz == 2) return scalarpol(x,varn(y));
    1117    18330536 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1118    18330513 :   gel(z,2) = gadd(gel(y,2),x);
    1119    83540656 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1120             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1121             :    * an inexact 0 */
    1122    18330630 :   return normalizepol_lg(z,lz);
    1123             : }
    1124             : GEN
    1125       64685 : RgX_Rg_add_shallow(GEN y, GEN x)
    1126             : {
    1127             :   GEN z;
    1128       64685 :   long lz = lg(y), i;
    1129       64685 :   if (lz == 2) return scalarpol(x,varn(y));
    1130       64685 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1131       64685 :   gel(z,2) = gadd(gel(y,2),x);
    1132      129489 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1133       64683 :   return z = normalizepol_lg(z,lz);
    1134             : }
    1135             : GEN
    1136      175905 : RgX_Rg_sub(GEN y, GEN x)
    1137             : {
    1138             :   GEN z;
    1139      175905 :   long lz = lg(y), i;
    1140      175905 :   if (lz == 2)
    1141             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1142        2548 :     long v = varn(y);
    1143        2548 :     if (isrationalzero(x)) return pol_0(v);
    1144          14 :     z = cgetg(3,t_POL);
    1145          14 :     z[1] = gequal0(x)? evalvarn(v)
    1146          14 :                    : evalvarn(v) | evalsigne(1);
    1147          14 :     gel(z,2) = gneg(x); return z;
    1148             :   }
    1149      173357 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1150      173356 :   gel(z,2) = gsub(gel(y,2),x);
    1151      453576 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1152      173357 :   return z = normalizepol_lg(z,lz);
    1153             : }
    1154             : GEN
    1155     3683338 : Rg_RgX_sub(GEN x, GEN y)
    1156             : {
    1157             :   GEN z;
    1158     3683338 :   long lz = lg(y), i;
    1159     3683338 :   if (lz == 2) return scalarpol(x,varn(y));
    1160     3682260 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1161     3682243 :   gel(z,2) = gsub(x, gel(y,2));
    1162     6732022 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1163     3682308 :   return z = normalizepol_lg(z,lz);
    1164             : }
    1165             : /*******************************************************************/
    1166             : /*                                                                 */
    1167             : /*                  KARATSUBA MULTIPLICATION                       */
    1168             : /*                                                                 */
    1169             : /*******************************************************************/
    1170             : #if 0
    1171             : /* to debug Karatsuba-like routines */
    1172             : GEN
    1173             : zx_debug_spec(GEN x, long nx)
    1174             : {
    1175             :   GEN z = cgetg(nx+2,t_POL);
    1176             :   long i;
    1177             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1178             :   z[1] = evalsigne(1); return z;
    1179             : }
    1180             : 
    1181             : GEN
    1182             : RgX_debug_spec(GEN x, long nx)
    1183             : {
    1184             :   GEN z = cgetg(nx+2,t_POL);
    1185             :   long i;
    1186             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1187             :   z[1] = evalsigne(1); return z;
    1188             : }
    1189             : #endif
    1190             : 
    1191             : /* generic multiplication */
    1192             : GEN
    1193     7831548 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1194             : {
    1195             :   GEN z, t;
    1196             :   long i;
    1197     7831548 :   if (nx == ny) {
    1198     1081254 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1199     3883920 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1200     1081248 :     return normalizepol_lg(z, nx+2);
    1201             :   }
    1202     6750294 :   if (ny < nx) {
    1203     6556210 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1204    26287518 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1205    15988093 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1206     6555561 :     return normalizepol_lg(z, nx+2);
    1207             :   } else {
    1208      194084 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1209     3708661 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1210      477764 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1211      194114 :     return normalizepol_lg(z, ny+2);
    1212             :   }
    1213             : }
    1214             : GEN
    1215      236988 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1216             : {
    1217             :   GEN z, t;
    1218             :   long i;
    1219      236988 :   if (nx == ny) {
    1220       12649 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1221     2143645 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1222       12649 :     return normalizepol_lg(z, nx+2);
    1223             :   }
    1224      224339 :   if (ny < nx) {
    1225      222638 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1226     3770543 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1227     2566055 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1228      222638 :     return normalizepol_lg(z, nx+2);
    1229             :   } else {
    1230        1701 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1231      327040 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1232       12040 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1233        1701 :     return normalizepol_lg(z, ny+2);
    1234             :   }
    1235             : }
    1236             : 
    1237             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1238             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1239             :  * t_VECSMALL, to hold this, which can be left on stack: gerepile
    1240             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1241             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1242             : static GEN
    1243    43581592 : RgXspec_kill0(GEN x, long lx)
    1244             : {
    1245    43581592 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
    1246             :   long i;
    1247   160376298 :   for (i=0; i <lx; i++)
    1248             :   {
    1249   116794768 :     GEN c = gel(x,i);
    1250   116794768 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1251             :   }
    1252    43581530 :   return z;
    1253             : }
    1254             : 
    1255             : INLINE GEN
    1256    83410117 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1257             : {
    1258    83410117 :   pari_sp av = avma;
    1259    83410117 :   GEN s = NULL;
    1260             :   long i;
    1261             : 
    1262   341696674 :   for (i=a; i<b; i++)
    1263   258299229 :     if (gel(y,i) && gel(x,-i))
    1264             :     {
    1265   175608068 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1266   175600576 :       s = s? gadd(s, t): t;
    1267             :     }
    1268    83397445 :   return s? gerepileupto(av, s): gen_0;
    1269             : }
    1270             : 
    1271             : /* assume nx >= ny > 0, return x * y * t^v */
    1272             : static GEN
    1273    16925751 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1274             : {
    1275             :   long i, lz, nz;
    1276             :   GEN z;
    1277             : 
    1278    16925751 :   x = RgXspec_kill0(x,nx);
    1279    16925734 :   y = RgXspec_kill0(y,ny);
    1280    16925730 :   lz = nx + ny + 1; nz = lz-2;
    1281    16925730 :   lz += v;
    1282    16925730 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1283    28653730 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1284    48887251 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1285    34275810 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1286    31961488 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1287    16925114 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1288             : }
    1289             : 
    1290             : /* return (x * X^d) + y. Assume d > 0 */
    1291             : GEN
    1292     7962211 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1293             : {
    1294             :   GEN x, y, xd, yd, zd;
    1295             :   long a, lz, nx, ny;
    1296             : 
    1297     7962211 :   if (!signe(x0)) return y0;
    1298     7862045 :   ny = lgpol(y0);
    1299     7862044 :   nx = lgpol(x0);
    1300     7862164 :   zd = (GEN)avma;
    1301     7862164 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1302     7862164 :   if (a <= 0)
    1303             :   {
    1304     1017079 :     lz = nx+d+2;
    1305     1017079 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1306     2472466 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1307     1017082 :     x = zd + a;
    1308     1029771 :     while (zd > x) gel(--zd,0) = gen_0;
    1309             :   }
    1310             :   else
    1311             :   {
    1312     6845085 :     xd = new_chunk(d); yd = y+d;
    1313     6845075 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1314     6845038 :     lz = (a>nx)? ny+2: lg(x)+d;
    1315    37680180 :     x += 2; while (xd > x) *--zd = *--xd;
    1316             :   }
    1317    20031309 :   while (yd > y) *--zd = *--yd;
    1318     7862120 :   *--zd = x0[1];
    1319     7862120 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1320             : }
    1321             : GEN
    1322      529526 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1323             : {
    1324             :   GEN x, y, xd, yd, zd;
    1325             :   long a, lz, nx, ny;
    1326             : 
    1327      529526 :   if (!signe(x0)) return RgX_copy(y0);
    1328      528756 :   nx = lgpol(x0);
    1329      528756 :   ny = lgpol(y0);
    1330      528756 :   zd = (GEN)avma;
    1331      528756 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1332      528756 :   if (a <= 0)
    1333             :   {
    1334      291768 :     lz = nx+d+2;
    1335      291768 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1336     4287579 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1337      291768 :     x = zd + a;
    1338      757913 :     while (zd > x) gel(--zd,0) = gen_0;
    1339             :   }
    1340             :   else
    1341             :   {
    1342      236988 :     xd = new_chunk(d); yd = y+d;
    1343      236988 :     x = RgX_addspec(x,yd, nx,a);
    1344      236988 :     lz = (a>nx)? ny+2: lg(x)+d;
    1345     8594984 :     x += 2; while (xd > x) *--zd = *--xd;
    1346             :   }
    1347     2721800 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1348      528756 :   *--zd = x0[1];
    1349      528756 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1350             : }
    1351             : 
    1352             : /* return x * y mod t^n */
    1353             : static GEN
    1354     4856477 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1355             : {
    1356     4856477 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1357             :   GEN z;
    1358     4856477 :   if (lx < 0) return pol_0(varn(x));
    1359     4856477 :   if (ly < 0) return pol_0(varn(x));
    1360     4856477 :   z = cgetg(lz, t_POL) + 2;
    1361     4856477 :   x+=2; if (lx > n) lx = n;
    1362     4856477 :   y+=2; if (ly > n) ly = n;
    1363     4856477 :   z[-1] = x[-1];
    1364     4856477 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1365     4856477 :   x = RgXspec_kill0(x, lx);
    1366     4856477 :   y = RgXspec_kill0(y, ly);
    1367             :   /* x:y:z [i] = term of degree i */
    1368    21316397 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1369     7408831 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1370     4906715 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1371     4856477 :   return normalizepol_lg(z - 2, lz);
    1372             : }
    1373             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1374             : static GEN
    1375     6029030 : RgXn_mul2(GEN f, GEN g, long n)
    1376             : {
    1377     6029030 :   pari_sp av = avma;
    1378             :   GEN fe,fo, ge,go, l,h,m;
    1379             :   long n0, n1;
    1380     6029030 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1381     4891407 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1382       34930 :   n0 = n>>1; n1 = n-n0;
    1383       34930 :   RgX_even_odd(f, &fe, &fo);
    1384       34930 :   RgX_even_odd(g, &ge, &go);
    1385       34930 :   l = RgXn_mul2(fe,ge,n1);
    1386       34930 :   h = RgXn_mul2(fo,go,n0);
    1387       34930 :   m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1388             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1389             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1390       34930 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1391             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1392       34930 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1393       34930 :   m = RgX_inflate(m,2);
    1394             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1395       34930 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1396       34930 :   h = RgX_inflate(h,2);
    1397       34930 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1398       34930 :   return gerepileupto(av, h);
    1399             : }
    1400             : /* (f*g) \/ x^n */
    1401             : static GEN
    1402     1047569 : RgX_mulhigh_i2(GEN f, GEN g, long n)
    1403             : {
    1404     1047569 :   long d = degpol(f)+degpol(g) + 1 - n;
    1405             :   GEN h;
    1406     1047569 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1407       21310 :   h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
    1408             :                                   RgX_recip_i(g), d));
    1409       21310 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1410             : }
    1411             : 
    1412             : /* (f*g) \/ x^n */
    1413             : static GEN
    1414           0 : RgX_sqrhigh_i2(GEN f, long n)
    1415             : {
    1416           0 :   long d = 2*degpol(f)+ 1 - n;
    1417             :   GEN h;
    1418           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1419           0 :   h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
    1420           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1421             : }
    1422             : 
    1423             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1424             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1425             :  * Only c, c0, c1, c2 are genuine GEN.
    1426             :  */
    1427             : GEN
    1428    17663939 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1429             : {
    1430             :   GEN a0, c, c0;
    1431    17663939 :   long n0, n0a, i, v = 0;
    1432             :   pari_sp av;
    1433             : 
    1434    24045839 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1435    23544753 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1436    17663891 :   if (na < nb) swapspec(a,b, na,nb);
    1437    17663891 :   if (!nb) return pol_0(0);
    1438             : 
    1439    17420294 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1440      494537 :   RgX_shift_inplace_init(v);
    1441      494557 :   i = (na>>1); n0 = na-i; na = i;
    1442      494557 :   av = avma; a0 = a+n0; n0a = n0;
    1443     1350057 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1444             : 
    1445      494557 :   if (nb > n0)
    1446             :   {
    1447             :     GEN b0,c1,c2;
    1448             :     long n0b;
    1449             : 
    1450      493235 :     nb -= n0; b0 = b+n0; n0b = n0;
    1451     1461198 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1452      493235 :     c = RgX_mulspec(a,b,n0a,n0b);
    1453      493235 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1454             : 
    1455      493235 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1456      493235 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1457             : 
    1458      493235 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1459      493235 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1460      493235 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1461             :   }
    1462             :   else
    1463             :   {
    1464        1322 :     c = RgX_mulspec(a,b,n0a,nb);
    1465        1322 :     c0 = RgX_mulspec(a0,b,na,nb);
    1466             :   }
    1467      494557 :   c0 = RgX_addmulXn(c0,c,n0);
    1468      494557 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1469             : }
    1470             : 
    1471             : INLINE GEN
    1472       78025 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1473             : {
    1474       78025 :   pari_sp av = avma;
    1475       78025 :   GEN s = NULL;
    1476       78025 :   long j, l = (i+1)>>1;
    1477      176330 :   for (j=a; j<l; j++)
    1478             :   {
    1479       98306 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1480       98306 :     if (xj && xx)
    1481             :     {
    1482       91040 :       GEN t = gmul(xj, xx);
    1483       91039 :       s = s? gadd(s, t): t;
    1484             :     }
    1485             :   }
    1486       78024 :   if (s) s = gshift(s,1);
    1487       78026 :   if ((i&1) == 0)
    1488             :   {
    1489       47695 :     GEN t = gel(x, i>>1);
    1490       47695 :     if (t) {
    1491       45399 :       t = gsqr(t);
    1492       45392 :       s = s? gadd(s, t): t;
    1493             :     }
    1494             :   }
    1495       78019 :   return s? gerepileupto(av,s): gen_0;
    1496             : }
    1497             : static GEN
    1498       17359 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1499             : {
    1500             :   long i, lz, nz;
    1501             :   GEN z;
    1502             : 
    1503       17359 :   if (!nx) return pol_0(0);
    1504       17359 :   x = RgXspec_kill0(x,nx);
    1505       17362 :   lz = (nx << 1) + 1, nz = lz-2;
    1506       17362 :   lz += v;
    1507       17362 :   z = cgetg(lz,t_POL) + 2;
    1508       21030 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1509       65052 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1510       47693 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1511       17360 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1512             : }
    1513             : /* return x^2 mod t^n */
    1514             : static GEN
    1515           0 : RgXn_sqr_basecase(GEN x, long n)
    1516             : {
    1517           0 :   long i, lz = n+2, lx = lgpol(x);
    1518             :   GEN z;
    1519           0 :   if (lx < 0) return pol_0(varn(x));
    1520           0 :   z = cgetg(lz, t_POL);
    1521           0 :   z[1] = x[1];
    1522           0 :   x+=2; if (lx > n) lx = n;
    1523           0 :   x = RgXspec_kill0(x,lx);
    1524           0 :   z+=2;/* x:z [i] = term of degree i */
    1525           0 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1526           0 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1527           0 :   z -= 2; return normalizepol_lg(z, lz);
    1528             : }
    1529             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1530             : static GEN
    1531         315 : RgXn_sqr2(GEN f, long n)
    1532             : {
    1533         315 :   pari_sp av = avma;
    1534             :   GEN fe,fo, l,h,m;
    1535             :   long n0, n1;
    1536         315 :   if (2*degpol(f) < n) return RgX_sqr_i(f);
    1537           0 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1538           0 :   n0 = n>>1; n1 = n-n0;
    1539           0 :   RgX_even_odd(f, &fe, &fo);
    1540           0 :   l = RgXn_sqr(fe,n1);
    1541           0 :   h = RgXn_sqr(fo,n0);
    1542           0 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1543             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1544             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1545           0 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1546             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1547           0 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1548           0 :   m = RgX_inflate(m,2);
    1549             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1550           0 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1551           0 :   h = RgX_inflate(h,2);
    1552           0 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1553           0 :   return gerepileupto(av, h);
    1554             : }
    1555             : GEN
    1556       17398 : RgX_sqrspec(GEN a, long na)
    1557             : {
    1558             :   GEN a0, c, c0, c1;
    1559       17398 :   long n0, n0a, i, v = 0;
    1560             :   pari_sp av;
    1561             : 
    1562       19232 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1563       17398 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1564          39 :   RgX_shift_inplace_init(v);
    1565          39 :   i = (na>>1); n0 = na-i; na = i;
    1566          39 :   av = avma; a0 = a+n0; n0a = n0;
    1567          39 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1568             : 
    1569          39 :   c = RgX_sqrspec(a,n0a);
    1570          39 :   c0 = RgX_sqrspec(a0,na);
    1571          39 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1572          39 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1573          39 :   c0 = RgX_addmulXn(c0,c,n0);
    1574          39 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1575             : }
    1576             : 
    1577             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1578             : GEN
    1579     1456976 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1580             : {
    1581     1456976 :   GEN z = RgX_mul(A, B);
    1582     1456965 :   if (a < b)
    1583       10732 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1584     1446233 :   else if (a > b)
    1585      830421 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1586             :   else
    1587      615812 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1588     1456963 :   return z;
    1589             : }
    1590             : 
    1591             : GEN
    1592    16181650 : RgX_mul_i(GEN x, GEN y)
    1593             : {
    1594    16181650 :   GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
    1595    16181110 :   setvarn(z, varn(x)); return z;
    1596             : }
    1597             : 
    1598             : GEN
    1599       17319 : RgX_sqr_i(GEN x)
    1600             : {
    1601       17319 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1602       17320 :   setvarn(z,varn(x)); return z;
    1603             : }
    1604             : 
    1605             : /*******************************************************************/
    1606             : /*                                                                 */
    1607             : /*                               DIVISION                          */
    1608             : /*                                                                 */
    1609             : /*******************************************************************/
    1610             : GEN
    1611     4020258 : RgX_Rg_divexact(GEN x, GEN y) {
    1612     4020258 :   long i, lx = lg(x);
    1613             :   GEN z;
    1614     4020258 :   if (lx == 2) return gcopy(x);
    1615     3970657 :   switch(typ(y))
    1616             :   {
    1617     3941637 :     case t_INT:
    1618     3941637 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1619     3778105 :       break;
    1620        3969 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1621             :   }
    1622     3803156 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1623    24432791 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1624     3791077 :   return z;
    1625             : }
    1626             : GEN
    1627    27900797 : RgX_Rg_div(GEN x, GEN y) {
    1628    27900797 :   long i, lx = lg(x);
    1629             :   GEN z;
    1630    27900797 :   if (lx == 2) return gcopy(x);
    1631    27675797 :   switch(typ(y))
    1632             :   {
    1633    18515973 :     case t_INT:
    1634    18515973 :       if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1635     3721976 :       break;
    1636        4193 :     case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
    1637             :   }
    1638    12877607 :   z = cgetg(lx, t_POL); z[1] = x[1];
    1639    49508047 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1640    12877268 :   return normalizepol_lg(z, lx);
    1641             : }
    1642             : GEN
    1643       13510 : RgX_normalize(GEN x)
    1644             : {
    1645       13510 :   GEN z, d = NULL;
    1646       13510 :   long i, n = lg(x)-1;
    1647       13510 :   for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
    1648       13510 :   if (i == 1) return pol_0(varn(x));
    1649       13510 :   if (i == n && isint1(d)) return x;
    1650        5236 :   n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
    1651       12509 :   for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
    1652        5236 :   gel(z,n) = Rg_get_1(d); return z;
    1653             : }
    1654             : GEN
    1655        2625 : RgX_divs(GEN x, long y) {
    1656             :   long i, lx;
    1657        2625 :   GEN z = cgetg_copy(x, &lx); z[1] = x[1];
    1658       10318 :   for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),y);
    1659        2625 :   return normalizepol_lg(z, lx);
    1660             : }
    1661             : GEN
    1662      230284 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1663             : {
    1664      230284 :   long l = lg(a), i;
    1665             :   GEN a0, z0, z;
    1666             : 
    1667      230284 :   if (l <= 3)
    1668             :   {
    1669           0 :     if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
    1670           0 :     return pol_0(0);
    1671             :   }
    1672      230284 :   z = cgetg(l-1, t_POL);
    1673      230292 :   z[1] = a[1];
    1674      230292 :   a0 = a + l-1;
    1675      230292 :   z0 = z + l-2; *z0 = *a0--;
    1676     2841646 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1677             :   {
    1678     2611362 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1679     2611354 :     gel(z0,0) = t;
    1680             :   }
    1681      230284 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1682      230284 :   return z;
    1683             : }
    1684             : /* Polynomial division x / y:
    1685             :  *   if pr = ONLY_REM return remainder, otherwise return quotient
    1686             :  *   if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
    1687             :  *   if pr != NULL set *pr to remainder, as the last object on stack */
    1688             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1689             : static GEN
    1690    20174795 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
    1691             : {
    1692             :   pari_sp avy, av, av1;
    1693             :   long dx,dy,dz,i,j,sx,lr;
    1694             :   GEN z,p1,p2,rem,y_lead,mod,p;
    1695             :   GEN (*f)(GEN,GEN);
    1696             : 
    1697    20174795 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1698             : 
    1699    20174795 :   dy = degpol(y);
    1700    20174742 :   y_lead = gel(y,dy+2);
    1701    20174742 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1702             :   {
    1703           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1704           0 :     for (dy--; dy>=0; dy--)
    1705             :     {
    1706           0 :       y_lead = gel(y,dy+2);
    1707           0 :       if (!gequal0(y_lead)) break;
    1708             :     }
    1709             :   }
    1710    20174706 :   if (!dy) /* y is constant */
    1711             :   {
    1712        3051 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1713        3051 :     z = RgX_Rg_div(x, y_lead);
    1714        3051 :     if (pr == ONLY_DIVIDES) return z;
    1715        1168 :     if (pr) *pr = pol_0(varn(x));
    1716        1168 :     return z;
    1717             :   }
    1718    20171655 :   dx = degpol(x);
    1719    20171592 :   if (dx < dy)
    1720             :   {
    1721     2288084 :     if (pr == ONLY_REM) return RgX_copy(x);
    1722      373861 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1723      373840 :     z = pol_0(varn(x));
    1724      373840 :     if (pr) *pr = RgX_copy(x);
    1725      373840 :     return z;
    1726             :   }
    1727             : 
    1728             :   /* x,y in R[X], y non constant */
    1729    17883508 :   av = avma;
    1730    17883508 :   p = NULL;
    1731    17883508 :   if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1732             :   {
    1733      225841 :     z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
    1734      225841 :     if (!z) return gc_NULL(av);
    1735      225841 :     z = FpX_to_mod(z, p);
    1736      225841 :     if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
    1737      121681 :       return gerepileupto(av, z);
    1738      104160 :     *pr = FpX_to_mod(*pr, p);
    1739      104160 :     gerepileall(av, 2, pr, &z);
    1740      104160 :     return z;
    1741             :   }
    1742    17657820 :   switch(typ(y_lead))
    1743             :   {
    1744           0 :     case t_REAL:
    1745           0 :       y_lead = ginv(y_lead);
    1746           0 :       f = gmul; mod = NULL;
    1747           0 :       break;
    1748        1954 :     case t_INTMOD:
    1749        1954 :     case t_POLMOD: y_lead = ginv(y_lead);
    1750        1954 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1751        1954 :       break;
    1752    17655866 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1753    17655847 :       f = gdiv; mod = NULL;
    1754             :   }
    1755             : 
    1756    17657801 :   if (y_lead == NULL)
    1757    15841833 :     p2 = gel(x,dx+2);
    1758             :   else {
    1759             :     for(;;) {
    1760     1815968 :       p2 = f(gel(x,dx+2),y_lead);
    1761     1815969 :       p2 = simplify_shallow(p2);
    1762     1815969 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1763             :     }
    1764     1815969 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1765             :     {
    1766           0 :       if (pr == ONLY_DIVIDES) {
    1767           0 :         set_avma(av);
    1768           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1769             :       }
    1770           0 :       if (pr == ONLY_REM)
    1771             :       {
    1772           0 :         if (dx < 0)
    1773           0 :           return gerepilecopy(av, scalarpol(p2, varn(x)));
    1774             :         else
    1775             :         {
    1776             :           GEN t;
    1777           0 :           set_avma(av);
    1778           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1779           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1780           0 :           return t;
    1781             :         }
    1782             :       }
    1783           0 :       if (pr) /* cf ONLY_REM above */
    1784             :       {
    1785           0 :         if (dx < 0)
    1786             :         {
    1787           0 :           p2 = gclone(p2);
    1788           0 :           set_avma(av);
    1789           0 :           z = pol_0(varn(x));
    1790           0 :           x = scalarpol(p2, varn(x));
    1791           0 :           gunclone(p2);
    1792             :         }
    1793             :         else
    1794             :         {
    1795             :           GEN t;
    1796           0 :           set_avma(av);
    1797           0 :           z = pol_0(varn(x));
    1798           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1799           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1800           0 :           x = t;
    1801             :         }
    1802           0 :         *pr = x;
    1803             :       }
    1804             :       else
    1805             :       {
    1806           0 :         set_avma(av);
    1807           0 :         z = pol_0(varn(x));
    1808             :       }
    1809           0 :       return z;
    1810             :     }
    1811             :   }
    1812             :   /* dx >= dy */
    1813    17657802 :   avy = avma;
    1814    17657802 :   dz = dx-dy;
    1815    17657802 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1816    17657767 :   x += 2;
    1817    17657767 :   z += 2;
    1818    17657767 :   y += 2;
    1819    17657767 :   gel(z,dz) = gcopy(p2);
    1820             : 
    1821    49702392 :   for (i=dx-1; i>=dy; i--)
    1822             :   {
    1823    32044691 :     av1=avma; p1=gel(x,i);
    1824  1101733912 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1825    32038466 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1826             : 
    1827    32038464 :     if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
    1828             :     else
    1829    22483221 :       p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
    1830    32044076 :     gel(z,i-dy) = p1;
    1831             :   }
    1832    17657701 :   if (!pr) return gerepileupto(av,z-2);
    1833             : 
    1834    10947109 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1835    10947069 :   for (sx=0; ; i--)
    1836             :   {
    1837      349430 :     p1 = gel(x,i);
    1838             :     /* we always enter this loop at least once */
    1839    30294318 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1840    11295176 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1841    11295176 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
    1842     3500212 :     if (!isexactzero(p1)) break;
    1843     3405570 :     if (!i) break;
    1844      349427 :     set_avma(av1);
    1845             :   }
    1846    10946197 :   if (pr == ONLY_DIVIDES)
    1847             :   {
    1848        1687 :     if (sx) return gc_NULL(av);
    1849        1680 :     set_avma((pari_sp)rem); return gerepileupto(av,z-2);
    1850             :   }
    1851    10944510 :   lr=i+3; rem -= lr;
    1852    10944510 :   if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
    1853    10900912 :   else p1 = gerepileupto((pari_sp)rem,p1);
    1854    10945385 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1855    10945360 :   rem[1] = z[-1];
    1856    10945360 :   rem += 2;
    1857    10945360 :   gel(rem,i) = p1;
    1858    17137762 :   for (i--; i>=0; i--)
    1859             :   {
    1860     6192433 :     av1=avma; p1 = gel(x,i);
    1861    17080244 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1862     6191790 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1863     6192059 :     gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
    1864             :   }
    1865    10945329 :   rem -= 2;
    1866    10945329 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1867    10945552 :   if (pr == ONLY_REM) return gerepileupto(av,rem);
    1868     5924180 :   z -= 2;
    1869             :   {
    1870     5924180 :     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
    1871     5924180 :     gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
    1872             :   }
    1873             : }
    1874             : 
    1875             : GEN
    1876    13239253 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1877             : {
    1878    13239253 :   if (pr == ONLY_REM) return RgX_rem(x, y);
    1879    13239253 :   return RgX_divrem_i(x, y, pr);
    1880             : }
    1881             : 
    1882             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1883             : GEN
    1884      132778 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1885             : {
    1886             :   long vx, dx, dy, dz, i, j, sx, lr;
    1887             :   pari_sp av0, av, tetpil;
    1888             :   GEN z,p1,rem,lead;
    1889             : 
    1890      132778 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1891      132778 :   vx = varn(x);
    1892      132778 :   dx = degpol(x);
    1893      132778 :   dy = degpol(y);
    1894      132778 :   if (dx < dy)
    1895             :   {
    1896       25955 :     if (pr)
    1897             :     {
    1898       25955 :       av0 = avma; x = RgXQX_red(x, T);
    1899       25955 :       if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
    1900       25955 :       if (pr == ONLY_REM) return x;
    1901           0 :       *pr = x;
    1902             :     }
    1903           0 :     return pol_0(vx);
    1904             :   }
    1905      106823 :   lead = leading_coeff(y);
    1906      106823 :   if (!dy) /* y is constant */
    1907             :   {
    1908         546 :     if (pr && pr != ONLY_DIVIDES)
    1909             :     {
    1910           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1911           0 :       *pr = pol_0(vx);
    1912             :     }
    1913         546 :     if (gequal1(lead)) return RgX_copy(x);
    1914           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1915           0 :     return gerepile(av0,tetpil,RgXQX_red(x,T));
    1916             :   }
    1917      106277 :   av0 = avma; dz = dx-dy;
    1918      106277 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1919      106277 :   set_avma(av0);
    1920      106277 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1921      106277 :   x += 2; y += 2; z += 2;
    1922             : 
    1923      106277 :   p1 = gel(x,dx); av = avma;
    1924      106277 :   gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    1925      532381 :   for (i=dx-1; i>=dy; i--)
    1926             :   {
    1927      426104 :     av=avma; p1=gel(x,i);
    1928     2230596 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1929      426079 :     if (lead) p1 = gmul(grem(p1, T), lead);
    1930      426076 :     tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
    1931             :   }
    1932      106277 :   if (!pr) { guncloneNULL(lead); return z-2; }
    1933             : 
    1934      105031 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    1935      105031 :   for (sx=0; ; i--)
    1936             :   {
    1937       70675 :     p1 = gel(x,i);
    1938      603116 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1939      175704 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    1940       97827 :     if (!i) break;
    1941       70675 :     set_avma(av);
    1942             :   }
    1943      105030 :   if (pr == ONLY_DIVIDES)
    1944             :   {
    1945       26305 :     guncloneNULL(lead);
    1946       26306 :     if (sx) return gc_NULL(av0);
    1947       23793 :     return gc_const((pari_sp)rem, z-2);
    1948             :   }
    1949       78725 :   lr=i+3; rem -= lr;
    1950       78725 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1951       78725 :   rem[1] = z[-1];
    1952       78725 :   p1 = gerepile((pari_sp)rem,tetpil,p1);
    1953       78725 :   rem += 2; gel(rem,i) = p1;
    1954      167838 :   for (i--; i>=0; i--)
    1955             :   {
    1956       89113 :     av=avma; p1 = gel(x,i);
    1957      303249 :     for (j=0; j<=i && j<=dz; j++)
    1958      214136 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1959       89113 :     tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
    1960             :   }
    1961       78725 :   rem -= 2;
    1962       78725 :   guncloneNULL(lead);
    1963       78725 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1964       78725 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
    1965         154 :   *pr = rem; return z-2;
    1966             : }
    1967             : 
    1968             : /*******************************************************************/
    1969             : /*                                                                 */
    1970             : /*                        PSEUDO-DIVISION                          */
    1971             : /*                                                                 */
    1972             : /*******************************************************************/
    1973             : INLINE GEN
    1974     1314594 : rem(GEN c, GEN T)
    1975             : {
    1976     1314594 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    1977     1314596 :   return c;
    1978             : }
    1979             : 
    1980             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    1981             : int
    1982       15188 : ZXQX_dvd(GEN x, GEN y, GEN T)
    1983             : {
    1984             :   long dx, dy, dz, i, p, T_ismonic;
    1985       15188 :   pari_sp av = avma, av2;
    1986             :   GEN y_lead;
    1987             : 
    1988       15188 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    1989       15188 :   dy = degpol(y); y_lead = gel(y,dy+2);
    1990       15188 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    1991             :   /* if monic, no point in using pseudo-division */
    1992       15188 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    1993       11856 :   T_ismonic = gequal1(leading_coeff(T));
    1994       11856 :   dx = degpol(x);
    1995       11856 :   if (dx < dy) return !signe(x);
    1996       11856 :   (void)new_chunk(2);
    1997       11856 :   x = RgX_recip_i(x)+2;
    1998       11856 :   y = RgX_recip_i(y)+2;
    1999             :   /* pay attention to sparse divisors */
    2000       24232 :   for (i = 1; i <= dy; i++)
    2001       12376 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    2002       11856 :   dz = dx-dy; p = dz+1;
    2003       11856 :   av2 = avma;
    2004             :   for (;;)
    2005      111383 :   {
    2006      123239 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    2007      123239 :     x0 = gneg(x0); p--;
    2008      123239 :     m = gcdii(cx, y0);
    2009      123239 :     if (!equali1(m))
    2010             :     {
    2011      121380 :       x0 = gdiv(x0, m);
    2012      121380 :       y0 = diviiexact(y0, m);
    2013      121379 :       if (equali1(y0)) y0 = NULL;
    2014             :     }
    2015      250418 :     for (i=1; i<=dy; i++)
    2016             :     {
    2017      127179 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2018      127179 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    2019      127180 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2020      127180 :       gel(x,i) = c;
    2021             :     }
    2022     1249303 :     for (   ; i<=dx; i++)
    2023             :     {
    2024     1126064 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    2025     1126064 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    2026     1126064 :       gel(x,i) = c;
    2027             :     }
    2028      135160 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    2029      123239 :     if (dx < dy) break;
    2030      111383 :     if (gc_needed(av2,1))
    2031             :     {
    2032          28 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    2033          28 :       gerepilecoeffs(av2,x,dx+1);
    2034             :     }
    2035             :   }
    2036       11856 :   return gc_bool(av, dx < 0);
    2037             : }
    2038             : 
    2039             : /* T either NULL or a t_POL. */
    2040             : GEN
    2041      134860 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    2042             : {
    2043      134860 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    2044      134860 :   pari_sp av = avma, av2;
    2045             :   GEN y_lead;
    2046             : 
    2047      134860 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    2048      134860 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2049             :   /* if monic, no point in using pseudo-division */
    2050      134860 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    2051      129261 :   dx = degpol(x);
    2052      129262 :   if (dx < dy) return RgX_copy(x);
    2053      129255 :   (void)new_chunk(2);
    2054      129255 :   x = RgX_recip_i(x)+2;
    2055      129254 :   y = RgX_recip_i(y)+2;
    2056             :   /* pay attention to sparse divisors */
    2057      369382 :   for (i = 1; i <= dy; i++)
    2058      240128 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2059      129254 :   dz = dx-dy; p = dz+1;
    2060      129254 :   av2 = avma;
    2061             :   for (;;)
    2062             :   {
    2063      360491 :     gel(x,0) = gneg(gel(x,0)); p--;
    2064     1070490 :     for (i=1; i<=dy; i++)
    2065             :     {
    2066      710007 :       GEN c = gmul(y_lead, gel(x,i));
    2067      709957 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    2068      709985 :       gel(x,i) = rem(c, T);
    2069             :     }
    2070      848966 :     for (   ; i<=dx; i++)
    2071             :     {
    2072      488478 :       GEN c = gmul(y_lead, gel(x,i));
    2073      488468 :       gel(x,i) = rem(c, T);
    2074             :     }
    2075      387987 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    2076      360491 :     if (dx < dy) break;
    2077      231237 :     if (gc_needed(av2,1))
    2078             :     {
    2079           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    2080           0 :       gerepilecoeffs(av2,x,dx+1);
    2081             :     }
    2082             :   }
    2083      129254 :   if (dx < 0) return pol_0(vx);
    2084      127721 :   lx = dx+3; x -= 2;
    2085      127721 :   x[0] = evaltyp(t_POL) | evallg(lx);
    2086      127720 :   x[1] = evalsigne(1) | evalvarn(vx);
    2087      127720 :   x = RgX_recip_i(x);
    2088      127720 :   if (p)
    2089             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    2090       12239 :     GEN t = y_lead;
    2091       12239 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    2092           0 :       t = RgXQ_powu(t, p, T);
    2093             :     else
    2094       12239 :       t = gpowgs(t, p);
    2095       32833 :     for (i=2; i<lx; i++)
    2096             :     {
    2097       20594 :       GEN c = gmul(gel(x,i), t);
    2098       20594 :       gel(x,i) = rem(c,T);
    2099             :     }
    2100       12239 :     if (!T) return gerepileupto(av, x);
    2101             :   }
    2102      115481 :   return gerepilecopy(av, x);
    2103             : }
    2104             : 
    2105             : GEN
    2106      134860 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2107             : 
    2108             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2109             : GEN
    2110        7897 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2111             : {
    2112        7897 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2113        7897 :   pari_sp av = avma, av2;
    2114             :   GEN z, r, ypow, y_lead;
    2115             : 
    2116        7897 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2117        7897 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2118        7897 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2119        7323 :   dx = degpol(x);
    2120        7323 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2121        7323 :   if (dx == dy)
    2122             :   {
    2123          70 :     GEN x_lead = gel(x,lg(x)-1);
    2124          70 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2125          70 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2126          70 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2127          70 :     *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
    2128             :   }
    2129        7253 :   (void)new_chunk(2);
    2130        7253 :   x = RgX_recip_i(x)+2;
    2131        7253 :   y = RgX_recip_i(y)+2;
    2132             :   /* pay attention to sparse divisors */
    2133       36948 :   for (i = 1; i <= dy; i++)
    2134       29695 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2135        7253 :   dz = dx-dy; p = dz+1;
    2136        7253 :   lz = dz+3;
    2137        7253 :   z = cgetg(lz, t_POL);
    2138        7253 :   z[1] = evalsigne(1) | evalvarn(vx);
    2139       27759 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2140        7253 :   ypow = new_chunk(dz+1);
    2141        7253 :   gel(ypow,0) = gen_1;
    2142        7253 :   gel(ypow,1) = y_lead;
    2143       13253 :   for (i=2; i<=dz; i++)
    2144             :   {
    2145        6000 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2146        6000 :     gel(ypow,i) = rem(c,T);
    2147             :   }
    2148        7253 :   av2 = avma;
    2149        7253 :   for (iz=2;;)
    2150             :   {
    2151        7869 :     p--;
    2152       15122 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2153       65385 :     for (i=1; i<=dy; i++)
    2154             :     {
    2155       50263 :       GEN c = gmul(y_lead, gel(x,i));
    2156       50263 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2157       50263 :       gel(x,i) = rem(c, T);
    2158             :     }
    2159       39344 :     for (   ; i<=dx; i++)
    2160             :     {
    2161       24222 :       GEN c = gmul(y_lead, gel(x,i));
    2162       24222 :       gel(x,i) = rem(c,T);
    2163             :     }
    2164       15122 :     x++; dx--;
    2165       20506 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2166       15122 :     if (dx < dy) break;
    2167        7869 :     if (gc_needed(av2,1))
    2168             :     {
    2169           0 :       GEN X = x-2;
    2170           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2171           0 :       X[0] = evaltyp(t_POL)|evallg(dx+3); X[1] = z[1]; /* hack */
    2172           0 :       gerepileall(av2,2, &X, &z); x = X+2;
    2173             :     }
    2174             :   }
    2175       12993 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2176        7253 :   if (dx < 0)
    2177          98 :     x = pol_0(vx);
    2178             :   else
    2179             :   {
    2180        7155 :     lx = dx+3; x -= 2;
    2181        7155 :     x[0] = evaltyp(t_POL) | evallg(lx);
    2182        7155 :     x[1] = evalsigne(1) | evalvarn(vx);
    2183        7155 :     x = RgX_recip_i(x);
    2184             :   }
    2185        7253 :   z = RgX_recip_i(z);
    2186        7253 :   r = x;
    2187        7253 :   if (p)
    2188             :   {
    2189        3241 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2190        3241 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2191             :   }
    2192        7253 :   gerepileall(av, 2, &z, &r);
    2193        7253 :   *ptr = r; return z;
    2194             : }
    2195             : GEN
    2196        7666 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2197        7666 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2198             : 
    2199             : GEN
    2200           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2201             : {
    2202           0 :   return RgXQX_red(RgX_mul(x,y), T);
    2203             : }
    2204             : GEN
    2205   133950905 : RgX_Rg_mul(GEN y, GEN x) {
    2206             :   long i, ly;
    2207   133950905 :   GEN z = cgetg_copy(y, &ly); z[1] = y[1];
    2208   133952153 :   if (ly == 2) return z;
    2209   546458100 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2210   133760271 :   return normalizepol_lg(z,ly);
    2211             : }
    2212             : GEN
    2213        2835 : RgX_muls(GEN y, long x) {
    2214             :   long i, ly;
    2215        2835 :   GEN z = cgetg_copy(y, &ly); z[1] = y[1];
    2216        2835 :   if (ly == 2) return z;
    2217       18739 :   for (i = 2; i < ly; i++) gel(z,i) = gmulsg(x,gel(y,i));
    2218        2380 :   return normalizepol_lg(z,ly);
    2219             : }
    2220             : GEN
    2221         938 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T)
    2222             : {
    2223         938 :   return RgXQX_red(RgX_Rg_mul(x,y), T);
    2224             : }
    2225             : GEN
    2226         133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T)
    2227             : {
    2228         133 :   return RgXQV_red(RgV_Rg_mul(v,x), T);
    2229             : }
    2230             : 
    2231             : GEN
    2232           0 : RgXQX_sqr(GEN x, GEN T)
    2233             : {
    2234           0 :   return RgXQX_red(RgX_sqr(x), T);
    2235             : }
    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      596419 : _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       70066 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2255             : static GEN
    2256       31598 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
    2257             : static GEN
    2258      216497 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2259             : static GEN
    2260      994268 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2261             : static GEN
    2262      969779 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2263             : static GEN
    2264         826 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2265             : static GEN
    2266      696504 : _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      397708 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2279             : {
    2280      397708 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2281      397707 :   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       11851 : _sqrXn(void *data, GEN x) {
    2291       11851 :   struct modXn *S = (struct modXn*)data;
    2292       11851 :   return RgXn_sqr(x, S->n);
    2293             : }
    2294             : static GEN
    2295        4493 : _mulXn(void *data, GEN x, GEN y) {
    2296        4493 :   struct modXn *S = (struct modXn*)data;
    2297        4493 :   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        2279 : RgXn_powu_i(GEN x, ulong m, long n)
    2324             : {
    2325             :   struct modXn S;
    2326             :   long v;
    2327        2279 :   if (n == 1) return x;
    2328        2237 :   v = RgX_valrem(x, &x);
    2329        2237 :   if (v) n -= m * v;
    2330        2237 :   S.v = varn(x); S.n = n;
    2331        2237 :   x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2332        2237 :   if (v) x = RgX_shift_shallow(x, m * v);
    2333        2237 :   return x;
    2334             : }
    2335             : GEN
    2336           0 : RgXn_powu(GEN x, ulong m, long n)
    2337             : {
    2338             :   pari_sp av;
    2339           0 :   if (n == 1) 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        4361 : RgXn_eval(GEN Q, GEN x, long n)
    2363             : {
    2364        4361 :   long d = degpol(x);
    2365             :   int use_sqr;
    2366             :   struct modXn S;
    2367        4361 :   if (d == 1 && isrationalzero(gel(x,2)))
    2368             :   {
    2369        4354 :     GEN y = RgX_unscale(Q, gel(x,3));
    2370        4354 :     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     1676886 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2381             : {
    2382     1676886 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2383     1676886 :   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           0 : RgXn_sqrhigh(GEN f, long n2, long n)
    2389             : {
    2390           0 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2391           0 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2392             : }
    2393             : 
    2394             : GEN
    2395     1004898 : RgXn_inv_i(GEN f, long e)
    2396             : {
    2397             :   pari_sp av;
    2398             :   ulong mask;
    2399             :   GEN W, a;
    2400     1004898 :   long v = varn(f), n = 1;
    2401             : 
    2402     1004898 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2403     1004891 :   a = ginv(gel(f,2));
    2404     1004888 :   if (e == 1) return scalarpol(a, v);
    2405     1000702 :   else if (e == 2)
    2406             :   {
    2407             :     GEN b;
    2408      224064 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2409      201652 :     b = gneg(b);
    2410      201651 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2411      201651 :     return deg1pol(b, a, v);
    2412             :   }
    2413      776638 :   av = avma;
    2414      776638 :   W = scalarpol_shallow(a,v);
    2415      776638 :   mask = quadratic_prec_mask(e);
    2416     2448414 :   while (mask > 1)
    2417             :   {
    2418             :     GEN u, fr;
    2419     1671776 :     long n2 = n;
    2420     1671776 :     n<<=1; if (mask & 1) n--;
    2421     1671776 :     mask >>= 1;
    2422     1671776 :     fr = RgXn_red_shallow(f, n);
    2423     1671776 :     u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2424     1671776 :     W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2425     1671776 :     if (gc_needed(av,2))
    2426             :     {
    2427           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2428           0 :       W = gerepileupto(av, W);
    2429             :     }
    2430             :   }
    2431      776638 :   return W;
    2432             : }
    2433             : 
    2434             : static GEN
    2435           0 : RgXn_inv_FpX(GEN x, long e, GEN p)
    2436             : {
    2437           0 :   pari_sp av = avma;
    2438             :   GEN r;
    2439           0 :   if (lgefint(p) == 3)
    2440             :   {
    2441           0 :     ulong pp = uel(p, 2);
    2442           0 :     if (pp == 2)
    2443           0 :       r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
    2444             :     else
    2445           0 :       r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
    2446             :   }
    2447             :   else
    2448           0 :     r = FpXn_inv(RgX_to_FpX(x, p), e, p);
    2449           0 :   return gerepileupto(av, FpX_to_mod(r, p));
    2450             : }
    2451             : 
    2452             : static GEN
    2453           0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
    2454             : {
    2455           0 :   pari_sp av = avma;
    2456           0 :   GEN r, T = RgX_to_FpX(pol, p);
    2457           0 :   if (signe(T) == 0) pari_err_OP("/", gen_1, x);
    2458           0 :   r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
    2459           0 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    2460             : }
    2461             : 
    2462             : #define code(t1,t2) ((t1 << 6) | t2)
    2463             : 
    2464             : static GEN
    2465       79216 : RgXn_inv_fast(GEN x, long e)
    2466             : {
    2467             :   GEN p, pol;
    2468             :   long pa;
    2469       79216 :   long t = RgX_type(x,&p,&pol,&pa);
    2470       79215 :   switch(t)
    2471             :   {
    2472           0 :     case t_INTMOD: return RgXn_inv_FpX(x, e, p);
    2473           0 :     case code(t_POLMOD, t_INTMOD):
    2474           0 :                    return RgXn_inv_FpXQX(x, e, pol, p);
    2475       79215 :     default:       return NULL;
    2476             :   }
    2477             : }
    2478             : #undef code
    2479             : 
    2480             : GEN
    2481       79216 : RgXn_inv(GEN f, long e)
    2482             : {
    2483       79216 :   pari_sp av = avma;
    2484       79216 :   GEN h = RgXn_inv_fast(f, e);
    2485       79215 :   if (h) return h;
    2486       79215 :   return gerepileupto(av, RgXn_inv_i(f, e));
    2487             : }
    2488             : 
    2489             : /* Compute intformal(x^n*S)/x^(n+1) */
    2490             : static GEN
    2491       21021 : RgX_integXn(GEN x, long n)
    2492             : {
    2493       21021 :   long i, lx = lg(x);
    2494             :   GEN y;
    2495       21021 :   if (lx == 2) return RgX_copy(x);
    2496       17717 :   y = cgetg(lx, t_POL); y[1] = x[1];
    2497       40432 :   for (i=2; i<lx; i++)
    2498       22715 :     gel(y,i) = gdivgs(gel(x,i), n+i-1);
    2499       17717 :   return RgX_renormalize_lg(y, lx);;
    2500             : }
    2501             : 
    2502             : GEN
    2503       15911 : RgXn_expint(GEN h, long e)
    2504             : {
    2505       15911 :   pari_sp av = avma, av2;
    2506       15911 :   long v = varn(h), n;
    2507       15911 :   GEN f = pol_1(v), g;
    2508             :   ulong mask;
    2509             : 
    2510       15911 :   if (!signe(h)) return f;
    2511       15911 :   g = pol_1(v);
    2512       15911 :   n = 1; mask = quadratic_prec_mask(e);
    2513       15911 :   av2 = avma;
    2514       21021 :   for (;mask>1;)
    2515             :   {
    2516             :     GEN u, w;
    2517       21021 :     long n2 = n;
    2518       21021 :     n<<=1; if (mask & 1) n--;
    2519       21021 :     mask >>= 1;
    2520       21021 :     u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
    2521       21021 :     u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
    2522       21021 :     w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
    2523       21021 :     f = RgX_add(f, RgX_shift_shallow(w, n2));
    2524       21021 :     if (mask<=1) break;
    2525        5110 :     u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
    2526        5110 :     g = RgX_sub(g, RgX_shift_shallow(u, n2));
    2527        5110 :     if (gc_needed(av2,2))
    2528             :     {
    2529           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
    2530           0 :       gerepileall(av2, 2, &f, &g);
    2531             :     }
    2532             :   }
    2533       15911 :   return gerepileupto(av, f);
    2534             : }
    2535             : 
    2536             : GEN
    2537           0 : RgXn_exp(GEN h, long e)
    2538             : {
    2539           0 :   long d = degpol(h);
    2540           0 :   if (d < 0) return pol_1(varn(h));
    2541           0 :   if (!d || !gequal0(gel(h,2)))
    2542           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2543           0 :   return RgXn_expint(RgX_deriv(h), e);
    2544             : }
    2545             : 
    2546             : GEN
    2547         154 : RgXn_reverse(GEN f, long e)
    2548             : {
    2549         154 :   pari_sp av = avma, av2;
    2550             :   ulong mask;
    2551             :   GEN fi, a, df, W, an;
    2552         154 :   long v = varn(f), n=1;
    2553         154 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2554           0 :     pari_err_INV("serreverse",f);
    2555         154 :   fi = ginv(gel(f,3));
    2556         154 :   a = deg1pol_shallow(fi,gen_0,v);
    2557         154 :   if (e <= 2) return gerepilecopy(av, a);
    2558         133 :   W = scalarpol(fi,v);
    2559         133 :   df = RgX_deriv(f);
    2560         133 :   mask = quadratic_prec_mask(e);
    2561         133 :   av2 = avma;
    2562         616 :   for (;mask>1;)
    2563             :   {
    2564             :     GEN u, fa, fr;
    2565         483 :     long n2 = n, rt;
    2566         483 :     n<<=1; if (mask & 1) n--;
    2567         483 :     mask >>= 1;
    2568         483 :     fr = RgXn_red_shallow(f, n);
    2569         483 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2570         483 :     an = RgXn_powers(a, rt, n);
    2571         483 :     if (n>1)
    2572             :     {
    2573         483 :       long n4 = (n2+1)>>1;
    2574         483 :       GEN dfr = RgXn_red_shallow(df, n2);
    2575         483 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2576         483 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2577         483 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2578             :     }
    2579         483 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2580         483 :     fa = RgX_shift(fa, -n2);
    2581         483 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2582         483 :     if (gc_needed(av2,2))
    2583             :     {
    2584           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2585           0 :       gerepileall(av2, 2, &a, &W);
    2586             :     }
    2587             :   }
    2588         133 :   return gerepileupto(av, a);
    2589             : }
    2590             : 
    2591             : GEN
    2592           0 : RgXn_sqrt(GEN h, long e)
    2593             : {
    2594           0 :   pari_sp av = avma, av2;
    2595           0 :   long v = varn(h), n = 1;
    2596           0 :   GEN f = scalarpol(gen_1, v), df = f;
    2597           0 :   ulong mask = quadratic_prec_mask(e);
    2598           0 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2599           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2600           0 :   av2 = avma;
    2601             :   while(1)
    2602           0 :   {
    2603           0 :     long n2 = n, m;
    2604             :     GEN g;
    2605           0 :     n<<=1; if (mask & 1) n--;
    2606           0 :     mask >>= 1;
    2607           0 :     m = n-n2;
    2608           0 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2609           0 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2610           0 :     if (mask==1) return gerepileupto(av, f);
    2611           0 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2612           0 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2613           0 :     if (gc_needed(av2,2))
    2614             :     {
    2615           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2616           0 :       gerepileall(av2, 2, &f, &df);
    2617             :     }
    2618             :   }
    2619             : }
    2620             : 
    2621             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2622             : GEN
    2623      112735 : RgXQ_powu(GEN x, ulong n, GEN T)
    2624             : {
    2625      112735 :   pari_sp av = avma;
    2626             : 
    2627      112735 :   if (!n) return pol_1(varn(x));
    2628       65481 :   if (n == 1) return RgX_copy(x);
    2629       17926 :   x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
    2630       17927 :   return gerepilecopy(av, x);
    2631             : }
    2632             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2633             : GEN
    2634       52653 : RgXQ_pow(GEN x, GEN n, GEN T)
    2635             : {
    2636             :   pari_sp av;
    2637       52653 :   long s = signe(n);
    2638             : 
    2639       52653 :   if (!s) return pol_1(varn(x));
    2640       52653 :   if (is_pm1(n) == 1)
    2641       31598 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2642       21055 :   av = avma;
    2643       21055 :   if (s < 0) x = RgXQ_inv(x, T);
    2644       21055 :   x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
    2645       21055 :   return gerepilecopy(av, x);
    2646             : }
    2647             : static GEN
    2648      148541 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
    2649             : static GEN
    2650       86595 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
    2651             : 
    2652             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2653             : GEN
    2654       10542 : ZXQ_powers(GEN x, long l, GEN T)
    2655             : {
    2656       10542 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2657       10542 :   return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
    2658             : }
    2659             : 
    2660             : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
    2661             : GEN
    2662      138433 : ZXQ_powu(GEN x, ulong n, GEN T)
    2663             : {
    2664      138433 :   pari_sp av = avma;
    2665             : 
    2666      138433 :   if (!n) return pol_1(varn(x));
    2667      138433 :   if (n == 1) return ZX_copy(x);
    2668       90051 :   x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
    2669       90059 :   return gerepilecopy(av, x);
    2670             : }
    2671             : 
    2672             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2673             : GEN
    2674        4186 : RgXQ_powers(GEN x, long l, GEN T)
    2675             : {
    2676        4186 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2677        4186 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2678             : }
    2679             : 
    2680             : GEN
    2681        2954 : RgXQV_factorback(GEN L, GEN e, GEN T)
    2682             : {
    2683        2954 :   return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
    2684             : }
    2685             : 
    2686             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2687             : GEN
    2688        8029 : QXQ_powers(GEN a, long n, GEN T)
    2689             : {
    2690        8029 :   GEN den, v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
    2691             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2692        8029 :   if (den)
    2693             :   { /* restore denominators */
    2694        5089 :     GEN d = den;
    2695             :     long i;
    2696        5089 :     gel(v,2) = a;
    2697       25753 :     for (i=3; i<=n+1; i++) {
    2698       20664 :       d = mulii(d,den);
    2699       20663 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2700             :     }
    2701             :   }
    2702        8029 :   return v;
    2703             : }
    2704             : 
    2705             : static GEN
    2706        2317 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2707             : {
    2708        2317 :   long l, i, m = 0;
    2709             :   GEN dz, z;
    2710        2317 :   GEN V = cgetg_copy(v, &l);
    2711        8330 :   for (i = imin; i < l; i++)
    2712             :   {
    2713        6013 :     GEN c = gel(v, i);
    2714        6013 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2715             :   }
    2716        2317 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2717        2611 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2718        8330 :   for (i = imin; i < l; i++)
    2719             :   {
    2720        6013 :     GEN c = gel(v,i);
    2721        6013 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2722        6013 :     gel(V,i) = c;
    2723             :   }
    2724        2317 :   return V;
    2725             : }
    2726             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2727             : GEN
    2728        2023 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2729        2023 : { return do_QXQ_eval(v, 1, a, T); }
    2730             : 
    2731             : GEN
    2732         294 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
    2733         294 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2734             : 
    2735             : GEN
    2736        2884 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2737             : {
    2738        2884 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2739             : }
    2740             : 
    2741             : GEN
    2742        4998 : RgXQ_norm(GEN x, GEN T)
    2743             : {
    2744             :   pari_sp av;
    2745        4998 :   long dx = degpol(x);
    2746             :   GEN L, y;
    2747             : 
    2748        4998 :   av = avma; y = resultant(T, x);
    2749        4998 :   L = leading_coeff(T);
    2750        4998 :   if (gequal1(L) || !signe(x)) return y;
    2751           0 :   return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
    2752             : }
    2753             : 
    2754             : GEN
    2755         406 : RgXQ_trace(GEN x, GEN T)
    2756             : {
    2757         406 :   pari_sp av = avma;
    2758         406 :   GEN dT = RgX_deriv(T);
    2759         406 :   long n = degpol(dT);
    2760         406 :   GEN z = RgXQ_mul(x, dT, T);
    2761         406 :   if (degpol(z)<n) return gc_const(av, gen_0);
    2762         357 :   return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
    2763             : }
    2764             : 
    2765             : GEN
    2766     1823922 : RgX_blocks(GEN P, long n, long m)
    2767             : {
    2768     1823922 :   GEN z = cgetg(m+1,t_VEC);
    2769     1823922 :   long i,j, k=2, l = lg(P);
    2770     5664812 :   for(i=1; i<=m; i++)
    2771             :   {
    2772     3840890 :     GEN zi = cgetg(n+2,t_POL);
    2773     3840890 :     zi[1] = P[1];
    2774     3840890 :     gel(z,i) = zi;
    2775    12234566 :     for(j=2; j<n+2; j++)
    2776     8393676 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2777     3840890 :     zi = RgX_renormalize_lg(zi, n+2);
    2778             :   }
    2779     1823922 :   return z;
    2780             : }
    2781             : 
    2782             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2783             : void
    2784    36027653 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2785             : {
    2786    36027653 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2787             :   GEN p0, p1;
    2788             : 
    2789    36027251 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2790             : 
    2791    36026964 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2792    36026964 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2793    36029204 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2794   117160335 :   for (i=0; i<n1; i++)
    2795             :   {
    2796    81128633 :     p0[2+i] = p[2+(i<<1)];
    2797    81128633 :     p1[2+i] = p[3+(i<<1)];
    2798             :   }
    2799    36031702 :   if (n1 != n0)
    2800    16972628 :     p0[2+i] = p[2+(i<<1)];
    2801    36031702 :   *pe = normalizepol(p0);
    2802    36025784 :   *po = normalizepol(p1);
    2803             : }
    2804             : 
    2805             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2806             : GEN
    2807       40922 : RgX_splitting(GEN p, long k)
    2808             : {
    2809       40922 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2810             :   GEN r;
    2811             : 
    2812       40922 :   m = n/k;
    2813       40922 :   r = cgetg(k+1,t_VEC);
    2814      225666 :   for(i=1; i<=k; i++)
    2815             :   {
    2816      184744 :     gel(r,i) = cgetg(m+3, t_POL);
    2817      184744 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2818             :   }
    2819      554981 :   for (j=1, i=0, l=2; i<=n; i++)
    2820             :   {
    2821      514059 :     gmael(r,j,l) = gel(p,2+i);
    2822      514059 :     if (j==k) { j=1; l++; } else j++;
    2823             :   }
    2824      225666 :   for(i=1; i<=k; i++)
    2825      184744 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2826       40922 :   return r;
    2827             : }
    2828             : 
    2829             : /*******************************************************************/
    2830             : /*                                                                 */
    2831             : /*                        Kronecker form                           */
    2832             : /*                                                                 */
    2833             : /*******************************************************************/
    2834             : 
    2835             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2836             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2837             :  * normalized coefficients */
    2838             : GEN
    2839      176285 : Kronecker_to_mod(GEN z, GEN T)
    2840             : {
    2841      176285 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2842      176285 :   GEN x, t = cgetg(N,t_POL);
    2843      176285 :   t[1] = T[1];
    2844      176285 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2845      176285 :   x[1] = z[1];
    2846      176285 :   T = RgX_copy(T);
    2847     1355497 :   for (i=2; i<lx+2; i++, z+= N-2)
    2848             :   {
    2849     5600108 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2850     1179212 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2851             :   }
    2852      176285 :   N = (l-2) % (N-2) + 2;
    2853      454673 :   for (j=2; j<N; j++) t[j] = z[j];
    2854      176285 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2855      176285 :   return normalizepol_lg(x, i+1);
    2856             : }
    2857             : 
    2858             : /*******************************************************************/
    2859             : /*                                                                 */
    2860             : /*                        Domain detection                         */
    2861             : /*                                                                 */
    2862             : /*******************************************************************/
    2863             : 
    2864             : static GEN
    2865      151642 : zero_FpX_mod(GEN p, long v)
    2866             : {
    2867      151642 :   GEN r = cgetg(3,t_POL);
    2868      151642 :   r[1] = evalvarn(v);
    2869      151642 :   gel(r,2) = mkintmod(gen_0, icopy(p));
    2870      151642 :   return r;
    2871             : }
    2872             : 
    2873             : static GEN
    2874      500345 : RgX_mul_FpX(GEN x, GEN y, GEN p)
    2875             : {
    2876      500345 :   pari_sp av = avma;
    2877             :   GEN r;
    2878      500345 :   if (lgefint(p) == 3)
    2879             :   {
    2880      470213 :     ulong pp = uel(p, 2);
    2881      470213 :     r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
    2882             :                                   RgX_to_Flx(y, pp), pp));
    2883             :   }
    2884             :   else
    2885       30132 :     r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2886      500345 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2887      386412 :   return gerepileupto(av, FpX_to_mod(r, p));
    2888             : }
    2889             : 
    2890             : static GEN
    2891         532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
    2892             : {
    2893         532 :   GEN r = cgetg(3,t_POL);
    2894         532 :   r[1] = evalvarn(v);
    2895         532 :   gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
    2896         532 :   return r;
    2897             : }
    2898             : 
    2899             : static GEN
    2900        2037 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2901             : {
    2902        2037 :   pari_sp av = avma;
    2903             :   long dT;
    2904             :   GEN kx, ky, r;
    2905        2037 :   GEN T = RgX_to_FpX(pol, p);
    2906        2037 :   if (signe(T)==0) pari_err_OP("*", x, y);
    2907        2030 :   dT = degpol(T);
    2908        2030 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    2909        2030 :   ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
    2910        2030 :   r = FpX_mul(kx, ky, p);
    2911        2030 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    2912        1505 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    2913             : }
    2914             : 
    2915             : static GEN
    2916      409031 : RgX_liftred(GEN x, GEN T)
    2917      409031 : { return RgXQX_red(liftpol_shallow(x), T); }
    2918             : 
    2919             : static GEN
    2920      161104 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
    2921             : {
    2922      161104 :   pari_sp av = avma;
    2923      161104 :   long dT = degpol(T);
    2924      161104 :   GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
    2925             :                  RgXX_to_Kronecker(RgX_liftred(y, T), dT));
    2926      161104 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    2927             : }
    2928             : 
    2929             : static GEN
    2930        1358 : RgX_sqr_FpX(GEN x, GEN p)
    2931             : {
    2932        1358 :   pari_sp av = avma;
    2933             :   GEN r;
    2934        1358 :   if (lgefint(p) == 3)
    2935             :   {
    2936        1154 :     ulong pp = uel(p, 2);
    2937        1154 :     r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
    2938             :   }
    2939             :   else
    2940         204 :     r = FpX_sqr(RgX_to_FpX(x, p), p);
    2941        1358 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2942        1239 :   return gerepileupto(av, FpX_to_mod(r, p));
    2943             : }
    2944             : 
    2945             : static GEN
    2946         196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
    2947             : {
    2948         196 :   pari_sp av = avma;
    2949             :   long dT;
    2950         196 :   GEN kx, r, T = RgX_to_FpX(pol, p);
    2951         196 :   if (signe(T)==0) pari_err_OP("*",x,x);
    2952         189 :   dT = degpol(T);
    2953         189 :   kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
    2954         189 :   r = FpX_sqr(kx, p);
    2955         189 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    2956         189 :   return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
    2957             : }
    2958             : 
    2959             : static GEN
    2960       13487 : RgX_sqr_QXQX(GEN x, GEN T)
    2961             : {
    2962       13487 :   pari_sp av = avma;
    2963       13487 :   long dT = degpol(T);
    2964       13487 :   GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
    2965       13487 :   return gerepileupto(av, Kronecker_to_mod(r, T));
    2966             : }
    2967             : 
    2968             : static GEN
    2969       67466 : RgX_rem_FpX(GEN x, GEN y, GEN p)
    2970             : {
    2971       67466 :   pari_sp av = avma;
    2972             :   GEN r;
    2973       67466 :   if (lgefint(p) == 3)
    2974             :   {
    2975       50676 :     ulong pp = uel(p, 2);
    2976       50676 :     r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
    2977             :                                   RgX_to_Flx(y, pp), pp));
    2978             :   }
    2979             :   else
    2980       16790 :     r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
    2981       67466 :   if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
    2982       29876 :   return gerepileupto(av, FpX_to_mod(r, p));
    2983             : }
    2984             : 
    2985             : static GEN
    2986       36668 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
    2987             : {
    2988       36668 :   pari_sp av = avma;
    2989             :   GEN r;
    2990       36668 :   r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
    2991       36668 :   return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
    2992             : }
    2993             : static GEN
    2994          70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
    2995             : {
    2996          70 :   pari_sp av = avma;
    2997             :   GEN r;
    2998          70 :   GEN T = RgX_to_FpX(pol, p);
    2999          70 :   if (signe(T) == 0) pari_err_OP("%", x, y);
    3000          63 :   if (lgefint(p) == 3)
    3001             :   {
    3002          55 :     ulong pp = uel(p, 2);
    3003          55 :     GEN Tp = ZX_to_Flx(T, pp);
    3004          55 :     r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
    3005             :                               RgX_to_FlxqX(y, Tp, pp), Tp, pp));
    3006             :   }
    3007             :   else
    3008           8 :     r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
    3009          63 :   if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
    3010          56 :   return gerepileupto(av, FpXQX_to_mod(r, T, p));
    3011             : }
    3012             : 
    3013             : #define code(t1,t2) ((t1 << 6) | t2)
    3014             : static GEN
    3015    62724398 : RgX_mul_fast(GEN x, GEN y)
    3016             : {
    3017             :   GEN p, pol;
    3018             :   long pa;
    3019    62724398 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3020    62724558 :   switch(t)
    3021             :   {
    3022    36313478 :     case t_INT:    return ZX_mul(x,y);
    3023     2510779 :     case t_FRAC:   return QX_mul(x,y);
    3024      104665 :     case t_FFELT:  return FFX_mul(x, y, pol);
    3025      500345 :     case t_INTMOD: return RgX_mul_FpX(x, y, p);
    3026      161104 :     case code(t_POLMOD, t_INT):
    3027             :     case code(t_POLMOD, t_FRAC):
    3028      161104 :                    return RgX_mul_QXQX(x, y, pol);
    3029        2025 :     case code(t_POLMOD, t_INTMOD):
    3030        2025 :                    return RgX_mul_FpXQX(x, y, pol, p);
    3031    23132162 :     default:       return NULL;
    3032             :   }
    3033             : }
    3034             : static GEN
    3035      185657 : RgX_sqr_fast(GEN x)
    3036             : {
    3037             :   GEN p, pol;
    3038             :   long pa;
    3039      185657 :   long t = RgX_type(x,&p,&pol,&pa);
    3040      185656 :   switch(t)
    3041             :   {
    3042      138337 :     case t_INT:    return ZX_sqr(x);
    3043       12663 :     case t_FRAC:   return QX_sqr(x);
    3044        2296 :     case t_FFELT:  return FFX_sqr(x, pol);
    3045        1358 :     case t_INTMOD: return RgX_sqr_FpX(x, p);
    3046       13487 :     case code(t_POLMOD, t_INT):
    3047             :     case code(t_POLMOD, t_FRAC):
    3048       13487 :                    return RgX_sqr_QXQX(x, pol);
    3049         194 :     case code(t_POLMOD, t_INTMOD):
    3050         194 :                    return RgX_sqr_FpXQX(x, pol, p);
    3051       17321 :     default:       return NULL;
    3052             :   }
    3053             : }
    3054             : 
    3055             : static GEN
    3056    11907556 : RgX_rem_fast(GEN x, GEN y)
    3057             : {
    3058             :   GEN p, pol;
    3059             :   long pa;
    3060    11907556 :   long t = RgX_type2(x,y, &p,&pol,&pa);
    3061    11907810 :   switch(t)
    3062             :   {
    3063     3222661 :     case t_INT:    return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
    3064     1660500 :     case t_FRAC:   return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
    3065          84 :     case t_FFELT:  return FFX_rem(x, y, pol);
    3066       67466 :     case t_INTMOD: return RgX_rem_FpX(x, y, p);
    3067       36668 :     case code(t_POLMOD, t_INT):
    3068             :     case code(t_POLMOD, t_FRAC):
    3069       36668 :                    return RgX_rem_QXQX(x, y, pol);
    3070          91 :     case code(t_POLMOD, t_INTMOD):
    3071          91 :                    return RgX_rem_FpXQX(x, y, pol, p);
    3072     6920340 :     default:       return NULL;
    3073             :   }
    3074             : }
    3075             : 
    3076             : #undef code
    3077             : 
    3078             : GEN
    3079    53781989 : RgX_mul(GEN x, GEN y)
    3080             : {
    3081    53781989 :   GEN z = RgX_mul_fast(x,y);
    3082    53782019 :   if (!z) z = RgX_mul_i(x,y);
    3083    53781454 :   return z;
    3084             : }
    3085             : 
    3086             : GEN
    3087      173561 : RgX_sqr(GEN x)
    3088             : {
    3089      173561 :   GEN z = RgX_sqr_fast(x);
    3090      173546 :   if (!z) z = RgX_sqr_i(x);
    3091      173549 :   return z;
    3092             : }
    3093             : 
    3094             : GEN
    3095    11907577 : RgX_rem(GEN x, GEN y)
    3096             : {
    3097    11907577 :   GEN z = RgX_rem_fast(x, y);
    3098    11907751 :   if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
    3099    11907597 :   return z;
    3100             : }
    3101             : 
    3102             : GEN
    3103     7244500 : RgXn_mul(GEN f, GEN g, long n)
    3104             : {
    3105     7244500 :   pari_sp av = avma;
    3106     7244500 :   GEN h = RgX_mul_fast(f,g);
    3107     7244500 :   if (!h) return RgXn_mul2(f,g,n);
    3108     1341570 :   if (degpol(h) < n) return h;
    3109      514639 :   return gerepilecopy(av, RgXn_red_shallow(h, n));
    3110             : }
    3111             : 
    3112             : GEN
    3113       12096 : RgXn_sqr(GEN f, long n)
    3114             : {
    3115       12096 :   pari_sp av = avma;
    3116       12096 :   GEN g = RgX_sqr_fast(f);
    3117       12096 :   if (!g) return RgXn_sqr2(f,n);
    3118       11781 :   if (degpol(g) < n) return g;
    3119       10654 :   return gerepilecopy(av, RgXn_red_shallow(g, n));
    3120             : }
    3121             : 
    3122             : /* (f*g) \/ x^n */
    3123             : GEN
    3124     1697907 : RgX_mulhigh_i(GEN f, GEN g, long n)
    3125             : {
    3126     1697907 :   GEN h = RgX_mul_fast(f,g);
    3127     1697907 :   return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
    3128             : }
    3129             : 
    3130             : /* (f*g) \/ x^n */
    3131             : GEN
    3132           0 : RgX_sqrhigh_i(GEN f, long n)
    3133             : {
    3134           0 :   GEN h = RgX_sqr_fast(f);
    3135           0 :   return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
    3136             : }

Generated by: LCOV version 1.13