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.1 lcov report (development 25819-e703fe1174) Lines: 1837 1936 94.9 %
Date: 2020-09-18 06:10:04 Functions: 94 94 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          28 : kro_quad(GEN x, GEN y)
      42          28 : { 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   105741697 : addRc(GEN x, GEN y) {
      58   105741697 :   GEN z = cgetg(3,t_COMPLEX);
      59   105741760 :   gel(z,1) = gadd(x,gel(y,1));
      60   105740848 :   gel(z,2) = gcopy(gel(y,2)); return z;
      61             : }
      62             : static GEN
      63   116573250 : mulRc(GEN x, GEN y) {
      64   116573250 :   GEN z = cgetg(3,t_COMPLEX);
      65   116573040 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      66   116573410 :   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      725046 : divRc(GEN x, GEN y) {
      77      725046 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      78      725046 :   GEN z = cgetg(3,t_COMPLEX);
      79      725046 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      80      725046 :   gel(z,2) = gmul(mt, gel(y,2));
      81      725046 :   return z;
      82             : }
      83             : static GEN
      84     8519407 : divcR(GEN x, GEN y) {
      85     8519407 :   GEN z = cgetg(3,t_COMPLEX);
      86     8519407 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      87     8519407 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      88             : }
      89             : static GEN
      90        1120 : addRq(GEN x, GEN y) {
      91        1120 :   GEN z = cgetg(4,t_QUAD);
      92        1120 :   gel(z,1) = ZX_copy(gel(y,1));
      93        1120 :   gel(z,2) = gadd(x, gel(y,2));
      94        1120 :   gel(z,3) = gcopy(gel(y,3)); return z;
      95             : }
      96             : static GEN
      97        2191 : mulRq(GEN x, GEN y) {
      98        2191 :   GEN z = cgetg(4,t_QUAD);
      99        2191 :   gel(z,1) = ZX_copy(gel(y,1));
     100        2191 :   gel(z,2) = gmul(x,gel(y,2));
     101        2191 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     102             : }
     103             : static GEN
     104          63 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     105          63 :   long i = gexpo(x) - gexpo(y);
     106          63 :   if (i > 0) prec += nbits2extraprec( i );
     107          63 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     108             : }
     109             : static GEN
     110    12365490 : mulrfrac(GEN x, GEN y)
     111             : {
     112             :   pari_sp av;
     113    12365490 :   GEN z, a = gel(y,1), b = gel(y,2);
     114    12365490 :   if (is_pm1(a)) /* frequent special case */
     115             :   {
     116     3423989 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     117     3423912 :     return z;
     118             :   }
     119     8941501 :   av = avma;
     120     8941501 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     121             : }
     122             : static GEN
     123          21 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     124          21 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     125             : }
     126             : static GEN
     127          63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     128          63 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     129             : }
     130             : static GEN
     131          42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     132          42 :   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      894547 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     148      894547 :   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        2184 : divTp(GEN x, GEN y) { pari_sp av = avma;
     153        2184 :   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       11410 : divpT(GEN x, GEN y) { pari_sp av = avma;
     159       11410 :   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     1886656 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     166     1886656 :   if (lgefint(X) == 3) {
     167     1796676 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     168     1796676 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     169             :   }
     170             :   else {
     171       89980 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     172       89964 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     173             :   }
     174     1886653 :   gel(z,1) = icopy(X); return z;
     175             : }
     176             : static GEN
     177     1156725 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     178     1156725 :   if (lgefint(X) == 3) {
     179      779302 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     180      779302 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     181             :   }
     182             :   else {
     183      377423 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     184      377419 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     185             :   }
     186     1156725 :   gel(z,1) = icopy(X); return z;
     187             : }
     188             : /* cf add_intmod_same */
     189             : static GEN
     190     2731238 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     191     2731238 :   if (lgefint(X) == 3) {
     192     2614874 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     193     2614874 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     194             :   }
     195             :   else
     196      116364 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     197     2731248 :   gel(z,1) = icopy(X); return z;
     198             : }
     199             : /* cf add_intmod_same */
     200             : static GEN
     201       43436 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     202             : {
     203       43436 :   if (lgefint(X) == 3) {
     204       20621 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     205       20614 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     206             :   }
     207             :   else
     208       22815 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     209       43430 :   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     9824551 : rfrac_denom_mul_scal(GEN d, GEN y)
     222             : {
     223     9824551 :   GEN D = RgX_Rg_mul(d, y);
     224     9824551 :   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     9824551 :   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    57774522 : gred_rfrac_simple(GEN n, GEN d)
     235             : {
     236             :   GEN c, cn, cd, z;
     237    57774522 :   long dd = degpol(d);
     238             : 
     239    57774522 :   if (dd <= 0)
     240             :   {
     241        4417 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     242        4417 :     n = gdiv(n, gel(d,2));
     243        4417 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     244        4417 :     return n;
     245             :   }
     246             : 
     247    57770105 :   cd = content(d);
     248    59636020 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     249    57770105 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     250    57770105 :   if (!gequal1(cd)) {
     251     6574556 :     d = RgX_Rg_div(d,cd);
     252     6574556 :     if (!gequal1(cn))
     253             :     {
     254     1302444 :       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     1302395 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     260     1302395 :         c = gdiv(cn,cd);
     261             :       }
     262             :     }
     263             :     else
     264     5272112 :       c = ginv(cd);
     265             :   } else {
     266    51195549 :     if (!gequal1(cn))
     267             :     {
     268     3220759 :       if (gequal0(cn)) {
     269         910 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     270          21 :         c = gen_1;
     271             :       } else {
     272     3219849 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     273     3219849 :         c = cn;
     274             :       }
     275             :     } else {
     276    47974790 :       GEN y = cgetg(3,t_RFRAC);
     277    47974790 :       gel(y,1) = gcopy(n);
     278    47974790 :       gel(y,2) = RgX_copy(d); return y;
     279             :     }
     280             :   }
     281             : 
     282     9794377 :   if (typ(c) == t_POL)
     283             :   {
     284      904048 :     z = c;
     285      943185 :     do { z = content(z); } while (typ(z) == t_POL);
     286      904048 :     cd = denom_i(z);
     287      904048 :     cn = gmul(c, cd);
     288             :   }
     289             :   else
     290             :   {
     291     8890329 :     cn = numer_i(c);
     292     8890329 :     cd = denom_i(c);
     293             :   }
     294     9794377 :   z = cgetg(3,t_RFRAC);
     295     9794377 :   gel(z,1) = gmul(n, cn);
     296     9794377 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     297     9794377 :   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      248955 : fix_rfrac(GEN x, long d)
     303             : {
     304             :   GEN z, N, D;
     305      248955 :   if (!d || typ(x) == t_POL) return x;
     306      175469 :   z = cgetg(3, t_RFRAC);
     307      175469 :   N = gel(x,1);
     308      175469 :   D = gel(x,2);
     309      175469 :   if (d > 0) {
     310       11683 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     311      177793 :                                                    : monomialcopy(N,d,varn(D));
     312      166110 :     gel(z, 2) = RgX_copy(D);
     313             :   } else {
     314        9359 :     gel(z, 1) = gcopy(N);
     315        9359 :     gel(z, 2) = RgX_shift(D, -d);
     316             :   }
     317      175469 :   return z;
     318             : }
     319             : 
     320             : /* assume d != 0 */
     321             : static GEN
     322    44752477 : gred_rfrac2(GEN n, GEN d)
     323             : {
     324             :   GEN y, z;
     325             :   long v, vd, vn;
     326             : 
     327    44752477 :   n = simplify_shallow(n);
     328    44752477 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     329    37922031 :   d = simplify_shallow(d);
     330    37922031 :   if (typ(d) != t_POL) return gdiv(n,d);
     331    36742340 :   vd = varn(d);
     332    36742340 :   if (typ(n) != t_POL)
     333             :   {
     334    20616115 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     335    20614708 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     336           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     337             :   }
     338    16126225 :   vn = varn(n);
     339    16126225 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     340    15988579 :   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    15828827 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     344    15828827 :   if (!degpol(d))
     345             :   {
     346    12756290 :     n = RgX_Rg_div(n,gel(d,2));
     347    12756290 :     return v? RgX_mulXn(n,v): n;
     348             :   }
     349             : 
     350             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     351     3072537 :   if (!isinexact(n) && !isinexact(d))
     352             :   {
     353     3072264 :     y = RgX_divrem(n, d, &z);
     354     3072264 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     355      248682 :     z = RgX_gcd(d, z);
     356      248682 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     357             :   }
     358      248955 :   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    74829719 : Qdivii(GEN x, GEN y)
     364             : {
     365    74829719 :   pari_sp av = avma;
     366             :   GEN r, q;
     367             : 
     368    74829719 :   if (lgefint(y) == 3)
     369             :   {
     370    62171187 :     q = Qdiviu(x, y[2]);
     371    62170874 :     if (signe(y) > 0) return q;
     372    10631807 :     if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
     373    10632023 :     return q;
     374             :   }
     375    12658532 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     376    12658752 :   if (equali1(x))
     377             :   {
     378     5173356 :     if (!signe(y)) pari_err_INV("gdiv",y);
     379     5173251 :     retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
     380             :   }
     381     7485411 :   q = dvmdii(x,y,&r);
     382     7485246 :   if (r == gen_0) return q; /* gen_0 intended */
     383     4672572 :   r = gcdii(y, r);
     384     4672489 :   if (lgefint(r) == 3)
     385             :   {
     386     4296692 :     ulong t = r[2];
     387     4296692 :     set_avma(av);
     388     4296694 :     if (t == 1) q = mkfraccopy(x,y);
     389             :     else
     390             :     {
     391     1830635 :       q = cgetg(3,t_FRAC);
     392     1830850 :       gel(q,1) = diviuexact(x,t);
     393     1830755 :       gel(q,2) = diviuexact(y,t);
     394             :     }
     395             :   }
     396             :   else
     397             :   { /* rare: r and q left on stack for efficiency */
     398      375797 :     q = cgetg(3,t_FRAC);
     399      375801 :     gel(q,1) = diviiexact(x,r);
     400      375802 :     gel(q,2) = diviiexact(y,r);
     401             :   }
     402     4672496 :   normalize_frac(q); return q;
     403             : }
     404             : 
     405             : /* x t_INT, return x/y in reduced form */
     406             : GEN
     407    70075005 : Qdiviu(GEN x, ulong y)
     408             : {
     409    70075005 :   pari_sp av = avma;
     410             :   ulong r, t;
     411             :   GEN q;
     412             : 
     413    70075005 :   if (y == 1) return icopy(x);
     414    53664067 :   if (!y) pari_err_INV("gdiv",gen_0);
     415    53664067 :   if (equali1(x)) retmkfrac(gen_1, utoipos(y));
     416    49680222 :   q = absdiviu_rem(x,y,&r);
     417    49679581 :   if (!r)
     418             :   {
     419    22202553 :     if (signe(x) < 0) togglesign(q);
     420    22202637 :     return q;
     421             :   }
     422    27477028 :   t = ugcd(y, r); set_avma(av);
     423    27477389 :   if (t == 1) retmkfrac(icopy(x), utoipos(y));
     424     9735254 :   retmkfrac(diviuexact(x,t), utoipos(y / t));
     425             : }
     426             : 
     427             : /* x t_INT, return x/y in reduced form */
     428             : GEN
     429     7982866 : Qdivis(GEN x, long y)
     430             : {
     431     7982866 :   pari_sp av = avma;
     432             :   ulong r, t;
     433             :   long s;
     434             :   GEN q;
     435             : 
     436     7982866 :   if (y > 0) return Qdiviu(x, y);
     437       79507 :   if (!y) pari_err_INV("gdiv",gen_0);
     438       79511 :   s = signe(x);
     439       79511 :   if (!s) return gen_0;
     440       61934 :   if (y < 0) { y = -y; s = -s; }
     441       61934 :   if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
     442       61654 :   if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
     443       60443 :   q = absdiviu_rem(x,y,&r);
     444       60443 :   if (!r)
     445             :   {
     446       37707 :     if (s < 0) togglesign(q);
     447       37707 :     return q;
     448             :   }
     449       22736 :   t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
     450       22736 :   if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
     451       22736 :   gel(q,1) = x; setsigne(x, s);
     452       22736 :   gel(q,2) = utoipos(y); return q;
     453             : }
     454             : 
     455             : /*******************************************************************/
     456             : /*                                                                 */
     457             : /*                          CONJUGATION                            */
     458             : /*                                                                 */
     459             : /*******************************************************************/
     460             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     461             : static GEN
     462       16366 : quad_polmod_conj(GEN x, GEN y)
     463             : {
     464             :   GEN z, u, v, a, b;
     465       16366 :   if (typ(x) != t_POL) return gcopy(x);
     466       16366 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     467       16366 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     468       16366 :   b = gel(y,3); v = gel(x,2);
     469       16366 :   z = cgetg(4, t_POL); z[1] = x[1];
     470       16366 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     471       16366 :   gel(z,3) = gneg(u); return z;
     472             : }
     473             : static GEN
     474       16366 : quad_polmod_norm(GEN x, GEN y)
     475             : {
     476             :   GEN z, u, v, a, b, c;
     477       16366 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     478       16366 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     479       16366 :   b = gel(y,3); v = gel(x,2);
     480       16366 :   c = gel(y,2);
     481       16366 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     482       16366 :   if (!gequal1(a)) z = gdiv(z, a);
     483       16366 :   return gadd(z, gsqr(v));
     484             : }
     485             : 
     486             : GEN
     487    11879764 : conj_i(GEN x)
     488             : {
     489             :   long lx, i;
     490             :   GEN y;
     491             : 
     492    11879764 :   switch(typ(x))
     493             :   {
     494     5332283 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     495     5332283 :       return x;
     496             : 
     497     6507721 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     498             : 
     499         812 :     case t_QUAD:
     500         812 :       y = cgetg(4,t_QUAD);
     501         812 :       gel(y,1) = gel(x,1);
     502         812 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     503         812 :                                       : gadd(gel(x,2), gel(x,3));
     504         812 :       gel(y,3) = gneg(gel(x,3)); return y;
     505             : 
     506        9205 :     case t_POL: case t_SER:
     507        9205 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     508       30317 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     509        9205 :       return y;
     510             : 
     511       29743 :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     512       29743 :       y = cgetg_copy(x, &lx);
     513      107583 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     514       29743 :       return y;
     515             : 
     516           0 :     case t_POLMOD:
     517             :     {
     518           0 :       GEN X = gel(x,1);
     519           0 :       long d = degpol(X);
     520           0 :       if (d < 2) return x;
     521           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     522             :     }
     523             :   }
     524           0 :   pari_err_TYPE("gconj",x);
     525             :   return NULL; /* LCOV_EXCL_LINE */
     526             : }
     527             : GEN
     528       69071 : gconj(GEN x)
     529       69071 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     530             : 
     531             : GEN
     532          84 : conjvec(GEN x,long prec)
     533             : {
     534             :   long lx, s, i;
     535             :   GEN z;
     536             : 
     537          84 :   switch(typ(x))
     538             :   {
     539           0 :     case t_INT: case t_INTMOD: case t_FRAC:
     540           0 :       return mkcolcopy(x);
     541             : 
     542           0 :     case t_COMPLEX: case t_QUAD:
     543           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     544             : 
     545          28 :     case t_FFELT:
     546          28 :       return FF_conjvec(x);
     547             : 
     548           0 :     case t_VEC: case t_COL:
     549           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     550           0 :       if (lx == 1) return z;
     551           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     552           0 :       s = lgcols(z);
     553           0 :       for (i=2; i<lx; i++)
     554             :       {
     555           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     556           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     557             :       }
     558           0 :       break;
     559             : 
     560          56 :     case t_POLMOD: {
     561          56 :       GEN T = gel(x,1), r;
     562             :       pari_sp av;
     563             : 
     564          56 :       lx = lg(T);
     565          56 :       if (lx <= 3) return cgetg(1,t_COL);
     566          56 :       x = gel(x,2);
     567         238 :       for (i=2; i<lx; i++)
     568             :       {
     569         189 :         GEN c = gel(T,i);
     570         189 :         switch(typ(c)) {
     571           7 :           case t_INTMOD: {
     572           7 :             GEN p = gel(c,1);
     573             :             pari_sp av;
     574           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     575           7 :             av = avma;
     576           7 :             T = RgX_to_FpX(T,p);
     577           7 :             x = RgX_to_FpX(x, p);
     578           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     579           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     580           7 :             return gerepileupto(av, z);
     581             :           }
     582         182 :           case t_INT:
     583         182 :           case t_FRAC: break;
     584           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     585             :         }
     586             :       }
     587          49 :       if (typ(x) != t_POL)
     588             :       {
     589           0 :         if (!is_rational_t(typ(x)))
     590           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     591           0 :         retconst_col(lx-3, gcopy(x));
     592             :       }
     593          49 :       RgX_check_QX(x,"conjvec");
     594          49 :       av = avma;
     595          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     596          49 :       r = cleanroots(T,prec);
     597          49 :       z = cgetg(lx-2,t_COL);
     598         182 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     599          49 :       return gerepileupto(av, z);
     600             :     }
     601             : 
     602           0 :     default:
     603           0 :       pari_err_TYPE("conjvec",x);
     604             :       return NULL; /* LCOV_EXCL_LINE */
     605             :   }
     606           0 :   return z;
     607             : }
     608             : 
     609             : /********************************************************************/
     610             : /**                                                                **/
     611             : /**                           ADDITION                             **/
     612             : /**                                                                **/
     613             : /********************************************************************/
     614             : /* x, y compatible PADIC, op = add or sub */
     615             : static GEN
     616     1693575 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     617             : {
     618     1693575 :   pari_sp av = avma;
     619             :   long d,e,r,rx,ry;
     620     1693575 :   GEN u, z, mod, p = gel(x,2);
     621             :   int swap;
     622             : 
     623     1693575 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     624     1693575 :   e = valp(x);
     625     1693575 :   r = valp(y); d = r-e;
     626     1693575 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     627     1693575 :   rx = precp(x);
     628     1693575 :   ry = precp(y);
     629     1693575 :   if (d) /* v(x) < v(y) */
     630             :   {
     631      861700 :     r = d+ry; z = powiu(p,d);
     632      861700 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     633      861700 :     z = mulii(z,gel(y,4));
     634      861700 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     635             :   }
     636             :   else
     637             :   {
     638             :     long c;
     639      831875 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     640      831875 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     641      831875 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     642             :     {
     643       66794 :       set_avma(av); return zeropadic(p, e+r);
     644             :     }
     645      765081 :     if (c)
     646             :     {
     647      175364 :       mod = diviiexact(mod, powiu(p,c));
     648      175364 :       r -= c;
     649      175364 :       e += c;
     650             :     }
     651             :   }
     652     1626781 :   u = modii(u, mod);
     653     1626781 :   set_avma(av); z = cgetg(5,t_PADIC);
     654     1626781 :   z[1] = evalprecp(r) | evalvalp(e);
     655     1626781 :   gel(z,2) = icopy(p);
     656     1626781 :   gel(z,3) = icopy(mod);
     657     1626781 :   gel(z,4) = icopy(u); return z;
     658             : }
     659             : /* Rg_to_Fp(t_FRAC) without GC */
     660             : static GEN
     661      934489 : Q_to_Fp(GEN x, GEN p)
     662      934489 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     663             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     664             : static GEN
     665     3462543 : addQp(GEN x, GEN y)
     666             : {
     667     3462543 :   pari_sp av = avma;
     668     3462543 :   long d, r, e, vy = valp(y), py = precp(y);
     669     3462543 :   GEN z, q, mod, u, p = gel(y,2);
     670             : 
     671     3462543 :   e = Q_pvalrem(x, p, &x);
     672     3462478 :   d = vy - e; r = d + py;
     673     3462478 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     674     3460539 :   mod = gel(y,3);
     675     3460539 :   u   = gel(y,4);
     676     3460539 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     677             : 
     678     3460597 :   if (d > 0)
     679             :   {
     680     1846994 :     q = powiu(p,d);
     681     1846953 :     mod = mulii(mod, q);
     682     1846955 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     683     1846955 :     u = addii(x,  mulii(u, q));
     684             :   }
     685     1613603 :   else if (d < 0)
     686             :   {
     687      207493 :     q = powiu(p,-d);
     688      207491 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     689      207491 :     u = addii(u, mulii(x, q));
     690      207489 :     r = py; e = vy;
     691             :   }
     692             :   else
     693             :   {
     694             :     long c;
     695     1406110 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     696     1406110 :     u = addii(u, x);
     697     1406106 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     698             :     {
     699         973 :       set_avma(av); return zeropadic(p,e+r);
     700             :     }
     701     1405133 :     if (c)
     702             :     {
     703      639389 :       mod = diviiexact(mod, powiu(p,c));
     704      639390 :       r -= c;
     705      639390 :       e += c;
     706             :     }
     707             :   }
     708     3459572 :   u = modii(u, mod); set_avma(av);
     709     3459564 :   z = cgetg(5,t_PADIC);
     710     3459609 :   z[1] = evalprecp(r) | evalvalp(e);
     711     3459619 :   gel(z,2) = icopy(p);
     712     3459627 :   gel(z,3) = icopy(mod);
     713     3459635 :   gel(z,4) = icopy(u); return z;
     714             : }
     715             : 
     716             : /* Mod(x,X) + Mod(y,X) */
     717             : #define addsub_polmod_same addsub_polmod_scal
     718             : /* Mod(x,X) +/- Mod(y,Y) */
     719             : static GEN
     720        2142 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     721             : {
     722        2142 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     723        2142 :   GEN z = cgetg(3,t_POLMOD);
     724        2142 :   long vx = varn(X), vy = varn(Y);
     725        2142 :   if (vx==vy) {
     726             :     pari_sp av;
     727          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     728          14 :     warn_coercion(X,Y,gel(z,1));
     729          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     730             :   }
     731        2128 :   if (varncmp(vx, vy) < 0)
     732        2128 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     733             :   else
     734           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     735        2128 :   gel(z,2) = op(x, y); return z;
     736             : }
     737             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     738             : static GEN
     739    11827433 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     740    11827433 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     741             : 
     742             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     743             : static GEN
     744      372178 : add_ser_scal(GEN y, GEN x)
     745             : {
     746             :   long i, l, ly, vy;
     747             :   GEN z;
     748             : 
     749      372178 :   if (isrationalzero(x)) return gcopy(y);
     750      346509 :   ly = lg(y);
     751      346509 :   l = valp(y);
     752      346509 :   if (l < 3-ly) return gcopy(y);
     753             :   /* l + ly >= 3 */
     754      345837 :   if (l < 0)
     755             :   {
     756        1127 :     z = cgetg(ly,t_SER); z[1] = y[1];
     757        3192 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     758        1127 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     759        3199 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     760        1127 :     return z;
     761             :   }
     762      344710 :   vy = varn(y);
     763      344710 :   if (l > 0)
     764             :   {
     765       16898 :     if (ser_isexactzero(y))
     766        8925 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     767        7973 :     y -= l; ly += l;
     768        7973 :     z = cgetg(ly,t_SER);
     769        7973 :     x = gcopy(x);
     770       17409 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     771             :   }
     772             :   else
     773             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     774      327812 :     z = cgetg(ly,t_SER);
     775      327812 :     x = gadd(x, gel(y,2));
     776      327812 :     i = 3;
     777             :   }
     778     1369325 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     779      335785 :   gel(z,2) = x;
     780      335785 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     781      335785 :   return gequal0(x)? normalize(z): z;
     782             : }
     783             : static long
     784     4794386 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     785             : /* x,y t_SER in the same variable: x+y */
     786             : static GEN
     787     2397655 : ser_add(GEN x, GEN y)
     788             : {
     789     2397655 :   long i, lx,ly, n = valp(y) - valp(x);
     790             :   GEN z;
     791     2397655 :   if (n < 0) { n = -n; swap(x,y); }
     792             :   /* valp(x) <= valp(y) */
     793     2397655 :   lx = _serprec(x);
     794     2397655 :   if (lx == 2) /* don't lose type information */
     795         924 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
     796     2396731 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     797     2396731 :   if (n)
     798             :   {
     799       95196 :     if (n+2 > lx) return gcopy(x);
     800       94202 :     z = cgetg(ly,t_SER);
     801      642996 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     802      437751 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     803             :   } else {
     804     2301535 :     z = cgetg(ly,t_SER);
     805    12411279 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     806             :   }
     807     2395737 :   z[1] = x[1]; return normalize(z);
     808             : }
     809             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     810             : static GEN
     811     8823623 : add_rfrac_scal(GEN y, GEN x)
     812             : {
     813             :   pari_sp av;
     814             :   GEN n;
     815             : 
     816     8823623 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     817     5150040 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     818     5150040 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     819             : }
     820             : 
     821             : /* x "scalar", ty != t_MAT and non-scalar */
     822             : static GEN
     823    17181849 : add_scal(GEN y, GEN x, long ty)
     824             : {
     825    17181849 :   switch(ty)
     826             :   {
     827    12330095 :     case t_POL: return RgX_Rg_add(y, x);
     828      372157 :     case t_SER: return add_ser_scal(y, x);
     829     4468337 :     case t_RFRAC: return add_rfrac_scal(y, x);
     830           0 :     case t_COL: return RgC_Rg_add(y, x);
     831       11253 :     case t_VEC:
     832       11253 :       if (isintzero(x)) return gcopy(y);
     833         175 :       break;
     834             :   }
     835         182 :   pari_err_TYPE2("+",x,y);
     836             :   return NULL; /* LCOV_EXCL_LINE */
     837             : }
     838             : 
     839             : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
     840             : static GEN
     841    10540387 : setfrac(GEN z, GEN a, GEN b)
     842             : {
     843    10540387 :   gel(z,1) = icopy_avma(a, (pari_sp)z);
     844    10540387 :   gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
     845    10540387 :   set_avma((pari_sp)gel(z,2)); return z;
     846             : }
     847             : /* z <- a / (b*Q), (Q,a) = 1 */
     848             : static GEN
     849     9496019 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
     850             : {
     851     9496019 :   GEN q = Qdivii(a, b); /* != 0 */
     852     9496019 :   if (typ(q) == t_INT)
     853             :   {
     854      888085 :     gel(z,1) = gerepileuptoint((pari_sp)Q, q);
     855      888085 :     gel(z,2) = Q; return z;
     856             :   }
     857     8607934 :   return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
     858             : }
     859             : static GEN
     860    19061619 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     861             : {
     862    19061619 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     863    19061619 :   GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
     864    19061619 :   int s = cmpii(x2, y2);
     865             : 
     866             :   /* common denominator: (x1 op y1) / x2 */
     867    19061619 :   if (!s)
     868             :   {
     869     6347038 :     pari_sp av = avma;
     870     6347038 :     return gerepileupto(av, Qdivii(op(x1, y1), x2));
     871             :   }
     872    12714581 :   z = cgetg(3, t_FRAC);
     873    12714581 :   if (s < 0)
     874             :   {
     875     6815163 :     Q = dvmdii(y2, x2, &r);
     876             :     /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     877     6815163 :     if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
     878     2138486 :     delta = gcdii(x2,r);
     879             :   }
     880             :   else
     881             :   {
     882     5899418 :     Q = dvmdii(x2, y2, &r);
     883             :     /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     884     5899418 :     if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
     885     1080076 :     delta = gcdii(y2,r);
     886             :   }
     887             :   /* delta = gcd(x2,y2) */
     888     3218562 :   if (equali1(delta))
     889             :   { /* numerator is non-zero */
     890     1286109 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     891     1286109 :     gel(z,2) = mulii(x2,y2); return z;
     892             :   }
     893     1932453 :   x2 = diviiexact(x2,delta);
     894     1932453 :   y2 = diviiexact(y2,delta);
     895     1932453 :   n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
     896     1932453 :   q = dvmdii(n, delta, &r);
     897     1932453 :   if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
     898     1760781 :   r = gcdii(delta, r);
     899     1760781 :   if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
     900     1760781 :   return setfrac(z, n, mulii(mulii(x2, y2), delta));
     901             : }
     902             : 
     903             : /* assume x2, y2 are t_POLs in the same variable */
     904             : static GEN
     905     3037465 : add_rfrac(GEN x, GEN y)
     906             : {
     907     3037465 :   pari_sp av = avma;
     908     3037465 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     909     3037465 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     910             : 
     911     3037465 :   delta = RgX_gcd(x2,y2);
     912     3037465 :   if (!degpol(delta))
     913             :   {
     914         665 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     915         665 :     d = RgX_mul(x2, y2);
     916         665 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     917             :   }
     918     3036800 :   x2 = RgX_div(x2,delta);
     919     3036800 :   y2 = RgX_div(y2,delta);
     920     3036800 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     921     3036800 :   if (!signe(n))
     922             :   {
     923      721152 :     n = simplify_shallow(n);
     924      721152 :     if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
     925           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     926             :   }
     927     2315648 :   if (degpol(n) == 0)
     928     1149728 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     929     1165920 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     930     1165920 :   if (isexactzero(r))
     931             :   {
     932             :     GEN z;
     933      227996 :     d = RgX_mul(x2, y2);
     934             :     /* "constant" denominator ? */
     935      227996 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     936      227996 :     return gerepileupto(av, z);
     937             :   }
     938      937924 :   r = RgX_gcd(delta, r);
     939      937924 :   if (degpol(r))
     940             :   {
     941      160705 :     n = RgX_div(n, r);
     942      160705 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     943             :   }
     944             :   else
     945      777219 :     d = RgX_mul(gel(x,2), y2);
     946      937924 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     947             : }
     948             : 
     949             : GEN
     950  2088819808 : gadd(GEN x, GEN y)
     951             : {
     952  2088819808 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     953             :   pari_sp av, tetpil;
     954             :   GEN z, p1;
     955             : 
     956  2088819808 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     957             :   {
     958   839126822 :     case t_INT: return addii(x,y);
     959   698221476 :     case t_REAL: return addrr(x,y);
     960     1591065 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     961     1591065 :       z = cgetg(3,t_INTMOD);
     962     1591065 :       if (X==Y || equalii(X,Y))
     963     1591051 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     964          14 :       gel(z,1) = gcdii(X,Y);
     965          14 :       warn_coercion(X,Y,gel(z,1));
     966          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     967          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     968             :     }
     969    14842284 :     case t_FRAC: return addsub_frac(x,y,addii);
     970   125161848 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     971   125161183 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     972   125154856 :       if (isintzero(gel(z,2)))
     973             :       {
     974      359646 :         set_avma((pari_sp)(z+3));
     975      359646 :         return gadd(gel(x,1),gel(y,1));
     976             :       }
     977   124795555 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     978   124796299 :       return z;
     979     1171976 :     case t_PADIC:
     980     1171976 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     981     1171976 :       return addsub_pp(x,y, addii);
     982         672 :     case t_QUAD: z = cgetg(4,t_QUAD);
     983         672 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     984         672 :       gel(z,1) = ZX_copy(gel(x,1));
     985         672 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     986         672 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     987    10271692 :     case t_POLMOD:
     988    10271692 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     989    10269585 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     990        2107 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     991    14130204 :     case t_FFELT: return FF_add(x,y);
     992    32932027 :     case t_POL:
     993    32932027 :       vx = varn(x);
     994    32932027 :       vy = varn(y);
     995    32932027 :       if (vx != vy) {
     996      870604 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     997       64685 :         else                     return RgX_Rg_add(y, x);
     998             :       }
     999    32061423 :       return RgX_add(x, y);
    1000     2392468 :     case t_SER:
    1001     2392468 :       vx = varn(x);
    1002     2392468 :       vy = varn(y);
    1003     2392468 :       if (vx != vy) {
    1004          21 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1005          21 :         else                     return add_ser_scal(y, x);
    1006             :       }
    1007     2392447 :       return ser_add(x, y);
    1008     4329352 :     case t_RFRAC:
    1009     4329352 :       vx = varn(gel(x,2));
    1010     4329352 :       vy = varn(gel(y,2));
    1011     4329352 :       if (vx != vy) {
    1012     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1013      538397 :         else                     return add_rfrac_scal(y, x);
    1014             :       }
    1015     3037465 :       return add_rfrac(x,y);
    1016      144693 :     case t_VEC:
    1017      144693 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1018      144693 :       return RgV_add(x,y);
    1019      858283 :     case t_COL:
    1020      858283 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1021      858283 :       return RgC_add(x,y);
    1022      671335 :     case t_MAT:
    1023      671335 :       lx = lg(x);
    1024      671335 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1025      671335 :       if (lx == 1) return cgetg(1, t_MAT);
    1026      671335 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1027      671328 :       return RgM_add(x,y);
    1028             : 
    1029           0 :     default: pari_err_TYPE2("+",x,y);
    1030             :   }
    1031             :   /* tx != ty */
    1032   342997334 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1033             : 
    1034   342997334 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1035             :   {
    1036   282196446 :     case t_INT:
    1037             :       switch(ty)
    1038             :       {
    1039   172823256 :         case t_REAL: return addir(x,y);
    1040      278028 :         case t_INTMOD:
    1041      278028 :           z = cgetg(3, t_INTMOD);
    1042      278028 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1043    25395329 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1044    25395313 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1045    25395362 :           gel(z,2) = icopy(gel(y,2)); return z;
    1046    79645544 :         case t_COMPLEX: return addRc(x, y);
    1047     2502285 :         case t_PADIC:
    1048     2502285 :           if (!signe(x)) return gcopy(y);
    1049     2366165 :           return addQp(x,y);
    1050         966 :         case t_QUAD: return addRq(x, y);
    1051     1551163 :         case t_FFELT: return FF_Z_add(y,x);
    1052             :       }
    1053             : 
    1054             :     case t_REAL:
    1055    28371347 :       switch(ty)
    1056             :       {
    1057    12506490 :         case t_FRAC:
    1058    12506490 :           if (!signe(gel(y,1))) return rcopy(x);
    1059    12506490 :           if (!signe(x))
    1060             :           {
    1061        7653 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1062        7653 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1063             :           }
    1064    12498837 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1065    12497748 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1066    15864801 :         case t_COMPLEX: return addRc(x, y);
    1067          56 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1068             : 
    1069           0 :         default: pari_err_TYPE2("+",x,y);
    1070             :       }
    1071             : 
    1072       17640 :     case t_INTMOD:
    1073             :       switch(ty)
    1074             :       {
    1075       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1076       17507 :           z = cgetg(3, t_INTMOD);
    1077       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1078       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1079             :         }
    1080          14 :         case t_FFELT:
    1081          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1082           0 :             pari_err_OP("+",x,y);
    1083          14 :           return FF_Z_add(y,gel(x,2));
    1084          91 :         case t_COMPLEX: return addRc(x, y);
    1085           0 :         case t_PADIC: { GEN X = gel(x,1);
    1086           0 :           z = cgetg(3, t_INTMOD);
    1087           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1088             :         }
    1089          28 :         case t_QUAD: return addRq(x, y);
    1090             :       }
    1091             : 
    1092             :     case t_FRAC:
    1093             :       switch (ty)
    1094             :       {
    1095    10231516 :         case t_COMPLEX: return addRc(x, y);
    1096     1096381 :         case t_PADIC:
    1097     1096381 :           if (!signe(gel(x,1))) return gcopy(y);
    1098     1096381 :           return addQp(x,y);
    1099         126 :         case t_QUAD: return addRq(x, y);
    1100        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1101             :       }
    1102             : 
    1103             :     case t_FFELT:
    1104           0 :       pari_err_TYPE2("+",x,y);
    1105             : 
    1106          21 :     case t_COMPLEX:
    1107             :       switch(ty)
    1108             :       {
    1109          14 :         case t_PADIC:
    1110          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1111           7 :         case t_QUAD:
    1112           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1113           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1114             :       }
    1115             : 
    1116             :     case t_PADIC: /* ty == t_QUAD */
    1117           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1118             :   }
    1119             :   /* tx < ty, !is_const_t(y) */
    1120    21166931 :   switch(ty)
    1121             :   {
    1122        5607 :     case t_MAT:
    1123        5607 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1124        5607 :       if (isrationalzero(x)) return gcopy(y);
    1125        5530 :       return RgM_Rg_add(y, x);
    1126      196203 :     case t_COL:
    1127      196203 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1128      196203 :       return RgC_Rg_add(y, x);
    1129      676458 :     case t_POLMOD: /* is_const_t(tx) in this case */
    1130      676458 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1131             :   }
    1132    20288663 :   if (is_scalar_t(tx))  {
    1133    17197002 :     if (tx == t_POLMOD)
    1134             :     {
    1135       71979 :       vx = varn(gel(x,1));
    1136       71979 :       vy = gvar(y);
    1137       71979 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1138             :       else
    1139       49236 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1140       37772 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1141             :     }
    1142    17125023 :     return add_scal(y, x, ty);
    1143             :   }
    1144             :   /* x and y are not scalars, ty != t_MAT */
    1145     3091660 :   vx = gvar(x);
    1146     3091660 :   vy = gvar(y);
    1147     3091660 :   if (vx != vy) { /* x or y is treated as a scalar */
    1148       22626 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1149       32151 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1150       32151 :                                 : add_scal(y, x, ty);
    1151             :   }
    1152             :   /* vx = vy */
    1153     3069034 :   switch(tx)
    1154             :   {
    1155     3068628 :     case t_POL:
    1156             :       switch (ty)
    1157             :       {
    1158        5229 :         case t_SER:
    1159        5229 :           if (lg(x) == 2) return gcopy(y);
    1160        5215 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1161        5215 :           i = lg(y) + valp(y) - i;
    1162        5215 :           if (i < 3) return gcopy(y);
    1163        5208 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1164        5208 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1165        5208 :           return y;
    1166             : 
    1167     3063399 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1168             :       }
    1169           0 :       break;
    1170             : 
    1171         406 :     case t_SER:
    1172         406 :       if (ty == t_RFRAC)
    1173             :       {
    1174             :         GEN n, d;
    1175             :         long vn, vd;
    1176         406 :         av = avma;
    1177         406 :         n = gel(y,1); vn = gval(n, vy);
    1178         406 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1179             : 
    1180         406 :         l = lg(x) + valp(x) - (vn - vd);
    1181         406 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1182             : 
    1183             :         /* take advantage of y = t^n ! */
    1184         406 :         if (degpol(d))
    1185          70 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1186             :         else {
    1187         336 :           y = gdiv(n, gel(d,2));
    1188         336 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1189             :         }
    1190         406 :         setvalp(y, valp(y) - vd);
    1191         406 :         return gerepileupto(av, gadd(y, x));
    1192             :       }
    1193           0 :       break;
    1194             :   }
    1195           0 :   pari_err_TYPE2("+",x,y);
    1196             :   return NULL; /* LCOV_EXCL_LINE */
    1197             : }
    1198             : 
    1199             : GEN
    1200    32650476 : gaddsg(long x, GEN y)
    1201             : {
    1202    32650476 :   long ty = typ(y);
    1203             :   GEN z;
    1204             : 
    1205    32650476 :   switch(ty)
    1206             :   {
    1207    12284202 :     case t_INT:  return addsi(x,y);
    1208     7269096 :     case t_REAL: return addsr(x,y);
    1209          14 :     case t_INTMOD:
    1210          14 :       z = cgetg(3, t_INTMOD);
    1211          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1212     5006681 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1213     5006681 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1214     5006681 :       gel(z,2) = icopy(gel(y,2)); return z;
    1215     7597941 :     case t_COMPLEX:
    1216     7597941 :       z = cgetg(3, t_COMPLEX);
    1217     7597941 :       gel(z,1) = gaddsg(x, gel(y,1));
    1218     7597941 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1219             : 
    1220      492542 :     default: return gadd(stoi(x), y);
    1221             :   }
    1222             : }
    1223             : 
    1224             : GEN
    1225     1184203 : gsubsg(long x, GEN y)
    1226             : {
    1227             :   GEN z, a, b;
    1228             :   pari_sp av;
    1229             : 
    1230     1184203 :   switch(typ(y))
    1231             :   {
    1232       39760 :     case t_INT:  return subsi(x,y);
    1233      744994 :     case t_REAL: return subsr(x,y);
    1234          56 :     case t_INTMOD:
    1235          56 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1236          56 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1237       25627 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1238       25627 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1239       25627 :       gel(z,2) = icopy(gel(y,2)); return z;
    1240      344268 :     case t_COMPLEX:
    1241      344268 :       z = cgetg(3, t_COMPLEX);
    1242      344268 :       gel(z,1) = gsubsg(x, gel(y,1));
    1243      344268 :       gel(z,2) = gneg(gel(y,2)); return z;
    1244             :   }
    1245       29498 :   av = avma;
    1246       29498 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1247             : }
    1248             : 
    1249             : /********************************************************************/
    1250             : /**                                                                **/
    1251             : /**                          SUBTRACTION                           **/
    1252             : /**                                                                **/
    1253             : /********************************************************************/
    1254             : 
    1255             : GEN
    1256  2142885010 : gsub(GEN x, GEN y)
    1257             : {
    1258  2142885010 :   long tx = typ(x), ty = typ(y);
    1259             :   pari_sp av;
    1260             :   GEN z;
    1261  2142885010 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1262             :   {
    1263  1802464812 :     case t_INT: return subii(x,y);
    1264   253163318 :     case t_REAL: return subrr(x,y);
    1265     1156739 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1266     1156739 :       z = cgetg(3,t_INTMOD);
    1267     1156739 :       if (X==Y || equalii(X,Y))
    1268     1156725 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1269          14 :       gel(z,1) = gcdii(X,Y);
    1270          14 :       warn_coercion(X,Y,gel(z,1));
    1271          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1272          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1273             :     }
    1274     4219335 :     case t_FRAC: return addsub_frac(x,y, subii);
    1275    40846222 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1276    40846222 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1277    40846222 :       if (isintzero(gel(z,2)))
    1278             :       {
    1279       68837 :         set_avma((pari_sp)(z+3));
    1280       68837 :         return gsub(gel(x,1),gel(y,1));
    1281             :       }
    1282    40777385 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1283    40777385 :       return z;
    1284      521592 :     case t_PADIC:
    1285      521592 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1286      521592 :       return addsub_pp(x,y, subii);
    1287         168 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1288         168 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1289         168 :       gel(z,1) = ZX_copy(gel(x,1));
    1290         168 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1291         168 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1292      843653 :     case t_POLMOD:
    1293      843653 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1294      843618 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1295          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1296      282778 :     case t_FFELT: return FF_sub(x,y);
    1297     8105008 :     case t_POL: {
    1298     8105008 :       long vx = varn(x);
    1299     8105008 :       long vy = varn(y);
    1300     8105008 :       if (vx != vy) {
    1301       18915 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1302        5579 :         else                     return Rg_RgX_sub(x, y);
    1303             :       }
    1304     8086093 :       return RgX_sub(x, y);
    1305             :     }
    1306      271413 :     case t_VEC:
    1307      271413 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1308      271413 :       return RgV_sub(x,y);
    1309     5384711 :     case t_COL:
    1310     5384711 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1311     5384711 :       return RgC_sub(x,y);
    1312       43434 :     case t_MAT: {
    1313       43434 :       long lx = lg(x);
    1314       43434 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1315       43434 :       if (lx == 1) return cgetg(1, t_MAT);
    1316       35174 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1317       35174 :       return RgM_sub(x,y);
    1318             :     }
    1319     2220792 :     case t_RFRAC: case t_SER: break;
    1320             : 
    1321           0 :     default: pari_err_TYPE2("+",x,y);
    1322             :   }
    1323    23361042 :   av = avma;
    1324    25581834 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1325             : }
    1326             : 
    1327             : /********************************************************************/
    1328             : /**                                                                **/
    1329             : /**                        MULTIPLICATION                          **/
    1330             : /**                                                                **/
    1331             : /********************************************************************/
    1332             : static GEN
    1333      303682 : mul_ser_scal(GEN y, GEN x) {
    1334             :   long l, i;
    1335             :   GEN z;
    1336      303682 :   if (isrationalzero(x)) return gmul(Rg_get_0(y), x);
    1337      301008 :   if (ser_isexactzero(y))
    1338         448 :     return scalarser(lg(y) == 2? Rg_get_0(x)
    1339         448 :                                : gmul(gel(y,2), x), varn(y), valp(y));
    1340      300560 :   z = cgetg_copy(y, &l); z[1] = y[1];
    1341     3834539 :   for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
    1342      300560 :   return normalize(z);
    1343             : }
    1344             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1345             :  * [n/d a valid RFRAC]  */
    1346             : static GEN
    1347    10444976 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1348             : {
    1349    10444976 :   pari_sp av = avma;
    1350             :   GEN z;
    1351             : 
    1352    10444976 :   switch(typ(x))
    1353             :   {
    1354          14 :     case t_PADIC:
    1355          14 :       n = gmul(n, x);
    1356          14 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1357          14 :       return gerepileupto(av, gdiv(n,d));
    1358             : 
    1359         434 :     case t_INTMOD: case t_POLMOD:
    1360         434 :       n = gmul(n, x);
    1361         434 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1362         434 :       return gerepileupto(av, gdiv(n,d));
    1363             :   }
    1364    10444528 :   z = gred_rfrac2(x, d);
    1365    10444528 :   n = simplify_shallow(n);
    1366    10444528 :   if (typ(z) == t_RFRAC)
    1367             :   {
    1368     7929000 :     n = gmul(gel(z,1), n);
    1369     7929000 :     d = gel(z,2);
    1370     7929000 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1371           0 :       z = RgX_Rg_div(n, d);
    1372             :     else
    1373     7929000 :       z = gred_rfrac_simple(n, d);
    1374             :   }
    1375             :   else
    1376     2515528 :     z = gmul(z, n);
    1377    10444528 :   return gerepileupto(av, z);
    1378             : }
    1379             : static GEN
    1380    73401736 : mul_scal(GEN y, GEN x, long ty)
    1381             : {
    1382    73401736 :   switch(ty)
    1383             :   {
    1384    64524441 :     case t_POL:
    1385    64524441 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1386    61674235 :       return RgX_Rg_mul(y, x);
    1387      295604 :     case t_SER: return mul_ser_scal(y, x);
    1388     8581664 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1389          14 :     case t_QFI: case t_QFR:
    1390          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1391             :   }
    1392          13 :   pari_err_TYPE2("*",x,y);
    1393             :   return NULL; /* LCOV_EXCL_LINE */
    1394             : }
    1395             : 
    1396             : static GEN
    1397      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1398             : {
    1399      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1400      160830 :   long vx = gvar(X), vy = varn(y2);
    1401      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1402        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1403             : }
    1404             : /* (x1/x2) * (y1/y2) */
    1405             : static GEN
    1406     7904440 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1407             : {
    1408             :   GEN z, X, Y;
    1409     7904440 :   pari_sp av = avma;
    1410             : 
    1411     7904440 :   X = gred_rfrac2(x1, y2);
    1412     7904440 :   Y = gred_rfrac2(y1, x2);
    1413     7904440 :   if (typ(X) == t_RFRAC)
    1414             :   {
    1415     6625477 :     if (typ(Y) == t_RFRAC) {
    1416     6559901 :       x1 = gel(X,1);
    1417     6559901 :       x2 = gel(X,2);
    1418     6559901 :       y1 = gel(Y,1);
    1419     6559901 :       y2 = gel(Y,2);
    1420     6559901 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1421             :     } else
    1422       65576 :       z = mul_gen_rfrac(Y, X);
    1423             :   }
    1424     1278963 :   else if (typ(Y) == t_RFRAC)
    1425       95254 :     z = mul_gen_rfrac(X, Y);
    1426             :   else
    1427     1183709 :     z = gmul(X, Y);
    1428     7904440 :   return gerepileupto(av, z);
    1429             : }
    1430             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1431             : static GEN
    1432      269701 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1433             : {
    1434      269701 :   pari_sp av = avma;
    1435      269701 :   GEN X = gred_rfrac2(x1, y2);
    1436      269701 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1437             :   {
    1438      263149 :     x2 = RgX_mul(gel(X,2), x2);
    1439      263149 :     x1 = gel(X,1);
    1440             :   }
    1441             :   else
    1442        6552 :     x1 = X;
    1443      269701 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1444             : }
    1445             : 
    1446             : /* Mod(y, Y) * x,  assuming x scalar */
    1447             : static GEN
    1448     3261446 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1449     3261446 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1450             : 
    1451             : /* cf mulqq */
    1452             : static GEN
    1453     5562335 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1454             : {
    1455     5562335 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1456     5562335 :   pari_sp tetpil, av = avma;
    1457     5562335 :   T[1] = x[1];
    1458     5562335 :   p2 = gmul(gel(x,2), gel(y,2));
    1459     5562335 :   p3 = gmul(gel(x,3), gel(y,3));
    1460     5562335 :   p1 = gmul(gneg_i(c),p3);
    1461             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1462     5562335 :   if (typ(b) == t_INT)
    1463             :   {
    1464     5562314 :     if (signe(b))
    1465             :     {
    1466     4957423 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1467     4957423 :       if (is_pm1(b))
    1468             :       {
    1469     4956765 :         if (signe(b) > 0) p3 = gneg(p3);
    1470             :       }
    1471             :       else
    1472         658 :         p3 = gmul(negi(b), p3);
    1473             :     }
    1474             :     else
    1475             :     {
    1476      604891 :       p3 = gmul(gel(x,2),gel(y,3));
    1477      604891 :       p4 = gmul(gel(x,3),gel(y,2));
    1478             :     }
    1479             :   }
    1480             :   else
    1481             :   {
    1482          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1483          21 :     p3 = gmul(gneg_i(b), p3);
    1484             :   }
    1485     5562335 :   tetpil = avma;
    1486     5562335 :   gel(T,2) = gadd(p2, p1);
    1487     5562335 :   gel(T,3) = gadd(p4, p3);
    1488     5562335 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1489     5562335 :   return normalizepol_lg(T,4);
    1490             : }
    1491             : /* Mod(x,T) * Mod(y,T) */
    1492             : static GEN
    1493     7103136 : mul_polmod_same(GEN T, GEN x, GEN y)
    1494             : {
    1495     7103136 :   GEN z = cgetg(3,t_POLMOD), a;
    1496     7103136 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1497     7103136 :   gel(z,1) = RgX_copy(T);
    1498             :   /* x * y mod T optimised */
    1499     7103136 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1500     6467443 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1501     1103468 :     a = gmul(x, y);
    1502             :   else
    1503             :   {
    1504     5999668 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1505     5558240 :       a = quad_polmod_mul(T, x, y);
    1506             :     else
    1507      441428 :       a = RgXQ_mul(x, y, gel(z,1));
    1508             :   }
    1509     7103136 :   gel(z,2) = a; return z;
    1510             : }
    1511             : static GEN
    1512       43882 : sqr_polmod(GEN T, GEN x)
    1513             : {
    1514       43882 :   GEN a, z = cgetg(3,t_POLMOD);
    1515       43882 :   gel(z,1) = RgX_copy(T);
    1516       43882 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1517       13330 :     a = gsqr(x);
    1518             :   else
    1519             :   {
    1520       30552 :     pari_sp av = avma;
    1521       30552 :     a = RgXQ_sqr(x, gel(z,1));
    1522       30552 :     a = gerepileupto(av, a);
    1523             :   }
    1524       43882 :   gel(z,2) = a; return z;
    1525             : }
    1526             : /* Mod(x,X) * Mod(y,Y) */
    1527             : static GEN
    1528     2661848 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1529             : {
    1530     2661848 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1531     2661848 :   long vx = varn(X), vy = varn(Y);
    1532     2661848 :   GEN z = cgetg(3,t_POLMOD);
    1533             : 
    1534     2661848 :   if (vx==vy) {
    1535             :     pari_sp av;
    1536          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1537          14 :     warn_coercion(X,Y,gel(z,1));
    1538          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1539          14 :     return z;
    1540             :   }
    1541     2661834 :   if (varncmp(vx, vy) < 0)
    1542      124509 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1543             :   else
    1544     2537325 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1545     2661834 :   gel(z,2) = gmul(x, y); return z;
    1546             : }
    1547             : 
    1548             : #if 0 /* used by 3M only */
    1549             : /* set z = x+y and return 1 if x,y have the same sign
    1550             :  * set z = x-y and return 0 otherwise */
    1551             : static int
    1552             : did_add(GEN x, GEN y, GEN *z)
    1553             : {
    1554             :   long tx = typ(x), ty = typ(y);
    1555             :   if (tx == ty) switch(tx)
    1556             :   {
    1557             :     case t_INT: *z = addii(x,y); return 1;
    1558             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1559             :     case t_REAL:
    1560             :       if (signe(x) == -signe(y))
    1561             :       { *z = subrr(x,y); return 0; }
    1562             :       else
    1563             :       { *z = addrr(x,y); return 1; }
    1564             :   }
    1565             :   if (tx == t_REAL) switch(ty)
    1566             :   {
    1567             :     case t_INT:
    1568             :       if (signe(x) == -signe(y))
    1569             :       { *z = subri(x,y); return 0; }
    1570             :       else
    1571             :       { *z = addri(x,y); return 1; }
    1572             :     case t_FRAC:
    1573             :       if (signe(x) == -signe(gel(y,1)))
    1574             :       { *z = gsub(x,y); return 0; }
    1575             :       else
    1576             :       { *z = gadd(x,y); return 1; }
    1577             :   }
    1578             :   else if (ty == t_REAL) switch(tx)
    1579             :   {
    1580             :     case t_INT:
    1581             :       if (signe(x) == -signe(y))
    1582             :       { *z = subir(x,y); return 0; }
    1583             :       else
    1584             :       { *z = addir(x,y); return 1; }
    1585             :     case t_FRAC:
    1586             :       if (signe(gel(x,1)) == -signe(y))
    1587             :       { *z = gsub(x,y); return 0; }
    1588             :       else
    1589             :       { *z = gadd(x,y); return 1; }
    1590             :   }
    1591             :   *z = gadd(x,y); return 1;
    1592             : }
    1593             : #endif
    1594             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1595             : static GEN
    1596     5371795 : mulcIR(GEN x, GEN y)
    1597             : {
    1598     5371795 :   GEN z = cgetg(3,t_COMPLEX);
    1599     5371795 :   pari_sp av = avma;
    1600     5371795 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1601     5371795 :   gel(z,2) = gmul(y, gel(x,1));
    1602     5371795 :   return z;
    1603             : 
    1604             : }
    1605             : /* x,y COMPLEX */
    1606             : static GEN
    1607   135307108 : mulcc(GEN x, GEN y)
    1608             : {
    1609   135307108 :   GEN xr = gel(x,1), xi = gel(x,2);
    1610   135307108 :   GEN yr = gel(y,1), yi = gel(y,2);
    1611             :   GEN p1, p2, p3, p4, z;
    1612             :   pari_sp tetpil, av;
    1613             : 
    1614   135307108 :   if (isintzero(xr))
    1615             :   {
    1616     4695270 :     if (isintzero(yr)) {
    1617      710529 :       av = avma;
    1618      710529 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1619             :     }
    1620     3984741 :     return mulcIR(y, xi);
    1621             :   }
    1622   130611838 :   if (isintzero(yr)) return mulcIR(x, yi);
    1623             : 
    1624   129224784 :   z = cgetg(3,t_COMPLEX); av = avma;
    1625             : #if 0
    1626             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1627             :    * e.g. xr + xi if exponents differ */
    1628             :   if (did_add(xr, xi, &p3))
    1629             :   {
    1630             :     if (did_add(yr, yi, &p4)) {
    1631             :     /* R = xr*yr - xi*yi
    1632             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1633             :       p1 = gmul(xr,yr);
    1634             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1635             :       p3 = gmul(p3, p4);
    1636             :       p4 = gsub(p2, p1);
    1637             :     } else {
    1638             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1639             :      * I = xr*yi + xi*yr */
    1640             :       p1 = gmul(p3,p4);
    1641             :       p3 = gmul(xr,yi);
    1642             :       p4 = gmul(xi,yr);
    1643             :       p2 = gsub(p3, p4);
    1644             :     }
    1645             :   } else {
    1646             :     if (did_add(yr, yi, &p4)) {
    1647             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1648             :       * I = xr*yi +xi*yr */
    1649             :       p1 = gmul(p3,p4);
    1650             :       p3 = gmul(xr,yi);
    1651             :       p4 = gmul(xi,yr);
    1652             :       p2 = gsub(p4, p3);
    1653             :     } else {
    1654             :     /* R = xr*yr - xi*yi
    1655             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1656             :       p3 = gneg( gmul(p3, p4) );
    1657             :       p1 = gmul(xr,yr);
    1658             :       p2 = gmul(xi,yi);
    1659             :       p4 = gadd(p1, p2);
    1660             : 
    1661             :       p2 = gneg(p2);
    1662             :     }
    1663             :   }
    1664             :   tetpil = avma;
    1665             :   gel(z,1) = gadd(p1,p2);
    1666             :   gel(z,2) = gadd(p3,p4);
    1667             : #else
    1668   129224784 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1669             :   { /* 3M formula */
    1670      615990 :     p3 = addii(xr,xi);
    1671      615990 :     p4 = addii(yr,yi);
    1672      615990 :     p1 = mulii(xr,yr);
    1673      615990 :     p2 = mulii(xi,yi);
    1674      615990 :     p3 = mulii(p3,p4);
    1675      615990 :     p4 = addii(p2,p1);
    1676      615990 :     tetpil = avma;
    1677      615990 :     gel(z,1) = subii(p1,p2);
    1678      615990 :     gel(z,2) = subii(p3,p4);
    1679      615990 :     if (!signe(gel(z,2)))
    1680      113015 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1681             :   }
    1682             :   else
    1683             :   { /* naive 4M formula: avoid all loss of accuracy */
    1684   128608794 :     p1 = gmul(xr,yr);
    1685   128608794 :     p2 = gmul(xi,yi);
    1686   128608794 :     p3 = gmul(xr,yi);
    1687   128608794 :     p4 = gmul(xi,yr);
    1688   128608794 :     tetpil = avma;
    1689   128608794 :     gel(z,1) = gsub(p1,p2);
    1690   128608794 :     gel(z,2) = gadd(p3,p4);
    1691   128608794 :     if (isintzero(gel(z,2)))
    1692             :     {
    1693         441 :       cgiv(gel(z,2));
    1694         441 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1695             :     }
    1696             :   }
    1697             : #endif
    1698             : 
    1699   129111328 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1700             : }
    1701             : /* x,y PADIC */
    1702             : static GEN
    1703     5543116 : mulpp(GEN x, GEN y) {
    1704     5543116 :   long l = valp(x) + valp(y);
    1705             :   pari_sp av;
    1706             :   GEN z, t;
    1707     5543116 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1708     5543111 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1709     5344640 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1710             : 
    1711     5266485 :   t = (precp(x) > precp(y))? y: x;
    1712     5266485 :   z = cgetp(t); setvalp(z,l); av = avma;
    1713     5266485 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1714     5266450 :   set_avma(av); return z;
    1715             : }
    1716             : /* x,y QUAD */
    1717             : static GEN
    1718        1071 : mulqq(GEN x, GEN y) {
    1719        1071 :   GEN z = cgetg(4,t_QUAD);
    1720        1071 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1721             :   pari_sp av, tetpil;
    1722        1071 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1723             : 
    1724        1071 :   gel(z,1) = ZX_copy(P); av = avma;
    1725        1071 :   p2 = gmul(gel(x,2),gel(y,2));
    1726        1071 :   p3 = gmul(gel(x,3),gel(y,3));
    1727        1071 :   p1 = gmul(gneg_i(c),p3);
    1728             : 
    1729        1071 :   if (signe(b))
    1730         952 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1731             :   else
    1732             :   {
    1733         119 :     p3 = gmul(gel(x,2),gel(y,3));
    1734         119 :     p4 = gmul(gel(x,3),gel(y,2));
    1735             :   }
    1736        1071 :   tetpil = avma;
    1737        1071 :   gel(z,2) = gadd(p2,p1);
    1738        1071 :   gel(z,3) = gadd(p4,p3);
    1739        1071 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1740             : }
    1741             : 
    1742             : GEN
    1743     3584408 : mulcxI(GEN x)
    1744             : {
    1745             :   GEN z;
    1746     3584408 :   switch(typ(x))
    1747             :   {
    1748      418788 :     case t_INT: case t_REAL: case t_FRAC:
    1749      418788 :       return mkcomplex(gen_0, x);
    1750     3165564 :     case t_COMPLEX:
    1751     3165564 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1752     3165361 :       z = cgetg(3,t_COMPLEX);
    1753     3165361 :       gel(z,1) = gneg(gel(x,2));
    1754     3165361 :       gel(z,2) = gel(x,1); return z;
    1755          56 :     default:
    1756          56 :       return gmul(gen_I(), x);
    1757             :   }
    1758             : }
    1759             : GEN
    1760      109606 : mulcxmI(GEN x)
    1761             : {
    1762             :   GEN z;
    1763      109606 :   switch(typ(x))
    1764             :   {
    1765        3213 :     case t_INT: case t_REAL: case t_FRAC:
    1766        3213 :       return mkcomplex(gen_0, gneg(x));
    1767      106330 :     case t_COMPLEX:
    1768      106330 :       if (isintzero(gel(x,1))) return gel(x,2);
    1769      104643 :       z = cgetg(3,t_COMPLEX);
    1770      104643 :       gel(z,1) = gel(x,2);
    1771      104643 :       gel(z,2) = gneg(gel(x,1)); return z;
    1772          63 :     default:
    1773          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1774             :   }
    1775             : }
    1776             : /* x * I^k */
    1777             : GEN
    1778      265529 : mulcxpowIs(GEN x, long k)
    1779             : {
    1780      265529 :   switch (k & 3)
    1781             :   {
    1782       47404 :     case 1: return mulcxI(x);
    1783       42469 :     case 2: return gneg(x);
    1784       46291 :     case 3: return mulcxmI(x);
    1785             :   }
    1786      129365 :   return x;
    1787             : }
    1788             : 
    1789             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1790             : static GEN
    1791     3806099 : fill_ser(GEN z, GEN y)
    1792             : {
    1793     3806099 :   long i, lx = lg(z), ly = lg(y);
    1794     3806099 :   if (ly >= lx) {
    1795    17561751 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1796             :   } else {
    1797      313757 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1798      183650 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1799             :   }
    1800     3806099 :   return normalize(z);
    1801             : }
    1802             : 
    1803             : GEN
    1804  4139977815 : gmul(GEN x, GEN y)
    1805             : {
    1806             :   long tx, ty, lx, ly, vx, vy, i, l;
    1807             :   pari_sp av, tetpil;
    1808             :   GEN z, p1, p2;
    1809             : 
    1810  4139977815 :   if (x == y) return gsqr(x);
    1811  3260400201 :   tx = typ(x); ty = typ(y);
    1812  3260400201 :   if (tx == ty) switch(tx)
    1813             :   {
    1814  1495304215 :     case t_INT: return mulii(x,y);
    1815   859347035 :     case t_REAL: return mulrr(x,y);
    1816     1759720 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1817     1759720 :       z = cgetg(3,t_INTMOD);
    1818     1759720 :       if (X==Y || equalii(X,Y))
    1819     1759685 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1820          35 :       gel(z,1) = gcdii(X,Y);
    1821          35 :       warn_coercion(X,Y,gel(z,1));
    1822          35 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1823          35 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1824             :     }
    1825     9667530 :     case t_FRAC:
    1826             :     {
    1827     9667530 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1828     9667530 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1829     9667530 :       z=cgetg(3,t_FRAC);
    1830     9667530 :       p1 = gcdii(x1, y2);
    1831     9667530 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1832     9667530 :       p1 = gcdii(x2, y1);
    1833     9667530 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1834     9667530 :       tetpil = avma;
    1835     9667530 :       gel(z,2) = mulii(x2,y2);
    1836     9667530 :       gel(z,1) = mulii(x1,y1);
    1837     9667530 :       fix_frac_if_int_GC(z,tetpil); return z;
    1838             :     }
    1839   131725938 :     case t_COMPLEX: return mulcc(x, y);
    1840     4647002 :     case t_PADIC: return mulpp(x, y);
    1841         840 :     case t_QUAD: return mulqq(x, y);
    1842    14563716 :     case t_FFELT: return FF_mul(x, y);
    1843     9345579 :     case t_POLMOD:
    1844     9345579 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1845     6683731 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1846     2661848 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1847    40913577 :     case t_POL:
    1848    40913577 :       vx = varn(x);
    1849    40913577 :       vy = varn(y);
    1850    40913577 :       if (vx != vy) {
    1851     4843672 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1852     2058961 :         else                     return RgX_Rg_mul(y, x);
    1853             :       }
    1854    36069905 :       return RgX_mul(x, y);
    1855             : 
    1856     2969139 :     case t_SER: {
    1857     2969139 :       vx = varn(x);
    1858     2969139 :       vy = varn(y);
    1859     2969139 :       if (vx != vy) {
    1860        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1861        3675 :         return mul_ser_scal(y, x);
    1862             :       }
    1863     2965464 :       lx = minss(lg(x), lg(y));
    1864     2965464 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1865     2965436 :       av = avma; z = cgetg(lx,t_SER);
    1866     2965436 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1867     2965436 :       x = ser2pol_i(x, lx);
    1868     2965436 :       y = ser2pol_i(y, lx);
    1869     2965436 :       y = RgXn_mul(x, y, lx-2);
    1870     2965436 :       return gerepilecopy(av, fill_ser(z,y));
    1871             :     }
    1872    16116800 :     case t_QFI: return qficomp(x,y);
    1873          14 :     case t_QFR: return qfrcomp(x,y);
    1874     6717896 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1875     2638782 :     case t_MAT: return RgM_mul(x, y);
    1876             : 
    1877        1631 :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1878        1631 :       z = cgetg_copy(x, &l);
    1879        1631 :       if (l != lg(y)) break;
    1880       17325 :       for (i=1; i<l; i++)
    1881             :       {
    1882       15694 :         long yi = y[i];
    1883       15694 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1884       15694 :         z[i] = x[yi];
    1885             :       }
    1886        1631 :       return z;
    1887             : 
    1888        2420 :     default:
    1889        2420 :       pari_err_TYPE2("*",x,y);
    1890             :   }
    1891             :   /* tx != ty */
    1892   664678367 :   if (is_const_t(ty) && is_const_t(tx))  {
    1893   577936353 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1894   577936353 :     switch(tx) {
    1895   528993458 :     case t_INT:
    1896             :       switch(ty)
    1897             :       {
    1898   355631940 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1899      971316 :         case t_INTMOD:
    1900      971316 :           z = cgetg(3, t_INTMOD);
    1901      971314 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1902    43072771 :         case t_FRAC:
    1903    43072771 :           if (!signe(x)) return gen_0;
    1904    28828875 :           z=cgetg(3,t_FRAC);
    1905    28828875 :           p1 = gcdii(x,gel(y,2));
    1906    28828872 :           if (equali1(p1))
    1907             :           {
    1908    15880194 :             set_avma((pari_sp)z);
    1909    15880196 :             gel(z,2) = icopy(gel(y,2));
    1910    15880198 :             gel(z,1) = mulii(gel(y,1), x);
    1911             :           }
    1912             :           else
    1913             :           {
    1914    12948676 :             x = diviiexact(x,p1); tetpil = avma;
    1915    12948676 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1916    12948677 :             gel(z,1) = mulii(gel(y,1), x);
    1917    12948677 :             fix_frac_if_int_GC(z,tetpil);
    1918             :           }
    1919    28828870 :           return z;
    1920   126625304 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1921      866826 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1922        1554 :         case t_QUAD: return mulRq(x,y);
    1923     1824267 :         case t_FFELT: return FF_Z_mul(y,x);
    1924             :       }
    1925             : 
    1926             :     case t_REAL:
    1927    46511690 :       switch(ty)
    1928             :       {
    1929    12365493 :         case t_FRAC: return mulrfrac(x, y);
    1930    34146155 :         case t_COMPLEX: return mulRc(x, y);
    1931          14 :         case t_QUAD: return mulqf(y, x, realprec(x));
    1932          28 :         default: pari_err_TYPE2("*",x,y);
    1933             :       }
    1934             : 
    1935        8057 :     case t_INTMOD:
    1936             :       switch(ty)
    1937             :       {
    1938        7350 :         case t_FRAC: { GEN X = gel(x,1);
    1939        7350 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1940        7350 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1941             :         }
    1942          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1943         259 :         case t_PADIC: { GEN X = gel(x,1);
    1944         259 :           z = cgetg(3, t_INTMOD);
    1945         259 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1946             :         }
    1947          49 :         case t_QUAD: return mulRq(x, y);
    1948         350 :         case t_FFELT:
    1949         350 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1950           0 :             pari_err_OP("*",x,y);
    1951         350 :           return FF_Z_mul(y,gel(x,2));
    1952             :       }
    1953             : 
    1954             :     case t_FRAC:
    1955             :       switch(ty)
    1956             :       {
    1957     2210446 :         case t_COMPLEX: return mulRc(x, y);
    1958      223973 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1959         343 :         case t_QUAD:  return mulRq(x, y);
    1960        1995 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1961             :       }
    1962             : 
    1963             :     case t_FFELT:
    1964          35 :       pari_err_TYPE2("*",x,y);
    1965             : 
    1966          21 :     case t_COMPLEX:
    1967             :       switch(ty)
    1968             :       {
    1969          14 :         case t_PADIC:
    1970          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1971           7 :         case t_QUAD:
    1972           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1973           7 :           return mulqf(y, x, lx);
    1974             :       }
    1975             : 
    1976             :     case t_PADIC: /* ty == t_QUAD */
    1977          28 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1978             :     }
    1979             :   }
    1980             : 
    1981    86816666 :   if (is_matvec_t(ty))
    1982             :   {
    1983     7506478 :     if (!is_matvec_t(tx))
    1984             :     {
    1985     7361582 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1986     7361582 :       z = cgetg_copy(y, &ly);
    1987   624940802 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1988     7361568 :       return z;
    1989             :     }
    1990      144896 :     switch(tx)
    1991             :     {
    1992       40652 :       case t_VEC:
    1993             :         switch(ty) {
    1994       14258 :           case t_COL: return RgV_RgC_mul(x,y);
    1995       26394 :           case t_MAT: return RgV_RgM_mul(x,y);
    1996             :         }
    1997           0 :         break;
    1998        1687 :       case t_COL:
    1999             :         switch(ty) {
    2000        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    2001           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    2002             :         }
    2003           0 :         break;
    2004      102557 :       case t_MAT:
    2005             :         switch(ty) {
    2006           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    2007      102557 :           case t_COL: return RgM_RgC_mul(x,y);
    2008             :         }
    2009             :     }
    2010             :   }
    2011    79325224 :   if (is_matvec_t(tx))
    2012             :   {
    2013      528199 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2014      528204 :     z = cgetg_copy(x, &lx);
    2015     2826915 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2016      527907 :     return z;
    2017             :   }
    2018    78797037 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2019             :   /* tx < ty, !ismatvec(x and y) */
    2020             : 
    2021    78797037 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2022     3214567 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2023    75582470 :   if (is_scalar_t(tx)) {
    2024    72293761 :     if (tx == t_POLMOD) {
    2025     3223268 :       vx = varn(gel(x,1));
    2026     3223268 :       vy = gvar(y);
    2027     3223268 :       if (vx != vy) {
    2028     2804248 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2029       46879 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2030             :       }
    2031             :       /* error if ty == t_SER */
    2032      419020 :       av = avma; y = gmod(y, gel(x,1));
    2033      419013 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2034             :     }
    2035    69070493 :     return mul_scal(y, x, ty);
    2036             :   }
    2037             : 
    2038             :   /* x and y are not scalars, nor matvec */
    2039     3288710 :   vx = gvar(x);
    2040     3288710 :   vy = gvar(y);
    2041     3288710 :   if (vx != vy) /* x or y is treated as a scalar */
    2042     2780040 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2043     2780040 :                                 : mul_scal(y, x, ty);
    2044             :   /* vx = vy */
    2045     1869927 :   switch(tx)
    2046             :   {
    2047     1869899 :     case t_POL:
    2048             :       switch (ty)
    2049             :       {
    2050        6622 :         case t_SER:
    2051             :         {
    2052             :           long vn;
    2053        6622 :           if (lg(x) == 2) return zeropol(vx);
    2054        6622 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2055        6391 :           av = avma;
    2056        6391 :           vn = RgX_valrem(x, &x);
    2057             :           /* take advantage of x = t^n ! */
    2058        6391 :           if (degpol(x)) {
    2059        1988 :             p1 = RgX_to_ser(x,lg(y));
    2060        1988 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2061        1988 :             p2 = gmul(p1,y);
    2062        1988 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2063             :           } else {
    2064        4403 :             set_avma(av);
    2065        4403 :             p2 = mul_ser_scal(y, gel(x,2));
    2066             :           }
    2067        6391 :           setvalp(p2, valp(p2) + vn);
    2068        6391 :           return p2;
    2069             :         }
    2070             : 
    2071     1863277 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2072             :       }
    2073           0 :       break;
    2074             : 
    2075          28 :     case t_SER:
    2076             :       switch (ty)
    2077             :       {
    2078          28 :         case t_RFRAC:
    2079          28 :           av = avma;
    2080          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2081             :       }
    2082           0 :       break;
    2083             :   }
    2084           0 :   pari_err_TYPE2("*",x,y);
    2085             :   return NULL; /* LCOV_EXCL_LINE */
    2086             : }
    2087             : 
    2088             : /* return a non-normalized result */
    2089             : GEN
    2090       93110 : sqr_ser_part(GEN x, long l1, long l2)
    2091             : {
    2092             :   long i, j, l;
    2093             :   pari_sp av;
    2094             :   GEN Z, z, p1, p2;
    2095             :   long mi;
    2096       93110 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2097       93096 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2098       93096 :   Z = cgetg(l2-l1+3, t_SER);
    2099       93096 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2100       93096 :   z = Z + 2-l1;
    2101       93096 :   x += 2; mi = 0;
    2102      348502 :   for (i=0; i<l1; i++)
    2103             :   {
    2104      255406 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2105             :   }
    2106             : 
    2107      624333 :   for (i=l1; i<=l2; i++)
    2108             :   {
    2109      531237 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2110      531237 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2111     1099313 :     for (j=i-mi; j<=minss(l,mi); j++)
    2112      568076 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2113      531237 :     p1 = gshift(p1,1);
    2114      531237 :     if ((i&1) == 0 && p2[i>>1])
    2115       57395 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2116      531237 :     gel(z,i) = gerepileupto(av,p1);
    2117             :   }
    2118       93096 :   return Z;
    2119             : }
    2120             : 
    2121             : GEN
    2122   984259532 : gsqr(GEN x)
    2123             : {
    2124             :   long i, lx;
    2125             :   pari_sp av, tetpil;
    2126             :   GEN z, p1, p2, p3, p4;
    2127             : 
    2128   984259532 :   switch(typ(x))
    2129             :   {
    2130   890058143 :     case t_INT: return sqri(x);
    2131    46874112 :     case t_REAL: return sqrr(x);
    2132      145476 :     case t_INTMOD: { GEN X = gel(x,1);
    2133      145476 :       z = cgetg(3,t_INTMOD);
    2134      145477 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2135      145478 :       gel(z,1) = icopy(X); return z;
    2136             :     }
    2137     1643217 :     case t_FRAC: return sqrfrac(x);
    2138             : 
    2139     3820458 :     case t_COMPLEX:
    2140     3820458 :       if (isintzero(gel(x,1))) {
    2141       51949 :         av = avma;
    2142       51949 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2143             :       }
    2144     3768509 :       z = cgetg(3,t_COMPLEX); av = avma;
    2145     3768509 :       p1 = gadd(gel(x,1),gel(x,2));
    2146     3768509 :       p2 = gsub(gel(x,1), gel(x,2));
    2147     3768509 :       p3 = gmul(gel(x,1),gel(x,2));
    2148     3768509 :       tetpil = avma;
    2149     3768509 :       gel(z,1) = gmul(p1,p2);
    2150     3768509 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2151             : 
    2152      186766 :     case t_PADIC:
    2153      186766 :       z = cgetg(5,t_PADIC);
    2154      186766 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2155      186766 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2156      186766 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2157      186766 :       gel(z,2) = icopy(gel(x,2));
    2158      186766 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2159      186766 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2160      186766 :       return z;
    2161             : 
    2162          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2163          28 :       p1 = gel(x,1);
    2164          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2165          28 :       p2 = gsqr(gel(x,2));
    2166          28 :       p3 = gsqr(gel(x,3));
    2167          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2168             : 
    2169          28 :       if (gequal0(gel(p1,3)))
    2170             :       {
    2171           7 :         tetpil = avma;
    2172           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2173           7 :         av = avma;
    2174           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2175           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2176             :       }
    2177             : 
    2178          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2179          21 :       tetpil = avma;
    2180          21 :       gel(z,2) = gadd(p2,p4);
    2181          21 :       gel(z,3) = gadd(p1,p3);
    2182          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2183             : 
    2184       43882 :     case t_POLMOD:
    2185       43882 :       return sqr_polmod(gel(x,1), gel(x,2));
    2186             : 
    2187     2894969 :     case t_FFELT: return FF_sqr(x);
    2188             : 
    2189      109364 :     case t_POL: return RgX_sqr(x);
    2190             : 
    2191       27426 :     case t_SER:
    2192       27426 :       lx = lg(x);
    2193       27426 :       if (ser_isexactzero(x)) {
    2194           7 :         GEN z = gcopy(x);
    2195           7 :         setvalp(z, 2*valp(x));
    2196           7 :         return z;
    2197             :       }
    2198       27419 :       if (lx < 40)
    2199       27174 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2200             :       else
    2201             :       {
    2202         245 :         pari_sp av = avma;
    2203         245 :         GEN z = cgetg(lx,t_SER);
    2204         245 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2205         245 :         x = ser2pol_i(x,lx);
    2206         245 :         x = RgXn_sqr(x, lx-2);
    2207         245 :         return gerepilecopy(av, fill_ser(z,x));
    2208             :       }
    2209             : 
    2210          14 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2211          14 :       gel(z,1) = gsqr(gel(x,1));
    2212          14 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2213             : 
    2214         980 :     case t_MAT: return RgM_sqr(x);
    2215          14 :     case t_QFR: return qfrsqr(x);
    2216    38454032 :     case t_QFI: return qfisqr(x);
    2217         658 :     case t_VECSMALL:
    2218         658 :       z = cgetg_copy(x, &lx);
    2219       16289 :       for (i=1; i<lx; i++)
    2220             :       {
    2221       15631 :         long xi = x[i];
    2222       15631 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2223       15631 :         z[i] = x[xi];
    2224             :       }
    2225         658 :       return z;
    2226             :   }
    2227           0 :   pari_err_TYPE2("*",x,x);
    2228             :   return NULL; /* LCOV_EXCL_LINE */
    2229             : }
    2230             : 
    2231             : /********************************************************************/
    2232             : /**                                                                **/
    2233             : /**                           DIVISION                             **/
    2234             : /**                                                                **/
    2235             : /********************************************************************/
    2236             : static GEN
    2237       30174 : div_rfrac_scal(GEN x, GEN y)
    2238             : {
    2239       30174 :   pari_sp av = avma;
    2240       30174 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2241       30174 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2242             : }
    2243             : static GEN
    2244       37452 : div_scal_rfrac(GEN x, GEN y)
    2245             : {
    2246       37452 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2247       37452 :   pari_sp av = avma;
    2248       37452 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2249             :   {
    2250           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2251           0 :     y1 = gel(y1,2);
    2252             :   }
    2253       37452 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2254             : }
    2255             : static GEN
    2256     1186544 : div_rfrac(GEN x, GEN y)
    2257     1186544 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2258             : 
    2259             : /* x != 0 */
    2260             : static GEN
    2261      870306 : div_ser_scal(GEN y, GEN x) {
    2262             :   long i, l;
    2263             :   GEN z;
    2264      870306 :   if (ser_isexactzero(y))
    2265           7 :     return scalarser(lg(y) == 2? Rg_get_0(x)
    2266           7 :                                : gdiv(gel(y,2), x), varn(y), valp(y));
    2267      870299 :   z = cgetg_copy(y, &l); z[1] = y[1];
    2268     4125705 :   for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
    2269      870299 :   return normalize(z);
    2270             : }
    2271             : GEN
    2272        7525 : ser_normalize(GEN x)
    2273             : {
    2274        7525 :   long i, lx = lg(x);
    2275             :   GEN c, z;
    2276        7525 :   if (lx == 2) return x;
    2277        7525 :   c = gel(x,2); if (gequal1(c)) return x;
    2278        7448 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2279      107499 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2280        7448 :   return z;
    2281             : }
    2282             : 
    2283             : /* y != 0 */
    2284             : static GEN
    2285     2522155 : div_T_scal(GEN x, GEN y, long tx) {
    2286     2522155 :   switch(tx)
    2287             :   {
    2288     1625546 :     case t_POL: return RgX_Rg_div(x, y);
    2289      870299 :     case t_SER: return div_ser_scal(x, y);
    2290       26310 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2291             :   }
    2292           0 :   pari_err_TYPE2("/",x,y);
    2293             :   return NULL; /* LCOV_EXCL_LINE */
    2294             : }
    2295             : 
    2296             : static GEN
    2297     9215554 : div_scal_pol(GEN x, GEN y) {
    2298     9215554 :   long ly = lg(y);
    2299             :   pari_sp av;
    2300     9215554 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2301     9137503 :   if (isrationalzero(x)) return zeropol(varn(y));
    2302     7090552 :   av = avma;
    2303     7090552 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2304             : }
    2305             : static GEN
    2306       18872 : div_scal_ser(GEN x, GEN y)
    2307       18872 : { pari_sp av = avma; return gerepileupto(av, gmul(x, ser_inv(y))); }
    2308             : static GEN
    2309     9224754 : div_scal_T(GEN x, GEN y, long ty) {
    2310     9224754 :   switch(ty)
    2311             :   {
    2312     9169193 :     case t_POL: return div_scal_pol(x, y);
    2313       18865 :     case t_SER: return div_scal_ser(x, y);
    2314       36696 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2315             :   }
    2316           0 :   pari_err_TYPE2("/",x,y);
    2317             :   return NULL; /* LCOV_EXCL_LINE */
    2318             : }
    2319             : 
    2320             : static GEN
    2321      900186 : ser2pol_ii(GEN x, long lx, long y1)
    2322             : {
    2323      900186 :   GEN y = cgetg(lx, t_POL);
    2324             :   long i;
    2325     3915637 :   for (i = lx-1; i > 1; i--) gel(y,i) = gel(x,i);
    2326      900186 :   y[1] = y1; return y;
    2327             : }
    2328             : 
    2329             : /* assume tx = ty = t_SER, same variable vx */
    2330             : static GEN
    2331      960120 : div_ser(GEN x, GEN y, long vx)
    2332             : {
    2333      960120 :   long v = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2334             :   GEN y_lead, z;
    2335      960120 :   pari_sp av = avma;
    2336             : 
    2337      960120 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2338      960113 :   if (ser_isexactzero(x))
    2339             :   {
    2340       59920 :     if (lx == 2) return zeroser(vx, v);
    2341         154 :     return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), v);
    2342             :   }
    2343      900193 :   y_lead = gel(y,2);
    2344      900193 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2345             :   {
    2346          14 :     GEN y0 = y;
    2347          14 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2348          35 :     for (v--, ly--,y++; ly > 2; v--, ly--, y++)
    2349             :     {
    2350          28 :       y_lead = gel(y,2);
    2351          28 :       if (!gequal0(y_lead)) break;
    2352             :     }
    2353          14 :     if (ly <= 2) pari_err_INV("div_ser", y0);
    2354             :   }
    2355      900186 :   if (ly < lx) lx = ly;
    2356      900186 :   z = cgetg(lx,t_SER); z[1] = evalvalp(v) | evalvarn(vx) | evalsigne(1);
    2357      900186 :   x = ser2pol_i(x, lx);
    2358      900186 :   y = ser2pol_ii(y, lx, z[1] & ~VALPBITS);
    2359      900186 :   if (lx == 3) gel(z,2) = gdiv(gel(x,2), gel(y,2));
    2360             :   else
    2361             :   { /* FIXME: use Karp/Markstein */
    2362      840418 :     y = RgXn_mul(x, RgXn_inv_i(y, lx-2), lx-2);
    2363      840418 :     z = fill_ser(z,y);
    2364             :   }
    2365      900186 :   return gerepilecopy(av, z);
    2366             : }
    2367             : /* x,y compatible PADIC */
    2368             : static GEN
    2369      428500 : divpp(GEN x, GEN y) {
    2370             :   pari_sp av;
    2371             :   long a, b;
    2372             :   GEN z, M;
    2373             : 
    2374      428500 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2375      428497 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2376      428196 :   a = precp(x);
    2377      428196 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2378      428196 :   z = cgetg(5, t_PADIC);
    2379      428263 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2380      428259 :   gel(z,2) = icopy(gel(x,2));
    2381      428263 :   gel(z,3) = icopy(M); av = avma;
    2382      428264 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2383      428259 :   return z;
    2384             : }
    2385             : static GEN
    2386       26061 : div_polmod_same(GEN T, GEN x, GEN y)
    2387             : {
    2388       26061 :   long v = varn(T);
    2389       26061 :   GEN a, z = cgetg(3, t_POLMOD);
    2390       26061 :   gel(z,1) = RgX_copy(T);
    2391       26061 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2392        9478 :     a = gdiv(x, y);
    2393       16583 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2394         119 :   {
    2395         119 :     pari_sp av = avma;
    2396         119 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2397             :   }
    2398       16464 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2399        4095 :   {
    2400        4095 :     pari_sp av = avma;
    2401        4095 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2402        4095 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2403        4095 :     a = gerepileupto(av, a);
    2404             :   }
    2405             :   else
    2406             :   {
    2407       12369 :     pari_sp av = avma;
    2408       12369 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2409       12369 :     a = gerepileupto(av, a);
    2410             :   }
    2411       26061 :   gel(z,2) = a; return z;
    2412             : }
    2413             : GEN
    2414   158664433 : gdiv(GEN x, GEN y)
    2415             : {
    2416   158664433 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2417             :   pari_sp av, tetpil;
    2418             :   GEN z, p1, p2;
    2419             : 
    2420   158664433 :   if (tx == ty) switch(tx)
    2421             :   {
    2422    50637553 :     case t_INT:
    2423    50637553 :       return Qdivii(x,y);
    2424             : 
    2425    46070433 :     case t_REAL: return divrr(x,y);
    2426       21029 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2427       21029 :       z = cgetg(3,t_INTMOD);
    2428       21029 :       if (X==Y || equalii(X,Y))
    2429       21022 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2430           7 :       gel(z,1) = gcdii(X,Y);
    2431           7 :       warn_coercion(X,Y,gel(z,1));
    2432           7 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2433           7 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2434             :     }
    2435     1667267 :     case t_FRAC: {
    2436     1667267 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2437     1667267 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2438     1667267 :       z = cgetg(3, t_FRAC);
    2439     1667267 :       p1 = gcdii(x1, y1);
    2440     1667267 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2441     1667267 :       p1 = gcdii(x2, y2);
    2442     1667267 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2443     1667267 :       tetpil = avma;
    2444     1667267 :       gel(z,2) = mulii(x2,y1);
    2445     1667267 :       gel(z,1) = mulii(x1,y2);
    2446     1667267 :       normalize_frac(z);
    2447     1667267 :       fix_frac_if_int_GC(z,tetpil);
    2448     1667267 :       return z;
    2449             :     }
    2450     3639053 :     case t_COMPLEX:
    2451     3639053 :       if (isintzero(gel(y,1)))
    2452             :       {
    2453       57883 :         y = gel(y,2);
    2454       57883 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2455       40992 :         z = cgetg(3,t_COMPLEX);
    2456       40992 :         gel(z,1) = gdiv(gel(x,2), y);
    2457       40992 :         av = avma;
    2458       40992 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2459       40992 :         return z;
    2460             :       }
    2461     3581170 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2462     3581170 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2463             : 
    2464      317506 :     case t_PADIC:
    2465      317506 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2466      317506 :       return divpp(x, y);
    2467             : 
    2468         245 :     case t_QUAD:
    2469         245 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2470         231 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2471         231 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2472             : 
    2473      236488 :     case t_FFELT: return FF_div(x,y);
    2474             : 
    2475       30436 :     case t_POLMOD:
    2476       30436 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2477       26061 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2478             :       else {
    2479        4375 :         av = avma;
    2480        4375 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2481             :       }
    2482       30436 :       return z;
    2483             : 
    2484    18442352 :     case t_POL:
    2485    18442352 :       vx = varn(x);
    2486    18442352 :       vy = varn(y);
    2487    18442352 :       if (vx != vy) {
    2488       95221 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2489       46361 :                             else return div_scal_pol(x, y);
    2490             :       }
    2491    18347131 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2492    18347131 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2493    18229368 :       av = avma;
    2494    18229368 :       return gerepileupto(av, gred_rfrac2(x,y));
    2495             : 
    2496      523055 :     case t_SER:
    2497      523055 :       vx = varn(x);
    2498      523055 :       vy = varn(y);
    2499      523055 :       if (vx != vy) {
    2500          21 :         if (varncmp(vx, vy) < 0)
    2501             :         {
    2502          14 :           if (!signe(y)) pari_err_INV("gdiv",y);
    2503           7 :           return div_ser_scal(x, y);
    2504             :         }
    2505           7 :         return div_scal_ser(x, y);
    2506             :       }
    2507      523034 :       return div_ser(x, y, vx);
    2508     1191024 :     case t_RFRAC:
    2509     1191024 :       vx = varn(gel(x,2));
    2510     1191024 :       vy = varn(gel(y,2));
    2511     1191024 :       if (vx != vy) {
    2512        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2513         756 :                             else return div_scal_rfrac(x, y);
    2514             :       }
    2515     1186544 :       return div_rfrac(x,y);
    2516             : 
    2517           7 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2518           7 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2519             : 
    2520          14 :     case t_MAT:
    2521          14 :       av = avma; p1 = RgM_inv(y);
    2522          14 :       if (!p1) pari_err_INV("gdiv",y);
    2523          14 :       return gerepileupto(av, RgM_mul(x, p1));
    2524             : 
    2525           0 :     default: pari_err_TYPE2("/",x,y);
    2526             :   }
    2527             : 
    2528    35888270 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2529             :   {
    2530     2459367 :     long s = signe(x);
    2531     2459367 :     if (!s) {
    2532      425868 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2533      425868 :       switch (ty)
    2534             :       {
    2535      422711 :       default: return gen_0;
    2536          21 :       case t_INTMOD:
    2537          21 :         z = cgetg(3,t_INTMOD);
    2538          21 :         gel(z,1) = icopy(gel(y,1));
    2539          21 :         gel(z,2) = gen_0; return z;
    2540        3136 :       case t_FFELT: return FF_zero(y);
    2541             :       }
    2542             :     }
    2543     2033499 :     if (is_pm1(x)) {
    2544      578401 :       if (s > 0) return ginv(y);
    2545       52169 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2546             :     }
    2547     1455098 :     switch(ty)
    2548             :     {
    2549      160890 :       case t_REAL: return divir(x,y);
    2550          21 :       case t_INTMOD:
    2551          21 :         z = cgetg(3, t_INTMOD);
    2552          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2553     1183608 :       case t_FRAC:
    2554     1183608 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2555     1183608 :         if (equali1(p1))
    2556             :         {
    2557      807753 :           set_avma((pari_sp)z);
    2558      807753 :           gel(z,2) = icopy(gel(y,1));
    2559      807753 :           gel(z,1) = mulii(gel(y,2), x);
    2560      807753 :           normalize_frac(z);
    2561      807753 :           fix_frac_if_int(z);
    2562             :         }
    2563             :         else
    2564             :         {
    2565      375855 :           x = diviiexact(x,p1); tetpil = avma;
    2566      375855 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2567      375855 :           gel(z,1) = mulii(gel(y,2), x);
    2568      375855 :           normalize_frac(z);
    2569      375855 :           fix_frac_if_int_GC(z,tetpil);
    2570             :         }
    2571     1183608 :         return z;
    2572             : 
    2573         266 :       case t_FFELT: return Z_FF_div(x,y);
    2574      109991 :       case t_COMPLEX: return divRc(x,y);
    2575          91 :       case t_PADIC: return divTp(x, y);
    2576         231 :       case t_QUAD:
    2577         231 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2578         231 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2579             :     }
    2580             :   }
    2581    33428904 :   if (gequal0(y))
    2582             :   {
    2583          49 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2584          28 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2585             :   }
    2586             : 
    2587    33429382 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2588             :   {
    2589    12587066 :     case t_REAL:
    2590    12587066 :       switch(ty)
    2591             :       {
    2592    11468811 :         case t_INT: return divri(x,y);
    2593      503599 :         case t_FRAC:
    2594      503599 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2595      503599 :           return gerepileuptoleaf(av, z);
    2596      614593 :         case t_COMPLEX: return divRc(x, y);
    2597          42 :         case t_QUAD: return divfq(x, y, realprec(x));
    2598          21 :         default: pari_err_TYPE2("/",x,y);
    2599             :       }
    2600             : 
    2601         287 :     case t_INTMOD:
    2602             :       switch(ty)
    2603             :       {
    2604         203 :         case t_INT:
    2605         203 :           z = cgetg(3, t_INTMOD);
    2606         203 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2607          28 :         case t_FRAC: { GEN X = gel(x,1);
    2608          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2609          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2610             :         }
    2611           7 :         case t_FFELT:
    2612           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2613           0 :             pari_err_OP("/",x,y);
    2614           7 :           return Z_FF_div(gel(x,2),y);
    2615             : 
    2616           0 :         case t_COMPLEX:
    2617           0 :           av = avma;
    2618           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2619             : 
    2620          14 :         case t_QUAD:
    2621          14 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2622          14 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2623             : 
    2624           7 :         case t_PADIC: { GEN X = gel(x,1);
    2625           7 :           z = cgetg(3, t_INTMOD);
    2626           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2627             :         }
    2628          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2629             :       }
    2630             : 
    2631             :     case t_FRAC:
    2632             :       switch(ty)
    2633             :       {
    2634      575055 :         case t_INT: z = cgetg(3, t_FRAC);
    2635      575055 :         p1 = gcdii(y,gel(x,1));
    2636      575055 :         if (equali1(p1))
    2637             :         {
    2638      286988 :           set_avma((pari_sp)z); tetpil = 0;
    2639      286988 :           gel(z,1) = icopy(gel(x,1));
    2640             :         }
    2641             :         else
    2642             :         {
    2643      288067 :           y = diviiexact(y,p1); tetpil = avma;
    2644      288067 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2645             :         }
    2646      575055 :         gel(z,2) = mulii(gel(x,2),y);
    2647      575055 :         normalize_frac(z);
    2648      575055 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2649      575055 :         return z;
    2650             : 
    2651       10115 :         case t_REAL:
    2652       10115 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2653       10115 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2654             : 
    2655           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2656           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2657           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2658             :         }
    2659             : 
    2660          28 :         case t_FFELT: av=avma;
    2661          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2662             : 
    2663         462 :         case t_COMPLEX: return divRc(x, y);
    2664             : 
    2665        2079 :         case t_PADIC:
    2666        2079 :           if (!signe(gel(x,1))) return gen_0;
    2667        2079 :           return divTp(x, y);
    2668             : 
    2669          42 :         case t_QUAD:
    2670          42 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2671          42 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2672             :       }
    2673             : 
    2674             :     case t_FFELT:
    2675         133 :       switch (ty)
    2676             :       {
    2677          49 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2678          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2679           7 :         case t_INTMOD:
    2680           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2681           0 :             pari_err_OP("/",x,y);
    2682           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2683          49 :         default:
    2684          49 :         pari_err_TYPE2("/",x,y);
    2685             :       }
    2686           0 :       break;
    2687             : 
    2688     5543430 :     case t_COMPLEX:
    2689             :       switch(ty)
    2690             :       {
    2691     5543430 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2692           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2693           0 :         case t_PADIC:
    2694           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2695           0 :         case t_QUAD:
    2696           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2697           0 :           return divfq(x, y, lx);
    2698             :       }
    2699             : 
    2700             :     case t_PADIC:
    2701             :       switch(ty)
    2702             :       {
    2703       11501 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2704       11410 :           return signe(gel(x,4))? divpT(x, y)
    2705       22911 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2706             :         }
    2707           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2708           7 :           z = cgetg(3, t_INTMOD);
    2709           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2710             :         }
    2711          14 :         case t_COMPLEX: case t_QUAD:
    2712          14 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2713          14 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2714             : 
    2715          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2716             :       }
    2717             : 
    2718             :     case t_QUAD:
    2719             :       switch (ty)
    2720             :       {
    2721        1197 :         case t_INT: case t_INTMOD: case t_FRAC:
    2722        1197 :           z = cgetg(4,t_QUAD);
    2723        1197 :           gel(z,1) = ZX_copy(gel(x,1));
    2724        1197 :           gel(z,2) = gdiv(gel(x,2), y);
    2725        1197 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2726          63 :         case t_REAL: return divqf(x, y, realprec(y));
    2727          14 :         case t_PADIC: return divTp(x, y);
    2728           0 :         case t_COMPLEX:
    2729           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2730           0 :           return divqf(x, y, ly);
    2731             :       }
    2732             :   }
    2733    14697853 :   switch(ty) {
    2734       41939 :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2735       41939 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2736          42 :     case t_MAT:
    2737          42 :       av = avma; p1 = RgM_inv(y);
    2738          35 :       if (!p1) pari_err_INV("gdiv",y);
    2739          28 :       return gerepileupto(av, gmul(x, p1));
    2740           0 :     case t_VEC: case t_COL:
    2741             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2742           0 :       pari_err_TYPE2("/",x,y);
    2743             :   }
    2744    14655872 :   switch(tx) {
    2745      955873 :     case t_VEC: case t_COL: case t_MAT:
    2746      955873 :       z = cgetg_copy(x, &lx);
    2747     4154471 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2748      955873 :       return z;
    2749           0 :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2750           0 :       pari_err_TYPE2("/",x,y);
    2751             :   }
    2752             : 
    2753    13699999 :   vy = gvar(y);
    2754    13699999 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2755      200435 :     vx = varn(X);
    2756      200435 :     if (vx != vy) {
    2757      200036 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2758      199707 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2759             :     }
    2760             :     /* y is POL, SER or RFRAC */
    2761         399 :     av = avma;
    2762         399 :     switch(ty)
    2763             :     {
    2764           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2765         399 :       default: y = ginvmod(gmod(y,X), X);
    2766             :     }
    2767         392 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2768             :   }
    2769             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2770             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2771             :    * variable NO_VARIABLE, then the operation is incorrect */
    2772    13499564 :   vx = gvar(x);
    2773    13499563 :   if (vx != vy) { /* includes cases where one is scalar */
    2774    11746584 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2775     9224425 :                         else return div_scal_T(x, y, ty);
    2776             :   }
    2777     1752979 :   switch(tx)
    2778             :   {
    2779     1046143 :     case t_POL:
    2780             :       switch(ty)
    2781             :       {
    2782          91 :         case t_SER:
    2783          91 :           if (lg(y) == 2)
    2784           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2785          91 :           p1 = RgX_to_ser(x,lg(y));
    2786          91 :           p2 = div_ser(p1, y, vx);
    2787          84 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2788          84 :           return p2;
    2789             : 
    2790     1046052 :         case t_RFRAC:
    2791             :         {
    2792     1046052 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2793     1046052 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2794           0 :             return mul_rfrac_scal(y2, y1, x);
    2795     1046052 :           av = avma;
    2796     1046052 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2797             :         }
    2798             :       }
    2799           0 :       break;
    2800             : 
    2801      437114 :     case t_SER:
    2802             :       switch(ty)
    2803             :       {
    2804      437100 :         case t_POL:
    2805      437100 :           if (lg(x) == 2)
    2806         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2807      436995 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2808      436995 :           p2 = div_ser(x, p1, vx);
    2809      436995 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2810      436995 :           return p2;
    2811          14 :         case t_RFRAC:
    2812          14 :           av = avma;
    2813          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2814             :       }
    2815           0 :       break;
    2816             : 
    2817      269701 :     case t_RFRAC:
    2818             :       switch(ty)
    2819             :       {
    2820      269701 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2821           0 :         case t_SER:
    2822           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2823           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2824             :       }
    2825           0 :       break;
    2826             :   }
    2827          21 :   pari_err_TYPE2("/",x,y);
    2828             :   return NULL; /* LCOV_EXCL_LINE */
    2829             : }
    2830             : 
    2831             : /********************************************************************/
    2832             : /**                                                                **/
    2833             : /**                     SIMPLE MULTIPLICATION                      **/
    2834             : /**                                                                **/
    2835             : /********************************************************************/
    2836             : GEN
    2837   155063653 : gmulsg(long s, GEN y)
    2838             : {
    2839             :   long ly, i;
    2840             :   pari_sp av;
    2841             :   GEN z;
    2842             : 
    2843   155063653 :   switch(typ(y))
    2844             :   {
    2845   102505586 :     case t_INT:  return mulsi(s,y);
    2846    42281943 :     case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
    2847      176451 :     case t_INTMOD: { GEN p = gel(y,1);
    2848      176451 :       z = cgetg(3,t_INTMOD);
    2849      176452 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2850      176453 :       gel(z,1) = icopy(p); return z;
    2851             :     }
    2852      627975 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2853     1161165 :     case t_FRAC:
    2854     1161165 :       if (!s) return gen_0;
    2855     1156594 :       z = cgetg(3,t_FRAC);
    2856     1156594 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2857     1156594 :       if (i == 1)
    2858             :       {
    2859      699397 :         gel(z,2) = icopy(gel(y,2));
    2860      699397 :         gel(z,1) = mulis(gel(y,1), s);
    2861             :       }
    2862             :       else
    2863             :       {
    2864      457197 :         gel(z,2) = divis(gel(y,2), i);
    2865      457197 :         gel(z,1) = mulis(gel(y,1), s/i);
    2866      457197 :         fix_frac_if_int(z);
    2867             :       }
    2868     1156594 :       return z;
    2869             : 
    2870     7075944 :     case t_COMPLEX:
    2871     7075944 :       if (!s) return gen_0;
    2872     7075930 :       z = cgetg(3, t_COMPLEX);
    2873     7075930 :       gel(z,1) = gmulsg(s,gel(y,1));
    2874     7075930 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2875             : 
    2876        1575 :     case t_PADIC:
    2877        1575 :       if (!s) return gen_0;
    2878        1575 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2879             : 
    2880           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2881           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2882           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2883           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2884             : 
    2885       78274 :     case t_POLMOD:
    2886       78274 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2887             : 
    2888      569785 :     case t_POL:
    2889      569785 :       if (!signe(y)) return RgX_copy(y);
    2890      554581 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2891      553125 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2892     1854177 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2893      553125 :       return normalizepol_lg(z, ly);
    2894             : 
    2895         182 :     case t_SER:
    2896         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2897         182 :       if (!s) return Rg_get_0(y);
    2898         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2899        3864 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2900         182 :       return normalize(z);
    2901             : 
    2902          14 :     case t_RFRAC:
    2903          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2904          14 :       if (s == 1) return gcopy(y);
    2905          14 :       if (s == -1) return gneg(y);
    2906          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2907             : 
    2908      589565 :     case t_VEC: case t_COL: case t_MAT:
    2909      589565 :       z = cgetg_copy(y, &ly);
    2910     1624050 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2911      589565 :       return z;
    2912             :   }
    2913           0 :   pari_err_TYPE("gmulsg",y);
    2914             :   return NULL; /* LCOV_EXCL_LINE */
    2915             : }
    2916             : 
    2917             : /********************************************************************/
    2918             : /**                                                                **/
    2919             : /**                       SIMPLE DIVISION                          **/
    2920             : /**                                                                **/
    2921             : /********************************************************************/
    2922             : 
    2923             : GEN
    2924    20293055 : gdivgs(GEN x, long s)
    2925             : {
    2926    20293055 :   long tx = typ(x), lx, i;
    2927             :   pari_sp av;
    2928             :   GEN z;
    2929             : 
    2930    20293055 :   if (!s)
    2931             :   {
    2932           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2933           0 :     pari_err_INV("gdivgs",gen_0);
    2934             :   }
    2935    20293063 :   switch(tx)
    2936             :   {
    2937     7982896 :     case t_INT: return Qdivis(x, s);
    2938     8556134 :     case t_REAL: return divrs(x,s);
    2939             : 
    2940       14651 :     case t_INTMOD:
    2941       14651 :       z = cgetg(3, t_INTMOD);
    2942       14651 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    2943             : 
    2944        1043 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    2945             : 
    2946      951238 :     case t_FRAC: z = cgetg(3, t_FRAC);
    2947      951238 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    2948      951238 :       if (i == 1)
    2949             :       {
    2950      627392 :         gel(z,2) = mulsi(s, gel(x,2));
    2951      627392 :         gel(z,1) = icopy(gel(x,1));
    2952             :       }
    2953             :       else
    2954             :       {
    2955      323846 :         gel(z,2) = mulsi(s/i, gel(x,2));
    2956      323846 :         gel(z,1) = divis(gel(x,1), i);
    2957             :       }
    2958      951238 :       normalize_frac(z);
    2959      951238 :       fix_frac_if_int(z); return z;
    2960             : 
    2961     2597933 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2962     2597933 :       gel(z,1) = gdivgs(gel(x,1),s);
    2963     2597933 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    2964             : 
    2965       98189 :     case t_PADIC: /* divpT */
    2966             :     {
    2967       98189 :       GEN p = gel(x,2);
    2968       98189 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    2969       97468 :       av = avma;
    2970       97468 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    2971             :     }
    2972             : 
    2973          28 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2974          28 :       gel(z,1) = ZX_copy(gel(x,1));
    2975          28 :       gel(z,2) = gdivgs(gel(x,2),s);
    2976          28 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    2977             : 
    2978       23135 :     case t_POLMOD:
    2979       23135 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    2980             : 
    2981         147 :     case t_RFRAC:
    2982         147 :       if (s == 1) return gcopy(x);
    2983         140 :       else if (s == -1) return gneg(x);
    2984         140 :       return div_rfrac_scal(x, stoi(s));
    2985             : 
    2986       66892 :     case t_POL: case t_SER:
    2987       66892 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    2988      328824 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2989       66892 :       return z;
    2990         777 :     case t_VEC: case t_COL: case t_MAT:
    2991         777 :       z = cgetg_copy(x, &lx);
    2992       39844 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2993         777 :       return z;
    2994             : 
    2995             :   }
    2996           0 :   pari_err_TYPE2("/",x, stoi(s));
    2997             :   return NULL; /* LCOV_EXCL_LINE */
    2998             : }
    2999             : 
    3000             : /* True shift (exact multiplication by 2^n) */
    3001             : GEN
    3002    23378267 : gmul2n(GEN x, long n)
    3003             : {
    3004             :   long lx, i, k, l;
    3005             :   GEN z, a, b;
    3006             : 
    3007    23378267 :   switch(typ(x))
    3008             :   {
    3009    12021758 :     case t_INT:
    3010    12021758 :       if (n>=0) return shifti(x,n);
    3011     2387370 :       if (!signe(x)) return gen_0;
    3012     2313425 :       l = vali(x); n = -n;
    3013     2313425 :       if (n<=l) return shifti(x,-n);
    3014      317063 :       z = cgetg(3,t_FRAC);
    3015      317063 :       gel(z,1) = shifti(x,-l);
    3016      317063 :       gel(z,2) = int2n(n-l); return z;
    3017             : 
    3018     7901278 :     case t_REAL:
    3019     7901278 :       return shiftr(x,n);
    3020             : 
    3021       76737 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3022       76737 :       z = cgetg(3,t_INTMOD);
    3023       76738 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3024       76598 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3025       76598 :       gel(z,1) = icopy(b); return z;
    3026             : 
    3027      216727 :     case t_FFELT: return FF_mul2n(x,n);
    3028             : 
    3029      603309 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3030      603309 :       l = vali(a);
    3031      603309 :       k = vali(b);
    3032      603309 :       if (n+l >= k)
    3033             :       {
    3034      170191 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3035      116109 :         l = n-k; k = -k;
    3036             :       }
    3037             :       else
    3038             :       {
    3039      433118 :         k = -(l+n); l = -l;
    3040             :       }
    3041      549227 :       z = cgetg(3,t_FRAC);
    3042      549227 :       gel(z,1) = shifti(a,l);
    3043      549227 :       gel(z,2) = shifti(b,k); return z;
    3044             : 
    3045      864469 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3046      864469 :       gel(z,1) = gmul2n(gel(x,1),n);
    3047      864468 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3048             : 
    3049         105 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3050         105 :       gel(z,1) = ZX_copy(gel(x,1));
    3051         105 :       gel(z,2) = gmul2n(gel(x,2),n);
    3052         105 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3053             : 
    3054       42994 :     case t_POLMOD:
    3055       42994 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3056             : 
    3057      478161 :     case t_POL:
    3058      478161 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3059     3177755 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3060      478161 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3061       83807 :     case t_SER:
    3062       83807 :       if (ser_isexactzero(x)) return gcopy(x);
    3063       83786 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3064      438256 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3065       83786 :       return normalize(z); /* needed if char = 2 */
    3066     1027039 :     case t_VEC: case t_COL: case t_MAT:
    3067     1027039 :       z = cgetg_copy(x, &lx);
    3068     4633019 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3069     1027039 :       return z;
    3070             : 
    3071          21 :     case t_RFRAC: /* int2n wrong if n < 0 */
    3072          21 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3073             : 
    3074       61867 :     case t_PADIC: /* int2n wrong if n < 0 */
    3075       61867 :       return gmul(gmul2n(gen_1,n),x);
    3076             :   }
    3077           0 :   pari_err_TYPE("gmul2n",x);
    3078             :   return NULL; /* LCOV_EXCL_LINE */
    3079             : }
    3080             : 
    3081             : /*******************************************************************/
    3082             : /*                                                                 */
    3083             : /*                              INVERSE                            */
    3084             : /*                                                                 */
    3085             : /*******************************************************************/
    3086             : static GEN
    3087       78851 : inv_polmod(GEN T, GEN x)
    3088             : {
    3089       78851 :   GEN z = cgetg(3,t_POLMOD), a;
    3090       78848 :   gel(z,1) = RgX_copy(T);
    3091       78848 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3092       62441 :     a = ginv(x);
    3093             :   else
    3094             :   {
    3095       16407 :     if (lg(T) == 5) /* quadratic fields */
    3096       12271 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3097             :     else
    3098        4136 :       a = RgXQ_inv(x, T);
    3099             :   }
    3100       78847 :   gel(z,2) = a; return z;
    3101             : }
    3102             : 
    3103             : GEN
    3104    18894566 : ginv(GEN x)
    3105             : {
    3106             :   long s;
    3107             :   pari_sp av, tetpil;
    3108             :   GEN z, y, p1, p2;
    3109             : 
    3110    18894566 :   switch(typ(x))
    3111             :   {
    3112     2377220 :     case t_INT:
    3113     2377220 :       if (is_pm1(x)) return icopy(x);
    3114     1941758 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3115     1941744 :       z = cgetg(3,t_FRAC);
    3116     1941744 :       gel(z,1) = s<0? gen_m1: gen_1;
    3117     1941744 :       gel(z,2) = absi(x); return z;
    3118             : 
    3119     2953751 :     case t_REAL: return invr(x);
    3120             : 
    3121        3948 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3122        3948 :       gel(z,1) = icopy(gel(x,1));
    3123        3948 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3124             : 
    3125      314778 :     case t_FRAC: {
    3126      314778 :       GEN a = gel(x,1), b = gel(x,2);
    3127      314778 :       s = signe(a);
    3128      314778 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3129      170541 :       z = cgetg(3,t_FRAC);
    3130      170541 :       gel(z,1) = icopy(b);
    3131      170541 :       gel(z,2) = icopy(a);
    3132      170541 :       normalize_frac(z); return z;
    3133             :     }
    3134     2975977 :     case t_COMPLEX:
    3135     2975977 :       av=avma;
    3136     2975977 :       p1=cxnorm(x);
    3137     2975977 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3138     2975977 :       tetpil=avma;
    3139     2975977 :       return gerepile(av,tetpil,divcR(p2,p1));
    3140             : 
    3141         273 :     case t_QUAD:
    3142         273 :       av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
    3143         273 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3144             : 
    3145        4522 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3146        4522 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3147        4515 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3148        4515 :       gel(z,2) = icopy(gel(x,2));
    3149        4515 :       gel(z,3) = icopy(gel(x,3));
    3150        4515 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3151             : 
    3152       78851 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3153       14308 :     case t_FFELT: return FF_inv(x);
    3154     5320846 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3155       15883 :     case t_SER: return ser_inv(x);
    3156        2912 :     case t_RFRAC:
    3157             :     {
    3158        2912 :       GEN n = gel(x,1), d = gel(x,2);
    3159        2912 :       pari_sp av = avma, ltop;
    3160        2912 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3161             : 
    3162        2912 :       n = simplify_shallow(n);
    3163        2912 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3164             :       {
    3165        2912 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3166         679 :         ltop = avma;
    3167         679 :         z = RgX_Rg_div(d,n);
    3168             :       } else {
    3169           0 :         ltop = avma;
    3170           0 :         z = cgetg(3,t_RFRAC);
    3171           0 :         gel(z,1) = RgX_copy(d);
    3172           0 :         gel(z,2) = RgX_copy(n);
    3173             :       }
    3174         679 :       stackdummy(av, ltop);
    3175         679 :       return z;
    3176             :     }
    3177             : 
    3178          14 :     case t_QFR:
    3179          14 :       av = avma; z = cgetg(5, t_QFR);
    3180          14 :       gel(z,1) = gel(x,1);
    3181          14 :       gel(z,2) = negi( gel(x,2) );
    3182          14 :       gel(z,3) = gel(x,3);
    3183          14 :       gel(z,4) = negr( gel(x,4) );
    3184          14 :       return gerepileupto(av, redreal(z));
    3185             : 
    3186     4810511 :     case t_QFI:
    3187     4810511 :       y = gcopy(x);
    3188     4810511 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3189     4695187 :         togglesign(gel(y,2));
    3190     4810511 :       return y;
    3191       20762 :     case t_MAT:
    3192       20762 :       y = RgM_inv(x);
    3193       20755 :       if (!y) pari_err_INV("ginv",x);
    3194       20685 :       return y;
    3195          28 :     case t_VECSMALL:
    3196             :     {
    3197          28 :       long i, lx = lg(x)-1;
    3198          28 :       y = zero_zv(lx);
    3199         112 :       for (i=1; i<=lx; i++)
    3200             :       {
    3201          84 :         long xi = x[i];
    3202          84 :         if (xi<1 || xi>lx || y[xi])
    3203           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3204          84 :         y[xi] = i;
    3205             :       }
    3206          28 :       return y;
    3207             :     }
    3208             :   }
    3209           0 :   pari_err_TYPE("inverse",x);
    3210             :   return NULL; /* LCOV_EXCL_LINE */
    3211             : }

Generated by: LCOV version 1.13