Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - gen1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23344-f0cc1b3f6) Lines: 1708 1810 94.4 %
Date: 2018-12-12 05:41:43 Functions: 89 89 100.0 %
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             : /**                         (first part)                           **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /* assume z[1] was created last */
      24             : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
      25             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));
      26             : 
      27             : /* assume z[1] was created last */
      28             : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
      29             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));\
      30             : else\
      31             :   gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
      32             : 
      33             : static void
      34          98 : warn_coercion(GEN x, GEN y, GEN z)
      35             : {
      36          98 :   if (DEBUGLEVEL)
      37          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      38          98 : }
      39             : 
      40             : static long
      41          21 : kro_quad(GEN x, GEN y)
      42          21 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
      43             : 
      44             : /* is -1 not a square in Zp, assume p prime */
      45             : INLINE int
      46          28 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      47             : 
      48             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      49             : static GEN addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN));
      50             : static GEN mulpp(GEN x, GEN y);
      51             : static GEN divpp(GEN x, GEN y);
      52             : /* Argument codes for inline routines
      53             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      54             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      55             :  * T: some type (to be converted to PADIC)
      56             :  */
      57             : static GEN
      58    75238890 : addRc(GEN x, GEN y) {
      59    75238890 :   GEN z = cgetg(3,t_COMPLEX);
      60    75238890 :   gel(z,1) = gadd(x,gel(y,1));
      61    75238890 :   gel(z,2) = gcopy(gel(y,2)); return z;
      62             : }
      63             : static GEN
      64    75174574 : mulRc(GEN x, GEN y) {
      65    75174574 :   GEN z = cgetg(3,t_COMPLEX);
      66    75174574 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      67    75174574 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      68             : }
      69             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      70             : static GEN
      71          49 : mulRc_direct(GEN x, GEN y) {
      72          49 :   GEN z = cgetg(3,t_COMPLEX);
      73          49 :   gel(z,1) = gmul(x,gel(y,1));
      74          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      75             : }
      76             : static GEN
      77      567871 : divRc(GEN x, GEN y) {
      78      567871 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      79      567871 :   GEN z = cgetg(3,t_COMPLEX);
      80      567871 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      81      567871 :   gel(z,2) = gmul(mt, gel(y,2));
      82      567871 :   return z;
      83             : }
      84             : static GEN
      85     6689207 : divcR(GEN x, GEN y) {
      86     6689207 :   GEN z = cgetg(3,t_COMPLEX);
      87     6689207 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      88     6689207 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      89             : }
      90             : static GEN
      91         763 : addRq(GEN x, GEN y) {
      92         763 :   GEN z = cgetg(4,t_QUAD);
      93         763 :   gel(z,1) = ZX_copy(gel(y,1));
      94         763 :   gel(z,2) = gadd(x, gel(y,2));
      95         763 :   gel(z,3) = gcopy(gel(y,3)); return z;
      96             : }
      97             : static GEN
      98        1736 : mulRq(GEN x, GEN y) {
      99        1736 :   GEN z = cgetg(4,t_QUAD);
     100        1736 :   gel(z,1) = ZX_copy(gel(y,1));
     101        1736 :   gel(z,2) = gmul(x,gel(y,2));
     102        1736 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     103             : }
     104             : static GEN
     105          28 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     106          28 :   long i = gexpo(x) - gexpo(y);
     107          28 :   if (i > 0) prec += nbits2extraprec( i );
     108          28 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     109             : }
     110             : static GEN
     111     9975336 : mulrfrac(GEN x, GEN y)
     112             : {
     113             :   pari_sp av;
     114     9975336 :   GEN z, a = gel(y,1), b = gel(y,2);
     115     9975336 :   if (is_pm1(a)) /* frequent special case */
     116             :   {
     117     2149656 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     118     2149656 :     return z;
     119             :   }
     120     7825680 :   av = avma;
     121     7825680 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     122             : }
     123             : static GEN
     124          14 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     125          14 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     126             : }
     127             : static GEN
     128          28 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     129          28 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     130             : }
     131             : static GEN
     132           7 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     133           7 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     134             : }
     135             : /* y PADIC, x + y by converting x to padic */
     136             : static GEN
     137           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     138             : 
     139           7 :   if (!valp(y)) z = cvtop2(x,y);
     140             :   else {
     141           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     142           7 :     z  = cvtop(x, gel(y,2), l);
     143             :   }
     144           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     145             : }
     146             : /* y PADIC, x * y by converting x to padic */
     147             : static GEN
     148      167055 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     149      167055 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     150             : }
     151             : /* y PADIC, non zero x / y by converting x to padic */
     152             : static GEN
     153       14119 : divTp(GEN x, GEN y) { pari_sp av = avma;
     154       14119 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     155             : }
     156             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     157             :  * converted to O(p^e) and division by 0 */
     158             : static GEN
     159        1820 : divpT(GEN x, GEN y) { pari_sp av = avma;
     160        1820 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     161             : }
     162             : 
     163             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     164             :  * clean memory from z on */
     165             : static GEN
     166     2258142 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     167     2258142 :   if (lgefint(X) == 3) {
     168     2188226 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     169     2188226 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     170             :   }
     171             :   else {
     172       69916 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     173       69877 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     174             :   }
     175     2258148 :   gel(z,1) = icopy(X); return z;
     176             : }
     177             : static GEN
     178     1109349 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     179     1109349 :   if (lgefint(X) == 3) {
     180      751299 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     181      751299 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     182             :   }
     183             :   else {
     184      358050 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     185      358037 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     186             :   }
     187     1109350 :   gel(z,1) = icopy(X); return z;
     188             : }
     189             : /* cf add_intmod_same */
     190             : static GEN
     191     2707268 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     192     2707268 :   if (lgefint(X) == 3) {
     193     2593926 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     194     2593926 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     195             :   }
     196             :   else
     197      113342 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     198     2707338 :   gel(z,1) = icopy(X); return z;
     199             : }
     200             : /* cf add_intmod_same */
     201             : static GEN
     202       73771 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     203             : {
     204       73771 :   if (lgefint(X) == 3) {
     205       51541 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     206       51534 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     207             :   }
     208             :   else
     209       22230 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     210       73764 :   gel(z,1) = icopy(X); return z;
     211             : }
     212             : 
     213             : /*******************************************************************/
     214             : /*                                                                 */
     215             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     216             : /*                                                                 */
     217             : /* (static routines are not memory clean, but OK for gerepileupto) */
     218             : /*******************************************************************/
     219             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     220             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     221             : static GEN
     222    10070353 : rfrac_denom_mul_scal(GEN d, GEN y)
     223             : {
     224    10070353 :   GEN D = RgX_Rg_mul(d, y);
     225    10070353 :   if (lg(D) != lg(d))
     226             :   { /* try to generate a meaningful diagnostic */
     227           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     228           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     229             :   }
     230    10070353 :   return D;
     231             : }
     232             : 
     233             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     234             : GEN
     235    58105479 : gred_rfrac_simple(GEN n, GEN d)
     236             : {
     237             :   GEN c, cn, cd, z;
     238    58105479 :   long dd = degpol(d);
     239             : 
     240    58105479 :   if (dd <= 0)
     241             :   {
     242       90466 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     243       90466 :     n = gdiv(n, gel(d,2));
     244       90465 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     245       90466 :     return n;
     246             :   }
     247             : 
     248    58015013 :   cd = content(d);
     249    58015013 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     250    58015013 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     251    58015013 :   if (!gequal1(cd)) {
     252     6576656 :     d = RgX_Rg_div(d,cd);
     253     6576656 :     if (!gequal1(cn))
     254             :     {
     255     1294328 :       if (gequal0(cn)) {
     256          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     257           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     258           0 :         c = gen_1;
     259             :       } else {
     260     1294279 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     261     1294279 :         c = gdiv(cn,cd);
     262             :       }
     263             :     }
     264             :     else
     265     5282328 :       c = ginv(cd);
     266             :   } else {
     267    51438357 :     if (!gequal1(cn))
     268             :     {
     269     3223613 :       if (gequal0(cn)) {
     270         910 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     271          21 :         c = gen_1;
     272             :       } else {
     273     3222703 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     274     3222703 :         c = cn;
     275             :       }
     276             :     } else {
     277    48214744 :       GEN y = cgetg(3,t_RFRAC);
     278    48214744 :       gel(y,1) = gcopy(n);
     279    48214744 :       gel(y,2) = RgX_copy(d); return y;
     280             :     }
     281             :   }
     282             : 
     283     9799331 :   if (typ(c) == t_POL)
     284             :   {
     285      905546 :     z = c;
     286      944984 :     do { z = content(z); } while (typ(z) == t_POL);
     287      905546 :     cd = denom_i(z);
     288      905546 :     cn = gmul(c, cd);
     289             :   }
     290             :   else
     291             :   {
     292     8893785 :     cn = numer_i(c);
     293     8893785 :     cd = denom_i(c);
     294             :   }
     295     9799331 :   z = cgetg(3,t_RFRAC);
     296     9799331 :   gel(z,1) = gmul(n, cn);
     297     9799331 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     298     9799331 :   return z;
     299             : }
     300             : 
     301             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     302             : static GEN
     303      249074 : fix_rfrac(GEN x, long d)
     304             : {
     305             :   GEN z, N, D;
     306      249074 :   if (!d || typ(x) == t_POL) return x;
     307      175630 :   z = cgetg(3, t_RFRAC);
     308      175630 :   N = gel(x,1);
     309      175630 :   D = gel(x,2);
     310      175630 :   if (d > 0) {
     311      344267 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     312      177975 :                                                    : monomialcopy(N,d,varn(D));
     313      166243 :     gel(z, 2) = RgX_copy(D);
     314             :   } else {
     315        9387 :     gel(z, 1) = gcopy(N);
     316        9387 :     gel(z, 2) = RgX_shift(D, -d);
     317             :   }
     318      175630 :   return z;
     319             : }
     320             : 
     321             : /* assume d != 0 */
     322             : static GEN
     323    44743745 : gred_rfrac2(GEN n, GEN d)
     324             : {
     325             :   GEN y, z;
     326             :   long v, vd, vn;
     327             : 
     328    44743745 :   n = simplify_shallow(n);
     329    44743745 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     330    37913117 :   d = simplify_shallow(d);
     331    37913117 :   if (typ(d) != t_POL) return gdiv(n,d);
     332    36733461 :   vd = varn(d);
     333    36733461 :   if (typ(n) != t_POL)
     334             :   {
     335    20606760 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     336    20605353 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     337           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     338             :   }
     339    16126701 :   vn = varn(n);
     340    16126701 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     341    15989083 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     342             : 
     343             :   /* now n and d are t_POLs in the same variable */
     344    15829331 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     345    15829331 :   if (!degpol(d))
     346             :   {
     347    12756633 :     n = RgX_Rg_div(n,gel(d,2));
     348    12756633 :     return v? RgX_mulXn(n,v): n;
     349             :   }
     350             : 
     351             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     352     3072698 :   if (!isinexact(n) && !isinexact(d))
     353             :   {
     354     3072488 :     y = RgX_divrem(n, d, &z);
     355     3072488 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     356      248864 :     z = RgX_gcd(d, z);
     357      248864 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     358             :   }
     359      249074 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     360             : }
     361             : 
     362             : /* x,y t_INT, return x/y in reduced form */
     363             : GEN
     364    71420263 : Qdivii(GEN x, GEN y)
     365             : {
     366    71420263 :   pari_sp av = avma;
     367             :   GEN r, q;
     368             : 
     369    71420263 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     370    33991122 :   if (equali1(x)) {
     371     5455145 :     long s = signe(y);
     372     5455145 :     if (!s) pari_err_INV("gdiv",y);
     373     5455040 :     if (signe(x) < 0) s = -s;
     374     5455040 :     q = cgetg(3, t_FRAC);
     375     5455040 :     gel(q,1) = s<0? gen_m1: gen_1;
     376     5455040 :     gel(q,2) = absi(y); return q;
     377             :   }
     378    28535977 :   q = dvmdii(x,y,&r);
     379    28535977 :   if (r == gen_0) return q; /* gen_0 intended */
     380    12436587 :   r = gcdii(y,r);
     381    12436587 :   if (lgefint(r) == 3)
     382             :   {
     383    12193810 :     ulong t = r[2];
     384    12193810 :     set_avma(av);
     385    12193810 :     if (t == 1) q = mkfraccopy(x,y);
     386             :     else
     387             :     {
     388     5518547 :       q = cgetg(3,t_FRAC);
     389     5518547 :       gel(q,1) = diviuexact(x,t);
     390     5518547 :       gel(q,2) = diviuexact(y,t);
     391             :     }
     392             :   }
     393             :   else
     394             :   { /* rare: r and q left on stack for efficiency */
     395      242777 :     q = cgetg(3,t_FRAC);
     396      242777 :     gel(q,1) = diviiexact(x,r);
     397      242777 :     gel(q,2) = diviiexact(y,r);
     398             :   }
     399    12436587 :   normalize_frac(q); return q;
     400             : }
     401             : 
     402             : /*******************************************************************/
     403             : /*                                                                 */
     404             : /*                          CONJUGATION                            */
     405             : /*                                                                 */
     406             : /*******************************************************************/
     407             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     408             : static GEN
     409       28007 : quad_polmod_conj(GEN x, GEN y)
     410             : {
     411             :   GEN z, u, v, a, b;
     412       28007 :   if (typ(x) != t_POL) return gcopy(x);
     413       28007 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     414       28007 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     415       28007 :   b = gel(y,3); v = gel(x,2);
     416       28007 :   z = cgetg(4, t_POL); z[1] = x[1];
     417       28007 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     418       28007 :   gel(z,3) = gneg(u); return z;
     419             : }
     420             : static GEN
     421       28007 : quad_polmod_norm(GEN x, GEN y)
     422             : {
     423             :   GEN z, u, v, a, b, c;
     424       28007 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     425       28007 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     426       28007 :   b = gel(y,3); v = gel(x,2);
     427       28007 :   c = gel(y,2);
     428       28007 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     429       28007 :   if (!gequal1(a)) z = gdiv(z, a);
     430       28007 :   return gadd(z, gsqr(v));
     431             : }
     432             : 
     433             : GEN
     434     8819001 : conj_i(GEN x)
     435             : {
     436             :   long lx, i;
     437             :   GEN y;
     438             : 
     439     8819001 :   switch(typ(x))
     440             :   {
     441             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     442     4418479 :       return x;
     443             : 
     444     4366016 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     445             : 
     446             :     case t_QUAD:
     447         595 :       y = cgetg(4,t_QUAD);
     448         595 :       gel(y,1) = gel(x,1);
     449        1190 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     450         595 :                                       : gadd(gel(x,2), gel(x,3));
     451         595 :       gel(y,3) = gneg(gel(x,3)); return y;
     452             : 
     453             :     case t_POL: case t_SER:
     454       10556 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     455       10556 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     456       10556 :       return y;
     457             : 
     458             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     459       23355 :       y = cgetg_copy(x, &lx);
     460       23355 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     461       23355 :       return y;
     462             : 
     463             :     case t_POLMOD:
     464             :     {
     465           0 :       GEN X = gel(x,1);
     466           0 :       long d = degpol(X);
     467           0 :       if (d < 2) return x;
     468           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     469             :     }
     470             :   }
     471           0 :   pari_err_TYPE("gconj",x);
     472             :   return NULL; /* LCOV_EXCL_LINE */
     473             : }
     474             : GEN
     475       30291 : gconj(GEN x)
     476       30291 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     477             : 
     478             : GEN
     479          84 : conjvec(GEN x,long prec)
     480             : {
     481             :   long lx, s, i;
     482             :   GEN z;
     483             : 
     484          84 :   switch(typ(x))
     485             :   {
     486             :     case t_INT: case t_INTMOD: case t_FRAC:
     487           0 :       return mkcolcopy(x);
     488             : 
     489             :     case t_COMPLEX: case t_QUAD:
     490           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     491             : 
     492             :     case t_FFELT:
     493          28 :       return FF_conjvec(x);
     494             : 
     495             :     case t_VEC: case t_COL:
     496           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     497           0 :       if (lx == 1) return z;
     498           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     499           0 :       s = lgcols(z);
     500           0 :       for (i=2; i<lx; i++)
     501             :       {
     502           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     503           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     504             :       }
     505           0 :       break;
     506             : 
     507             :     case t_POLMOD: {
     508          56 :       GEN T = gel(x,1), r;
     509             :       pari_sp av;
     510             : 
     511          56 :       lx = lg(T);
     512          56 :       if (lx <= 3) return cgetg(1,t_COL);
     513          56 :       x = gel(x,2);
     514         238 :       for (i=2; i<lx; i++)
     515             :       {
     516         189 :         GEN c = gel(T,i);
     517         189 :         switch(typ(c)) {
     518             :           case t_INTMOD: {
     519           7 :             GEN p = gel(c,1);
     520             :             pari_sp av;
     521           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     522           7 :             av = avma;
     523           7 :             T = RgX_to_FpX(T,p);
     524           7 :             x = RgX_to_FpX(x, p);
     525           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     526           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     527           7 :             return gerepileupto(av, z);
     528             :           }
     529             :           case t_INT:
     530         182 :           case t_FRAC: break;
     531           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     532             :         }
     533             :       }
     534          49 :       if (typ(x) != t_POL)
     535             :       {
     536           0 :         if (!is_rational_t(typ(x)))
     537           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     538           0 :         retconst_col(lx-3, gcopy(x));
     539             :       }
     540          49 :       RgX_check_QX(x,"conjvec");
     541          49 :       av = avma;
     542          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     543          49 :       r = cleanroots(T,prec);
     544          49 :       z = cgetg(lx-2,t_COL);
     545          49 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     546          49 :       return gerepileupto(av, z);
     547             :     }
     548             : 
     549             :     default:
     550           0 :       pari_err_TYPE("conjvec",x);
     551             :       return NULL; /* LCOV_EXCL_LINE */
     552             :   }
     553           0 :   return z;
     554             : }
     555             : 
     556             : 
     557             : /********************************************************************/
     558             : /**                                                                **/
     559             : /**                           ADDITION                             **/
     560             : /**                                                                **/
     561             : /********************************************************************/
     562             : /* x, y compatible PADIC, op = add or sub */
     563             : static GEN
     564     1248762 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     565             : {
     566     1248762 :   pari_sp av = avma;
     567             :   long d,e,r,rx,ry;
     568     1248762 :   GEN u, z, mod, p = gel(x,2);
     569             :   int swap;
     570             : 
     571     1248762 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     572     1248762 :   e = valp(x);
     573     1248762 :   r = valp(y); d = r-e;
     574     1248762 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     575     1248762 :   rx = precp(x);
     576     1248762 :   ry = precp(y);
     577     1248762 :   if (d) /* v(x) < v(y) */
     578             :   {
     579      602899 :     r = d+ry; z = powiu(p,d);
     580      602899 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     581      602899 :     z = mulii(z,gel(y,4));
     582      602899 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     583             :   }
     584             :   else
     585             :   {
     586             :     long c;
     587      645863 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     588      645863 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     589      645863 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     590             :     {
     591       67508 :       set_avma(av); return zeropadic(p, e+r);
     592             :     }
     593      578355 :     if (c)
     594             :     {
     595       99225 :       mod = diviiexact(mod, powiu(p,c));
     596       99225 :       r -= c;
     597       99225 :       e += c;
     598             :     }
     599             :   }
     600     1181254 :   u = modii(u, mod);
     601     1181254 :   set_avma(av); z = cgetg(5,t_PADIC);
     602     1181254 :   z[1] = evalprecp(r) | evalvalp(e);
     603     1181254 :   gel(z,2) = icopy(p);
     604     1181254 :   gel(z,3) = icopy(mod);
     605     1181254 :   gel(z,4) = icopy(u); return z;
     606             : }
     607             : /* Rg_to_Fp(t_FRAC) without GC */
     608             : static GEN
     609       56483 : Q_to_Fp(GEN x, GEN p)
     610       56483 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     611             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     612             : static GEN
     613      374701 : addQp(GEN x, GEN y)
     614             : {
     615      374701 :   pari_sp av = avma;
     616      374701 :   long d, r, e, vy = valp(y), py = precp(y);
     617      374701 :   GEN z, q, mod, u, p = gel(y,2);
     618             : 
     619      374701 :   e = Q_pvalrem(x, p, &x);
     620      374651 :   d = vy - e; r = d + py;
     621      374651 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     622      373069 :   mod = gel(y,3);
     623      373069 :   u   = gel(y,4);
     624      373069 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     625             : 
     626      373125 :   if (d > 0)
     627             :   {
     628      236827 :     q = powiu(p,d);
     629      236775 :     mod = mulii(mod, q);
     630      236781 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     631      236781 :     u = addii(x,  mulii(u, q));
     632             :   }
     633      136298 :   else if (d < 0)
     634             :   {
     635       25199 :     q = powiu(p,-d);
     636       25197 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     637       25197 :     u = addii(u, mulii(x, q));
     638       25198 :     r = py; e = vy;
     639             :   }
     640             :   else
     641             :   {
     642             :     long c;
     643      111099 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     644      111099 :     u = addii(u, x);
     645      111096 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     646             :     {
     647         966 :       set_avma(av); return zeropadic(p,e+r);
     648             :     }
     649      110131 :     if (c)
     650             :     {
     651       59865 :       mod = diviiexact(mod, powiu(p,c));
     652       59865 :       r -= c;
     653       59865 :       e += c;
     654             :     }
     655             :   }
     656      372117 :   u = modii(u, mod); set_avma(av);
     657      372097 :   z = cgetg(5,t_PADIC);
     658      372144 :   z[1] = evalprecp(r) | evalvalp(e);
     659      372150 :   gel(z,2) = icopy(p);
     660      372162 :   gel(z,3) = icopy(mod);
     661      372159 :   gel(z,4) = icopy(u); return z;
     662             : }
     663             : 
     664             : /* Mod(x,X) + Mod(y,X) */
     665             : #define addsub_polmod_same addsub_polmod_scal
     666             : /* Mod(x,X) +/- Mod(y,Y) */
     667             : static GEN
     668        1750 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     669             : {
     670        1750 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     671        1750 :   GEN z = cgetg(3,t_POLMOD);
     672        1750 :   long vx = varn(X), vy = varn(Y);
     673        1750 :   if (vx==vy) {
     674             :     pari_sp av;
     675          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     676          14 :     warn_coercion(X,Y,gel(z,1));
     677          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     678             :   }
     679        1736 :   if (varncmp(vx, vy) < 0)
     680        1736 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     681             :   else
     682           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     683        1736 :   gel(z,2) = op(x, y); return z;
     684             : }
     685             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     686             : static GEN
     687    12322416 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     688    12322416 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     689             : 
     690             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     691             : static GEN
     692      307734 : add_ser_scal(GEN y, GEN x)
     693             : {
     694             :   long i, l, ly, vy;
     695             :   GEN z;
     696             : 
     697      307734 :   if (isrationalzero(x)) return gcopy(y);
     698      281372 :   ly = lg(y);
     699      281372 :   l = valp(y);
     700      281372 :   if (l < 3-ly) return gcopy(y);
     701             :   /* l + ly >= 3 */
     702      281204 :   if (l < 0)
     703             :   {
     704         735 :     z = cgetg(ly,t_SER); z[1] = y[1];
     705         735 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     706         735 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     707         735 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     708         735 :     return z;
     709             :   }
     710      280469 :   vy = varn(y);
     711      280469 :   if (l > 0)
     712             :   {
     713       21238 :     if (ser_isexactzero(y))
     714       13650 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     715        7588 :     y -= l; ly += l;
     716        7588 :     z = cgetg(ly,t_SER);
     717        7588 :     x = gcopy(x);
     718        7588 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     719             :   }
     720             :   else
     721             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     722      259231 :     z = cgetg(ly,t_SER);
     723      259231 :     x = gadd(x, gel(y,2));
     724      259231 :     i = 3;
     725             :   }
     726      266819 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     727      266819 :   gel(z,2) = x;
     728      266819 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     729      266819 :   return gequal0(x)? normalize(z): z;
     730             : }
     731             : static long
     732     3542600 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     733             : /* x,y t_SER in the same variable: x+y */
     734             : static GEN
     735     1771692 : ser_add(GEN x, GEN y)
     736             : {
     737     1771692 :   long i, lx,ly, n = valp(y) - valp(x);
     738             :   GEN z;
     739     1771692 :   if (n < 0) { n = -n; swap(x,y); }
     740             :   /* valp(x) <= valp(y) */
     741     1771692 :   lx = _serprec(x);
     742     1771692 :   if (lx == 2) /* don't lose type information */
     743         784 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
     744     1770908 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     745     1770908 :   if (n)
     746             :   {
     747       91010 :     if (n+2 > lx) return gcopy(x);
     748       90457 :     z = cgetg(ly,t_SER);
     749       90457 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     750       90457 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     751             :   } else {
     752     1679898 :     z = cgetg(ly,t_SER);
     753     1679898 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     754             :   }
     755     1770355 :   z[1] = x[1]; return normalize(z);
     756             : }
     757             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     758             : static GEN
     759     8823518 : add_rfrac_scal(GEN y, GEN x)
     760             : {
     761             :   pari_sp av;
     762             :   GEN n;
     763             : 
     764     8823518 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     765     5149242 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     766     5149242 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     767             : }
     768             : 
     769             : /* x "scalar", ty != t_MAT and non-scalar */
     770             : static GEN
     771    19789172 : add_scal(GEN y, GEN x, long ty)
     772             : {
     773    19789172 :   switch(ty)
     774             :   {
     775    12153216 :     case t_POL: return RgX_Rg_add(y, x);
     776      307720 :     case t_SER: return add_ser_scal(y, x);
     777     4470248 :     case t_RFRAC: return add_rfrac_scal(y, x);
     778           0 :     case t_COL: return RgC_Rg_add(y, x);
     779             :     case t_VEC:
     780     2857981 :       if (isintzero(x)) return gcopy(y);
     781         168 :       break;
     782             :   }
     783         175 :   pari_err_TYPE2("+",x,y);
     784             :   return NULL; /* LCOV_EXCL_LINE */
     785             : }
     786             : 
     787             : static GEN
     788    18736552 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     789             : {
     790    18736552 :   pari_sp av0 = avma;
     791    18736552 :   GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
     792    18736552 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     793    18736552 :   int s = cmpii(x2, y2);
     794             : 
     795    18736552 :   if (!s)
     796             :   { /* common denominator: (x1 op y1) / x2 */
     797     6396511 :     n = op(x1, y1);
     798     6396511 :     if (!signe(n)) { set_avma(av0); return gen_0; }
     799     5764979 :     d = x2;
     800     5764979 :     q = dvmdii(n, d, &r);
     801     5764979 :     if (r == gen_0) { set_avma(av0); return icopy(q); }
     802     3878539 :     r = gcdii(d, r);
     803     3878539 :     if (!equali1(r)) { n = diviiexact(n, r); d = diviiexact(d, r); }
     804     3878539 :     set_avma((pari_sp)z);
     805     3878539 :     gel(z,1) = icopy(n);
     806     3878539 :     gel(z,2) = icopy(d); return z;
     807             :   }
     808    12340041 :   if (s < 0)
     809             :   {
     810     6327089 :     GEN Q = dvmdii(y2, x2, &r);
     811     6327089 :     if (r == gen_0)
     812             :     { /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     813     4321474 :       pari_sp av = avma;
     814     4321474 :       n = op(mulii(Q,x1), y1);
     815     4321474 :       q = dvmdii(n, x2, &r);
     816     4321474 :       if (r == gen_0)
     817             :       {
     818      611766 :         gel(z,1) = gerepileuptoint(av, q);
     819      611766 :         gel(z,2) = Q; return z;
     820             :       }
     821     3709708 :       r = gcdii(x2, r);
     822     3709708 :       if (!equali1(r)) { n = diviiexact(n, r); x2 = diviiexact(x2, r); }
     823     3709708 :       d = mulii(x2,Q); set_avma((pari_sp)z);
     824     3709708 :       gel(z,1) = icopy(n);
     825     3709708 :       gel(z,2) = icopy(d); return z;
     826             :     }
     827     2005615 :     delta = gcdii(x2,r);
     828             :   }
     829             :   else
     830             :   {
     831     6012952 :     GEN Q = dvmdii(x2, y2, &r);
     832     6012952 :     if (r == gen_0)
     833             :     { /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     834     4966489 :       pari_sp av = avma;
     835     4966489 :       n = op(x1, mulii(Q,y1));
     836     4966489 :       q = dvmdii(n, y2, &r);
     837     4966489 :       if (r == gen_0)
     838             :       {
     839      144384 :         gel(z,1) = gerepileuptoint(av, q);
     840      144384 :         gel(z,2) = Q; return z;
     841             :       }
     842     4822105 :       r = gcdii(y2, r);
     843     4822105 :       if (!equali1(r)) { n = diviiexact(n, r); y2 = diviiexact(y2, r); }
     844     4822105 :       d = mulii(y2,Q); set_avma((pari_sp)z);
     845     4822105 :       gel(z,1) = icopy(n);
     846     4822105 :       gel(z,2) = icopy(d); return z;
     847             :     }
     848     1046463 :     delta = gcdii(y2,r);
     849             :   }
     850             :   /* delta = gcd(x2,y2) */
     851     3052078 :   if (equali1(delta))
     852             :   { /* numerator is non-zero */
     853     1222024 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     854     1222024 :     gel(z,2) = mulii(x2,y2); return z;
     855             :   }
     856     1830054 :   x2 = diviiexact(x2,delta);
     857     1830054 :   y2 = diviiexact(y2,delta);
     858     1830054 :   n = op(mulii(x1,y2), mulii(y1,x2));
     859     1830054 :   if (!signe(n)) { set_avma(av0); return gen_0; }
     860     1830054 :   d = mulii(x2, y2);
     861     1830054 :   q = dvmdii(n, delta, &r);
     862     1830054 :   if (r == gen_0)
     863             :   {
     864      142226 :     if (equali1(d)) { set_avma(av0); return icopy(q); }
     865      142226 :     set_avma((pari_sp)z);
     866      142226 :     gel(z,2) = icopy(d);
     867      142226 :     gel(z,1) = icopy(q); return z;
     868             :   }
     869     1687828 :   r = gcdii(delta, r);
     870     1687828 :   if (!equali1(r))
     871             :   {
     872      402896 :     n     = diviiexact(n, r);
     873      402896 :     delta = diviiexact(delta, r);
     874             :   }
     875     1687828 :   d = mulii(d,delta); set_avma((pari_sp)z);
     876     1687828 :   gel(z,1) = icopy(n);
     877     1687828 :   gel(z,2) = icopy(d); return z;
     878             : }
     879             : 
     880             : /* assume x2, y2 are t_POLs in the same variable */
     881             : static GEN
     882     3037311 : add_rfrac(GEN x, GEN y)
     883             : {
     884     3037311 :   pari_sp av = avma;
     885     3037311 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     886     3037311 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     887             : 
     888     3037311 :   delta = RgX_gcd(x2,y2);
     889     3037311 :   if (!degpol(delta))
     890             :   {
     891         658 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     892         658 :     d = RgX_mul(x2, y2);
     893         658 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     894             :   }
     895     3036653 :   x2 = RgX_div(x2,delta);
     896     3036653 :   y2 = RgX_div(y2,delta);
     897     3036653 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     898     3036653 :   if (!signe(n))
     899             :   {
     900      721054 :     n = simplify_shallow(n);
     901      721054 :     if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
     902           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     903             :   }
     904     2315599 :   if (degpol(n) == 0)
     905     1149721 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     906     1165878 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     907     1165878 :   if (isexactzero(r))
     908             :   {
     909             :     GEN z;
     910      228045 :     d = RgX_mul(x2, y2);
     911             :     /* "constant" denominator ? */
     912      228045 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     913      228045 :     return gerepileupto(av, z);
     914             :   }
     915      937833 :   r = RgX_gcd(delta, r);
     916      937833 :   if (degpol(r))
     917             :   {
     918      160705 :     n = RgX_div(n, r);
     919      160705 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     920             :   }
     921             :   else
     922      777128 :     d = RgX_mul(gel(x,2), y2);
     923      937833 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     924             : }
     925             : 
     926             : GEN
     927  2130299099 : gadd(GEN x, GEN y)
     928             : {
     929  2130299099 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     930             :   pari_sp av, tetpil;
     931             :   GEN z, p1;
     932             : 
     933  2130299099 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     934             :   {
     935  1203570927 :     case t_INT: return addii(x,y);
     936   538755690 :     case t_REAL: return addrr(x,y);
     937     1525837 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     938     1525837 :       z = cgetg(3,t_INTMOD);
     939     1525800 :       if (X==Y || equalii(X,Y))
     940     1525777 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     941          14 :       gel(z,1) = gcdii(X,Y);
     942          14 :       warn_coercion(X,Y,gel(z,1));
     943          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     944          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     945             :     }
     946    13983289 :     case t_FRAC: return addsub_frac(x,y,addii);
     947    77852338 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     948    77852338 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     949    77852338 :       if (isintzero(gel(z,2)))
     950             :       {
     951      328937 :         set_avma((pari_sp)(z+3));
     952      328937 :         return gadd(gel(x,1),gel(y,1));
     953             :       }
     954    77523401 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     955    77523401 :       return z;
     956             :     case t_PADIC:
     957     1019900 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     958     1019900 :       return addsub_pp(x,y, addii);
     959         476 :     case t_QUAD: z = cgetg(4,t_QUAD);
     960         476 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     961         476 :       gel(z,1) = ZX_copy(gel(x,1));
     962         476 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     963         476 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     964             :     case t_POLMOD:
     965     2332729 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     966     2331014 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     967        1715 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     968     9202141 :     case t_FFELT: return FF_add(x,y);
     969             :     case t_POL:
     970    19280842 :       vx = varn(x);
     971    19280842 :       vy = varn(y);
     972    19280842 :       if (vx != vy) {
     973      806637 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     974       18401 :         else                     return RgX_Rg_add(y, x);
     975             :       }
     976    18474205 :       return RgX_add(x, y);
     977             :     case t_SER:
     978     1766876 :       vx = varn(x);
     979     1766876 :       vy = varn(y);
     980     1766876 :       if (vx != vy) {
     981          14 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
     982          14 :         else                     return add_ser_scal(y, x);
     983             :       }
     984     1766862 :       return ser_add(x, y);
     985             :     case t_RFRAC:
     986     4329198 :       vx = varn(gel(x,2));
     987     4329198 :       vy = varn(gel(y,2));
     988     4329198 :       if (vx != vy) {
     989     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
     990      538397 :         else                     return add_rfrac_scal(y, x);
     991             :       }
     992     3037311 :       return add_rfrac(x,y);
     993             :     case t_VEC:
     994      606757 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
     995      606757 :       return RgV_add(x,y);
     996             :     case t_COL:
     997      851043 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
     998      851043 :       return RgC_add(x,y);
     999             :     case t_MAT:
    1000      671902 :       lx = lg(x);
    1001      671902 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1002      671902 :       if (lx == 1) return cgetg(1, t_MAT);
    1003      671902 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1004      671895 :       return RgM_add(x,y);
    1005             : 
    1006           0 :     default: pari_err_TYPE2("+",x,y);
    1007             :   }
    1008             :   /* tx != ty */
    1009   254549171 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1010             : 
    1011   254549171 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1012             :   {
    1013             :     case t_INT:
    1014   198292999 :       switch(ty)
    1015             :       {
    1016   113294396 :         case t_REAL: return addir(x,y);
    1017             :         case t_INTMOD:
    1018      714779 :           z = cgetg(3, t_INTMOD);
    1019      714779 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1020    26937539 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1021    26937539 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1022    26937539 :           gel(z,2) = icopy(gel(y,2)); return z;
    1023    55351097 :         case t_COMPLEX: return addRc(x, y);
    1024             :         case t_PADIC:
    1025      336133 :           if (!signe(x)) return gcopy(y);
    1026      303944 :           return addQp(x,y);
    1027         721 :         case t_QUAD: return addRq(x, y);
    1028     1658334 :         case t_FFELT: return FF_Z_add(y,x);
    1029             :       }
    1030             : 
    1031             :     case t_REAL:
    1032    22623963 :       switch(ty)
    1033             :       {
    1034             :         case t_FRAC:
    1035    12648993 :           if (!signe(gel(y,1))) return rcopy(x);
    1036    12648993 :           if (!signe(x))
    1037             :           {
    1038        4825 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1039        4825 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1040             :           }
    1041    12644168 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1042    12644168 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1043     9974949 :         case t_COMPLEX: return addRc(x, y);
    1044          21 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1045             : 
    1046           0 :         default: pari_err_TYPE2("+",x,y);
    1047             :       }
    1048             : 
    1049             :     case t_INTMOD:
    1050       17640 :       switch(ty)
    1051             :       {
    1052       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1053       17507 :           z = cgetg(3, t_INTMOD);
    1054       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1055       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1056             :         }
    1057             :         case t_FFELT:
    1058          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1059           0 :             pari_err_OP("+",x,y);
    1060          14 :           return FF_Z_add(y,gel(x,2));
    1061          91 :         case t_COMPLEX: return addRc(x, y);
    1062           0 :         case t_PADIC: { GEN X = gel(x,1);
    1063           0 :           z = cgetg(3, t_INTMOD);
    1064           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1065             :         }
    1066          28 :         case t_QUAD: return addRq(x, y);
    1067             :       }
    1068             : 
    1069             :     case t_FRAC:
    1070     9984853 :       switch (ty)
    1071             :       {
    1072     9912746 :         case t_COMPLEX: return addRc(x, y);
    1073             :         case t_PADIC:
    1074       70756 :           if (!signe(gel(x,1))) return gcopy(y);
    1075       70756 :           return addQp(x,y);
    1076          14 :         case t_QUAD: return addRq(x, y);
    1077        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1078             :       }
    1079             : 
    1080             :     case t_FFELT:
    1081           0 :       pari_err_TYPE2("+",x,y);
    1082             : 
    1083             :     case t_COMPLEX:
    1084          21 :       switch(ty)
    1085             :       {
    1086             :         case t_PADIC:
    1087          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1088             :         case t_QUAD:
    1089           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1090           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1091             :       }
    1092             : 
    1093             :     case t_PADIC: /* ty == t_QUAD */
    1094           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1095             :   }
    1096             :   /* tx < ty, !is_const_t(y) */
    1097    23629693 :   switch(ty)
    1098             :   {
    1099             :     case t_MAT:
    1100       65548 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1101       65548 :       if (isrationalzero(x)) return gcopy(y);
    1102        5474 :       return RgM_Rg_add(y, x);
    1103             :     case t_COL:
    1104      330823 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1105      330823 :       return RgC_Rg_add(y, x);
    1106             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1107      332269 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1108             :   }
    1109    22901053 :   if (is_scalar_t(tx))  {
    1110    19811815 :     if (tx == t_POLMOD)
    1111             :     {
    1112       82705 :       vx = varn(gel(x,1));
    1113       82705 :       vy = gvar(y);
    1114       82705 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1115             :       else
    1116       48062 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1117       45234 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1118             :     }
    1119    19729110 :     return add_scal(y, x, ty);
    1120             :   }
    1121             :   /* x and y are not scalars, ty != t_MAT */
    1122     3089238 :   vx = gvar(x);
    1123     3089238 :   vy = gvar(y);
    1124     3089238 :   if (vx != vy) { /* x or y is treated as a scalar */
    1125       22598 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1126       22591 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1127       22591 :                                 : add_scal(y, x, ty);
    1128             :   }
    1129             :   /* vx = vy */
    1130     3066640 :   switch(tx)
    1131             :   {
    1132             :     case t_POL:
    1133     3066234 :       switch (ty)
    1134             :       {
    1135             :         case t_SER:
    1136        4851 :           if (lg(x) == 2) return gcopy(y);
    1137        4837 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1138        4837 :           i = lg(y) + valp(y) - i;
    1139        4837 :           if (i < 3) return gcopy(y);
    1140        4830 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1141        4830 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1142        4830 :           return y;
    1143             : 
    1144     3061383 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1145             :       }
    1146           0 :       break;
    1147             : 
    1148             :     case t_SER:
    1149         406 :       if (ty == t_RFRAC)
    1150             :       {
    1151             :         GEN n, d;
    1152             :         long vn, vd;
    1153         406 :         av = avma;
    1154         406 :         n = gel(y,1); vn = gval(n, vy);
    1155         406 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1156             : 
    1157         406 :         l = lg(x) + valp(x) - (vn - vd);
    1158         406 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1159             : 
    1160             :         /* take advantage of y = t^n ! */
    1161         406 :         if (degpol(d))
    1162          70 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1163             :         else {
    1164         336 :           y = gdiv(n, gel(d,2));
    1165         336 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1166             :         }
    1167         406 :         setvalp(y, valp(y) - vd);
    1168         406 :         return gerepileupto(av, gadd(y, x));
    1169             :       }
    1170           0 :       break;
    1171             :   }
    1172           0 :   pari_err_TYPE2("+",x,y);
    1173             :   return NULL; /* LCOV_EXCL_LINE */
    1174             : }
    1175             : 
    1176             : GEN
    1177    33625714 : gaddsg(long x, GEN y)
    1178             : {
    1179    33625714 :   long ty = typ(y);
    1180             :   GEN z;
    1181             : 
    1182    33625714 :   switch(ty)
    1183             :   {
    1184    12371638 :     case t_INT:  return addsi(x,y);
    1185     8118992 :     case t_REAL: return addsr(x,y);
    1186             :     case t_INTMOD:
    1187          14 :       z = cgetg(3, t_INTMOD);
    1188          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1189     5001316 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1190     5001316 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1191     5001316 :       gel(z,2) = icopy(gel(y,2)); return z;
    1192             :     case t_COMPLEX:
    1193     7538323 :       z = cgetg(3, t_COMPLEX);
    1194     7538323 :       gel(z,1) = gaddsg(x, gel(y,1));
    1195     7538323 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1196             : 
    1197      595431 :     default: return gadd(stoi(x), y);
    1198             :   }
    1199             : }
    1200             : 
    1201             : GEN
    1202     1050701 : gsubsg(long x, GEN y)
    1203             : {
    1204             :   GEN z, a, b;
    1205             :   pari_sp av;
    1206             : 
    1207     1050701 :   switch(typ(y))
    1208             :   {
    1209      101294 :     case t_INT:  return subsi(x,y);
    1210      560594 :     case t_REAL: return subsr(x,y);
    1211             :     case t_INTMOD:
    1212          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1213          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1214       42006 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1215       42006 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1216       42006 :       gel(z,2) = icopy(gel(y,2)); return z;
    1217             :     case t_COMPLEX:
    1218      318618 :       z = cgetg(3, t_COMPLEX);
    1219      318618 :       gel(z,1) = gsubsg(x, gel(y,1));
    1220      318618 :       gel(z,2) = gneg(gel(y,2)); return z;
    1221             :   }
    1222       28140 :   av = avma;
    1223       28140 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1224             : }
    1225             : 
    1226             : /********************************************************************/
    1227             : /**                                                                **/
    1228             : /**                          SUBTRACTION                           **/
    1229             : /**                                                                **/
    1230             : /********************************************************************/
    1231             : 
    1232             : GEN
    1233  1472066580 : gsub(GEN x, GEN y)
    1234             : {
    1235  1472066580 :   long tx = typ(x), ty = typ(y);
    1236             :   pari_sp av;
    1237             :   GEN z;
    1238  1472066580 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1239             :   {
    1240  1194563703 :     case t_INT: return subii(x,y);
    1241   198038743 :     case t_REAL: return subrr(x,y);
    1242     1109371 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1243     1109371 :       z = cgetg(3,t_INTMOD);
    1244     1109366 :       if (X==Y || equalii(X,Y))
    1245     1109350 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1246          14 :       gel(z,1) = gcdii(X,Y);
    1247          14 :       warn_coercion(X,Y,gel(z,1));
    1248          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1249          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1250             :     }
    1251     4753263 :     case t_FRAC: return addsub_frac(x,y, subii);
    1252    26255562 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1253    26255562 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1254    26255562 :       if (isintzero(gel(z,2)))
    1255             :       {
    1256       80919 :         set_avma((pari_sp)(z+3));
    1257       80919 :         return gsub(gel(x,1),gel(y,1));
    1258             :       }
    1259    26174643 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1260    26174643 :       return z;
    1261             :     case t_PADIC:
    1262      228855 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1263      228855 :       return addsub_pp(x,y, subii);
    1264         595 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1265         595 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1266         595 :       gel(z,1) = ZX_copy(gel(x,1));
    1267         595 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1268         595 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1269             :     case t_POLMOD:
    1270     9613934 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1271     9613899 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1272          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1273      283085 :     case t_FFELT: return FF_sub(x,y);
    1274             :     case t_POL: {
    1275    16364770 :       long vx = varn(x);
    1276    16364770 :       long vy = varn(y);
    1277    16364770 :       if (vx != vy) {
    1278       12496 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1279        2072 :         else                     return Rg_RgX_sub(x, y);
    1280             :       }
    1281    16352274 :       return RgX_sub(x, y);
    1282             :     }
    1283             :     case t_VEC:
    1284      271041 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1285      271041 :       return RgV_sub(x,y);
    1286             :     case t_COL:
    1287     3266585 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1288     3266585 :       return RgC_sub(x,y);
    1289             :     case t_MAT: {
    1290       12881 :       long lx = lg(x);
    1291       12881 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1292       12881 :       if (lx == 1) return cgetg(1, t_MAT);
    1293       12881 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1294       12881 :       return RgM_sub(x,y);
    1295             :     }
    1296     2095983 :     case t_RFRAC: case t_SER: break;
    1297             : 
    1298           0 :     default: pari_err_TYPE2("+",x,y);
    1299             :   }
    1300    17304192 :   av = avma;
    1301    17304192 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1302             : }
    1303             : 
    1304             : /********************************************************************/
    1305             : /**                                                                **/
    1306             : /**                        MULTIPLICATION                          **/
    1307             : /**                                                                **/
    1308             : /********************************************************************/
    1309             : static GEN
    1310      516194 : mul_ser_scal(GEN y, GEN x) {
    1311             :   long ly, i;
    1312             :   GEN z;
    1313      516194 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1314      512946 :   if (ser_isexactzero(y))
    1315             :   {
    1316         413 :     if (lg(y) == 2) return gcopy(y);
    1317           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1318             :   }
    1319      512533 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1320      512533 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1321      512533 :   return normalize(z);
    1322             : }
    1323             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1324             :  * [n/d a valid RFRAC]  */
    1325             : static GEN
    1326    10444535 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1327             : {
    1328    10444535 :   pari_sp av = avma;
    1329             :   GEN z;
    1330             : 
    1331    10444535 :   switch(typ(x))
    1332             :   {
    1333             :     case t_PADIC:
    1334          14 :       n = gmul(n, x);
    1335          14 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1336          14 :       return gerepileupto(av, gdiv(n,d));
    1337             : 
    1338             :     case t_INTMOD: case t_POLMOD:
    1339         434 :       n = gmul(n, x);
    1340         434 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1341         434 :       return gerepileupto(av, gdiv(n,d));
    1342             :   }
    1343    10444087 :   z = gred_rfrac2(x, d);
    1344    10444087 :   n = simplify_shallow(n);
    1345    10444087 :   if (typ(z) == t_RFRAC)
    1346             :   {
    1347     7928181 :     n = gmul(gel(z,1), n);
    1348     7928181 :     d = gel(z,2);
    1349     7928181 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1350           0 :       z = RgX_Rg_div(n, d);
    1351             :     else
    1352     7928181 :       z = gred_rfrac_simple(n, d);
    1353             :   }
    1354             :   else
    1355     2515906 :     z = gmul(z, n);
    1356    10444087 :   return gerepileupto(av, z);
    1357             : }
    1358             : static GEN
    1359    67383051 : mul_scal(GEN y, GEN x, long ty)
    1360             : {
    1361    67383051 :   switch(ty)
    1362             :   {
    1363             :     case t_POL:
    1364    58421620 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1365    54174519 :       return RgX_Rg_mul(y, x);
    1366      381468 :     case t_SER: return mul_ser_scal(y, x);
    1367     8579921 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1368             :     case t_QFI: case t_QFR:
    1369          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1370             :   }
    1371          28 :   pari_err_TYPE2("*",x,y);
    1372             :   return NULL; /* LCOV_EXCL_LINE */
    1373             : }
    1374             : 
    1375             : static GEN
    1376      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1377             : {
    1378      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1379      160830 :   long vx = gvar(X), vy = varn(y2);
    1380      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1381        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1382             : }
    1383             : /* (x1/x2) * (y1/y2) */
    1384             : static GEN
    1385     7904132 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1386             : {
    1387             :   GEN z, X, Y;
    1388     7904132 :   pari_sp av = avma;
    1389             : 
    1390     7904132 :   X = gred_rfrac2(x1, y2);
    1391     7904132 :   Y = gred_rfrac2(y1, x2);
    1392     7904132 :   if (typ(X) == t_RFRAC)
    1393             :   {
    1394     6625204 :     if (typ(Y) == t_RFRAC) {
    1395     6559628 :       x1 = gel(X,1);
    1396     6559628 :       x2 = gel(X,2);
    1397     6559628 :       y1 = gel(Y,1);
    1398     6559628 :       y2 = gel(Y,2);
    1399     6559628 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1400             :     } else
    1401       65576 :       z = mul_gen_rfrac(Y, X);
    1402             :   }
    1403     1278928 :   else if (typ(Y) == t_RFRAC)
    1404       95254 :     z = mul_gen_rfrac(X, Y);
    1405             :   else
    1406     1183674 :     z = gmul(X, Y);
    1407     7904132 :   return gerepileupto(av, z);
    1408             : }
    1409             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1410             : static GEN
    1411      269855 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1412             : {
    1413      269855 :   pari_sp av = avma;
    1414      269855 :   GEN X = gred_rfrac2(x1, y2);
    1415      269855 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1416             :   {
    1417      263303 :     x2 = RgX_mul(gel(X,2), x2);
    1418      263303 :     x1 = gel(X,1);
    1419             :   }
    1420             :   else
    1421        6552 :     x1 = X;
    1422      269855 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1423             : }
    1424             : 
    1425             : /* Mod(y, Y) * x,  assuming x scalar */
    1426             : static GEN
    1427     4490838 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1428     4490838 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1429             : 
    1430             : /* cf mulqq */
    1431             : static GEN
    1432     2679719 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1433             : {
    1434     2679719 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1435     2679719 :   pari_sp tetpil, av = avma;
    1436     2679719 :   T[1] = x[1];
    1437     2679719 :   p2 = gmul(gel(x,2), gel(y,2));
    1438     2679719 :   p3 = gmul(gel(x,3), gel(y,3));
    1439     2679719 :   p1 = gmul(gneg_i(c),p3);
    1440             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1441     2679719 :   if (typ(b) == t_INT)
    1442             :   {
    1443     2679698 :     if (signe(b))
    1444             :     {
    1445     2068402 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1446     2068402 :       if (is_pm1(b))
    1447             :       {
    1448     2067744 :         if (signe(b) > 0) p3 = gneg(p3);
    1449             :       }
    1450             :       else
    1451         658 :         p3 = gmul(negi(b), p3);
    1452             :     }
    1453             :     else
    1454             :     {
    1455      611296 :       p3 = gmul(gel(x,2),gel(y,3));
    1456      611296 :       p4 = gmul(gel(x,3),gel(y,2));
    1457             :     }
    1458             :   }
    1459             :   else
    1460             :   {
    1461          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1462          21 :     p3 = gmul(gneg_i(b), p3);
    1463             :   }
    1464     2679719 :   tetpil = avma;
    1465     2679719 :   gel(T,2) = gadd(p2, p1);
    1466     2679719 :   gel(T,3) = gadd(p4, p3);
    1467     2679719 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1468     2679719 :   return normalizepol_lg(T,4);
    1469             : }
    1470             : /* Mod(x,T) * Mod(y,T) */
    1471             : static GEN
    1472     7031181 : mul_polmod_same(GEN T, GEN x, GEN y)
    1473             : {
    1474     7031181 :   GEN z = cgetg(3,t_POLMOD), a;
    1475     7031181 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1476     7031181 :   gel(z,1) = RgX_copy(T);
    1477             :   /* x * y mod T optimised */
    1478     7031181 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1479     5086956 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1480     3942581 :     a = gmul(x, y);
    1481             :   else
    1482             :   {
    1483     3088600 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1484     2673622 :       a = quad_polmod_mul(T, x, y);
    1485             :     else
    1486      414978 :       a = RgXQ_mul(x, y, gel(z,1));
    1487             :   }
    1488     7031181 :   gel(z,2) = a; return z;
    1489             : }
    1490             : static GEN
    1491       43246 : sqr_polmod(GEN T, GEN x)
    1492             : {
    1493       43246 :   GEN a, z = cgetg(3,t_POLMOD);
    1494       43246 :   gel(z,1) = RgX_copy(T);
    1495       43246 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1496       14098 :     a = gsqr(x);
    1497             :   else
    1498             :   {
    1499       29148 :     pari_sp av = avma;
    1500       29148 :     a = RgXQ_sqr(x, gel(z,1));
    1501       29148 :     a = gerepileupto(av, a);
    1502             :   }
    1503       43246 :   gel(z,2) = a; return z;
    1504             : }
    1505             : /* Mod(x,X) * Mod(y,Y) */
    1506             : static GEN
    1507     4301353 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1508             : {
    1509     4301353 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1510     4301353 :   long vx = varn(X), vy = varn(Y);
    1511     4301353 :   GEN z = cgetg(3,t_POLMOD);
    1512             : 
    1513     4301353 :   if (vx==vy) {
    1514             :     pari_sp av;
    1515          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1516          14 :     warn_coercion(X,Y,gel(z,1));
    1517          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1518          14 :     return z;
    1519             :   }
    1520     4301339 :   if (varncmp(vx, vy) < 0)
    1521     4251527 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1522             :   else
    1523       49812 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1524     4301339 :   gel(z,2) = gmul(x, y); return z;
    1525             : }
    1526             : 
    1527             : #if 0 /* used by 3M only */
    1528             : /* set z = x+y and return 1 if x,y have the same sign
    1529             :  * set z = x-y and return 0 otherwise */
    1530             : static int
    1531             : did_add(GEN x, GEN y, GEN *z)
    1532             : {
    1533             :   long tx = typ(x), ty = typ(y);
    1534             :   if (tx == ty) switch(tx)
    1535             :   {
    1536             :     case t_INT: *z = addii(x,y); return 1;
    1537             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1538             :     case t_REAL:
    1539             :       if (signe(x) == -signe(y))
    1540             :       { *z = subrr(x,y); return 0; }
    1541             :       else
    1542             :       { *z = addrr(x,y); return 1; }
    1543             :   }
    1544             :   if (tx == t_REAL) switch(ty)
    1545             :   {
    1546             :     case t_INT:
    1547             :       if (signe(x) == -signe(y))
    1548             :       { *z = subri(x,y); return 0; }
    1549             :       else
    1550             :       { *z = addri(x,y); return 1; }
    1551             :     case t_FRAC:
    1552             :       if (signe(x) == -signe(gel(y,1)))
    1553             :       { *z = gsub(x,y); return 0; }
    1554             :       else
    1555             :       { *z = gadd(x,y); return 1; }
    1556             :   }
    1557             :   else if (ty == t_REAL) switch(tx)
    1558             :   {
    1559             :     case t_INT:
    1560             :       if (signe(x) == -signe(y))
    1561             :       { *z = subir(x,y); return 0; }
    1562             :       else
    1563             :       { *z = addir(x,y); return 1; }
    1564             :     case t_FRAC:
    1565             :       if (signe(gel(x,1)) == -signe(y))
    1566             :       { *z = gsub(x,y); return 0; }
    1567             :       else
    1568             :       { *z = gadd(x,y); return 1; }
    1569             :   }
    1570             :   *z = gadd(x,y); return 1;
    1571             : }
    1572             : #endif
    1573             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1574             : static GEN
    1575     4286258 : mulcIR(GEN x, GEN y)
    1576             : {
    1577     4286258 :   GEN z = cgetg(3,t_COMPLEX);
    1578     4286258 :   pari_sp av = avma;
    1579     4286258 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1580     4286258 :   gel(z,2) = gmul(y, gel(x,1));
    1581     4286258 :   return z;
    1582             : 
    1583             : }
    1584             : /* x,y COMPLEX */
    1585             : static GEN
    1586   109103968 : mulcc(GEN x, GEN y)
    1587             : {
    1588   109103968 :   GEN xr = gel(x,1), xi = gel(x,2);
    1589   109103968 :   GEN yr = gel(y,1), yi = gel(y,2);
    1590             :   GEN p1, p2, p3, p4, z;
    1591             :   pari_sp tetpil, av;
    1592             : 
    1593   109103968 :   if (isintzero(xr))
    1594             :   {
    1595     4167114 :     if (isintzero(yr)) {
    1596      610903 :       av = avma;
    1597      610903 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1598             :     }
    1599     3556211 :     return mulcIR(y, xi);
    1600             :   }
    1601   104936854 :   if (isintzero(yr)) return mulcIR(x, yi);
    1602             : 
    1603   104206807 :   z = cgetg(3,t_COMPLEX); av = avma;
    1604             : #if 0
    1605             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1606             :    * e.g. xr + xi if exponents differ */
    1607             :   if (did_add(xr, xi, &p3))
    1608             :   {
    1609             :     if (did_add(yr, yi, &p4)) {
    1610             :     /* R = xr*yr - xi*yi
    1611             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1612             :       p1 = gmul(xr,yr);
    1613             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1614             :       p3 = gmul(p3, p4);
    1615             :       p4 = gsub(p2, p1);
    1616             :     } else {
    1617             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1618             :      * I = xr*yi + xi*yr */
    1619             :       p1 = gmul(p3,p4);
    1620             :       p3 = gmul(xr,yi);
    1621             :       p4 = gmul(xi,yr);
    1622             :       p2 = gsub(p3, p4);
    1623             :     }
    1624             :   } else {
    1625             :     if (did_add(yr, yi, &p4)) {
    1626             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1627             :       * I = xr*yi +xi*yr */
    1628             :       p1 = gmul(p3,p4);
    1629             :       p3 = gmul(xr,yi);
    1630             :       p4 = gmul(xi,yr);
    1631             :       p2 = gsub(p4, p3);
    1632             :     } else {
    1633             :     /* R = xr*yr - xi*yi
    1634             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1635             :       p3 = gneg( gmul(p3, p4) );
    1636             :       p1 = gmul(xr,yr);
    1637             :       p2 = gmul(xi,yi);
    1638             :       p4 = gadd(p1, p2);
    1639             : 
    1640             :       p2 = gneg(p2);
    1641             :     }
    1642             :   }
    1643             :   tetpil = avma;
    1644             :   gel(z,1) = gadd(p1,p2);
    1645             :   gel(z,2) = gadd(p3,p4);
    1646             : #else
    1647   104206807 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1648             :   { /* 3M formula */
    1649      615927 :     p3 = addii(xr,xi);
    1650      615927 :     p4 = addii(yr,yi);
    1651      615927 :     p1 = mulii(xr,yr);
    1652      615927 :     p2 = mulii(xi,yi);
    1653      615927 :     p3 = mulii(p3,p4);
    1654      615927 :     p4 = addii(p2,p1);
    1655      615927 :     tetpil = avma;
    1656      615927 :     gel(z,1) = subii(p1,p2);
    1657      615927 :     gel(z,2) = subii(p3,p4);
    1658     1118888 :     if (!signe(gel(z,2)))
    1659      112966 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1660             :   }
    1661             :   else
    1662             :   { /* naive 4M formula: avoid all loss of accuracy */
    1663   103590880 :     p1 = gmul(xr,yr);
    1664   103590880 :     p2 = gmul(xi,yi);
    1665   103590880 :     p3 = gmul(xr,yi);
    1666   103590880 :     p4 = gmul(xi,yr);
    1667   103590880 :     tetpil = avma;
    1668   103590880 :     gel(z,1) = gsub(p1,p2);
    1669   103590880 :     gel(z,2) = gadd(p3,p4);
    1670   103590880 :     if (isintzero(gel(z,2)))
    1671             :     {
    1672       14312 :       cgiv(gel(z,2));
    1673       14312 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1674             :     }
    1675             :   }
    1676             : #endif
    1677             : 
    1678   104079529 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1679             : }
    1680             : /* x,y PADIC */
    1681             : static GEN
    1682     1857631 : mulpp(GEN x, GEN y) {
    1683     1857631 :   long l = valp(x) + valp(y);
    1684             :   pari_sp av;
    1685             :   GEN z, t;
    1686     1857631 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1687     1857642 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1688     1659009 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1689             : 
    1690     1579703 :   t = (precp(x) > precp(y))? y: x;
    1691     1579703 :   z = cgetp(t); setvalp(z,l); av = avma;
    1692     1579694 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1693     1579656 :   set_avma(av); return z;
    1694             : }
    1695             : /* x,y QUAD */
    1696             : static GEN
    1697        1316 : mulqq(GEN x, GEN y) {
    1698        1316 :   GEN z = cgetg(4,t_QUAD);
    1699        1316 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1700             :   pari_sp av, tetpil;
    1701        1316 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1702             : 
    1703        1316 :   gel(z,1) = ZX_copy(P); av = avma;
    1704        1316 :   p2 = gmul(gel(x,2),gel(y,2));
    1705        1316 :   p3 = gmul(gel(x,3),gel(y,3));
    1706        1316 :   p1 = gmul(gneg_i(c),p3);
    1707             : 
    1708        1316 :   if (signe(b))
    1709        1204 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1710             :   else
    1711             :   {
    1712         112 :     p3 = gmul(gel(x,2),gel(y,3));
    1713         112 :     p4 = gmul(gel(x,3),gel(y,2));
    1714             :   }
    1715        1316 :   tetpil = avma;
    1716        1316 :   gel(z,2) = gadd(p2,p1);
    1717        1316 :   gel(z,3) = gadd(p4,p3);
    1718        1316 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1719             : }
    1720             : 
    1721             : GEN
    1722     2328328 : mulcxI(GEN x)
    1723             : {
    1724             :   GEN z;
    1725     2328328 :   switch(typ(x))
    1726             :   {
    1727             :     case t_INT: case t_REAL: case t_FRAC:
    1728      287777 :       return mkcomplex(gen_0, x);
    1729             :     case t_COMPLEX:
    1730     2040495 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1731     2038178 :       z = cgetg(3,t_COMPLEX);
    1732     2038178 :       gel(z,1) = gneg(gel(x,2));
    1733     2038178 :       gel(z,2) = gel(x,1); return z;
    1734             :     default:
    1735          56 :       return gmul(gen_I(), x);
    1736             :   }
    1737             : }
    1738             : GEN
    1739       49161 : mulcxmI(GEN x)
    1740             : {
    1741             :   GEN z;
    1742       49161 :   switch(typ(x))
    1743             :   {
    1744             :     case t_INT: case t_REAL: case t_FRAC:
    1745        1834 :       return mkcomplex(gen_0, gneg(x));
    1746             :     case t_COMPLEX:
    1747       47264 :       if (isintzero(gel(x,1))) return gel(x,2);
    1748       46165 :       z = cgetg(3,t_COMPLEX);
    1749       46165 :       gel(z,1) = gel(x,2);
    1750       46165 :       gel(z,2) = gneg(gel(x,1)); return z;
    1751             :     default:
    1752          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1753             :   }
    1754             : }
    1755             : /* x * I^k */
    1756             : GEN
    1757        1477 : mulcxpowIs(GEN x, long k)
    1758             : {
    1759        1477 :   switch (k & 3)
    1760             :   {
    1761         280 :     case 1: return mulcxI(x);
    1762         483 :     case 2: return gneg(x);
    1763         364 :     case 3: return mulcxmI(x);
    1764             :   }
    1765         350 :   return x;
    1766             : }
    1767             : 
    1768             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1769             : static GEN
    1770     2233737 : fill_ser(GEN z, GEN y)
    1771             : {
    1772     2233737 :   long i, lx = lg(z), ly = lg(y);
    1773     2233737 :   if (ly >= lx) {
    1774     2216461 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1775             :   } else {
    1776       17276 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1777       17276 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1778             :   }
    1779     2233737 :   return normalize(z);
    1780             : }
    1781             : 
    1782             : GEN
    1783  3588840487 : gmul(GEN x, GEN y)
    1784             : {
    1785             :   long tx, ty, lx, ly, vx, vy, i, l;
    1786             :   pari_sp av, tetpil;
    1787             :   GEN z, p1, p2;
    1788             : 
    1789  3588840487 :   if (x == y) return gsqr(x);
    1790  2654487447 :   tx = typ(x); ty = typ(y);
    1791  2654487447 :   if (tx == ty) switch(tx)
    1792             :   {
    1793  1237826196 :     case t_INT: return mulii(x,y);
    1794   724903927 :     case t_REAL: return mulrr(x,y);
    1795     1733936 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1796     1733936 :       z = cgetg(3,t_INTMOD);
    1797     1733924 :       if (X==Y || equalii(X,Y))
    1798     1733889 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1799          14 :       gel(z,1) = gcdii(X,Y);
    1800          14 :       warn_coercion(X,Y,gel(z,1));
    1801          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1802          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1803             :     }
    1804             :     case t_FRAC:
    1805             :     {
    1806     9659156 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1807     9659156 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1808     9659156 :       z=cgetg(3,t_FRAC);
    1809     9659156 :       p1 = gcdii(x1, y2);
    1810     9659156 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1811     9659156 :       p1 = gcdii(x2, y1);
    1812     9659156 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1813     9659156 :       tetpil = avma;
    1814     9659156 :       gel(z,2) = mulii(x2,y2);
    1815     9659156 :       gel(z,1) = mulii(x1,y1);
    1816     9659156 :       fix_frac_if_int_GC(z,tetpil); return z;
    1817             :     }
    1818   106518915 :     case t_COMPLEX: return mulcc(x, y);
    1819     1681373 :     case t_PADIC: return mulpp(x, y);
    1820         994 :     case t_QUAD: return mulqq(x, y);
    1821     9584562 :     case t_FFELT: return FF_mul(x, y);
    1822             :     case t_POLMOD:
    1823    10902202 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1824     6600849 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1825     4301353 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1826             :     case t_POL:
    1827    43353721 :       vx = varn(x);
    1828    43353721 :       vy = varn(y);
    1829    43353721 :       if (vx != vy) {
    1830     4780302 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1831     2063855 :         else                     return RgX_Rg_mul(y, x);
    1832             :       }
    1833    38573419 :       return RgX_mul(x, y);
    1834             : 
    1835             :     case t_SER: {
    1836     2237251 :       vx = varn(x);
    1837     2237251 :       vy = varn(y);
    1838     2237251 :       if (vx != vy) {
    1839        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1840        3675 :         else                     return mul_ser_scal(y, x);
    1841             :       }
    1842     2233576 :       lx = lg(x);
    1843     2233576 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1844     2233576 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1845     2233548 :       av = avma; z = cgetg(lx,t_SER);
    1846     2233548 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1847     2233548 :       x = ser2pol_i(x, lx);
    1848     2233548 :       y = ser2pol_i(y, lx);
    1849     2233548 :       y = RgXn_mul(x, y, lx-2);
    1850     2233548 :       return gerepilecopy(av, fill_ser(z,y));
    1851             :     }
    1852    16035591 :     case t_QFI: return qficomp(x,y);
    1853          14 :     case t_QFR: return qfrcomp(x,y);
    1854     6717623 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1855     1778358 :     case t_MAT: return RgM_mul(x, y);
    1856             : 
    1857             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1858        1533 :       z = cgetg_copy(x, &l);
    1859        1533 :       if (l != lg(y)) break;
    1860       16933 :       for (i=1; i<l; i++)
    1861             :       {
    1862       15400 :         long yi = y[i];
    1863       15400 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1864       15400 :         z[i] = x[yi];
    1865             :       }
    1866        1533 :       return z;
    1867             : 
    1868             : 
    1869             :     default:
    1870           0 :       pari_err_TYPE2("*",x,y);
    1871             :   }
    1872             :   /* tx != ty */
    1873   481552095 :   if (is_const_t(ty) && is_const_t(tx))  {
    1874   401602226 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1875   401602226 :     switch(tx) {
    1876             :     case t_INT:
    1877   359660618 :       switch(ty)
    1878             :       {
    1879   233282073 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1880             :         case t_INTMOD:
    1881      973175 :           z = cgetg(3, t_INTMOD);
    1882      973149 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1883             :         case t_FRAC:
    1884    45380529 :           if (!signe(x)) return gen_0;
    1885    28275841 :           z=cgetg(3,t_FRAC);
    1886    28275841 :           p1 = gcdii(x,gel(y,2));
    1887    28275841 :           if (equali1(p1))
    1888             :           {
    1889    16424392 :             set_avma((pari_sp)z);
    1890    16424392 :             gel(z,2) = icopy(gel(y,2));
    1891    16424392 :             gel(z,1) = mulii(gel(y,1), x);
    1892             :           }
    1893             :           else
    1894             :           {
    1895    11851449 :             x = diviiexact(x,p1); tetpil = avma;
    1896    11851449 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1897    11851449 :             gel(z,1) = mulii(gel(y,1), x);
    1898    11851449 :             fix_frac_if_int_GC(z,tetpil);
    1899             :           }
    1900    28275841 :           return z;
    1901    77996142 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1902      412562 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1903        1330 :         case t_QUAD: return mulRq(x,y);
    1904     1614814 :         case t_FFELT: return FF_Z_mul(y,x);
    1905             :       }
    1906             : 
    1907             :     case t_REAL:
    1908    39975883 :       switch(ty)
    1909             :       {
    1910     9975336 :         case t_FRAC: return mulrfrac(x, y);
    1911    30000498 :         case t_COMPLEX: return mulRc(x, y);
    1912           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1913          42 :         default: pari_err_TYPE2("*",x,y);
    1914             :       }
    1915             : 
    1916             :     case t_INTMOD:
    1917        7819 :       switch(ty)
    1918             :       {
    1919        7364 :         case t_FRAC: { GEN X = gel(x,1);
    1920        7364 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1921        7364 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1922             :         }
    1923          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1924         301 :         case t_PADIC: { GEN X = gel(x,1);
    1925         301 :           z = cgetg(3, t_INTMOD);
    1926         301 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1927             :         }
    1928          42 :         case t_QUAD: return mulRq(x, y);
    1929             :         case t_FFELT:
    1930          63 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1931           0 :             pari_err_OP("*",x,y);
    1932          63 :           return FF_Z_mul(y,gel(x,2));
    1933             :       }
    1934             : 
    1935             :     case t_FRAC:
    1936     1957840 :       switch(ty)
    1937             :       {
    1938     1924660 :         case t_COMPLEX: return mulRc(x, y);
    1939       30877 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1940         308 :         case t_QUAD:  return mulRq(x, y);
    1941        1995 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1942             :       }
    1943             : 
    1944             :     case t_FFELT:
    1945          35 :       pari_err_TYPE2("*",x,y);
    1946             : 
    1947             :     case t_COMPLEX:
    1948          21 :       switch(ty)
    1949             :       {
    1950             :         case t_PADIC:
    1951          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1952             :         case t_QUAD:
    1953           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1954           7 :           return mulqf(y, x, lx);
    1955             :       }
    1956             : 
    1957             :     case t_PADIC: /* ty == t_QUAD */
    1958          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1959             :     }
    1960             :   }
    1961             : 
    1962    79950040 :   if (is_matvec_t(ty))
    1963             :   {
    1964     5429725 :     if (!is_matvec_t(tx))
    1965             :     {
    1966     5278155 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1967     5278155 :       z = cgetg_copy(y, &ly);
    1968     5278155 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1969     5278141 :       return z;
    1970             :     }
    1971      151570 :     switch(tx)
    1972             :     {
    1973             :       case t_VEC:
    1974       41290 :         switch(ty) {
    1975       15015 :           case t_COL: return RgV_RgC_mul(x,y);
    1976       26275 :           case t_MAT: return RgV_RgM_mul(x,y);
    1977             :         }
    1978           0 :         break;
    1979             :       case t_COL:
    1980        1687 :         switch(ty) {
    1981        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    1982           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    1983             :         }
    1984           0 :         break;
    1985             :       case t_MAT:
    1986      108593 :         switch(ty) {
    1987           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    1988      108593 :           case t_COL: return RgM_RgC_mul(x,y);
    1989             :         }
    1990             :     }
    1991             :   }
    1992    74520349 :   if (is_matvec_t(tx))
    1993             :   {
    1994      373718 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    1995      373718 :     z = cgetg_copy(x, &lx);
    1996      373723 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    1997      373719 :     return z;
    1998             :   }
    1999    74146630 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2000             :   /* tx < ty, !ismatvec(x and y) */
    2001             : 
    2002    74146630 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2003     4452303 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2004    69694327 :   if (is_scalar_t(tx)) {
    2005    66276372 :     if (tx == t_POLMOD) {
    2006     4861589 :       vx = varn(gel(x,1));
    2007     4861589 :       vy = gvar(y);
    2008     4861589 :       if (vx != vy) {
    2009     4431642 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2010       38535 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2011             :       }
    2012             :       /* error if ty == t_SER */
    2013      429947 :       av = avma; y = gmod(y, gel(x,1));
    2014      429940 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2015             :     }
    2016    61414783 :     return mul_scal(y, x, ty);
    2017             :   }
    2018             : 
    2019             :   /* x and y are not scalars, nor matvec */
    2020     3417955 :   vx = gvar(x);
    2021     3417955 :   vy = gvar(y);
    2022     3417955 :   if (vx != vy) /* x or y is treated as a scalar */
    2023     2781482 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2024     2781482 :                                 : mul_scal(y, x, ty);
    2025             :   /* vx = vy */
    2026     1997884 :   switch(tx)
    2027             :   {
    2028             :     case t_POL:
    2029     1997856 :       switch (ty)
    2030             :       {
    2031             :         case t_SER:
    2032             :         {
    2033             :           long vn;
    2034      133270 :           if (lg(x) == 2) return zeropol(vx);
    2035      133270 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2036      133039 :           av = avma;
    2037      133039 :           vn = RgX_valrem(x, &x);
    2038             :           /* take advantage of x = t^n ! */
    2039      133039 :           if (degpol(x)) {
    2040        1988 :             p1 = RgX_to_ser(x,lg(y));
    2041        1988 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2042        1988 :             p2 = gmul(p1,y);
    2043        1988 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2044             :           } else {
    2045      131051 :             set_avma(av);
    2046      131051 :             p2 = mul_ser_scal(y, gel(x,2));
    2047             :           }
    2048      133039 :           setvalp(p2, valp(p2) + vn);
    2049      133039 :           return p2;
    2050             :         }
    2051             : 
    2052     1864586 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2053             :       }
    2054           0 :       break;
    2055             : 
    2056             :     case t_SER:
    2057          28 :       switch (ty)
    2058             :       {
    2059             :         case t_RFRAC:
    2060          28 :           av = avma;
    2061          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2062             :       }
    2063           0 :       break;
    2064             :   }
    2065           0 :   pari_err_TYPE2("*",x,y);
    2066             :   return NULL; /* LCOV_EXCL_LINE */
    2067             : }
    2068             : 
    2069             : /* return a non-normalized result */
    2070             : GEN
    2071       88763 : sqr_ser_part(GEN x, long l1, long l2)
    2072             : {
    2073             :   long i, j, l;
    2074             :   pari_sp av;
    2075             :   GEN Z, z, p1, p2;
    2076             :   long mi;
    2077       88763 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2078       88756 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2079       88756 :   Z = cgetg(l2-l1+3, t_SER);
    2080       88756 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2081       88756 :   z = Z + 2-l1;
    2082       88756 :   x += 2; mi = 0;
    2083      340921 :   for (i=0; i<l1; i++)
    2084             :   {
    2085      252165 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2086             :   }
    2087             : 
    2088      608352 :   for (i=l1; i<=l2; i++)
    2089             :   {
    2090      519596 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2091      519596 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2092     1064075 :     for (j=i-mi; j<=minss(l,mi); j++)
    2093      544479 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2094      519596 :     p1 = gshift(p1,1);
    2095      519596 :     if ((i&1) == 0 && p2[i>>1])
    2096       50500 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2097      519596 :     gel(z,i) = gerepileupto(av,p1);
    2098             :   }
    2099       88756 :   return Z;
    2100             : }
    2101             : 
    2102             : GEN
    2103  1017025521 : gsqr(GEN x)
    2104             : {
    2105             :   long i, lx;
    2106             :   pari_sp av, tetpil;
    2107             :   GEN z, p1, p2, p3, p4;
    2108             : 
    2109  1017025521 :   switch(typ(x))
    2110             :   {
    2111   944213555 :     case t_INT: return sqri(x);
    2112    26900050 :     case t_REAL: return sqrr(x);
    2113      141958 :     case t_INTMOD: { GEN X = gel(x,1);
    2114      141958 :       z = cgetg(3,t_INTMOD);
    2115      141943 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2116      141949 :       gel(z,1) = icopy(X); return z;
    2117             :     }
    2118     1283258 :     case t_FRAC: return sqrfrac(x);
    2119             : 
    2120             :     case t_COMPLEX:
    2121     3121587 :       if (isintzero(gel(x,1))) {
    2122       31509 :         av = avma;
    2123       31509 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2124             :       }
    2125     3090078 :       z = cgetg(3,t_COMPLEX); av = avma;
    2126     3090078 :       p1 = gadd(gel(x,1),gel(x,2));
    2127     3090078 :       p2 = gsub(gel(x,1), gel(x,2));
    2128     3090078 :       p3 = gmul(gel(x,1),gel(x,2));
    2129     3090078 :       tetpil = avma;
    2130     3090078 :       gel(z,1) = gmul(p1,p2);
    2131     3090078 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2132             : 
    2133             :     case t_PADIC:
    2134       28777 :       z = cgetg(5,t_PADIC);
    2135       28777 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2136       28777 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2137       28777 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2138       28777 :       gel(z,2) = icopy(gel(x,2));
    2139       28777 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2140       28777 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2141       28777 :       return z;
    2142             : 
    2143          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2144          28 :       p1 = gel(x,1);
    2145          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2146          28 :       p2 = gsqr(gel(x,2));
    2147          28 :       p3 = gsqr(gel(x,3));
    2148          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2149             : 
    2150          28 :       if (gequal0(gel(p1,3)))
    2151             :       {
    2152           7 :         tetpil = avma;
    2153           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2154           7 :         av = avma;
    2155           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2156           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2157             :       }
    2158             : 
    2159          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2160          21 :       tetpil = avma;
    2161          21 :       gel(z,2) = gadd(p2,p4);
    2162          21 :       gel(z,3) = gadd(p1,p3);
    2163          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2164             : 
    2165             :     case t_POLMOD:
    2166       43246 :       return sqr_polmod(gel(x,1), gel(x,2));
    2167             : 
    2168     2893686 :     case t_FFELT: return FF_sqr(x);
    2169             : 
    2170      108824 :     case t_POL: return RgX_sqr(x);
    2171             : 
    2172             :     case t_SER:
    2173       24143 :       lx = lg(x);
    2174       24143 :       if (ser_isexactzero(x)) {
    2175           7 :         GEN z = gcopy(x);
    2176           7 :         setvalp(z, 2*valp(x));
    2177           7 :         return z;
    2178             :       }
    2179       24136 :       if (lx < 40)
    2180       23947 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2181             :       else
    2182             :       {
    2183         189 :         pari_sp av = avma;
    2184         189 :         GEN z = cgetg(lx,t_SER);
    2185         189 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2186         189 :         x = ser2pol_i(x,lx);
    2187         189 :         x = RgXn_sqr(x, lx-2);
    2188         189 :         return gerepilecopy(av, fill_ser(z,x));
    2189             :       }
    2190             : 
    2191           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2192           7 :       gel(z,1) = gsqr(gel(x,1));
    2193           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2194             : 
    2195         980 :     case t_MAT: return RgM_sqr(x);
    2196          14 :     case t_QFR: return qfrsqr(x);
    2197    38264792 :     case t_QFI: return qfisqr(x);
    2198             :     case t_VECSMALL:
    2199         616 :       z = cgetg_copy(x, &lx);
    2200        6720 :       for (i=1; i<lx; i++)
    2201             :       {
    2202        6104 :         long xi = x[i];
    2203        6104 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2204        6104 :         z[i] = x[xi];
    2205             :       }
    2206         616 :       return z;
    2207             :   }
    2208           0 :   pari_err_TYPE2("*",x,x);
    2209             :   return NULL; /* LCOV_EXCL_LINE */
    2210             : }
    2211             : 
    2212             : /********************************************************************/
    2213             : /**                                                                **/
    2214             : /**                           DIVISION                             **/
    2215             : /**                                                                **/
    2216             : /********************************************************************/
    2217             : static GEN
    2218      271022 : div_rfrac_scal(GEN x, GEN y)
    2219             : {
    2220      271022 :   pari_sp av = avma;
    2221      271022 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2222      271022 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2223             : }
    2224             : static GEN
    2225       37445 : div_scal_rfrac(GEN x, GEN y)
    2226             : {
    2227       37445 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2228       37445 :   pari_sp av = avma;
    2229       37445 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2230             :   {
    2231           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2232           0 :     y1 = gel(y1,2);
    2233             :   }
    2234       37445 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2235             : }
    2236             : static GEN
    2237     1186509 : div_rfrac(GEN x, GEN y)
    2238     1186509 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2239             : 
    2240             : static GEN
    2241      622417 : div_ser_scal(GEN x, GEN y) {
    2242             :   long i, lx;
    2243             :   GEN z;
    2244      622417 :   if (ser_isexactzero(x))
    2245             :   {
    2246           7 :     if (lg(x) == 2) return gcopy(x);
    2247           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2248             :   }
    2249      622410 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2250      622410 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2251      622410 :   return normalize(z);
    2252             : }
    2253             : GEN
    2254        7525 : ser_normalize(GEN x)
    2255             : {
    2256        7525 :   long i, lx = lg(x);
    2257             :   GEN c, z;
    2258        7525 :   if (lx == 2) return x;
    2259        7525 :   c = gel(x,2); if (gequal1(c)) return x;
    2260        7448 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2261        7448 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2262        7448 :   return z;
    2263             : }
    2264             : 
    2265             : static GEN
    2266     1633634 : div_T_scal(GEN x, GEN y, long tx) {
    2267     1633634 :   switch(tx)
    2268             :   {
    2269      744059 :     case t_POL: return RgX_Rg_div(x, y);
    2270      622410 :     case t_SER: return div_ser_scal(x, y);
    2271      267165 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2272             :   }
    2273           0 :   pari_err_TYPE2("/",x,y);
    2274             :   return NULL; /* LCOV_EXCL_LINE */
    2275             : }
    2276             : 
    2277             : static GEN
    2278     9216975 : div_scal_pol(GEN x, GEN y) {
    2279     9216975 :   long ly = lg(y);
    2280             :   pari_sp av;
    2281     9216975 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2282     9139029 :   if (isrationalzero(x)) return zeropol(varn(y));
    2283     7091630 :   av = avma;
    2284     7091630 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2285             : }
    2286             : static GEN
    2287       32816 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2288             :   GEN z;
    2289             :   long ly, i;
    2290       32816 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2291       32816 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2292       32816 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2293       32816 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2294       32816 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2295       32816 :   y = gdiv(z,y); pari_free(z); return y;
    2296             : }
    2297             : static GEN
    2298     9240147 : div_scal_T(GEN x, GEN y, long ty) {
    2299     9240147 :   switch(ty)
    2300             :   {
    2301     9170642 :     case t_POL: return div_scal_pol(x, y);
    2302       32816 :     case t_SER: return div_scal_ser(x, y);
    2303       36689 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2304             :   }
    2305           0 :   pari_err_TYPE2("/",x,y);
    2306             :   return NULL; /* LCOV_EXCL_LINE */
    2307             : }
    2308             : 
    2309             : /* assume tx = ty = t_SER, same variable vx */
    2310             : static GEN
    2311      877445 : div_ser(GEN x, GEN y, long vx)
    2312             : {
    2313      877445 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2314             :   GEN y_lead, p1, p2, z;
    2315             : 
    2316      877445 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2317      877438 :   if (ser_isexactzero(x))
    2318             :   {
    2319       59843 :     if (lx == 2) return zeroser(vx, l);
    2320         154 :     return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), l);
    2321             :   }
    2322      817595 :   y_lead = gel(y,2);
    2323      817595 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2324             :   {
    2325          14 :     GEN y0 = y;
    2326          14 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2327          35 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2328             :     {
    2329          28 :       y_lead = gel(y,2);
    2330          28 :       if (!gequal0(y_lead)) break;
    2331             :     }
    2332          14 :     if (ly <= 2) pari_err_INV("div_ser", y0);
    2333             :   }
    2334      817588 :   if (ly < lx) lx = ly;
    2335      817588 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2336     3024841 :   for (i=3; i<lx; i++)
    2337             :   {
    2338     2207253 :     p1 = gel(y,i);
    2339     2207253 :     if (isrationalzero(p1)) p1 = NULL;
    2340     2207253 :     gel(p2,i) = p1;
    2341             :   }
    2342      817588 :   z = cgetg(lx,t_SER);
    2343      817588 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2344      817588 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2345     3024841 :   for (i=3; i<lx; i++)
    2346             :   {
    2347     2207253 :     pari_sp av = avma;
    2348     2207253 :     p1 = gel(x,i);
    2349    18591892 :     for (j=2, l=i; j<i; j++, l--)
    2350    16384639 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2351     2207253 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2352             :   }
    2353      817588 :   return normalize(z);
    2354             : }
    2355             : /* x,y compatible PADIC */
    2356             : static GEN
    2357      225693 : divpp(GEN x, GEN y) {
    2358             :   pari_sp av;
    2359             :   long a, b;
    2360             :   GEN z, M;
    2361             : 
    2362      225693 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2363      225689 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2364      225031 :   a = precp(x);
    2365      225031 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2366      225031 :   z = cgetg(5, t_PADIC);
    2367      225066 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2368      225068 :   gel(z,2) = icopy(gel(x,2));
    2369      225077 :   gel(z,3) = icopy(M); av = avma;
    2370      225078 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2371      225075 :   return z;
    2372             : }
    2373             : static GEN
    2374       43840 : div_polmod_same(GEN T, GEN x, GEN y)
    2375             : {
    2376       43840 :   long v = varn(T);
    2377       43840 :   GEN a, z = cgetg(3, t_POLMOD);
    2378       43840 :   gel(z,1) = RgX_copy(T);
    2379       43840 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2380       25102 :     a = gdiv(x, y);
    2381       18738 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2382         273 :   {
    2383         273 :     pari_sp av = avma;
    2384         273 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2385             :   }
    2386       18465 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2387        6097 :   {
    2388        6097 :     pari_sp av = avma;
    2389        6097 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2390        6097 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2391        6097 :     a = gerepileupto(av, a);
    2392             :   }
    2393             :   else
    2394             :   {
    2395       12368 :     pari_sp av = avma;
    2396       12368 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2397       12368 :     a = gerepileupto(av, a);
    2398             :   }
    2399       43840 :   gel(z,2) = a; return z;
    2400             : }
    2401             : GEN
    2402   171703833 : gdiv(GEN x, GEN y)
    2403             : {
    2404   171703833 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2405             :   pari_sp av, tetpil;
    2406             :   GEN z, p1, p2;
    2407             : 
    2408   171703833 :   if (tx == ty) switch(tx)
    2409             :   {
    2410             :     case t_INT:
    2411    66801342 :       return Qdivii(x,y);
    2412             : 
    2413    49390220 :     case t_REAL: return divrr(x,y);
    2414       38045 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2415       38045 :       z = cgetg(3,t_INTMOD);
    2416       38045 :       if (X==Y || equalii(X,Y))
    2417       38014 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2418          28 :       gel(z,1) = gcdii(X,Y);
    2419          28 :       warn_coercion(X,Y,gel(z,1));
    2420          28 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2421          28 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2422             :     }
    2423             :     case t_FRAC: {
    2424     1669768 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2425     1669768 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2426     1669768 :       z = cgetg(3, t_FRAC);
    2427     1669768 :       p1 = gcdii(x1, y1);
    2428     1669768 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2429     1669768 :       p1 = gcdii(x2, y2);
    2430     1669768 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2431     1669768 :       tetpil = avma;
    2432     1669768 :       gel(z,2) = mulii(x2,y1);
    2433     1669768 :       gel(z,1) = mulii(x1,y2);
    2434     1669768 :       normalize_frac(z);
    2435     1669768 :       fix_frac_if_int_GC(z,tetpil);
    2436     1669768 :       return z;
    2437             :     }
    2438             :     case t_COMPLEX:
    2439     2643020 :       if (isintzero(gel(y,1)))
    2440             :       {
    2441       57967 :         y = gel(y,2);
    2442       57967 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2443       40768 :         z = cgetg(3,t_COMPLEX);
    2444       40768 :         gel(z,1) = gdiv(gel(x,2), y);
    2445       40768 :         av = avma;
    2446       40768 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2447       40768 :         return z;
    2448             :       }
    2449     2585053 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2450     2585053 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2451             : 
    2452             :     case t_PADIC:
    2453      123180 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2454      123180 :       return divpp(x, y);
    2455             : 
    2456             :     case t_QUAD:
    2457         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2458         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2459         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2460             : 
    2461      236789 :     case t_FFELT: return FF_div(x,y);
    2462             : 
    2463             :     case t_POLMOD:
    2464       87415 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2465       43840 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2466             :       else {
    2467       43575 :         av = avma;
    2468       43575 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2469             :       }
    2470       87415 :       return z;
    2471             : 
    2472             :     case t_POL:
    2473    18502676 :       vx = varn(x);
    2474    18502676 :       vy = varn(y);
    2475    18502676 :       if (vx != vy) {
    2476       95046 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2477       46333 :                             else return div_scal_pol(x, y);
    2478             :       }
    2479    18407630 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2480    18407630 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2481    18221539 :       av = avma;
    2482    18221539 :       return gerepileupto(av, gred_rfrac2(x,y));
    2483             : 
    2484             :     case t_SER:
    2485      431713 :       vx = varn(x);
    2486      431713 :       vy = varn(y);
    2487      431713 :       if (vx != vy) {
    2488           7 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2489           0 :                             else return div_scal_ser(x, y);
    2490             :       }
    2491      431706 :       return div_ser(x, y, vx);
    2492             :     case t_RFRAC:
    2493     1190989 :       vx = varn(gel(x,2));
    2494     1190989 :       vy = varn(gel(y,2));
    2495     1190989 :       if (vx != vy) {
    2496        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2497         756 :                             else return div_scal_rfrac(x, y);
    2498             :       }
    2499     1186509 :       return div_rfrac(x,y);
    2500             : 
    2501           7 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2502           7 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2503             : 
    2504             :     case t_MAT:
    2505           7 :       av = avma; p1 = RgM_inv(y);
    2506           7 :       if (!p1) pari_err_INV("gdiv",y);
    2507           7 :       return gerepileupto(av, RgM_mul(x, p1));
    2508             : 
    2509           0 :     default: pari_err_TYPE2("/",x,y);
    2510             :   }
    2511             : 
    2512    30588334 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2513             :   {
    2514     2694111 :     long s = signe(x);
    2515     2694111 :     if (!s) {
    2516      681888 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2517      681888 :       switch (ty)
    2518             :       {
    2519      678696 :       default: return gen_0;
    2520             :       case t_INTMOD:
    2521          56 :         z = cgetg(3,t_INTMOD);
    2522          56 :         gel(z,1) = icopy(gel(y,1));
    2523          56 :         gel(z,2) = gen_0; return z;
    2524        3136 :       case t_FFELT: return FF_zero(y);
    2525             :       }
    2526             :     }
    2527     2012223 :     if (is_pm1(x)) {
    2528      553665 :       if (s > 0) return ginv(y);
    2529       41258 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2530             :     }
    2531     1458558 :     switch(ty)
    2532             :     {
    2533      153703 :       case t_REAL: return divir(x,y);
    2534             :       case t_INTMOD:
    2535          21 :         z = cgetg(3, t_INTMOD);
    2536          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2537             :       case t_FRAC:
    2538     1187501 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2539     1187501 :         if (equali1(p1))
    2540             :         {
    2541      810684 :           set_avma((pari_sp)z);
    2542      810684 :           gel(z,2) = icopy(gel(y,1));
    2543      810684 :           gel(z,1) = mulii(gel(y,2), x);
    2544      810684 :           normalize_frac(z);
    2545      810684 :           fix_frac_if_int(z);
    2546             :         }
    2547             :         else
    2548             :         {
    2549      376817 :           x = diviiexact(x,p1); tetpil = avma;
    2550      376817 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2551      376817 :           gel(z,1) = mulii(gel(y,2), x);
    2552      376817 :           normalize_frac(z);
    2553      376817 :           fix_frac_if_int_GC(z,tetpil);
    2554             :         }
    2555     1187501 :         return z;
    2556             : 
    2557         266 :       case t_FFELT: return Z_FF_div(x,y);
    2558      104824 :       case t_COMPLEX: return divRc(x,y);
    2559       12201 :       case t_PADIC: return divTp(x, y);
    2560             :       case t_QUAD:
    2561          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2562          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2563             :     }
    2564             :   }
    2565    27894223 :   if (gequal0(y))
    2566             :   {
    2567          42 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2568          21 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2569             :   }
    2570             : 
    2571    27894195 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2572             :   {
    2573             :     case t_REAL:
    2574     8863426 :       switch(ty)
    2575             :       {
    2576     7666705 :         case t_INT: return divri(x,y);
    2577             :         case t_FRAC:
    2578      734073 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2579      734073 :           return gerepileuptoleaf(av, z);
    2580      462620 :         case t_COMPLEX: return divRc(x, y);
    2581           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2582          21 :         default: pari_err_TYPE2("/",x,y);
    2583             :       }
    2584             : 
    2585             :     case t_INTMOD:
    2586       13482 :       switch(ty)
    2587             :       {
    2588             :         case t_INT:
    2589       13405 :           z = cgetg(3, t_INTMOD);
    2590       13405 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2591          28 :         case t_FRAC: { GEN X = gel(x,1);
    2592          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2593          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2594             :         }
    2595             :         case t_FFELT:
    2596           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2597           0 :             pari_err_OP("/",x,y);
    2598           7 :           return Z_FF_div(gel(x,2),y);
    2599             : 
    2600             :         case t_COMPLEX:
    2601           0 :           av = avma;
    2602           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2603             : 
    2604             :         case t_QUAD:
    2605           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2606           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2607             : 
    2608           7 :         case t_PADIC: { GEN X = gel(x,1);
    2609           7 :           z = cgetg(3, t_INTMOD);
    2610           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2611             :         }
    2612          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2613             :       }
    2614             : 
    2615             :     case t_FRAC:
    2616      685519 :       switch(ty)
    2617             :       {
    2618      642755 :         case t_INT: z = cgetg(3, t_FRAC);
    2619      642755 :         p1 = gcdii(y,gel(x,1));
    2620      642755 :         if (equali1(p1))
    2621             :         {
    2622      355486 :           set_avma((pari_sp)z); tetpil = 0;
    2623      355486 :           gel(z,1) = icopy(gel(x,1));
    2624             :         }
    2625             :         else
    2626             :         {
    2627      287269 :           y = diviiexact(y,p1); tetpil = avma;
    2628      287269 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2629             :         }
    2630      642755 :         gel(z,2) = mulii(gel(x,2),y);
    2631      642755 :         normalize_frac(z);
    2632      642755 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2633      642755 :         return z;
    2634             : 
    2635             :         case t_REAL:
    2636       40384 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2637       40384 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2638             : 
    2639           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2640           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2641           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2642             :         }
    2643             : 
    2644          28 :         case t_FFELT: av=avma;
    2645          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2646             : 
    2647         427 :         case t_COMPLEX: return divRc(x, y);
    2648             : 
    2649             :         case t_PADIC:
    2650        1911 :           if (!signe(gel(x,1))) return gen_0;
    2651        1911 :           return divTp(x, y);
    2652             : 
    2653             :         case t_QUAD:
    2654           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2655           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2656             :       }
    2657             : 
    2658             :     case t_FFELT:
    2659        1533 :       switch (ty)
    2660             :       {
    2661        1456 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2662          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2663             :         case t_INTMOD:
    2664           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2665           0 :             pari_err_OP("/",x,y);
    2666           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2667             :         default:
    2668          42 :         pari_err_TYPE2("/",x,y);
    2669             :       }
    2670           0 :       break;
    2671             : 
    2672             :     case t_COMPLEX:
    2673     4418291 :       switch(ty)
    2674             :       {
    2675     4418291 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2676           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2677             :         case t_PADIC:
    2678           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2679             :         case t_QUAD:
    2680           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2681           0 :           return divfq(x, y, lx);
    2682             :       }
    2683             : 
    2684             :     case t_PADIC:
    2685        1988 :       switch(ty)
    2686             :       {
    2687        1939 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2688        1939 :           return signe(gel(x,4))? divpT(x, y)
    2689        1939 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2690             :         }
    2691           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2692           7 :           z = cgetg(3, t_INTMOD);
    2693           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2694             :         }
    2695             :         case t_COMPLEX: case t_QUAD:
    2696           7 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2697           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2698             : 
    2699          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2700             :       }
    2701             : 
    2702             :     case t_QUAD:
    2703         896 :       switch (ty)
    2704             :       {
    2705             :         case t_INT: case t_INTMOD: case t_FRAC:
    2706         847 :           z = cgetg(4,t_QUAD);
    2707         847 :           gel(z,1) = ZX_copy(gel(x,1));
    2708         847 :           gel(z,2) = gdiv(gel(x,2), y);
    2709         847 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2710          28 :         case t_REAL: return divqf(x, y, lg(y));
    2711           7 :         case t_PADIC: return divTp(x, y);
    2712             :         case t_COMPLEX:
    2713           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2714           0 :           return divqf(x, y, ly);
    2715             :       }
    2716             :   }
    2717    13909081 :   switch(ty) {
    2718             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2719       31926 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2720             :     case t_MAT:
    2721          42 :       av = avma; p1 = RgM_inv(y);
    2722          35 :       if (!p1) pari_err_INV("gdiv",y);
    2723          28 :       return gerepileupto(av, gmul(x, p1));
    2724             :     case t_VEC: case t_COL:
    2725             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2726           0 :       pari_err_TYPE2("/",x,y);
    2727             :   }
    2728    13877113 :   switch(tx) {
    2729             :     case t_VEC: case t_COL: case t_MAT:
    2730     1026632 :       z = cgetg_copy(x, &lx);
    2731     1026632 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2732     1026632 :       return z;
    2733             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2734           0 :       pari_err_TYPE2("/",x,y);
    2735             :   }
    2736             : 
    2737    12850481 :   vy = gvar(y);
    2738    12850481 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2739      215264 :     vx = varn(X);
    2740      215264 :     if (vx != vy) {
    2741      214865 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2742      214536 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2743             :     }
    2744             :     /* y is POL, SER or RFRAC */
    2745         399 :     av = avma;
    2746         399 :     switch(ty)
    2747             :     {
    2748           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2749         399 :       default: y = ginvmod(gmod(y,X), X);
    2750             :     }
    2751         392 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2752             :   }
    2753             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2754             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2755             :    * variable NO_VARIABLE, then the operation is incorrect */
    2756    12635217 :   vx = gvar(x);
    2757    12635217 :   if (vx != vy) { /* includes cases where one is scalar */
    2758    10873452 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2759     9239818 :                         else return div_scal_T(x, y, ty);
    2760             :   }
    2761     1761765 :   switch(tx)
    2762             :   {
    2763             :     case t_POL:
    2764     1046108 :       switch(ty)
    2765             :       {
    2766             :         case t_SER:
    2767          70 :           if (lg(y) == 2)
    2768           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2769          70 :           p1 = RgX_to_ser(x,lg(y));
    2770          70 :           p2 = div_ser(p1, y, vx);
    2771          63 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2772          63 :           return p2;
    2773             : 
    2774             :         case t_RFRAC:
    2775             :         {
    2776     1046038 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2777     1046038 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2778           0 :             return mul_rfrac_scal(y2, y1, x);
    2779     1046038 :           av = avma;
    2780     1046038 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2781             :         }
    2782             :       }
    2783           0 :       break;
    2784             : 
    2785             :     case t_SER:
    2786      445788 :       switch(ty)
    2787             :       {
    2788             :         case t_POL:
    2789      445774 :           if (lg(x) == 2)
    2790         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2791      445669 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2792      445669 :           p2 = div_ser(x, p1, vx);
    2793      445669 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2794      445669 :           return p2;
    2795             :         case t_RFRAC:
    2796          14 :           av = avma;
    2797          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2798             :       }
    2799           0 :       break;
    2800             : 
    2801             :     case t_RFRAC:
    2802      269855 :       switch(ty)
    2803             :       {
    2804      269855 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2805             :         case t_SER:
    2806           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2807           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2808             :       }
    2809           0 :       break;
    2810             :   }
    2811          14 :   pari_err_TYPE2("/",x,y);
    2812             :   return NULL; /* LCOV_EXCL_LINE */
    2813             : }
    2814             : 
    2815             : /********************************************************************/
    2816             : /**                                                                **/
    2817             : /**                     SIMPLE MULTIPLICATION                      **/
    2818             : /**                                                                **/
    2819             : /********************************************************************/
    2820             : GEN
    2821   124377069 : gmulsg(long s, GEN y)
    2822             : {
    2823             :   long ly, i;
    2824             :   pari_sp av;
    2825             :   GEN z;
    2826             : 
    2827   124377069 :   switch(typ(y))
    2828             :   {
    2829    88948087 :     case t_INT:  return mulsi(s,y);
    2830    26556509 :     case t_REAL: return mulsr(s,y);
    2831      172368 :     case t_INTMOD: { GEN p = gel(y,1);
    2832      172368 :       z = cgetg(3,t_INTMOD);
    2833      172355 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2834      172380 :       gel(z,1) = icopy(p); return z;
    2835             :     }
    2836      627240 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2837             :     case t_FRAC:
    2838      842390 :       if (!s) return gen_0;
    2839      837987 :       z = cgetg(3,t_FRAC);
    2840      837987 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2841      837987 :       if (i == 1)
    2842             :       {
    2843      386851 :         gel(z,2) = icopy(gel(y,2));
    2844      386851 :         gel(z,1) = mulis(gel(y,1), s);
    2845             :       }
    2846             :       else
    2847             :       {
    2848      451136 :         gel(z,2) = divis(gel(y,2), i);
    2849      451136 :         gel(z,1) = mulis(gel(y,1), s/i);
    2850      451136 :         fix_frac_if_int(z);
    2851             :       }
    2852      837987 :       return z;
    2853             : 
    2854     5963683 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2855     5963683 :       gel(z,1) = gmulsg(s,gel(y,1));
    2856     5963683 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2857             : 
    2858             :     case t_PADIC:
    2859        9212 :       if (!s) return gen_0;
    2860        9212 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2861             : 
    2862           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2863           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2864           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2865           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2866             : 
    2867             :     case t_POLMOD:
    2868       83594 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2869             : 
    2870             :     case t_POL:
    2871      588711 :       if (!signe(y)) return RgX_copy(y);
    2872      573255 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2873      571897 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2874      571907 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2875      571888 :       return normalizepol_lg(z, ly);
    2876             : 
    2877             :     case t_SER:
    2878         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2879         182 :       if (!s) return Rg_get_0(y);
    2880         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2881         182 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2882         182 :       return normalize(z);
    2883             : 
    2884             :     case t_RFRAC:
    2885          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2886          14 :       if (s == 1) return gcopy(y);
    2887          14 :       if (s == -1) return gneg(y);
    2888          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2889             : 
    2890             :     case t_VEC: case t_COL: case t_MAT:
    2891      585072 :       z = cgetg_copy(y, &ly);
    2892      585072 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2893      585072 :       return z;
    2894             :   }
    2895           0 :   pari_err_TYPE("gmulsg",y);
    2896             :   return NULL; /* LCOV_EXCL_LINE */
    2897             : }
    2898             : 
    2899             : /********************************************************************/
    2900             : /**                                                                **/
    2901             : /**                       SIMPLE DIVISION                          **/
    2902             : /**                                                                **/
    2903             : /********************************************************************/
    2904             : 
    2905             : GEN
    2906    18442949 : gdivgs(GEN x, long s)
    2907             : {
    2908    18442949 :   long tx = typ(x), lx, i;
    2909             :   pari_sp av;
    2910             :   GEN z, y;
    2911             : 
    2912    18442949 :   if (!s)
    2913             :   {
    2914           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2915           0 :     pari_err_INV("gdivgs",gen_0);
    2916             :   }
    2917    18442947 :   switch(tx)
    2918             :   {
    2919             :     case t_INT:
    2920     7210427 :       av = avma; z = divis_rem(x,s,&i);
    2921     7210427 :       if (!i) return z;
    2922             : 
    2923     4876424 :       i = cgcd(s, i);
    2924     4876424 :       set_avma(av); z = cgetg(3,t_FRAC);
    2925     4876424 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    2926     4876424 :       gel(z,1) = y;
    2927     4876424 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    2928             : 
    2929             :     case t_REAL:
    2930     7740146 :       return divrs(x,s);
    2931             : 
    2932             :     case t_INTMOD:
    2933       14651 :       z = cgetg(3, t_INTMOD);
    2934       14651 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    2935             : 
    2936        1043 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    2937             : 
    2938      882739 :     case t_FRAC: z = cgetg(3, t_FRAC);
    2939      882739 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    2940      882739 :       if (i == 1)
    2941             :       {
    2942      402083 :         gel(z,2) = mulsi(s, gel(x,2));
    2943      402083 :         gel(z,1) = icopy(gel(x,1));
    2944             :       }
    2945             :       else
    2946             :       {
    2947      480656 :         gel(z,2) = mulsi(s/i, gel(x,2));
    2948      480656 :         gel(z,1) = divis(gel(x,1), i);
    2949             :       }
    2950      882739 :       normalize_frac(z);
    2951      882739 :       fix_frac_if_int(z); return z;
    2952             : 
    2953     2401020 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2954     2401020 :       gel(z,1) = gdivgs(gel(x,1),s);
    2955     2401020 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    2956             : 
    2957             :     case t_PADIC: /* divpT */
    2958             :     {
    2959       87333 :       GEN p = gel(x,2);
    2960       87333 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    2961       86612 :       av = avma;
    2962       86612 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    2963             :     }
    2964             : 
    2965           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2966           0 :       gel(z,1) = ZX_copy(gel(x,1));
    2967           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    2968           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    2969             : 
    2970             :     case t_POLMOD:
    2971       22722 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    2972             : 
    2973             :     case t_RFRAC:
    2974         140 :       if (s == 1) return gcopy(x);
    2975         133 :       else if (s == -1) return gneg(x);
    2976         133 :       return div_rfrac_scal(x, stoi(s));
    2977             : 
    2978             :     case t_POL: case t_SER:
    2979       82425 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    2980       82425 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2981       82425 :       return z;
    2982             :     case t_VEC: case t_COL: case t_MAT:
    2983         301 :       z = cgetg_copy(x, &lx);
    2984         301 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2985         301 :       return z;
    2986             : 
    2987             :   }
    2988           0 :   pari_err_TYPE2("/",x, stoi(s));
    2989             :   return NULL; /* LCOV_EXCL_LINE */
    2990             : }
    2991             : 
    2992             : /* True shift (exact multiplication by 2^n) */
    2993             : GEN
    2994    51648409 : gmul2n(GEN x, long n)
    2995             : {
    2996             :   long lx, i, k, l;
    2997             :   GEN z, a, b;
    2998             : 
    2999    51648409 :   switch(typ(x))
    3000             :   {
    3001             :     case t_INT:
    3002    37027852 :       if (n>=0) return shifti(x,n);
    3003     2394910 :       if (!signe(x)) return gen_0;
    3004     2336778 :       l = vali(x); n = -n;
    3005     2336778 :       if (n<=l) return shifti(x,-n);
    3006      304145 :       z = cgetg(3,t_FRAC);
    3007      304145 :       gel(z,1) = shifti(x,-l);
    3008      304145 :       gel(z,2) = int2n(n-l); return z;
    3009             : 
    3010             :     case t_REAL:
    3011     7136992 :       return shiftr(x,n);
    3012             : 
    3013       75345 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3014       75345 :       z = cgetg(3,t_INTMOD);
    3015       75340 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3016       75074 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3017       75080 :       gel(z,1) = icopy(b); return z;
    3018             : 
    3019      216860 :     case t_FFELT: return FF_mul2n(x,n);
    3020             : 
    3021      525790 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3022      525790 :       l = vali(a);
    3023      525790 :       k = vali(b);
    3024      525790 :       if (n+l >= k)
    3025             :       {
    3026      139608 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3027       99890 :         l = n-k; k = -k;
    3028             :       }
    3029             :       else
    3030             :       {
    3031      386182 :         k = -(l+n); l = -l;
    3032             :       }
    3033      486072 :       z = cgetg(3,t_FRAC);
    3034      486072 :       gel(z,1) = shifti(a,l);
    3035      486072 :       gel(z,2) = shifti(b,k); return z;
    3036             : 
    3037     5156555 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3038     5156555 :       gel(z,1) = gmul2n(gel(x,1),n);
    3039     5156555 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3040             : 
    3041          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3042          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3043          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3044          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3045             : 
    3046             :     case t_POLMOD:
    3047       18466 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3048             : 
    3049             :     case t_POL:
    3050      433226 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3051      433226 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3052      433226 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3053             :     case t_SER:
    3054      122062 :       if (ser_isexactzero(x)) return gcopy(x);
    3055      122048 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3056      122048 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3057      122048 :       return normalize(z); /* needed if char = 2 */
    3058             :     case t_VEC: case t_COL: case t_MAT:
    3059      917191 :       z = cgetg_copy(x, &lx);
    3060      917191 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3061      917191 :       return z;
    3062             : 
    3063             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3064          14 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3065             : 
    3066             :     case t_PADIC: /* int2n wrong if n < 0 */
    3067       18042 :       return gmul(gmul2n(gen_1,n),x);
    3068             :   }
    3069           0 :   pari_err_TYPE("gmul2n",x);
    3070             :   return NULL; /* LCOV_EXCL_LINE */
    3071             : }
    3072             : 
    3073             : /*******************************************************************/
    3074             : /*                                                                 */
    3075             : /*                              INVERSE                            */
    3076             : /*                                                                 */
    3077             : /*******************************************************************/
    3078             : static GEN
    3079      131767 : inv_polmod(GEN T, GEN x)
    3080             : {
    3081      131767 :   GEN z = cgetg(3,t_POLMOD), a;
    3082      131767 :   gel(z,1) = RgX_copy(T);
    3083      131772 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3084      104458 :     a = ginv(x);
    3085             :   else
    3086             :   {
    3087       27314 :     if (lg(T) == 5) /* quadratic fields */
    3088       21910 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3089             :     else
    3090        5404 :       a = RgXQ_inv(x, T);
    3091             :   }
    3092      131772 :   gel(z,2) = a; return z;
    3093             : }
    3094             : 
    3095             : GEN
    3096    18058549 : ginv(GEN x)
    3097             : {
    3098             :   long s;
    3099             :   pari_sp av, tetpil;
    3100             :   GEN z, y, p1, p2;
    3101             : 
    3102    18058549 :   switch(typ(x))
    3103             :   {
    3104             :     case t_INT:
    3105     1953928 :       if (is_pm1(x)) return icopy(x);
    3106     1580273 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3107     1580273 :       z = cgetg(3,t_FRAC);
    3108     1580273 :       gel(z,1) = s<0? gen_m1: gen_1;
    3109     1580273 :       gel(z,2) = absi(x); return z;
    3110             : 
    3111     3199461 :     case t_REAL: return invr(x);
    3112             : 
    3113        4641 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3114        4641 :       gel(z,1) = icopy(gel(x,1));
    3115        4641 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3116             : 
    3117             :     case t_FRAC: {
    3118      249631 :       GEN a = gel(x,1), b = gel(x,2);
    3119      249631 :       s = signe(a);
    3120      249631 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3121      111765 :       z = cgetg(3,t_FRAC);
    3122      111765 :       gel(z,1) = icopy(b);
    3123      111765 :       gel(z,2) = icopy(a);
    3124      111765 :       normalize_frac(z); return z;
    3125             :     }
    3126             :     case t_COMPLEX:
    3127     2270916 :       av=avma;
    3128     2270916 :       p1=cxnorm(x);
    3129     2270916 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3130     2270916 :       tetpil=avma;
    3131     2270916 :       return gerepile(av,tetpil,divcR(p2,p1));
    3132             : 
    3133             :     case t_QUAD:
    3134         203 :       av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
    3135         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3136             : 
    3137        5460 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3138        5460 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3139        5453 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3140        5453 :       gel(z,2) = icopy(gel(x,2));
    3141        5453 :       gel(z,3) = icopy(gel(x,3));
    3142        5453 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3143             : 
    3144      131767 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3145       14077 :     case t_FFELT: return FF_inv(x);
    3146     5421318 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3147       15330 :     case t_SER: return gdiv(gen_1,x);
    3148             : 
    3149             :     case t_RFRAC:
    3150             :     {
    3151        2912 :       GEN n = gel(x,1), d = gel(x,2);
    3152        2912 :       pari_sp av = avma, ltop;
    3153        2912 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3154             : 
    3155        2912 :       n = simplify_shallow(n);
    3156        2912 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3157             :       {
    3158        2912 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3159         679 :         ltop = avma;
    3160         679 :         z = RgX_Rg_div(d,n);
    3161             :       } else {
    3162           0 :         ltop = avma;
    3163           0 :         z = cgetg(3,t_RFRAC);
    3164           0 :         gel(z,1) = RgX_copy(d);
    3165           0 :         gel(z,2) = RgX_copy(n);
    3166             :       }
    3167         679 :       stackdummy(av, ltop);
    3168         679 :       return z;
    3169             :     }
    3170             : 
    3171             :     case t_QFR:
    3172          14 :       av = avma; z = cgetg(5, t_QFR);
    3173          14 :       gel(z,1) = gel(x,1);
    3174          14 :       gel(z,2) = negi( gel(x,2) );
    3175          14 :       gel(z,3) = gel(x,3);
    3176          14 :       gel(z,4) = negr( gel(x,4) );
    3177          14 :       return gerepileupto(av, redreal(z));
    3178             : 
    3179             :     case t_QFI:
    3180     4784971 :       y = gcopy(x);
    3181     4784971 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3182     4670039 :         togglesign(gel(y,2));
    3183     4784971 :       return y;
    3184             :     case t_MAT:
    3185        3920 :       y = RgM_inv(x);
    3186        3913 :       if (!y) pari_err_INV("ginv",x);
    3187        3843 :       return y;
    3188             :     case t_VECSMALL:
    3189             :     {
    3190           0 :       long i, lx = lg(x)-1;
    3191           0 :       y = zero_zv(lx);
    3192           0 :       for (i=1; i<=lx; i++)
    3193             :       {
    3194           0 :         long xi = x[i];
    3195           0 :         if (xi<1 || xi>lx || y[xi])
    3196           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3197           0 :         y[xi] = i;
    3198             :       }
    3199           0 :       return y;
    3200             :     }
    3201             :   }
    3202           0 :   pari_err_TYPE("inverse",x);
    3203             :   return NULL; /* LCOV_EXCL_LINE */
    3204             : }

Generated by: LCOV version 1.13