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

Generated by: LCOV version 1.13