"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/RgX.c" (14 Jan 2021, 77499 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "RgX.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

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