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 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.10.0 lcov report (development 20777-d2a9243) Lines: 1306 1456 89.7 %
Date: 2017-06-25 05:59:24 Functions: 140 157 89.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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*******************************************************************/
      18             : /*                                                                 */
      19             : /*                         GENERIC                                 */
      20             : /*                                                                 */
      21             : /*******************************************************************/
      22             : 
      23             : /* Return optimal parameter l for the evaluation of n/m polynomials of degree d
      24             :    Fractional values can be used if the evaluations are done with different
      25             :    accuracies, and thus have different weights.
      26             :  */
      27             : long
      28     2275792 : brent_kung_optpow(long d, long n, long m)
      29             : {
      30             :   long p, r;
      31     2275792 :   long pold=1, rold=n*(d-1);
      32    12585092 :   for(p=2; p<=d; p++)
      33             :   {
      34    10309300 :     r = m*(p-1) + n*((d-1)/p);
      35    10309300 :     if (r<rold) { pold=p; rold=r; }
      36             :   }
      37     2275792 :   return pold;
      38             : }
      39             : 
      40             : static GEN
      41     9894554 : gen_RgXQ_eval_powers(GEN P, GEN V, long a, long n, void *E, const struct bb_algebra *ff,
      42             :                                            GEN cmul(void *E, GEN P, long a, GEN x))
      43             : {
      44     9894554 :   pari_sp av = avma;
      45             :   long i;
      46     9894554 :   GEN z = cmul(E,P,a,ff->one(E));
      47     9894441 :   if (!z) z = gen_0;
      48    60917269 :   for (i=1; i<=n; i++)
      49             :   {
      50    51022721 :     GEN t = cmul(E,P,a+i,gel(V,i+1));
      51    51022557 :     if (t) {
      52    49685412 :       z = ff->add(E, z, t);
      53    49685492 :       if (gc_needed(av,2)) z = gerepileupto(av, z);
      54             :     }
      55             :   }
      56     9894548 :   return ff->red(E,z);
      57             : }
      58             : 
      59             : /* Brent & Kung
      60             :  * (Fast algorithms for manipulating formal power series, JACM 25:581-595, 1978)
      61             :  *
      62             :  * V as output by FpXQ_powers(x,l,T,p). For optimal performance, l is as given
      63             :  * by brent_kung_optpow */
      64             : GEN
      65     6397856 : gen_bkeval_powers(GEN P, long d, GEN V, void *E, const struct bb_algebra *ff,
      66             :                                      GEN cmul(void *E, GEN P, long a, GEN x))
      67             : {
      68     6397856 :   pari_sp av = avma;
      69     6397856 :   long l = lg(V)-1;
      70             :   GEN z, u;
      71             : 
      72     6397856 :   if (d < 0) return ff->zero(E);
      73     5891101 :   if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
      74     2377194 :   if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
      75     2377194 :   if (DEBUGLEVEL>=8)
      76             :   {
      77           0 :     long cnt = 1 + (d - l) / (l-1);
      78           0 :     err_printf("RgX_RgXQV_eval(%ld/%ld): %ld RgXQ_mul\n", d, l-1, cnt);
      79             :   }
      80     2377194 :   d -= l;
      81     2377194 :   z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
      82     6380649 :   while (d >= l-1)
      83             :   {
      84     1626262 :     d -= l-1;
      85     1626262 :     u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
      86     1626247 :     z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
      87     1626261 :     if (gc_needed(av,2))
      88          61 :       z = gerepileupto(av, z);
      89             :   }
      90     2377194 :   u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
      91     2377195 :   z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
      92     2377196 :   return gerepileupto(av, ff->red(E,z));
      93             : }
      94             : 
      95             : GEN
      96     1059912 : gen_bkeval(GEN Q, long d, GEN x, int use_sqr, void *E, const struct bb_algebra *ff,
      97             :                                       GEN cmul(void *E, GEN P, long a, GEN x))
      98             : {
      99     1059912 :   pari_sp av = avma;
     100             :   GEN z, V;
     101             :   long rtd;
     102     1059912 :   if (d < 0) return ff->zero(E);
     103     1059807 :   rtd = (long) sqrt((double)d);
     104     1059807 :   V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
     105     1059804 :   z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
     106     1059806 :   return gerepileupto(av, z);
     107             : }
     108             : 
     109             : static GEN
     110      728360 : _gen_nored(void *E, GEN x) { (void)E; return x; }
     111             : static GEN
     112    20663605 : _gen_add(void *E, GEN x, GEN y) { (void)E; return gadd(x, y); }
     113             : static GEN
     114           0 : _gen_sub(void *E, GEN x, GEN y) { (void)E; return gsub(x, y); }
     115             : static GEN
     116      605719 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
     117             : static GEN
     118      193741 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
     119             : static GEN
     120      743088 : _gen_one(void *E) { (void)E; return gen_1; }
     121             : static GEN
     122       14483 : _gen_zero(void *E) { (void)E; return gen_0; }
     123             : 
     124             : static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
     125             :               _gen_mul, _gen_sqr,_gen_one,_gen_zero };
     126             : 
     127             : static GEN
     128      351407 : _gen_cmul(void *E, GEN P, long a, GEN x)
     129      351407 : {(void)E; return gmul(gel(P,a+2), x);}
     130             : 
     131             : GEN
     132      124047 : RgX_RgV_eval(GEN Q, GEN x)
     133             : {
     134      124047 :   return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
     135             : }
     136             : 
     137             : GEN
     138           0 : RgX_Rg_eval_bk(GEN Q, GEN x)
     139             : {
     140           0 :   return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
     141             : }
     142             : 
     143             : GEN
     144        2814 : RgXV_RgV_eval(GEN Q, GEN x)
     145             : {
     146        2814 :   long i, l = lg(Q), vQ = gvar(Q);
     147        2814 :   GEN v = cgetg(l, t_VEC);
     148      237538 :   for (i = 1; i < l; i++)
     149             :   {
     150      234724 :     GEN Qi = gel(Q, i);
     151      234724 :     gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
     152             :   }
     153        2814 :   return v;
     154             : }
     155             : 
     156             : const struct bb_algebra *
     157       92614 : get_Rg_algebra(void)
     158             : {
     159       92614 :   return &Rg_algebra;
     160             : }
     161             : 
     162             : static struct bb_ring Rg_ring = {  _gen_add, _gen_mul, _gen_sqr };
     163             : 
     164             : static GEN
     165           0 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
     166             : {
     167             :   (void) E;
     168           0 :   return RgX_divrem(x, y, r);
     169             : }
     170             : 
     171             : GEN
     172           0 : RgX_digits(GEN x, GEN T)
     173             : {
     174           0 :   pari_sp av = avma;
     175           0 :   long d = degpol(T), n = (lgpol(x)+d-1)/d;
     176           0 :   GEN z = gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
     177           0 :   return gerepileupto(av, z);
     178             : }
     179             : 
     180             : /*******************************************************************/
     181             : /*                                                                 */
     182             : /*                         RgX                                     */
     183             : /*                                                                 */
     184             : /*******************************************************************/
     185             : 
     186             : long
     187    21914694 : RgX_equal(GEN x, GEN y)
     188             : {
     189    21914694 :   long i = lg(x);
     190             : 
     191    21914694 :   if (i != lg(y)) return 0;
     192    93245909 :   for (i--; i > 1; i--)
     193    71392570 :     if (!gequal(gel(x,i),gel(y,i))) return 0;
     194    21853339 :   return 1;
     195             : }
     196             : 
     197             : /* Returns 1 in the base ring over which x is defined */
     198             : /* HACK: this also works for t_SER */
     199             : GEN
     200      673387 : RgX_get_1(GEN x)
     201             : {
     202             :   GEN p, T;
     203      673387 :   long i, lx, tx = RgX_type(x, &p, &T, &lx);
     204      673387 :   if (RgX_type_is_composite(tx))
     205        1253 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     206      673387 :   switch(tx)
     207             :   {
     208          77 :     case t_INTMOD: retmkintmod(gen_1, icopy(p));
     209           7 :     case t_PADIC: return cvtop(gen_1, p, lx);
     210          35 :     case t_FFELT: return FF_1(T);
     211      673268 :     default: return gen_1;
     212             :   }
     213             : }
     214             : /* Returns 0 in the base ring over which x is defined */
     215             : /* HACK: this also works for t_SER */
     216             : GEN
     217      130998 : RgX_get_0(GEN x)
     218             : {
     219             :   GEN p, T;
     220      130998 :   long i, lx, tx = RgX_type(x, &p, &T, &lx);
     221      130998 :   if (RgX_type_is_composite(tx))
     222       13447 :     RgX_type_decode(tx, &i /*junk*/, &tx);
     223      130998 :   switch(tx)
     224             :   {
     225         259 :     case t_INTMOD: retmkintmod(gen_0, icopy(p));
     226           0 :     case t_PADIC: return cvtop(gen_0, p, lx);
     227         210 :     case t_FFELT: return FF_zero(T);
     228      130529 :     default: return gen_0;
     229             :   }
     230             : }
     231             : 
     232             : GEN
     233        2450 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
     234             : {
     235        2450 :   long i, n = degpol(P);
     236             :   GEN z, dz, dP;
     237        2450 :   if (n < 0) return gen_0;
     238        2450 :   P = Q_remove_denom(P, &dP);
     239        2450 :   z = gel(P,2); if (n == 0) return icopy(z);
     240        1337 :   if (dV) z = mulii(dV, z); /* V[1] = dV */
     241        1337 :   z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
     242        1337 :   for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
     243        1337 :   dz = mul_denom(dP, dV);
     244        1337 :   return dz? RgX_Rg_div(z, dz): z;
     245             : }
     246             : 
     247             : /* Return P(h * x), not memory clean */
     248             : GEN
     249        3878 : RgX_unscale(GEN P, GEN h)
     250             : {
     251        3878 :   long i, l = lg(P);
     252        3878 :   GEN hi = gen_1, Q = cgetg(l, t_POL);
     253        3878 :   Q[1] = P[1];
     254        3878 :   if (l == 2) return Q;
     255        3878 :   gel(Q,2) = gcopy(gel(P,2));
     256       21245 :   for (i=3; i<l; i++)
     257             :   {
     258       17367 :     hi = gmul(hi,h);
     259       17367 :     gel(Q,i) = gmul(gel(P,i), hi);
     260             :   }
     261        3878 :   return Q;
     262             : }
     263             : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
     264             : GEN
     265      904015 : ZX_z_unscale(GEN P, long h)
     266             : {
     267      904015 :   long i, l = lg(P);
     268      904015 :   GEN Q = cgetg(l, t_POL);
     269      904015 :   Q[1] = P[1];
     270      904015 :   if (l == 2) return Q;
     271      904015 :   gel(Q,2) = gel(P,2);
     272      904015 :   if (l == 3) return Q;
     273      904015 :   if (h == -1)
     274      306383 :     for (i = 3; i < l; i++)
     275             :     {
     276      299403 :       gel(Q,i) = negi(gel(P,i));
     277      299403 :       if (++i == l) break;
     278      296142 :       gel(Q,i) = gel(P,i);
     279             :     }
     280             :   else
     281             :   {
     282             :     GEN hi;
     283      893774 :     gel(Q,3) = mulis(gel(P,3), h);
     284      893774 :     hi = sqrs(h);
     285     1852802 :     for (i = 4; i < l; i++)
     286             :     {
     287      959028 :       gel(Q,i) = mulii(gel(P,i), hi);
     288      959028 :       if (i != l-1) hi = mulis(hi,h);
     289             :     }
     290             :   }
     291      904015 :   return Q;
     292             : }
     293             : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
     294             : GEN
     295        7336 : ZX_unscale(GEN P, GEN h)
     296             : {
     297             :   long i, l;
     298             :   GEN Q, hi;
     299        7336 :   i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
     300          14 :   l = lg(P); Q = cgetg(l, t_POL);
     301          14 :   Q[1] = P[1];
     302          14 :   if (l == 2) return Q;
     303          14 :   gel(Q,2) = gel(P,2);
     304          14 :   if (l == 3) return Q;
     305          14 :   hi = h;
     306          14 :   gel(Q,3) = mulii(gel(P,3), hi);
     307          91 :   for (i = 4; i < l; i++)
     308             :   {
     309          77 :     hi = mulii(hi,h);
     310          77 :     gel(Q,i) = mulii(gel(P,i), hi);
     311             :   }
     312          14 :   return Q;
     313             : }
     314             : /* P a ZX. Return P(x << n), not memory clean */
     315             : GEN
     316       18133 : ZX_unscale2n(GEN P, long n)
     317             : {
     318       18133 :   long i, ni = n, l = lg(P);
     319       18133 :   GEN Q = cgetg(l, t_POL);
     320       18133 :   Q[1] = P[1];
     321       18133 :   if (l == 2) return Q;
     322       18133 :   gel(Q,2) = gel(P,2);
     323       18133 :   if (l == 3) return Q;
     324       18133 :   gel(Q,3) = shifti(gel(P,3), ni);
     325       73511 :   for (i=4; i<l; i++)
     326             :   {
     327       55378 :     ni += n;
     328       55378 :     gel(Q,i) = shifti(gel(P,i), ni);
     329             :   }
     330       18133 :   return Q;
     331             : }
     332             : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
     333             : GEN
     334        1176 : ZX_unscale_div(GEN P, GEN h)
     335             : {
     336        1176 :   long i, l = lg(P);
     337        1176 :   GEN hi, Q = cgetg(l, t_POL);
     338        1176 :   Q[1] = P[1];
     339        1176 :   if (l == 2) return Q;
     340        1176 :   gel(Q,2) = diviiexact(gel(P,2), h);
     341        1176 :   if (l == 3) return Q;
     342        1176 :   gel(Q,3) = gel(P,3);
     343        1176 :   if (l == 4) return Q;
     344        1176 :   hi = h;
     345        1176 :   gel(Q,4) = mulii(gel(P,4), hi);
     346        5194 :   for (i=5; i<l; i++)
     347             :   {
     348        4018 :     hi = mulii(hi,h);
     349        4018 :     gel(Q,i) = mulii(gel(P,i), hi);
     350             :   }
     351        1176 :   return Q;
     352             : }
     353             : 
     354             : GEN
     355         224 : RgXV_unscale(GEN v, GEN h)
     356             : {
     357             :   long i, l;
     358             :   GEN w;
     359         224 :   if (!h || isint1(h)) return v;
     360         168 :   w = cgetg_copy(v, &l);
     361         168 :   for (i=1; i<l; i++) gel(w,i) = RgX_unscale(gel(v,i), h);
     362         168 :   return w;
     363             : }
     364             : 
     365             : /* Return h^degpol(P) P(x / h), not memory clean */
     366             : GEN
     367        1540 : RgX_rescale(GEN P, GEN h)
     368             : {
     369        1540 :   long i, l = lg(P);
     370        1540 :   GEN Q = cgetg(l,t_POL), hi = h;
     371        1540 :   Q[l-1] = P[l-1];
     372        8099 :   for (i=l-2; i>=2; i--)
     373             :   {
     374        8099 :     gel(Q,i) = gmul(gel(P,i), hi);
     375        8099 :     if (i == 2) break;
     376        6559 :     hi = gmul(hi,h);
     377             :   }
     378        1540 :   Q[1] = P[1]; return Q;
     379             : }
     380             : 
     381             : /* A(X^d) --> A(X) */
     382             : GEN
     383      109528 : RgX_deflate(GEN x0, long d)
     384             : {
     385             :   GEN z, y, x;
     386      109528 :   long i,id, dy, dx = degpol(x0);
     387      109528 :   if (d == 1 || dx <= 0) return leafcopy(x0);
     388       63025 :   dy = dx/d;
     389       63025 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     390       63025 :   z = y + 2;
     391       63025 :   x = x0+ 2;
     392       63025 :   for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
     393       63025 :   return y;
     394             : }
     395             : 
     396             : /* return x0(X^d) */
     397             : GEN
     398      250673 : RgX_inflate(GEN x0, long d)
     399             : {
     400      250673 :   long i, id, dy, dx = degpol(x0);
     401      250673 :   GEN x = x0 + 2, z, y;
     402      250673 :   if (dx <= 0) return leafcopy(x0);
     403      247803 :   dy = dx*d;
     404      247803 :   y = cgetg(dy+3, t_POL); y[1] = x0[1];
     405      247804 :   z = y + 2;
     406      247804 :   for (i=0; i<=dy; i++) gel(z,i) = gen_0;
     407      247804 :   for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
     408      247804 :   return y;
     409             : }
     410             : 
     411             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     412             : GEN
     413     1017193 : RgX_translate(GEN P, GEN c)
     414             : {
     415     1017193 :   pari_sp av = avma;
     416             :   GEN Q, *R;
     417             :   long i, k, n;
     418             : 
     419     1017193 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     420     1014189 :   Q = leafcopy(P);
     421     1014189 :   R = (GEN*)(Q+2); n = degpol(P);
     422     1014189 :   if (gequal1(c))
     423             :   {
     424        2072 :     for (i=1; i<=n; i++)
     425             :     {
     426        1799 :       for (k=n-i; k<n; k++) R[k] = gadd(R[k], R[k+1]);
     427        1799 :       if (gc_needed(av,2))
     428             :       {
     429           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"TR_POL(1), i = %ld/%ld", i,n);
     430           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     431             :       }
     432             :     }
     433             :   }
     434     1013916 :   else if (gequalm1(c))
     435             :   {
     436      134169 :     for (i=1; i<=n; i++)
     437             :     {
     438      114800 :       for (k=n-i; k<n; k++) R[k] = gsub(R[k], R[k+1]);
     439      114800 :       if (gc_needed(av,2))
     440             :       {
     441           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"TR_POL(-1), i = %ld/%ld", i,n);
     442           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     443             :       }
     444             :     }
     445             :   }
     446             :   else
     447             :   {
     448     3420625 :     for (i=1; i<=n; i++)
     449             :     {
     450     2426078 :       for (k=n-i; k<n; k++) R[k] = gadd(R[k], gmul(c, R[k+1]));
     451     2426078 :       if (gc_needed(av,2))
     452             :       {
     453           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"TR_POL, i = %ld/%ld", i,n);
     454           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     455             :       }
     456             :     }
     457             :   }
     458     1014189 :   return gerepilecopy(av, Q);
     459             : }
     460             : 
     461             : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
     462             : GEN
     463      376147 : ZX_translate(GEN P, GEN c)
     464             : {
     465      376147 :   pari_sp av = avma;
     466             :   GEN Q, *R;
     467             :   long i, k, n;
     468             : 
     469      376147 :   if (!signe(P) || !signe(c)) return ZX_copy(P);
     470      376112 :   Q = leafcopy(P);
     471      376112 :   R = (GEN*)(Q+2); n = degpol(P);
     472      376112 :   if (equali1(c))
     473             :   {
     474     2356216 :     for (i=1; i<=n; i++)
     475             :     {
     476     2083272 :       for (k=n-i; k<n; k++) R[k] = addii(R[k], R[k+1]);
     477     2083272 :       if (gc_needed(av,2))
     478             :       {
     479           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(1), i = %ld/%ld", i,n);
     480           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     481             :       }
     482             :     }
     483             :   }
     484      103168 :   else if (equalim1(c))
     485             :   {
     486          70 :     for (i=1; i<=n; i++)
     487             :     {
     488          49 :       for (k=n-i; k<n; k++) R[k] = subii(R[k], R[k+1]);
     489          49 :       if (gc_needed(av,2))
     490             :       {
     491           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(-1), i = %ld/%ld", i,n);
     492           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     493             :       }
     494             :     }
     495             :   }
     496             :   else
     497             :   {
     498      768881 :     for (i=1; i<=n; i++)
     499             :     {
     500      665734 :       for (k=n-i; k<n; k++) R[k] = addii(R[k], mulii(c, R[k+1]));
     501      665734 :       if (gc_needed(av,2))
     502             :       {
     503           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate, i = %ld/%ld", i,n);
     504           0 :         Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     505             :       }
     506             :     }
     507             :   }
     508      376112 :   return gerepilecopy(av, Q);
     509             : }
     510             : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
     511             : GEN
     512        5866 : RgXQX_translate(GEN P, GEN c, GEN T)
     513             : {
     514        5866 :   pari_sp av = avma;
     515             :   GEN Q, *R;
     516             :   long i, k, n;
     517             : 
     518        5866 :   if (!signe(P) || gequal0(c)) return RgX_copy(P);
     519        5866 :   Q = leafcopy(P);
     520        5866 :   R = (GEN*)(Q+2); n = degpol(P);
     521       34181 :   for (i=1; i<=n; i++)
     522             :   {
     523      140077 :     for (k=n-i; k<n; k++)
     524             :     {
     525      111762 :       pari_sp av2 = avma;
     526      111762 :       R[k] = gerepileupto(av2, RgX_rem(gadd(R[k], gmul(c, R[k+1])), T));
     527             :     }
     528       28315 :     if (gc_needed(av,2))
     529             :     {
     530           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
     531           0 :       Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
     532             :     }
     533             :   }
     534        5866 :   return gerepilecopy(av, Q);
     535             : }
     536             : 
     537             : /********************************************************************/
     538             : /**                                                                **/
     539             : /**                          CONVERSIONS                           **/
     540             : /**                       (not memory clean)                       **/
     541             : /**                                                                **/
     542             : /********************************************************************/
     543             : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
     544             :  * but everything else is */
     545             : static GEN
     546       14648 : QXQ_to_mod_copy(GEN x, GEN T)
     547             : {
     548             :   long d;
     549       14648 :   switch(typ(x))
     550             :   {
     551        5201 :     case t_INT:  return icopy(x);
     552         399 :     case t_FRAC: return gcopy(x);
     553             :     case t_POL:
     554        9048 :       d = degpol(x);
     555        9048 :       if (d < 0) return gen_0;
     556        8768 :       if (d == 0) return gcopy(gel(x,2));
     557        8600 :       return mkpolmod(RgX_copy(x), T);
     558           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     559             :              return NULL;/* LCOV_EXCL_LINE */
     560             :   }
     561             : }
     562             : /* pure shallow version */
     563             : static GEN
     564      408335 : QXQ_to_mod(GEN x, GEN T)
     565             : {
     566             :   long d;
     567      408335 :   switch(typ(x))
     568             :   {
     569             :     case t_INT:
     570      355201 :     case t_FRAC: return x;
     571             :     case t_POL:
     572       53134 :       d = degpol(x);
     573       53134 :       if (d < 0) return gen_0;
     574       48902 :       if (d == 0) return gel(x,2);
     575       45045 :       return mkpolmod(x, T);
     576           0 :     default: pari_err_TYPE("QXQ_to_mod",x);
     577             :              return NULL;/* LCOV_EXCL_LINE */
     578             :   }
     579             : }
     580             : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod_copy to all coeffs.
     581             :  * Not memory clean because T not copied, but everything else is */
     582             : static GEN
     583        1848 : QXQX_to_mod(GEN z, GEN T)
     584             : {
     585        1848 :   long i,l = lg(z);
     586        1848 :   GEN x = cgetg(l,t_POL);
     587        1848 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_copy(gel(z,i), T);
     588        1848 :   x[1] = z[1]; return normalizepol_lg(x,l);
     589             : }
     590             : /* pure shallow version */
     591             : GEN
     592       83118 : QXQX_to_mod_shallow(GEN z, GEN T)
     593             : {
     594       83118 :   long i,l = lg(z);
     595       83118 :   GEN x = cgetg(l,t_POL);
     596       83118 :   for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
     597       83118 :   x[1] = z[1]; return normalizepol_lg(x,l);
     598             : }
     599             : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
     600             : GEN
     601         546 : QXQXV_to_mod(GEN V, GEN T)
     602             : {
     603         546 :   long i, l = lg(V);
     604         546 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     605         546 :   for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
     606         546 :   return z;
     607             : }
     608             : /* Apply QXQ_to_mod_copy to all entries. Memory-clean ! */
     609             : GEN
     610        1464 : QXQV_to_mod(GEN V, GEN T)
     611             : {
     612        1464 :   long i, l = lg(V);
     613        1464 :   GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
     614        1464 :   for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_copy(gel(V,i), T);
     615        1464 :   return z;
     616             : }
     617             : 
     618             : GEN
     619      720291 : RgX_renormalize_lg(GEN x, long lx)
     620             : {
     621             :   long i;
     622     1981265 :   for (i = lx-1; i>1; i--)
     623     1866094 :     if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
     624      720291 :   stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
     625      720291 :   setlg(x, i+1); setsigne(x, i != 1); return x;
     626             : }
     627             : 
     628             : GEN
     629      481330 : RgV_to_RgX(GEN x, long v)
     630             : {
     631      481330 :   long i, k = lg(x);
     632             :   GEN p;
     633             : 
     634      481330 :   while (--k && gequal0(gel(x,k)));
     635      481332 :   if (!k) return pol_0(v);
     636      463111 :   i = k+2; p = cgetg(i,t_POL);
     637      463112 :   p[1] = evalsigne(1) | evalvarn(v);
     638      463112 :   x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
     639      463112 :   return p;
     640             : }
     641             : GEN
     642      149566 : RgV_to_RgX_reverse(GEN x, long v)
     643             : {
     644      149566 :   long j, k, l = lg(x);
     645             :   GEN p;
     646             : 
     647      149566 :   for (k = 1; k < l; k++)
     648      149566 :     if (!gequal0(gel(x,k))) break;
     649      149566 :   if (k == l) return pol_0(v);
     650      149566 :   k -= 1;
     651      149566 :   l -= k;
     652      149566 :   x += k;
     653      149566 :   p = cgetg(l+1,t_POL);
     654      149566 :   p[1] = evalsigne(1) | evalvarn(v);
     655      149566 :   for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
     656      149566 :   return p;
     657             : }
     658             : 
     659             : /* return the (N-dimensional) vector of coeffs of p */
     660             : GEN
     661     4079650 : RgX_to_RgC(GEN x, long N)
     662             : {
     663             :   long i, l;
     664             :   GEN z;
     665     4079650 :   l = lg(x)-1; x++;
     666     4079650 :   if (l > N+1) l = N+1; /* truncate higher degree terms */
     667     4079650 :   z = cgetg(N+1,t_COL);
     668     4079650 :   for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
     669     4079650 :   for (   ; i<=N; i++) gel(z,i) = gen_0;
     670     4079650 :   return z;
     671             : }
     672             : GEN
     673      679022 : Rg_to_RgC(GEN x, long N)
     674             : {
     675      679022 :   return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
     676             : }
     677             : 
     678             : /* vector of polynomials (in v) whose coeffs are given by the columns of x */
     679             : GEN
     680       47724 : RgM_to_RgXV(GEN x, long v)
     681             : {
     682       47724 :   long j, lx = lg(x);
     683       47724 :   GEN y = cgetg(lx, t_VEC);
     684       47724 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), v);
     685       47724 :   return y;
     686             : }
     687             : 
     688             : /* matrix whose entries are given by the coeffs of the polynomials in
     689             :  * vector v (considered as degree n-1 polynomials) */
     690             : GEN
     691      103050 : RgV_to_RgM(GEN v, long n)
     692             : {
     693      103050 :   long j, N = lg(v);
     694      103050 :   GEN y = cgetg(N, t_MAT);
     695      103050 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j), n);
     696      103050 :   return y;
     697             : }
     698             : GEN
     699        2247 : RgXV_to_RgM(GEN v, long n)
     700             : {
     701        2247 :   long j, N = lg(v);
     702        2247 :   GEN y = cgetg(N, t_MAT);
     703        2247 :   for (j=1; j<N; j++) gel(y,j) = RgX_to_RgC(gel(v,j), n);
     704        2247 :   return y;
     705             : }
     706             : 
     707             : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
     708             : GEN
     709       15204 : RgM_to_RgXX(GEN x, long v,long w)
     710             : {
     711       15204 :   long j, lx = lg(x);
     712       15204 :   GEN y = cgetg(lx+1, t_POL);
     713       15204 :   y[1] = evalsigne(1) | evalvarn(v);
     714       15204 :   y++;
     715       15204 :   for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
     716       15204 :   return normalizepol_lg(--y, lx+1);
     717             : }
     718             : 
     719             : /* matrix whose entries are given by the coeffs of the polynomial v in
     720             :  * two variables (considered as degree n-1 polynomials) */
     721             : GEN
     722          21 : RgXX_to_RgM(GEN v, long n)
     723             : {
     724          21 :   long j, N = lg(v)-1;
     725          21 :   GEN y = cgetg(N, t_MAT);
     726          21 :   for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
     727          21 :   return y;
     728             : }
     729             : 
     730             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     731             : GEN
     732       11554 : RgXY_swapspec(GEN x, long n, long w, long nx)
     733             : {
     734       11554 :   long j, ly = n+3;
     735       11554 :   GEN y = cgetg(ly, t_POL);
     736       11554 :   y[1] = evalsigne(1);
     737      183010 :   for (j=2; j<ly; j++)
     738             :   {
     739             :     long k;
     740      171456 :     GEN a = cgetg(nx+2,t_POL);
     741      171456 :     a[1] = evalsigne(1) | evalvarn(w);
     742      931798 :     for (k=0; k<nx; k++)
     743             :     {
     744      760342 :       GEN xk = gel(x,k);
     745      760342 :       if (typ(xk)==t_POL)
     746      674767 :         gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
     747             :       else
     748       85575 :         gel(a,k+2) = j==2 ? xk: gen_0;
     749             :     }
     750      171456 :     gel(y,j) = normalizepol_lg(a, nx+2);
     751             :   }
     752       11554 :   return normalizepol_lg(y,ly);
     753             : }
     754             : 
     755             : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
     756             : GEN
     757         224 : RgXY_swap(GEN x, long n, long w)
     758             : {
     759         224 :   GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
     760         224 :   setvarn(z, varn(x)); return z;
     761             : }
     762             : 
     763             : long
     764           1 : RgXY_degreex(GEN b)
     765             : {
     766           1 :   long deg = -1, i;
     767           1 :   if (!signe(b)) return -1;
     768           3 :   for (i = 2; i < lg(b); ++i)
     769             :   {
     770           2 :     GEN bi = gel(b, i);
     771           2 :     if (typ(bi) == t_POL)
     772           1 :       deg = maxss(deg, degpol(bi));
     773             :   }
     774           1 :   return deg;
     775             : }
     776             : 
     777             : /* return (x % X^n). Shallow */
     778             : GEN
     779       81013 : RgXn_red_shallow(GEN a, long n)
     780             : {
     781       81013 :   long i, L = n+2, l = lg(a);
     782             :   GEN  b;
     783       81013 :   if (L >= l) return a; /* deg(x) < n */
     784       44998 :   b = cgetg(L, t_POL); b[1] = a[1];
     785       44998 :   for (i=2; i<L; i++) gel(b,i) = gel(a,i);
     786       44998 :   return normalizepol_lg(b,L);
     787             : }
     788             : 
     789             : GEN
     790         336 : RgXnV_red_shallow(GEN P, long n)
     791             : {
     792         336 :   long i, l = lg(P);
     793         336 :   GEN Q = cgetg(l, t_VEC);
     794         336 :   for (i=1; i<l; i++) gel(Q,i) = RgXn_red_shallow(gel(P,i), n);
     795         336 :   return Q;
     796             : }
     797             : 
     798             : /* return (x * X^n). Shallow */
     799             : GEN
     800    57140113 : RgX_shift_shallow(GEN a, long n)
     801             : {
     802    57140113 :   long i, l = lg(a);
     803             :   GEN  b;
     804    57140113 :   if (l == 2 || !n) return a;
     805    42257295 :   l += n;
     806    42257295 :   if (n < 0)
     807             :   {
     808    37521026 :     if (l <= 2) return pol_0(varn(a));
     809    37502518 :     b = cgetg(l, t_POL); b[1] = a[1];
     810    37502519 :     a -= n;
     811    37502519 :     for (i=2; i<l; i++) gel(b,i) = gel(a,i);
     812             :   } else {
     813     4736269 :     b = cgetg(l, t_POL); b[1] = a[1];
     814     4736268 :     a -= n; n += 2;
     815     4736268 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     816     4736268 :     for (   ; i<l; i++) gel(b,i) = gel(a,i);
     817             :   }
     818    42238787 :   return b;
     819             : }
     820             : /* return (x * X^n). */
     821             : GEN
     822     3393242 : RgX_shift(GEN a, long n)
     823             : {
     824     3393242 :   long i, l = lg(a);
     825             :   GEN  b;
     826     3393242 :   if (l == 2 || !n) return RgX_copy(a);
     827     3393011 :   l += n;
     828     3393011 :   if (n < 0)
     829             :   {
     830         595 :     if (l <= 2) return pol_0(varn(a));
     831         553 :     b = cgetg(l, t_POL); b[1] = a[1];
     832         553 :     a -= n;
     833         553 :     for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     834             :   } else {
     835     3392416 :     b = cgetg(l, t_POL); b[1] = a[1];
     836     3392416 :     a -= n; n += 2;
     837     3392416 :     for (i=2; i<n; i++) gel(b,i) = gen_0;
     838     3392416 :     for (   ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
     839             :   }
     840     3392969 :   return b;
     841             : }
     842             : 
     843             : GEN
     844      315273 : RgX_rotate_shallow(GEN P, long k, long p)
     845             : {
     846      315273 :   long i, l = lgpol(P);
     847             :   GEN r;
     848      315273 :   if (signe(P)==0)
     849         427 :     return pol_0(varn(P));
     850      314846 :   r = cgetg(p+2,t_POL); r[1] = P[1];
     851     2091390 :   for(i=0; i<p; i++)
     852             :   {
     853     1776544 :     long s = 2+(i+k)%p;
     854     1776544 :     gel(r,s) = i<l? gel(P,2+i): gen_0;
     855             :   }
     856      314846 :   return RgX_renormalize(r);
     857             : }
     858             : 
     859             : GEN
     860     2857944 : RgX_mulXn(GEN x, long d)
     861             : {
     862             :   pari_sp av;
     863             :   GEN z;
     864             :   long v;
     865     2857944 :   if (d >= 0) return RgX_shift(x, d);
     866     1364589 :   d = -d;
     867     1364589 :   v = RgX_val(x);
     868     1364589 :   if (v >= d) return RgX_shift(x, -d);
     869     1364582 :   av = avma;
     870     1364582 :   z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
     871     1364582 :   return gerepileupto(av, z);
     872             : }
     873             : 
     874             : long
     875     2092892 : RgX_val(GEN x)
     876             : {
     877     2092892 :   long i, lx = lg(x);
     878     2092892 :   if (lx == 2) return LONG_MAX;
     879     2118281 :   for (i = 2; i < lx; i++)
     880     2118267 :     if (!isexactzero(gel(x,i))) break;
     881     2092850 :   if (i == lx) return LONG_MAX;/* possible with non-rational zeros */
     882     2092836 :   return i - 2;
     883             : }
     884             : long
     885    42140073 : RgX_valrem(GEN x, GEN *Z)
     886             : {
     887    42140073 :   long v, i, lx = lg(x);
     888    42140073 :   if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
     889    81417220 :   for (i = 2; i < lx; i++)
     890    81417220 :     if (!isexactzero(gel(x,i))) break;
     891             :   /* possible with non-rational zeros */
     892    42140073 :   if (i == lx) { *Z = pol_0(varn(x)); return LONG_MAX; }
     893    42140073 :   v = i - 2;
     894    42140073 :   *Z = RgX_shift_shallow(x, -v);
     895    42140073 :   return v;
     896             : }
     897             : long
     898        6883 : RgX_valrem_inexact(GEN x, GEN *Z)
     899             : {
     900             :   long v;
     901        6883 :   if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
     902        7170 :   for (v = 0;; v++)
     903        7170 :     if (!gequal0(gel(x,2+v))) break;
     904         294 :   if (Z) *Z = RgX_shift_shallow(x, -v);
     905        6876 :   return v;
     906             : }
     907             : 
     908             : GEN
     909           0 : RgXQC_red(GEN P, GEN T)
     910             : {
     911           0 :   long i, l = lg(P);
     912           0 :   GEN Q = cgetg(l, t_COL);
     913           0 :   for (i=1; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
     914           0 :   return Q;
     915             : }
     916             : 
     917             : GEN
     918          56 : RgXQV_red(GEN P, GEN T)
     919             : {
     920          56 :   long i, l = lg(P);
     921          56 :   GEN Q = cgetg(l, t_VEC);
     922          56 :   for (i=1; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
     923          56 :   return Q;
     924             : }
     925             : 
     926             : GEN
     927           0 : RgXQM_red(GEN P, GEN T)
     928             : {
     929           0 :   long i, l = lg(P);
     930           0 :   GEN Q = cgetg(l, t_MAT);
     931           0 :   for (i=1; i<l; i++) gel(Q,i) = RgXQC_red(gel(P,i), T);
     932           0 :   return Q;
     933             : }
     934             : 
     935             : GEN
     936           0 : RgXQM_mul(GEN P, GEN Q, GEN T)
     937             : {
     938           0 :   return RgXQM_red(RgM_mul(P, Q), T);
     939             : }
     940             : 
     941             : GEN
     942        3927 : RgXQX_red(GEN P, GEN T)
     943             : {
     944        3927 :   long i, l = lg(P);
     945        3927 :   GEN Q = cgetg(l, t_POL);
     946        3927 :   Q[1] = P[1];
     947        3927 :   for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
     948        3927 :   return normalizepol_lg(Q, l);
     949             : }
     950             : 
     951             : GEN
     952      210819 : RgX_deriv(GEN x)
     953             : {
     954      210819 :   long i,lx = lg(x)-1;
     955             :   GEN y;
     956             : 
     957      210819 :   if (lx<3) return pol_0(varn(x));
     958      184142 :   y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
     959      184142 :   for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
     960      184142 :   y[1] = x[1]; return normalizepol_lg(y,i);
     961             : }
     962             : 
     963             : GEN
     964      295553 : RgX_recipspec_shallow(GEN x, long l, long n)
     965             : {
     966             :   long i;
     967      295553 :   GEN z=cgetg(n+2,t_POL)+2;
     968    14672111 :   for(i=0; i<l; i++)
     969    14376555 :     gel(z,n-i-1) = gel(x,i);
     970      383280 :   for(   ; i<n; i++)
     971       87724 :     gel(z, n-i-1) = gen_0;
     972      295556 :   return normalizepol_lg(z-2,n+2);
     973             : }
     974             : 
     975             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
     976             : GEN
     977         504 : RgX_recip(GEN x)
     978             : {
     979             :   long lx, i, j;
     980         504 :   GEN y = cgetg_copy(x, &lx);
     981         504 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
     982         504 :   return normalizepol_lg(y,lx);
     983             : }
     984             : /* shallow version */
     985             : GEN
     986      406891 : RgX_recip_shallow(GEN x)
     987             : {
     988             :   long lx, i, j;
     989      406891 :   GEN y = cgetg_copy(x, &lx);
     990      406893 :   y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
     991      406893 :   return y;
     992             : }
     993             : /*******************************************************************/
     994             : /*                                                                 */
     995             : /*                      ADDITION / SUBTRACTION                     */
     996             : /*                                                                 */
     997             : /*******************************************************************/
     998             : /* same variable */
     999             : GEN
    1000    26959144 : RgX_add(GEN x, GEN y)
    1001             : {
    1002    26959144 :   long i, lx = lg(x), ly = lg(y);
    1003             :   GEN z;
    1004    26959144 :   if (ly <= lx) {
    1005    25023145 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1006    25023153 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1007    25023145 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1008    25023145 :     z = normalizepol_lg(z, lx);
    1009             :   } else {
    1010     1935999 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1011     1936011 :     for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
    1012     1936002 :     for (   ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
    1013     1936002 :     z = normalizepol_lg(z, ly);
    1014             :   }
    1015    26959147 :   return z;
    1016             : }
    1017             : GEN
    1018    10862082 : RgX_sub(GEN x, GEN y)
    1019             : {
    1020    10862082 :   long i, lx = lg(x), ly = lg(y);
    1021             :   GEN z;
    1022    10862082 :   if (ly <= lx) {
    1023     8532100 :     z = cgetg(lx,t_POL); z[1] = x[1];
    1024     8532141 :     for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1025     8532095 :     for (   ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
    1026     8532095 :     z = normalizepol_lg(z, lx);
    1027             :   } else {
    1028     2329982 :     z = cgetg(ly,t_POL); z[1] = y[1];
    1029     2329982 :     for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
    1030     2329982 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    1031     2329982 :     z = normalizepol_lg(z, ly);
    1032             :   }
    1033    10862079 :   return z;
    1034             : }
    1035             : GEN
    1036      993469 : RgX_neg(GEN x)
    1037             : {
    1038      993469 :   long i, lx = lg(x);
    1039      993469 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    1040      993469 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    1041      993469 :   return y;
    1042             : }
    1043             : 
    1044             : GEN
    1045    11556429 : RgX_Rg_add(GEN y, GEN x)
    1046             : {
    1047             :   GEN z;
    1048    11556429 :   long lz = lg(y), i;
    1049    11556429 :   if (lz == 2) return scalarpol(x,varn(y));
    1050     9863936 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1051     9863936 :   gel(z,2) = gadd(gel(y,2),x);
    1052     9863936 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1053             :   /* probably useless unless lz = 3, but cannot be skipped if y is
    1054             :    * an inexact 0 */
    1055     9863936 :   return normalizepol_lg(z,lz);
    1056             : }
    1057             : GEN
    1058        2422 : RgX_Rg_add_shallow(GEN y, GEN x)
    1059             : {
    1060             :   GEN z;
    1061        2422 :   long lz = lg(y), i;
    1062        2422 :   if (lz == 2) return scalarpol(x,varn(y));
    1063        2422 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1064        2422 :   gel(z,2) = gadd(gel(y,2),x);
    1065        2422 :   for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
    1066        2422 :   return z = normalizepol_lg(z,lz);
    1067             : }
    1068             : GEN
    1069       32488 : RgX_Rg_sub(GEN y, GEN x)
    1070             : {
    1071             :   GEN z;
    1072       32488 :   long lz = lg(y), i;
    1073       32488 :   if (lz == 2)
    1074             :   { /* scalarpol(gneg(x),varn(y)) optimized */
    1075        3864 :     long v = varn(y);
    1076        3864 :     if (isrationalzero(x)) return pol_0(v);
    1077          14 :     z = cgetg(3,t_POL);
    1078          28 :     z[1] = gequal0(x)? evalvarn(v)
    1079          14 :                    : evalvarn(v) | evalsigne(1);
    1080          14 :     gel(z,2) = gneg(x); return z;
    1081             :   }
    1082       28624 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1083       28624 :   gel(z,2) = gsub(gel(y,2),x);
    1084       28624 :   for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
    1085       28624 :   return z = normalizepol_lg(z,lz);
    1086             : }
    1087             : GEN
    1088      380774 : Rg_RgX_sub(GEN x, GEN y)
    1089             : {
    1090             :   GEN z;
    1091      380774 :   long lz = lg(y), i;
    1092      380774 :   if (lz == 2) return scalarpol(x,varn(y));
    1093      379759 :   z = cgetg(lz,t_POL); z[1] = y[1];
    1094      379759 :   gel(z,2) = gsub(x, gel(y,2));
    1095      379759 :   for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
    1096      379759 :   return z = normalizepol_lg(z,lz);
    1097             : }
    1098             : /*******************************************************************/
    1099             : /*                                                                 */
    1100             : /*                  KARATSUBA MULTIPLICATION                       */
    1101             : /*                                                                 */
    1102             : /*******************************************************************/
    1103             : #if 0
    1104             : /* to debug Karatsuba-like routines */
    1105             : GEN
    1106             : zx_debug_spec(GEN x, long nx)
    1107             : {
    1108             :   GEN z = cgetg(nx+2,t_POL);
    1109             :   long i;
    1110             :   for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
    1111             :   z[1] = evalsigne(1); return z;
    1112             : }
    1113             : 
    1114             : GEN
    1115             : RgX_debug_spec(GEN x, long nx)
    1116             : {
    1117             :   GEN z = cgetg(nx+2,t_POL);
    1118             :   long i;
    1119             :   for (i=0; i<nx; i++) z[i+2] = x[i];
    1120             :   z[1] = evalsigne(1); return z;
    1121             : }
    1122             : #endif
    1123             : 
    1124             : /* generic multiplication */
    1125             : GEN
    1126     3156023 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
    1127             : {
    1128             :   GEN z, t;
    1129             :   long i;
    1130     3156023 :   if (nx == ny) {
    1131      748649 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1132      748666 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1133      748646 :     return normalizepol_lg(z, nx+2);
    1134             :   }
    1135     2407374 :   if (ny < nx) {
    1136     2106538 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1137     2106547 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1138     2106537 :     for (   ; i < nx; i++) gel(t,i) = gel(x,i);
    1139     2106537 :     return normalizepol_lg(z, nx+2);
    1140             :   } else {
    1141      300836 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1142      300858 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1143      300834 :     for (   ; i < ny; i++) gel(t,i) = gel(y,i);
    1144      300834 :     return normalizepol_lg(z, ny+2);
    1145             :   }
    1146             : }
    1147             : GEN
    1148      371290 : RgX_addspec(GEN x, GEN y, long nx, long ny)
    1149             : {
    1150             :   GEN z, t;
    1151             :   long i;
    1152      371290 :   if (nx == ny) {
    1153         601 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1154         601 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1155         601 :     return normalizepol_lg(z, nx+2);
    1156             :   }
    1157      370689 :   if (ny < nx) {
    1158      370689 :     z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1159      370709 :     for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1160      370691 :     for (   ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
    1161      370692 :     return normalizepol_lg(z, nx+2);
    1162             :   } else {
    1163           0 :     z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
    1164           0 :     for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
    1165           0 :     for (   ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
    1166           0 :     return normalizepol_lg(z, ny+2);
    1167             :   }
    1168             : }
    1169             : 
    1170             : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
    1171             :  * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
    1172             :  * t_VECSMALL, to hold this, which can be left on stack: gerepile
    1173             :  * will not crash on it. The returned vector itself is not a proper GEN,
    1174             :  * we access the coefficients as x[i], i = 0..deg(x) */
    1175             : static GEN
    1176    32163893 : RgXspec_kill0(GEN x, long lx)
    1177             : {
    1178    32163893 :   GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
    1179             :   long i;
    1180   136711673 :   for (i=0; i <lx; i++)
    1181             :   {
    1182   104547871 :     GEN c = gel(x,i);
    1183   104547871 :     z[i] = (long)(isrationalzero(c)? NULL: c);
    1184             :   }
    1185    32163802 :   return z;
    1186             : }
    1187             : 
    1188             : INLINE GEN
    1189    74727612 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
    1190             : {
    1191    74727612 :   pari_sp av = avma;
    1192    74727612 :   GEN s = NULL;
    1193             :   long i;
    1194             : 
    1195   298648865 :   for (i=a; i<b; i++)
    1196   223930713 :     if (gel(y,i) && gel(x,-i))
    1197             :     {
    1198   171915562 :       GEN t = gmul(gel(y,i), gel(x,-i));
    1199   171925383 :       s = s? gadd(s, t): t;
    1200             :     }
    1201    74718152 :   return s? gerepileupto(av, s): gen_0;
    1202             : }
    1203             : 
    1204             : /* assume nx >= ny > 0, return x * y * t^v */
    1205             : static GEN
    1206    12811713 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
    1207             : {
    1208             :   long i, lz, nz;
    1209             :   GEN z;
    1210             : 
    1211    12811713 :   x = RgXspec_kill0(x,nx);
    1212    12811687 :   y = RgXspec_kill0(y,ny);
    1213    12811697 :   lz = nx + ny + 1; nz = lz-2;
    1214    12811697 :   lz += v;
    1215    12811697 :   z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
    1216    12811808 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1217    12811808 :   for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
    1218    12811606 :   for (  ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
    1219    12811597 :   for (  ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
    1220    12811700 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1221             : }
    1222             : 
    1223             : /* return (x * X^d) + y. Assume d > 0 */
    1224             : GEN
    1225     1908622 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
    1226             : {
    1227             :   GEN x, y, xd, yd, zd;
    1228             :   long a, lz, nx, ny;
    1229             : 
    1230     1908622 :   if (!signe(x0)) return y0;
    1231     1887711 :   ny = lgpol(y0);
    1232     1887711 :   nx = lgpol(x0);
    1233     1887710 :   zd = (GEN)avma;
    1234     1887710 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1235     1887710 :   if (a <= 0)
    1236             :   {
    1237      159443 :     lz = nx+d+2;
    1238      159443 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1239      159447 :     while (xd > x) gel(--zd,0) = gel(--xd,0);
    1240      159447 :     x = zd + a;
    1241      159447 :     while (zd > x) gel(--zd,0) = gen_0;
    1242             :   }
    1243             :   else
    1244             :   {
    1245     1728267 :     xd = new_chunk(d); yd = y+d;
    1246     1728267 :     x = RgX_addspec_shallow(x,yd, nx,a);
    1247     1728267 :     lz = (a>nx)? ny+2: lg(x)+d;
    1248     1728267 :     x += 2; while (xd > x) *--zd = *--xd;
    1249             :   }
    1250     1887714 :   while (yd > y) *--zd = *--yd;
    1251     1887714 :   *--zd = x0[1];
    1252     1887714 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1253             : }
    1254             : GEN
    1255      665295 : RgX_addmulXn(GEN x0, GEN y0, long d)
    1256             : {
    1257             :   GEN x, y, xd, yd, zd;
    1258             :   long a, lz, nx, ny;
    1259             : 
    1260      665295 :   if (!signe(x0)) return RgX_copy(y0);
    1261      665002 :   nx = lgpol(x0);
    1262      665001 :   ny = lgpol(y0);
    1263      665000 :   zd = (GEN)avma;
    1264      665000 :   x = x0 + 2; y = y0 + 2; a = ny-d;
    1265      665000 :   if (a <= 0)
    1266             :   {
    1267      293711 :     lz = nx+d+2;
    1268      293711 :     (void)new_chunk(lz); xd = x+nx; yd = y+ny;
    1269      293727 :     while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
    1270      293727 :     x = zd + a;
    1271      293727 :     while (zd > x) gel(--zd,0) = gen_0;
    1272             :   }
    1273             :   else
    1274             :   {
    1275      371289 :     xd = new_chunk(d); yd = y+d;
    1276      371290 :     x = RgX_addspec(x,yd, nx,a);
    1277      371293 :     lz = (a>nx)? ny+2: lg(x)+d;
    1278      371293 :     x += 2; while (xd > x) *--zd = *--xd;
    1279             :   }
    1280      665020 :   while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
    1281      665003 :   *--zd = x0[1];
    1282      665003 :   *--zd = evaltyp(t_POL) | evallg(lz); return zd;
    1283             : }
    1284             : 
    1285             : /* return x * y mod t^n */
    1286             : static GEN
    1287     3105709 : RgXn_mul_basecase(GEN x, GEN y, long n)
    1288             : {
    1289     3105709 :   long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
    1290             :   GEN z;
    1291     3105709 :   if (lx < 0) return pol_0(varn(x));
    1292     3105709 :   if (ly < 0) return pol_0(varn(x));
    1293     3105709 :   z = cgetg(lz, t_POL) + 2;
    1294     3105709 :   x+=2; if (lx > n) lx = n;
    1295     3105709 :   y+=2; if (ly > n) ly = n;
    1296     3105709 :   z[-1] = x[-1];
    1297     3105709 :   if (ly > lx) { swap(x,y); lswap(lx,ly); }
    1298     3105709 :   x = RgXspec_kill0(x, lx);
    1299     3105709 :   y = RgXspec_kill0(y, ly);
    1300             :   /* x:y:z [i] = term of degree i */
    1301     3105709 :   for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
    1302     3105709 :   for (  ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
    1303     3105709 :   for (  ; i<n; i++)  gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
    1304     3105709 :   return normalizepol_lg(z - 2, lz);
    1305             : }
    1306             : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
    1307             : GEN
    1308     3667697 : RgXn_mul(GEN f, GEN g, long n)
    1309             : {
    1310     3667697 :   pari_sp av = avma;
    1311             :   GEN fe,fo, ge,go, l,h,m;
    1312             :   long n0, n1;
    1313     3667697 :   if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
    1314     3108061 :   if (n < 80) return RgXn_mul_basecase(f,g,n);
    1315        2352 :   n0 = n>>1; n1 = n-n0;
    1316        2352 :   RgX_even_odd(f, &fe, &fo);
    1317        2352 :   RgX_even_odd(g, &ge, &go);
    1318        2352 :   l = RgXn_mul(fe,ge,n1);
    1319        2352 :   h = RgXn_mul(fo,go,n0);
    1320        2352 :   m = RgX_sub(RgXn_mul(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
    1321             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1322             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1323        2352 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1324             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1325        2352 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1326        2352 :   m = RgX_inflate(m,2);
    1327             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1328        2352 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1329        2352 :   h = RgX_inflate(h,2);
    1330        2352 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1331        2352 :   return gerepileupto(av, h);
    1332             : }
    1333             : /* (f*g) \/ x^n */
    1334             : GEN
    1335       33774 : RgX_mulhigh_i(GEN f, GEN g, long n)
    1336             : {
    1337       33774 :   long d = degpol(f)+degpol(g) + 1 - n;
    1338             :   GEN h;
    1339       33774 :   if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
    1340        1795 :   h = RgX_recip_shallow(RgXn_mul(RgX_recip_shallow(f),
    1341             :                                  RgX_recip_shallow(g), d));
    1342        1795 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1343             : }
    1344             : 
    1345             : /* (f*g) \/ x^n */
    1346             : GEN
    1347           0 : RgX_sqrhigh_i(GEN f, long n)
    1348             : {
    1349           0 :   long d = 2*degpol(f)+ 1 - n;
    1350             :   GEN h;
    1351           0 :   if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
    1352           0 :   h = RgX_recip_shallow(RgXn_sqr(RgX_recip_shallow(f), d));
    1353           0 :   return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
    1354             : }
    1355             : 
    1356             : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
    1357             :  * b+2 were sent instead. na, nb = number of terms of a, b.
    1358             :  * Only c, c0, c1, c2 are genuine GEN.
    1359             :  */
    1360             : GEN
    1361    13681840 : RgX_mulspec(GEN a, GEN b, long na, long nb)
    1362             : {
    1363             :   GEN a0, c, c0;
    1364    13681840 :   long n0, n0a, i, v = 0;
    1365             :   pari_sp av;
    1366             : 
    1367    13681840 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
    1368    13681838 :   while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
    1369    13681838 :   if (na < nb) swapspec(a,b, na,nb);
    1370    13681838 :   if (!nb) return pol_0(0);
    1371             : 
    1372    13473835 :   if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
    1373      662110 :   RgX_shift_inplace_init(v);
    1374      662111 :   i = (na>>1); n0 = na-i; na = i;
    1375      662111 :   av = avma; a0 = a+n0; n0a = n0;
    1376      662111 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1377             : 
    1378      662110 :   if (nb > n0)
    1379             :   {
    1380             :     GEN b0,c1,c2;
    1381             :     long n0b;
    1382             : 
    1383      653359 :     nb -= n0; b0 = b+n0; n0b = n0;
    1384      653359 :     while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
    1385      653359 :     c = RgX_mulspec(a,b,n0a,n0b);
    1386      653357 :     c0 = RgX_mulspec(a0,b0, na,nb);
    1387             : 
    1388      653358 :     c2 = RgX_addspec_shallow(a0,a, na,n0a);
    1389      653360 :     c1 = RgX_addspec_shallow(b0,b, nb,n0b);
    1390             : 
    1391      653357 :     c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
    1392      653357 :     c2 = RgX_sub(c1, RgX_add(c0,c));
    1393      653358 :     c0 = RgX_addmulXn_shallow(c0, c2, n0);
    1394             :   }
    1395             :   else
    1396             :   {
    1397        8751 :     c = RgX_mulspec(a,b,n0a,nb);
    1398        8751 :     c0 = RgX_mulspec(a0,b,na,nb);
    1399             :   }
    1400      662110 :   c0 = RgX_addmulXn(c0,c,n0);
    1401      662111 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1402             : }
    1403             : 
    1404             : INLINE GEN
    1405     4142282 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
    1406             : {
    1407     4142282 :   pari_sp av = avma;
    1408     4142282 :   GEN s = NULL;
    1409     4142282 :   long j, l = (i+1)>>1;
    1410    17908686 :   for (j=a; j<l; j++)
    1411             :   {
    1412    13769580 :     GEN xj = gel(x,j), xx = gel(x,i-j);
    1413    13769580 :     if (xj && xx)
    1414             :     {
    1415     8527962 :       GEN t = gmul(xj, xx);
    1416     8531717 :       s = s? gadd(s, t): t;
    1417             :     }
    1418             :   }
    1419     4139106 :   if (s) s = gshift(s,1);
    1420     4139219 :   if ((i&1) == 0)
    1421             :   {
    1422     2234603 :     GEN t = gel(x, i>>1);
    1423     2234603 :     if (t) {
    1424     1843372 :       t = gsqr(t);
    1425     1843305 :       s = s? gadd(s, t): t;
    1426             :     }
    1427             :   }
    1428     4138871 :   return s? gerepileupto(av,s): gen_0;
    1429             : }
    1430             : static GEN
    1431      328420 : RgX_sqrspec_basecase(GEN x, long nx, long v)
    1432             : {
    1433             :   long i, lz, nz;
    1434             :   GEN z;
    1435             : 
    1436      328420 :   if (!nx) return pol_0(0);
    1437      328413 :   x = RgXspec_kill0(x,nx);
    1438      328425 :   lz = (nx << 1) + 1, nz = lz-2;
    1439      328425 :   lz += v;
    1440      328425 :   z = cgetg(lz,t_POL) + 2;
    1441      328529 :   for (i=0; i<v; i++) gel(z++, 0) = gen_0;
    1442      328529 :   for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1443      328401 :   for (  ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
    1444      328422 :   z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
    1445             : }
    1446             : /* return x^2 mod t^n */
    1447             : static GEN
    1448         686 : RgXn_sqr_basecase(GEN x, long n)
    1449             : {
    1450         686 :   long i, lz = n+2, lx = lgpol(x);
    1451             :   GEN z;
    1452         686 :   if (lx < 0) return pol_0(varn(x));
    1453         686 :   z = cgetg(lz, t_POL);
    1454         686 :   z[1] = x[1];
    1455         686 :   x+=2; if (lx > n) lx = n;
    1456         686 :   x = RgXspec_kill0(x,lx);
    1457         686 :   z+=2;/* x:z [i] = term of degree i */
    1458         686 :   for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
    1459         686 :   for (  ; i<n; i++)  gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
    1460         686 :   z -= 2; return normalizepol_lg(z, lz);
    1461             : }
    1462             : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
    1463             : GEN
    1464       15666 : RgXn_sqr(GEN f, long n)
    1465             : {
    1466       15666 :   pari_sp av = avma;
    1467             :   GEN fe,fo, l,h,m;
    1468             :   long n0, n1;
    1469       15666 :   if (2*degpol(f) < n) return RgX_sqr(f);
    1470         714 :   if (n < 80) return RgXn_sqr_basecase(f,n);
    1471          28 :   n0 = n>>1; n1 = n-n0;
    1472          28 :   RgX_even_odd(f, &fe, &fo);
    1473          28 :   l = RgXn_sqr(fe,n1);
    1474          28 :   h = RgXn_sqr(fo,n0);
    1475          28 :   m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
    1476             :   /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
    1477             :    * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
    1478          28 :   l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
    1479             :   /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
    1480          28 :   if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
    1481          28 :   m = RgX_inflate(m,2);
    1482             :   /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
    1483          28 :   if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
    1484          28 :   h = RgX_inflate(h,2);
    1485          28 :   h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
    1486          28 :   return gerepileupto(av, h);
    1487             : }
    1488             : 
    1489             : GEN
    1490      329236 : RgX_sqrspec(GEN a, long na)
    1491             : {
    1492             :   GEN a0, c, c0, c1;
    1493      329236 :   long n0, n0a, i, v = 0;
    1494             :   pari_sp av;
    1495             : 
    1496      329236 :   while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
    1497      329238 :   if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
    1498         805 :   RgX_shift_inplace_init(v);
    1499         805 :   i = (na>>1); n0 = na-i; na = i;
    1500         805 :   av = avma; a0 = a+n0; n0a = n0;
    1501         805 :   while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
    1502             : 
    1503         805 :   c = RgX_sqrspec(a,n0a);
    1504         805 :   c0 = RgX_sqrspec(a0,na);
    1505         805 :   c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
    1506         805 :   c0 = RgX_addmulXn_shallow(c0,c1, n0);
    1507         805 :   c0 = RgX_addmulXn(c0,c,n0);
    1508         805 :   return RgX_shift_inplace(gerepileupto(av,c0), v);
    1509             : }
    1510             : 
    1511             : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
    1512             : GEN
    1513      414604 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
    1514             : {
    1515      414604 :   GEN z = RgX_mul(A, B);
    1516      414604 :   if (a < b)
    1517        5222 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
    1518      409382 :   else if (a > b)
    1519      258125 :     z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
    1520             :   else
    1521      151257 :     z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
    1522      414604 :   return z;
    1523             : }
    1524             : 
    1525             : GEN
    1526    11703466 : RgX_mul(GEN x, GEN y)
    1527             : {
    1528    11703466 :   GEN z = RgX_mulspec(y+2, x+2, lgpol(y), lgpol(x));
    1529    11703461 :   setvarn(z,varn(x)); return z;
    1530             : }
    1531             : 
    1532             : GEN
    1533      327618 : RgX_sqr(GEN x)
    1534             : {
    1535      327618 :   GEN z = RgX_sqrspec(x+2, lgpol(x));
    1536      327619 :   setvarn(z,varn(x)); return z;
    1537             : }
    1538             : 
    1539             : /*******************************************************************/
    1540             : /*                                                                 */
    1541             : /*                               DIVISION                          */
    1542             : /*                                                                 */
    1543             : /*******************************************************************/
    1544             : GEN
    1545      534361 : RgX_Rg_divexact(GEN x, GEN y) {
    1546             :   long i, lx;
    1547             :   GEN z;
    1548      534361 :   if (typ(y) == t_INT && is_pm1(y))
    1549      101748 :     return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
    1550      432613 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    1551      432613 :   for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1552      432613 :   return z;
    1553             : }
    1554             : GEN
    1555    22535881 : RgX_Rg_div(GEN x, GEN y) {
    1556             :   long i, lx;
    1557    22535881 :   GEN z = cgetg_copy(x, &lx); z[1] = x[1];
    1558    22535881 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    1559    22535881 :   return normalizepol_lg(z, lx);
    1560             : }
    1561             : GEN
    1562        1526 : RgX_normalize(GEN x)
    1563             : {
    1564        1526 :   GEN d = NULL;
    1565        1526 :   long i, n = lg(x)-1;
    1566        1526 :   for (i = n; i > 1; i--)
    1567             :   {
    1568        1526 :     d = gel(x,i);
    1569        1526 :     if (!gequal0(d)) break;
    1570             :   }
    1571        1526 :   if (i == 1) return pol_0(varn(x));
    1572        1526 :   if (i == n && isint1(d)) return x;
    1573         315 :   return normalizepol_lg(RgX_Rg_div(x, d), i+1);
    1574             : }
    1575             : GEN
    1576        1925 : RgX_divs(GEN x, long y) {
    1577             :   long i, lx;
    1578        1925 :   GEN z = cgetg_copy(x, &lx); z[1] = x[1];
    1579        1925 :   for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),y);
    1580        1925 :   return normalizepol_lg(z, lx);
    1581             : }
    1582             : GEN
    1583       35101 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
    1584             : {
    1585       35101 :   long l = lg(a), i;
    1586       35101 :   GEN a0, z0, z = cgetg(l-1, t_POL);
    1587       35101 :   z[1] = a[1];
    1588       35101 :   a0 = a + l-1;
    1589       35101 :   z0 = z + l-2; *z0 = *a0--;
    1590      768370 :   for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
    1591             :   {
    1592      733269 :     GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
    1593      733269 :     gel(z0,0) = t;
    1594             :   }
    1595       35101 :   if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
    1596       35101 :   return z;
    1597             : }
    1598             : /* Polynomial division x / y:
    1599             :  *   if z = ONLY_REM  return remainder, otherwise return quotient
    1600             :  *   if z != NULL set *z to remainder
    1601             :  *   *z is the last object on stack (and thus can be disposed of with cgiv
    1602             :  *   instead of gerepile) */
    1603             : /* assume, typ(x) = typ(y) = t_POL, same variable */
    1604             : GEN
    1605    14411824 : RgX_divrem(GEN x, GEN y, GEN *pr)
    1606             : {
    1607             :   pari_sp avy, av, av1;
    1608             :   long dx,dy,dz,i,j,sx,lr;
    1609             :   GEN z,p1,p2,rem,y_lead,mod;
    1610             :   GEN (*f)(GEN,GEN);
    1611             : 
    1612    14411824 :   if (!signe(y)) pari_err_INV("RgX_divrem",y);
    1613             : 
    1614    14411824 :   dy = degpol(y);
    1615    14411747 :   y_lead = gel(y,dy+2);
    1616    14411747 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    1617             :   {
    1618           0 :     pari_warn(warner,"normalizing a polynomial with 0 leading term");
    1619           0 :     for (dy--; dy>=0; dy--)
    1620             :     {
    1621           0 :       y_lead = gel(y,dy+2);
    1622           0 :       if (!gequal0(y_lead)) break;
    1623             :     }
    1624             :   }
    1625    14411729 :   if (!dy) /* y is constant */
    1626             :   {
    1627       13451 :     if (pr == ONLY_REM) return pol_0(varn(x));
    1628       12842 :     z = RgX_Rg_div(x, y_lead);
    1629       12842 :     if (pr == ONLY_DIVIDES) return z;
    1630       12135 :     if (pr) *pr = pol_0(varn(x));
    1631       12135 :     return z;
    1632             :   }
    1633    14398278 :   dx = degpol(x);
    1634    14398242 :   if (dx < dy)
    1635             :   {
    1636     1349922 :     if (pr == ONLY_REM) return RgX_copy(x);
    1637      309268 :     if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
    1638      309247 :     z = pol_0(varn(x));
    1639      309247 :     if (pr) *pr = RgX_copy(x);
    1640      309247 :     return z;
    1641             :   }
    1642             : 
    1643             :   /* x,y in R[X], y non constant */
    1644    13048320 :   av = avma;
    1645    13048320 :   switch(typ(y_lead))
    1646             :   {
    1647             :     case t_REAL:
    1648           0 :       y_lead = ginv(y_lead);
    1649           0 :       f = gmul; mod = NULL;
    1650           0 :       break;
    1651             :     case t_INTMOD:
    1652        4619 :     case t_POLMOD: y_lead = ginv(y_lead);
    1653        4619 :       f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
    1654        4619 :       break;
    1655    13043701 :     default: if (gequal1(y_lead)) y_lead = NULL;
    1656    13043670 :       f = gdiv; mod = NULL;
    1657             :   }
    1658             : 
    1659    13048289 :   if (y_lead == NULL)
    1660    11275804 :     p2 = gel(x,dx+2);
    1661             :   else {
    1662             :     for(;;) {
    1663     1772485 :       p2 = f(gel(x,dx+2),y_lead);
    1664     1772535 :       p2 = simplify_shallow(p2);
    1665     1772535 :       if (!isexactzero(p2) || (--dx < 0)) break;
    1666           0 :     }
    1667     1772535 :     if (dx < dy) /* leading coeff of x was in fact zero */
    1668             :     {
    1669           0 :       if (pr == ONLY_DIVIDES) {
    1670           0 :         avma = av;
    1671           0 :         return (dx < 0)? pol_0(varn(x)) : NULL;
    1672             :       }
    1673           0 :       if (pr == ONLY_REM)
    1674             :       {
    1675           0 :         if (dx < 0)
    1676           0 :           return gerepilecopy(av, scalarpol(p2, varn(x)));
    1677             :         else
    1678             :         {
    1679             :           GEN t;
    1680           0 :           avma = av;
    1681           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1682           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1683           0 :           return t;
    1684             :         }
    1685             :       }
    1686           0 :       if (pr) /* cf ONLY_REM above */
    1687             :       {
    1688           0 :         if (dx < 0)
    1689             :         {
    1690           0 :           p2 = gclone(p2);
    1691           0 :           avma = av;
    1692           0 :           z = pol_0(varn(x));
    1693           0 :           x = scalarpol(p2, varn(x));
    1694           0 :           gunclone(p2);
    1695             :         }
    1696             :         else
    1697             :         {
    1698             :           GEN t;
    1699           0 :           avma = av;
    1700           0 :           z = pol_0(varn(x));
    1701           0 :           t = cgetg(dx + 3, t_POL); t[1] = x[1];
    1702           0 :           for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
    1703           0 :           x = t;
    1704             :         }
    1705           0 :         *pr = x;
    1706             :       }
    1707             :       else
    1708             :       {
    1709           0 :         avma = av;
    1710           0 :         z = pol_0(varn(x));
    1711             :       }
    1712           0 :       return z;
    1713             :     }
    1714             :   }
    1715             :   /* dx >= dy */
    1716    13048339 :   avy = avma;
    1717    13048339 :   dz = dx-dy;
    1718    13048339 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1719    13048169 :   x += 2;
    1720    13048169 :   z += 2;
    1721    13048169 :   y += 2;
    1722    13048169 :   gel(z,dz) = gcopy(p2);
    1723             : 
    1724    41537433 :   for (i=dx-1; i>=dy; i--)
    1725             :   {
    1726    28488964 :     av1=avma; p1=gel(x,i);
    1727    28488964 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1728    28451641 :     if (y_lead) p1 = simplify(f(p1,y_lead));
    1729             : 
    1730    28451641 :     if (isrationalzero(p1)) { avma=av1; p1 = gen_0; }
    1731             :     else
    1732    18341236 :       p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
    1733    28488296 :     gel(z,i-dy) = p1;
    1734             :   }
    1735    13048469 :   if (!pr) return gerepileupto(av,z-2);
    1736             : 
    1737     6638073 :   rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
    1738     7584588 :   for (sx=0; ; i--)
    1739             :   {
    1740     7584588 :     p1 = gel(x,i);
    1741             :     /* we always enter this loop at least once */
    1742     7584588 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1743     7582767 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1744     7582767 :     if (!gequal0(p1)) { sx = 1; break; } /* remainder is non-zero */
    1745     4071353 :     if (!isexactzero(p1)) break;
    1746     4062007 :     if (!i) break;
    1747      946585 :     avma=av1;
    1748      946585 :   }
    1749     6637799 :   if (pr == ONLY_DIVIDES)
    1750             :   {
    1751         966 :     if (sx) { avma=av; return NULL; }
    1752         959 :     avma = (pari_sp)rem;
    1753         959 :     return gerepileupto(av,z-2);
    1754             :   }
    1755     6636833 :   lr=i+3; rem -= lr;
    1756     6636833 :   if (avma==av1) { avma = (pari_sp)rem; p1 = gcopy(p1); }
    1757     6534658 :   else p1 = gerepileupto((pari_sp)rem,p1);
    1758     6637003 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1759     6636308 :   rem[1] = z[-1];
    1760     6636308 :   rem += 2;
    1761     6636308 :   gel(rem,i) = p1;
    1762    21324755 :   for (i--; i>=0; i--)
    1763             :   {
    1764    14687700 :     av1=avma; p1 = gel(x,i);
    1765    14687700 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1766    14656613 :     if (mod && avma==av1) p1 = gmul(p1,mod);
    1767    14680828 :     gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
    1768             :   }
    1769     6637055 :   rem -= 2;
    1770     6637055 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1771     6637061 :   if (pr == ONLY_REM) return gerepileupto(av,rem);
    1772     4020114 :   z -= 2;
    1773             :   {
    1774     4020114 :     GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
    1775     4020114 :     gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
    1776             :   }
    1777             : }
    1778             : 
    1779             : /* x and y in (R[Y]/T)[X]  (lifted), T in R[Y]. y preferably monic */
    1780             : GEN
    1781       19503 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
    1782             : {
    1783             :   long vx, dx, dy, dz, i, j, sx, lr;
    1784             :   pari_sp av0, av, tetpil;
    1785             :   GEN z,p1,rem,lead;
    1786             : 
    1787       19503 :   if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
    1788       19503 :   vx = varn(x);
    1789       19503 :   dx = degpol(x);
    1790       19503 :   dy = degpol(y);
    1791       19503 :   if (dx < dy)
    1792             :   {
    1793        1470 :     if (pr)
    1794             :     {
    1795        1470 :       av0 = avma; x = RgXQX_red(x, T);
    1796        1470 :       if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gen_0; }
    1797        1470 :       if (pr == ONLY_REM) return x;
    1798           0 :       *pr = x;
    1799             :     }
    1800           0 :     return pol_0(vx);
    1801             :   }
    1802       18033 :   lead = leading_coeff(y);
    1803       18033 :   if (!dy) /* y is constant */
    1804             :   {
    1805           7 :     if (pr && pr != ONLY_DIVIDES)
    1806             :     {
    1807           0 :       if (pr == ONLY_REM) return pol_0(vx);
    1808           0 :       *pr = pol_0(vx);
    1809             :     }
    1810           7 :     if (gequal1(lead)) return RgX_copy(x);
    1811           0 :     av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
    1812           0 :     return gerepile(av0,tetpil,RgXQX_red(x,T));
    1813             :   }
    1814       18026 :   av0 = avma; dz = dx-dy;
    1815       18026 :   lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
    1816       18026 :   avma = av0;
    1817       18026 :   z = cgetg(dz+3,t_POL); z[1] = x[1];
    1818       18026 :   x += 2; y += 2; z += 2;
    1819             : 
    1820       18026 :   p1 = gel(x,dx); av = avma;
    1821       18026 :   gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
    1822       84342 :   for (i=dx-1; i>=dy; i--)
    1823             :   {
    1824       66316 :     av=avma; p1=gel(x,i);
    1825       66316 :     for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1826       66316 :     if (lead) p1 = gmul(grem(p1, T), lead);
    1827       66316 :     tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
    1828             :   }
    1829       18026 :   if (!pr) { if (lead) gunclone(lead); return z-2; }
    1830             : 
    1831       18026 :   rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
    1832       28235 :   for (sx=0; ; i--)
    1833             :   {
    1834       28235 :     p1 = gel(x,i);
    1835       28235 :     for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1836       28235 :     tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
    1837       12793 :     if (!i) break;
    1838       10209 :     avma=av;
    1839       10209 :   }
    1840       18026 :   if (pr == ONLY_DIVIDES)
    1841             :   {
    1842        1660 :     if (lead) gunclone(lead);
    1843        1660 :     if (sx) { avma=av0; return NULL; }
    1844        1590 :     avma = (pari_sp)rem; return z-2;
    1845             :   }
    1846       16366 :   lr=i+3; rem -= lr;
    1847       16366 :   rem[0] = evaltyp(t_POL) | evallg(lr);
    1848       16366 :   rem[1] = z[-1];
    1849       16366 :   p1 = gerepile((pari_sp)rem,tetpil,p1);
    1850       16366 :   rem += 2; gel(rem,i) = p1;
    1851       31133 :   for (i--; i>=0; i--)
    1852             :   {
    1853       14767 :     av=avma; p1 = gel(x,i);
    1854       36916 :     for (j=0; j<=i && j<=dz; j++)
    1855       22149 :       p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
    1856       14767 :     tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
    1857             :   }
    1858       16366 :   rem -= 2;
    1859       16366 :   if (lead) gunclone(lead);
    1860       16366 :   if (!sx) (void)normalizepol_lg(rem, lr);
    1861       16366 :   if (pr == ONLY_REM) return gerepileupto(av0,rem);
    1862          70 :   *pr = rem; return z-2;
    1863             : }
    1864             : 
    1865             : /*******************************************************************/
    1866             : /*                                                                 */
    1867             : /*                        PSEUDO-DIVISION                          */
    1868             : /*                                                                 */
    1869             : /*******************************************************************/
    1870             : INLINE GEN
    1871      711772 : rem(GEN c, GEN T)
    1872             : {
    1873      711772 :   if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
    1874      711772 :   return c;
    1875             : }
    1876             : 
    1877             : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
    1878             : int
    1879        3567 : ZXQX_dvd(GEN x, GEN y, GEN T)
    1880             : {
    1881             :   long dx, dy, dz, i, p, T_ismonic;
    1882        3567 :   pari_sp av = avma, av2;
    1883             :   GEN y_lead;
    1884             : 
    1885        3567 :   if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
    1886        3567 :   dy = degpol(y); y_lead = gel(y,dy+2);
    1887        3567 :   if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
    1888             :   /* if monic, no point in using pseudo-division */
    1889        3567 :   if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
    1890        2594 :   T_ismonic = gequal1(leading_coeff(T));
    1891        2594 :   dx = degpol(x);
    1892        2594 :   if (dx < dy) return !signe(x);
    1893        2594 :   (void)new_chunk(2);
    1894        2594 :   x = RgX_recip_shallow(x)+2;
    1895        2594 :   y = RgX_recip_shallow(y)+2;
    1896             :   /* pay attention to sparse divisors */
    1897        5314 :   for (i = 1; i <= dy; i++)
    1898        2720 :     if (!signe(gel(y,i))) gel(y,i) = NULL;
    1899        2594 :   dz = dx-dy; p = dz+1;
    1900        2594 :   av2 = avma;
    1901             :   for (;;)
    1902             :   {
    1903       28685 :     GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
    1904       28685 :     x0 = gneg(x0); p--;
    1905       28685 :     m = gcdii(cx, y0);
    1906       28685 :     if (!equali1(m))
    1907             :     {
    1908       27754 :       x0 = gdiv(x0, m);
    1909       27754 :       y0 = diviiexact(y0, m);
    1910       27754 :       if (equali1(y0)) y0 = NULL;
    1911             :     }
    1912       58196 :     for (i=1; i<=dy; i++)
    1913             :     {
    1914       29511 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    1915       29511 :       if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
    1916       29511 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    1917       29511 :       gel(x,i) = c;
    1918             :     }
    1919      371893 :     for (   ; i<=dx; i++)
    1920             :     {
    1921      343208 :       GEN c = gel(x,i); if (y0) c = gmul(y0, c);
    1922      343208 :       if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
    1923      343208 :       gel(x,i) = c;
    1924             :     }
    1925       31328 :     do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
    1926       28685 :     if (dx < dy) break;
    1927       26091 :     if (gc_needed(av2,1))
    1928             :     {
    1929           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
    1930           0 :       gerepilecoeffs(av2,x,dx+1);
    1931             :     }
    1932       26091 :   }
    1933        2594 :   avma = av; return (dx < 0);
    1934             : }
    1935             : 
    1936             : /* T either NULL or a t_POL. */
    1937             : GEN
    1938       24597 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
    1939             : {
    1940       24597 :   long vx = varn(x), dx, dy, dz, i, lx, p;
    1941       24597 :   pari_sp av = avma, av2;
    1942             :   GEN y_lead;
    1943             : 
    1944       24597 :   if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
    1945       24597 :   dy = degpol(y); y_lead = gel(y,dy+2);
    1946             :   /* if monic, no point in using pseudo-division */
    1947       24597 :   if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
    1948       20845 :   dx = degpol(x);
    1949       20845 :   if (dx < dy) return RgX_copy(x);
    1950       20845 :   (void)new_chunk(2);
    1951       20845 :   x = RgX_recip_shallow(x)+2;
    1952       20845 :   y = RgX_recip_shallow(y)+2;
    1953             :   /* pay attention to sparse divisors */
    1954       64780 :   for (i = 1; i <= dy; i++)
    1955       43935 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    1956       20845 :   dz = dx-dy; p = dz+1;
    1957       20845 :   av2 = avma;
    1958             :   for (;;)
    1959             :   {
    1960       85083 :     gel(x,0) = gneg(gel(x,0)); p--;
    1961      267400 :     for (i=1; i<=dy; i++)
    1962             :     {
    1963      182317 :       GEN c = gmul(y_lead, gel(x,i));
    1964      182317 :       if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
    1965      182317 :       gel(x,i) = rem(c, T);
    1966             :     }
    1967      314097 :     for (   ; i<=dx; i++)
    1968             :     {
    1969      229014 :       GEN c = gmul(y_lead, gel(x,i));
    1970      229014 :       gel(x,i) = rem(c, T);
    1971             :     }
    1972       92013 :     do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
    1973       85083 :     if (dx < dy) break;
    1974       64238 :     if (gc_needed(av2,1))
    1975             :     {
    1976           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
    1977           0 :       gerepilecoeffs(av2,x,dx+1);
    1978             :     }
    1979       64238 :   }
    1980       20845 :   if (dx < 0) return pol_0(vx);
    1981       19102 :   lx = dx+3; x -= 2;
    1982       19102 :   x[0] = evaltyp(t_POL) | evallg(lx);
    1983       19102 :   x[1] = evalsigne(1) | evalvarn(vx);
    1984       19102 :   x = RgX_recip_shallow(x);
    1985       19102 :   if (p)
    1986             :   { /* multiply by y[0]^p   [beware dummy vars from FpX_FpXY_resultant] */
    1987        1127 :     GEN t = y_lead;
    1988        1127 :     if (T && typ(t) == t_POL && varn(t) == varn(T))
    1989           0 :       t = RgXQ_powu(t, p, T);
    1990             :     else
    1991        1127 :       t = gpowgs(t, p);
    1992        3682 :     for (i=2; i<lx; i++)
    1993             :     {
    1994        2555 :       GEN c = gmul(gel(x,i), t);
    1995        2555 :       gel(x,i) = rem(c,T);
    1996             :     }
    1997        1127 :     if (!T) return gerepileupto(av, x);
    1998             :   }
    1999       17975 :   return gerepilecopy(av, x);
    2000             : }
    2001             : 
    2002             : GEN
    2003       24597 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
    2004             : 
    2005             : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
    2006             : GEN
    2007       53845 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
    2008             : {
    2009       53845 :   long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
    2010       53845 :   pari_sp av = avma, av2;
    2011             :   GEN z, r, ypow, y_lead;
    2012             : 
    2013       53845 :   if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
    2014       53845 :   dy = degpol(y); y_lead = gel(y,dy+2);
    2015       53845 :   if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
    2016       27048 :   dx = degpol(x);
    2017       27048 :   if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
    2018       27048 :   if (dx == dy)
    2019             :   {
    2020          28 :     GEN x_lead = gel(x,lg(x)-1);
    2021          28 :     x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
    2022          28 :     y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
    2023          28 :     r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
    2024          28 :     *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
    2025             :   }
    2026       27020 :   (void)new_chunk(2);
    2027       27020 :   x = RgX_recip_shallow(x)+2;
    2028       27020 :   y = RgX_recip_shallow(y)+2;
    2029             :   /* pay attention to sparse divisors */
    2030      118055 :   for (i = 1; i <= dy; i++)
    2031       91035 :     if (isexactzero(gel(y,i))) gel(y,i) = NULL;
    2032       27020 :   dz = dx-dy; p = dz+1;
    2033       27020 :   lz = dz+3;
    2034       27020 :   z = cgetg(lz, t_POL);
    2035       27020 :   z[1] = evalsigne(1) | evalvarn(vx);
    2036       27020 :   for (i = 2; i < lz; i++) gel(z,i) = gen_0;
    2037       27020 :   ypow = new_chunk(dz+1);
    2038       27020 :   gel(ypow,0) = gen_1;
    2039       27020 :   gel(ypow,1) = y_lead;
    2040       34971 :   for (i=2; i<=dz; i++)
    2041             :   {
    2042        7951 :     GEN c = gmul(gel(ypow,i-1), y_lead);
    2043        7951 :     gel(ypow,i) = rem(c,T);
    2044             :   }
    2045       27020 :   av2 = avma;
    2046       27020 :   for (iz=2;;)
    2047             :   {
    2048       56043 :     p--;
    2049       56043 :     gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
    2050      240777 :     for (i=1; i<=dy; i++)
    2051             :     {
    2052      184734 :       GEN c = gmul(y_lead, gel(x,i));
    2053      184734 :       if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
    2054      184734 :       gel(x,i) = rem(c, T);
    2055             :     }
    2056      105201 :     for (   ; i<=dx; i++)
    2057             :     {
    2058       49158 :       GEN c = gmul(y_lead, gel(x,i));
    2059       49158 :       gel(x,i) = rem(c,T);
    2060             :     }
    2061       56043 :     x++; dx--;
    2062       56043 :     while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
    2063       56043 :     if (dx < dy) break;
    2064       29023 :     if (gc_needed(av2,1))
    2065             :     {
    2066           0 :       GEN X = x-2;
    2067           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
    2068           0 :       X[0] = evaltyp(t_POL)|evallg(dx+3); X[1] = z[1]; /* hack */
    2069           0 :       gerepileall(av2,2, &X, &z); x = X+2;
    2070             :     }
    2071       29023 :   }
    2072       27020 :   while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
    2073       27020 :   if (dx < 0)
    2074          98 :     x = pol_0(vx);
    2075             :   else
    2076             :   {
    2077       26922 :     lx = dx+3; x -= 2;
    2078       26922 :     x[0] = evaltyp(t_POL) | evallg(lx);
    2079       26922 :     x[1] = evalsigne(1) | evalvarn(vx);
    2080       26922 :     x = RgX_recip_shallow(x);
    2081             :   }
    2082       27020 :   z = RgX_recip_shallow(z);
    2083       27020 :   r = x;
    2084       27020 :   if (p)
    2085             :   {
    2086        4106 :     GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
    2087        4106 :     if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
    2088             :   }
    2089       27020 :   gerepileall(av, 2, &z, &r);
    2090       27020 :   *ptr = r; return z;
    2091             : }
    2092             : GEN
    2093       53698 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
    2094       53698 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
    2095             : 
    2096             : GEN
    2097           0 : RgXQX_mul(GEN x, GEN y, GEN T)
    2098             : {
    2099           0 :   return RgXQX_red(RgX_mul(x,y), T);
    2100             : }
    2101             : GEN
    2102    68928210 : RgX_Rg_mul(GEN y, GEN x) {
    2103             :   long i, ly;
    2104    68928210 :   GEN z = cgetg_copy(y, &ly); z[1] = y[1];
    2105    68928210 :   if (ly == 2) return z;
    2106    68870600 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2107    68870593 :   return normalizepol_lg(z,ly);
    2108             : }
    2109             : GEN
    2110       13881 : RgX_muls(GEN y, long x) {
    2111             :   long i, ly;
    2112       13881 :   GEN z = cgetg_copy(y, &ly); z[1] = y[1];
    2113       13881 :   if (ly == 2) return z;
    2114       13846 :   for (i = 2; i < ly; i++) gel(z,i) = gmulsg(x,gel(y,i));
    2115       13846 :   return normalizepol_lg(z,ly);
    2116             : }
    2117             : GEN
    2118          28 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T)
    2119             : {
    2120          28 :   return RgXQX_red(RgX_Rg_mul(x,y), T);
    2121             : }
    2122             : GEN
    2123          56 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T)
    2124             : {
    2125          56 :   return RgXQV_red(RgV_Rg_mul(v,x), T);
    2126             : }
    2127             : 
    2128             : GEN
    2129           0 : RgXQX_sqr(GEN x, GEN T)
    2130             : {
    2131           0 :   return RgXQX_red(RgX_sqr(x), T);
    2132             : }
    2133             : 
    2134             : static GEN
    2135       65093 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
    2136             : static GEN
    2137           0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
    2138             : static GEN
    2139      248259 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
    2140             : static GEN
    2141      101911 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
    2142             : static GEN
    2143      110593 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
    2144             : static GEN
    2145      105448 : _one(void *data) { return pol_1(varn((GEN)data)); }
    2146             : static GEN
    2147         105 : _zero(void *data) { return pol_0(varn((GEN)data)); }
    2148             : static GEN
    2149       70609 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
    2150             : 
    2151             : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
    2152             :               _mul, _sqr, _one, _zero };
    2153             : 
    2154             : GEN
    2155           0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
    2156             : {
    2157           0 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
    2158             : }
    2159             : 
    2160             : GEN
    2161       44926 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
    2162             : {
    2163       44926 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2164       44926 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
    2165             : }
    2166             : 
    2167             : /* mod X^n */
    2168             : struct modXn {
    2169             :   long v; /* varn(X) */
    2170             :   long n;
    2171             : } ;
    2172             : static GEN
    2173        1806 : _sqrXn(void *data, GEN x) {
    2174        1806 :   struct modXn *S = (struct modXn*)data;
    2175        1806 :   return RgXn_sqr(x, S->n);
    2176             : }
    2177             : static GEN
    2178        1190 : _mulXn(void *data, GEN x, GEN y) {
    2179        1190 :   struct modXn *S = (struct modXn*)data;
    2180        1190 :   return RgXn_mul(x,y, S->n);
    2181             : }
    2182             : static GEN
    2183        1407 : _oneXn(void *data) {
    2184        1407 :   struct modXn *S = (struct modXn*)data;
    2185        1407 :   return pol_1(S->v);
    2186             : }
    2187             : static GEN
    2188           0 : _zeroXn(void *data) {
    2189           0 :   struct modXn *S = (struct modXn*)data;
    2190           0 :   return pol_0(S->v);
    2191             : }
    2192             : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
    2193             :                                           _oneXn, _zeroXn };
    2194             : 
    2195             : GEN
    2196         336 : RgXn_powers(GEN x, long m, long n)
    2197             : {
    2198         336 :   long d = degpol(x);
    2199         336 :   int use_sqr = (d<<1) >= n;
    2200             :   struct modXn S;
    2201         336 :   S.v = varn(x); S.n = n;
    2202         336 :   return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
    2203             : }
    2204             : 
    2205             : GEN
    2206        1526 : RgXn_powu_i(GEN x, ulong m, long n)
    2207             : {
    2208             :   struct modXn S;
    2209        1526 :   S.v = varn(x); S.n = n;
    2210        1526 :   return gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
    2211             : }
    2212             : GEN
    2213           0 : RgXn_powu(GEN x, ulong m, long n)
    2214             : {
    2215             :   struct modXn S;
    2216           0 :   S.v = varn(x); S.n = n;
    2217           0 :   return gen_powu(x, m, (void*)&S,_sqrXn,_mulXn);
    2218             : }
    2219             : 
    2220             : GEN
    2221         672 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
    2222             : {
    2223             :   struct modXn S;
    2224         672 :   S.v = varn(gel(x,2)); S.n = n;
    2225         672 :   return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
    2226             : }
    2227             : 
    2228             : GEN
    2229           0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
    2230             : {
    2231           0 :   int use_sqr = 2*degpol(x) >= n;
    2232             :   struct modXn S;
    2233           0 :   S.v = varn(x); S.n = n;
    2234           0 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2235             : }
    2236             : 
    2237             : /* Q(x) mod t^n, x in R[t], n >= 1 */
    2238             : GEN
    2239        1799 : RgXn_eval(GEN Q, GEN x, long n)
    2240             : {
    2241        1799 :   long d = degpol(x);
    2242             :   int use_sqr;
    2243             :   struct modXn S;
    2244        1799 :   if (d == 1 && isrationalzero(gel(x,2)))
    2245             :   {
    2246        1792 :     GEN y = RgX_unscale(Q, gel(x,3));
    2247        1792 :     setvarn(y, varn(x)); return y;
    2248             :   }
    2249           7 :   S.v = varn(x);
    2250           7 :   S.n = n;
    2251           7 :   use_sqr = (d<<1) >= n;
    2252           7 :   return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
    2253             : }
    2254             : 
    2255             : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
    2256             : static GEN
    2257       33774 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
    2258             : {
    2259       33774 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2260       33774 :   return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
    2261             : }
    2262             : 
    2263             : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
    2264             : static GEN
    2265           0 : RgXn_sqrhigh(GEN f, long n2, long n)
    2266             : {
    2267           0 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
    2268           0 :   return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
    2269             : }
    2270             : 
    2271             : GEN
    2272      112595 : RgXn_inv(GEN f, long e)
    2273             : {
    2274      112595 :   pari_sp av = avma, av2;
    2275             :   ulong mask;
    2276             :   GEN W, a;
    2277      112595 :   long v = varn(f), n = 1;
    2278             : 
    2279      112595 :   if (!signe(f)) pari_err_INV("RgXn_inv",f);
    2280      112595 :   a = ginv(gel(f,2));
    2281      112595 :   if (e == 1) return scalarpol(a, v);
    2282      112595 :   else if (e == 2)
    2283             :   {
    2284             :     GEN b;
    2285       98994 :     if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
    2286       80185 :     b = gneg(b);
    2287       80185 :     if (!gequal1(a)) b = gmul(b, gsqr(a));
    2288       80185 :     W = deg1pol_shallow(b, a, v);
    2289       80185 :     return gerepilecopy(av, W);
    2290             :   }
    2291       13601 :   W = scalarpol_shallow(ginv(gel(f,2)),v);
    2292       13601 :   mask = quadratic_prec_mask(e);
    2293       13601 :   av2 = avma;
    2294       60976 :   for (;mask>1;)
    2295             :   {
    2296             :     GEN u, fr;
    2297       33774 :     long n2 = n;
    2298       33774 :     n<<=1; if (mask & 1) n--;
    2299       33774 :     mask >>= 1;
    2300       33774 :     fr = RgXn_red_shallow(f, n);
    2301       33774 :     u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
    2302       33774 :     W = RgX_sub(W, RgX_shift_shallow(u, n2));
    2303       33774 :     if (gc_needed(av2,2))
    2304             :     {
    2305           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
    2306           0 :       W = gerepileupto(av2, W);
    2307             :     }
    2308             :   }
    2309       13601 :   return gerepileupto(av, W);
    2310             : }
    2311             : 
    2312             : GEN
    2313       12586 : RgXn_exp(GEN h, long e)
    2314             : {
    2315       12586 :   pari_sp av = avma, av2;
    2316       12586 :   long v = varn(h), n=1;
    2317       12586 :   GEN f = pol_1(v), g = pol_1(v);
    2318       12586 :   ulong mask = quadratic_prec_mask(e);
    2319       12586 :   av2 = avma;
    2320       12586 :   if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
    2321           0 :     pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
    2322       38822 :   for (;mask>1;)
    2323             :   {
    2324             :     GEN q, w;
    2325       13650 :     long n2 = n;
    2326       13650 :     n<<=1; if (mask & 1) n--;
    2327       13650 :     mask >>= 1;
    2328       13650 :     g = RgX_sub(RgX_muls(g,2),RgXn_mul(f,RgXn_sqr(g,n2),n2));
    2329       13650 :     q = RgX_deriv(RgXn_red_shallow(h,n2));
    2330       13650 :     w = RgX_add(q, RgXn_mul(g, RgX_sub(RgX_deriv(f), RgXn_mul(f,q,n-1)),n-1));
    2331       13650 :     f = RgX_add(f, RgXn_mul(f, RgX_sub(RgXn_red_shallow(h, n), RgX_integ(w)), n));
    2332       13650 :     if (gc_needed(av2,2))
    2333             :     {
    2334           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_exp, e = %ld", n);
    2335           0 :       gerepileall(av2, 2, &f, &g);
    2336             :     }
    2337             :   }
    2338       12586 :   return gerepileupto(av, f);
    2339             : }
    2340             : 
    2341             : GEN
    2342          84 : RgXn_reverse(GEN f, long e)
    2343             : {
    2344          84 :   pari_sp av = avma, av2;
    2345             :   ulong mask;
    2346             :   GEN fi, a, df, W, an;
    2347          84 :   long v = varn(f), n=1;
    2348          84 :   if (degpol(f)<1 || !gequal0(gel(f,2)))
    2349           0 :     pari_err_INV("serreverse",f);
    2350          84 :   fi = ginv(gel(f,3));
    2351          84 :   a = deg1pol_shallow(fi,gen_0,v);
    2352          84 :   if (e <= 2) return gerepilecopy(av, a);
    2353          84 :   W = scalarpol(fi,v);
    2354          84 :   df = RgX_deriv(f);
    2355          84 :   mask = quadratic_prec_mask(e);
    2356          84 :   av2 = avma;
    2357         504 :   for (;mask>1;)
    2358             :   {
    2359             :     GEN u, fa, fr;
    2360         336 :     long n2 = n, rt;
    2361         336 :     n<<=1; if (mask & 1) n--;
    2362         336 :     mask >>= 1;
    2363         336 :     fr = RgXn_red_shallow(f, n);
    2364         336 :     rt = brent_kung_optpow(degpol(fr), 4, 3);
    2365         336 :     an = RgXn_powers(a, rt, n);
    2366         336 :     if (n>1)
    2367             :     {
    2368         336 :       long n4 = (n2+1)>>1;
    2369         336 :       GEN dfr = RgXn_red_shallow(df, n2);
    2370         336 :       dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
    2371         336 :       u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
    2372         336 :       W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
    2373             :     }
    2374         336 :     fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
    2375         336 :     fa = RgX_shift(fa, -n2);
    2376         336 :     a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
    2377         336 :     if (gc_needed(av2,2))
    2378             :     {
    2379           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
    2380           0 :       gerepileall(av2, 2, &a, &W);
    2381             :     }
    2382             :   }
    2383          84 :   return gerepileupto(av, a);
    2384             : }
    2385             : 
    2386             : GEN
    2387           0 : RgXn_sqrt(GEN h, long e)
    2388             : {
    2389           0 :   pari_sp av = avma, av2;
    2390           0 :   long v = varn(h), n = 1;
    2391           0 :   GEN f = scalarpol(gen_1, v), df = f;
    2392           0 :   ulong mask = quadratic_prec_mask(e);
    2393           0 :   if (degpol(h)<0 || !gequal1(gel(h,2)))
    2394           0 :     pari_err_SQRTN("RgXn_sqrt",h);
    2395           0 :   av2 = avma;
    2396             :   while(1)
    2397             :   {
    2398           0 :     long n2 = n, m;
    2399             :     GEN g;
    2400           0 :     n<<=1; if (mask & 1) n--;
    2401           0 :     mask >>= 1;
    2402           0 :     m = n-n2;
    2403           0 :     g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
    2404           0 :     f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
    2405           0 :     if (mask==1) return gerepileupto(av, f);
    2406           0 :     g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
    2407           0 :     df = RgX_sub(df, RgX_shift_shallow(g, n2));
    2408           0 :     if (gc_needed(av2,2))
    2409             :     {
    2410           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
    2411           0 :       gerepileall(av2, 2, &f, &df);
    2412             :     }
    2413           0 :   }
    2414             : }
    2415             : 
    2416             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2417             : GEN
    2418      205612 : RgXQ_powu(GEN x, ulong n, GEN T)
    2419             : {
    2420             :   pari_sp av;
    2421             :   GEN y;
    2422             : 
    2423      205612 :   if (!n) return pol_1(varn(x));
    2424      204058 :   if (n == 1) return RgX_copy(x);
    2425      139961 :   av = avma;
    2426      139961 :   y = gen_powu(x, n, (void*)T, &_sqr, &_mul);
    2427      139959 :   return gerepileupto(av, y);
    2428             : }
    2429             : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
    2430             : GEN
    2431       17611 : RgXQ_pow(GEN x, GEN n, GEN T)
    2432             : {
    2433             :   pari_sp av;
    2434       17611 :   long s = signe(n);
    2435             :   GEN y;
    2436             : 
    2437       17611 :   if (!s) return pol_1(varn(x));
    2438       17611 :   if (is_pm1(n) == 1)
    2439           0 :     return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
    2440       17611 :   av = avma;
    2441       17611 :   if (s < 0) x = RgXQ_inv(x, T);
    2442       17611 :   y = gen_pow(x, n, (void*)T, &_sqr, &_mul);
    2443       17611 :   return gerepileupto(av, y);
    2444             : }
    2445             : 
    2446             : /* generates the list of powers of x of degree 0,1,2,...,l*/
    2447             : GEN
    2448        2219 : RgXQ_powers(GEN x, long l, GEN T)
    2449             : {
    2450        2219 :   int use_sqr = 2*degpol(x) >= degpol(T);
    2451        2219 :   return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
    2452             : }
    2453             : 
    2454             : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
    2455             : GEN
    2456        1834 : QXQ_powers(GEN a, long n, GEN T)
    2457             : {
    2458        1834 :   GEN den, v = RgXQ_powers(Q_remove_denom(a, &den), n, T);
    2459             :   /* den*a integral; v[i+1] = (den*a)^i in K */
    2460        1834 :   if (den)
    2461             :   { /* restore denominators */
    2462        1246 :     GEN d = den;
    2463             :     long i;
    2464        1246 :     gel(v,2) = a;
    2465        3829 :     for (i=3; i<=n+1; i++) {
    2466        2583 :       d = mulii(d,den);
    2467        2583 :       gel(v,i) = RgX_Rg_div(gel(v,i), d);
    2468             :     }
    2469             :   }
    2470        1834 :   return v;
    2471             : }
    2472             : 
    2473             : static GEN
    2474        1176 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
    2475             : {
    2476        1176 :   long l, i, m = 0;
    2477             :   GEN dz, z;
    2478        1176 :   GEN V = cgetg_copy(v, &l);
    2479        3808 :   for (i = imin; i < l; i++)
    2480             :   {
    2481        2632 :     GEN c = gel(v, i);
    2482        2632 :     if (typ(c) == t_POL) m = maxss(m, degpol(c));
    2483             :   }
    2484        1176 :   z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
    2485        1176 :   for (i = 1; i < imin; i++) V[i] = v[i];
    2486        3808 :   for (i = imin; i < l; i++)
    2487             :   {
    2488        2632 :     GEN c = gel(v,i);
    2489        2632 :     if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
    2490        2632 :     gel(V,i) = c;
    2491             :   }
    2492        1176 :   return V;
    2493             : }
    2494             : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
    2495             : GEN
    2496        1113 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
    2497        1113 : { return do_QXQ_eval(v, 1, a, T); }
    2498             : GEN
    2499          63 : QXX_QXQ_eval(GEN v, GEN a, GEN T)
    2500          63 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
    2501             : 
    2502             : GEN
    2503         287 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
    2504             : {
    2505         287 :   return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
    2506             : }
    2507             : 
    2508             : GEN
    2509          56 : RgXQ_minpoly_naive(GEN y, GEN P)
    2510             : {
    2511          56 :   pari_sp ltop=avma;
    2512          56 :   long n=lgpol(P);
    2513          56 :   GEN M=ker(RgXQ_matrix_pow(y,n,n,P));
    2514          56 :   M=content(RgM_to_RgXV(M,varn(P)));
    2515          56 :   return gerepileupto(ltop,M);
    2516             : }
    2517             : 
    2518             : GEN
    2519       34668 : RgXQ_norm(GEN x, GEN T)
    2520             : {
    2521             :   pari_sp av;
    2522       34668 :   long dx = degpol(x);
    2523             :   GEN L, y;
    2524             : 
    2525       34668 :   av = avma; y = resultant(T, x);
    2526       34668 :   L = leading_coeff(T);
    2527       34668 :   if (gequal1(L) || !signe(x)) return y;
    2528           0 :   return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
    2529             : }
    2530             : 
    2531             : GEN
    2532      103886 : RgX_blocks(GEN P, long n, long m)
    2533             : {
    2534      103886 :   GEN z = cgetg(m+1,t_VEC);
    2535      103886 :   long i,j, k=2, l = lg(P);
    2536      504830 :   for(i=1; i<=m; i++)
    2537             :   {
    2538      400944 :     GEN zi = cgetg(n+2,t_POL);
    2539      400944 :     zi[1] = P[1];
    2540      400944 :     gel(z,i) = zi;
    2541     2372681 :     for(j=2; j<n+2; j++)
    2542     1971737 :       gel(zi, j) = k==l ? gen_0 : gel(P,k++);
    2543      400944 :     zi = RgX_renormalize_lg(zi, n+2);
    2544             :   }
    2545      103886 :   return z;
    2546             : }
    2547             : 
    2548             : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
    2549             : void
    2550       30211 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
    2551             : {
    2552       30211 :   long n = degpol(p), v = varn(p), n0, n1, i;
    2553             :   GEN p0, p1;
    2554             : 
    2555       60424 :   if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
    2556             : 
    2557       30212 :   n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
    2558       30212 :   p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
    2559       30212 :   p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
    2560      928712 :   for (i=0; i<n1; i++)
    2561             :   {
    2562      898500 :     p0[2+i] = p[2+(i<<1)];
    2563      898500 :     p1[2+i] = p[3+(i<<1)];
    2564             :   }
    2565       30212 :   if (n1 != n0)
    2566       21715 :     p0[2+i] = p[2+(i<<1)];
    2567       30212 :   *pe = normalizepol(p0);
    2568       30211 :   *po = normalizepol(p1);
    2569             : }
    2570             : 
    2571             : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
    2572             : GEN
    2573       40614 : RgX_splitting(GEN p, long k)
    2574             : {
    2575       40614 :   long n = degpol(p), v = varn(p), m, i, j, l;
    2576             :   GEN r;
    2577             : 
    2578       40614 :   m = n/k;
    2579       40614 :   r = cgetg(k+1,t_VEC);
    2580      223930 :   for(i=1; i<=k; i++)
    2581             :   {
    2582      183316 :     gel(r,i) = cgetg(m+3, t_POL);
    2583      183316 :     mael(r,i,1) = evalvarn(v)|evalsigne(1);
    2584             :   }
    2585      552426 :   for (j=1, i=0, l=2; i<=n; i++)
    2586             :   {
    2587      511812 :     gmael(r,j,l) = gel(p,2+i);
    2588      511812 :     if (j==k) { j=1; l++; } else j++;
    2589             :   }
    2590      223930 :   for(i=1; i<=k; i++)
    2591      183316 :     gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
    2592       40614 :   return r;
    2593             : }
    2594             : 
    2595             : /*******************************************************************/
    2596             : /*                                                                 */
    2597             : /*                        Kronecker form                           */
    2598             : /*                                                                 */
    2599             : /*******************************************************************/
    2600             : 
    2601             : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
    2602             :  * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
    2603             :  * normalized coefficients */
    2604             : GEN
    2605         189 : Kronecker_to_mod(GEN z, GEN T)
    2606             : {
    2607         189 :   long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
    2608         189 :   GEN x, t = cgetg(N,t_POL);
    2609         189 :   t[1] = T[1];
    2610         189 :   lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
    2611         189 :   x[1] = z[1];
    2612         189 :   T = RgX_copy(T);
    2613        4389 :   for (i=2; i<lx+2; i++, z+= N-2)
    2614             :   {
    2615        4200 :     for (j=2; j<N; j++) gel(t,j) = gel(z,j);
    2616        4200 :     gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2617             :   }
    2618         189 :   N = (l-2) % (N-2) + 2;
    2619         189 :   for (j=2; j<N; j++) t[j] = z[j];
    2620         189 :   gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
    2621         189 :   return normalizepol_lg(x, i+1);
    2622             : }

Generated by: LCOV version 1.11