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.14.0 lcov report (development 27775-aca467eab2) Lines: 2410 2605 92.5 %
Date: 2022-07-03 07:33:15 Functions: 231 241 95.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13