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 23690-5d6e28857) Lines: 1679 1780 94.3 %
Date: 2019-03-18 05:43:21 Functions: 91 91 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 mulpp(GEN x, GEN y);
      50             : static GEN divpp(GEN x, GEN y);
      51             : /* Argument codes for inline routines
      52             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      53             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      54             :  * T: some type (to be converted to PADIC)
      55             :  */
      56             : static GEN
      57    81448255 : addRc(GEN x, GEN y) {
      58    81448255 :   GEN z = cgetg(3,t_COMPLEX);
      59    81448255 :   gel(z,1) = gadd(x,gel(y,1));
      60    81448255 :   gel(z,2) = gcopy(gel(y,2)); return z;
      61             : }
      62             : static GEN
      63    79652783 : mulRc(GEN x, GEN y) {
      64    79652783 :   GEN z = cgetg(3,t_COMPLEX);
      65    79652783 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      66    79652783 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      67             : }
      68             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      69             : static GEN
      70          49 : mulRc_direct(GEN x, GEN y) {
      71          49 :   GEN z = cgetg(3,t_COMPLEX);
      72          49 :   gel(z,1) = gmul(x,gel(y,1));
      73          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      74             : }
      75             : static GEN
      76      582668 : divRc(GEN x, GEN y) {
      77      582668 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      78      582668 :   GEN z = cgetg(3,t_COMPLEX);
      79      582668 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      80      582668 :   gel(z,2) = gmul(mt, gel(y,2));
      81      582668 :   return z;
      82             : }
      83             : static GEN
      84     6948853 : divcR(GEN x, GEN y) {
      85     6948853 :   GEN z = cgetg(3,t_COMPLEX);
      86     6948853 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      87     6948853 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      88             : }
      89             : static GEN
      90         763 : addRq(GEN x, GEN y) {
      91         763 :   GEN z = cgetg(4,t_QUAD);
      92         763 :   gel(z,1) = ZX_copy(gel(y,1));
      93         763 :   gel(z,2) = gadd(x, gel(y,2));
      94         763 :   gel(z,3) = gcopy(gel(y,3)); return z;
      95             : }
      96             : static GEN
      97        1736 : mulRq(GEN x, GEN y) {
      98        1736 :   GEN z = cgetg(4,t_QUAD);
      99        1736 :   gel(z,1) = ZX_copy(gel(y,1));
     100        1736 :   gel(z,2) = gmul(x,gel(y,2));
     101        1736 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     102             : }
     103             : static GEN
     104          28 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     105          28 :   long i = gexpo(x) - gexpo(y);
     106          28 :   if (i > 0) prec += nbits2extraprec( i );
     107          28 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     108             : }
     109             : static GEN
     110    10057470 : mulrfrac(GEN x, GEN y)
     111             : {
     112             :   pari_sp av;
     113    10057470 :   GEN z, a = gel(y,1), b = gel(y,2);
     114    10057470 :   if (is_pm1(a)) /* frequent special case */
     115             :   {
     116     2182265 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     117     2182265 :     return z;
     118             :   }
     119     7875205 :   av = avma;
     120     7875205 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     121             : }
     122             : static GEN
     123          14 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     124          14 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     125             : }
     126             : static GEN
     127          28 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     128          28 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     129             : }
     130             : static GEN
     131           7 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     132           7 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     133             : }
     134             : /* y PADIC, x + y by converting x to padic */
     135             : static GEN
     136           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     137             : 
     138           7 :   if (!valp(y)) z = cvtop2(x,y);
     139             :   else {
     140           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     141           7 :     z  = cvtop(x, gel(y,2), l);
     142             :   }
     143           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     144             : }
     145             : /* y PADIC, x * y by converting x to padic */
     146             : static GEN
     147      167063 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     148      167063 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     149             : }
     150             : /* y PADIC, non zero x / y by converting x to padic */
     151             : static GEN
     152       14119 : divTp(GEN x, GEN y) { pari_sp av = avma;
     153       14119 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     154             : }
     155             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     156             :  * converted to O(p^e) and division by 0 */
     157             : static GEN
     158        1820 : divpT(GEN x, GEN y) { pari_sp av = avma;
     159        1820 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     160             : }
     161             : 
     162             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     163             :  * clean memory from z on */
     164             : static GEN
     165     2257881 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     166     2257881 :   if (lgefint(X) == 3) {
     167     2188184 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     168     2188184 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     169             :   }
     170             :   else {
     171       69697 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     172       69695 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     173             :   }
     174     2257859 :   gel(z,1) = icopy(X); return z;
     175             : }
     176             : static GEN
     177     1109290 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     178     1109290 :   if (lgefint(X) == 3) {
     179      751299 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     180      751299 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     181             :   }
     182             :   else {
     183      357991 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     184      357988 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     185             :   }
     186     1109294 :   gel(z,1) = icopy(X); return z;
     187             : }
     188             : /* cf add_intmod_same */
     189             : static GEN
     190     2705580 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     191     2705580 :   if (lgefint(X) == 3) {
     192     2594220 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     193     2594220 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     194             :   }
     195             :   else
     196      111360 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     197     2705577 :   gel(z,1) = icopy(X); return z;
     198             : }
     199             : /* cf add_intmod_same */
     200             : static GEN
     201       73740 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     202             : {
     203       73740 :   if (lgefint(X) == 3) {
     204       51527 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     205       51520 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     206             :   }
     207             :   else
     208       22213 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     209       73743 :   gel(z,1) = icopy(X); return z;
     210             : }
     211             : 
     212             : /*******************************************************************/
     213             : /*                                                                 */
     214             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     215             : /*                                                                 */
     216             : /* (static routines are not memory clean, but OK for gerepileupto) */
     217             : /*******************************************************************/
     218             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     219             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     220             : static GEN
     221    10070738 : rfrac_denom_mul_scal(GEN d, GEN y)
     222             : {
     223    10070738 :   GEN D = RgX_Rg_mul(d, y);
     224    10070738 :   if (lg(D) != lg(d))
     225             :   { /* try to generate a meaningful diagnostic */
     226           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     227           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     228             :   }
     229    10070738 :   return D;
     230             : }
     231             : 
     232             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     233             : GEN
     234    58115557 : gred_rfrac_simple(GEN n, GEN d)
     235             : {
     236             :   GEN c, cn, cd, z;
     237    58115557 :   long dd = degpol(d);
     238             : 
     239    58115557 :   if (dd <= 0)
     240             :   {
     241       99788 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     242       99788 :     n = gdiv(n, gel(d,2));
     243       99789 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     244       99786 :     return n;
     245             :   }
     246             : 
     247    58015769 :   cd = content(d);
     248    58015769 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     249    58015769 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     250    58015769 :   if (!gequal1(cd)) {
     251     6576915 :     d = RgX_Rg_div(d,cd);
     252     6576915 :     if (!gequal1(cn))
     253             :     {
     254     1294335 :       if (gequal0(cn)) {
     255          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     256           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     257           0 :         c = gen_1;
     258             :       } else {
     259     1294286 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     260     1294286 :         c = gdiv(cn,cd);
     261             :       }
     262             :     }
     263             :     else
     264     5282580 :       c = ginv(cd);
     265             :   } else {
     266    51438854 :     if (!gequal1(cn))
     267             :     {
     268     3223725 :       if (gequal0(cn)) {
     269         910 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     270          21 :         c = gen_1;
     271             :       } else {
     272     3222815 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     273     3222815 :         c = cn;
     274             :       }
     275             :     } else {
     276    48215129 :       GEN y = cgetg(3,t_RFRAC);
     277    48215129 :       gel(y,1) = gcopy(n);
     278    48215129 :       gel(y,2) = RgX_copy(d); return y;
     279             :     }
     280             :   }
     281             : 
     282     9799702 :   if (typ(c) == t_POL)
     283             :   {
     284      905602 :     z = c;
     285      945040 :     do { z = content(z); } while (typ(z) == t_POL);
     286      905602 :     cd = denom_i(z);
     287      905602 :     cn = gmul(c, cd);
     288             :   }
     289             :   else
     290             :   {
     291     8894100 :     cn = numer_i(c);
     292     8894100 :     cd = denom_i(c);
     293             :   }
     294     9799702 :   z = cgetg(3,t_RFRAC);
     295     9799702 :   gel(z,1) = gmul(n, cn);
     296     9799702 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     297     9799702 :   return z;
     298             : }
     299             : 
     300             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     301             : static GEN
     302      249095 : fix_rfrac(GEN x, long d)
     303             : {
     304             :   GEN z, N, D;
     305      249095 :   if (!d || typ(x) == t_POL) return x;
     306      175637 :   z = cgetg(3, t_RFRAC);
     307      175637 :   N = gel(x,1);
     308      175637 :   D = gel(x,2);
     309      175637 :   if (d > 0) {
     310      344281 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     311      177982 :                                                    : monomialcopy(N,d,varn(D));
     312      166250 :     gel(z, 2) = RgX_copy(D);
     313             :   } else {
     314        9387 :     gel(z, 1) = gcopy(N);
     315        9387 :     gel(z, 2) = RgX_shift(D, -d);
     316             :   }
     317      175637 :   return z;
     318             : }
     319             : 
     320             : /* assume d != 0 */
     321             : static GEN
     322    44744431 : gred_rfrac2(GEN n, GEN d)
     323             : {
     324             :   GEN y, z;
     325             :   long v, vd, vn;
     326             : 
     327    44744431 :   n = simplify_shallow(n);
     328    44744431 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     329    37913768 :   d = simplify_shallow(d);
     330    37913768 :   if (typ(d) != t_POL) return gdiv(n,d);
     331    36734077 :   vd = varn(d);
     332    36734077 :   if (typ(n) != t_POL)
     333             :   {
     334    20606830 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     335    20605423 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     336           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     337             :   }
     338    16127247 :   vn = varn(n);
     339    16127247 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     340    15989601 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     341             : 
     342             :   /* now n and d are t_POLs in the same variable */
     343    15829849 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     344    15829849 :   if (!degpol(d))
     345             :   {
     346    12756766 :     n = RgX_Rg_div(n,gel(d,2));
     347    12756766 :     return v? RgX_mulXn(n,v): n;
     348             :   }
     349             : 
     350             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     351     3073083 :   if (!isinexact(n) && !isinexact(d))
     352             :   {
     353     3072873 :     y = RgX_divrem(n, d, &z);
     354     3072873 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     355      248885 :     z = RgX_gcd(d, z);
     356      248885 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     357             :   }
     358      249095 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     359             : }
     360             : 
     361             : /* x,y t_INT, return x/y in reduced form */
     362             : GEN
     363    88210593 : Qdivii(GEN x, GEN y)
     364             : {
     365    88210593 :   pari_sp av = avma;
     366             :   GEN r, q;
     367             : 
     368    88210593 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     369    50999861 :   if (equali1(x)) {
     370     5953949 :     long s = signe(y);
     371     5953949 :     if (!s) pari_err_INV("gdiv",y);
     372     5953844 :     if (signe(x) < 0) s = -s;
     373     5953844 :     q = cgetg(3, t_FRAC);
     374     5953844 :     gel(q,1) = s<0? gen_m1: gen_1;
     375     5953844 :     gel(q,2) = absi(y); return q;
     376             :   }
     377    45045912 :   q = dvmdii(x,y,&r);
     378    45045911 :   if (r == gen_0) return q; /* gen_0 intended */
     379    25115453 :   r = gcdii(y,r);
     380    25115453 :   if (lgefint(r) == 3)
     381             :   {
     382    24780966 :     ulong t = r[2];
     383    24780966 :     set_avma(av);
     384    24780966 :     if (t == 1) q = mkfraccopy(x,y);
     385             :     else
     386             :     {
     387     8800196 :       q = cgetg(3,t_FRAC);
     388     8800196 :       gel(q,1) = diviuexact(x,t);
     389     8800196 :       gel(q,2) = diviuexact(y,t);
     390             :     }
     391             :   }
     392             :   else
     393             :   { /* rare: r and q left on stack for efficiency */
     394      334487 :     q = cgetg(3,t_FRAC);
     395      334487 :     gel(q,1) = diviiexact(x,r);
     396      334487 :     gel(q,2) = diviiexact(y,r);
     397             :   }
     398    25115455 :   normalize_frac(q); return q;
     399             : }
     400             : 
     401             : /*******************************************************************/
     402             : /*                                                                 */
     403             : /*                          CONJUGATION                            */
     404             : /*                                                                 */
     405             : /*******************************************************************/
     406             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     407             : static GEN
     408       28700 : quad_polmod_conj(GEN x, GEN y)
     409             : {
     410             :   GEN z, u, v, a, b;
     411       28700 :   if (typ(x) != t_POL) return gcopy(x);
     412       28700 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     413       28700 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     414       28700 :   b = gel(y,3); v = gel(x,2);
     415       28700 :   z = cgetg(4, t_POL); z[1] = x[1];
     416       28700 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     417       28700 :   gel(z,3) = gneg(u); return z;
     418             : }
     419             : static GEN
     420       28700 : quad_polmod_norm(GEN x, GEN y)
     421             : {
     422             :   GEN z, u, v, a, b, c;
     423       28700 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     424       28700 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     425       28700 :   b = gel(y,3); v = gel(x,2);
     426       28700 :   c = gel(y,2);
     427       28700 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     428       28700 :   if (!gequal1(a)) z = gdiv(z, a);
     429       28700 :   return gadd(z, gsqr(v));
     430             : }
     431             : 
     432             : GEN
     433    10712682 : conj_i(GEN x)
     434             : {
     435             :   long lx, i;
     436             :   GEN y;
     437             : 
     438    10712682 :   switch(typ(x))
     439             :   {
     440             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     441     5208149 :       return x;
     442             : 
     443     5468032 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     444             : 
     445             :     case t_QUAD:
     446         595 :       y = cgetg(4,t_QUAD);
     447         595 :       gel(y,1) = gel(x,1);
     448        1190 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     449         595 :                                       : gadd(gel(x,2), gel(x,3));
     450         595 :       gel(y,3) = gneg(gel(x,3)); return y;
     451             : 
     452             :     case t_POL: case t_SER:
     453       10647 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     454       10647 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     455       10647 :       return y;
     456             : 
     457             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     458       25259 :       y = cgetg_copy(x, &lx);
     459       25259 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     460       25259 :       return y;
     461             : 
     462             :     case t_POLMOD:
     463             :     {
     464           0 :       GEN X = gel(x,1);
     465           0 :       long d = degpol(X);
     466           0 :       if (d < 2) return x;
     467           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     468             :     }
     469             :   }
     470           0 :   pari_err_TYPE("gconj",x);
     471             :   return NULL; /* LCOV_EXCL_LINE */
     472             : }
     473             : GEN
     474       37086 : gconj(GEN x)
     475       37086 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     476             : 
     477             : GEN
     478          84 : conjvec(GEN x,long prec)
     479             : {
     480             :   long lx, s, i;
     481             :   GEN z;
     482             : 
     483          84 :   switch(typ(x))
     484             :   {
     485             :     case t_INT: case t_INTMOD: case t_FRAC:
     486           0 :       return mkcolcopy(x);
     487             : 
     488             :     case t_COMPLEX: case t_QUAD:
     489           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     490             : 
     491             :     case t_FFELT:
     492          28 :       return FF_conjvec(x);
     493             : 
     494             :     case t_VEC: case t_COL:
     495           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     496           0 :       if (lx == 1) return z;
     497           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     498           0 :       s = lgcols(z);
     499           0 :       for (i=2; i<lx; i++)
     500             :       {
     501           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     502           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     503             :       }
     504           0 :       break;
     505             : 
     506             :     case t_POLMOD: {
     507          56 :       GEN T = gel(x,1), r;
     508             :       pari_sp av;
     509             : 
     510          56 :       lx = lg(T);
     511          56 :       if (lx <= 3) return cgetg(1,t_COL);
     512          56 :       x = gel(x,2);
     513         238 :       for (i=2; i<lx; i++)
     514             :       {
     515         189 :         GEN c = gel(T,i);
     516         189 :         switch(typ(c)) {
     517             :           case t_INTMOD: {
     518           7 :             GEN p = gel(c,1);
     519             :             pari_sp av;
     520           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     521           7 :             av = avma;
     522           7 :             T = RgX_to_FpX(T,p);
     523           7 :             x = RgX_to_FpX(x, p);
     524           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     525           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     526           7 :             return gerepileupto(av, z);
     527             :           }
     528             :           case t_INT:
     529         182 :           case t_FRAC: break;
     530           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     531             :         }
     532             :       }
     533          49 :       if (typ(x) != t_POL)
     534             :       {
     535           0 :         if (!is_rational_t(typ(x)))
     536           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     537           0 :         retconst_col(lx-3, gcopy(x));
     538             :       }
     539          49 :       RgX_check_QX(x,"conjvec");
     540          49 :       av = avma;
     541          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     542          49 :       r = cleanroots(T,prec);
     543          49 :       z = cgetg(lx-2,t_COL);
     544          49 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     545          49 :       return gerepileupto(av, z);
     546             :     }
     547             : 
     548             :     default:
     549           0 :       pari_err_TYPE("conjvec",x);
     550             :       return NULL; /* LCOV_EXCL_LINE */
     551             :   }
     552           0 :   return z;
     553             : }
     554             : 
     555             : 
     556             : /********************************************************************/
     557             : /**                                                                **/
     558             : /**                           ADDITION                             **/
     559             : /**                                                                **/
     560             : /********************************************************************/
     561             : /* x, y compatible PADIC, op = add or sub */
     562             : static GEN
     563     1248762 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     564             : {
     565     1248762 :   pari_sp av = avma;
     566             :   long d,e,r,rx,ry;
     567     1248762 :   GEN u, z, mod, p = gel(x,2);
     568             :   int swap;
     569             : 
     570     1248762 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     571     1248762 :   e = valp(x);
     572     1248762 :   r = valp(y); d = r-e;
     573     1248762 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     574     1248762 :   rx = precp(x);
     575     1248762 :   ry = precp(y);
     576     1248762 :   if (d) /* v(x) < v(y) */
     577             :   {
     578      602899 :     r = d+ry; z = powiu(p,d);
     579      602899 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     580      602899 :     z = mulii(z,gel(y,4));
     581      602899 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     582             :   }
     583             :   else
     584             :   {
     585             :     long c;
     586      645863 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     587      645863 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     588      645863 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     589             :     {
     590       67508 :       set_avma(av); return zeropadic(p, e+r);
     591             :     }
     592      578355 :     if (c)
     593             :     {
     594       99225 :       mod = diviiexact(mod, powiu(p,c));
     595       99225 :       r -= c;
     596       99225 :       e += c;
     597             :     }
     598             :   }
     599     1181254 :   u = modii(u, mod);
     600     1181254 :   set_avma(av); z = cgetg(5,t_PADIC);
     601     1181254 :   z[1] = evalprecp(r) | evalvalp(e);
     602     1181254 :   gel(z,2) = icopy(p);
     603     1181254 :   gel(z,3) = icopy(mod);
     604     1181254 :   gel(z,4) = icopy(u); return z;
     605             : }
     606             : /* Rg_to_Fp(t_FRAC) without GC */
     607             : static GEN
     608       56483 : Q_to_Fp(GEN x, GEN p)
     609       56483 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     610             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     611             : static GEN
     612      375002 : addQp(GEN x, GEN y)
     613             : {
     614      375002 :   pari_sp av = avma;
     615      375002 :   long d, r, e, vy = valp(y), py = precp(y);
     616      375002 :   GEN z, q, mod, u, p = gel(y,2);
     617             : 
     618      375002 :   e = Q_pvalrem(x, p, &x);
     619      374983 :   d = vy - e; r = d + py;
     620      374983 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     621      373401 :   mod = gel(y,3);
     622      373401 :   u   = gel(y,4);
     623      373401 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     624             : 
     625      373420 :   if (d > 0)
     626             :   {
     627      237035 :     q = powiu(p,d);
     628      236999 :     mod = mulii(mod, q);
     629      236990 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     630      236990 :     u = addii(x,  mulii(u, q));
     631             :   }
     632      136385 :   else if (d < 0)
     633             :   {
     634       25243 :     q = powiu(p,-d);
     635       25240 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     636       25240 :     u = addii(u, mulii(x, q));
     637       25242 :     r = py; e = vy;
     638             :   }
     639             :   else
     640             :   {
     641             :     long c;
     642      111142 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     643      111142 :     u = addii(u, x);
     644      111140 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     645             :     {
     646         966 :       set_avma(av); return zeropadic(p,e+r);
     647             :     }
     648      110177 :     if (c)
     649             :     {
     650       59911 :       mod = diviiexact(mod, powiu(p,c));
     651       59910 :       r -= c;
     652       59910 :       e += c;
     653             :     }
     654             :   }
     655      372411 :   u = modii(u, mod); set_avma(av);
     656      372412 :   z = cgetg(5,t_PADIC);
     657      372444 :   z[1] = evalprecp(r) | evalvalp(e);
     658      372442 :   gel(z,2) = icopy(p);
     659      372452 :   gel(z,3) = icopy(mod);
     660      372455 :   gel(z,4) = icopy(u); return z;
     661             : }
     662             : 
     663             : /* Mod(x,X) + Mod(y,X) */
     664             : #define addsub_polmod_same addsub_polmod_scal
     665             : /* Mod(x,X) +/- Mod(y,Y) */
     666             : static GEN
     667        1750 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     668             : {
     669        1750 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     670        1750 :   GEN z = cgetg(3,t_POLMOD);
     671        1750 :   long vx = varn(X), vy = varn(Y);
     672        1750 :   if (vx==vy) {
     673             :     pari_sp av;
     674          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     675          14 :     warn_coercion(X,Y,gel(z,1));
     676          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     677             :   }
     678        1736 :   if (varncmp(vx, vy) < 0)
     679        1736 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     680             :   else
     681           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     682        1736 :   gel(z,2) = op(x, y); return z;
     683             : }
     684             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     685             : static GEN
     686    12343864 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     687    12343864 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     688             : 
     689             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     690             : static GEN
     691      308126 : add_ser_scal(GEN y, GEN x)
     692             : {
     693             :   long i, l, ly, vy;
     694             :   GEN z;
     695             : 
     696      308126 :   if (isrationalzero(x)) return gcopy(y);
     697      281666 :   ly = lg(y);
     698      281666 :   l = valp(y);
     699      281666 :   if (l < 3-ly) return gcopy(y);
     700             :   /* l + ly >= 3 */
     701      281498 :   if (l < 0)
     702             :   {
     703         735 :     z = cgetg(ly,t_SER); z[1] = y[1];
     704         735 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     705         735 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     706         735 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     707         735 :     return z;
     708             :   }
     709      280763 :   vy = varn(y);
     710      280763 :   if (l > 0)
     711             :   {
     712       21455 :     if (ser_isexactzero(y))
     713       13867 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     714        7588 :     y -= l; ly += l;
     715        7588 :     z = cgetg(ly,t_SER);
     716        7588 :     x = gcopy(x);
     717        7588 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     718             :   }
     719             :   else
     720             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     721      259308 :     z = cgetg(ly,t_SER);
     722      259308 :     x = gadd(x, gel(y,2));
     723      259308 :     i = 3;
     724             :   }
     725      266896 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     726      266896 :   gel(z,2) = x;
     727      266896 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     728      266896 :   return gequal0(x)? normalize(z): z;
     729             : }
     730             : static long
     731     3542740 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     732             : /* x,y t_SER in the same variable: x+y */
     733             : static GEN
     734     1771762 : ser_add(GEN x, GEN y)
     735             : {
     736     1771762 :   long i, lx,ly, n = valp(y) - valp(x);
     737             :   GEN z;
     738     1771762 :   if (n < 0) { n = -n; swap(x,y); }
     739             :   /* valp(x) <= valp(y) */
     740     1771762 :   lx = _serprec(x);
     741     1771762 :   if (lx == 2) /* don't lose type information */
     742         784 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
     743     1770978 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     744     1770978 :   if (n)
     745             :   {
     746       91080 :     if (n+2 > lx) return gcopy(x);
     747       90527 :     z = cgetg(ly,t_SER);
     748       90527 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     749       90527 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     750             :   } else {
     751     1679898 :     z = cgetg(ly,t_SER);
     752     1679898 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     753             :   }
     754     1770425 :   z[1] = x[1]; return normalize(z);
     755             : }
     756             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     757             : static GEN
     758     8823567 : add_rfrac_scal(GEN y, GEN x)
     759             : {
     760             :   pari_sp av;
     761             :   GEN n;
     762             : 
     763     8823567 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     764     5149291 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     765     5149291 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     766             : }
     767             : 
     768             : /* x "scalar", ty != t_MAT and non-scalar */
     769             : static GEN
     770    19827219 : add_scal(GEN y, GEN x, long ty)
     771             : {
     772    19827219 :   switch(ty)
     773             :   {
     774    12190851 :     case t_POL: return RgX_Rg_add(y, x);
     775      308112 :     case t_SER: return add_ser_scal(y, x);
     776     4470269 :     case t_RFRAC: return add_rfrac_scal(y, x);
     777           0 :     case t_COL: return RgC_Rg_add(y, x);
     778             :     case t_VEC:
     779     2857981 :       if (isintzero(x)) return gcopy(y);
     780         168 :       break;
     781             :   }
     782         174 :   pari_err_TYPE2("+",x,y);
     783             :   return NULL; /* LCOV_EXCL_LINE */
     784             : }
     785             : 
     786             : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
     787             : static GEN
     788    10571834 : setfrac(GEN z, GEN a, GEN b)
     789             : {
     790    10571834 :   gel(z,1) = icopy_avma(a, (pari_sp)z);
     791    10571834 :   gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
     792    10571834 :   set_avma((pari_sp)gel(z,2)); return z;
     793             : }
     794             : /* z <- a / (b*Q), (Q,a) = 1 */
     795             : static GEN
     796     9450536 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
     797             : {
     798     9450536 :   GEN q = Qdivii(a, b); /* != 0 */
     799     9450536 :   if (typ(q) == t_INT)
     800             :   {
     801      766934 :     gel(z,1) = gerepileuptoint((pari_sp)Q, q);
     802      766934 :     gel(z,2) = Q; return z;
     803             :   }
     804     8683602 :   return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
     805             : }
     806             : static GEN
     807    19172890 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     808             : {
     809    19172890 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     810    19172890 :   GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
     811    19172890 :   int s = cmpii(x2, y2);
     812             : 
     813             :   /* common denominator: (x1 op y1) / x2 */
     814    19172890 :   if (!s)
     815             :   {
     816     6601692 :     pari_sp av = avma;
     817     6601692 :     return gerepileupto(av, Qdivii(op(x1, y1), x2));
     818             :   }
     819    12571198 :   z = cgetg(3, t_FRAC);
     820    12571198 :   if (s < 0)
     821             :   {
     822     6456519 :     Q = dvmdii(y2, x2, &r);
     823             :     /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     824     6456519 :     if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
     825     2052289 :     delta = gcdii(x2,r);
     826             :   }
     827             :   else
     828             :   {
     829     6114679 :     Q = dvmdii(x2, y2, &r);
     830             :     /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     831     6114679 :     if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
     832     1068373 :     delta = gcdii(y2,r);
     833             :   }
     834             :   /* delta = gcd(x2,y2) */
     835     3120662 :   if (equali1(delta))
     836             :   { /* numerator is non-zero */
     837     1232430 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     838     1232430 :     gel(z,2) = mulii(x2,y2); return z;
     839             :   }
     840     1888232 :   x2 = diviiexact(x2,delta);
     841     1888232 :   y2 = diviiexact(y2,delta);
     842     1888232 :   n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
     843     1888232 :   q = dvmdii(n, delta, &r);
     844     1888232 :   if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
     845     1746014 :   r = gcdii(delta, r);
     846     1746014 :   if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
     847     1746014 :   return setfrac(z, n, mulii(mulii(x2, y2), delta));
     848             : }
     849             : 
     850             : /* assume x2, y2 are t_POLs in the same variable */
     851             : static GEN
     852     3037318 : add_rfrac(GEN x, GEN y)
     853             : {
     854     3037318 :   pari_sp av = avma;
     855     3037318 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     856     3037318 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     857             : 
     858     3037318 :   delta = RgX_gcd(x2,y2);
     859     3037318 :   if (!degpol(delta))
     860             :   {
     861         658 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     862         658 :     d = RgX_mul(x2, y2);
     863         658 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     864             :   }
     865     3036660 :   x2 = RgX_div(x2,delta);
     866     3036660 :   y2 = RgX_div(y2,delta);
     867     3036660 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     868     3036660 :   if (!signe(n))
     869             :   {
     870      721054 :     n = simplify_shallow(n);
     871      721054 :     if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
     872           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     873             :   }
     874     2315606 :   if (degpol(n) == 0)
     875     1149721 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     876     1165885 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     877     1165885 :   if (isexactzero(r))
     878             :   {
     879             :     GEN z;
     880      228052 :     d = RgX_mul(x2, y2);
     881             :     /* "constant" denominator ? */
     882      228052 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     883      228052 :     return gerepileupto(av, z);
     884             :   }
     885      937833 :   r = RgX_gcd(delta, r);
     886      937833 :   if (degpol(r))
     887             :   {
     888      160705 :     n = RgX_div(n, r);
     889      160705 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     890             :   }
     891             :   else
     892      777128 :     d = RgX_mul(gel(x,2), y2);
     893      937833 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     894             : }
     895             : 
     896             : GEN
     897  2227870657 : gadd(GEN x, GEN y)
     898             : {
     899  2227870657 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     900             :   pari_sp av, tetpil;
     901             :   GEN z, p1;
     902             : 
     903  2227870657 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     904             :   {
     905  1210351861 :     case t_INT: return addii(x,y);
     906   593650592 :     case t_REAL: return addrr(x,y);
     907     1525664 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     908     1525664 :       z = cgetg(3,t_INTMOD);
     909     1525662 :       if (X==Y || equalii(X,Y))
     910     1525611 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     911          14 :       gel(z,1) = gcdii(X,Y);
     912          14 :       warn_coercion(X,Y,gel(z,1));
     913          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     914          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     915             :     }
     916    14326397 :     case t_FRAC: return addsub_frac(x,y,addii);
     917    95127775 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     918    95127775 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     919    95127775 :       if (isintzero(gel(z,2)))
     920             :       {
     921      360500 :         set_avma((pari_sp)(z+3));
     922      360500 :         return gadd(gel(x,1),gel(y,1));
     923             :       }
     924    94767275 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     925    94767275 :       return z;
     926             :     case t_PADIC:
     927     1019900 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     928     1019900 :       return addsub_pp(x,y, addii);
     929         476 :     case t_QUAD: z = cgetg(4,t_QUAD);
     930         476 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     931         476 :       gel(z,1) = ZX_copy(gel(x,1));
     932         476 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     933         476 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     934             :     case t_POLMOD:
     935     2339302 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     936     2337587 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     937        1715 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     938    14133655 :     case t_FFELT: return FF_add(x,y);
     939             :     case t_POL:
     940    19526511 :       vx = varn(x);
     941    19526511 :       vy = varn(y);
     942    19526511 :       if (vx != vy) {
     943      806833 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     944       18548 :         else                     return RgX_Rg_add(y, x);
     945             :       }
     946    18719678 :       return RgX_add(x, y);
     947             :     case t_SER:
     948     1766946 :       vx = varn(x);
     949     1766946 :       vy = varn(y);
     950     1766946 :       if (vx != vy) {
     951          14 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
     952          14 :         else                     return add_ser_scal(y, x);
     953             :       }
     954     1766932 :       return ser_add(x, y);
     955             :     case t_RFRAC:
     956     4329205 :       vx = varn(gel(x,2));
     957     4329205 :       vy = varn(gel(y,2));
     958     4329205 :       if (vx != vy) {
     959     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
     960      538397 :         else                     return add_rfrac_scal(y, x);
     961             :       }
     962     3037318 :       return add_rfrac(x,y);
     963             :     case t_VEC:
     964      614058 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
     965      614058 :       return RgV_add(x,y);
     966             :     case t_COL:
     967      871602 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
     968      871602 :       return RgC_add(x,y);
     969             :     case t_MAT:
     970      671902 :       lx = lg(x);
     971      671902 :       if (lg(y) != lx) pari_err_OP("+",x,y);
     972      671902 :       if (lx == 1) return cgetg(1, t_MAT);
     973      671902 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
     974      671895 :       return RgM_add(x,y);
     975             : 
     976           0 :     default: pari_err_TYPE2("+",x,y);
     977             :   }
     978             :   /* tx != ty */
     979   267614979 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
     980             : 
     981   267614979 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
     982             :   {
     983             :     case t_INT:
     984   211153434 :       switch(ty)
     985             :       {
     986   119744219 :         case t_REAL: return addir(x,y);
     987             :         case t_INTMOD:
     988      714702 :           z = cgetg(3, t_INTMOD);
     989      714702 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
     990    27292902 :         case t_FRAC: z = cgetg(3,t_FRAC);
     991    27292902 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
     992    27292902 :           gel(z,2) = icopy(gel(y,2)); return z;
     993    61405980 :         case t_COMPLEX: return addRc(x, y);
     994             :         case t_PADIC:
     995      336436 :           if (!signe(x)) return gcopy(y);
     996      304247 :           return addQp(x,y);
     997         721 :         case t_QUAD: return addRq(x, y);
     998     1658474 :         case t_FFELT: return FF_Z_add(y,x);
     999             :       }
    1000             : 
    1001             :     case t_REAL:
    1002    22780450 :       switch(ty)
    1003             :       {
    1004             :         case t_FRAC:
    1005    12653035 :           if (!signe(gel(y,1))) return rcopy(x);
    1006    12653035 :           if (!signe(x))
    1007             :           {
    1008        4839 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1009        4839 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1010             :           }
    1011    12648196 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1012    12648196 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1013    10127394 :         case t_COMPLEX: return addRc(x, y);
    1014          21 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1015             : 
    1016           0 :         default: pari_err_TYPE2("+",x,y);
    1017             :       }
    1018             : 
    1019             :     case t_INTMOD:
    1020       17640 :       switch(ty)
    1021             :       {
    1022       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1023       17507 :           z = cgetg(3, t_INTMOD);
    1024       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1025       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1026             :         }
    1027             :         case t_FFELT:
    1028          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1029           0 :             pari_err_OP("+",x,y);
    1030          14 :           return FF_Z_add(y,gel(x,2));
    1031          91 :         case t_COMPLEX: return addRc(x, y);
    1032           0 :         case t_PADIC: { GEN X = gel(x,1);
    1033           0 :           z = cgetg(3, t_INTMOD);
    1034           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1035             :         }
    1036          28 :         case t_QUAD: return addRq(x, y);
    1037             :       }
    1038             : 
    1039             :     case t_FRAC:
    1040     9986890 :       switch (ty)
    1041             :       {
    1042     9914783 :         case t_COMPLEX: return addRc(x, y);
    1043             :         case t_PADIC:
    1044       70756 :           if (!signe(gel(x,1))) return gcopy(y);
    1045       70756 :           return addQp(x,y);
    1046          14 :         case t_QUAD: return addRq(x, y);
    1047        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1048             :       }
    1049             : 
    1050             :     case t_FFELT:
    1051           2 :       pari_err_TYPE2("+",x,y);
    1052             : 
    1053             :     case t_COMPLEX:
    1054          21 :       switch(ty)
    1055             :       {
    1056             :         case t_PADIC:
    1057          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1058             :         case t_QUAD:
    1059           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1060           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1061             :       }
    1062             : 
    1063             :     case t_PADIC: /* ty == t_QUAD */
    1064           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1065             :   }
    1066             :   /* tx < ty, !is_const_t(y) */
    1067    23676539 :   switch(ty)
    1068             :   {
    1069             :     case t_MAT:
    1070       65604 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1071       65604 :       if (isrationalzero(x)) return gcopy(y);
    1072        5488 :       return RgM_Rg_add(y, x);
    1073             :     case t_COL:
    1074      336444 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1075      336444 :       return RgC_Rg_add(y, x);
    1076             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1077      335146 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1078             :   }
    1079    22939345 :   if (is_scalar_t(tx))  {
    1080    19850073 :     if (tx == t_POLMOD)
    1081             :     {
    1082       84287 :       vx = varn(gel(x,1));
    1083       84287 :       vy = gvar(y);
    1084       84287 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1085             :       else
    1086       49427 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1087       45451 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1088             :     }
    1089    19765786 :     return add_scal(y, x, ty);
    1090             :   }
    1091             :   /* x and y are not scalars, ty != t_MAT */
    1092     3089273 :   vx = gvar(x);
    1093     3089273 :   vy = gvar(y);
    1094     3089273 :   if (vx != vy) { /* x or y is treated as a scalar */
    1095       22605 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1096       22598 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1097       22598 :                                 : add_scal(y, x, ty);
    1098             :   }
    1099             :   /* vx = vy */
    1100     3066668 :   switch(tx)
    1101             :   {
    1102             :     case t_POL:
    1103     3066262 :       switch (ty)
    1104             :       {
    1105             :         case t_SER:
    1106        4851 :           if (lg(x) == 2) return gcopy(y);
    1107        4837 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1108        4837 :           i = lg(y) + valp(y) - i;
    1109        4837 :           if (i < 3) return gcopy(y);
    1110        4830 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1111        4830 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1112        4830 :           return y;
    1113             : 
    1114     3061411 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1115             :       }
    1116           0 :       break;
    1117             : 
    1118             :     case t_SER:
    1119         406 :       if (ty == t_RFRAC)
    1120             :       {
    1121             :         GEN n, d;
    1122             :         long vn, vd;
    1123         406 :         av = avma;
    1124         406 :         n = gel(y,1); vn = gval(n, vy);
    1125         406 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1126             : 
    1127         406 :         l = lg(x) + valp(x) - (vn - vd);
    1128         406 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1129             : 
    1130             :         /* take advantage of y = t^n ! */
    1131         406 :         if (degpol(d))
    1132          70 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1133             :         else {
    1134         336 :           y = gdiv(n, gel(d,2));
    1135         336 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1136             :         }
    1137         406 :         setvalp(y, valp(y) - vd);
    1138         406 :         return gerepileupto(av, gadd(y, x));
    1139             :       }
    1140           0 :       break;
    1141             :   }
    1142           0 :   pari_err_TYPE2("+",x,y);
    1143             :   return NULL; /* LCOV_EXCL_LINE */
    1144             : }
    1145             : 
    1146             : GEN
    1147    33652116 : gaddsg(long x, GEN y)
    1148             : {
    1149    33652116 :   long ty = typ(y);
    1150             :   GEN z;
    1151             : 
    1152    33652116 :   switch(ty)
    1153             :   {
    1154    12387759 :     case t_INT:  return addsi(x,y);
    1155     8124409 :     case t_REAL: return addsr(x,y);
    1156             :     case t_INTMOD:
    1157          14 :       z = cgetg(3, t_INTMOD);
    1158          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1159     5002823 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1160     5002823 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1161     5002823 :       gel(z,2) = icopy(gel(y,2)); return z;
    1162             :     case t_COMPLEX:
    1163     7541653 :       z = cgetg(3, t_COMPLEX);
    1164     7541653 :       gel(z,1) = gaddsg(x, gel(y,1));
    1165     7541653 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1166             : 
    1167      595458 :     default: return gadd(stoi(x), y);
    1168             :   }
    1169             : }
    1170             : 
    1171             : GEN
    1172     1054748 : gsubsg(long x, GEN y)
    1173             : {
    1174             :   GEN z, a, b;
    1175             :   pari_sp av;
    1176             : 
    1177     1054748 :   switch(typ(y))
    1178             :   {
    1179      101399 :     case t_INT:  return subsi(x,y);
    1180      562488 :     case t_REAL: return subsr(x,y);
    1181             :     case t_INTMOD:
    1182          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1183          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1184       42160 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1185       42160 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1186       42160 :       gel(z,2) = icopy(gel(y,2)); return z;
    1187             :     case t_COMPLEX:
    1188      320512 :       z = cgetg(3, t_COMPLEX);
    1189      320512 :       gel(z,1) = gsubsg(x, gel(y,1));
    1190      320512 :       gel(z,2) = gneg(gel(y,2)); return z;
    1191             :   }
    1192       28140 :   av = avma;
    1193       28140 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1194             : }
    1195             : 
    1196             : /********************************************************************/
    1197             : /**                                                                **/
    1198             : /**                          SUBTRACTION                           **/
    1199             : /**                                                                **/
    1200             : /********************************************************************/
    1201             : 
    1202             : GEN
    1203  1499166309 : gsub(GEN x, GEN y)
    1204             : {
    1205  1499166309 :   long tx = typ(x), ty = typ(y);
    1206             :   pari_sp av;
    1207             :   GEN z;
    1208  1499166309 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1209             :   {
    1210  1197893858 :     case t_INT: return subii(x,y);
    1211   218781563 :     case t_REAL: return subrr(x,y);
    1212     1109318 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1213     1109318 :       z = cgetg(3,t_INTMOD);
    1214     1109316 :       if (X==Y || equalii(X,Y))
    1215     1109289 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1216          14 :       gel(z,1) = gcdii(X,Y);
    1217          14 :       warn_coercion(X,Y,gel(z,1));
    1218          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1219          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1220             :     }
    1221     4846493 :     case t_FRAC: return addsub_frac(x,y, subii);
    1222    28573736 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1223    28573736 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1224    28573736 :       if (isintzero(gel(z,2)))
    1225             :       {
    1226       81128 :         set_avma((pari_sp)(z+3));
    1227       81128 :         return gsub(gel(x,1),gel(y,1));
    1228             :       }
    1229    28492608 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1230    28492608 :       return z;
    1231             :     case t_PADIC:
    1232      228855 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1233      228855 :       return addsub_pp(x,y, subii);
    1234         595 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1235         595 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1236         595 :       gel(z,1) = ZX_copy(gel(x,1));
    1237         595 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1238         595 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1239             :     case t_POLMOD:
    1240     9625715 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1241     9625680 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1242          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1243      283094 :     case t_FFELT: return FF_sub(x,y);
    1244             :     case t_POL: {
    1245    16386511 :       long vx = varn(x);
    1246    16386511 :       long vy = varn(y);
    1247    16386511 :       if (vx != vy) {
    1248       12699 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1249        2072 :         else                     return Rg_RgX_sub(x, y);
    1250             :       }
    1251    16373812 :       return RgX_sub(x, y);
    1252             :     }
    1253             :     case t_VEC:
    1254      271055 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1255      271055 :       return RgV_sub(x,y);
    1256             :     case t_COL:
    1257     3319575 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1258     3319575 :       return RgC_sub(x,y);
    1259             :     case t_MAT: {
    1260       13994 :       long lx = lg(x);
    1261       13994 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1262       13994 :       if (lx == 1) return cgetg(1, t_MAT);
    1263       13994 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1264       13994 :       return RgM_sub(x,y);
    1265             :     }
    1266     2096060 :     case t_RFRAC: case t_SER: break;
    1267             : 
    1268           0 :     default: pari_err_TYPE2("+",x,y);
    1269             :   }
    1270    17831947 :   av = avma;
    1271    17831947 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1272             : }
    1273             : 
    1274             : /********************************************************************/
    1275             : /**                                                                **/
    1276             : /**                        MULTIPLICATION                          **/
    1277             : /**                                                                **/
    1278             : /********************************************************************/
    1279             : static GEN
    1280      516810 : mul_ser_scal(GEN y, GEN x) {
    1281             :   long ly, i;
    1282             :   GEN z;
    1283      516810 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1284      513562 :   if (ser_isexactzero(y))
    1285             :   {
    1286         413 :     if (lg(y) == 2) return gcopy(y);
    1287           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1288             :   }
    1289      513149 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1290      513149 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1291      513149 :   return normalize(z);
    1292             : }
    1293             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1294             :  * [n/d a valid RFRAC]  */
    1295             : static GEN
    1296    10445004 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1297             : {
    1298    10445004 :   pari_sp av = avma;
    1299             :   GEN z;
    1300             : 
    1301    10445004 :   switch(typ(x))
    1302             :   {
    1303             :     case t_PADIC:
    1304          14 :       n = gmul(n, x);
    1305          14 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1306          14 :       return gerepileupto(av, gdiv(n,d));
    1307             : 
    1308             :     case t_INTMOD: case t_POLMOD:
    1309         434 :       n = gmul(n, x);
    1310         434 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1311         434 :       return gerepileupto(av, gdiv(n,d));
    1312             :   }
    1313    10444556 :   z = gred_rfrac2(x, d);
    1314    10444556 :   n = simplify_shallow(n);
    1315    10444556 :   if (typ(z) == t_RFRAC)
    1316             :   {
    1317     7928244 :     n = gmul(gel(z,1), n);
    1318     7928244 :     d = gel(z,2);
    1319     7928244 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1320           0 :       z = RgX_Rg_div(n, d);
    1321             :     else
    1322     7928244 :       z = gred_rfrac_simple(n, d);
    1323             :   }
    1324             :   else
    1325     2516312 :     z = gmul(z, n);
    1326    10444556 :   return gerepileupto(av, z);
    1327             : }
    1328             : static GEN
    1329    67586585 : mul_scal(GEN y, GEN x, long ty)
    1330             : {
    1331    67586585 :   switch(ty)
    1332             :   {
    1333             :     case t_POL:
    1334    58624517 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1335    54373811 :       return RgX_Rg_mul(y, x);
    1336      382077 :     case t_SER: return mul_ser_scal(y, x);
    1337     8579963 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1338             :     case t_QFI: case t_QFR:
    1339          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1340             :   }
    1341          14 :   pari_err_TYPE2("*",x,y);
    1342             :   return NULL; /* LCOV_EXCL_LINE */
    1343             : }
    1344             : 
    1345             : static GEN
    1346      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1347             : {
    1348      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1349      160830 :   long vx = gvar(X), vy = varn(y2);
    1350      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1351        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1352             : }
    1353             : /* (x1/x2) * (y1/y2) */
    1354             : static GEN
    1355     7904181 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1356             : {
    1357             :   GEN z, X, Y;
    1358     7904181 :   pari_sp av = avma;
    1359             : 
    1360     7904181 :   X = gred_rfrac2(x1, y2);
    1361     7904181 :   Y = gred_rfrac2(y1, x2);
    1362     7904181 :   if (typ(X) == t_RFRAC)
    1363             :   {
    1364     6625218 :     if (typ(Y) == t_RFRAC) {
    1365     6559642 :       x1 = gel(X,1);
    1366     6559642 :       x2 = gel(X,2);
    1367     6559642 :       y1 = gel(Y,1);
    1368     6559642 :       y2 = gel(Y,2);
    1369     6559642 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1370             :     } else
    1371       65576 :       z = mul_gen_rfrac(Y, X);
    1372             :   }
    1373     1278963 :   else if (typ(Y) == t_RFRAC)
    1374       95254 :     z = mul_gen_rfrac(X, Y);
    1375             :   else
    1376     1183709 :     z = gmul(X, Y);
    1377     7904181 :   return gerepileupto(av, z);
    1378             : }
    1379             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1380             : static GEN
    1381      269855 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1382             : {
    1383      269855 :   pari_sp av = avma;
    1384      269855 :   GEN X = gred_rfrac2(x1, y2);
    1385      269855 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1386             :   {
    1387      263303 :     x2 = RgX_mul(gel(X,2), x2);
    1388      263303 :     x1 = gel(X,1);
    1389             :   }
    1390             :   else
    1391        6552 :     x1 = X;
    1392      269855 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1393             : }
    1394             : 
    1395             : /* Mod(y, Y) * x,  assuming x scalar */
    1396             : static GEN
    1397     4497271 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1398     4497271 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1399             : 
    1400             : /* cf mulqq */
    1401             : static GEN
    1402     2684206 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1403             : {
    1404     2684206 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1405     2684206 :   pari_sp tetpil, av = avma;
    1406     2684206 :   T[1] = x[1];
    1407     2684206 :   p2 = gmul(gel(x,2), gel(y,2));
    1408     2684206 :   p3 = gmul(gel(x,3), gel(y,3));
    1409     2684206 :   p1 = gmul(gneg_i(c),p3);
    1410             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1411     2684206 :   if (typ(b) == t_INT)
    1412             :   {
    1413     2684185 :     if (signe(b))
    1414             :     {
    1415     2072777 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1416     2072777 :       if (is_pm1(b))
    1417             :       {
    1418     2072119 :         if (signe(b) > 0) p3 = gneg(p3);
    1419             :       }
    1420             :       else
    1421         658 :         p3 = gmul(negi(b), p3);
    1422             :     }
    1423             :     else
    1424             :     {
    1425      611408 :       p3 = gmul(gel(x,2),gel(y,3));
    1426      611408 :       p4 = gmul(gel(x,3),gel(y,2));
    1427             :     }
    1428             :   }
    1429             :   else
    1430             :   {
    1431          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1432          21 :     p3 = gmul(gneg_i(b), p3);
    1433             :   }
    1434     2684206 :   tetpil = avma;
    1435     2684206 :   gel(T,2) = gadd(p2, p1);
    1436     2684206 :   gel(T,3) = gadd(p4, p3);
    1437     2684206 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1438     2684206 :   return normalizepol_lg(T,4);
    1439             : }
    1440             : /* Mod(x,T) * Mod(y,T) */
    1441             : static GEN
    1442     7051145 : mul_polmod_same(GEN T, GEN x, GEN y)
    1443             : {
    1444     7051145 :   GEN z = cgetg(3,t_POLMOD), a;
    1445     7051145 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1446     7051145 :   gel(z,1) = RgX_copy(T);
    1447             :   /* x * y mod T optimised */
    1448     7051145 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1449     5093676 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1450     3956700 :     a = gmul(x, y);
    1451             :   else
    1452             :   {
    1453     3094445 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1454     2678088 :       a = quad_polmod_mul(T, x, y);
    1455             :     else
    1456      416357 :       a = RgXQ_mul(x, y, gel(z,1));
    1457             :   }
    1458     7051145 :   gel(z,2) = a; return z;
    1459             : }
    1460             : static GEN
    1461       43246 : sqr_polmod(GEN T, GEN x)
    1462             : {
    1463       43246 :   GEN a, z = cgetg(3,t_POLMOD);
    1464       43246 :   gel(z,1) = RgX_copy(T);
    1465       43246 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1466       14098 :     a = gsqr(x);
    1467             :   else
    1468             :   {
    1469       29148 :     pari_sp av = avma;
    1470       29148 :     a = RgXQ_sqr(x, gel(z,1));
    1471       29148 :     a = gerepileupto(av, a);
    1472             :   }
    1473       43246 :   gel(z,2) = a; return z;
    1474             : }
    1475             : /* Mod(x,X) * Mod(y,Y) */
    1476             : static GEN
    1477     4301353 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1478             : {
    1479     4301353 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1480     4301353 :   long vx = varn(X), vy = varn(Y);
    1481     4301353 :   GEN z = cgetg(3,t_POLMOD);
    1482             : 
    1483     4301353 :   if (vx==vy) {
    1484             :     pari_sp av;
    1485          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1486          14 :     warn_coercion(X,Y,gel(z,1));
    1487          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1488          14 :     return z;
    1489             :   }
    1490     4301339 :   if (varncmp(vx, vy) < 0)
    1491     4251527 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1492             :   else
    1493       49812 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1494     4301339 :   gel(z,2) = gmul(x, y); return z;
    1495             : }
    1496             : 
    1497             : #if 0 /* used by 3M only */
    1498             : /* set z = x+y and return 1 if x,y have the same sign
    1499             :  * set z = x-y and return 0 otherwise */
    1500             : static int
    1501             : did_add(GEN x, GEN y, GEN *z)
    1502             : {
    1503             :   long tx = typ(x), ty = typ(y);
    1504             :   if (tx == ty) switch(tx)
    1505             :   {
    1506             :     case t_INT: *z = addii(x,y); return 1;
    1507             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1508             :     case t_REAL:
    1509             :       if (signe(x) == -signe(y))
    1510             :       { *z = subrr(x,y); return 0; }
    1511             :       else
    1512             :       { *z = addrr(x,y); return 1; }
    1513             :   }
    1514             :   if (tx == t_REAL) switch(ty)
    1515             :   {
    1516             :     case t_INT:
    1517             :       if (signe(x) == -signe(y))
    1518             :       { *z = subri(x,y); return 0; }
    1519             :       else
    1520             :       { *z = addri(x,y); return 1; }
    1521             :     case t_FRAC:
    1522             :       if (signe(x) == -signe(gel(y,1)))
    1523             :       { *z = gsub(x,y); return 0; }
    1524             :       else
    1525             :       { *z = gadd(x,y); return 1; }
    1526             :   }
    1527             :   else if (ty == t_REAL) switch(tx)
    1528             :   {
    1529             :     case t_INT:
    1530             :       if (signe(x) == -signe(y))
    1531             :       { *z = subir(x,y); return 0; }
    1532             :       else
    1533             :       { *z = addir(x,y); return 1; }
    1534             :     case t_FRAC:
    1535             :       if (signe(gel(x,1)) == -signe(y))
    1536             :       { *z = gsub(x,y); return 0; }
    1537             :       else
    1538             :       { *z = gadd(x,y); return 1; }
    1539             :   }
    1540             :   *z = gadd(x,y); return 1;
    1541             : }
    1542             : #endif
    1543             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1544             : static GEN
    1545     4759322 : mulcIR(GEN x, GEN y)
    1546             : {
    1547     4759322 :   GEN z = cgetg(3,t_COMPLEX);
    1548     4759322 :   pari_sp av = avma;
    1549     4759322 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1550     4759322 :   gel(z,2) = gmul(y, gel(x,1));
    1551     4759322 :   return z;
    1552             : 
    1553             : }
    1554             : /* x,y COMPLEX */
    1555             : static GEN
    1556   125268194 : mulcc(GEN x, GEN y)
    1557             : {
    1558   125268194 :   GEN xr = gel(x,1), xi = gel(x,2);
    1559   125268194 :   GEN yr = gel(y,1), yi = gel(y,2);
    1560             :   GEN p1, p2, p3, p4, z;
    1561             :   pari_sp tetpil, av;
    1562             : 
    1563   125268194 :   if (isintzero(xr))
    1564             :   {
    1565     4633321 :     if (isintzero(yr)) {
    1566      634015 :       av = avma;
    1567      634015 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1568             :     }
    1569     3999306 :     return mulcIR(y, xi);
    1570             :   }
    1571   120634873 :   if (isintzero(yr)) return mulcIR(x, yi);
    1572             : 
    1573   119874857 :   z = cgetg(3,t_COMPLEX); av = avma;
    1574             : #if 0
    1575             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1576             :    * e.g. xr + xi if exponents differ */
    1577             :   if (did_add(xr, xi, &p3))
    1578             :   {
    1579             :     if (did_add(yr, yi, &p4)) {
    1580             :     /* R = xr*yr - xi*yi
    1581             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1582             :       p1 = gmul(xr,yr);
    1583             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1584             :       p3 = gmul(p3, p4);
    1585             :       p4 = gsub(p2, p1);
    1586             :     } else {
    1587             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1588             :      * I = xr*yi + xi*yr */
    1589             :       p1 = gmul(p3,p4);
    1590             :       p3 = gmul(xr,yi);
    1591             :       p4 = gmul(xi,yr);
    1592             :       p2 = gsub(p3, p4);
    1593             :     }
    1594             :   } else {
    1595             :     if (did_add(yr, yi, &p4)) {
    1596             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1597             :       * I = xr*yi +xi*yr */
    1598             :       p1 = gmul(p3,p4);
    1599             :       p3 = gmul(xr,yi);
    1600             :       p4 = gmul(xi,yr);
    1601             :       p2 = gsub(p4, p3);
    1602             :     } else {
    1603             :     /* R = xr*yr - xi*yi
    1604             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1605             :       p3 = gneg( gmul(p3, p4) );
    1606             :       p1 = gmul(xr,yr);
    1607             :       p2 = gmul(xi,yi);
    1608             :       p4 = gadd(p1, p2);
    1609             : 
    1610             :       p2 = gneg(p2);
    1611             :     }
    1612             :   }
    1613             :   tetpil = avma;
    1614             :   gel(z,1) = gadd(p1,p2);
    1615             :   gel(z,2) = gadd(p3,p4);
    1616             : #else
    1617   119874857 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1618             :   { /* 3M formula */
    1619      615927 :     p3 = addii(xr,xi);
    1620      615927 :     p4 = addii(yr,yi);
    1621      615927 :     p1 = mulii(xr,yr);
    1622      615927 :     p2 = mulii(xi,yi);
    1623      615927 :     p3 = mulii(p3,p4);
    1624      615927 :     p4 = addii(p2,p1);
    1625      615927 :     tetpil = avma;
    1626      615927 :     gel(z,1) = subii(p1,p2);
    1627      615927 :     gel(z,2) = subii(p3,p4);
    1628     1118888 :     if (!signe(gel(z,2)))
    1629      112966 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1630             :   }
    1631             :   else
    1632             :   { /* naive 4M formula: avoid all loss of accuracy */
    1633   119258930 :     p1 = gmul(xr,yr);
    1634   119258930 :     p2 = gmul(xi,yi);
    1635   119258930 :     p3 = gmul(xr,yi);
    1636   119258930 :     p4 = gmul(xi,yr);
    1637   119258930 :     tetpil = avma;
    1638   119258930 :     gel(z,1) = gsub(p1,p2);
    1639   119258930 :     gel(z,2) = gadd(p3,p4);
    1640   119258930 :     if (isintzero(gel(z,2)))
    1641             :     {
    1642       16661 :       cgiv(gel(z,2));
    1643       16661 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1644             :     }
    1645             :   }
    1646             : #endif
    1647             : 
    1648   119745230 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1649             : }
    1650             : /* x,y PADIC */
    1651             : static GEN
    1652     1857930 : mulpp(GEN x, GEN y) {
    1653     1857930 :   long l = valp(x) + valp(y);
    1654             :   pari_sp av;
    1655             :   GEN z, t;
    1656     1857930 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1657     1857935 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1658     1659302 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1659             : 
    1660     1579996 :   t = (precp(x) > precp(y))? y: x;
    1661     1579996 :   z = cgetp(t); setvalp(z,l); av = avma;
    1662     1579986 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1663     1579964 :   set_avma(av); return z;
    1664             : }
    1665             : /* x,y QUAD */
    1666             : static GEN
    1667        1316 : mulqq(GEN x, GEN y) {
    1668        1316 :   GEN z = cgetg(4,t_QUAD);
    1669        1316 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1670             :   pari_sp av, tetpil;
    1671        1316 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1672             : 
    1673        1316 :   gel(z,1) = ZX_copy(P); av = avma;
    1674        1316 :   p2 = gmul(gel(x,2),gel(y,2));
    1675        1316 :   p3 = gmul(gel(x,3),gel(y,3));
    1676        1316 :   p1 = gmul(gneg_i(c),p3);
    1677             : 
    1678        1316 :   if (signe(b))
    1679        1204 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1680             :   else
    1681             :   {
    1682         112 :     p3 = gmul(gel(x,2),gel(y,3));
    1683         112 :     p4 = gmul(gel(x,3),gel(y,2));
    1684             :   }
    1685        1316 :   tetpil = avma;
    1686        1316 :   gel(z,2) = gadd(p2,p1);
    1687        1316 :   gel(z,3) = gadd(p4,p3);
    1688        1316 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1689             : }
    1690             : 
    1691             : GEN
    1692     2503828 : mulcxI(GEN x)
    1693             : {
    1694             :   GEN z;
    1695     2503828 :   switch(typ(x))
    1696             :   {
    1697             :     case t_INT: case t_REAL: case t_FRAC:
    1698      308193 :       return mkcomplex(gen_0, x);
    1699             :     case t_COMPLEX:
    1700     2195579 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1701     2193262 :       z = cgetg(3,t_COMPLEX);
    1702     2193262 :       gel(z,1) = gneg(gel(x,2));
    1703     2193262 :       gel(z,2) = gel(x,1); return z;
    1704             :     default:
    1705          56 :       return gmul(gen_I(), x);
    1706             :   }
    1707             : }
    1708             : GEN
    1709       50106 : mulcxmI(GEN x)
    1710             : {
    1711             :   GEN z;
    1712       50106 :   switch(typ(x))
    1713             :   {
    1714             :     case t_INT: case t_REAL: case t_FRAC:
    1715        2100 :       return mkcomplex(gen_0, gneg(x));
    1716             :     case t_COMPLEX:
    1717       47943 :       if (isintzero(gel(x,1))) return gel(x,2);
    1718       46620 :       z = cgetg(3,t_COMPLEX);
    1719       46620 :       gel(z,1) = gel(x,2);
    1720       46620 :       gel(z,2) = gneg(gel(x,1)); return z;
    1721             :     default:
    1722          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1723             :   }
    1724             : }
    1725             : /* x * I^k */
    1726             : GEN
    1727        1491 : mulcxpowIs(GEN x, long k)
    1728             : {
    1729        1491 :   switch (k & 3)
    1730             :   {
    1731         280 :     case 1: return mulcxI(x);
    1732         483 :     case 2: return gneg(x);
    1733         378 :     case 3: return mulcxmI(x);
    1734             :   }
    1735         350 :   return x;
    1736             : }
    1737             : 
    1738             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1739             : static GEN
    1740     2234017 : fill_ser(GEN z, GEN y)
    1741             : {
    1742     2234017 :   long i, lx = lg(z), ly = lg(y);
    1743     2234017 :   if (ly >= lx) {
    1744     2216741 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1745             :   } else {
    1746       17276 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1747       17276 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1748             :   }
    1749     2234017 :   return normalize(z);
    1750             : }
    1751             : 
    1752             : GEN
    1753  3702187410 : gmul(GEN x, GEN y)
    1754             : {
    1755             :   long tx, ty, lx, ly, vx, vy, i, l;
    1756             :   pari_sp av, tetpil;
    1757             :   GEN z, p1, p2;
    1758             : 
    1759  3702187410 :   if (x == y) return gsqr(x);
    1760  2766618739 :   tx = typ(x); ty = typ(y);
    1761  2766618739 :   if (tx == ty) switch(tx)
    1762             :   {
    1763  1242258473 :     case t_INT: return mulii(x,y);
    1764   793874571 :     case t_REAL: return mulrr(x,y);
    1765     1734138 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1766     1734138 :       z = cgetg(3,t_INTMOD);
    1767     1734140 :       if (X==Y || equalii(X,Y))
    1768     1734087 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1769          14 :       gel(z,1) = gcdii(X,Y);
    1770          14 :       warn_coercion(X,Y,gel(z,1));
    1771          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1772          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1773             :     }
    1774             :     case t_FRAC:
    1775             :     {
    1776     9696040 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1777     9696040 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1778     9696040 :       z=cgetg(3,t_FRAC);
    1779     9696040 :       p1 = gcdii(x1, y2);
    1780     9696040 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1781     9696040 :       p1 = gcdii(x2, y1);
    1782     9696040 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1783     9696040 :       tetpil = avma;
    1784     9696040 :       gel(z,2) = mulii(x2,y2);
    1785     9696040 :       gel(z,1) = mulii(x1,y1);
    1786     9696040 :       fix_frac_if_int_GC(z,tetpil); return z;
    1787             :     }
    1788   122545984 :     case t_COMPLEX: return mulcc(x, y);
    1789     1681665 :     case t_PADIC: return mulpp(x, y);
    1790         994 :     case t_QUAD: return mulqq(x, y);
    1791    14566643 :     case t_FFELT: return FF_mul(x, y);
    1792             :     case t_POLMOD:
    1793    10917091 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1794     6615738 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1795     4301353 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1796             :     case t_POL:
    1797    43456748 :       vx = varn(x);
    1798    43456748 :       vy = varn(y);
    1799    43456748 :       if (vx != vy) {
    1800     4782024 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1801     2065507 :         else                     return RgX_Rg_mul(y, x);
    1802             :       }
    1803    38674724 :       return RgX_mul(x, y);
    1804             : 
    1805             :     case t_SER: {
    1806     2237531 :       vx = varn(x);
    1807     2237531 :       vy = varn(y);
    1808     2237531 :       if (vx != vy) {
    1809        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1810        3675 :         else                     return mul_ser_scal(y, x);
    1811             :       }
    1812     2233856 :       lx = lg(x);
    1813     2233856 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1814     2233856 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1815     2233828 :       av = avma; z = cgetg(lx,t_SER);
    1816     2233828 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1817     2233828 :       x = ser2pol_i(x, lx);
    1818     2233828 :       y = ser2pol_i(y, lx);
    1819     2233828 :       y = RgXn_mul(x, y, lx-2);
    1820     2233828 :       return gerepilecopy(av, fill_ser(z,y));
    1821             :     }
    1822    16035591 :     case t_QFI: return qficomp(x,y);
    1823          14 :     case t_QFR: return qfrcomp(x,y);
    1824     6717637 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1825     1780110 :     case t_MAT: return RgM_mul(x, y);
    1826             : 
    1827             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1828        1533 :       z = cgetg_copy(x, &l);
    1829        1533 :       if (l != lg(y)) break;
    1830       16933 :       for (i=1; i<l; i++)
    1831             :       {
    1832       15400 :         long yi = y[i];
    1833       15400 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1834       15400 :         z[i] = x[yi];
    1835             :       }
    1836        1533 :       return z;
    1837             : 
    1838             : 
    1839             :     default:
    1840           0 :       pari_err_TYPE2("*",x,y);
    1841             :   }
    1842             :   /* tx != ty */
    1843   499113976 :   if (is_const_t(ty) && is_const_t(tx))  {
    1844   418900051 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1845   418900051 :     switch(tx) {
    1846             :     case t_INT:
    1847   376059889 :       switch(ty)
    1848             :       {
    1849   242637257 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1850             :         case t_INTMOD:
    1851      971368 :           z = cgetg(3, t_INTMOD);
    1852      971334 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1853             :         case t_FRAC:
    1854    45867819 :           if (!signe(x)) return gen_0;
    1855    28700147 :           z=cgetg(3,t_FRAC);
    1856    28700147 :           p1 = gcdii(x,gel(y,2));
    1857    28700147 :           if (equali1(p1))
    1858             :           {
    1859    16670994 :             set_avma((pari_sp)z);
    1860    16670994 :             gel(z,2) = icopy(gel(y,2));
    1861    16670994 :             gel(z,1) = mulii(gel(y,1), x);
    1862             :           }
    1863             :           else
    1864             :           {
    1865    12029153 :             x = diviiexact(x,p1); tetpil = avma;
    1866    12029153 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1867    12029153 :             gel(z,1) = mulii(gel(y,1), x);
    1868    12029153 :             fix_frac_if_int_GC(z,tetpil);
    1869             :           }
    1870    28700147 :           return z;
    1871    84346402 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1872      412570 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1873        1330 :         case t_QUAD: return mulRq(x,y);
    1874     1823142 :         case t_FFELT: return FF_Z_mul(y,x);
    1875             :       }
    1876             : 
    1877             :     case t_REAL:
    1878    40839760 :       switch(ty)
    1879             :       {
    1880    10057470 :         case t_FRAC: return mulrfrac(x, y);
    1881    30782255 :         case t_COMPLEX: return mulRc(x, y);
    1882           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1883          28 :         default: pari_err_TYPE2("*",x,y);
    1884             :       }
    1885             : 
    1886             :     case t_INTMOD:
    1887        8092 :       switch(ty)
    1888             :       {
    1889        7350 :         case t_FRAC: { GEN X = gel(x,1);
    1890        7350 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1891        7350 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1892             :         }
    1893          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1894         301 :         case t_PADIC: { GEN X = gel(x,1);
    1895         301 :           z = cgetg(3, t_INTMOD);
    1896         301 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1897             :         }
    1898          42 :         case t_QUAD: return mulRq(x, y);
    1899             :         case t_FFELT:
    1900         350 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1901           0 :             pari_err_OP("*",x,y);
    1902         350 :           return FF_Z_mul(y,gel(x,2));
    1903             :       }
    1904             : 
    1905             :     case t_FRAC:
    1906     1992275 :       switch(ty)
    1907             :       {
    1908     1959095 :         case t_COMPLEX: return mulRc(x, y);
    1909       30877 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1910         308 :         case t_QUAD:  return mulRq(x, y);
    1911        1995 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1912             :       }
    1913             : 
    1914             :     case t_FFELT:
    1915          35 :       pari_err_TYPE2("*",x,y);
    1916             : 
    1917             :     case t_COMPLEX:
    1918          21 :       switch(ty)
    1919             :       {
    1920             :         case t_PADIC:
    1921          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1922             :         case t_QUAD:
    1923           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1924           7 :           return mulqf(y, x, lx);
    1925             :       }
    1926             : 
    1927             :     case t_PADIC: /* ty == t_QUAD */
    1928          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1929             :     }
    1930             :   }
    1931             : 
    1932    80214308 :   if (is_matvec_t(ty))
    1933             :   {
    1934     5467600 :     if (!is_matvec_t(tx))
    1935             :     {
    1936     5312434 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1937     5312434 :       z = cgetg_copy(y, &ly);
    1938     5312434 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1939     5312420 :       return z;
    1940             :     }
    1941      155166 :     switch(tx)
    1942             :     {
    1943             :       case t_VEC:
    1944       41708 :         switch(ty) {
    1945       15434 :           case t_COL: return RgV_RgC_mul(x,y);
    1946       26274 :           case t_MAT: return RgV_RgM_mul(x,y);
    1947             :         }
    1948           0 :         break;
    1949             :       case t_COL:
    1950        1687 :         switch(ty) {
    1951        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    1952           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    1953             :         }
    1954           0 :         break;
    1955             :       case t_MAT:
    1956      111771 :         switch(ty) {
    1957           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    1958      111771 :           case t_COL: return RgM_RgC_mul(x,y);
    1959             :         }
    1960             :     }
    1961             :   }
    1962    74746783 :   if (is_matvec_t(tx))
    1963             :   {
    1964      384684 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    1965      384684 :     z = cgetg_copy(x, &lx);
    1966      384685 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    1967      384682 :     return z;
    1968             :   }
    1969    74362099 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1970             :   /* tx < ty, !ismatvec(x and y) */
    1971             : 
    1972    74362099 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    1973     4457210 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    1974    69904889 :   if (is_scalar_t(tx)) {
    1975    66486465 :     if (tx == t_POLMOD) {
    1976     4869373 :       vx = varn(gel(x,1));
    1977     4869373 :       vy = gvar(y);
    1978     4869373 :       if (vx != vy) {
    1979     4434351 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    1980       40061 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    1981             :       }
    1982             :       /* error if ty == t_SER */
    1983      435022 :       av = avma; y = gmod(y, gel(x,1));
    1984      435015 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    1985             :     }
    1986    61617092 :     return mul_scal(y, x, ty);
    1987             :   }
    1988             : 
    1989             :   /* x and y are not scalars, nor matvec */
    1990     3418424 :   vx = gvar(x);
    1991     3418424 :   vy = gvar(y);
    1992     3418424 :   if (vx != vy) /* x or y is treated as a scalar */
    1993     2781531 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    1994     2781531 :                                 : mul_scal(y, x, ty);
    1995             :   /* vx = vy */
    1996     1998311 :   switch(tx)
    1997             :   {
    1998             :     case t_POL:
    1999     1998283 :       switch (ty)
    2000             :       {
    2001             :         case t_SER:
    2002             :         {
    2003             :           long vn;
    2004      133277 :           if (lg(x) == 2) return zeropol(vx);
    2005      133277 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2006      133046 :           av = avma;
    2007      133046 :           vn = RgX_valrem(x, &x);
    2008             :           /* take advantage of x = t^n ! */
    2009      133046 :           if (degpol(x)) {
    2010        1988 :             p1 = RgX_to_ser(x,lg(y));
    2011        1988 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2012        1988 :             p2 = gmul(p1,y);
    2013        1988 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2014             :           } else {
    2015      131058 :             set_avma(av);
    2016      131058 :             p2 = mul_ser_scal(y, gel(x,2));
    2017             :           }
    2018      133046 :           setvalp(p2, valp(p2) + vn);
    2019      133046 :           return p2;
    2020             :         }
    2021             : 
    2022     1865006 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2023             :       }
    2024           0 :       break;
    2025             : 
    2026             :     case t_SER:
    2027          28 :       switch (ty)
    2028             :       {
    2029             :         case t_RFRAC:
    2030          28 :           av = avma;
    2031          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2032             :       }
    2033           0 :       break;
    2034             :   }
    2035           0 :   pari_err_TYPE2("*",x,y);
    2036             :   return NULL; /* LCOV_EXCL_LINE */
    2037             : }
    2038             : 
    2039             : /* return a non-normalized result */
    2040             : GEN
    2041       88833 : sqr_ser_part(GEN x, long l1, long l2)
    2042             : {
    2043             :   long i, j, l;
    2044             :   pari_sp av;
    2045             :   GEN Z, z, p1, p2;
    2046             :   long mi;
    2047       88833 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2048       88826 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2049       88826 :   Z = cgetg(l2-l1+3, t_SER);
    2050       88826 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2051       88826 :   z = Z + 2-l1;
    2052       88826 :   x += 2; mi = 0;
    2053      341138 :   for (i=0; i<l1; i++)
    2054             :   {
    2055      252312 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2056             :   }
    2057             : 
    2058      608562 :   for (i=l1; i<=l2; i++)
    2059             :   {
    2060      519736 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2061      519736 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2062     1064250 :     for (j=i-mi; j<=minss(l,mi); j++)
    2063      544514 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2064      519736 :     p1 = gshift(p1,1);
    2065      519736 :     if ((i&1) == 0 && p2[i>>1])
    2066       50556 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2067      519736 :     gel(z,i) = gerepileupto(av,p1);
    2068             :   }
    2069       88826 :   return Z;
    2070             : }
    2071             : 
    2072             : GEN
    2073  1020956793 : gsqr(GEN x)
    2074             : {
    2075             :   long i, lx;
    2076             :   pari_sp av, tetpil;
    2077             :   GEN z, p1, p2, p3, p4;
    2078             : 
    2079  1020956793 :   switch(typ(x))
    2080             :   {
    2081   945708434 :     case t_INT: return sqri(x);
    2082    29275786 :     case t_REAL: return sqrr(x);
    2083      141849 :     case t_INTMOD: { GEN X = gel(x,1);
    2084      141849 :       z = cgetg(3,t_INTMOD);
    2085      141831 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2086      141813 :       gel(z,1) = icopy(X); return z;
    2087             :     }
    2088     1285617 :     case t_FRAC: return sqrfrac(x);
    2089             : 
    2090             :     case t_COMPLEX:
    2091     3178339 :       if (isintzero(gel(x,1))) {
    2092       36693 :         av = avma;
    2093       36693 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2094             :       }
    2095     3141646 :       z = cgetg(3,t_COMPLEX); av = avma;
    2096     3141646 :       p1 = gadd(gel(x,1),gel(x,2));
    2097     3141646 :       p2 = gsub(gel(x,1), gel(x,2));
    2098     3141646 :       p3 = gmul(gel(x,1),gel(x,2));
    2099     3141646 :       tetpil = avma;
    2100     3141646 :       gel(z,1) = gmul(p1,p2);
    2101     3141646 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2102             : 
    2103             :     case t_PADIC:
    2104       28777 :       z = cgetg(5,t_PADIC);
    2105       28777 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2106       28777 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2107       28777 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2108       28777 :       gel(z,2) = icopy(gel(x,2));
    2109       28777 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2110       28777 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2111       28777 :       return z;
    2112             : 
    2113          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2114          28 :       p1 = gel(x,1);
    2115          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2116          28 :       p2 = gsqr(gel(x,2));
    2117          28 :       p3 = gsqr(gel(x,3));
    2118          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2119             : 
    2120          28 :       if (gequal0(gel(p1,3)))
    2121             :       {
    2122           7 :         tetpil = avma;
    2123           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2124           7 :         av = avma;
    2125           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2126           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2127             :       }
    2128             : 
    2129          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2130          21 :       tetpil = avma;
    2131          21 :       gel(z,2) = gadd(p2,p4);
    2132          21 :       gel(z,3) = gadd(p1,p3);
    2133          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2134             : 
    2135             :     case t_POLMOD:
    2136       43246 :       return sqr_polmod(gel(x,1), gel(x,2));
    2137             : 
    2138     2893766 :     case t_FFELT: return FF_sqr(x);
    2139             : 
    2140      110399 :     case t_POL: return RgX_sqr(x);
    2141             : 
    2142             :     case t_SER:
    2143       24143 :       lx = lg(x);
    2144       24143 :       if (ser_isexactzero(x)) {
    2145           7 :         GEN z = gcopy(x);
    2146           7 :         setvalp(z, 2*valp(x));
    2147           7 :         return z;
    2148             :       }
    2149       24136 :       if (lx < 40)
    2150       23947 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2151             :       else
    2152             :       {
    2153         189 :         pari_sp av = avma;
    2154         189 :         GEN z = cgetg(lx,t_SER);
    2155         189 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2156         189 :         x = ser2pol_i(x,lx);
    2157         189 :         x = RgXn_sqr(x, lx-2);
    2158         189 :         return gerepilecopy(av, fill_ser(z,x));
    2159             :       }
    2160             : 
    2161           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2162           7 :       gel(z,1) = gsqr(gel(x,1));
    2163           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2164             : 
    2165         980 :     case t_MAT: return RgM_sqr(x);
    2166          14 :     case t_QFR: return qfrsqr(x);
    2167    38264792 :     case t_QFI: return qfisqr(x);
    2168             :     case t_VECSMALL:
    2169         616 :       z = cgetg_copy(x, &lx);
    2170        6720 :       for (i=1; i<lx; i++)
    2171             :       {
    2172        6104 :         long xi = x[i];
    2173        6104 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2174        6104 :         z[i] = x[xi];
    2175             :       }
    2176         616 :       return z;
    2177             :   }
    2178           0 :   pari_err_TYPE2("*",x,x);
    2179             :   return NULL; /* LCOV_EXCL_LINE */
    2180             : }
    2181             : 
    2182             : /********************************************************************/
    2183             : /**                                                                **/
    2184             : /**                           DIVISION                             **/
    2185             : /**                                                                **/
    2186             : /********************************************************************/
    2187             : static GEN
    2188      271036 : div_rfrac_scal(GEN x, GEN y)
    2189             : {
    2190      271036 :   pari_sp av = avma;
    2191      271036 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2192      271036 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2193             : }
    2194             : static GEN
    2195       37452 : div_scal_rfrac(GEN x, GEN y)
    2196             : {
    2197       37452 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2198       37452 :   pari_sp av = avma;
    2199       37452 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2200             :   {
    2201           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2202           0 :     y1 = gel(y1,2);
    2203             :   }
    2204       37452 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2205             : }
    2206             : static GEN
    2207     1186544 : div_rfrac(GEN x, GEN y)
    2208     1186544 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2209             : 
    2210             : static GEN
    2211      622417 : div_ser_scal(GEN x, GEN y) {
    2212             :   long i, lx;
    2213             :   GEN z;
    2214      622417 :   if (ser_isexactzero(x))
    2215             :   {
    2216           7 :     if (lg(x) == 2) return gcopy(x);
    2217           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2218             :   }
    2219      622410 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2220      622410 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2221      622410 :   return normalize(z);
    2222             : }
    2223             : GEN
    2224        7525 : ser_normalize(GEN x)
    2225             : {
    2226        7525 :   long i, lx = lg(x);
    2227             :   GEN c, z;
    2228        7525 :   if (lx == 2) return x;
    2229        7525 :   c = gel(x,2); if (gequal1(c)) return x;
    2230        7448 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2231        7448 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2232        7448 :   return z;
    2233             : }
    2234             : 
    2235             : static GEN
    2236     1644993 : div_T_scal(GEN x, GEN y, long tx) {
    2237     1644993 :   switch(tx)
    2238             :   {
    2239      755404 :     case t_POL: return RgX_Rg_div(x, y);
    2240      622410 :     case t_SER: return div_ser_scal(x, y);
    2241      267179 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2242             :   }
    2243           0 :   pari_err_TYPE2("/",x,y);
    2244             :   return NULL; /* LCOV_EXCL_LINE */
    2245             : }
    2246             : 
    2247             : static GEN
    2248     9217101 : div_scal_pol(GEN x, GEN y) {
    2249     9217101 :   long ly = lg(y);
    2250             :   pari_sp av;
    2251     9217101 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2252     9139155 :   if (isrationalzero(x)) return zeropol(varn(y));
    2253     7091735 :   av = avma;
    2254     7091735 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2255             : }
    2256             : static GEN
    2257       32907 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2258             :   GEN z;
    2259             :   long ly, i;
    2260       32907 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2261       32907 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2262       32907 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2263       32907 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2264       32907 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2265       32907 :   y = gdiv(z,y); pari_free(z); return y;
    2266             : }
    2267             : static GEN
    2268     9240357 : div_scal_T(GEN x, GEN y, long ty) {
    2269     9240357 :   switch(ty)
    2270             :   {
    2271     9170754 :     case t_POL: return div_scal_pol(x, y);
    2272       32907 :     case t_SER: return div_scal_ser(x, y);
    2273       36696 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2274             :   }
    2275           0 :   pari_err_TYPE2("/",x,y);
    2276             :   return NULL; /* LCOV_EXCL_LINE */
    2277             : }
    2278             : 
    2279             : /* assume tx = ty = t_SER, same variable vx */
    2280             : static GEN
    2281      877704 : div_ser(GEN x, GEN y, long vx)
    2282             : {
    2283      877704 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2284             :   GEN y_lead, p1, p2, z;
    2285             : 
    2286      877704 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2287      877697 :   if (ser_isexactzero(x))
    2288             :   {
    2289       59843 :     if (lx == 2) return zeroser(vx, l);
    2290         154 :     return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), l);
    2291             :   }
    2292      817854 :   y_lead = gel(y,2);
    2293      817854 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2294             :   {
    2295          14 :     GEN y0 = y;
    2296          14 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2297          35 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2298             :     {
    2299          28 :       y_lead = gel(y,2);
    2300          28 :       if (!gequal0(y_lead)) break;
    2301             :     }
    2302          14 :     if (ly <= 2) pari_err_INV("div_ser", y0);
    2303             :   }
    2304      817847 :   if (ly < lx) lx = ly;
    2305      817847 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2306     3022314 :   for (i=3; i<lx; i++)
    2307             :   {
    2308     2204467 :     p1 = gel(y,i);
    2309     2204467 :     if (isrationalzero(p1)) p1 = NULL;
    2310     2204467 :     gel(p2,i) = p1;
    2311             :   }
    2312      817847 :   z = cgetg(lx,t_SER);
    2313      817847 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2314      817847 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2315     3022314 :   for (i=3; i<lx; i++)
    2316             :   {
    2317     2204467 :     pari_sp av = avma;
    2318     2204467 :     p1 = gel(x,i);
    2319    18358379 :     for (j=2, l=i; j<i; j++, l--)
    2320    16153912 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2321     2204467 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2322             :   }
    2323      817847 :   return normalize(z);
    2324             : }
    2325             : /* x,y compatible PADIC */
    2326             : static GEN
    2327      225998 : divpp(GEN x, GEN y) {
    2328             :   pari_sp av;
    2329             :   long a, b;
    2330             :   GEN z, M;
    2331             : 
    2332      225998 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2333      225996 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2334      225338 :   a = precp(x);
    2335      225338 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2336      225338 :   z = cgetg(5, t_PADIC);
    2337      225365 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2338      225365 :   gel(z,2) = icopy(gel(x,2));
    2339      225372 :   gel(z,3) = icopy(M); av = avma;
    2340      225379 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2341      225374 :   return z;
    2342             : }
    2343             : static GEN
    2344       45534 : div_polmod_same(GEN T, GEN x, GEN y)
    2345             : {
    2346       45534 :   long v = varn(T);
    2347       45534 :   GEN a, z = cgetg(3, t_POLMOD);
    2348       45534 :   gel(z,1) = RgX_copy(T);
    2349       45534 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2350       26775 :     a = gdiv(x, y);
    2351       18759 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2352         273 :   {
    2353         273 :     pari_sp av = avma;
    2354         273 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2355             :   }
    2356       18486 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2357        6118 :   {
    2358        6118 :     pari_sp av = avma;
    2359        6118 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2360        6118 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2361        6118 :     a = gerepileupto(av, a);
    2362             :   }
    2363             :   else
    2364             :   {
    2365       12368 :     pari_sp av = avma;
    2366       12368 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2367       12368 :     a = gerepileupto(av, a);
    2368             :   }
    2369       45534 :   gel(z,2) = a; return z;
    2370             : }
    2371             : GEN
    2372   173359798 : gdiv(GEN x, GEN y)
    2373             : {
    2374   173359798 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2375             :   pari_sp av, tetpil;
    2376             :   GEN z, p1, p2;
    2377             : 
    2378   173359798 :   if (tx == ty) switch(tx)
    2379             :   {
    2380             :     case t_INT:
    2381    67189159 :       return Qdivii(x,y);
    2382             : 
    2383    50016226 :     case t_REAL: return divrr(x,y);
    2384       38030 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2385       38030 :       z = cgetg(3,t_INTMOD);
    2386       38029 :       if (X==Y || equalii(X,Y))
    2387       37997 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2388          28 :       gel(z,1) = gcdii(X,Y);
    2389          28 :       warn_coercion(X,Y,gel(z,1));
    2390          28 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2391          28 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2392             :     }
    2393             :     case t_FRAC: {
    2394     1671494 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2395     1671494 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2396     1671494 :       z = cgetg(3, t_FRAC);
    2397     1671494 :       p1 = gcdii(x1, y1);
    2398     1671494 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2399     1671494 :       p1 = gcdii(x2, y2);
    2400     1671494 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2401     1671494 :       tetpil = avma;
    2402     1671494 :       gel(z,2) = mulii(x2,y1);
    2403     1671494 :       gel(z,1) = mulii(x1,y2);
    2404     1671494 :       normalize_frac(z);
    2405     1671494 :       fix_frac_if_int_GC(z,tetpil);
    2406     1671494 :       return z;
    2407             :     }
    2408             :     case t_COMPLEX:
    2409     2780373 :       if (isintzero(gel(y,1)))
    2410             :       {
    2411       58163 :         y = gel(y,2);
    2412       58163 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2413       40964 :         z = cgetg(3,t_COMPLEX);
    2414       40964 :         gel(z,1) = gdiv(gel(x,2), y);
    2415       40964 :         av = avma;
    2416       40964 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2417       40964 :         return z;
    2418             :       }
    2419     2722210 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2420     2722210 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2421             : 
    2422             :     case t_PADIC:
    2423      123180 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2424      123180 :       return divpp(x, y);
    2425             : 
    2426             :     case t_QUAD:
    2427         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2428         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2429         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2430             : 
    2431      236797 :     case t_FFELT: return FF_div(x,y);
    2432             : 
    2433             :     case t_POLMOD:
    2434       89109 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2435       45534 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2436             :       else {
    2437       43575 :         av = avma;
    2438       43575 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2439             :       }
    2440       89109 :       return z;
    2441             : 
    2442             :     case t_POL:
    2443    18503845 :       vx = varn(x);
    2444    18503845 :       vy = varn(y);
    2445    18503845 :       if (vx != vy) {
    2446       95088 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2447       46347 :                             else return div_scal_pol(x, y);
    2448             :       }
    2449    18408757 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2450    18408757 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2451    18221658 :       av = avma;
    2452    18221658 :       return gerepileupto(av, gred_rfrac2(x,y));
    2453             : 
    2454             :     case t_SER:
    2455      431979 :       vx = varn(x);
    2456      431979 :       vy = varn(y);
    2457      431979 :       if (vx != vy) {
    2458           7 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2459           0 :                             else return div_scal_ser(x, y);
    2460             :       }
    2461      431972 :       return div_ser(x, y, vx);
    2462             :     case t_RFRAC:
    2463     1191024 :       vx = varn(gel(x,2));
    2464     1191024 :       vy = varn(gel(y,2));
    2465     1191024 :       if (vx != vy) {
    2466        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2467         756 :                             else return div_scal_rfrac(x, y);
    2468             :       }
    2469     1186544 :       return div_rfrac(x,y);
    2470             : 
    2471           7 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2472           7 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2473             : 
    2474             :     case t_MAT:
    2475           7 :       av = avma; p1 = RgM_inv(y);
    2476           7 :       if (!p1) pari_err_INV("gdiv",y);
    2477           7 :       return gerepileupto(av, RgM_mul(x, p1));
    2478             : 
    2479           0 :     default: pari_err_TYPE2("/",x,y);
    2480             :   }
    2481             : 
    2482    31088239 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2483             :   {
    2484     2716461 :     long s = signe(x);
    2485     2716461 :     if (!s) {
    2486      698309 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2487      698309 :       switch (ty)
    2488             :       {
    2489      695117 :       default: return gen_0;
    2490             :       case t_INTMOD:
    2491          56 :         z = cgetg(3,t_INTMOD);
    2492          56 :         gel(z,1) = icopy(gel(y,1));
    2493          56 :         gel(z,2) = gen_0; return z;
    2494        3136 :       case t_FFELT: return FF_zero(y);
    2495             :       }
    2496             :     }
    2497     2018152 :     if (is_pm1(x)) {
    2498      557744 :       if (s > 0) return ginv(y);
    2499       43369 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2500             :     }
    2501     1460408 :     switch(ty)
    2502             :     {
    2503      154263 :       case t_REAL: return divir(x,y);
    2504             :       case t_INTMOD:
    2505          21 :         z = cgetg(3, t_INTMOD);
    2506          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2507             :       case t_FRAC:
    2508     1188238 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2509     1188238 :         if (equali1(p1))
    2510             :         {
    2511      811415 :           set_avma((pari_sp)z);
    2512      811415 :           gel(z,2) = icopy(gel(y,1));
    2513      811415 :           gel(z,1) = mulii(gel(y,2), x);
    2514      811415 :           normalize_frac(z);
    2515      811415 :           fix_frac_if_int(z);
    2516             :         }
    2517             :         else
    2518             :         {
    2519      376823 :           x = diviiexact(x,p1); tetpil = avma;
    2520      376823 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2521      376823 :           gel(z,1) = mulii(gel(y,2), x);
    2522      376823 :           normalize_frac(z);
    2523      376823 :           fix_frac_if_int_GC(z,tetpil);
    2524             :         }
    2525     1188238 :         return z;
    2526             : 
    2527         266 :       case t_FFELT: return Z_FF_div(x,y);
    2528      105377 :       case t_COMPLEX: return divRc(x,y);
    2529       12201 :       case t_PADIC: return divTp(x, y);
    2530             :       case t_QUAD:
    2531          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2532          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2533             :     }
    2534             :   }
    2535    28371778 :   if (gequal0(y))
    2536             :   {
    2537          42 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2538          21 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2539             :   }
    2540             : 
    2541    28371750 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2542             :   {
    2543             :     case t_REAL:
    2544     9147110 :       switch(ty)
    2545             :       {
    2546     7935859 :         case t_INT: return divri(x,y);
    2547             :         case t_FRAC:
    2548      734359 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2549      734359 :           return gerepileuptoleaf(av, z);
    2550      476864 :         case t_COMPLEX: return divRc(x, y);
    2551           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2552          21 :         default: pari_err_TYPE2("/",x,y);
    2553             :       }
    2554             : 
    2555             :     case t_INTMOD:
    2556       13482 :       switch(ty)
    2557             :       {
    2558             :         case t_INT:
    2559       13405 :           z = cgetg(3, t_INTMOD);
    2560       13405 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2561          28 :         case t_FRAC: { GEN X = gel(x,1);
    2562          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2563          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2564             :         }
    2565             :         case t_FFELT:
    2566           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2567           0 :             pari_err_OP("/",x,y);
    2568           7 :           return Z_FF_div(gel(x,2),y);
    2569             : 
    2570             :         case t_COMPLEX:
    2571           0 :           av = avma;
    2572           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2573             : 
    2574             :         case t_QUAD:
    2575           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2576           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2577             : 
    2578           7 :         case t_PADIC: { GEN X = gel(x,1);
    2579           7 :           z = cgetg(3, t_INTMOD);
    2580           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2581             :         }
    2582          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2583             :       }
    2584             : 
    2585             :     case t_FRAC:
    2586      685858 :       switch(ty)
    2587             :       {
    2588      643080 :         case t_INT: z = cgetg(3, t_FRAC);
    2589      643080 :         p1 = gcdii(y,gel(x,1));
    2590      643080 :         if (equali1(p1))
    2591             :         {
    2592      355827 :           set_avma((pari_sp)z); tetpil = 0;
    2593      355827 :           gel(z,1) = icopy(gel(x,1));
    2594             :         }
    2595             :         else
    2596             :         {
    2597      287253 :           y = diviiexact(y,p1); tetpil = avma;
    2598      287253 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2599             :         }
    2600      643080 :         gel(z,2) = mulii(gel(x,2),y);
    2601      643080 :         normalize_frac(z);
    2602      643080 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2603      643080 :         return z;
    2604             : 
    2605             :         case t_REAL:
    2606       40398 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2607       40398 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2608             : 
    2609           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2610           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2611           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2612             :         }
    2613             : 
    2614          28 :         case t_FFELT: av=avma;
    2615          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2616             : 
    2617         427 :         case t_COMPLEX: return divRc(x, y);
    2618             : 
    2619             :         case t_PADIC:
    2620        1911 :           if (!signe(gel(x,1))) return gen_0;
    2621        1911 :           return divTp(x, y);
    2622             : 
    2623             :         case t_QUAD:
    2624           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2625           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2626             :       }
    2627             : 
    2628             :     case t_FFELT:
    2629        1533 :       switch (ty)
    2630             :       {
    2631        1456 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2632          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2633             :         case t_INTMOD:
    2634           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2635           0 :             pari_err_OP("/",x,y);
    2636           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2637             :         default:
    2638          42 :         pari_err_TYPE2("/",x,y);
    2639             :       }
    2640           0 :       break;
    2641             : 
    2642             :     case t_COMPLEX:
    2643     4591882 :       switch(ty)
    2644             :       {
    2645     4591882 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2646           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2647             :         case t_PADIC:
    2648           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2649             :         case t_QUAD:
    2650           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2651           0 :           return divfq(x, y, lx);
    2652             :       }
    2653             : 
    2654             :     case t_PADIC:
    2655        1988 :       switch(ty)
    2656             :       {
    2657        1939 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2658        1939 :           return signe(gel(x,4))? divpT(x, y)
    2659        1939 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2660             :         }
    2661           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2662           7 :           z = cgetg(3, t_INTMOD);
    2663           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2664             :         }
    2665             :         case t_COMPLEX: case t_QUAD:
    2666           7 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2667           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2668             : 
    2669          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2670             :       }
    2671             : 
    2672             :     case t_QUAD:
    2673         896 :       switch (ty)
    2674             :       {
    2675             :         case t_INT: case t_INTMOD: case t_FRAC:
    2676         847 :           z = cgetg(4,t_QUAD);
    2677         847 :           gel(z,1) = ZX_copy(gel(x,1));
    2678         847 :           gel(z,2) = gdiv(gel(x,2), y);
    2679         847 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2680          28 :         case t_REAL: return divqf(x, y, lg(y));
    2681           7 :         case t_PADIC: return divTp(x, y);
    2682             :         case t_COMPLEX:
    2683           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2684           0 :           return divqf(x, y, ly);
    2685             :       }
    2686             :   }
    2687    13929022 :   switch(ty) {
    2688             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2689       33592 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2690             :     case t_MAT:
    2691          42 :       av = avma; p1 = RgM_inv(y);
    2692          35 :       if (!p1) pari_err_INV("gdiv",y);
    2693          28 :       return gerepileupto(av, gmul(x, p1));
    2694             :     case t_VEC: case t_COL:
    2695             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2696           0 :       pari_err_TYPE2("/",x,y);
    2697             :   }
    2698    13895388 :   switch(tx) {
    2699             :     case t_VEC: case t_COL: case t_MAT:
    2700     1033226 :       z = cgetg_copy(x, &lx);
    2701     1033226 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2702     1033226 :       return z;
    2703             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2704           0 :       pari_err_TYPE2("/",x,y);
    2705             :   }
    2706             : 
    2707    12862162 :   vy = gvar(y);
    2708    12862162 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2709      215369 :     vx = varn(X);
    2710      215369 :     if (vx != vy) {
    2711      214970 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2712      214641 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2713             :     }
    2714             :     /* y is POL, SER or RFRAC */
    2715         399 :     av = avma;
    2716         399 :     switch(ty)
    2717             :     {
    2718           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2719         399 :       default: y = ginvmod(gmod(y,X), X);
    2720             :     }
    2721         392 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2722             :   }
    2723             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2724             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2725             :    * variable NO_VARIABLE, then the operation is incorrect */
    2726    12646793 :   vx = gvar(x);
    2727    12646793 :   if (vx != vy) { /* includes cases where one is scalar */
    2728    10885021 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2729     9240028 :                         else return div_scal_T(x, y, ty);
    2730             :   }
    2731     1761772 :   switch(tx)
    2732             :   {
    2733             :     case t_POL:
    2734     1046122 :       switch(ty)
    2735             :       {
    2736             :         case t_SER:
    2737          70 :           if (lg(y) == 2)
    2738           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2739          70 :           p1 = RgX_to_ser(x,lg(y));
    2740          70 :           p2 = div_ser(p1, y, vx);
    2741          63 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2742          63 :           return p2;
    2743             : 
    2744             :         case t_RFRAC:
    2745             :         {
    2746     1046052 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2747     1046052 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2748           0 :             return mul_rfrac_scal(y2, y1, x);
    2749     1046052 :           av = avma;
    2750     1046052 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2751             :         }
    2752             :       }
    2753           0 :       break;
    2754             : 
    2755             :     case t_SER:
    2756      445781 :       switch(ty)
    2757             :       {
    2758             :         case t_POL:
    2759      445767 :           if (lg(x) == 2)
    2760         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2761      445662 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2762      445662 :           p2 = div_ser(x, p1, vx);
    2763      445662 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2764      445662 :           return p2;
    2765             :         case t_RFRAC:
    2766          14 :           av = avma;
    2767          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2768             :       }
    2769           0 :       break;
    2770             : 
    2771             :     case t_RFRAC:
    2772      269855 :       switch(ty)
    2773             :       {
    2774      269855 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2775             :         case t_SER:
    2776           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2777           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2778             :       }
    2779           0 :       break;
    2780             :   }
    2781          14 :   pari_err_TYPE2("/",x,y);
    2782             :   return NULL; /* LCOV_EXCL_LINE */
    2783             : }
    2784             : 
    2785             : /********************************************************************/
    2786             : /**                                                                **/
    2787             : /**                     SIMPLE MULTIPLICATION                      **/
    2788             : /**                                                                **/
    2789             : /********************************************************************/
    2790             : GEN
    2791   125360720 : gmulsg(long s, GEN y)
    2792             : {
    2793             :   long ly, i;
    2794             :   pari_sp av;
    2795             :   GEN z;
    2796             : 
    2797   125360720 :   switch(typ(y))
    2798             :   {
    2799    88981359 :     case t_INT:  return mulsi(s,y);
    2800    27321036 :     case t_REAL: return mulsr(s,y);
    2801      172253 :     case t_INTMOD: { GEN p = gel(y,1);
    2802      172253 :       z = cgetg(3,t_INTMOD);
    2803      172247 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2804      172236 :       gel(z,1) = icopy(p); return z;
    2805             :     }
    2806      627261 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2807             :     case t_FRAC:
    2808      844534 :       if (!s) return gen_0;
    2809      840131 :       z = cgetg(3,t_FRAC);
    2810      840131 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2811      840131 :       if (i == 1)
    2812             :       {
    2813      388828 :         gel(z,2) = icopy(gel(y,2));
    2814      388828 :         gel(z,1) = mulis(gel(y,1), s);
    2815             :       }
    2816             :       else
    2817             :       {
    2818      451303 :         gel(z,2) = divis(gel(y,2), i);
    2819      451303 :         gel(z,1) = mulis(gel(y,1), s/i);
    2820      451303 :         fix_frac_if_int(z);
    2821             :       }
    2822      840131 :       return z;
    2823             : 
    2824     6137179 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2825     6137179 :       gel(z,1) = gmulsg(s,gel(y,1));
    2826     6137179 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2827             : 
    2828             :     case t_PADIC:
    2829        9212 :       if (!s) return gen_0;
    2830        9212 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2831             : 
    2832           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2833           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2834           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2835           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2836             : 
    2837             :     case t_POLMOD:
    2838       85596 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2839             : 
    2840             :     case t_POL:
    2841      590050 :       if (!signe(y)) return RgX_copy(y);
    2842      574566 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2843      573208 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2844      573213 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2845      573182 :       return normalizepol_lg(z, ly);
    2846             : 
    2847             :     case t_SER:
    2848         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2849         182 :       if (!s) return Rg_get_0(y);
    2850         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2851         182 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2852         182 :       return normalize(z);
    2853             : 
    2854             :     case t_RFRAC:
    2855          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2856          14 :       if (s == 1) return gcopy(y);
    2857          14 :       if (s == -1) return gneg(y);
    2858          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2859             : 
    2860             :     case t_VEC: case t_COL: case t_MAT:
    2861      592037 :       z = cgetg_copy(y, &ly);
    2862      592037 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2863      592037 :       return z;
    2864             :   }
    2865           0 :   pari_err_TYPE("gmulsg",y);
    2866             :   return NULL; /* LCOV_EXCL_LINE */
    2867             : }
    2868             : 
    2869             : /********************************************************************/
    2870             : /**                                                                **/
    2871             : /**                       SIMPLE DIVISION                          **/
    2872             : /**                                                                **/
    2873             : /********************************************************************/
    2874             : 
    2875             : GEN
    2876    18427548 : gdivgs(GEN x, long s)
    2877             : {
    2878    18427548 :   long tx = typ(x), lx, i;
    2879             :   pari_sp av;
    2880             :   GEN z, y;
    2881             : 
    2882    18427548 :   if (!s)
    2883             :   {
    2884           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2885           0 :     pari_err_INV("gdivgs",gen_0);
    2886             :   }
    2887    18427548 :   switch(tx)
    2888             :   {
    2889             :     case t_INT:
    2890     7069615 :       av = avma; z = divis_rem(x,s,&i);
    2891     7069615 :       if (!i) return z;
    2892             : 
    2893     4860198 :       i = cgcd(s, i);
    2894     4860198 :       set_avma(av); z = cgetg(3,t_FRAC);
    2895     4860198 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    2896     4860198 :       gel(z,1) = y;
    2897     4860198 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    2898             : 
    2899             :     case t_REAL:
    2900     7833768 :       return divrs(x,s);
    2901             : 
    2902             :     case t_INTMOD:
    2903       14651 :       z = cgetg(3, t_INTMOD);
    2904       14651 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    2905             : 
    2906        1043 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    2907             : 
    2908      890243 :     case t_FRAC: z = cgetg(3, t_FRAC);
    2909      890243 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    2910      890243 :       if (i == 1)
    2911             :       {
    2912      408279 :         gel(z,2) = mulsi(s, gel(x,2));
    2913      408279 :         gel(z,1) = icopy(gel(x,1));
    2914             :       }
    2915             :       else
    2916             :       {
    2917      481964 :         gel(z,2) = mulsi(s/i, gel(x,2));
    2918      481964 :         gel(z,1) = divis(gel(x,1), i);
    2919             :       }
    2920      890243 :       normalize_frac(z);
    2921      890243 :       fix_frac_if_int(z); return z;
    2922             : 
    2923     2434953 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2924     2434953 :       gel(z,1) = gdivgs(gel(x,1),s);
    2925     2434953 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    2926             : 
    2927             :     case t_PADIC: /* divpT */
    2928             :     {
    2929       87648 :       GEN p = gel(x,2);
    2930       87648 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    2931       86927 :       av = avma;
    2932       86927 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    2933             :     }
    2934             : 
    2935           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2936           0 :       gel(z,1) = ZX_copy(gel(x,1));
    2937           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    2938           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    2939             : 
    2940             :     case t_POLMOD:
    2941       24822 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    2942             : 
    2943             :     case t_RFRAC:
    2944         140 :       if (s == 1) return gcopy(x);
    2945         133 :       else if (s == -1) return gneg(x);
    2946         133 :       return div_rfrac_scal(x, stoi(s));
    2947             : 
    2948             :     case t_POL: case t_SER:
    2949       70364 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    2950       70364 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2951       70364 :       return z;
    2952             :     case t_VEC: case t_COL: case t_MAT:
    2953         301 :       z = cgetg_copy(x, &lx);
    2954         301 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2955         301 :       return z;
    2956             : 
    2957             :   }
    2958           0 :   pari_err_TYPE2("/",x, stoi(s));
    2959             :   return NULL; /* LCOV_EXCL_LINE */
    2960             : }
    2961             : 
    2962             : /* True shift (exact multiplication by 2^n) */
    2963             : GEN
    2964    55165141 : gmul2n(GEN x, long n)
    2965             : {
    2966             :   long lx, i, k, l;
    2967             :   GEN z, a, b;
    2968             : 
    2969    55165141 :   switch(typ(x))
    2970             :   {
    2971             :     case t_INT:
    2972    39684653 :       if (n>=0) return shifti(x,n);
    2973     2399153 :       if (!signe(x)) return gen_0;
    2974     2341409 :       l = vali(x); n = -n;
    2975     2341409 :       if (n<=l) return shifti(x,-n);
    2976      304306 :       z = cgetg(3,t_FRAC);
    2977      304306 :       gel(z,1) = shifti(x,-l);
    2978      304306 :       gel(z,2) = int2n(n-l); return z;
    2979             : 
    2980             :     case t_REAL:
    2981     7396317 :       return shiftr(x,n);
    2982             : 
    2983       75296 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    2984       75296 :       z = cgetg(3,t_INTMOD);
    2985       75288 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    2986       75022 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    2987       75015 :       gel(z,1) = icopy(b); return z;
    2988             : 
    2989      216867 :     case t_FFELT: return FF_mul2n(x,n);
    2990             : 
    2991      528268 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    2992      528268 :       l = vali(a);
    2993      528268 :       k = vali(b);
    2994      528268 :       if (n+l >= k)
    2995             :       {
    2996      141736 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    2997      101969 :         l = n-k; k = -k;
    2998             :       }
    2999             :       else
    3000             :       {
    3001      386532 :         k = -(l+n); l = -l;
    3002             :       }
    3003      488501 :       z = cgetg(3,t_FRAC);
    3004      488501 :       gel(z,1) = shifti(a,l);
    3005      488501 :       gel(z,2) = shifti(b,k); return z;
    3006             : 
    3007     5670963 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3008     5670963 :       gel(z,1) = gmul2n(gel(x,1),n);
    3009     5670963 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3010             : 
    3011          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3012          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3013          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3014          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3015             : 
    3016             :     case t_POLMOD:
    3017       18508 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3018             : 
    3019             :     case t_POL:
    3020      492194 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3021      492194 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3022      492194 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3023             :     case t_SER:
    3024      122153 :       if (ser_isexactzero(x)) return gcopy(x);
    3025      122139 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3026      122139 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3027      122139 :       return normalize(z); /* needed if char = 2 */
    3028             :     case t_VEC: case t_COL: case t_MAT:
    3029      941845 :       z = cgetg_copy(x, &lx);
    3030      941845 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3031      941845 :       return z;
    3032             : 
    3033             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3034          21 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3035             : 
    3036             :     case t_PADIC: /* int2n wrong if n < 0 */
    3037       18042 :       return gmul(gmul2n(gen_1,n),x);
    3038             :   }
    3039           0 :   pari_err_TYPE("gmul2n",x);
    3040             :   return NULL; /* LCOV_EXCL_LINE */
    3041             : }
    3042             : 
    3043             : /*******************************************************************/
    3044             : /*                                                                 */
    3045             : /*                              INVERSE                            */
    3046             : /*                                                                 */
    3047             : /*******************************************************************/
    3048             : static GEN
    3049      133800 : inv_polmod(GEN T, GEN x)
    3050             : {
    3051      133800 :   GEN z = cgetg(3,t_POLMOD), a;
    3052      133800 :   gel(z,1) = RgX_copy(T);
    3053      133800 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3054      105751 :     a = ginv(x);
    3055             :   else
    3056             :   {
    3057       28049 :     if (lg(T) == 5) /* quadratic fields */
    3058       22582 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3059             :     else
    3060        5467 :       a = RgXQ_inv(x, T);
    3061             :   }
    3062      133798 :   gel(z,2) = a; return z;
    3063             : }
    3064             : 
    3065             : GEN
    3066    18173915 : ginv(GEN x)
    3067             : {
    3068             :   long s;
    3069             :   pari_sp av, tetpil;
    3070             :   GEN z, y, p1, p2;
    3071             : 
    3072    18173915 :   switch(typ(x))
    3073             :   {
    3074             :     case t_INT:
    3075     1955026 :       if (is_pm1(x)) return icopy(x);
    3076     1587350 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3077     1587350 :       z = cgetg(3,t_FRAC);
    3078     1587350 :       gel(z,1) = s<0? gen_m1: gen_1;
    3079     1587350 :       gel(z,2) = absi(x); return z;
    3080             : 
    3081     3214221 :     case t_REAL: return invr(x);
    3082             : 
    3083        4641 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3084        4641 :       gel(z,1) = icopy(gel(x,1));
    3085        4641 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3086             : 
    3087             :     case t_FRAC: {
    3088      251112 :       GEN a = gel(x,1), b = gel(x,2);
    3089      251112 :       s = signe(a);
    3090      251112 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3091      111956 :       z = cgetg(3,t_FRAC);
    3092      111956 :       gel(z,1) = icopy(b);
    3093      111956 :       gel(z,2) = icopy(a);
    3094      111956 :       normalize_frac(z); return z;
    3095             :     }
    3096             :     case t_COMPLEX:
    3097     2356971 :       av=avma;
    3098     2356971 :       p1=cxnorm(x);
    3099     2356971 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3100     2356971 :       tetpil=avma;
    3101     2356971 :       return gerepile(av,tetpil,divcR(p2,p1));
    3102             : 
    3103             :     case t_QUAD:
    3104         203 :       av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
    3105         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3106             : 
    3107        5460 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3108        5460 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3109        5453 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3110        5453 :       gel(z,2) = icopy(gel(x,2));
    3111        5453 :       gel(z,3) = icopy(gel(x,3));
    3112        5453 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3113             : 
    3114      133800 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3115       14085 :     case t_FFELT: return FF_inv(x);
    3116     5430983 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3117       15421 :     case t_SER: return gdiv(gen_1,x);
    3118             : 
    3119             :     case t_RFRAC:
    3120             :     {
    3121        2912 :       GEN n = gel(x,1), d = gel(x,2);
    3122        2912 :       pari_sp av = avma, ltop;
    3123        2912 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3124             : 
    3125        2912 :       n = simplify_shallow(n);
    3126        2912 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3127             :       {
    3128        2912 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3129         679 :         ltop = avma;
    3130         679 :         z = RgX_Rg_div(d,n);
    3131             :       } else {
    3132           0 :         ltop = avma;
    3133           0 :         z = cgetg(3,t_RFRAC);
    3134           0 :         gel(z,1) = RgX_copy(d);
    3135           0 :         gel(z,2) = RgX_copy(n);
    3136             :       }
    3137         679 :       stackdummy(av, ltop);
    3138         679 :       return z;
    3139             :     }
    3140             : 
    3141             :     case t_QFR:
    3142          14 :       av = avma; z = cgetg(5, t_QFR);
    3143          14 :       gel(z,1) = gel(x,1);
    3144          14 :       gel(z,2) = negi( gel(x,2) );
    3145          14 :       gel(z,3) = gel(x,3);
    3146          14 :       gel(z,4) = negr( gel(x,4) );
    3147          14 :       return gerepileupto(av, redreal(z));
    3148             : 
    3149             :     case t_QFI:
    3150     4784971 :       y = gcopy(x);
    3151     4784971 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3152     4670039 :         togglesign(gel(y,2));
    3153     4784971 :       return y;
    3154             :     case t_MAT:
    3155        4095 :       y = RgM_inv(x);
    3156        4088 :       if (!y) pari_err_INV("ginv",x);
    3157        4018 :       return y;
    3158             :     case t_VECSMALL:
    3159             :     {
    3160           0 :       long i, lx = lg(x)-1;
    3161           0 :       y = zero_zv(lx);
    3162           0 :       for (i=1; i<=lx; i++)
    3163             :       {
    3164           0 :         long xi = x[i];
    3165           0 :         if (xi<1 || xi>lx || y[xi])
    3166           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3167           0 :         y[xi] = i;
    3168             :       }
    3169           0 :       return y;
    3170             :     }
    3171             :   }
    3172           0 :   pari_err_TYPE("inverse",x);
    3173             :   return NULL; /* LCOV_EXCL_LINE */
    3174             : }

Generated by: LCOV version 1.13