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.12.1 lcov report (development 25819-e703fe1174) Lines: 2436 2633 92.5 %
Date: 2020-09-18 06:10:04 Functions: 226 236 95.8 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13