"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/gen3.c" (14 Jan 2021, 113362 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 "gen3.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 /********************************************************************/
   16 /**                                                                **/
   17 /**                      GENERIC OPERATIONS                        **/
   18 /**                         (third part)                           **/
   19 /**                                                                **/
   20 /********************************************************************/
   21 #include "pari.h"
   22 #include "paripriv.h"
   23 
   24 /********************************************************************/
   25 /**                                                                **/
   26 /**                 PRINCIPAL VARIABLE NUMBER                      **/
   27 /**                                                                **/
   28 /********************************************************************/
   29 static void
   30 recvar(hashtable *h, GEN x)
   31 {
   32   long i = 1, lx = lg(x);
   33   void *v;
   34   switch(typ(x))
   35   {
   36     case t_POL: case t_SER:
   37       v = (void*)varn(x);
   38       if (!hash_search(h, v)) hash_insert(h, v, NULL);
   39       i = 2; break;
   40     case t_POLMOD: case t_RFRAC:
   41     case t_VEC: case t_COL: case t_MAT: break;
   42     case t_LIST:
   43       x = list_data(x);
   44       lx = x? lg(x): 1; break;
   45     default:
   46       return;
   47   }
   48   for (; i < lx; i++) recvar(h, gel(x,i));
   49 }
   50 
   51 GEN
   52 variables_vecsmall(GEN x)
   53 {
   54   hashtable *h = hash_create_ulong(100, 1);
   55   recvar(h, x);
   56   return vars_sort_inplace(hash_keys(h));
   57 }
   58 
   59 GEN
   60 variables_vec(GEN x)
   61 { return x? vars_to_RgXV(variables_vecsmall(x)): gpolvar(NULL); }
   62 
   63 long
   64 gvar(GEN x)
   65 {
   66   long i, v, w, lx;
   67   switch(typ(x))
   68   {
   69     case t_POL: case t_SER: return varn(x);
   70     case t_POLMOD: return varn(gel(x,1));
   71     case t_RFRAC:  return varn(gel(x,2));
   72     case t_VEC: case t_COL: case t_MAT:
   73       lx = lg(x); break;
   74     case t_LIST:
   75       x = list_data(x);
   76       lx = x? lg(x): 1; break;
   77     default:
   78       return NO_VARIABLE;
   79   }
   80   v = NO_VARIABLE;
   81   for (i=1; i < lx; i++) { w = gvar(gel(x,i)); if (varncmp(w,v) < 0) v = w; }
   82   return v;
   83 }
   84 /* T main polynomial in R[X], A auxiliary in R[X] (possibly degree 0).
   85  * Guess and return the main variable of R */
   86 static long
   87 var2_aux(GEN T, GEN A)
   88 {
   89   long a = gvar2(T);
   90   long b = (typ(A) == t_POL && varn(A) == varn(T))? gvar2(A): gvar(A);
   91   if (varncmp(a, b) > 0) a = b;
   92   return a;
   93 }
   94 static long
   95 var2_rfrac(GEN x)  { return var2_aux(gel(x,2), gel(x,1)); }
   96 static long
   97 var2_polmod(GEN x) { return var2_aux(gel(x,1), gel(x,2)); }
   98 
   99 /* main variable of x, with the convention that the "natural" main
  100  * variable of a POLMOD is mute, so we want the next one. */
  101 static long
  102 gvar9(GEN x)
  103 { return (typ(x) == t_POLMOD)? var2_polmod(x): gvar(x); }
  104 
  105 /* main variable of the ring over wich x is defined */
  106 long
  107 gvar2(GEN x)
  108 {
  109   long i, v, w;
  110   switch(typ(x))
  111   {
  112     case t_POLMOD:
  113       return var2_polmod(x);
  114     case t_POL: case t_SER:
  115       v = NO_VARIABLE;
  116       for (i=2; i < lg(x); i++) {
  117         w = gvar9(gel(x,i));
  118         if (varncmp(w,v) < 0) v=w;
  119       }
  120       return v;
  121     case t_RFRAC:
  122       return var2_rfrac(x);
  123     case t_VEC: case t_COL: case t_MAT:
  124       v = NO_VARIABLE;
  125       for (i=1; i < lg(x); i++) {
  126         w = gvar2(gel(x,i));
  127         if (varncmp(w,v)<0) v=w;
  128       }
  129       return v;
  130   }
  131   return NO_VARIABLE;
  132 }
  133 
  134 /*******************************************************************/
  135 /*                                                                 */
  136 /*                    PRECISION OF SCALAR OBJECTS                  */
  137 /*                                                                 */
  138 /*******************************************************************/
  139 static long
  140 prec0(long e) { return (e < 0)? nbits2prec(-e): 2; }
  141 static long
  142 precREAL(GEN x) { return signe(x) ? realprec(x): prec0(expo(x)); }
  143 /* x t_REAL, y an exact noncomplex type. Return precision(|x| + |y|) */
  144 static long
  145 precrealexact(GEN x, GEN y)
  146 {
  147   long lx, ey = gexpo(y), ex, e;
  148   if (ey == -(long)HIGHEXPOBIT) return precREAL(x);
  149   ex = expo(x);
  150   e = ey - ex;
  151   if (!signe(x)) return prec0((e >= 0)? -e: ex);
  152   lx = realprec(x);
  153   return (e > 0)? lx + nbits2extraprec(e): lx;
  154 }
  155 static long
  156 precCOMPLEX(GEN z)
  157 { /* ~ precision(|x| + |y|) */
  158   GEN x = gel(z,1), y = gel(z,2);
  159   long e, ex, ey, lz, lx, ly;
  160   if (typ(x) != t_REAL) {
  161     if (typ(y) != t_REAL) return 0;
  162     return precrealexact(y, x);
  163   }
  164   if (typ(y) != t_REAL) return precrealexact(x, y);
  165   /* x, y are t_REALs, cf addrr_sign */
  166   ex = expo(x);
  167   ey = expo(y);
  168   e = ey - ex;
  169   if (!signe(x)) {
  170     if (!signe(y)) return prec0( minss(ex,ey) );
  171     if (e <= 0) return prec0(ex);
  172     lz = nbits2prec(e);
  173     ly = realprec(y); if (lz > ly) lz = ly;
  174     return lz;
  175   }
  176   if (!signe(y)) {
  177     if (e >= 0) return prec0(ey);
  178     lz = nbits2prec(-e);
  179     lx = realprec(x); if (lz > lx) lz = lx;
  180     return lz;
  181   }
  182   if (e < 0) { swap(x, y); e = -e; }
  183   lx = realprec(x);
  184   ly = realprec(y);
  185   if (e) {
  186     long d = nbits2extraprec(e), l = ly-d;
  187     return (l > lx)? lx + d: ly;
  188   }
  189   return minss(lx, ly);
  190 }
  191 long
  192 precision(GEN z)
  193 {
  194   switch(typ(z))
  195   {
  196     case t_REAL: return precREAL(z);
  197     case t_COMPLEX: return precCOMPLEX(z);
  198   }
  199   return 0;
  200 }
  201 
  202 long
  203 gprecision(GEN x)
  204 {
  205   long i, k, l;
  206 
  207   switch(typ(x))
  208   {
  209     case t_REAL: return precREAL(x);
  210     case t_COMPLEX: return precCOMPLEX(x);
  211     case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT:
  212     case t_PADIC: case t_QUAD: case t_POLMOD:
  213       return 0;
  214 
  215     case t_POL: case t_SER:
  216       k = LONG_MAX;
  217       for (i=lg(x)-1; i>1; i--)
  218       {
  219         l = gprecision(gel(x,i));
  220         if (l && l<k) k = l;
  221       }
  222       return (k==LONG_MAX)? 0: k;
  223     case t_VEC: case t_COL: case t_MAT:
  224       k = LONG_MAX;
  225       for (i=lg(x)-1; i>0; i--)
  226       {
  227         l = gprecision(gel(x,i));
  228         if (l && l<k) k = l;
  229       }
  230       return (k==LONG_MAX)? 0: k;
  231 
  232     case t_RFRAC:
  233     {
  234       k=gprecision(gel(x,1));
  235       l=gprecision(gel(x,2)); if (l && (!k || l<k)) k=l;
  236       return k;
  237     }
  238     case t_QFR:
  239       return gprecision(gel(x,4));
  240   }
  241   return 0;
  242 }
  243 
  244 static long
  245 vec_padicprec_relative(GEN x, long imin)
  246 {
  247   long s, t, i;
  248   for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
  249   {
  250     t = padicprec_relative(gel(x,i)); if (t<s) s = t;
  251   }
  252   return s;
  253 }
  254 /* RELATIVE padic precision. Only accept decent types: don't try to make sense
  255  * of everything like padicprec */
  256 long
  257 padicprec_relative(GEN x)
  258 {
  259   switch(typ(x))
  260   {
  261     case t_INT: case t_FRAC:
  262       return LONG_MAX;
  263     case t_PADIC:
  264       return signe(gel(x,4))? precp(x): 0;
  265     case t_POLMOD: case t_VEC: case t_COL: case t_MAT:
  266       return vec_padicprec_relative(x, 1);
  267     case t_POL: case t_SER:
  268       return vec_padicprec_relative(x, 2);
  269   }
  270   pari_err_TYPE("padicprec_relative",x);
  271   return 0; /* LCOV_EXCL_LINE */
  272 }
  273 
  274 static long
  275 vec_padicprec(GEN x, GEN p, long imin)
  276 {
  277   long s, t, i;
  278   for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
  279   {
  280     t = padicprec(gel(x,i),p); if (t<s) s = t;
  281   }
  282   return s;
  283 }
  284 static long
  285 vec_serprec(GEN x, long v, long imin)
  286 {
  287   long s, t, i;
  288   for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
  289   {
  290     t = serprec(gel(x,i),v); if (t<s) s = t;
  291   }
  292   return s;
  293 }
  294 
  295 /* ABSOLUTE padic precision */
  296 long
  297 padicprec(GEN x, GEN p)
  298 {
  299   if (typ(p) != t_INT) pari_err_TYPE("padicprec",p);
  300   switch(typ(x))
  301   {
  302     case t_INT: case t_FRAC:
  303       return LONG_MAX;
  304 
  305     case t_INTMOD:
  306       return Z_pval(gel(x,1),p);
  307 
  308     case t_PADIC:
  309       if (!equalii(gel(x,2),p)) pari_err_MODULUS("padicprec", gel(x,2), p);
  310       return precp(x)+valp(x);
  311 
  312     case t_POL: case t_SER:
  313       return vec_padicprec(x, p, 2);
  314     case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_RFRAC:
  315     case t_VEC: case t_COL: case t_MAT:
  316       return vec_padicprec(x, p, 1);
  317   }
  318   pari_err_TYPE("padicprec",x);
  319   return 0; /* LCOV_EXCL_LINE */
  320 }
  321 GEN
  322 gppadicprec(GEN x, GEN p)
  323 {
  324   long v = padicprec(x,p);
  325   return v == LONG_MAX? mkoo(): stoi(v);
  326 }
  327 
  328 /* ABSOLUTE padic precision */
  329 long
  330 serprec(GEN x, long v)
  331 {
  332   long w;
  333   switch(typ(x))
  334   {
  335     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
  336     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_QFR:
  337       return LONG_MAX;
  338 
  339     case t_POL:
  340       w = varn(x);
  341       if (varncmp(v,w) <= 0) return LONG_MAX;
  342       return vec_serprec(x, v, 2);
  343     case t_SER:
  344       w = varn(x);
  345       if (w == v) return lg(x)-2+valp(x);
  346       if (varncmp(v,w) < 0) return LONG_MAX;
  347       return vec_serprec(x, v, 2);
  348     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
  349       return vec_serprec(x, v, 1);
  350   }
  351   pari_err_TYPE("serprec",x);
  352   return 0; /* LCOV_EXCL_LINE */
  353 }
  354 GEN
  355 gpserprec(GEN x, long v)
  356 {
  357   long p = serprec(x,v);
  358   return p == LONG_MAX? mkoo(): stoi(p);
  359 }
  360 
  361 /* Degree of x (scalar, t_POL, t_RFRAC) wrt variable v if v >= 0,
  362  * wrt to main variable if v < 0. */
  363 long
  364 poldegree(GEN x, long v)
  365 {
  366   const long DEGREE0 = -LONG_MAX;
  367   long tx = typ(x), lx,w,i,d;
  368 
  369   if (is_scalar_t(tx)) return gequal0(x)? DEGREE0: 0;
  370   switch(tx)
  371   {
  372     case t_POL:
  373       if (!signe(x)) return DEGREE0;
  374       w = varn(x);
  375       if (v < 0 || v == w) return degpol(x);
  376       if (varncmp(v, w) < 0) return 0;
  377       lx = lg(x); d = DEGREE0;
  378       for (i=2; i<lx; i++)
  379       {
  380         long e = poldegree(gel(x,i), v);
  381         if (e > d) d = e;
  382       }
  383       return d;
  384 
  385     case t_RFRAC:
  386     {
  387       GEN a = gel(x,1), b = gel(x,2);
  388       if (gequal0(a)) return DEGREE0;
  389       if (v < 0)
  390       {
  391         v = varn(b); d = -degpol(b);
  392         if (typ(a) == t_POL && varn(a) == v) d += degpol(a);
  393         return d;
  394       }
  395       return poldegree(a,v) - poldegree(b,v);
  396     }
  397   }
  398   pari_err_TYPE("degree",x);
  399   return 0; /* LCOV_EXCL_LINE  */
  400 }
  401 GEN
  402 gppoldegree(GEN x, long v)
  403 {
  404   long d = poldegree(x,v);
  405   return d == -LONG_MAX? mkmoo(): stoi(d);
  406 }
  407 
  408 /* assume v >= 0 and x is a POLYNOMIAL in v, return deg_v(x) */
  409 long
  410 RgX_degree(GEN x, long v)
  411 {
  412   long tx = typ(x), lx, w, i, d;
  413 
  414   if (is_scalar_t(tx)) return gequal0(x)? -1: 0;
  415   switch(tx)
  416   {
  417     case t_POL:
  418       if (!signe(x)) return -1;
  419       w = varn(x);
  420       if (v == w) return degpol(x);
  421       if (varncmp(v, w) < 0) return 0;
  422       lx = lg(x); d = -1;
  423       for (i=2; i<lx; i++)
  424       {
  425         long e = RgX_degree(gel(x,i), v);
  426         if (e > d) d = e;
  427       }
  428       return d;
  429 
  430     case t_RFRAC:
  431       w = varn(gel(x,2));
  432       if (varncmp(v, w) < 0) return 0;
  433       if (RgX_degree(gel(x,2),v)) pari_err_TYPE("RgX_degree", x);
  434       return RgX_degree(gel(x,1),v);
  435   }
  436   pari_err_TYPE("RgX_degree",x);
  437   return 0; /* LCOV_EXCL_LINE  */
  438 }
  439 
  440 long
  441 degree(GEN x)
  442 {
  443   return poldegree(x,-1);
  444 }
  445 
  446 /* If v<0, leading coeff with respect to the main variable, otherwise wrt v. */
  447 GEN
  448 pollead(GEN x, long v)
  449 {
  450   long tx = typ(x), w;
  451   pari_sp av;
  452 
  453   if (is_scalar_t(tx)) return gcopy(x);
  454   w = varn(x);
  455   switch(tx)
  456   {
  457     case t_POL:
  458       if (v < 0 || v == w)
  459       {
  460         long l = lg(x);
  461         return (l==2)? gen_0: gcopy(gel(x,l-1));
  462       }
  463       break;
  464 
  465     case t_SER:
  466       if (v < 0 || v == w) return signe(x)? gcopy(gel(x,2)): gen_0;
  467       if (varncmp(v, w) > 0) x = polcoef_i(x, valp(x), v);
  468       break;
  469 
  470     default:
  471       pari_err_TYPE("pollead",x);
  472       return NULL; /* LCOV_EXCL_LINE */
  473   }
  474   if (varncmp(v, w) < 0) return gcopy(x);
  475   av = avma; w = fetch_var_higher();
  476   x = gsubst(x, v, pol_x(w));
  477   x = pollead(x, w);
  478   delete_var(); return gerepileupto(av, x);
  479 }
  480 
  481 /* returns 1 if there's a real component in the structure, 0 otherwise */
  482 int
  483 isinexactreal(GEN x)
  484 {
  485   long i;
  486   switch(typ(x))
  487   {
  488     case t_REAL: return 1;
  489     case t_COMPLEX: return (typ(gel(x,1))==t_REAL || typ(gel(x,2))==t_REAL);
  490 
  491     case t_INT: case t_INTMOD: case t_FRAC:
  492     case t_FFELT: case t_PADIC: case t_QUAD:
  493     case t_QFR: case t_QFI: return 0;
  494 
  495     case t_RFRAC: case t_POLMOD:
  496       return isinexactreal(gel(x,1)) || isinexactreal(gel(x,2));
  497 
  498     case t_POL: case t_SER:
  499       for (i=lg(x)-1; i>1; i--)
  500         if (isinexactreal(gel(x,i))) return 1;
  501       return 0;
  502 
  503     case t_VEC: case t_COL: case t_MAT:
  504       for (i=lg(x)-1; i>0; i--)
  505         if (isinexactreal(gel(x,i))) return 1;
  506       return 0;
  507     default: return 0;
  508   }
  509 }
  510 /* Check if x is approximately real with precision e */
  511 int
  512 isrealappr(GEN x, long e)
  513 {
  514   long i;
  515   switch(typ(x))
  516   {
  517     case t_INT: case t_REAL: case t_FRAC:
  518       return 1;
  519     case t_COMPLEX:
  520       return (gexpo(gel(x,2)) < e);
  521 
  522     case t_POL: case t_SER:
  523       for (i=lg(x)-1; i>1; i--)
  524         if (! isrealappr(gel(x,i),e)) return 0;
  525       return 1;
  526 
  527     case t_RFRAC: case t_POLMOD:
  528       return isrealappr(gel(x,1),e) && isrealappr(gel(x,2),e);
  529 
  530     case t_VEC: case t_COL: case t_MAT:
  531       for (i=lg(x)-1; i>0; i--)
  532         if (! isrealappr(gel(x,i),e)) return 0;
  533       return 1;
  534     default: pari_err_TYPE("isrealappr",x); return 0;
  535   }
  536 }
  537 
  538 /* returns 1 if there's an inexact component in the structure, and
  539  * 0 otherwise. */
  540 int
  541 isinexact(GEN x)
  542 {
  543   long lx, i;
  544 
  545   switch(typ(x))
  546   {
  547     case t_REAL: case t_PADIC: case t_SER:
  548       return 1;
  549     case t_INT: case t_INTMOD: case t_FFELT: case t_FRAC:
  550     case t_QFR: case t_QFI:
  551       return 0;
  552     case t_COMPLEX: case t_QUAD: case t_RFRAC: case t_POLMOD:
  553       return isinexact(gel(x,1)) || isinexact(gel(x,2));
  554     case t_POL:
  555       for (i=lg(x)-1; i>1; i--)
  556         if (isinexact(gel(x,i))) return 1;
  557       return 0;
  558     case t_VEC: case t_COL: case t_MAT:
  559       for (i=lg(x)-1; i>0; i--)
  560         if (isinexact(gel(x,i))) return 1;
  561       return 0;
  562     case t_LIST:
  563       x = list_data(x); lx = x? lg(x): 1;
  564       for (i=1; i<lx; i++)
  565         if (isinexact(gel(x,i))) return 1;
  566       return 0;
  567   }
  568   return 0;
  569 }
  570 
  571 int
  572 isrationalzeroscalar(GEN g)
  573 {
  574   switch (typ(g))
  575   {
  576     case t_INT:     return !signe(g);
  577     case t_COMPLEX: return isintzero(gel(g,1)) && isintzero(gel(g,2));
  578     case t_QUAD:    return isintzero(gel(g,2)) && isintzero(gel(g,3));
  579   }
  580   return 0;
  581 }
  582 
  583 int
  584 iscomplex(GEN x)
  585 {
  586   switch(typ(x))
  587   {
  588     case t_INT: case t_REAL: case t_FRAC:
  589       return 0;
  590     case t_COMPLEX:
  591       return !gequal0(gel(x,2));
  592     case t_QUAD:
  593       return signe(gmael(x,1,2)) > 0;
  594   }
  595   pari_err_TYPE("iscomplex",x);
  596   return 0; /* LCOV_EXCL_LINE */
  597 }
  598 
  599 /*******************************************************************/
  600 /*                                                                 */
  601 /*                    GENERIC REMAINDER                            */
  602 /*                                                                 */
  603 /*******************************************************************/
  604 static int
  605 is_realquad(GEN x) { GEN Q = gel(x,1); return signe(gel(Q,2)) < 0; }
  606 static int
  607 is_realext(GEN x)
  608 { long t = typ(x);
  609   return (t == t_QUAD)? is_realquad(x): is_real_t(t);
  610 }
  611 /* euclidean quotient for scalars of admissible types */
  612 static GEN
  613 _quot(GEN x, GEN y)
  614 {
  615   GEN q = gdiv(x,y), f = gfloor(q);
  616   if (gsigne(y) < 0 && !gequal(f,q)) f = addiu(f, 1);
  617   return f;
  618 }
  619 /* y t_REAL, x \ y */
  620 static GEN
  621 _quotsr(long x, GEN y)
  622 {
  623   GEN q, f;
  624   if (!x) return gen_0;
  625   q = divsr(x,y); f = floorr(q);
  626   if (signe(y) < 0 && signe(subir(f,q))) f = addiu(f, 1);
  627   return f;
  628 }
  629 /* x t_REAL, x \ y */
  630 static GEN
  631 _quotrs(GEN x, long y)
  632 {
  633   GEN q = divrs(x,y), f = floorr(q);
  634   if (y < 0 && signe(subir(f,q))) f = addiu(f, 1);
  635   return f;
  636 }
  637 static GEN
  638 _quotri(GEN x, GEN y)
  639 {
  640   GEN q = divri(x,y), f = floorr(q);
  641   if (signe(y) < 0 && signe(subir(f,q))) f = addiu(f, 1);
  642   return f;
  643 }
  644 static GEN
  645 _quotsq(long x, GEN y)
  646 {
  647   GEN f = gfloor(gdivsg(x,y));
  648   if (gsigne(y) < 0) f = gaddgs(f, 1);
  649   return f;
  650 }
  651 static GEN
  652 _quotqs(GEN x, long y)
  653 {
  654   GEN f = gfloor(gdivgs(x,y));
  655   if (y < 0) f = addiu(f, 1);
  656   return f;
  657 }
  658 
  659 /* y t_FRAC, x \ y */
  660 static GEN
  661 _quotsf(long x, GEN y)
  662 { return truedivii(mulis(gel(y,2),x), gel(y,1)); }
  663 /* x t_FRAC, x \ y */
  664 static GEN
  665 _quotfs(GEN x, long y)
  666 { return truedivii(gel(x,1),mulis(gel(x,2),y)); }
  667 /* x t_FRAC, y t_INT, x \ y */
  668 static GEN
  669 _quotfi(GEN x, GEN y)
  670 { return truedivii(gel(x,1),mulii(gel(x,2),y)); }
  671 
  672 static GEN
  673 quot(GEN x, GEN y)
  674 { pari_sp av = avma; return gerepileupto(av, _quot(x, y)); }
  675 static GEN
  676 quotrs(GEN x, long y)
  677 { pari_sp av = avma; return gerepileuptoleaf(av, _quotrs(x,y)); }
  678 static GEN
  679 quotfs(GEN x, long s)
  680 { pari_sp av = avma; return gerepileuptoleaf(av, _quotfs(x,s)); }
  681 static GEN
  682 quotsr(long x, GEN y)
  683 { pari_sp av = avma; return gerepileuptoleaf(av, _quotsr(x, y)); }
  684 static GEN
  685 quotsf(long x, GEN y)
  686 { pari_sp av = avma; return gerepileuptoleaf(av, _quotsf(x, y)); }
  687 static GEN
  688 quotsq(long x, GEN y)
  689 { pari_sp av = avma; return gerepileuptoleaf(av, _quotsq(x, y)); }
  690 static GEN
  691 quotqs(GEN x, long y)
  692 { pari_sp av = avma; return gerepileuptoleaf(av, _quotqs(x, y)); }
  693 static GEN
  694 quotfi(GEN x, GEN y)
  695 { pari_sp av = avma; return gerepileuptoleaf(av, _quotfi(x, y)); }
  696 static GEN
  697 quotri(GEN x, GEN y)
  698 { pari_sp av = avma; return gerepileuptoleaf(av, _quotri(x, y)); }
  699 
  700 static GEN
  701 modrs(GEN x, long y)
  702 {
  703   pari_sp av = avma;
  704   GEN q = _quotrs(x,y);
  705   if (!signe(q)) { set_avma(av); return rcopy(x); }
  706   return gerepileuptoleaf(av, subri(x, mulis(q,y)));
  707 }
  708 static GEN
  709 modsr(long x, GEN y)
  710 {
  711   pari_sp av = avma;
  712   GEN q = _quotsr(x,y);
  713   if (!signe(q)) { set_avma(av); return stoi(x); }
  714   return gerepileuptoleaf(av, subsr(x, mulir(q,y)));
  715 }
  716 static GEN
  717 modsf(long x, GEN y)
  718 {
  719   pari_sp av = avma;
  720   return gerepileupto(av, Qdivii(modii(mulis(gel(y,2),x), gel(y,1)), gel(y,2)));
  721 }
  722 
  723 /* assume y a t_REAL, x a t_INT, t_FRAC or t_REAL.
  724  * Return x mod y or NULL if accuracy error */
  725 GEN
  726 modRr_safe(GEN x, GEN y)
  727 {
  728   GEN q, f;
  729   long e;
  730   if (isintzero(x)) return gen_0;
  731   q = gdiv(x,y); /* t_REAL */
  732 
  733   e = expo(q);
  734   if (e >= 0 && nbits2prec(e+1) > realprec(q)) return NULL;
  735   f = floorr(q);
  736   if (signe(y) < 0 && signe(subri(q,f))) f = addiu(f, 1);
  737   return signe(f)? gsub(x, mulir(f,y)): x;
  738 }
  739 
  740 GEN
  741 gmod(GEN x, GEN y)
  742 {
  743   pari_sp av;
  744   long i, lx, ty, tx;
  745   GEN z;
  746 
  747   tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gmodsg(itos(x),y);
  748   ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gmodgs(x,itos(y));
  749   if (is_matvec_t(tx))
  750   {
  751     z = cgetg_copy(x, &lx);
  752     for (i=1; i<lx; i++) gel(z,i) = gmod(gel(x,i),y);
  753     return z;
  754   }
  755   if (tx == t_POL || ty == t_POL) return grem(x,y);
  756   if (!is_scalar_t(tx) || !is_scalar_t(ty)) pari_err_TYPE2("%",x,y);
  757   switch(ty)
  758   {
  759     case t_INT:
  760       switch(tx)
  761       {
  762         case t_INT: return modii(x,y);
  763         case t_INTMOD: z=cgetg(3, t_INTMOD);
  764           gel(z,1) = gcdii(gel(x,1),y);
  765           gel(z,2) = modii(gel(x,2),gel(z,1)); return z;
  766         case t_FRAC: return Fp_div(gel(x,1),gel(x,2),y);
  767         case t_PADIC: return padic_to_Fp(x, y);
  768         case t_QUAD: if (!is_realquad(x)) break;
  769         case t_REAL:
  770           av = avma; /* NB: conflicting semantic with lift(x * Mod(1,y)). */
  771           return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
  772       }
  773       break;
  774     case t_QUAD:
  775       if (!is_realquad(y)) break;
  776     case t_REAL: case t_FRAC:
  777       if (!is_realext(x)) break;
  778       av = avma;
  779       return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
  780   }
  781   pari_err_TYPE2("%",x,y);
  782   return NULL; /* LCOV_EXCL_LINE */
  783 }
  784 
  785 GEN
  786 gmodgs(GEN x, long y)
  787 {
  788   ulong u;
  789   long i, lx, tx = typ(x);
  790   GEN z;
  791   if (is_matvec_t(tx))
  792   {
  793     z = cgetg_copy(x, &lx);
  794     for (i=1; i<lx; i++) gel(z,i) = gmodgs(gel(x,i),y);
  795     return z;
  796   }
  797   if (!y) pari_err_INV("gmodgs",gen_0);
  798   switch(tx)
  799   {
  800     case t_INT: return modis(x,y);
  801     case t_REAL: return modrs(x,y);
  802 
  803     case t_INTMOD: z=cgetg(3, t_INTMOD);
  804       u = (ulong)labs(y);
  805       i = ugcdiu(gel(x,1), u);
  806       gel(z,1) = utoi(i);
  807       gel(z,2) = modis(gel(x,2), i); return z;
  808 
  809     case t_FRAC:
  810       u = (ulong)labs(y);
  811       return utoi( Fl_div(umodiu(gel(x,1), u),
  812                           umodiu(gel(x,2), u), u) );
  813     case t_QUAD:
  814     {
  815       pari_sp av = avma;
  816       if (!is_realquad(x)) break;
  817       return gerepileupto(av, gsub(x, gmulgs(_quotqs(x,y),y)));
  818     }
  819     case t_PADIC: return padic_to_Fp(x, stoi(y));
  820     case t_POL: return scalarpol(Rg_get_0(x), varn(x));
  821     case t_POLMOD: return gmul(gen_0,x);
  822   }
  823   pari_err_TYPE2("%",x,stoi(y));
  824   return NULL; /* LCOV_EXCL_LINE */
  825 }
  826 GEN
  827 gmodsg(long x, GEN y)
  828 {
  829   switch(typ(y))
  830   {
  831     case t_INT: return modsi(x,y);
  832     case t_REAL: return modsr(x,y);
  833     case t_FRAC: return modsf(x,y);
  834     case t_QUAD:
  835     {
  836       pari_sp av = avma;
  837       if (!is_realquad(y)) break;
  838       return gerepileupto(av, gsubsg(x, gmul(_quotsq(x,y),y)));
  839     }
  840     case t_POL:
  841       if (!signe(y)) pari_err_INV("gmodsg",y);
  842       return degpol(y)? gmulsg(x, Rg_get_1(y)): Rg_get_0(y);
  843   }
  844   pari_err_TYPE2("%",stoi(x),y);
  845   return NULL; /* LCOV_EXCL_LINE */
  846 }
  847 /* divisibility: return 1 if y | x, 0 otherwise */
  848 int
  849 gdvd(GEN x, GEN y)
  850 {
  851   pari_sp av = avma;
  852   return gc_bool(av, gequal0( gmod(x,y) ));
  853 }
  854 
  855 GEN
  856 gmodulss(long x, long y)
  857 {
  858   if (!y) pari_err_INV("%",gen_0);
  859   y = labs(y);
  860   retmkintmod(utoi(umodsu(x, y)), utoipos(y));
  861 }
  862 GEN
  863 gmodulsg(long x, GEN y)
  864 {
  865   switch(typ(y))
  866   {
  867     case t_INT:
  868       if (!is_bigint(y)) return gmodulss(x,itos(y));
  869       retmkintmod(modsi(x,y), absi(y));
  870     case t_POL:
  871       if (!signe(y)) pari_err_INV("%", y);
  872       retmkpolmod(degpol(y)? stoi(x): gen_0,RgX_copy(y));
  873   }
  874   pari_err_TYPE2("%",stoi(x),y);
  875   return NULL; /* LCOV_EXCL_LINE */
  876 }
  877 GEN
  878 gmodulo(GEN x,GEN y)
  879 {
  880   long tx = typ(x), vx, vy;
  881   if (tx == t_INT && !is_bigint(x)) return gmodulsg(itos(x), y);
  882   if (is_matvec_t(tx))
  883   {
  884     long l, i;
  885     GEN z = cgetg_copy(x, &l);
  886     for (i=1; i<l; i++) gel(z,i) = gmodulo(gel(x,i),y);
  887     return z;
  888   }
  889   switch(typ(y))
  890   {
  891     case t_INT:
  892       if (!is_const_t(tx)) return gmul(x, gmodulsg(1,y));
  893       if (tx == t_INTMOD) return gmod(x,y);
  894       retmkintmod(Rg_to_Fp(x,y), absi(y));
  895     case t_POL:
  896       vx = gvar(x); vy = varn(y);
  897       if (varncmp(vy, vx) > 0) return gmul(x, gmodulsg(1,y));
  898       if (vx == vy && tx == t_POLMOD) return grem(x,y);
  899       retmkpolmod(grem(x,y), RgX_copy(y));
  900   }
  901   pari_err_TYPE2("%",x,y);
  902   return NULL; /* LCOV_EXCL_LINE */
  903 }
  904 
  905 /*******************************************************************/
  906 /*                                                                 */
  907 /*                 GENERIC EUCLIDEAN DIVISION                      */
  908 /*                                                                 */
  909 /*******************************************************************/
  910 GEN
  911 gdivent(GEN x, GEN y)
  912 {
  913   long tx, ty;
  914   tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gdiventsg(itos(x),y);
  915   ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gdiventgs(x,itos(y));
  916   if (is_matvec_t(tx))
  917   {
  918     long i, lx;
  919     GEN z = cgetg_copy(x, &lx);
  920     for (i=1; i<lx; i++) gel(z,i) = gdivent(gel(x,i),y);
  921     return z;
  922   }
  923   if (tx == t_POL || ty == t_POL) return gdeuc(x,y);
  924   switch(ty)
  925   {
  926     case t_INT:
  927       switch(tx)
  928       {
  929         case t_INT: return truedivii(x,y);
  930         case t_REAL: return quotri(x,y);
  931         case t_FRAC: return quotfi(x,y);
  932         case t_QUAD:
  933           if (!is_realquad(x)) break;
  934           return quot(x,y);
  935       }
  936       break;
  937     case t_QUAD:
  938       if (!is_realext(x) || !is_realquad(y)) break;
  939     case t_REAL: case t_FRAC:
  940       return quot(x,y);
  941   }
  942   pari_err_TYPE2("\\",x,y);
  943   return NULL; /* LCOV_EXCL_LINE */
  944 }
  945 
  946 GEN
  947 gdiventgs(GEN x, long y)
  948 {
  949   long i, lx;
  950   GEN z;
  951   switch(typ(x))
  952   {
  953     case t_INT:  return truedivis(x,y);
  954     case t_REAL: return quotrs(x,y);
  955     case t_FRAC: return quotfs(x,y);
  956     case t_QUAD: if (!is_realquad(x)) break;
  957                  return quotqs(x,y);
  958     case t_POL:  return gdivgs(x,y);
  959     case t_VEC: case t_COL: case t_MAT:
  960       z = cgetg_copy(x, &lx);
  961       for (i=1; i<lx; i++) gel(z,i) = gdiventgs(gel(x,i),y);
  962       return z;
  963   }
  964   pari_err_TYPE2("\\",x,stoi(y));
  965   return NULL; /* LCOV_EXCL_LINE */
  966 }
  967 GEN
  968 gdiventsg(long x, GEN y)
  969 {
  970   switch(typ(y))
  971   {
  972     case t_INT:  return truedivsi(x,y);
  973     case t_REAL: return quotsr(x,y);
  974     case t_FRAC: return quotsf(x,y);
  975     case t_QUAD: if (!is_realquad(y)) break;
  976                  return quotsq(x,y);
  977     case t_POL:
  978       if (!signe(y)) pari_err_INV("gdiventsg",y);
  979       return degpol(y)? Rg_get_0(y): gdivsg(x,gel(y,2));
  980   }
  981   pari_err_TYPE2("\\",stoi(x),y);
  982   return NULL; /* LCOV_EXCL_LINE */
  983 }
  984 
  985 /* with remainder */
  986 static GEN
  987 quotrem(GEN x, GEN y, GEN *r)
  988 {
  989   GEN q = quot(x,y);
  990   pari_sp av = avma;
  991   *r = gerepileupto(av, gsub(x, gmul(q,y)));
  992   return q;
  993 }
  994 
  995 GEN
  996 gdiventres(GEN x, GEN y)
  997 {
  998   long tx = typ(x), ty = typ(y);
  999   GEN z;
 1000 
 1001   if (is_matvec_t(tx))
 1002   {
 1003     long i, lx;
 1004     z = cgetg_copy(x, &lx);
 1005     for (i=1; i<lx; i++) gel(z,i) = gdiventres(gel(x,i),y);
 1006     return z;
 1007   }
 1008   z = cgetg(3,t_COL);
 1009   if (tx == t_POL || ty == t_POL)
 1010   {
 1011     gel(z,1) = poldivrem(x,y,(GEN*)(z+2));
 1012     return z;
 1013   }
 1014   switch(ty)
 1015   {
 1016     case t_INT:
 1017       switch(tx)
 1018       { /* equal to, but more efficient than next case */
 1019         case t_INT:
 1020           gel(z,1) = truedvmdii(x,y,(GEN*)(z+2));
 1021           return z;
 1022         case t_QUAD:
 1023           if (!is_realquad(x)) break;
 1024         case t_REAL: case t_FRAC:
 1025           gel(z,1) = quotrem(x,y,&gel(z,2));
 1026           return z;
 1027       }
 1028       break;
 1029     case t_QUAD:
 1030       if (!is_realext(x) || !is_realquad(y)) break;
 1031     case t_REAL: case t_FRAC:
 1032       gel(z,1) = quotrem(x,y,&gel(z,2));
 1033       return z;
 1034   }
 1035   pari_err_TYPE2("\\",x,y);
 1036   return NULL; /* LCOV_EXCL_LINE */
 1037 }
 1038 
 1039 GEN
 1040 divrem(GEN x, GEN y, long v)
 1041 {
 1042   pari_sp av = avma;
 1043   long vx, vy;
 1044   GEN q, r;
 1045   if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
 1046   vx = varn(x); if (vx != v) x = swap_vars(x,v);
 1047   vy = varn(y); if (vy != v) y = swap_vars(y,v);
 1048   q = RgX_divrem(x,y, &r);
 1049   if (v && (vx != v || vy != v))
 1050   {
 1051     GEN X = pol_x(v);
 1052     q = gsubst(q, v, X); /* poleval broken for t_RFRAC, subst is safe */
 1053     r = gsubst(r, v, X);
 1054   }
 1055   return gerepilecopy(av, mkcol2(q, r));
 1056 }
 1057 
 1058 GEN
 1059 diviiround(GEN x, GEN y)
 1060 {
 1061   pari_sp av1, av = avma;
 1062   GEN q,r;
 1063   int fl;
 1064 
 1065   q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
 1066   if (r==gen_0) return q;
 1067   av1 = avma;
 1068   fl = abscmpii(shifti(r,1),y);
 1069   set_avma(av1); cgiv(r);
 1070   if (fl >= 0) /* If 2*|r| >= |y| */
 1071   {
 1072     long sz = signe(x)*signe(y);
 1073     if (fl || sz > 0) q = gerepileuptoint(av, addis(q,sz));
 1074   }
 1075   return q;
 1076 }
 1077 
 1078 static GEN
 1079 _abs(GEN x)
 1080 {
 1081   if (typ(x) == t_QUAD) return (gsigne(x) < 0)? gneg(x): x;
 1082   return R_abs_shallow(x);
 1083 }
 1084 
 1085 /* If x and y are not both scalars, same as gdivent.
 1086  * Otherwise, compute the quotient x/y, rounded to the nearest integer
 1087  * (towards +oo in case of tie). */
 1088 GEN
 1089 gdivround(GEN x, GEN y)
 1090 {
 1091   pari_sp av;
 1092   long tx = typ(x), ty = typ(y);
 1093   GEN q, r;
 1094 
 1095   if (tx == t_INT && ty == t_INT) return diviiround(x,y);
 1096   av = avma;
 1097   if (is_realext(x) && is_realext(y))
 1098   { /* same as diviiround, less efficient */
 1099     pari_sp av1;
 1100     int fl;
 1101     q = quotrem(x,y,&r); av1 = avma;
 1102     fl = gcmp(gmul2n(_abs(r),1), _abs(y));
 1103     set_avma(av1); cgiv(r);
 1104     if (fl >= 0) /* If 2*|r| >= |y| */
 1105     {
 1106       long sz = gsigne(y);
 1107       if (fl || sz > 0) q = gerepileupto(av, gaddgs(q, sz));
 1108     }
 1109     return q;
 1110   }
 1111   if (is_matvec_t(tx))
 1112   {
 1113     long i, lx;
 1114     GEN z = cgetg_copy(x, &lx);
 1115     for (i=1; i<lx; i++) gel(z,i) = gdivround(gel(x,i),y);
 1116     return z;
 1117   }
 1118   return gdivent(x,y);
 1119 }
 1120 
 1121 GEN
 1122 gdivmod(GEN x, GEN y, GEN *pr)
 1123 {
 1124   switch(typ(x))
 1125   {
 1126     case t_INT:
 1127       switch(typ(y))
 1128       {
 1129         case t_INT: return dvmdii(x,y,pr);
 1130         case t_POL: *pr=icopy(x); return gen_0;
 1131       }
 1132       break;
 1133     case t_POL: return poldivrem(x,y,pr);
 1134   }
 1135   pari_err_TYPE2("gdivmod",x,y);
 1136   return NULL;/*LCOV_EXCL_LINE*/
 1137 }
 1138 
 1139 /*******************************************************************/
 1140 /*                                                                 */
 1141 /*                               SHIFT                             */
 1142 /*                                                                 */
 1143 /*******************************************************************/
 1144 
 1145 /* Shift tronque si n<0 (multiplication tronquee par 2^n)  */
 1146 
 1147 GEN
 1148 gshift(GEN x, long n)
 1149 {
 1150   long i, lx;
 1151   GEN y;
 1152 
 1153   switch(typ(x))
 1154   {
 1155     case t_INT: return shifti(x,n);
 1156     case t_REAL:return shiftr(x,n);
 1157 
 1158     case t_VEC: case t_COL: case t_MAT:
 1159       y = cgetg_copy(x, &lx);
 1160       for (i=1; i<lx; i++) gel(y,i) = gshift(gel(x,i),n);
 1161       return y;
 1162   }
 1163   return gmul2n(x,n);
 1164 }
 1165 
 1166 /*******************************************************************/
 1167 /*                                                                 */
 1168 /*           SUBSTITUTION DANS UN POLYNOME OU UNE SERIE            */
 1169 /*                                                                 */
 1170 /*******************************************************************/
 1171 
 1172 /* Convert t_SER --> t_POL, ignoring valp. INTERNAL ! */
 1173 GEN
 1174 ser2pol_i(GEN x, long lx)
 1175 {
 1176   long i = lx-1;
 1177   GEN y;
 1178   while (i > 1 && isexactzero(gel(x,i))) i--;
 1179   y = cgetg(i+1, t_POL); y[1] = x[1] & ~VALPBITS;
 1180   for ( ; i > 1; i--) gel(y,i) = gel(x,i);
 1181   return y;
 1182 }
 1183 
 1184 GEN
 1185 ser_inv(GEN b)
 1186 {
 1187   pari_sp av = avma;
 1188   long l = lg(b), e = valp(b), prec = l-2;
 1189   GEN y = RgXn_inv_i(ser2pol_i(b, l), prec);
 1190   GEN x = RgX_to_ser(y, l);
 1191   setvalp(x, -e); return gerepilecopy(av, x);
 1192 }
 1193 
 1194 /* T t_POL in var v, mod out by T components of x which are t_POL/t_RFRAC in v.
 1195  * Recursively. Make sure that resulting polynomials of degree 0 in v are
 1196  * simplified (map K[X]_0 to K) */
 1197 static GEN
 1198 mod_r(GEN x, long v, GEN T)
 1199 {
 1200   long i, w, lx, tx = typ(x);
 1201   GEN y;
 1202 
 1203   if (is_const_t(tx)) return x;
 1204   switch(tx)
 1205   {
 1206     case t_POLMOD:
 1207       w = varn(gel(x,1));
 1208       if (w == v) pari_err_PRIORITY("subst", gel(x,1), "=", v);
 1209       if (varncmp(v, w) < 0) return x;
 1210       return gmodulo(mod_r(gel(x,2),v,T), mod_r(gel(x,1),v,T));
 1211     case t_SER:
 1212       w = varn(x);
 1213       if (w == v) break; /* fail */
 1214       if (varncmp(v, w) < 0 || ser_isexactzero(x)) return x;
 1215       y = cgetg_copy(x, &lx); y[1] = x[1];
 1216       for (i = 2; i < lx; i++) gel(y,i) = mod_r(gel(x,i),v,T);
 1217       return normalize(y);
 1218     case t_POL:
 1219       w = varn(x);
 1220       if (w == v)
 1221       {
 1222         x = RgX_rem(x, T);
 1223         if (!degpol(x)) x = gel(x,2);
 1224         return x;
 1225       }
 1226       if (varncmp(v, w) < 0) return x;
 1227       y = cgetg_copy(x, &lx); y[1] = x[1];
 1228       for (i = 2; i < lx; i++) gel(y,i) = mod_r(gel(x,i),v,T);
 1229       return normalizepol_lg(y, lx);
 1230     case t_RFRAC:
 1231       x = gdiv(mod_r(gel(x,1),v,T), mod_r(gel(x,2),v,T));
 1232       if (typ(x) == t_POL && varn(x) == v && lg(x) == 3) x = gel(x,2);
 1233       return x;
 1234     case t_VEC: case t_COL: case t_MAT:
 1235       y = cgetg_copy(x, &lx);
 1236       for (i = 1; i < lx; i++) gel(y,i) = mod_r(gel(x,i),v,T);
 1237       return y;
 1238     case t_LIST:
 1239       y = mklist();
 1240       list_data(y) = list_data(x)? mod_r(list_data(x),v,T): NULL;
 1241       return y;
 1242   }
 1243   pari_err_TYPE("substpol",x);
 1244   return NULL;/*LCOV_EXCL_LINE*/
 1245 }
 1246 GEN
 1247 gsubstpol(GEN x, GEN T, GEN y)
 1248 {
 1249   pari_sp av = avma;
 1250   long v;
 1251   GEN z;
 1252   if (typ(T) == t_POL && RgX_is_monomial(T) && gequal1(leading_coeff(T)))
 1253   { /* T = t^d */
 1254     long d = degpol(T);
 1255     v = varn(T); z = (d==1)? x: gdeflate(x, v, d);
 1256     if (z) return gerepileupto(av, gsubst(z, v, y));
 1257   }
 1258   v = fetch_var(); T = simplify_shallow(T);
 1259   if (typ(T) == t_RFRAC)
 1260     z = gsub(gel(T,1), gmul(pol_x(v), gel(T,2)));
 1261   else
 1262     z = gsub(T, pol_x(v));
 1263   z = mod_r(x, gvar(T), z);
 1264   z = gsubst(z, v, y); (void)delete_var();
 1265   return gerepileupto(av, z);
 1266 }
 1267 
 1268 long
 1269 RgX_deflate_order(GEN x)
 1270 {
 1271   ulong d = 0, i, lx = (ulong)lg(x);
 1272   for (i=3; i<lx; i++)
 1273     if (!gequal0(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
 1274   return d? (long)d: 1;
 1275 }
 1276 long
 1277 ZX_deflate_order(GEN x)
 1278 {
 1279   ulong d = 0, i, lx = (ulong)lg(x);
 1280   for (i=3; i<lx; i++)
 1281     if (signe(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
 1282   return d? (long)d: 1;
 1283 }
 1284 
 1285 /* deflate (non-leaf) x recursively */
 1286 static GEN
 1287 vdeflate(GEN x, long v, long d)
 1288 {
 1289   long i = lontyp[typ(x)], lx;
 1290   GEN z = cgetg_copy(x, &lx);
 1291   if (i == 2) z[1] = x[1];
 1292   for (; i<lx; i++)
 1293   {
 1294     gel(z,i) = gdeflate(gel(x,i),v,d);
 1295     if (!z[i]) return NULL;
 1296   }
 1297   return z;
 1298 }
 1299 
 1300 /* don't return NULL if substitution fails (fallback won't be able to handle
 1301  * t_SER anyway), fail with a meaningful message */
 1302 static GEN
 1303 serdeflate(GEN x, long v, long d)
 1304 {
 1305   long V, dy, lx, vx = varn(x);
 1306   pari_sp av;
 1307   GEN y;
 1308   if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
 1309   if (varncmp(vx, v) > 0) return gcopy(x);
 1310   av = avma;
 1311   V = valp(x);
 1312   lx = lg(x);
 1313   if (lx == 2) return zeroser(v, V / d);
 1314   y = ser2pol_i(x, lx);
 1315   dy = degpol(y);
 1316   if (V % d != 0 || (dy > 0 && RgX_deflate_order(y) % d != 0))
 1317   {
 1318     const char *s = stack_sprintf("valuation(x) %% %ld", d);
 1319     pari_err_DOMAIN("gdeflate", s, "!=", gen_0,x);
 1320   }
 1321   if (dy > 0) y = RgX_deflate(y, d);
 1322   y = RgX_to_ser(y, 3 + (lx-3)/d);
 1323   setvalp(y, V/d); return gerepilecopy(av, y);
 1324 }
 1325 static GEN
 1326 poldeflate(GEN x, long v, long d)
 1327 {
 1328   long vx = varn(x);
 1329   pari_sp av;
 1330   if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
 1331   if (varncmp(vx, v) > 0 || degpol(x) <= 0) return gcopy(x);
 1332   av = avma;
 1333   /* x nonconstant */
 1334   if (RgX_deflate_order(x) % d != 0) return NULL;
 1335   return gerepilecopy(av, RgX_deflate(x,d));
 1336 }
 1337 static GEN
 1338 listdeflate(GEN x, long v, long d)
 1339 {
 1340   GEN y = NULL, z = mklist();
 1341   if (list_data(x))
 1342   {
 1343     y = vdeflate(list_data(x),v,d);
 1344     if (!y) return NULL;
 1345   }
 1346   list_data(z) = y; return z;
 1347 }
 1348 /* return NULL if substitution fails */
 1349 GEN
 1350 gdeflate(GEN x, long v, long d)
 1351 {
 1352   if (d <= 0) pari_err_DOMAIN("gdeflate", "degree", "<=", gen_0,stoi(d));
 1353   switch(typ(x))
 1354   {
 1355     case t_INT:
 1356     case t_REAL:
 1357     case t_INTMOD:
 1358     case t_FRAC:
 1359     case t_FFELT:
 1360     case t_COMPLEX:
 1361     case t_PADIC:
 1362     case t_QUAD: return gcopy(x);
 1363     case t_POL: return poldeflate(x,v,d);
 1364     case t_SER: return serdeflate(x,v,d);
 1365     case t_POLMOD:
 1366       if (varncmp(varn(gel(x,1)), v) >= 0) return gcopy(x);
 1367       /* fall through */
 1368     case t_RFRAC:
 1369     case t_VEC:
 1370     case t_COL:
 1371     case t_MAT: return vdeflate(x,v,d);
 1372     case t_LIST: return listdeflate(x,v,d);
 1373   }
 1374   pari_err_TYPE("gdeflate",x);
 1375   return NULL; /* LCOV_EXCL_LINE */
 1376 }
 1377 
 1378 /* set *m to the largest d such that x0 = A(X^d); return A */
 1379 GEN
 1380 RgX_deflate_max(GEN x, long *m)
 1381 {
 1382   *m = RgX_deflate_order(x);
 1383   return RgX_deflate(x, *m);
 1384 }
 1385 GEN
 1386 ZX_deflate_max(GEN x, long *m)
 1387 {
 1388   *m = ZX_deflate_order(x);
 1389   return RgX_deflate(x, *m);
 1390 }
 1391 
 1392 static int
 1393 serequalXk(GEN x)
 1394 {
 1395   long i, l = lg(x);
 1396   if (l == 2 || !isint1(gel(x,2))) return 0;
 1397   for (i = 3; i < l; i++)
 1398     if (!isintzero(gel(x,i))) return 0;
 1399   return 1;
 1400 }
 1401 
 1402 static GEN
 1403 gsubst_v(GEN e, long v, GEN x)
 1404 { pari_APPLY_same(gsubst(e, v, gel(x,i))); }
 1405 
 1406 static GEN
 1407 constmat(GEN z, long n)
 1408 {
 1409   GEN y = cgetg(n+1, t_MAT), c = const_col(n, gcopy(z));
 1410   long i;
 1411   for (i = 1; i <= n; i++) gel(y, i) = c;
 1412   return y;
 1413 }
 1414 static GEN
 1415 scalarmat2(GEN o, GEN z, long n)
 1416 {
 1417   GEN y;
 1418   long i;
 1419   if (n == 0) return cgetg(1, t_MAT);
 1420   if (n == 1) retmkmat(mkcol(gcopy(o)));
 1421   y = cgetg(n+1, t_MAT); z = gcopy(z); o = gcopy(o);
 1422   for (i = 1; i <= n; i++) { gel(y, i) = const_col(n, z); gcoeff(y,i,i) = o; }
 1423   return y;
 1424 }
 1425 /* x * y^0, n = dim(y) if t_MAT, else -1 */
 1426 static GEN
 1427 subst_higher(GEN x, GEN y, long n)
 1428 {
 1429   GEN o = Rg_get_1(y);
 1430   if (o == gen_1) return n < 0? gcopy(x): scalarmat(x,n);
 1431   x = gmul(x,o); return n < 0? x: scalarmat2(x, Rg_get_0(y), n);
 1432 }
 1433 
 1434 /* x t_POLMOD, v strictly lower priority than var(x) */
 1435 static GEN
 1436 subst_polmod(GEN x, long v, GEN y)
 1437 {
 1438   long l, i;
 1439   GEN a = gsubst(gel(x,2),v,y), b = gsubst(gel(x,1),v,y), z;
 1440 
 1441   if (typ(b) != t_POL) pari_err_TYPE2("substitution",x,y);
 1442   if (typ(a) != t_POL || varncmp(varn(a), varn(b)) >= 0) return gmodulo(a, b);
 1443   l = lg(a); z = cgetg(l,t_POL); z[1] = a[1];
 1444   for (i = 2; i < l; i++) gel(z,i) = gmodulo(gel(a,i),b);
 1445   return normalizepol_lg(z, l);
 1446 }
 1447 GEN
 1448 gsubst(GEN x, long v, GEN y)
 1449 {
 1450   long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
 1451   long l, vx, vy, ex, ey, i, j, k, jb, matn;
 1452   pari_sp av, av2;
 1453   GEN X, t, z;
 1454 
 1455   switch(ty)
 1456   {
 1457     case t_VEC: case t_COL:
 1458       return gsubst_v(x, v, y);
 1459     case t_MAT:
 1460       if (ly==1) return cgetg(1,t_MAT);
 1461       if (ly == lgcols(y)) { matn = ly - 1; break; }
 1462       /* fall through */
 1463     case t_QFR: case t_QFI:
 1464       pari_err_TYPE2("substitution",x,y);
 1465     default: matn = -1;
 1466   }
 1467   if (is_scalar_t(tx))
 1468   {
 1469     if (tx == t_POLMOD && varncmp(v, varn(gel(x,1))) > 0)
 1470     {
 1471       av = avma;
 1472       return gerepileupto(av, subst_polmod(x, v, y));
 1473     }
 1474     return subst_higher(x, y, matn);
 1475   }
 1476 
 1477   switch(tx)
 1478   {
 1479     case t_POL:
 1480       vx = varn(x);
 1481       if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
 1482       if (varncmp(vx, v) < 0)
 1483       {
 1484         av = avma; z = cgetg(lx, t_POL); z[1] = x[1];
 1485         if (lx == 2) return z;
 1486         for (i = 2; i < lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
 1487         z = normalizepol_lg(z, lx); lx = lg(z);
 1488         if (lx == 2) { set_avma(av); return zeropol(vx); }
 1489         if (lx == 3) return gerepileupto(av, gmul(pol_1(vx), gel(z,2)));
 1490         return gerepileupto(av, poleval(z, pol_x(vx)));
 1491       }
 1492       /* v = vx */
 1493       if (lx == 2)
 1494       {
 1495         GEN z = Rg_get_0(y);
 1496         return matn >= 0? constmat(z, matn): z;
 1497       }
 1498       if (lx == 3)
 1499       {
 1500         x = subst_higher(gel(x,2), y, matn);
 1501         if (matn >= 0) return x;
 1502         vy = gvar(y);
 1503         return (vy == NO_VARIABLE)? x: gmul(x, pol_1(vy));
 1504       }
 1505       return matn >= 0? RgX_RgM_eval(x, y): poleval(x,y);
 1506 
 1507     case t_SER:
 1508       vx = varn(x);
 1509       if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
 1510       ex = valp(x);
 1511       if (varncmp(vx, v) < 0)
 1512       {
 1513         if (lx == 2) return matn >= 0? scalarmat(x, matn): gcopy(x);
 1514         av = avma; X = pol_x(vx);
 1515         av2 = avma;
 1516         z = gadd(gsubst(gel(x,lx-1),v,y), zeroser(vx,1));
 1517         for (i = lx-2; i>=2; i--)
 1518         {
 1519           z = gadd(gmul(z,X), gsubst(gel(x,i),v,y));
 1520           if (gc_needed(av2,1))
 1521           {
 1522             if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
 1523             z = gerepileupto(av2, z);
 1524           }
 1525         }
 1526         if (ex) z = gmul(z, pol_xnall(ex,vx));
 1527         return gerepileupto(av, z);
 1528       }
 1529       switch(ty) /* here vx == v */
 1530       {
 1531         case t_SER:
 1532           vy = varn(y); ey = valp(y);
 1533           if (ey < 1 || lx == 2) return zeroser(vy, ey*(ex+lx-2));
 1534           if (ey == 1 && serequalXk(y)
 1535                       && (varncmp(vx,vy) >= 0 || varncmp(gvar2(x), vy) >= 0))
 1536           { /* y = t + O(t^N) */
 1537             if (lx > ly)
 1538             { /* correct number of significant terms */
 1539               l = ly;
 1540               if (!ex)
 1541                 for (i = 3; i < lx; i++)
 1542                   if (++l >= lx || !gequal0(gel(x,i))) break;
 1543               lx = l;
 1544             }
 1545             z = cgetg(lx, t_SER); z[1] = x[1];
 1546             for (i = 2; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
 1547             if (vx != vy) setvarn(z,vy);
 1548             return z;
 1549           }
 1550           if (vy != vx)
 1551           {
 1552             av = avma; z = gel(x,lx-1);
 1553             for (i=lx-2; i>=2; i--)
 1554             {
 1555               z = gadd(gmul(y,z), gel(x,i));
 1556               if (gc_needed(av,1))
 1557               {
 1558                 if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
 1559                 z = gerepileupto(av, z);
 1560               }
 1561             }
 1562             if (ex) z = gmul(z, gpowgs(y,ex));
 1563             return gerepileupto(av,z);
 1564           }
 1565           l = (lx-2)*ey + 2;
 1566           if (ex) { if (l>ly) l = ly; }
 1567           else if (lx != 3)
 1568           {
 1569             for (i = 3; i < lx; i++)
 1570               if (!isexactzero(gel(x,i))) break;
 1571             l = minss(l, (i-2)*ey + (gequal0(y)? 2 : ly));
 1572           }
 1573           av = avma; t = leafcopy(y);
 1574           if (l < ly) setlg(t, l);
 1575           z = scalarser(gen_1, varn(y), l-2);
 1576           gel(z,2) = gel(x,2); /* ensure lg(z) = l even if x[2] = 0 */
 1577           for (i = 3, jb = ey; jb <= l-2; i++,jb += ey)
 1578           {
 1579             if (i < lx) {
 1580               for (j = jb+2; j < minss(l, jb+ly); j++)
 1581                 gel(z,j) = gadd(gel(z,j), gmul(gel(x,i),gel(t,j-jb)));
 1582             }
 1583             for (j = minss(ly-1, l-1-jb-ey); j > 1; j--)
 1584             {
 1585               GEN a = gmul(gel(t,2), gel(y,j));
 1586               for (k=2; k<j; k++) a = gadd(a, gmul(gel(t,j-k+2), gel(y,k)));
 1587               gel(t,j) = a;
 1588             }
 1589             if (gc_needed(av,1))
 1590             {
 1591               if(DEBUGMEM>1) pari_warn(warnmem,"gsubst");
 1592               gerepileall(av,2, &z,&t);
 1593             }
 1594           }
 1595           if (!ex) return gerepilecopy(av,z);
 1596           return gerepileupto(av, gmul(z, gpowgs(y, ex)));
 1597 
 1598         case t_POL: case t_RFRAC:
 1599         {
 1600           long N, n = lx-2;
 1601           vy = gvar(y); ey = gval(y,vy);
 1602           if (ey == LONG_MAX)
 1603           { /* y = 0 */
 1604             if (ex < 0) pari_err_INV("gsubst",y);
 1605             if (!n) return gcopy(x);
 1606             if (ex > 0) return Rg_get_0(ty == t_RFRAC? gel(y,2): y);
 1607             y = Rg_get_1(ty == t_RFRAC? gel(y,2): y);
 1608             return gmul(y, gel(x,2));
 1609           }
 1610           if (ey < 1 || n == 0) return zeroser(vy, ey*(ex+n));
 1611           av = avma;
 1612           n *= ey;
 1613           N = ex? n: maxss(n-ey,1);
 1614           y = (ty == t_RFRAC)? rfrac_to_ser(y, N+2): RgX_to_ser(y, N+2);
 1615           if (lg(y)-2 > n) setlg(y, n+2);
 1616           x = ser2pol_i(x, lx);
 1617           if (varncmp(vy,vx) > 0)
 1618             z = gadd(poleval(x, y), zeroser(vy,n));
 1619           else
 1620           {
 1621             z = RgXn_eval(x, ser2rfrac_i(y), n);
 1622             if (varn(z) == vy) z = RgX_to_ser(z, n+2);
 1623           }
 1624           switch(typ(z))
 1625           {
 1626             case t_SER:
 1627             case t_POL:
 1628               if (varncmp(varn(z),vy) <= 0) break;
 1629             default: z = scalarser(z, vy, n);
 1630           }
 1631           if (!ex) return gerepilecopy(av, z);
 1632           return gerepileupto(av,  gmul(z, gpowgs(y,ex)));
 1633         }
 1634 
 1635         default:
 1636           if (isexactzero(y))
 1637           {
 1638             if (ex < 0) pari_err_INV("gsubst",y);
 1639             if (ex > 0) return gcopy(y);
 1640             if (lx > 2) return gadd(gel(x,2), y); /*add maps to correct ring*/
 1641           }
 1642           pari_err_TYPE2("substitution",x,y);
 1643       }
 1644       break;
 1645 
 1646     case t_RFRAC:
 1647     {
 1648       GEN a = gel(x,1), b = gel(x,2);
 1649       av = avma;
 1650       a = gsubst(a, v, y);
 1651       b = gsubst(b, v, y); return gerepileupto(av, gdiv(a, b));
 1652     }
 1653 
 1654     case t_VEC: case t_COL: case t_MAT:
 1655       z = cgetg_copy(x, &lx);
 1656       for (i=1; i<lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
 1657       return z;
 1658     case t_LIST:
 1659       z = mklist();
 1660       list_data(z) = list_data(x)? gsubst(list_data(x),v,y): NULL;
 1661       return z;
 1662   }
 1663   return gcopy(x);
 1664 }
 1665 
 1666 /* Return P(x * h), not memory clean */
 1667 GEN
 1668 ser_unscale(GEN P, GEN h)
 1669 {
 1670   long l = lg(P);
 1671   GEN Q = cgetg(l,t_SER);
 1672   Q[1] = P[1];
 1673   if (l != 2)
 1674   {
 1675     long i = 2;
 1676     GEN hi = gpowgs(h, valp(P));
 1677     gel(Q,i) = gmul(gel(P,i), hi);
 1678     for (i++; i<l; i++) { hi = gmul(hi,h); gel(Q,i) = gmul(gel(P,i), hi); }
 1679   }
 1680   return Q;
 1681 }
 1682 
 1683 GEN
 1684 gsubstvec(GEN e, GEN v, GEN r)
 1685 {
 1686   pari_sp av = avma;
 1687   long i, j, l = lg(v);
 1688   GEN w, z, R;
 1689   if ( !is_vec_t(typ(v)) ) pari_err_TYPE("substvec",v);
 1690   if ( !is_vec_t(typ(r)) ) pari_err_TYPE("substvec",r);
 1691   if (lg(r)!=l) pari_err_DIM("substvec");
 1692   w = cgetg(l,t_VECSMALL);
 1693   z = cgetg(l,t_VECSMALL);
 1694   R = cgetg(l,t_VEC);
 1695   for(i=j=1;i<l;i++)
 1696   {
 1697     GEN T = gel(v,i), ri = gel(r,i);
 1698     if (!gequalX(T)) pari_err_TYPE("substvec [not a variable]", T);
 1699     if (gvar(ri) == NO_VARIABLE) /* no need to take precautions */
 1700       e = gsubst(e, varn(T), ri);
 1701     else
 1702     {
 1703       w[j] = varn(T);
 1704       z[j] = fetch_var();
 1705       gel(R,j) = ri;
 1706       j++;
 1707     }
 1708   }
 1709   for(i = 1; i < j; i++) e = gsubst(e,w[i],pol_x(z[i]));
 1710   for(i = 1; i < j; i++) e = gsubst(e,z[i],gel(R,i));
 1711   for(i = 1; i < j; i++) (void)delete_var();
 1712   return gerepileupto(av, e);
 1713 }
 1714 
 1715 /*******************************************************************/
 1716 /*                                                                 */
 1717 /*                SERIE RECIPROQUE D'UNE SERIE                     */
 1718 /*                                                                 */
 1719 /*******************************************************************/
 1720 
 1721 GEN
 1722 serreverse(GEN x)
 1723 {
 1724   long v=varn(x), lx = lg(x), i, mi;
 1725   pari_sp av0 = avma, av;
 1726   GEN a, y, u;
 1727 
 1728   if (typ(x)!=t_SER) pari_err_TYPE("serreverse",x);
 1729   if (valp(x)!=1) pari_err_DOMAIN("serreverse", "valuation", "!=", gen_1,x);
 1730   if (lx < 3) pari_err_DOMAIN("serreverse", "x", "=", gen_0,x);
 1731   y = ser_normalize(x);
 1732   if (y == x) a = NULL; else { a = gel(x,2); x = y; }
 1733   av = avma;
 1734   mi = lx-1; while (mi>=3 && gequal0(gel(x,mi))) mi--;
 1735   u = cgetg(lx,t_SER);
 1736   y = cgetg(lx,t_SER);
 1737   u[1] = y[1] = evalsigne(1) | _evalvalp(1) | evalvarn(v);
 1738   gel(u,2) = gel(y,2) = gen_1;
 1739   if (lx > 3)
 1740   {
 1741     gel(u,3) = gmulsg(-2,gel(x,3));
 1742     gel(y,3) = gneg(gel(x,3));
 1743   }
 1744   for (i=3; i<lx-1; )
 1745   {
 1746     pari_sp av2;
 1747     GEN p1;
 1748     long j, k, K = minss(i,mi);
 1749     for (j=3; j<i+1; j++)
 1750     {
 1751       av2 = avma; p1 = gel(x,j);
 1752       for (k = maxss(3,j+2-mi); k < j; k++)
 1753         p1 = gadd(p1, gmul(gel(u,k),gel(x,j-k+2)));
 1754       p1 = gneg(p1);
 1755       gel(u,j) = gerepileupto(av2, gadd(gel(u,j), p1));
 1756     }
 1757     av2 = avma;
 1758     p1 = gmulsg(i,gel(x,i+1));
 1759     for (k = 2; k < K; k++)
 1760     {
 1761       GEN p2 = gmul(gel(x,k+1),gel(u,i-k+2));
 1762       p1 = gadd(p1, gmulsg(k,p2));
 1763     }
 1764     i++;
 1765     gel(u,i) = gerepileupto(av2, gneg(p1));
 1766     gel(y,i) = gdivgs(gel(u,i), i-1);
 1767     if (gc_needed(av,2))
 1768     {
 1769       GEN dummy = cgetg(1,t_VEC);
 1770       if(DEBUGMEM>1) pari_warn(warnmem,"serreverse");
 1771       for(k=i+1; k<lx; k++) gel(u,k) = gel(y,k) = dummy;
 1772       gerepileall(av,2, &u,&y);
 1773     }
 1774   }
 1775   if (a) y = ser_unscale(y, ginv(a));
 1776   return gerepilecopy(av0,y);
 1777 }
 1778 
 1779 /*******************************************************************/
 1780 /*                                                                 */
 1781 /*                    DERIVATION ET INTEGRATION                    */
 1782 /*                                                                 */
 1783 /*******************************************************************/
 1784 GEN
 1785 derivser(GEN x)
 1786 {
 1787   long i, vx = varn(x), e = valp(x), lx = lg(x);
 1788   GEN y;
 1789   if (ser_isexactzero(x))
 1790   {
 1791     x = gcopy(x);
 1792     if (e) setvalp(x,e-1);
 1793     return x;
 1794   }
 1795   if (e)
 1796   {
 1797     y = cgetg(lx,t_SER); y[1] = evalsigne(1)|evalvalp(e-1) | evalvarn(vx);
 1798     for (i=2; i<lx; i++) gel(y,i) = gmulsg(i+e-2,gel(x,i));
 1799   } else {
 1800     if (lx == 3) return zeroser(vx, 0);
 1801     lx--;
 1802     y = cgetg(lx,t_SER); y[1] = evalsigne(1)|_evalvalp(0) | evalvarn(vx);
 1803     for (i=2; i<lx; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
 1804   }
 1805   return normalize(y);
 1806 }
 1807 
 1808 static GEN
 1809 rfrac_deriv(GEN x, long v)
 1810 {
 1811   pari_sp av = avma;
 1812   GEN y = cgetg(3,t_RFRAC), a = gel(x,1), b = gel(x,2), bp, b0, t, T;
 1813   long vx = varn(b);
 1814 
 1815   bp = deriv(b, v);
 1816   t = simplify_shallow(RgX_gcd(bp, b));
 1817   if (typ(t) != t_POL || varn(t) != vx)
 1818   {
 1819     if (gequal1(t)) b0 = b;
 1820     else
 1821     {
 1822       b0 = RgX_Rg_div(b, t);
 1823       bp = RgX_Rg_div(bp, t);
 1824     }
 1825     a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
 1826     if (isexactzero(a)) return gerepileupto(av, a);
 1827     if (b0 == b)
 1828     {
 1829       gel(y,1) = gerepileupto((pari_sp)y, a);
 1830       gel(y,2) = RgX_sqr(b);
 1831     }
 1832     else
 1833     {
 1834       gel(y,1) = a;
 1835       gel(y,2) = RgX_Rg_mul(RgX_sqr(b0), t);
 1836       y = gerepilecopy(av, y);
 1837     }
 1838     return y;
 1839   }
 1840   b0 = gdivexact(b, t);
 1841   bp = gdivexact(bp,t);
 1842   a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
 1843   if (isexactzero(a)) return gerepileupto(av, a);
 1844   T = RgX_gcd(a, t);
 1845   if (typ(T) != t_POL || varn(T) != vx)
 1846   {
 1847     a = gdiv(a, T);
 1848     t = gdiv(t, T);
 1849   }
 1850   else if (!gequal1(T))
 1851   {
 1852     a = gdivexact(a, T);
 1853     t = gdivexact(t, T);
 1854   }
 1855   gel(y,1) = a;
 1856   gel(y,2) = gmul(RgX_sqr(b0), t);
 1857   return gerepilecopy(av, y);
 1858 }
 1859 
 1860 GEN
 1861 deriv(GEN x, long v)
 1862 {
 1863   long lx, tx, i, j;
 1864   GEN y;
 1865 
 1866   tx = typ(x);
 1867   if (is_const_t(tx))
 1868     switch(tx)
 1869     {
 1870       case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
 1871       case t_FFELT: return FF_zero(x);
 1872       default: return gen_0;
 1873     }
 1874   if (v < 0)
 1875   {
 1876     if (tx == t_CLOSURE) return closure_deriv(x);
 1877     v = gvar9(x);
 1878   }
 1879   switch(tx)
 1880   {
 1881     case t_POLMOD:
 1882     {
 1883       GEN a = gel(x,2), b = gel(x,1);
 1884       if (v == varn(b)) return Rg_get_0(b);
 1885       retmkpolmod(deriv(a,v), RgX_copy(b));
 1886     }
 1887     case t_POL:
 1888       switch(varncmp(varn(x), v))
 1889       {
 1890         case 1: return Rg_get_0(x);
 1891         case 0: return RgX_deriv(x);
 1892       }
 1893       y = cgetg_copy(x, &lx); y[1] = x[1];
 1894       for (i=2; i<lx; i++) gel(y,i) = deriv(gel(x,i),v);
 1895       return normalizepol_lg(y,i);
 1896 
 1897     case t_SER:
 1898       switch(varncmp(varn(x), v))
 1899       {
 1900         case 1: return Rg_get_0(x);
 1901         case 0: return derivser(x);
 1902       }
 1903       if (ser_isexactzero(x)) return gcopy(x);
 1904       y = cgetg_copy(x, &lx); y[1] = x[1];
 1905       for (j=2; j<lx; j++) gel(y,j) = deriv(gel(x,j),v);
 1906       return normalize(y);
 1907 
 1908     case t_RFRAC:
 1909       return rfrac_deriv(x,v);
 1910 
 1911     case t_VEC: case t_COL: case t_MAT:
 1912       y = cgetg_copy(x, &lx);
 1913       for (i=1; i<lx; i++) gel(y,i) = deriv(gel(x,i),v);
 1914       return y;
 1915   }
 1916   pari_err_TYPE("deriv",x);
 1917   return NULL; /* LCOV_EXCL_LINE */
 1918 }
 1919 
 1920 /* n-th derivative of t_SER x, n > 0 */
 1921 static GEN
 1922 derivnser(GEN x, long n)
 1923 {
 1924   long i, vx = varn(x), e = valp(x), lx = lg(x);
 1925   GEN y;
 1926   if (ser_isexactzero(x))
 1927   {
 1928     x = gcopy(x);
 1929     if (e) setvalp(x,e-n);
 1930     return x;
 1931   }
 1932   if (e < 0 || e >= n)
 1933   {
 1934     y = cgetg(lx,t_SER);
 1935     y[1] = evalsigne(1)| evalvalp(e-n) | evalvarn(vx);
 1936     for (i=0; i<lx-2; i++)
 1937       gel(y,i+2) = gmul(muls_interval(i+e-n+1,i+e), gel(x,i+2));
 1938   } else {
 1939     if (lx <= n+2) return zeroser(vx, 0);
 1940     lx -= n;
 1941     y = cgetg(lx,t_SER);
 1942     y[1] = evalsigne(1)|_evalvalp(0) | evalvarn(vx);
 1943     for (i=0; i<lx-2; i++)
 1944       gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n-e));
 1945   }
 1946   return normalize(y);
 1947 }
 1948 
 1949 /* n-th derivative of t_POL x, n > 0 */
 1950 static GEN
 1951 RgX_derivn(GEN x, long n)
 1952 {
 1953   long i, vx = varn(x), lx = lg(x);
 1954   GEN y;
 1955   if (lx <= n+2) return pol_0(vx);
 1956   lx -= n;
 1957   y = cgetg(lx,t_POL);
 1958   y[1] = evalsigne(1)| evalvarn(vx);
 1959   for (i=0; i<lx-2; i++)
 1960     gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n));
 1961   return normalizepol_lg(y, lx);
 1962 }
 1963 
 1964 static GEN
 1965 rfrac_derivn(GEN x, long n, long vs)
 1966 {
 1967   pari_sp av = avma;
 1968   GEN u  = gel(x,1), v = gel(x,2);
 1969   GEN dv = deriv(v, vs);
 1970   long i;
 1971   for (i=1; i<=n; i++)
 1972   {
 1973     GEN du = deriv(u, vs);
 1974     u = gadd(gmul(du, v), gmulsg (-i, gmul(dv, u)));
 1975   }
 1976   v = gpowgs(v, n+1);
 1977   return gerepileupto(av, gdiv(u, v));
 1978 }
 1979 
 1980 GEN
 1981 derivn(GEN x, long n, long v)
 1982 {
 1983   long lx, tx, i, j;
 1984   GEN y;
 1985   if (n < 0)  pari_err_DOMAIN("derivn","n","<", gen_0, stoi(n));
 1986   if (n == 0) return gcopy(x);
 1987   tx = typ(x);
 1988   if (is_const_t(tx))
 1989     switch(tx)
 1990     {
 1991       case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
 1992       case t_FFELT: return FF_zero(x);
 1993       default: return gen_0;
 1994     }
 1995   if (v < 0)
 1996   {
 1997     if (tx == t_CLOSURE) return closure_derivn(x, n);
 1998     v = gvar9(x);
 1999   }
 2000   switch(tx)
 2001   {
 2002     case t_POLMOD:
 2003     {
 2004       GEN a = gel(x,2), b = gel(x,1);
 2005       if (v == varn(b)) return Rg_get_0(b);
 2006       retmkpolmod(derivn(a,n,v), RgX_copy(b));
 2007     }
 2008     case t_POL:
 2009       switch(varncmp(varn(x), v))
 2010       {
 2011         case 1: return Rg_get_0(x);
 2012         case 0: return RgX_derivn(x,n);
 2013       }
 2014       y = cgetg_copy(x, &lx); y[1] = x[1];
 2015       for (i=2; i<lx; i++) gel(y,i) = derivn(gel(x,i),n,v);
 2016       return normalizepol_lg(y,i);
 2017 
 2018     case t_SER:
 2019       switch(varncmp(varn(x), v))
 2020       {
 2021         case 1: return Rg_get_0(x);
 2022         case 0: return derivnser(x, n);
 2023       }
 2024       if (ser_isexactzero(x)) return gcopy(x);
 2025       y = cgetg_copy(x, &lx); y[1] = x[1];
 2026       for (j=2; j<lx; j++) gel(y,j) = derivn(gel(x,j),n,v);
 2027       return normalize(y);
 2028 
 2029     case t_RFRAC:
 2030       return rfrac_derivn(x, n, v);
 2031 
 2032     case t_VEC: case t_COL: case t_MAT:
 2033       y = cgetg_copy(x, &lx);
 2034       for (i=1; i<lx; i++) gel(y,i) = derivn(gel(x,i),n,v);
 2035       return y;
 2036   }
 2037   pari_err_TYPE("derivn",x);
 2038   return NULL; /* LCOV_EXCL_LINE */
 2039 }
 2040 
 2041 static long
 2042 lookup(GEN v, long vx)
 2043 {
 2044   long i ,l = lg(v);
 2045   for(i=1; i<l; i++)
 2046     if (varn(gel(v,i)) == vx) return i;
 2047   return 0;
 2048 }
 2049 
 2050 GEN
 2051 diffop(GEN x, GEN v, GEN dv)
 2052 {
 2053   pari_sp av;
 2054   long i, idx, lx, tx = typ(x), vx;
 2055   GEN y;
 2056   if (!is_vec_t(typ(v))) pari_err_TYPE("diffop",v);
 2057   if (!is_vec_t(typ(dv))) pari_err_TYPE("diffop",dv);
 2058   if (lg(v)!=lg(dv)) pari_err_DIM("diffop");
 2059   if (is_const_t(tx)) return gen_0;
 2060   switch(tx)
 2061   {
 2062     case t_POLMOD:
 2063       av = avma;
 2064       vx  = varn(gel(x,1)); idx = lookup(v,vx);
 2065       if (idx) /*Assume the users now what they are doing */
 2066         y = gmodulo(diffop(gel(x,2),v,dv), gel(x,1));
 2067       else
 2068       {
 2069         GEN m = gel(x,1), pol=gel(x,2);
 2070         GEN u = gneg(gdiv(diffop(m,v,dv),RgX_deriv(m)));
 2071         y = diffop(pol,v,dv);
 2072         if (typ(pol)==t_POL && varn(pol)==varn(m))
 2073           y = gadd(y, gmul(u,RgX_deriv(pol)));
 2074         y = gmodulo(y, gel(x,1));
 2075       }
 2076       return gerepileupto(av, y);
 2077     case t_POL:
 2078       if (signe(x)==0) return gen_0;
 2079       vx  = varn(x); idx = lookup(v,vx);
 2080       av = avma; lx = lg(x);
 2081       y = diffop(gel(x,lx-1),v,dv);
 2082       for (i=lx-2; i>=2; i--) y = gadd(gmul(y,pol_x(vx)),diffop(gel(x,i),v,dv));
 2083       if (idx) y = gadd(y, gmul(gel(dv,idx),RgX_deriv(x)));
 2084       return gerepileupto(av, y);
 2085 
 2086     case t_SER:
 2087       if (signe(x)==0) return gen_0;
 2088       vx  = varn(x); idx = lookup(v,vx);
 2089       if (!idx) return gen_0;
 2090       av = avma;
 2091       if (ser_isexactzero(x)) y = x;
 2092       else
 2093       {
 2094         y = cgetg_copy(x, &lx); y[1] = x[1];
 2095         for (i=2; i<lx; i++) gel(y,i) = diffop(gel(x,i),v,dv);
 2096         y = normalize(y); /* y is probably invalid */
 2097         y = gsubst(y, vx, pol_x(vx)); /* Fix that */
 2098       }
 2099       y = gadd(y, gmul(gel(dv,idx),derivser(x)));
 2100       return gerepileupto(av, y);
 2101 
 2102     case t_RFRAC: {
 2103       GEN a = gel(x,1), b = gel(x,2), ap, bp;
 2104       av = avma;
 2105       ap = diffop(a, v, dv); bp = diffop(b, v, dv);
 2106       y = gsub(gdiv(ap,b),gdiv(gmul(a,bp),gsqr(b)));
 2107       return gerepileupto(av, y);
 2108     }
 2109 
 2110     case t_VEC: case t_COL: case t_MAT:
 2111       y = cgetg_copy(x, &lx);
 2112       for (i=1; i<lx; i++) gel(y,i) = diffop(gel(x,i),v,dv);
 2113       return y;
 2114 
 2115   }
 2116   pari_err_TYPE("diffop",x);
 2117   return NULL; /* LCOV_EXCL_LINE */
 2118 }
 2119 
 2120 GEN
 2121 diffop0(GEN x, GEN v, GEN dv, long n)
 2122 {
 2123   pari_sp av=avma;
 2124   long i;
 2125   for(i=1; i<=n; i++)
 2126     x = gerepileupto(av, diffop(x,v,dv));
 2127   return x;
 2128 }
 2129 
 2130 /********************************************************************/
 2131 /**                                                                **/
 2132 /**                         TAYLOR SERIES                          **/
 2133 /**                                                                **/
 2134 /********************************************************************/
 2135 /* swap vars (vx,v) in x (assume vx < v, vx main variable in x), then call
 2136  * act(data, v, x). FIXME: use in other places */
 2137 static GEN
 2138 swapvar_act(GEN x, long vx, long v, GEN (*act)(void*, long, GEN), void *data)
 2139 {
 2140   long v0 = fetch_var();
 2141   GEN y = act(data, v, gsubst(x,vx,pol_x(v0)));
 2142   y = gsubst(y,v0,pol_x(vx));
 2143   (void)delete_var(); return y;
 2144 }
 2145 /* x + O(v^data) */
 2146 static GEN
 2147 tayl_act(void *data, long v, GEN x) { return gadd(zeroser(v, (long)data), x); }
 2148 static  GEN
 2149 integ_act(void *data, long v, GEN x) { (void)data; return integ(x,v); }
 2150 
 2151 GEN
 2152 tayl(GEN x, long v, long precS)
 2153 {
 2154   long vx = gvar9(x);
 2155   pari_sp av;
 2156 
 2157   if (varncmp(v, vx) <= 0) return gadd(zeroser(v,precS), x);
 2158   av = avma;
 2159   return gerepileupto(av, swapvar_act(x, vx, v, tayl_act, (void*)precS));
 2160 }
 2161 
 2162 GEN
 2163 ggrando(GEN x, long n)
 2164 {
 2165   long m, v;
 2166 
 2167   switch(typ(x))
 2168   {
 2169   case t_INT:/* bug 3 + O(1) */
 2170     if (signe(x) <= 0) pari_err_DOMAIN("O", "x", "<=", gen_0, x);
 2171     if (!is_pm1(x)) return zeropadic(x,n);
 2172     /* +/-1 = x^0 */
 2173     v = m = 0; break;
 2174   case t_POL:
 2175     if (!signe(x)) pari_err_DOMAIN("O", "x", "=", gen_0, x);
 2176     v = varn(x);
 2177     m = n * RgX_val(x); break;
 2178   case t_RFRAC:
 2179     if (gequal0(gel(x,1))) pari_err_DOMAIN("O", "x", "=", gen_0, x);
 2180     v = gvar(x);
 2181     m = n * gval(x,v); break;
 2182     default:  pari_err_TYPE("O", x);
 2183       v = m = 0; /* LCOV_EXCL_LINE */
 2184   }
 2185   return zeroser(v,m);
 2186 }
 2187 
 2188 /*******************************************************************/
 2189 /*                                                                 */
 2190 /*                    FORMAL INTEGRATION                           */
 2191 /*                                                                 */
 2192 /*******************************************************************/
 2193 
 2194 static GEN
 2195 triv_integ(GEN x, long v)
 2196 {
 2197   long i, lx;
 2198   GEN y = cgetg_copy(x, &lx); y[1] = x[1];
 2199   for (i=2; i<lx; i++) gel(y,i) = integ(gel(x,i),v);
 2200   return y;
 2201 }
 2202 
 2203 GEN
 2204 RgX_integ(GEN x)
 2205 {
 2206   long i, lx = lg(x);
 2207   GEN y;
 2208   if (lx == 2) return RgX_copy(x);
 2209   y = cgetg(lx+1, t_POL); y[1] = x[1]; gel(y,2) = gen_0;
 2210   for (i=3; i<=lx; i++) gel(y,i) = gdivgs(gel(x,i-1),i-2);
 2211   return y;
 2212 }
 2213 
 2214 static void
 2215 err_intformal(GEN x)
 2216 { pari_err_DOMAIN("intformal", "residue(series, pole)", "!=", gen_0, x); }
 2217 
 2218 GEN
 2219 integser(GEN x)
 2220 {
 2221   long i, lx = lg(x), vx = varn(x), e = valp(x);
 2222   GEN y;
 2223   if (lx == 2) return zeroser(vx, e+1);
 2224   y = cgetg(lx, t_SER);
 2225   for (i=2; i<lx; i++)
 2226   {
 2227     long j = i+e-1;
 2228     GEN c = gel(x,i);
 2229     if (j)
 2230       c = gdivgs(c, j);
 2231     else
 2232     { /* should be isexactzero, but try to avoid error */
 2233       if (!gequal0(c)) err_intformal(x);
 2234       c = gen_0;
 2235     }
 2236     gel(y,i) = c;
 2237   }
 2238   y[1] = evalsigne(1) | evalvarn(vx) | evalvalp(e+1); return y;
 2239 }
 2240 
 2241 GEN
 2242 integ(GEN x, long v)
 2243 {
 2244   long lx, tx, i, vx, n;
 2245   pari_sp av = avma;
 2246   GEN y,p1;
 2247 
 2248   tx = typ(x);
 2249   if (v < 0) { v = gvar9(x); if (v == NO_VARIABLE) v = 0; }
 2250   if (is_scalar_t(tx))
 2251   {
 2252     if (tx == t_POLMOD)
 2253     {
 2254       GEN a = gel(x,2), b = gel(x,1);
 2255       vx = varn(b);
 2256       if (varncmp(v, vx) > 0) retmkpolmod(integ(a,v), RgX_copy(b));
 2257       if (v == vx) pari_err_PRIORITY("intformal",x,"=",v);
 2258     }
 2259     return deg1pol(x, gen_0, v);
 2260   }
 2261 
 2262   switch(tx)
 2263   {
 2264     case t_POL:
 2265       vx = varn(x);
 2266       if (v == vx) return RgX_integ(x);
 2267       if (lg(x) == 2) {
 2268         if (varncmp(vx, v) < 0) v = vx;
 2269         return zeropol(v);
 2270       }
 2271       if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
 2272       return triv_integ(x,v);
 2273 
 2274     case t_SER:
 2275       vx = varn(x);
 2276       if (v == vx) return integser(x);
 2277       if (lg(x) == 2) {
 2278         if (varncmp(vx, v) < 0) v = vx;
 2279         return zeroser(v, valp(x));
 2280       }
 2281       if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
 2282       return triv_integ(x,v);
 2283 
 2284     case t_RFRAC:
 2285     {
 2286       GEN a = gel(x,1), b = gel(x,2), c, d, s;
 2287       vx = varn(b);
 2288       if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
 2289       if (varncmp(vx, v) < 0)
 2290         return gerepileupto(av, swapvar_act(x, vx, v, integ_act, NULL));
 2291 
 2292       n = degpol(b);
 2293       if (typ(a) == t_POL && varn(a) == vx) n += degpol(a);
 2294       y = integ(gadd(x, zeroser(v,n + 2)), v);
 2295       y = gdiv(gtrunc(gmul(b, y)), b);
 2296       if (typ(y) != t_RFRAC) pari_err_BUG("intformal(t_RFRAC)");
 2297       c = gel(y,1); d = gel(y,2);
 2298       s = gsub(gmul(deriv(c,v),d), gmul(c,deriv(d,v)));
 2299       /* (c'd-cd')/d^2 = y' = x = a/b ? */
 2300       if (!gequal(gmul(s,b), gmul(a,gsqr(d)))) err_intformal(x);
 2301       if (typ(y)==t_RFRAC && lg(gel(y,1)) == lg(gel(y,2)))
 2302       {
 2303         GEN p2 = leading_coeff(gel(y,2));
 2304         p1 = gel(y,1);
 2305         if (typ(p1) == t_POL && varn(p1) == vx) p1 = leading_coeff(p1);
 2306         y = gsub(y, gdiv(p1,p2));
 2307       }
 2308       return gerepileupto(av,y);
 2309     }
 2310 
 2311     case t_VEC: case t_COL: case t_MAT:
 2312       y = cgetg_copy(x, &lx);
 2313       for (i=1; i<lg(x); i++) gel(y,i) = integ(gel(x,i),v);
 2314       return y;
 2315   }
 2316   pari_err_TYPE("integ",x);
 2317   return NULL; /* LCOV_EXCL_LINE */
 2318 }
 2319 
 2320 /*******************************************************************/
 2321 /*                                                                 */
 2322 /*                             FLOOR                               */
 2323 /*                                                                 */
 2324 /*******************************************************************/
 2325 
 2326 GEN
 2327 gfloor(GEN x)
 2328 {
 2329   GEN y;
 2330   long i, lx;
 2331 
 2332   switch(typ(x))
 2333   {
 2334     case t_INT: return icopy(x);
 2335     case t_POL: return RgX_copy(x);
 2336     case t_REAL: return floorr(x);
 2337     case t_FRAC: return truedivii(gel(x,1),gel(x,2));
 2338     case t_QUAD:
 2339     {
 2340       pari_sp av = avma;
 2341       GEN Q = gel(x,1), D = quad_disc(x), u, v, b, d, z;
 2342       if (signe(D) < 0) break;
 2343       x = Q_remove_denom(x, &d);
 2344       u = gel(x,2);
 2345       v = gel(x,3); b = gel(Q,3);
 2346       /* x0 = (2u + v*(-b + sqrt(D))) / (2d) */
 2347       z = sqrti(mulii(D, sqri(v)));
 2348       if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
 2349       /* z = floor(v * sqrt(D)) */
 2350       z = addii(subii(shifti(u,1), mulii(v,b)), z);
 2351       d = d? shifti(d,1): gen_2;
 2352       return gerepileuptoint(av, truedivii(z, d));
 2353     }
 2354     case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
 2355     case t_VEC: case t_COL: case t_MAT:
 2356       y = cgetg_copy(x, &lx);
 2357       for (i=1; i<lx; i++) gel(y,i) = gfloor(gel(x,i));
 2358       return y;
 2359   }
 2360   pari_err_TYPE("gfloor",x);
 2361   return NULL; /* LCOV_EXCL_LINE */
 2362 }
 2363 
 2364 GEN
 2365 gfrac(GEN x)
 2366 {
 2367   pari_sp av = avma;
 2368   return gerepileupto(av, gsub(x,gfloor(x)));
 2369 }
 2370 
 2371 /* assume x t_REAL */
 2372 GEN
 2373 ceilr(GEN x) {
 2374   pari_sp av = avma;
 2375   GEN y = floorr(x);
 2376   if (cmpri(x, y)) return gerepileuptoint(av, addui(1,y));
 2377   return y;
 2378 }
 2379 
 2380 GEN
 2381 gceil(GEN x)
 2382 {
 2383   GEN y;
 2384   long i, lx;
 2385   pari_sp av;
 2386 
 2387   switch(typ(x))
 2388   {
 2389     case t_INT: return icopy(x);
 2390     case t_POL: return RgX_copy(x);
 2391     case t_REAL: return ceilr(x);
 2392     case t_FRAC:
 2393       av = avma; y = divii(gel(x,1),gel(x,2));
 2394       if (signe(gel(x,1)) > 0) y = gerepileuptoint(av, addui(1,y));
 2395       return y;
 2396     case t_QUAD:
 2397       if (!is_realquad(x)) break;
 2398       av = avma; return gerepileupto(av, addiu(gfloor(x), 1));
 2399     case t_RFRAC:
 2400       return gdeuc(gel(x,1),gel(x,2));
 2401 
 2402     case t_VEC: case t_COL: case t_MAT:
 2403       y = cgetg_copy(x, &lx);
 2404       for (i=1; i<lx; i++) gel(y,i) = gceil(gel(x,i));
 2405       return y;
 2406   }
 2407   pari_err_TYPE("gceil",x);
 2408   return NULL; /* LCOV_EXCL_LINE */
 2409 }
 2410 
 2411 GEN
 2412 round0(GEN x, GEN *pte)
 2413 {
 2414   if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
 2415   return ground(x);
 2416 }
 2417 
 2418 /* x t_REAL, return q=floor(x+1/2), set e = expo(x-q) */
 2419 static GEN
 2420 round_i(GEN x, long *pe)
 2421 {
 2422   long e;
 2423   GEN B, q,r, m = mantissa_real(x, &e); /* x = m/2^e */
 2424   if (e <= 0)
 2425   {
 2426     if (e) m = shifti(m,-e);
 2427     *pe = -e; return m;
 2428   }
 2429   B = int2n(e-1);
 2430   m = addii(m, B);
 2431   q = shifti(m, -e);
 2432   r = remi2n(m, e);
 2433   /* 2^e (x+1/2) = m = 2^e q + r, sgn(r)=sgn(m), |r|<2^e */
 2434   if (!signe(r))
 2435     *pe = -1;
 2436   else
 2437   {
 2438     if (signe(m) < 0)
 2439     {
 2440       q = subiu(q,1);
 2441       r = addii(r, B);
 2442     }
 2443     else
 2444       r = subii(r, B);
 2445     /* |x - q| = |r| / 2^e */
 2446     *pe = signe(r)? expi(r) - e: -e;
 2447     cgiv(r);
 2448   }
 2449   return q;
 2450 }
 2451 /* assume x a t_REAL */
 2452 GEN
 2453 roundr(GEN x)
 2454 {
 2455   long ex, s = signe(x);
 2456   pari_sp av;
 2457   if (!s || (ex=expo(x)) < -1) return gen_0;
 2458   if (ex == -1) return s>0? gen_1:
 2459                             absrnz_equal2n(x)? gen_0: gen_m1;
 2460   av = avma; x = round_i(x, &ex);
 2461   if (ex >= 0) pari_err_PREC( "roundr (precision loss in truncation)");
 2462   return gerepileuptoint(av, x);
 2463 }
 2464 GEN
 2465 roundr_safe(GEN x)
 2466 {
 2467   long ex, s = signe(x);
 2468   pari_sp av;
 2469 
 2470   if (!s || (ex = expo(x)) < -1) return gen_0;
 2471   if (ex == -1) return s>0? gen_1:
 2472                             absrnz_equal2n(x)? gen_0: gen_m1;
 2473   av = avma; x = round_i(x, &ex);
 2474   return gerepileuptoint(av, x);
 2475 }
 2476 
 2477 GEN
 2478 ground(GEN x)
 2479 {
 2480   GEN y;
 2481   long i, lx;
 2482   pari_sp av;
 2483 
 2484   switch(typ(x))
 2485   {
 2486     case t_INT: return icopy(x);
 2487     case t_INTMOD: return gcopy(x);
 2488     case t_REAL: return roundr(x);
 2489     case t_FRAC: return diviiround(gel(x,1), gel(x,2));
 2490     case t_QUAD:
 2491     {
 2492       GEN Q = gel(x,1), u, v, b, d, z;
 2493       av = avma;
 2494       if (is_realquad(x)) return gerepileupto(av, gfloor(gadd(x, ghalf)));
 2495       u = gel(x,2);
 2496       v = gel(x,3); b = gel(Q,3);
 2497       u = ground(gsub(u, gmul2n(gmul(v,b),-1)));
 2498       v = Q_remove_denom(v, &d);
 2499       if (!d) d = gen_1;
 2500       /* Im x = v sqrt(|D|) / (2d),
 2501        * Im(round(x)) = floor((d + v sqrt(|D|)) / (2d))
 2502        *              = floor(floor(d + v sqrt(|D|)) / (2d)) */
 2503       z = sqrti(mulii(sqri(v), quad_disc(x)));
 2504       if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
 2505       /* z = floor(v * sqrt(|D|)) */
 2506       v = truedivii(addii(z, d), shifti(d,1));
 2507       return gerepilecopy(av, signe(v)? mkcomplex(u,v): u);
 2508     }
 2509     case t_POLMOD:
 2510       retmkpolmod(ground(gel(x,2)), RgX_copy(gel(x,1)));
 2511 
 2512     case t_COMPLEX:
 2513       av = avma; y = cgetg(3, t_COMPLEX);
 2514       gel(y,2) = ground(gel(x,2));
 2515       if (!signe(gel(y,2))) { set_avma(av); return ground(gel(x,1)); }
 2516       gel(y,1) = ground(gel(x,1)); return y;
 2517 
 2518     case t_POL:
 2519       y = cgetg_copy(x, &lx); y[1] = x[1];
 2520       for (i=2; i<lx; i++) gel(y,i) = ground(gel(x,i));
 2521       return normalizepol_lg(y, lx);
 2522     case t_SER:
 2523       if (ser_isexactzero(x)) return gcopy(x);
 2524       y = cgetg_copy(x, &lx); y[1] = x[1];
 2525       for (i=2; i<lx; i++) gel(y,i) = ground(gel(x,i));
 2526       return normalize(y);
 2527     case t_RFRAC:
 2528       av = avma;
 2529       return gerepileupto(av, gdiv(ground(gel(x,1)), ground(gel(x,2))));
 2530     case t_VEC: case t_COL: case t_MAT:
 2531       y = cgetg_copy(x, &lx);
 2532       for (i=1; i<lx; i++) gel(y,i) = ground(gel(x,i));
 2533       return y;
 2534   }
 2535   pari_err_TYPE("ground",x);
 2536   return NULL; /* LCOV_EXCL_LINE */
 2537 }
 2538 
 2539 /* e = number of error bits on integral part */
 2540 GEN
 2541 grndtoi(GEN x, long *e)
 2542 {
 2543   GEN y;
 2544   long i, lx, e1;
 2545   pari_sp av;
 2546 
 2547   *e = -(long)HIGHEXPOBIT;
 2548   switch(typ(x))
 2549   {
 2550     case t_INT: return icopy(x);
 2551     case t_REAL: {
 2552       long ex = expo(x);
 2553       if (!signe(x) || ex < -1) { *e = ex; return gen_0; }
 2554       av = avma; x = round_i(x, e);
 2555       return gerepileuptoint(av, x);
 2556     }
 2557     case t_FRAC: return diviiround(gel(x,1), gel(x,2));
 2558     case t_INTMOD: return gcopy(x);
 2559     case t_QUAD:
 2560       y = ground(x); av = avma;
 2561       *e = gexpo(gsub(x,y)); set_avma(avma); return y;
 2562     case t_COMPLEX:
 2563       av = avma; y = cgetg(3, t_COMPLEX);
 2564       gel(y,2) = grndtoi(gel(x,2), e);
 2565       if (!signe(gel(y,2))) {
 2566         set_avma(av);
 2567         y = grndtoi(gel(x,1), &e1);
 2568       }
 2569       else
 2570         gel(y,1) = grndtoi(gel(x,1), &e1);
 2571       if (e1 > *e) *e = e1;
 2572       return y;
 2573 
 2574     case t_POLMOD:
 2575       retmkpolmod(grndtoi(gel(x,2), e), RgX_copy(gel(x,1)));
 2576 
 2577     case t_POL:
 2578       y = cgetg_copy(x, &lx); y[1] = x[1];
 2579       for (i=2; i<lx; i++)
 2580       {
 2581         gel(y,i) = grndtoi(gel(x,i),&e1);
 2582         if (e1 > *e) *e = e1;
 2583       }
 2584       return normalizepol_lg(y, lx);
 2585     case t_SER:
 2586       if (ser_isexactzero(x)) return gcopy(x);
 2587       y = cgetg_copy(x, &lx); y[1] = x[1];
 2588       for (i=2; i<lx; i++)
 2589       {
 2590         gel(y,i) = grndtoi(gel(x,i),&e1);
 2591         if (e1 > *e) *e = e1;
 2592       }
 2593       return normalize(y);
 2594     case t_RFRAC:
 2595       y = cgetg(3,t_RFRAC);
 2596       gel(y,1) = grndtoi(gel(x,1),&e1); if (e1 > *e) *e = e1;
 2597       gel(y,2) = grndtoi(gel(x,2),&e1); if (e1 > *e) *e = e1;
 2598       return y;
 2599     case t_VEC: case t_COL: case t_MAT:
 2600       y = cgetg_copy(x, &lx);
 2601       for (i=1; i<lx; i++)
 2602       {
 2603         gel(y,i) = grndtoi(gel(x,i),&e1);
 2604         if (e1 > *e) *e = e1;
 2605       }
 2606       return y;
 2607   }
 2608   pari_err_TYPE("grndtoi",x);
 2609   return NULL; /* LCOV_EXCL_LINE */
 2610 }
 2611 
 2612 /* trunc(x * 2^s). lindep() sanity checks rely on this function to return a
 2613  * t_INT or fail when fed a non-t_COMPLEX input; so do not make this one
 2614  * recursive [ or change the lindep call ] */
 2615 GEN
 2616 gtrunc2n(GEN x, long s)
 2617 {
 2618   GEN z;
 2619   switch(typ(x))
 2620   {
 2621     case t_INT:  return shifti(x, s);
 2622     case t_REAL: return trunc2nr(x, s);
 2623     case t_FRAC: {
 2624       pari_sp av;
 2625       GEN a = gel(x,1), b = gel(x,2), q;
 2626       if (s == 0) return divii(a, b);
 2627       av = avma;
 2628       if (s < 0) q = divii(shifti(a, s), b);
 2629       else {
 2630         GEN r;
 2631         q = dvmdii(a, b, &r);
 2632         q = addii(shifti(q,s), divii(shifti(r,s), b));
 2633       }
 2634       return gerepileuptoint(av, q);
 2635     }
 2636     case t_COMPLEX:
 2637       z = cgetg(3, t_COMPLEX);
 2638       gel(z,2) = gtrunc2n(gel(x,2), s);
 2639       if (!signe(gel(z,2))) {
 2640         set_avma((pari_sp)(z + 3));
 2641         return gtrunc2n(gel(x,1), s);
 2642       }
 2643       gel(z,1) = gtrunc2n(gel(x,1), s);
 2644       return z;
 2645     default: pari_err_TYPE("gtrunc2n",x);
 2646       return NULL; /* LCOV_EXCL_LINE */
 2647   }
 2648 }
 2649 
 2650 /* e = number of error bits on integral part */
 2651 GEN
 2652 gcvtoi(GEN x, long *e)
 2653 {
 2654   long tx = typ(x), lx, e1;
 2655   GEN y;
 2656 
 2657   if (tx == t_REAL)
 2658   {
 2659     long ex = expo(x); if (ex < 0) { *e = ex; return gen_0; }
 2660     e1 = ex - bit_prec(x) + 1;
 2661     y = mantissa2nr(x, e1);
 2662     if (e1 <= 0) { pari_sp av = avma; e1 = expo(subri(x,y)); set_avma(av); }
 2663     *e = e1; return y;
 2664   }
 2665   *e = -(long)HIGHEXPOBIT;
 2666   if (is_matvec_t(tx))
 2667   {
 2668     long i;
 2669     y = cgetg_copy(x, &lx);
 2670     for (i=1; i<lx; i++)
 2671     {
 2672       gel(y,i) = gcvtoi(gel(x,i),&e1);
 2673       if (e1 > *e) *e = e1;
 2674     }
 2675     return y;
 2676   }
 2677   return gtrunc(x);
 2678 }
 2679 
 2680 int
 2681 isint(GEN n, GEN *ptk)
 2682 {
 2683   switch(typ(n))
 2684   {
 2685     case t_INT: *ptk = n; return 1;
 2686     case t_REAL: {
 2687       pari_sp av0 = avma;
 2688       GEN z = floorr(n);
 2689       pari_sp av = avma;
 2690       long s = signe(subri(n, z));
 2691       if (s) return gc_bool(av0,0);
 2692       *ptk = z; return gc_bool(av,1);
 2693     }
 2694     case t_FRAC:    return 0;
 2695     case t_COMPLEX: return gequal0(gel(n,2)) && isint(gel(n,1),ptk);
 2696     case t_QUAD:    return gequal0(gel(n,3)) && isint(gel(n,2),ptk);
 2697     default: pari_err_TYPE("isint",n);
 2698              return 0; /* LCOV_EXCL_LINE */
 2699   }
 2700 }
 2701 
 2702 int
 2703 issmall(GEN n, long *ptk)
 2704 {
 2705   pari_sp av = avma;
 2706   GEN z;
 2707   long k;
 2708   if (!isint(n, &z)) return 0;
 2709   k = itos_or_0(z); set_avma(av);
 2710   if (k || lgefint(z) == 2) { *ptk = k; return 1; }
 2711   return 0;
 2712 }
 2713 
 2714 /* smallest integer greater than any incarnations of the real x
 2715  * Avoid "precision loss in truncation" */
 2716 GEN
 2717 ceil_safe(GEN x)
 2718 {
 2719   pari_sp av = avma;
 2720   long e, tx = typ(x);
 2721   GEN y;
 2722 
 2723   if (is_rational_t(tx)) return gceil(x);
 2724   y = gcvtoi(x,&e);
 2725   if (gsigne(x) >= 0)
 2726   {
 2727     if (e < 0) e = 0;
 2728     y = addii(y, int2n(e));
 2729   }
 2730   return gerepileuptoint(av, y);
 2731 }
 2732 /* largest integer smaller than any incarnations of the real x
 2733  * Avoid "precision loss in truncation" */
 2734 GEN
 2735 floor_safe(GEN x)
 2736 {
 2737   pari_sp av = avma;
 2738   long e, tx = typ(x);
 2739   GEN y;
 2740 
 2741   if (is_rational_t(tx)) return gfloor(x);
 2742   y = gcvtoi(x,&e);
 2743   if (gsigne(x) <= 0)
 2744   {
 2745     if (e < 0) e = 0;
 2746     y = subii(y, int2n(e));
 2747   }
 2748   return gerepileuptoint(av, y);
 2749 }
 2750 
 2751 GEN
 2752 ser2rfrac_i(GEN x)
 2753 {
 2754   long e = valp(x);
 2755   GEN a = ser2pol_i(x, lg(x));
 2756   if (e) {
 2757     if (e > 0) a = RgX_shift_shallow(a, e);
 2758     else a = gred_rfrac_simple(a, pol_xn(-e, varn(a)));
 2759   }
 2760   return a;
 2761 }
 2762 
 2763 static GEN
 2764 ser2rfrac(GEN x)
 2765 {
 2766   pari_sp av = avma;
 2767   return gerepilecopy(av, ser2rfrac_i(x));
 2768 }
 2769 
 2770 /* x t_PADIC, truncate to rational (t_INT/t_FRAC) */
 2771 GEN
 2772 padic_to_Q(GEN x)
 2773 {
 2774   GEN u = gel(x,4), p;
 2775   long v;
 2776   if (!signe(u)) return gen_0;
 2777   v = valp(x);
 2778   if (!v) return icopy(u);
 2779   p = gel(x,2);
 2780   if (v>0)
 2781   {
 2782     pari_sp av = avma;
 2783     return gerepileuptoint(av, mulii(u, powiu(p,v)));
 2784   }
 2785   retmkfrac(icopy(u), powiu(p,-v));
 2786 }
 2787 GEN
 2788 padic_to_Q_shallow(GEN x)
 2789 {
 2790   GEN u = gel(x,4), p, q, q2;
 2791   long v;
 2792   if (!signe(u)) return gen_0;
 2793   q = gel(x,3); q2 = shifti(q,-1);
 2794   v = valp(x);
 2795   u = Fp_center_i(u, q, q2);
 2796   if (!v) return u;
 2797   p = gel(x,2);
 2798   if (v>0) return mulii(powiu(p,v), u);
 2799   return mkfrac(u, powiu(p,-v));
 2800 }
 2801 GEN
 2802 QpV_to_QV(GEN v)
 2803 {
 2804   long i, l;
 2805   GEN w = cgetg_copy(v, &l);
 2806   for (i = 1; i < l; i++)
 2807   {
 2808     GEN c = gel(v,i);
 2809     switch(typ(c))
 2810     {
 2811       case t_INT:
 2812       case t_FRAC: break;
 2813       case t_PADIC: c = padic_to_Q_shallow(c); break;
 2814       default: pari_err_TYPE("padic_to_Q", v);
 2815     }
 2816     gel(w,i) = c;
 2817   }
 2818   return w;
 2819 }
 2820 
 2821 GEN
 2822 gtrunc(GEN x)
 2823 {
 2824   switch(typ(x))
 2825   {
 2826     case t_INT: return icopy(x);
 2827     case t_REAL: return truncr(x);
 2828     case t_FRAC: return divii(gel(x,1),gel(x,2));
 2829     case t_PADIC: return padic_to_Q(x);
 2830     case t_POL: return RgX_copy(x);
 2831     case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
 2832     case t_SER: return ser2rfrac(x);
 2833     case t_VEC: case t_COL: case t_MAT:
 2834     {
 2835       long i, lx;
 2836       GEN y = cgetg_copy(x, &lx);
 2837       for (i=1; i<lx; i++) gel(y,i) = gtrunc(gel(x,i));
 2838       return y;
 2839     }
 2840   }
 2841   pari_err_TYPE("gtrunc",x);
 2842   return NULL; /* LCOV_EXCL_LINE */
 2843 }
 2844 
 2845 GEN
 2846 trunc0(GEN x, GEN *pte)
 2847 {
 2848   if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
 2849   return gtrunc(x);
 2850 }
 2851 /*******************************************************************/
 2852 /*                                                                 */
 2853 /*                  CONVERSIONS -->  INT, POL & SER                */
 2854 /*                                                                 */
 2855 /*******************************************************************/
 2856 
 2857 /* return a_(n-1) B^(n-1) + ... + a_0, where B = 2^32.
 2858  * The a_i are 32bits integers */
 2859 GEN
 2860 mkintn(long n, ...)
 2861 {
 2862   va_list ap;
 2863   GEN x, y;
 2864   long i;
 2865 #ifdef LONG_IS_64BIT
 2866   long e = (n&1);
 2867   n = (n+1) >> 1;
 2868 #endif
 2869   va_start(ap,n);
 2870   x = cgetipos(n+2);
 2871   y = int_MSW(x);
 2872   for (i=0; i <n; i++)
 2873   {
 2874 #ifdef LONG_IS_64BIT
 2875     ulong a = (e && !i)? 0: (ulong) va_arg(ap, unsigned int);
 2876     ulong b = (ulong) va_arg(ap, unsigned int);
 2877     *y = (a << 32) | b;
 2878 #else
 2879     *y = (ulong) va_arg(ap, unsigned int);
 2880 #endif
 2881     y = int_precW(y);
 2882   }
 2883   va_end(ap);
 2884   return int_normalize(x, 0);
 2885 }
 2886 
 2887 /* 2^32 a + b */
 2888 GEN
 2889 uu32toi(ulong a, ulong b)
 2890 {
 2891 #ifdef LONG_IS_64BIT
 2892   return utoi((a<<32) | b);
 2893 #else
 2894   return uutoi(a, b);
 2895 #endif
 2896 }
 2897 /* - (2^32 a + b), assume a or b != 0 */
 2898 GEN
 2899 uu32toineg(ulong a, ulong b)
 2900 {
 2901 #ifdef LONG_IS_64BIT
 2902   return utoineg((a<<32) | b);
 2903 #else
 2904   return uutoineg(a, b);
 2905 #endif
 2906 }
 2907 
 2908 /* return a_(n-1) x^(n-1) + ... + a_0 */
 2909 GEN
 2910 mkpoln(long n, ...)
 2911 {
 2912   va_list ap;
 2913   GEN x, y;
 2914   long i, l = n+2;
 2915   va_start(ap,n);
 2916   x = cgetg(l, t_POL); y = x + 2;
 2917   x[1] = evalvarn(0);
 2918   for (i=n-1; i >= 0; i--) gel(y,i) = va_arg(ap, GEN);
 2919   va_end(ap); return normalizepol_lg(x, l);
 2920 }
 2921 
 2922 /* return [a_1, ..., a_n] */
 2923 GEN
 2924 mkvecn(long n, ...)
 2925 {
 2926   va_list ap;
 2927   GEN x;
 2928   long i;
 2929   va_start(ap,n);
 2930   x = cgetg(n+1, t_VEC);
 2931   for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
 2932   va_end(ap); return x;
 2933 }
 2934 
 2935 GEN
 2936 mkcoln(long n, ...)
 2937 {
 2938   va_list ap;
 2939   GEN x;
 2940   long i;
 2941   va_start(ap,n);
 2942   x = cgetg(n+1, t_COL);
 2943   for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
 2944   va_end(ap); return x;
 2945 }
 2946 
 2947 GEN
 2948 mkvecsmalln(long n, ...)
 2949 {
 2950   va_list ap;
 2951   GEN x;
 2952   long i;
 2953   va_start(ap,n);
 2954   x = cgetg(n+1, t_VECSMALL);
 2955   for (i=1; i <= n; i++) x[i] = va_arg(ap, long);
 2956   va_end(ap); return x;
 2957 }
 2958 
 2959 GEN
 2960 scalarpol(GEN x, long v)
 2961 {
 2962   GEN y;
 2963   if (isrationalzero(x)) return zeropol(v);
 2964   y = cgetg(3,t_POL);
 2965   y[1] = gequal0(x)? evalvarn(v)
 2966                    : evalvarn(v) | evalsigne(1);
 2967   gel(y,2) = gcopy(x); return y;
 2968 }
 2969 GEN
 2970 scalarpol_shallow(GEN x, long v)
 2971 {
 2972   GEN y;
 2973   if (isrationalzero(x)) return zeropol(v);
 2974   y = cgetg(3,t_POL);
 2975   y[1] = gequal0(x)? evalvarn(v)
 2976                    : evalvarn(v) | evalsigne(1);
 2977   gel(y,2) = x; return y;
 2978 }
 2979 
 2980 /* x0 + x1*T, do not assume x1 != 0 */
 2981 GEN
 2982 deg1pol(GEN x1, GEN x0,long v)
 2983 {
 2984   GEN x = cgetg(4,t_POL);
 2985   x[1] = evalsigne(1) | evalvarn(v);
 2986   gel(x,2) = x0 == gen_0? x0: gcopy(x0); /* gen_0 frequent */
 2987   gel(x,3) = gcopy(x1); return normalizepol_lg(x,4);
 2988 }
 2989 
 2990 /* same, no copy */
 2991 GEN
 2992 deg1pol_shallow(GEN x1, GEN x0,long v)
 2993 {
 2994   GEN x = cgetg(4,t_POL);
 2995   x[1] = evalsigne(1) | evalvarn(v);
 2996   gel(x,2) = x0;
 2997   gel(x,3) = x1; return normalizepol_lg(x,4);
 2998 }
 2999 
 3000 GEN
 3001 deg2pol_shallow(GEN x2, GEN x1, GEN x0, long v)
 3002 {
 3003   GEN x = cgetg(5,t_POL);
 3004   x[1] = evalsigne(1) | evalvarn(v);
 3005   gel(x,2) = x0;
 3006   gel(x,3) = x1;
 3007   gel(x,4) = x2;
 3008   return normalizepol_lg(x,5);
 3009 }
 3010 
 3011 static GEN
 3012 _gtopoly(GEN x, long v, int reverse)
 3013 {
 3014   long tx = typ(x);
 3015   GEN y;
 3016 
 3017   if (v<0) v = 0;
 3018   switch(tx)
 3019   {
 3020     case t_POL:
 3021       if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
 3022       y = RgX_copy(x); break;
 3023     case t_SER:
 3024       if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
 3025       y = ser2rfrac(x);
 3026       if (typ(y) != t_POL)
 3027         pari_err_DOMAIN("gtopoly", "valuation", "<", gen_0, x);
 3028       break;
 3029     case t_RFRAC:
 3030     {
 3031       GEN a = gel(x,1), b = gel(x,2);
 3032       long vb = varn(b);
 3033       if (varncmp(vb, v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
 3034       if (typ(a) != t_POL || varn(a) != vb) return zeropol(v);
 3035       y = RgX_div(a,b); break;
 3036     }
 3037     case t_VECSMALL: x = zv_to_ZV(x); /* fall through */
 3038     case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
 3039     {
 3040       long j, k, lx = lg(x);
 3041       GEN z;
 3042       if (tx == t_QFR) lx--;
 3043       if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("gtopoly", x, "<=", v);
 3044       y = cgetg(lx+1, t_POL);
 3045       y[1] = evalvarn(v);
 3046       if (reverse) {
 3047         x--;
 3048         for (j=2; j<=lx; j++) gel(y,j) = gel(x,j);
 3049       } else {
 3050         for (j=2, k=lx; k>=2; j++) gel(y,j) = gel(x,--k);
 3051       }
 3052       z = RgX_copy(normalizepol_lg(y,lx+1));
 3053       settyp(y, t_VECSMALL);/* left on stack */
 3054       return z;
 3055     }
 3056     default:
 3057       if (is_scalar_t(tx)) return scalarpol(x,v);
 3058       pari_err_TYPE("gtopoly",x);
 3059       return NULL; /* LCOV_EXCL_LINE */
 3060   }
 3061   setvarn(y,v); return y;
 3062 }
 3063 
 3064 GEN
 3065 gtopolyrev(GEN x, long v) { return _gtopoly(x,v,1); }
 3066 
 3067 GEN
 3068 gtopoly(GEN x, long v) { return _gtopoly(x,v,0); }
 3069 
 3070 static GEN
 3071 gtovecpost(GEN x, long n)
 3072 {
 3073   long i, imax, lx, tx = typ(x);
 3074   GEN y = zerovec(n);
 3075 
 3076   if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,1) = gcopy(x); return y; }
 3077   switch(tx)
 3078   {
 3079     case t_POL:
 3080       lx=lg(x); imax = minss(lx-2, n);
 3081       for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,lx-i));
 3082       return y;
 3083     case t_SER:
 3084       lx=lg(x); imax = minss(lx-2, n); x++;
 3085       for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
 3086       return y;
 3087     case t_QFR: case t_QFI: case t_VEC: case t_COL:
 3088       lx=lg(x); imax = minss(lx-1, n);
 3089       for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
 3090       return y;
 3091     case t_LIST:
 3092       if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
 3093       x = list_data(x); lx = x? lg(x): 1;
 3094       imax = minss(lx-1, n);
 3095       for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
 3096       return y;
 3097     case t_VECSMALL:
 3098       lx=lg(x);
 3099       imax = minss(lx-1, n);
 3100       for (i=1; i<=imax; i++) gel(y,i) = stoi(x[i]);
 3101       return y;
 3102     default: pari_err_TYPE("gtovec",x);
 3103       return NULL; /*LCOV_EXCL_LINE*/
 3104   }
 3105 }
 3106 
 3107 static GEN
 3108 init_vectopre(long a, long n, GEN y, long *imax)
 3109 {
 3110   *imax = minss(a, n);
 3111   return (n == *imax)?  y: y + n - a;
 3112 }
 3113 static GEN
 3114 gtovecpre(GEN x, long n)
 3115 {
 3116   long i, imax, lx, tx = typ(x);
 3117   GEN y = zerovec(n), y0;
 3118 
 3119   if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,n) = gcopy(x); return y; }
 3120   switch(tx)
 3121   {
 3122     case t_POL:
 3123       lx=lg(x);
 3124       y0 = init_vectopre(lx-2, n, y, &imax);
 3125       for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,lx-i));
 3126       return y;
 3127     case t_SER:
 3128       lx=lg(x); x++;
 3129       y0 = init_vectopre(lx-2, n, y, &imax);
 3130       for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
 3131       return y;
 3132     case t_QFR: case t_QFI: case t_VEC: case t_COL:
 3133       lx=lg(x);
 3134       y0 = init_vectopre(lx-1, n, y, &imax);
 3135       for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
 3136       return y;
 3137     case t_LIST:
 3138       if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
 3139       x = list_data(x); lx = x? lg(x): 1;
 3140       y0 = init_vectopre(lx-1, n, y, &imax);
 3141       for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
 3142       return y;
 3143     case t_VECSMALL:
 3144       lx=lg(x);
 3145       y0 = init_vectopre(lx-1, n, y, &imax);
 3146       for (i=1; i<=imax; i++) gel(y0,i) = stoi(x[i]);
 3147       return y;
 3148     default: pari_err_TYPE("gtovec",x);
 3149       return NULL; /*LCOV_EXCL_LINE*/
 3150   }
 3151 }
 3152 GEN
 3153 gtovec0(GEN x, long n)
 3154 {
 3155   if (!n) return gtovec(x);
 3156   if (n > 0) return gtovecpost(x, n);
 3157   return gtovecpre(x, -n);
 3158 }
 3159 
 3160 GEN
 3161 gtovec(GEN x)
 3162 {
 3163   long i, lx, tx = typ(x);
 3164   GEN y;
 3165 
 3166   if (is_scalar_t(tx)) return mkveccopy(x);
 3167   switch(tx)
 3168   {
 3169     case t_POL:
 3170       lx=lg(x); y=cgetg(lx-1,t_VEC);
 3171       for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,lx-i));
 3172       return y;
 3173     case t_SER:
 3174       lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
 3175       for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
 3176       return y;
 3177     case t_RFRAC: return mkveccopy(x);
 3178     case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
 3179       lx=lg(x); y=cgetg(lx,t_VEC);
 3180       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
 3181       return y;
 3182     case t_LIST:
 3183       if (list_typ(x) == t_LIST_MAP) return mapdomain(x);
 3184       x = list_data(x); lx = x? lg(x): 1;
 3185       y = cgetg(lx, t_VEC);
 3186       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
 3187       return y;
 3188     case t_STR:
 3189     {
 3190       char *s = GSTR(x);
 3191       lx = strlen(s)+1; y = cgetg(lx, t_VEC);
 3192       s--;
 3193       for (i=1; i<lx; i++) gel(y,i) = chartoGENstr(s[i]);
 3194       return y;
 3195     }
 3196     case t_VECSMALL:
 3197       return vecsmall_to_vec(x);
 3198     case t_ERROR:
 3199       lx=lg(x); y = cgetg(lx,t_VEC);
 3200       gel(y,1) = errname(x);
 3201       for (i=2; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
 3202       return y;
 3203     default: pari_err_TYPE("gtovec",x);
 3204       return NULL; /*LCOV_EXCL_LINE*/
 3205   }
 3206 }
 3207 
 3208 GEN
 3209 gtovecrev0(GEN x, long n)
 3210 {
 3211   GEN y = gtovec0(x, -n);
 3212   vecreverse_inplace(y);
 3213   return y;
 3214 }
 3215 GEN
 3216 gtovecrev(GEN x) { return gtovecrev0(x, 0); }
 3217 
 3218 GEN
 3219 gtocol0(GEN x, long n)
 3220 {
 3221   GEN y;
 3222   if (!n) return gtocol(x);
 3223   y = gtovec0(x, n);
 3224   settyp(y, t_COL); return y;
 3225 }
 3226 GEN
 3227 gtocol(GEN x)
 3228 {
 3229   long lx, tx, i, j, h;
 3230   GEN y;
 3231   tx = typ(x);
 3232   if (tx != t_MAT) { y = gtovec(x); settyp(y, t_COL); return y; }
 3233   lx = lg(x); if (lx == 1) return cgetg(1, t_COL);
 3234   h = lgcols(x); y = cgetg(h, t_COL);
 3235   for (i = 1 ; i < h; i++) {
 3236     gel(y,i) = cgetg(lx, t_VEC);
 3237     for (j = 1; j < lx; j++) gmael(y,i,j) = gcopy(gcoeff(x,i,j));
 3238   }
 3239   return y;
 3240 }
 3241 
 3242 GEN
 3243 gtocolrev0(GEN x, long n)
 3244 {
 3245   GEN y = gtocol0(x, -n);
 3246   long ly = lg(y), lim = ly>>1, i;
 3247   for (i = 1; i <= lim; i++) swap(gel(y,i), gel(y,ly-i));
 3248   return y;
 3249 }
 3250 GEN
 3251 gtocolrev(GEN x) { return gtocolrev0(x, 0); }
 3252 
 3253 static long
 3254 Itos(GEN x)
 3255 {
 3256    if (typ(x) != t_INT) pari_err_TYPE("vectosmall",x);
 3257    return itos(x);
 3258 }
 3259 
 3260 static GEN
 3261 gtovecsmallpost(GEN x, long n)
 3262 {
 3263   long i, imax, lx;
 3264   GEN y = zero_Flv(n);
 3265 
 3266   switch(typ(x))
 3267   {
 3268     case t_INT:
 3269       y[1] = itos(x); return y;
 3270     case t_POL:
 3271       lx=lg(x); imax = minss(lx-2, n);
 3272       for (i=1; i<=imax; i++) y[i] = Itos(gel(x,lx-i));
 3273       return y;
 3274     case t_SER:
 3275       lx=lg(x); imax = minss(lx-2, n); x++;
 3276       for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
 3277       return y;
 3278     case t_VEC: case t_COL:
 3279       lx=lg(x); imax = minss(lx-1, n);
 3280       for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
 3281       return y;
 3282     case t_LIST:
 3283       x = list_data(x); lx = x? lg(x): 1;
 3284       imax = minss(lx-1, n);
 3285       for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
 3286       return y;
 3287     case t_VECSMALL:
 3288       lx=lg(x);
 3289       imax = minss(lx-1, n);
 3290       for (i=1; i<=imax; i++) y[i] = x[i];
 3291       return y;
 3292     default: pari_err_TYPE("gtovecsmall",x);
 3293       return NULL; /*LCOV_EXCL_LINE*/
 3294   }
 3295 }
 3296 static GEN
 3297 gtovecsmallpre(GEN x, long n)
 3298 {
 3299   long i, imax, lx;
 3300   GEN y = zero_Flv(n), y0;
 3301 
 3302   switch(typ(x))
 3303   {
 3304     case t_INT:
 3305       y[n] = itos(x); return y;
 3306     case t_POL:
 3307       lx=lg(x);
 3308       y0 = init_vectopre(lx-2, n, y, &imax);
 3309       for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,lx-i));
 3310       return y;
 3311     case t_SER:
 3312       lx=lg(x); x++;
 3313       y0 = init_vectopre(lx-2, n, y, &imax);
 3314       for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
 3315       return y;
 3316     case t_VEC: case t_COL:
 3317       lx=lg(x);
 3318       y0 = init_vectopre(lx-1, n, y, &imax);
 3319       for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
 3320       return y;
 3321     case t_LIST:
 3322       x = list_data(x); lx = x? lg(x): 1;
 3323       y0 = init_vectopre(lx-1, n, y, &imax);
 3324       for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
 3325       return y;
 3326     case t_VECSMALL:
 3327       lx=lg(x);
 3328       y0 = init_vectopre(lx-1, n, y, &imax);
 3329       for (i=1; i<=imax; i++) y0[i] = x[i];
 3330       return y;
 3331     default: pari_err_TYPE("gtovecsmall",x);
 3332       return NULL; /*LCOV_EXCL_LINE*/
 3333   }
 3334 }
 3335 
 3336 GEN
 3337 gtovecsmall0(GEN x, long n)
 3338 {
 3339   if (!n) return gtovecsmall(x);
 3340   if (n > 0) return gtovecsmallpost(x, n);
 3341   return gtovecsmallpre(x, -n);
 3342 }
 3343 
 3344 GEN
 3345 gtovecsmall(GEN x)
 3346 {
 3347   GEN V;
 3348   long l, i;
 3349 
 3350   switch(typ(x))
 3351   {
 3352     case t_INT: return mkvecsmall(itos(x));
 3353     case t_STR:
 3354     {
 3355       unsigned char *s = (unsigned char*)GSTR(x);
 3356       l = strlen((const char *)s);
 3357       V = cgetg(l+1, t_VECSMALL);
 3358       s--;
 3359       for (i=1; i<=l; i++) V[i] = (long)s[i];
 3360       return V;
 3361     }
 3362     case t_VECSMALL: return leafcopy(x);
 3363     case t_LIST:
 3364       x = list_data(x);
 3365       if (!x) return cgetg(1, t_VECSMALL);
 3366       /* fall through */
 3367     case t_VEC: case t_COL:
 3368       l = lg(x); V = cgetg(l,t_VECSMALL);
 3369       for(i=1; i<l; i++) V[i] = Itos(gel(x,i));
 3370       return V;
 3371     case t_POL:
 3372       l = lg(x); V = cgetg(l-1,t_VECSMALL);
 3373       for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,l-i));
 3374       return V;
 3375     case t_SER:
 3376       l = lg(x); V = cgetg(l-1,t_VECSMALL); x++;
 3377       for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,i));
 3378       return V;
 3379     default:
 3380       pari_err_TYPE("vectosmall",x);
 3381       return NULL; /* LCOV_EXCL_LINE */
 3382   }
 3383 }
 3384 
 3385 GEN
 3386 compo(GEN x, long n)
 3387 {
 3388   long tx = typ(x);
 3389   ulong l, lx = (ulong)lg(x);
 3390 
 3391   if (!is_recursive_t(tx))
 3392   {
 3393     if (tx == t_VECSMALL)
 3394     {
 3395       if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
 3396       if ((ulong)n >= lx) pari_err_COMPONENT("", ">", utoi(lx-1), stoi(n));
 3397       return stoi(x[n]);
 3398     }
 3399     pari_err_TYPE("component [leaf]", x);
 3400   }
 3401   if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
 3402   if (tx == t_LIST) {
 3403     x = list_data(x); tx = t_VEC;
 3404     lx = (ulong)(x? lg(x): 1);
 3405   }
 3406   l = (ulong)lontyp[tx] + (ulong)n-1; /* beware overflow */
 3407   if (l >= lx) pari_err_COMPONENT("", ">", utoi(lx-lontyp[tx]), stoi(n));
 3408   return gcopy(gel(x,l));
 3409 }
 3410 
 3411 /* assume x a t_POL */
 3412 static GEN
 3413 _polcoef(GEN x, long n, long v)
 3414 {
 3415   long i, w, lx = lg(x), dx = lx-3;
 3416   GEN z;
 3417   if (dx < 0) return gen_0;
 3418   if (v < 0 || v == (w=varn(x)))
 3419     return (n < 0 || n > dx)? gen_0: gel(x,n+2);
 3420   if (varncmp(w,v) > 0) return n? gen_0: x;
 3421   /* w < v */
 3422   z = cgetg(lx, t_POL); z[1] = x[1];
 3423   for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
 3424   z = normalizepol_lg(z, lx);
 3425   switch(lg(z))
 3426   {
 3427     case 2: z = gen_0; break;
 3428     case 3: z = gel(z,2); break;
 3429   }
 3430   return z;
 3431 }
 3432 
 3433 /* assume x a t_SER */
 3434 static GEN
 3435 _sercoef(GEN x, long n, long v)
 3436 {
 3437   long i, w = varn(x), lx = lg(x), dx = lx-3, N;
 3438   GEN z;
 3439   if (v < 0) v = w;
 3440   N = v == w? n - valp(x): n;
 3441   if (dx < 0)
 3442   {
 3443     if (N >= 0) pari_err_DOMAIN("polcoef", "t_SER", "=", x, x);
 3444     return gen_0;
 3445   }
 3446   if (v == w)
 3447   {
 3448     if (N > dx)
 3449       pari_err_DOMAIN("polcoef", "degree", ">", stoi(dx+valp(x)), stoi(n));
 3450     return (N < 0)? gen_0: gel(x,N+2);
 3451   }
 3452   if (varncmp(w,v) > 0) return N? gen_0: x;
 3453   /* w < v */
 3454   z = cgetg(lx, t_SER); z[1] = x[1];
 3455   for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
 3456   return normalize(z);
 3457 }
 3458 
 3459 /* assume x a t_RFRAC(n) */
 3460 static GEN
 3461 _rfraccoef(GEN x, long n, long v)
 3462 {
 3463   GEN P,Q, p = gel(x,1), q = gel(x,2);
 3464   long vp = gvar(p), vq = gvar(q);
 3465   if (v < 0) v = varncmp(vp, vq) < 0? vp: vq;
 3466   P = (vp == v)? p: swap_vars(p, v);
 3467   Q = (vq == v)? q: swap_vars(q, v);
 3468   if (!RgX_is_monomial(Q)) pari_err_TYPE("polcoef", x);
 3469   n += degpol(Q);
 3470   return gdiv(_polcoef(P, n, v), leading_coeff(Q));
 3471 }
 3472 
 3473 GEN
 3474 polcoef_i(GEN x, long n, long v)
 3475 {
 3476   long tx = typ(x);
 3477   switch(tx)
 3478   {
 3479     case t_POL: return _polcoef(x,n,v);
 3480     case t_SER: return _sercoef(x,n,v);
 3481     case t_RFRAC: return _rfraccoef(x,n,v);
 3482   }
 3483   if (!is_scalar_t(tx)) pari_err_TYPE("polcoef", x);
 3484   return n? gen_0: x;
 3485 }
 3486 
 3487 /* with respect to the main variable if v<0, with respect to the variable v
 3488  * otherwise. v ignored if x is not a polynomial/series. */
 3489 GEN
 3490 polcoef(GEN x, long n, long v)
 3491 {
 3492   pari_sp av = avma;
 3493   x = polcoef_i(x,n,v);
 3494   if (x == gen_0) return x;
 3495   return (avma == av)? gcopy(x): gerepilecopy(av, x);
 3496 }
 3497 
 3498 static GEN
 3499 vecdenom(GEN v, long imin, long imax)
 3500 {
 3501   long i = imin;
 3502   GEN s;
 3503   if (imin > imax) return gen_1;
 3504   s = denom_i(gel(v,i));
 3505   for (i++; i<=imax; i++)
 3506   {
 3507     GEN t = denom_i(gel(v,i));
 3508     if (t != gen_1) s = glcm(s,t);
 3509   }
 3510   return s;
 3511 }
 3512 static GEN denompol(GEN x, long v);
 3513 static GEN
 3514 vecdenompol(GEN v, long imin, long imax, long vx)
 3515 {
 3516   long i = imin;
 3517   GEN s;
 3518   if (imin > imax) return gen_1;
 3519   s = denompol(gel(v,i), vx);
 3520   for (i++; i<=imax; i++)
 3521   {
 3522     GEN t = denompol(gel(v,i), vx);
 3523     if (t != gen_1) s = glcm(s,t);
 3524   }
 3525   return s;
 3526 }
 3527 GEN
 3528 denom_i(GEN x)
 3529 {
 3530   switch(typ(x))
 3531   {
 3532     case t_INT:
 3533     case t_REAL:
 3534     case t_INTMOD:
 3535     case t_FFELT:
 3536     case t_PADIC:
 3537     case t_SER:
 3538     case t_VECSMALL: return gen_1;
 3539     case t_FRAC: return gel(x,2);
 3540     case t_COMPLEX: return vecdenom(x,1,2);
 3541     case t_QUAD: return vecdenom(x,2,3);
 3542     case t_POLMOD: return denom_i(gel(x,2));
 3543     case t_RFRAC: return gel(x,2);
 3544     case t_POL: return pol_1(varn(x));
 3545     case t_VEC: case t_COL: case t_MAT: return vecdenom(x, 1, lg(x)-1);
 3546   }
 3547   pari_err_TYPE("denom",x);
 3548   return NULL; /* LCOV_EXCL_LINE */
 3549 }
 3550 /* v has lower (or equal) priority as x's main variable */
 3551 static GEN
 3552 denompol(GEN x, long v)
 3553 {
 3554   long vx, tx = typ(x);
 3555   if (is_scalar_t(tx)) return gen_1;
 3556   switch(typ(x))
 3557   {
 3558     case t_SER:
 3559       if (varn(x) != v) return x;
 3560       vx = valp(x); return vx < 0? pol_xn(-vx, v): pol_1(v);
 3561     case t_RFRAC: x = gel(x,2); return varn(x) == v? x: pol_1(v);
 3562     case t_POL: return pol_1(v);
 3563     case t_VEC: case t_COL: case t_MAT: return vecdenompol(x, 1, lg(x)-1, v);
 3564   }
 3565   pari_err_TYPE("denom",x);
 3566   return NULL; /* LCOV_EXCL_LINE */
 3567 }
 3568 GEN
 3569 denom(GEN x) { pari_sp av = avma; return gerepilecopy(av, denom_i(x)); }
 3570 
 3571 static GEN
 3572 denominator_v(GEN x, long v)
 3573 {
 3574   long v0 = gvar(x);
 3575   GEN d;
 3576   if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
 3577   if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
 3578   d = denompol(x, v0);
 3579   if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
 3580   return d;
 3581 }
 3582 GEN
 3583 denominator(GEN x, GEN D)
 3584 {
 3585   pari_sp av = avma;
 3586   GEN d;
 3587   if (!D) return denom(x);
 3588   if (isint1(D))
 3589   {
 3590     d = Q_denom_safe(x);
 3591     if (!d) { set_avma(av); return gen_1; }
 3592     return gerepilecopy(av, d);
 3593   }
 3594   if (!gequalX(D)) pari_err_TYPE("denominator", D);
 3595   return gerepileupto(av, denominator_v(x, varn(D)));
 3596 }
 3597 GEN
 3598 numerator(GEN x, GEN D)
 3599 {
 3600   pari_sp av = avma;
 3601   long v;
 3602   if (!D) return numer(x);
 3603   if (isint1(D)) return Q_remove_denom(x,NULL);
 3604   if (!gequalX(D)) pari_err_TYPE("numerator", D);
 3605   v = varn(D); /* optimization */
 3606   if (typ(x) == t_RFRAC && varn(gel(x,2)) == v) return gcopy(gel(x,2));
 3607   return gerepileupto(av, gmul(x, denominator_v(x,v)));
 3608 }
 3609 GEN
 3610 content0(GEN x, GEN D)
 3611 {
 3612   pari_sp av = avma;
 3613   long v, v0;
 3614   GEN d;
 3615   if (!D) return content(x);
 3616   if (isint1(D))
 3617   {
 3618     d = Q_content_safe(x);
 3619     return d? d: gen_1;
 3620   }
 3621   if (!gequalX(D)) pari_err_TYPE("content", D);
 3622   v = varn(D);
 3623   v0 = gvar(x); if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
 3624   if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
 3625   d = content(x);
 3626   /* gsubst is needed because of content([x]) = x */
 3627   if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
 3628   return gerepileupto(av, d);
 3629 }
 3630 
 3631 GEN
 3632 numer_i(GEN x)
 3633 {
 3634   switch(typ(x))
 3635   {
 3636     case t_INT:
 3637     case t_REAL:
 3638     case t_INTMOD:
 3639     case t_FFELT:
 3640     case t_PADIC:
 3641     case t_SER:
 3642     case t_VECSMALL:
 3643     case t_POL: return x;
 3644     case t_POLMOD: return mkpolmod(numer_i(gel(x,2)), gel(x,1));
 3645     case t_FRAC:
 3646     case t_RFRAC: return gel(x,1);
 3647     case t_COMPLEX:
 3648     case t_QUAD:
 3649     case t_VEC:
 3650     case t_COL:
 3651     case t_MAT: return gmul(denom_i(x),x);
 3652   }
 3653   pari_err_TYPE("numer",x);
 3654   return NULL; /* LCOV_EXCL_LINE */
 3655 }
 3656 GEN
 3657 numer(GEN x) { pari_sp av = avma; return gerepilecopy(av, numer_i(x)); }
 3658 
 3659 /* Lift only intmods if v does not occur in x, lift with respect to main
 3660  * variable of x if v < 0, with respect to variable v otherwise */
 3661 GEN
 3662 lift0(GEN x, long v)
 3663 {
 3664   long lx, i;
 3665   GEN y;
 3666 
 3667   switch(typ(x))
 3668   {
 3669     case t_INT: return icopy(x);
 3670     case t_INTMOD: return v < 0? icopy(gel(x,2)): gcopy(x);
 3671     case t_POLMOD:
 3672       if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
 3673       y = cgetg(3, t_POLMOD);
 3674       gel(y,1) = lift0(gel(x,1),v);
 3675       gel(y,2) = lift0(gel(x,2),v); return y;
 3676     case t_PADIC: return v < 0? padic_to_Q(x): gcopy(x);
 3677     case t_POL:
 3678       y = cgetg_copy(x, &lx); y[1] = x[1];
 3679       for (i=2; i<lx; i++) gel(y,i) = lift0(gel(x,i), v);
 3680       return normalizepol_lg(y,lx);
 3681     case t_SER:
 3682       if (ser_isexactzero(x))
 3683       {
 3684         if (lg(x) == 2) return gcopy(x);
 3685         return scalarser(lift0(gel(x,2),v), varn(x), valp(x));
 3686       }
 3687       y = cgetg_copy(x, &lx); y[1] = x[1];
 3688       for (i=2; i<lx; i++) gel(y,i) = lift0(gel(x,i), v);
 3689       return normalize(y);
 3690     case t_COMPLEX: case t_QUAD: case t_RFRAC:
 3691     case t_VEC: case t_COL: case t_MAT:
 3692       y = cgetg_copy(x, &lx);
 3693       for (i=1; i<lx; i++) gel(y,i) = lift0(gel(x,i), v);
 3694       return y;
 3695     default: return gcopy(x);
 3696   }
 3697 }
 3698 /* same as lift, shallow */
 3699 GEN
 3700 lift_shallow(GEN x)
 3701 {
 3702   long i, l;
 3703   GEN y;
 3704   switch(typ(x))
 3705   {
 3706     case t_INTMOD: case t_POLMOD: return gel(x,2);
 3707     case t_PADIC: return padic_to_Q(x);
 3708     case t_SER:
 3709       if (ser_isexactzero(x))
 3710       {
 3711         if (lg(x) == 2) return x;
 3712         return scalarser(lift_shallow(gel(x,2)), varn(x), valp(x));
 3713       }
 3714       y = cgetg_copy(x,&l); y[1] = x[1];
 3715       for (i = 2; i < l; i++) gel(y,i) = lift_shallow(gel(x,i));
 3716       return normalize(y);
 3717     case t_POL:
 3718       y = cgetg_copy(x,&l); y[1] = x[1];
 3719       for (i = 2; i < l; i++) gel(y,i) = lift_shallow(gel(x,i));
 3720       return normalizepol(y);
 3721     case t_COMPLEX: case t_QUAD: case t_RFRAC:
 3722     case t_VEC: case t_COL: case t_MAT:
 3723       y = cgetg_copy(x,&l);
 3724       for (i = 1; i < l; i++) gel(y,i) = lift_shallow(gel(x,i));
 3725       return y;
 3726     default: return x;
 3727   }
 3728 }
 3729 GEN
 3730 lift(GEN x) { return lift0(x,-1); }
 3731 
 3732 GEN
 3733 liftall_shallow(GEN x)
 3734 {
 3735   long lx, i;
 3736   GEN y;
 3737 
 3738   switch(typ(x))
 3739   {
 3740     case t_INTMOD: return gel(x,2);
 3741     case t_POLMOD:
 3742       return liftall_shallow(gel(x,2));
 3743     case t_PADIC: return padic_to_Q(x);
 3744     case t_POL:
 3745       y = cgetg_copy(x, &lx); y[1] = x[1];
 3746       for (i=2; i<lx; i++) gel(y,i) = liftall_shallow(gel(x,i));
 3747       return normalizepol_lg(y,lx);
 3748     case t_SER:
 3749       if (ser_isexactzero(x))
 3750       {
 3751         if (lg(x) == 2) return x;
 3752         return scalarser(liftall_shallow(gel(x,2)), varn(x), valp(x));
 3753       }
 3754       y = cgetg_copy(x, &lx); y[1] = x[1];
 3755       for (i=2; i<lx; i++) gel(y,i) = liftall_shallow(gel(x,i));
 3756       return normalize(y);
 3757     case t_COMPLEX: case t_QUAD: case t_RFRAC:
 3758     case t_VEC: case t_COL: case t_MAT:
 3759       y = cgetg_copy(x, &lx);
 3760       for (i=1; i<lx; i++) gel(y,i) = liftall_shallow(gel(x,i));
 3761       return y;
 3762     default: return x;
 3763   }
 3764 }
 3765 GEN
 3766 liftall(GEN x)
 3767 { pari_sp av = avma; return gerepilecopy(av, liftall_shallow(x)); }
 3768 
 3769 GEN
 3770 liftint_shallow(GEN x)
 3771 {
 3772   long lx, i;
 3773   GEN y;
 3774 
 3775   switch(typ(x))
 3776   {
 3777     case t_INTMOD: return gel(x,2);
 3778     case t_PADIC: return padic_to_Q(x);
 3779     case t_POL:
 3780       y = cgetg_copy(x, &lx); y[1] = x[1];
 3781       for (i=2; i<lx; i++) gel(y,i) = liftint_shallow(gel(x,i));
 3782       return normalizepol_lg(y,lx);
 3783     case t_SER:
 3784       if (ser_isexactzero(x))
 3785       {
 3786         if (lg(x) == 2) return x;
 3787         return scalarser(liftint_shallow(gel(x,2)), varn(x), valp(x));
 3788       }
 3789       y = cgetg_copy(x, &lx); y[1] = x[1];
 3790       for (i=2; i<lx; i++) gel(y,i) = liftint_shallow(gel(x,i));
 3791       return normalize(y);
 3792     case t_POLMOD: case t_COMPLEX: case t_QUAD: case t_RFRAC:
 3793     case t_VEC: case t_COL: case t_MAT:
 3794       y = cgetg_copy(x, &lx);
 3795       for (i=1; i<lx; i++) gel(y,i) = liftint_shallow(gel(x,i));
 3796       return y;
 3797     default: return x;
 3798   }
 3799 }
 3800 GEN
 3801 liftint(GEN x)
 3802 { pari_sp av = avma; return gerepilecopy(av, liftint_shallow(x)); }
 3803 
 3804 GEN
 3805 liftpol_shallow(GEN x)
 3806 {
 3807   long lx, i;
 3808   GEN y;
 3809 
 3810   switch(typ(x))
 3811   {
 3812     case t_POLMOD:
 3813       return liftpol_shallow(gel(x,2));
 3814     case t_POL:
 3815       y = cgetg_copy(x, &lx); y[1] = x[1];
 3816       for (i=2; i<lx; i++) gel(y,i) = liftpol_shallow(gel(x,i));
 3817       return normalizepol_lg(y,lx);
 3818     case t_SER:
 3819       if (ser_isexactzero(x))
 3820       {
 3821         if (lg(x) == 2) return x;
 3822         return scalarser(liftpol_shallow(gel(x,2)), varn(x), valp(x));
 3823       }
 3824       y = cgetg_copy(x, &lx); y[1] = x[1];
 3825       for (i=2; i<lx; i++) gel(y,i) = liftpol_shallow(gel(x,i));
 3826       return normalize(y);
 3827     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
 3828       y = cgetg_copy(x, &lx);
 3829       for (i=1; i<lx; i++) gel(y,i) = liftpol_shallow(gel(x,i));
 3830       return y;
 3831     default: return x;
 3832   }
 3833 }
 3834 GEN
 3835 liftpol(GEN x)
 3836 { pari_sp av = avma; return gerepilecopy(av, liftpol_shallow(x)); }
 3837 
 3838 static GEN
 3839 centerliftii(GEN x, GEN y)
 3840 {
 3841   pari_sp av = avma;
 3842   long i = cmpii(shifti(x,1), y);
 3843   set_avma(av); return (i > 0)? subii(x,y): icopy(x);
 3844 }
 3845 
 3846 /* see lift0 */
 3847 GEN
 3848 centerlift0(GEN x, long v)
 3849 { return v < 0? centerlift(x): lift0(x,v); }
 3850 GEN
 3851 centerlift(GEN x)
 3852 {
 3853   long i, v, lx;
 3854   GEN y;
 3855   switch(typ(x))
 3856   {
 3857     case t_INT: return icopy(x);
 3858     case t_INTMOD: return centerliftii(gel(x,2), gel(x,1));
 3859     case t_POLMOD: return gcopy(gel(x,2));
 3860    case t_POL:
 3861       y = cgetg_copy(x, &lx); y[1] = x[1];
 3862       for (i=2; i<lx; i++) gel(y,i) = centerlift(gel(x,i));
 3863       return normalizepol_lg(y,lx);
 3864    case t_SER:
 3865       if (ser_isexactzero(x)) return lift(x);
 3866       y = cgetg_copy(x, &lx); y[1] = x[1];
 3867       for (i=2; i<lx; i++) gel(y,i) = centerlift(gel(x,i));
 3868       return normalize(y);
 3869    case t_COMPLEX: case t_QUAD: case t_RFRAC:
 3870    case t_VEC: case t_COL: case t_MAT:
 3871       y = cgetg_copy(x, &lx);
 3872       for (i=1; i<lx; i++) gel(y,i) = centerlift(gel(x,i));
 3873       return y;
 3874     case t_PADIC:
 3875       if (!signe(gel(x,4))) return gen_0;
 3876       v = valp(x);
 3877       if (v>=0)
 3878       { /* here p^v is an integer */
 3879         GEN z =  centerliftii(gel(x,4), gel(x,3));
 3880         pari_sp av;
 3881         if (!v) return z;
 3882         av = avma; y = powiu(gel(x,2),v);
 3883         return gerepileuptoint(av, mulii(y,z));
 3884       }
 3885       y = cgetg(3,t_FRAC);
 3886       gel(y,1) = centerliftii(gel(x,4), gel(x,3));
 3887       gel(y,2) = powiu(gel(x,2),-v);
 3888       return y;
 3889     default: return gcopy(x);
 3890   }
 3891 }
 3892 
 3893 /*******************************************************************/
 3894 /*                                                                 */
 3895 /*                      REAL & IMAGINARY PARTS                     */
 3896 /*                                                                 */
 3897 /*******************************************************************/
 3898 
 3899 static GEN
 3900 op_ReIm(GEN f(GEN), GEN x)
 3901 {
 3902   long lx, i, j;
 3903   pari_sp av;
 3904   GEN z;
 3905 
 3906   switch(typ(x))
 3907   {
 3908     case t_POL:
 3909       z = cgetg_copy(x, &lx); z[1] = x[1];
 3910       for (j=2; j<lx; j++) gel(z,j) = f(gel(x,j));
 3911       return normalizepol_lg(z, lx);
 3912 
 3913     case t_SER:
 3914       z = cgetg_copy(x, &lx); z[1] = x[1];
 3915       for (j=2; j<lx; j++) gel(z,j) = f(gel(x,j));
 3916       return normalize(z);
 3917 
 3918     case t_RFRAC: {
 3919       GEN dxb, n, d;
 3920       av = avma; dxb = conj_i(gel(x,2));
 3921       n = gmul(gel(x,1), dxb);
 3922       d = gmul(gel(x,2), dxb);
 3923       return gerepileupto(av, gdiv(f(n), d));
 3924     }
 3925 
 3926     case t_VEC: case t_COL: case t_MAT:
 3927       z = cgetg_copy(x, &lx);
 3928       for (i=1; i<lx; i++) gel(z,i) = f(gel(x,i));
 3929       return z;
 3930   }
 3931   pari_err_TYPE("greal/gimag",x);
 3932   return NULL; /* LCOV_EXCL_LINE */
 3933 }
 3934 
 3935 GEN
 3936 real_i(GEN x)
 3937 {
 3938   switch(typ(x))
 3939   {
 3940     case t_INT: case t_REAL: case t_FRAC:
 3941       return x;
 3942     case t_COMPLEX:
 3943       return gel(x,1);
 3944     case t_QUAD:
 3945       return gel(x,2);
 3946   }
 3947   return op_ReIm(real_i,x);
 3948 }
 3949 GEN
 3950 imag_i(GEN x)
 3951 {
 3952   switch(typ(x))
 3953   {
 3954     case t_INT: case t_REAL: case t_FRAC:
 3955       return gen_0;
 3956     case t_COMPLEX:
 3957       return gel(x,2);
 3958     case t_QUAD:
 3959       return gel(x,3);
 3960   }
 3961   return op_ReIm(imag_i,x);
 3962 }
 3963 GEN
 3964 greal(GEN x)
 3965 {
 3966   switch(typ(x))
 3967   {
 3968     case t_INT: case t_REAL:
 3969       return mpcopy(x);
 3970 
 3971     case t_FRAC:
 3972       return gcopy(x);
 3973 
 3974     case t_COMPLEX:
 3975       return gcopy(gel(x,1));
 3976 
 3977     case t_QUAD:
 3978       return gcopy(gel(x,2));
 3979   }
 3980   return op_ReIm(greal,x);
 3981 }
 3982 GEN
 3983 gimag(GEN x)
 3984 {
 3985   switch(typ(x))
 3986   {
 3987     case t_INT: case t_REAL: case t_FRAC:
 3988       return gen_0;
 3989 
 3990     case t_COMPLEX:
 3991       return gcopy(gel(x,2));
 3992 
 3993     case t_QUAD:
 3994       return gcopy(gel(x,3));
 3995   }
 3996   return op_ReIm(gimag,x);
 3997 }
 3998 
 3999 /* return Re(x * y), x and y scalars */
 4000 GEN
 4001 mulreal(GEN x, GEN y)
 4002 {
 4003   if (typ(x) == t_COMPLEX)
 4004   {
 4005     if (typ(y) == t_COMPLEX)
 4006     {
 4007       pari_sp av = avma;
 4008       GEN a = gmul(gel(x,1), gel(y,1));
 4009       GEN b = gmul(gel(x,2), gel(y,2));
 4010       return gerepileupto(av, gsub(a, b));
 4011     }
 4012     x = gel(x,1);
 4013   }
 4014   else
 4015     if (typ(y) == t_COMPLEX) y = gel(y,1);
 4016   return gmul(x,y);
 4017 }
 4018 /* Compute Re(x * y), x and y matrices of compatible dimensions
 4019  * assume scalar entries */
 4020 GEN
 4021 RgM_mulreal(GEN x, GEN y)
 4022 {
 4023   long i, j, k, l, lx = lg(x), ly = lg(y);
 4024   GEN z = cgetg(ly,t_MAT);
 4025   l = (lx == 1)? 1: lgcols(x);
 4026   for (j=1; j<ly; j++)
 4027   {
 4028     GEN zj = cgetg(l,t_COL), yj = gel(y,j);
 4029     gel(z,j) = zj;
 4030     for (i=1; i<l; i++)
 4031     {
 4032       pari_sp av = avma;
 4033       GEN c = mulreal(gcoeff(x,i,1),gel(yj,1));
 4034       for (k=2; k<lx; k++) c = gadd(c, mulreal(gcoeff(x,i,k),gel(yj,k)));
 4035       gel(zj, i) = gerepileupto(av, c);
 4036     }
 4037   }
 4038   return z;
 4039 }
 4040 
 4041 /*******************************************************************/
 4042 /*                                                                 */
 4043 /*                     LOGICAL OPERATIONS                          */
 4044 /*                                                                 */
 4045 /*******************************************************************/
 4046 static long
 4047 _egal_i(GEN x, GEN y)
 4048 {
 4049   x = simplify_shallow(x);
 4050   y = simplify_shallow(y);
 4051   if (typ(y) == t_INT)
 4052   {
 4053     if (equali1(y)) return gequal1(x);
 4054     if (equalim1(y)) return gequalm1(x);
 4055   }
 4056   else if (typ(x) == t_INT)
 4057   {
 4058     if (equali1(x)) return gequal1(y);
 4059     if (equalim1(x)) return gequalm1(y);
 4060   }
 4061   return gequal(x, y);
 4062 }
 4063 static long
 4064 _egal(GEN x, GEN y)
 4065 { pari_sp av = avma; return gc_long(av, _egal_i(x, y)); }
 4066 
 4067 GEN
 4068 glt(GEN x, GEN y) { return gcmp(x,y)<0? gen_1: gen_0; }
 4069 
 4070 GEN
 4071 gle(GEN x, GEN y) { return gcmp(x,y)<=0? gen_1: gen_0; }
 4072 
 4073 GEN
 4074 gge(GEN x, GEN y) { return gcmp(x,y)>=0? gen_1: gen_0; }
 4075 
 4076 GEN
 4077 ggt(GEN x, GEN y) { return gcmp(x,y)>0? gen_1: gen_0; }
 4078 
 4079 GEN
 4080 geq(GEN x, GEN y) { return _egal(x,y)? gen_1: gen_0; }
 4081 
 4082 GEN
 4083 gne(GEN x, GEN y) { return _egal(x,y)? gen_0: gen_1; }
 4084 
 4085 GEN
 4086 gnot(GEN x) { return gequal0(x)? gen_1: gen_0; }
 4087 
 4088 /*******************************************************************/
 4089 /*                                                                 */
 4090 /*                      FORMAL SIMPLIFICATIONS                     */
 4091 /*                                                                 */
 4092 /*******************************************************************/
 4093 
 4094 GEN
 4095 geval_gp(GEN x, GEN t)
 4096 {
 4097   long lx, i, tx = typ(x);
 4098   pari_sp av;
 4099   GEN y, z;
 4100 
 4101   if (is_const_t(tx) || tx==t_VECSMALL) return gcopy(x);
 4102   switch(tx)
 4103   {
 4104     case t_STR:
 4105       return localvars_read_str(GSTR(x),t);
 4106 
 4107     case t_POLMOD:
 4108       av = avma;
 4109       return gerepileupto(av, gmodulo(geval_gp(gel(x,2),t),
 4110                                       geval_gp(gel(x,1),t)));
 4111 
 4112     case t_POL:
 4113       lx=lg(x); if (lx==2) return gen_0;
 4114       z = fetch_var_value(varn(x),t);
 4115       if (!z) return RgX_copy(x);
 4116       av = avma; y = geval_gp(gel(x,lx-1),t);
 4117       for (i=lx-2; i>1; i--)
 4118         y = gadd(geval_gp(gel(x,i),t), gmul(z,y));
 4119       return gerepileupto(av, y);
 4120 
 4121     case t_SER:
 4122       pari_err_IMPL( "evaluation of a power series");
 4123 
 4124     case t_RFRAC:
 4125       av = avma;
 4126       return gerepileupto(av, gdiv(geval_gp(gel(x,1),t), geval_gp(gel(x,2),t)));
 4127 
 4128     case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
 4129       y = cgetg_copy(x, &lx);
 4130       for (i=1; i<lx; i++) gel(y,i) = geval_gp(gel(x,i),t);
 4131       return y;
 4132 
 4133     case t_CLOSURE:
 4134       if (closure_arity(x) || closure_is_variadic(x))
 4135         pari_err_IMPL("eval on functions with parameters");
 4136       return closure_evalres(x);
 4137   }
 4138   pari_err_TYPE("geval",x);
 4139   return NULL; /* LCOV_EXCL_LINE */
 4140 }
 4141 GEN
 4142 geval(GEN x) { return geval_gp(x,NULL); }
 4143 
 4144 GEN
 4145 simplify_shallow(GEN x)
 4146 {
 4147   long i, lx;
 4148   GEN y, z;
 4149   if (!x) pari_err_BUG("simplify, NULL input");
 4150 
 4151   switch(typ(x))
 4152   {
 4153     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
 4154     case t_PADIC: case t_QFR: case t_QFI: case t_LIST: case t_STR: case t_VECSMALL:
 4155     case t_CLOSURE: case t_ERROR: case t_INFINITY:
 4156       return x;
 4157     case t_COMPLEX: return isintzero(gel(x,2))? gel(x,1): x;
 4158     case t_QUAD:    return isintzero(gel(x,3))? gel(x,2): x;
 4159 
 4160     case t_POLMOD: y = cgetg(3,t_POLMOD);
 4161       z = simplify_shallow(gel(x,1));
 4162       if (typ(z) != t_POL) z = scalarpol(z, varn(gel(x,1)));
 4163       gel(y,1) = z; /* z must be a t_POL: invalid object otherwise */
 4164       gel(y,2) = simplify_shallow(gel(x,2)); return y;
 4165 
 4166     case t_POL:
 4167       lx = lg(x);
 4168       if (lx==2) return gen_0;
 4169       if (lx==3) return simplify_shallow(gel(x,2));
 4170       y = cgetg(lx,t_POL); y[1] = x[1];
 4171       for (i=2; i<lx; i++) gel(y,i) = simplify_shallow(gel(x,i));
 4172       return y;
 4173 
 4174     case t_SER:
 4175       y = cgetg_copy(x, &lx); y[1] = x[1];
 4176       for (i=2; i<lx; i++) gel(y,i) = simplify_shallow(gel(x,i));
 4177       return y;
 4178 
 4179     case t_RFRAC: y = cgetg(3,t_RFRAC);
 4180       gel(y,1) = simplify_shallow(gel(x,1));
 4181       z = simplify_shallow(gel(x,2));
 4182       if (typ(z) != t_POL) return gdiv(gel(y,1), z);
 4183       gel(y,2) = z; return y;
 4184 
 4185     case t_VEC: case t_COL: case t_MAT:
 4186       y = cgetg_copy(x, &lx);
 4187       for (i=1; i<lx; i++) gel(y,i) = simplify_shallow(gel(x,i));
 4188       return y;
 4189   }
 4190   pari_err_BUG("simplify_shallow, type unknown");
 4191   return NULL; /* LCOV_EXCL_LINE */
 4192 }
 4193 
 4194 GEN
 4195 simplify(GEN x)
 4196 {
 4197   pari_sp av = avma;
 4198   GEN y = simplify_shallow(x);
 4199   return av == avma ? gcopy(y): gerepilecopy(av, y);
 4200 }
 4201 
 4202 /*******************************************************************/
 4203 /*                                                                 */
 4204 /*                EVALUATION OF SOME SIMPLE OBJECTS                */
 4205 /*                                                                 */
 4206 /*******************************************************************/
 4207 /* q is a real symetric matrix, x a RgV. Horner-type evaluation of q(x)
 4208  * using (n^2+3n-2)/2 mul */
 4209 GEN
 4210 qfeval(GEN q, GEN x)
 4211 {
 4212   pari_sp av = avma;
 4213   long i, j, l = lg(q);
 4214   GEN z;
 4215   if (lg(x) != l) pari_err_DIM("qfeval");
 4216   if (l==1) return gen_0;
 4217   if (lgcols(q) != l) pari_err_DIM("qfeval");
 4218   /* l = lg(x) = lg(q) > 1 */
 4219   z = gmul(gcoeff(q,1,1), gsqr(gel(x,1)));
 4220   for (i=2; i<l; i++)
 4221   {
 4222     GEN c = gel(q,i), s;
 4223     if (isintzero(gel(x,i))) continue;
 4224     s = gmul(gel(c,1), gel(x,1));
 4225     for (j=2; j<i; j++) s = gadd(s, gmul(gel(c,j),gel(x,j)));
 4226     s = gadd(gshift(s,1), gmul(gel(c,i),gel(x,i)));
 4227     z = gadd(z, gmul(gel(x,i), s));
 4228   }
 4229   return gerepileupto(av,z);
 4230 }
 4231 
 4232 static GEN
 4233 qfbeval(GEN q, GEN z)
 4234 {
 4235   GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3), x = gel(z,1), y = gel(z,2);
 4236   pari_sp av = avma;
 4237   A = gadd(gmul(x, gadd(gmul(a,x), gmul(b,y))), gmul(c, gsqr(y)));
 4238   return gerepileupto(av, A);
 4239 }
 4240 static GEN
 4241 qfbevalb(GEN q, GEN z, GEN z2)
 4242 {
 4243   GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3);
 4244   GEN x = gel(z,1), y = gel(z,2);
 4245   GEN X = gel(z2,1), Y = gel(z2,2);
 4246   GEN a2 = shifti(a,1), c2 = shifti(c,1);
 4247   pari_sp av = avma;
 4248   /* a2 x X + b (x Y + X y) + c2 y Y */
 4249   A = gadd(gmul(x, gadd(gmul(a2,X), gmul(b,Y))),
 4250            gmul(y, gadd(gmul(c2,Y), gmul(b,X))));
 4251   return gerepileupto(av, gmul2n(A, -1));
 4252 }
 4253 GEN
 4254 qfb_apply_ZM(GEN q, GEN M)
 4255 {
 4256   pari_sp av = avma;
 4257   GEN A, B, C, a = gel(q,1), b = gel(q,2), c = gel(q,3);
 4258   GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
 4259   GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
 4260   GEN by = mulii(b,y), bt = mulii(b,t), bz  = mulii(b,z);
 4261   GEN a2 = shifti(a,1), c2 = shifti(c,1);
 4262 
 4263   A = addii(mulii(x, addii(mulii(a,x), by)), mulii(c, sqri(y)));
 4264   B = addii(mulii(x, addii(mulii(a2,z), bt)),
 4265             mulii(y, addii(mulii(c2,t), bz)));
 4266   C = addii(mulii(z, addii(mulii(a,z), bt)), mulii(c