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

Generated by: LCOV version 1.13