Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - gen3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 29950-285c5b69ed) Lines: 2439 2641 92.4 %
Date: 2025-02-05 09:09:51 Functions: 234 245 95.5 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16