Code coverage tests

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

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

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

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

Generated by: LCOV version 1.13