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.14.0 lcov report (development 27775-aca467eab2) Lines: 1910 2031 94.0 %
Date: 2022-07-03 07:33:15 Functions: 98 98 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                      GENERIC OPERATIONS                        **/
      18             : /**                         (first part)                           **/
      19             : /**                                                                **/
      20             : /********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_mod
      25             : 
      26             : /* assume z[1] was created last */
      27             : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
      28             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));
      29             : 
      30             : /* assume z[1] was created last */
      31             : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
      32             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));\
      33             : else\
      34             :   gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
      35             : 
      36             : static void
      37          98 : warn_coercion(GEN x, GEN y, GEN z)
      38             : {
      39          98 :   if (DEBUGLEVEL)
      40          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      41          98 : }
      42             : 
      43             : static long
      44          28 : kro_quad(GEN x, GEN y)
      45          28 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
      46             : 
      47             : /* is -1 not a square in Zp, assume p prime */
      48             : INLINE int
      49          42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      50             : 
      51             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      52             : static GEN mulpp(GEN x, GEN y);
      53             : static GEN divpp(GEN x, GEN y);
      54             : /* Argument codes for inline routines
      55             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      56             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      57             :  * T: some type (to be converted to PADIC)
      58             :  */
      59             : static GEN
      60   216795396 : addRc(GEN x, GEN y) {
      61   216795396 :   GEN z = cgetg(3,t_COMPLEX);
      62   216794613 :   gel(z,1) = gadd(x,gel(y,1));
      63   216740524 :   gel(z,2) = gcopy(gel(y,2)); return z;
      64             : }
      65             : static GEN
      66   266442132 : mulRc(GEN x, GEN y) {
      67   266442132 :   GEN z = cgetg(3,t_COMPLEX);
      68   266447362 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      69   266415934 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      70             : }
      71             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      72             : static GEN
      73          49 : mulRc_direct(GEN x, GEN y) {
      74          49 :   GEN z = cgetg(3,t_COMPLEX);
      75          49 :   gel(z,1) = gmul(x,gel(y,1));
      76          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      77             : }
      78             : static GEN
      79     2512625 : divRc(GEN x, GEN y) {
      80     2512625 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      81     2512624 :   GEN z = cgetg(3,t_COMPLEX);
      82     2512624 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      83     2512600 :   gel(z,2) = gmul(mt, gel(y,2));
      84     2512601 :   return z;
      85             : }
      86             : static GEN
      87    18609208 : divcR(GEN x, GEN y) {
      88    18609208 :   GEN z = cgetg(3,t_COMPLEX);
      89    18609124 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      90    18608188 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      91             : }
      92             : static GEN
      93        1267 : addRq(GEN x, GEN y) {
      94        1267 :   GEN z = cgetg(4,t_QUAD);
      95        1267 :   gel(z,1) = ZX_copy(gel(y,1));
      96        1267 :   gel(z,2) = gadd(x, gel(y,2));
      97        1267 :   gel(z,3) = gcopy(gel(y,3)); return z;
      98             : }
      99             : static GEN
     100        2352 : mulRq(GEN x, GEN y) {
     101        2352 :   GEN z = cgetg(4,t_QUAD);
     102        2352 :   gel(z,1) = ZX_copy(gel(y,1));
     103        2352 :   gel(z,2) = gmul(x,gel(y,2));
     104        2352 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     105             : }
     106             : static GEN
     107          70 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     108          70 :   long i = gexpo(x) - gexpo(y);
     109          70 :   if (i > 0) prec += nbits2extraprec(i);
     110          70 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     111             : }
     112             : static GEN
     113    29262851 : mulrfrac(GEN x, GEN y)
     114             : {
     115             :   pari_sp av;
     116    29262851 :   GEN z, a = gel(y,1), b = gel(y,2);
     117    29262851 :   if (is_pm1(a)) /* frequent special case */
     118             :   {
     119     5242059 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     120     5241913 :     return z;
     121             :   }
     122    24020800 :   av = avma;
     123    24020800 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     124             : }
     125             : static GEN
     126          28 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     127          28 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     128             : }
     129             : static GEN
     130          63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     131          63 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     132             : }
     133             : static GEN
     134          42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     135          42 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     136             : }
     137             : /* y PADIC, x + y by converting x to padic */
     138             : static GEN
     139           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     140             : 
     141           7 :   if (!valp(y)) z = cvtop2(x,y);
     142             :   else {
     143           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     144           7 :     z  = cvtop(x, gel(y,2), l);
     145             :   }
     146           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     147             : }
     148             : /* y PADIC, x * y by converting x to padic */
     149             : static GEN
     150     6496826 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     151     6496826 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     152             : }
     153             : /* y PADIC, non zero x / y by converting x to padic */
     154             : static GEN
     155        3647 : divTp(GEN x, GEN y) { pari_sp av = avma;
     156        3647 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     157             : }
     158             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     159             :  * converted to O(p^e) and division by 0 */
     160             : static GEN
     161     1327200 : divpT(GEN x, GEN y) { pari_sp av = avma;
     162     1327200 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     163             : }
     164             : 
     165             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     166             :  * clean memory from z on */
     167             : static GEN
     168     1974870 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     169     1974870 :   if (lgefint(X) == 3) {
     170     1894040 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     171     1894040 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     172             :   }
     173             :   else {
     174       80830 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     175       80829 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     176             :   }
     177     1974871 :   gel(z,1) = icopy(X); return z;
     178             : }
     179             : static GEN
     180     1153254 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     181     1153254 :   if (lgefint(X) == 3) {
     182      779220 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     183      779220 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     184             :   }
     185             :   else {
     186      374034 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     187      374034 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     188             :   }
     189     1153254 :   gel(z,1) = icopy(X); return z;
     190             : }
     191             : /* cf add_intmod_same */
     192             : static GEN
     193     3218325 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     194     3218325 :   if (lgefint(X) == 3) {
     195     3109932 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     196     3109932 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     197             :   }
     198             :   else
     199      108393 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     200     3218334 :   gel(z,1) = icopy(X); return z;
     201             : }
     202             : /* cf add_intmod_same */
     203             : static GEN
     204      337532 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     205             : {
     206      337532 :   if (lgefint(X) == 3) {
     207      314719 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     208      314712 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     209             :   }
     210             :   else
     211       22813 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     212      337526 :   gel(z,1) = icopy(X); return z;
     213             : }
     214             : 
     215             : /*******************************************************************/
     216             : /*                                                                 */
     217             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     218             : /*                                                                 */
     219             : /* (static routines are not memory clean, but OK for gerepileupto) */
     220             : /*******************************************************************/
     221             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     222             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     223             : static GEN
     224     9884933 : rfrac_denom_mul_scal(GEN d, GEN y)
     225             : {
     226     9884933 :   GEN D = RgX_Rg_mul(d, y);
     227     9884933 :   if (lg(D) != lg(d))
     228             :   { /* try to generate a meaningful diagnostic */
     229           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     230           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     231             :   }
     232     9884933 :   return D;
     233             : }
     234             : 
     235             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     236             : GEN
     237    57837375 : gred_rfrac_simple(GEN n, GEN d)
     238             : {
     239             :   GEN c, cn, cd, z;
     240    57837375 :   long dd = degpol(d);
     241             : 
     242    57837375 :   if (dd <= 0)
     243             :   {
     244        4123 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     245        4123 :     n = gdiv(n, gel(d,2));
     246        4123 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     247        4123 :     return n;
     248             :   }
     249             : 
     250    57833252 :   cd = content(d);
     251    59699419 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     252    57833252 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     253    57833252 :   if (!gequal1(cd)) {
     254     6582165 :     d = RgX_Rg_div(d,cd);
     255     6582165 :     if (!gequal1(cn))
     256             :     {
     257     1309528 :       if (gequal0(cn)) {
     258          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     259           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     260           0 :         c = gen_1;
     261             :       } else {
     262     1309479 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     263     1309479 :         c = gdiv(cn,cd);
     264             :       }
     265             :     }
     266             :     else
     267     5272637 :       c = ginv(cd);
     268             :   } else {
     269    51251087 :     if (!gequal1(cn))
     270             :     {
     271     3272895 :       if (gequal0(cn)) {
     272         931 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     273          21 :         c = gen_1;
     274             :       } else {
     275     3271964 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     276     3271964 :         c = cn;
     277             :       }
     278             :     } else {
     279    47978192 :       GEN y = cgetg(3,t_RFRAC);
     280    47978192 :       gel(y,1) = gcopy(n);
     281    47978192 :       gel(y,2) = RgX_copy(d); return y;
     282             :     }
     283             :   }
     284             : 
     285     9854101 :   if (typ(c) == t_POL)
     286             :   {
     287      904069 :     z = c;
     288      943213 :     do { z = content(z); } while (typ(z) == t_POL);
     289      904069 :     cd = denom_i(z);
     290      904069 :     cn = gmul(c, cd);
     291             :   }
     292             :   else
     293             :   {
     294     8950032 :     cn = numer_i(c);
     295     8950032 :     cd = denom_i(c);
     296             :   }
     297     9854101 :   z = cgetg(3,t_RFRAC);
     298     9854101 :   gel(z,1) = gmul(n, cn);
     299     9854101 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     300     9854101 :   return z;
     301             : }
     302             : 
     303             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     304             : static GEN
     305      249123 : fix_rfrac(GEN x, long d)
     306             : {
     307             :   GEN z, N, D;
     308      249123 :   if (!d || typ(x) == t_POL) return x;
     309      175511 :   z = cgetg(3, t_RFRAC);
     310      175511 :   N = gel(x,1);
     311      175511 :   D = gel(x,2);
     312      175511 :   if (d > 0) {
     313       11690 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     314      177842 :                                                    : monomialcopy(N,d,varn(D));
     315      166152 :     gel(z, 2) = RgX_copy(D);
     316             :   } else {
     317        9359 :     gel(z, 1) = gcopy(N);
     318        9359 :     gel(z, 2) = RgX_shift(D, -d);
     319             :   }
     320      175511 :   return z;
     321             : }
     322             : 
     323             : /* assume d != 0 */
     324             : static GEN
     325    44767912 : gred_rfrac2(GEN n, GEN d)
     326             : {
     327             :   GEN y, z;
     328             :   long v, vd, vn;
     329             : 
     330    44767912 :   n = simplify_shallow(n);
     331    44767912 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     332    37936360 :   d = simplify_shallow(d);
     333    37936360 :   if (typ(d) != t_POL) return gdiv(n,d);
     334    36756564 :   vd = varn(d);
     335    36756564 :   if (typ(n) != t_POL)
     336             :   {
     337    20622506 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     338    20621099 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     339           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     340             :   }
     341    16134058 :   vn = varn(n);
     342    16134058 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     343    15996412 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     344             : 
     345             :   /* now n and d are t_POLs in the same variable */
     346    15836660 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     347    15836660 :   if (!degpol(d))
     348             :   {
     349    12760350 :     n = RgX_Rg_div(n,gel(d,2));
     350    12760350 :     return v? RgX_mulXn(n,v): n;
     351             :   }
     352             : 
     353             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     354     3076310 :   if (!isinexact(n) && !isinexact(d))
     355             :   {
     356     3076093 :     y = RgX_divrem(n, d, &z);
     357     3076093 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     358      248906 :     z = RgX_gcd(d, z);
     359      248906 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     360             :   }
     361      249123 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     362             : }
     363             : 
     364             : /* x,y t_INT, return x/y in reduced form */
     365             : GEN
     366   109262585 : Qdivii(GEN x, GEN y)
     367             : {
     368   109262585 :   pari_sp av = avma;
     369             :   GEN r, q;
     370             : 
     371   109262585 :   if (lgefint(y) == 3)
     372             :   {
     373    94169777 :     q = Qdiviu(x, y[2]);
     374    94166229 :     if (signe(y) > 0) return q;
     375    10672004 :     if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
     376    10672184 :     return q;
     377             :   }
     378    15092808 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     379    15094103 :   if (equali1(x))
     380             :   {
     381     5177462 :     if (!signe(y)) pari_err_INV("gdiv",y);
     382     5177350 :     retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
     383             :   }
     384     9916649 :   q = dvmdii(x,y,&r);
     385     9916402 :   if (r == gen_0) return q; /* gen_0 intended */
     386     7012052 :   r = gcdii(y, r);
     387     7011885 :   if (lgefint(r) == 3)
     388             :   {
     389     6184303 :     ulong t = r[2];
     390     6184303 :     set_avma(av);
     391     6184388 :     if (t == 1) q = mkfraccopy(x,y);
     392             :     else
     393             :     {
     394     2463879 :       q = cgetg(3,t_FRAC);
     395     2464262 :       gel(q,1) = diviuexact(x,t);
     396     2463957 :       gel(q,2) = diviuexact(y,t);
     397             :     }
     398             :   }
     399             :   else
     400             :   { /* rare: r and q left on stack for efficiency */
     401      827582 :     q = cgetg(3,t_FRAC);
     402      827620 :     gel(q,1) = diviiexact(x,r);
     403      827626 :     gel(q,2) = diviiexact(y,r);
     404             :   }
     405     7011813 :   normalize_frac(q); return q;
     406             : }
     407             : 
     408             : /* x t_INT, return x/y in reduced form */
     409             : GEN
     410   112342243 : Qdiviu(GEN x, ulong y)
     411             : {
     412   112342243 :   pari_sp av = avma;
     413             :   ulong r, t;
     414             :   GEN q;
     415             : 
     416   112342243 :   if (y == 1) return icopy(x);
     417    93151721 :   if (!y) pari_err_INV("gdiv",gen_0);
     418    93151721 :   if (equali1(x)) retmkfrac(gen_1, utoipos(y));
     419    87468311 :   q = absdiviu_rem(x,y,&r);
     420    87460479 :   if (!r)
     421             :   {
     422    48985498 :     if (signe(x) < 0) togglesign(q);
     423    48986523 :     return q;
     424             :   }
     425    38474981 :   t = ugcd(y, r); set_avma(av);
     426    38481479 :   if (t == 1) retmkfrac(icopy(x), utoipos(y));
     427    15758744 :   retmkfrac(diviuexact(x,t), utoipos(y / t));
     428             : }
     429             : 
     430             : /* x t_INT, return x/y in reduced form */
     431             : GEN
     432     1326121 : Qdivis(GEN x, long y)
     433             : {
     434     1326121 :   pari_sp av = avma;
     435             :   ulong r, t;
     436             :   long s;
     437             :   GEN q;
     438             : 
     439     1326121 :   if (y > 0) return Qdiviu(x, y);
     440       98434 :   if (!y) pari_err_INV("gdiv",gen_0);
     441       98434 :   s = signe(x);
     442       98434 :   if (!s) return gen_0;
     443       70140 :   if (y < 0) { y = -y; s = -s; }
     444       70140 :   if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
     445       69860 :   if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
     446       67830 :   q = absdiviu_rem(x,y,&r);
     447       67830 :   if (!r)
     448             :   {
     449       38402 :     if (s < 0) togglesign(q);
     450       38402 :     return q;
     451             :   }
     452       29428 :   t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
     453       29428 :   if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
     454       29428 :   gel(q,1) = x; setsigne(x, s);
     455       29428 :   gel(q,2) = utoipos(y); return q;
     456             : }
     457             : 
     458             : /*******************************************************************/
     459             : /*                                                                 */
     460             : /*                          CONJUGATION                            */
     461             : /*                                                                 */
     462             : /*******************************************************************/
     463             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     464             : static GEN
     465       17850 : quad_polmod_conj(GEN x, GEN y)
     466             : {
     467             :   GEN z, u, v, a, b;
     468       17850 :   if (typ(x) != t_POL) return gcopy(x);
     469       17850 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     470       17850 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     471       17850 :   b = gel(y,3); v = gel(x,2);
     472       17850 :   z = cgetg(4, t_POL); z[1] = x[1];
     473       17850 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     474       17850 :   gel(z,3) = gneg(u); return z;
     475             : }
     476             : static GEN
     477       17850 : quad_polmod_norm(GEN x, GEN y)
     478             : {
     479             :   GEN z, u, v, a, b, c;
     480       17850 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     481       17850 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     482       17850 :   b = gel(y,3); v = gel(x,2);
     483       17850 :   c = gel(y,2);
     484       17850 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     485       17850 :   if (!gequal1(a)) z = gdiv(z, a);
     486       17850 :   return gadd(z, gsqr(v));
     487             : }
     488             : 
     489             : GEN
     490    17163436 : conj_i(GEN x)
     491             : {
     492             :   long lx, i;
     493             :   GEN y;
     494             : 
     495    17163436 :   switch(typ(x))
     496             :   {
     497     5783852 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     498     5783852 :       return x;
     499             : 
     500    11220800 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     501             : 
     502         924 :     case t_QUAD:
     503         924 :       y = cgetg(4,t_QUAD);
     504         924 :       gel(y,1) = gel(x,1);
     505         924 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     506         924 :                                       : gadd(gel(x,2), gel(x,3));
     507         924 :       gel(y,3) = gneg(gel(x,3)); return y;
     508             : 
     509        8708 :     case t_POL: case t_SER:
     510        8708 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     511       29323 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     512        8708 :       return y;
     513             : 
     514      149184 :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     515      149184 :       y = cgetg_copy(x, &lx);
     516      538094 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     517      149183 :       return y;
     518             : 
     519           0 :     case t_POLMOD:
     520             :     {
     521           0 :       GEN X = gel(x,1);
     522           0 :       long d = degpol(X);
     523           0 :       if (d < 2) return x;
     524           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     525             :     }
     526             :   }
     527           0 :   pari_err_TYPE("gconj",x);
     528             :   return NULL; /* LCOV_EXCL_LINE */
     529             : }
     530             : GEN
     531       95388 : gconj(GEN x)
     532       95388 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     533             : 
     534             : GEN
     535          84 : conjvec(GEN x,long prec)
     536             : {
     537             :   long lx, s, i;
     538             :   GEN z;
     539             : 
     540          84 :   switch(typ(x))
     541             :   {
     542           0 :     case t_INT: case t_INTMOD: case t_FRAC:
     543           0 :       return mkcolcopy(x);
     544             : 
     545           0 :     case t_COMPLEX: case t_QUAD:
     546           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     547             : 
     548          28 :     case t_FFELT:
     549          28 :       return FF_conjvec(x);
     550             : 
     551           0 :     case t_VEC: case t_COL:
     552           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     553           0 :       if (lx == 1) return z;
     554           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     555           0 :       s = lgcols(z);
     556           0 :       for (i=2; i<lx; i++)
     557             :       {
     558           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     559           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     560             :       }
     561           0 :       break;
     562             : 
     563          56 :     case t_POLMOD: {
     564          56 :       GEN T = gel(x,1), r;
     565             :       pari_sp av;
     566             : 
     567          56 :       lx = lg(T);
     568          56 :       if (lx <= 3) return cgetg(1,t_COL);
     569          56 :       x = gel(x,2);
     570         238 :       for (i=2; i<lx; i++)
     571             :       {
     572         189 :         GEN c = gel(T,i);
     573         189 :         switch(typ(c)) {
     574           7 :           case t_INTMOD: {
     575           7 :             GEN p = gel(c,1);
     576             :             pari_sp av;
     577           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     578           7 :             av = avma;
     579           7 :             T = RgX_to_FpX(T,p);
     580           7 :             x = RgX_to_FpX(x, p);
     581           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     582           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     583           7 :             return gerepileupto(av, z);
     584             :           }
     585         182 :           case t_INT:
     586         182 :           case t_FRAC: break;
     587           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     588             :         }
     589             :       }
     590          49 :       if (typ(x) != t_POL)
     591             :       {
     592           0 :         if (!is_rational_t(typ(x)))
     593           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     594           0 :         retconst_col(lx-3, gcopy(x));
     595             :       }
     596          49 :       RgX_check_QX(x,"conjvec");
     597          49 :       av = avma;
     598          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     599          49 :       r = cleanroots(T,prec);
     600          49 :       z = cgetg(lx-2,t_COL);
     601         182 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     602          49 :       return gerepileupto(av, z);
     603             :     }
     604             : 
     605           0 :     default:
     606           0 :       pari_err_TYPE("conjvec",x);
     607             :       return NULL; /* LCOV_EXCL_LINE */
     608             :   }
     609           0 :   return z;
     610             : }
     611             : 
     612             : /********************************************************************/
     613             : /**                                                                **/
     614             : /**                           ADDITION                             **/
     615             : /**                                                                **/
     616             : /********************************************************************/
     617             : /* x, y compatible PADIC, op = add or sub */
     618             : static GEN
     619    10531647 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     620             : {
     621    10531647 :   pari_sp av = avma;
     622             :   long d,e,r,rx,ry;
     623    10531647 :   GEN u, z, mod, p = gel(x,2);
     624             :   int swap;
     625             : 
     626    10531647 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     627    10531385 :   e = valp(x);
     628    10531385 :   r = valp(y); d = r-e;
     629    10531385 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     630    10531385 :   rx = precp(x);
     631    10531385 :   ry = precp(y);
     632    10531385 :   if (d) /* v(x) < v(y) */
     633             :   {
     634     7368739 :     r = d+ry; z = powiu(p,d);
     635     7368823 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     636     7368840 :     z = mulii(z,gel(y,4));
     637     7368787 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     638             :   }
     639             :   else
     640             :   {
     641             :     long c;
     642     3162646 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     643     3162646 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     644     3163446 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     645             :     {
     646      138619 :       set_avma(av); return zeropadic(p, e+r);
     647             :     }
     648     3025098 :     if (c)
     649             :     {
     650     1870118 :       mod = diviiexact(mod, powiu(p,c));
     651     1870115 :       r -= c;
     652     1870115 :       e += c;
     653             :     }
     654             :   }
     655    10393850 :   u = modii(u, mod);
     656    10392873 :   set_avma(av); z = cgetg(5,t_PADIC);
     657    10392210 :   z[1] = evalprecp(r) | evalvalp(e);
     658    10392160 :   gel(z,2) = icopy(p);
     659    10391998 :   gel(z,3) = icopy(mod);
     660    10392057 :   gel(z,4) = icopy(u); return z;
     661             : }
     662             : /* Rg_to_Fp(t_FRAC) without GC */
     663             : static GEN
     664       23555 : Q_to_Fp(GEN x, GEN p)
     665       23555 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     666             : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
     667             : static GEN
     668     3513464 : addQp(GEN x, GEN y)
     669             : {
     670     3513464 :   pari_sp av = avma;
     671     3513464 :   long d, r, e, vy = valp(y), py = precp(y);
     672     3513464 :   GEN z, q, mod, u, p = gel(y,2);
     673             : 
     674     3513464 :   e = Q_pvalrem(x, p, &x);
     675     3513450 :   d = vy - e; r = d + py;
     676     3513450 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     677     3511497 :   mod = gel(y,3);
     678     3511497 :   u   = gel(y,4);
     679     3511497 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     680             : 
     681     3511491 :   if (d > 0)
     682             :   {
     683      893989 :     q = powiu(p,d);
     684      894075 :     mod = mulii(mod, q);
     685      894075 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     686      894075 :     u = addii(x,  mulii(u, q));
     687             :   }
     688     2617502 :   else if (d < 0)
     689             :   {
     690      397578 :     q = powiu(p,-d);
     691      397578 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     692      397578 :     u = addii(u, mulii(x, q));
     693      397578 :     r = py; e = vy;
     694             :   }
     695             :   else
     696             :   {
     697             :     long c;
     698     2219924 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     699     2219924 :     u = addii(u, x);
     700     2219909 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     701             :     {
     702         904 :       set_avma(av); return zeropadic(p,e+r);
     703             :     }
     704     2219012 :     if (c)
     705             :     {
     706      796797 :       mod = diviiexact(mod, powiu(p,c));
     707      796797 :       r -= c;
     708      796797 :       e += c;
     709             :     }
     710             :   }
     711     3510668 :   u = modii(u, mod); set_avma(av);
     712     3510609 :   z = cgetg(5,t_PADIC);
     713     3510560 :   z[1] = evalprecp(r) | evalvalp(e);
     714     3510555 :   gel(z,2) = icopy(p);
     715     3510536 :   gel(z,3) = icopy(mod);
     716     3510534 :   gel(z,4) = icopy(u); return z;
     717             : }
     718             : 
     719             : /* Mod(x,X) + Mod(y,X) */
     720             : #define addsub_polmod_same addsub_polmod_scal
     721             : /* Mod(x,X) +/- Mod(y,Y) */
     722             : static GEN
     723        2142 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     724             : {
     725        2142 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     726        2142 :   GEN z = cgetg(3,t_POLMOD);
     727        2142 :   long vx = varn(X), vy = varn(Y);
     728        2142 :   if (vx==vy) {
     729             :     pari_sp av;
     730          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     731          14 :     warn_coercion(X,Y,gel(z,1));
     732          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     733             :   }
     734        2128 :   if (varncmp(vx, vy) < 0)
     735        2128 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     736             :   else
     737           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     738        2128 :   gel(z,2) = op(x, y); return z;
     739             : }
     740             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     741             : static GEN
     742    12527896 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     743    12527896 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     744             : 
     745             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     746             : static GEN
     747      369328 : add_ser_scal(GEN y, GEN x)
     748             : {
     749             :   long i, v, ly, vy;
     750             :   GEN z;
     751             : 
     752      369328 :   if (isrationalzero(x)) return gcopy(y);
     753      344275 :   ly = lg(y);
     754      344275 :   v = valp(y);
     755      344275 :   if (v < 3-ly) return gcopy(y);
     756             :   /* v + ly >= 3 */
     757      344051 :   if (v < 0)
     758             :   {
     759        1169 :     z = cgetg(ly,t_SER); z[1] = y[1];
     760        3304 :     for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
     761        1169 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     762        3276 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     763        1169 :     return normalizeser(z);
     764             :   }
     765      342882 :   vy = varn(y);
     766      342882 :   if (v > 0)
     767             :   {
     768       15946 :     if (ser_isexactzero(y))
     769        7406 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
     770        8540 :     y -= v; ly += v;
     771        8540 :     z = cgetg(ly,t_SER);
     772        8540 :     x = gcopy(x);
     773       19138 :     for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
     774             :   }
     775      326936 :   else if (ser_isexactzero(y)) return gcopy(y);
     776             :   else
     777             :   { /* v = 0, ly >= 3 */
     778      326929 :     z = cgetg(ly,t_SER);
     779      326929 :     x = gadd(x, gel(y,2));
     780      326929 :     i = 3;
     781             :   }
     782     1508196 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     783      335469 :   gel(z,2) = x;
     784      335469 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     785      335469 :   return gequal0(x)? normalizeser(z): z;
     786             : }
     787             : static long
     788     6921726 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     789             : /* x,y t_SER in the same variable: x+y */
     790             : static GEN
     791     3461262 : ser_add(GEN x, GEN y)
     792             : {
     793     3461262 :   long i, lx,ly, n = valp(y) - valp(x);
     794             :   GEN z;
     795     3461262 :   if (n < 0) { n = -n; swap(x,y); }
     796             :   /* valp(x) <= valp(y) */
     797     3461262 :   lx = _serprec(x);
     798     3461262 :   if (lx == 2) /* don't lose type information */
     799             :   {
     800         798 :     z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
     801         798 :     setvalp(z, valp(x)); return z;
     802             :   }
     803     3460464 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     804     3460464 :   if (n)
     805             :   {
     806      114887 :     if (n+2 > lx) return gcopy(x);
     807      114313 :     z = cgetg(ly,t_SER);
     808      841684 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     809      529248 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     810             :   } else {
     811     3345577 :     z = cgetg(ly,t_SER);
     812    16837374 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     813             :   }
     814     3459890 :   z[1] = x[1]; return normalizeser(z);
     815             : }
     816             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     817             : static GEN
     818     8824491 : add_rfrac_scal(GEN y, GEN x)
     819             : {
     820             :   pari_sp av;
     821             :   GEN n;
     822             : 
     823     8824491 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     824     5150726 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     825     5150726 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     826             : }
     827             : 
     828             : /* x "scalar", ty != t_MAT and nonscalar */
     829             : static GEN
     830    21229791 : add_scal(GEN y, GEN x, long ty)
     831             : {
     832    21229791 :   switch(ty)
     833             :   {
     834    16377606 :     case t_POL: return RgX_Rg_add(y, x);
     835      369307 :     case t_SER: return add_ser_scal(y, x);
     836     4468757 :     case t_RFRAC: return add_rfrac_scal(y, x);
     837           0 :     case t_COL: return RgC_Rg_add(y, x);
     838       14120 :     case t_VEC:
     839       14120 :       if (isintzero(x)) return gcopy(y);
     840         182 :       break;
     841             :   }
     842         183 :   pari_err_TYPE2("+",x,y);
     843             :   return NULL; /* LCOV_EXCL_LINE */
     844             : }
     845             : 
     846             : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
     847             : static GEN
     848    12571103 : setfrac(GEN z, GEN a, GEN b)
     849             : {
     850    12571103 :   gel(z,1) = icopy_avma(a, (pari_sp)z);
     851    12571131 :   gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
     852    12571209 :   set_avma((pari_sp)gel(z,2)); return z;
     853             : }
     854             : /* z <- a / (b*Q), (Q,a) = 1 */
     855             : static GEN
     856    11829893 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
     857             : {
     858    11829893 :   GEN q = Qdivii(a, b); /* != 0 */
     859    11829989 :   if (typ(q) == t_INT)
     860             :   {
     861     1538905 :     gel(z,1) = gerepileuptoint((pari_sp)Q, q);
     862     1538905 :     gel(z,2) = Q; return z;
     863             :   }
     864    10291084 :   return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
     865             : }
     866             : static GEN
     867    24655885 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     868             : {
     869    24655885 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     870    24655885 :   GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
     871    24655885 :   int s = cmpii(x2, y2);
     872             : 
     873             :   /* common denominator: (x1 op y1) / x2 */
     874    24655886 :   if (!s)
     875             :   {
     876     9026700 :     pari_sp av = avma;
     877     9026700 :     return gerepileupto(av, Qdivii(op(x1, y1), x2));
     878             :   }
     879    15629186 :   z = cgetg(3, t_FRAC);
     880    15629198 :   if (s < 0)
     881             :   {
     882     8804523 :     Q = dvmdii(y2, x2, &r);
     883             :     /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     884     8804427 :     if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
     885     2493024 :     delta = gcdii(x2,r);
     886             :   }
     887             :   else
     888             :   {
     889     6824675 :     Q = dvmdii(x2, y2, &r);
     890             :     /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     891     6824696 :     if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
     892     1306130 :     delta = gcdii(y2,r);
     893             :   }
     894             :   /* delta = gcd(x2,y2) */
     895     3799212 :   if (equali1(delta))
     896             :   { /* numerator is nonzero */
     897     1519056 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     898     1519057 :     gel(z,2) = mulii(x2,y2); return z;
     899             :   }
     900     2280150 :   x2 = diviiexact(x2,delta);
     901     2280151 :   y2 = diviiexact(y2,delta);
     902     2280151 :   n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
     903     2280152 :   q = dvmdii(n, delta, &r);
     904     2280150 :   if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
     905     2084280 :   r = gcdii(delta, r);
     906     2084284 :   if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
     907     2084284 :   return setfrac(z, n, mulii(mulii(x2, y2), delta));
     908             : }
     909             : 
     910             : /* assume x2, y2 are t_POLs in the same variable */
     911             : static GEN
     912     3041742 : add_rfrac(GEN x, GEN y)
     913             : {
     914     3041742 :   pari_sp av = avma;
     915     3041742 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     916     3041742 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     917             : 
     918     3041742 :   delta = RgX_gcd(x2,y2);
     919     3041742 :   if (!degpol(delta))
     920             :   {
     921         665 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     922         665 :     d = RgX_mul(x2, y2);
     923         665 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     924             :   }
     925     3041077 :   x2 = RgX_div(x2,delta);
     926     3041077 :   y2 = RgX_div(y2,delta);
     927     3041077 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     928     3041077 :   if (!signe(n))
     929             :   {
     930      721418 :     n = simplify_shallow(n);
     931      721418 :     if (isexactzero(n))
     932             :     {
     933      721411 :       if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
     934           7 :       return gerepileupto(av, scalarpol(n, varn(x2)));
     935             :     }
     936           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     937             :   }
     938     2319659 :   if (degpol(n) == 0)
     939     1150554 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     940     1169105 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     941     1169105 :   if (isexactzero(r))
     942             :   {
     943             :     GEN z;
     944      228010 :     d = RgX_mul(x2, y2);
     945             :     /* "constant" denominator ? */
     946      228010 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     947      228010 :     return gerepileupto(av, z);
     948             :   }
     949      941095 :   r = RgX_gcd(delta, r);
     950      941095 :   if (degpol(r))
     951             :   {
     952      160726 :     n = RgX_div(n, r);
     953      160726 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     954             :   }
     955             :   else
     956      780369 :     d = RgX_mul(gel(x,2), y2);
     957      941095 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     958             : }
     959             : 
     960             : GEN
     961  4473397305 : gadd(GEN x, GEN y)
     962             : {
     963  4473397305 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     964             :   pari_sp av, tetpil;
     965             :   GEN z, p1;
     966             : 
     967  4473397305 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     968             :   {
     969  1981321952 :     case t_INT: return addii(x,y);
     970  1392709627 :     case t_REAL: return addrr(x,y);
     971     1679202 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     972     1679202 :       z = cgetg(3,t_INTMOD);
     973     1679202 :       if (X==Y || equalii(X,Y))
     974     1679188 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     975          14 :       gel(z,1) = gcdii(X,Y);
     976          14 :       warn_coercion(X,Y,gel(z,1));
     977          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     978          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     979             :     }
     980    19667547 :     case t_FRAC: return addsub_frac(x,y,addii);
     981   270502650 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     982   270511768 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     983   270387980 :       if (isintzero(gel(z,2)))
     984             :       {
     985      362131 :         set_avma((pari_sp)(z+3));
     986      362131 :         return gadd(gel(x,1),gel(y,1));
     987             :       }
     988   270026315 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     989   270024005 :       return z;
     990     2769673 :     case t_PADIC:
     991     2769673 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     992     2769463 :       return addsub_pp(x,y, addii);
     993         672 :     case t_QUAD: z = cgetg(4,t_QUAD);
     994         672 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     995         672 :       gel(z,1) = ZX_copy(gel(x,1));
     996         672 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     997         672 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     998     9857782 :     case t_POLMOD:
     999     9857782 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1000     9855675 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
    1001        2107 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
    1002    14006310 :     case t_FFELT: return FF_add(x,y);
    1003    66605126 :     case t_POL:
    1004    66605126 :       vx = varn(x);
    1005    66605126 :       vy = varn(y);
    1006    66605126 :       if (vx != vy) {
    1007      872067 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
    1008       65742 :         else                     return RgX_Rg_add(y, x);
    1009             :       }
    1010    65733059 :       return RgX_add(x, y);
    1011     3456005 :     case t_SER:
    1012     3456005 :       vx = varn(x);
    1013     3456005 :       vy = varn(y);
    1014     3456005 :       if (vx != vy) {
    1015          21 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1016          21 :         else                     return add_ser_scal(y, x);
    1017             :       }
    1018     3455984 :       return ser_add(x, y);
    1019     4333629 :     case t_RFRAC:
    1020     4333629 :       vx = varn(gel(x,2));
    1021     4333629 :       vy = varn(gel(y,2));
    1022     4333629 :       if (vx != vy) {
    1023     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1024      538397 :         else                     return add_rfrac_scal(y, x);
    1025             :       }
    1026     3041742 :       return add_rfrac(x,y);
    1027      288213 :     case t_VEC:
    1028      288213 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1029      288213 :       return RgV_add(x,y);
    1030     1074681 :     case t_COL:
    1031     1074681 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1032     1074681 :       return RgC_add(x,y);
    1033      671237 :     case t_MAT:
    1034      671237 :       lx = lg(x);
    1035      671237 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1036      671237 :       if (lx == 1) return cgetg(1, t_MAT);
    1037      671237 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1038      671230 :       return RgM_add(x,y);
    1039             : 
    1040           0 :     default: pari_err_TYPE2("+",x,y);
    1041             :   }
    1042             :   /* tx != ty */
    1043   706303478 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1044             : 
    1045   706303478 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1046             :   {
    1047   626251733 :     case t_INT:
    1048             :       switch(ty)
    1049             :       {
    1050   401842987 :         case t_REAL: return addir(x,y);
    1051      278105 :         case t_INTMOD:
    1052      278105 :           z = cgetg(3, t_INTMOD);
    1053      278105 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1054    43598417 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1055    43598403 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1056    43598393 :           gel(z,2) = icopy(gel(y,2)); return z;
    1057   174579305 :         case t_COMPLEX: return addRc(x, y);
    1058     4637450 :         case t_PADIC:
    1059     4637450 :           if (!signe(x)) return gcopy(y);
    1060     3487771 :           return addQp(x,y);
    1061        1113 :         case t_QUAD: return addRq(x, y);
    1062     1356296 :         case t_FFELT: return FF_Z_add(y,x);
    1063             :       }
    1064             : 
    1065             :     case t_REAL:
    1066    44896255 :       switch(ty)
    1067             :       {
    1068    13118398 :         case t_FRAC:
    1069    13118398 :           if (!signe(gel(y,1))) return rcopy(x);
    1070    13118398 :           if (!signe(x))
    1071             :           {
    1072        8473 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1073        8473 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1074             :           }
    1075    13109925 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1076    13107501 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1077    31777794 :         case t_COMPLEX: return addRc(x, y);
    1078          63 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1079             : 
    1080           0 :         default: pari_err_TYPE2("+",x,y);
    1081             :       }
    1082             : 
    1083       17640 :     case t_INTMOD:
    1084             :       switch(ty)
    1085             :       {
    1086       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1087       17507 :           z = cgetg(3, t_INTMOD);
    1088       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1089       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1090             :         }
    1091          14 :         case t_FFELT:
    1092          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1093           0 :             pari_err_OP("+",x,y);
    1094          14 :           return FF_Z_add(y,gel(x,2));
    1095          91 :         case t_COMPLEX: return addRc(x, y);
    1096           0 :         case t_PADIC: { GEN X = gel(x,1);
    1097           0 :           z = cgetg(3, t_INTMOD);
    1098           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1099             :         }
    1100          28 :         case t_QUAD: return addRq(x, y);
    1101             :       }
    1102             : 
    1103             :     case t_FRAC:
    1104             :       switch (ty)
    1105             :       {
    1106    10445443 :         case t_COMPLEX: return addRc(x, y);
    1107       25690 :         case t_PADIC:
    1108       25690 :           if (!signe(gel(x,1))) return gcopy(y);
    1109       25690 :           return addQp(x,y);
    1110         126 :         case t_QUAD: return addRq(x, y);
    1111        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1112             :       }
    1113             : 
    1114             :     case t_FFELT:
    1115           0 :       pari_err_TYPE2("+",x,y);
    1116             : 
    1117          35 :     case t_COMPLEX:
    1118             :       switch(ty)
    1119             :       {
    1120          28 :         case t_PADIC:
    1121          28 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1122           7 :         case t_QUAD:
    1123           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1124           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1125             :       }
    1126             : 
    1127             :     case t_PADIC: /* ty == t_QUAD */
    1128           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1129             :   }
    1130             :   /* tx < ty, !is_const_t(y) */
    1131    26388726 :   switch(ty)
    1132             :   {
    1133        7993 :     case t_MAT:
    1134        7993 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1135        7993 :       if (isrationalzero(x)) return gcopy(y);
    1136        7916 :       return RgM_Rg_add(y, x);
    1137      202306 :     case t_COL:
    1138      202306 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1139      202306 :       return RgC_Rg_add(y, x);
    1140     1838867 :     case t_POLMOD: /* is_const_t(tx) in this case */
    1141     1838867 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1142             :   }
    1143    24339560 :   if (is_scalar_t(tx))  {
    1144    21247240 :     if (tx == t_POLMOD)
    1145             :     {
    1146      106400 :       vx = varn(gel(x,1));
    1147      106400 :       vy = gvar(y);
    1148      106399 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1149             :       else
    1150       81304 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1151       40124 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1152             :     }
    1153    21140840 :     return add_scal(y, x, ty);
    1154             :   }
    1155             :   /* x and y are not scalars, ty != t_MAT */
    1156     3092290 :   vx = gvar(x);
    1157     3092290 :   vy = gvar(y);
    1158     3092290 :   if (vx != vy) { /* x or y is treated as a scalar */
    1159       22682 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1160       32263 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1161       32263 :                                 : add_scal(y, x, ty);
    1162             :   }
    1163             :   /* vx = vy */
    1164     3069608 :   switch(tx)
    1165             :   {
    1166     3069146 :     case t_POL:
    1167             :       switch (ty)
    1168             :       {
    1169        5299 :         case t_SER:
    1170        5299 :           if (lg(x) == 2) return gcopy(y);
    1171        5285 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1172        5285 :           i = lg(y) + valp(y) - i;
    1173        5285 :           if (i < 3) return gcopy(y);
    1174        5278 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1175        5278 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1176        5278 :           return y;
    1177             : 
    1178     3063847 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1179             :       }
    1180           0 :       break;
    1181             : 
    1182         462 :     case t_SER:
    1183         462 :       if (ty == t_RFRAC)
    1184             :       {
    1185             :         long vn, vd;
    1186         462 :         av = avma;
    1187         462 :         vn = gval(gel(y,1), vy);
    1188         462 :         vd = RgX_valrem_inexact(gel(y,2), NULL);
    1189         462 :         if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
    1190             : 
    1191         462 :         l = lg(x) + valp(x) - (vn - vd);
    1192         462 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1193         462 :         return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
    1194             :       }
    1195           0 :       break;
    1196             :   }
    1197           0 :   pari_err_TYPE2("+",x,y);
    1198             :   return NULL; /* LCOV_EXCL_LINE */
    1199             : }
    1200             : 
    1201             : GEN
    1202   283926422 : gaddsg(long x, GEN y)
    1203             : {
    1204   283926422 :   long ty = typ(y);
    1205             :   GEN z;
    1206             : 
    1207   283926422 :   switch(ty)
    1208             :   {
    1209   133240627 :     case t_INT:  return addsi(x,y);
    1210   126998568 :     case t_REAL: return addsr(x,y);
    1211          14 :     case t_INTMOD:
    1212          14 :       z = cgetg(3, t_INTMOD);
    1213          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1214    14710451 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1215    14710451 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1216    14710451 :       gel(z,2) = icopy(gel(y,2)); return z;
    1217     8235120 :     case t_COMPLEX:
    1218     8235120 :       z = cgetg(3, t_COMPLEX);
    1219     8235120 :       gel(z,1) = gaddsg(x, gel(y,1));
    1220     8235120 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1221             : 
    1222      741642 :     default: return gadd(stoi(x), y);
    1223             :   }
    1224             : }
    1225             : 
    1226             : GEN
    1227     3144202 : gsubsg(long x, GEN y)
    1228             : {
    1229             :   GEN z, a, b;
    1230             :   pari_sp av;
    1231             : 
    1232     3144202 :   switch(typ(y))
    1233             :   {
    1234      274342 :     case t_INT:  return subsi(x,y);
    1235     1212608 :     case t_REAL: return subsr(x,y);
    1236          56 :     case t_INTMOD:
    1237          56 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1238          56 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1239      729137 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1240      729137 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1241      729134 :       gel(z,2) = icopy(gel(y,2)); return z;
    1242      891337 :     case t_COMPLEX:
    1243      891337 :       z = cgetg(3, t_COMPLEX);
    1244      891337 :       gel(z,1) = gsubsg(x, gel(y,1));
    1245      891337 :       gel(z,2) = gneg(gel(y,2)); return z;
    1246             :   }
    1247       36722 :   av = avma;
    1248       36722 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1249             : }
    1250             : 
    1251             : /********************************************************************/
    1252             : /**                                                                **/
    1253             : /**                          SUBTRACTION                           **/
    1254             : /**                                                                **/
    1255             : /********************************************************************/
    1256             : 
    1257             : GEN
    1258  3427362920 : gsub(GEN x, GEN y)
    1259             : {
    1260  3427362920 :   long tx = typ(x), ty = typ(y);
    1261             :   pari_sp av;
    1262             :   GEN z;
    1263  3427362920 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1264             :   {
    1265  2743228627 :     case t_INT: return subii(x,y);
    1266   495630649 :     case t_REAL: return subrr(x,y);
    1267     1153268 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1268     1153268 :       z = cgetg(3,t_INTMOD);
    1269     1153268 :       if (X==Y || equalii(X,Y))
    1270     1153254 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1271          14 :       gel(z,1) = gcdii(X,Y);
    1272          14 :       warn_coercion(X,Y,gel(z,1));
    1273          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1274          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1275             :     }
    1276     4988368 :     case t_FRAC: return addsub_frac(x,y, subii);
    1277   107628951 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1278   107669373 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1279   107558093 :       if (isintzero(gel(z,2)))
    1280             :       {
    1281       21237 :         set_avma((pari_sp)(z+3));
    1282       21237 :         return gsub(gel(x,1),gel(y,1));
    1283             :       }
    1284   107540652 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1285   107526884 :       return z;
    1286     7762227 :     case t_PADIC:
    1287     7762227 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1288     7762227 :       return addsub_pp(x,y, subii);
    1289         168 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1290         168 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1291         168 :       gel(z,1) = ZX_copy(gel(x,1));
    1292         168 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1293         168 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1294      793265 :     case t_POLMOD:
    1295      793265 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1296      793230 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1297          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1298      203128 :     case t_FFELT: return FF_sub(x,y);
    1299    10219580 :     case t_POL: {
    1300    10219580 :       long vx = varn(x);
    1301    10219580 :       long vy = varn(y);
    1302    10219580 :       if (vx != vy) {
    1303       30219 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1304        5614 :         else                     return Rg_RgX_sub(x, y);
    1305             :       }
    1306    10189361 :       return RgX_sub(x, y);
    1307             :     }
    1308      297384 :     case t_VEC:
    1309      297384 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1310      297384 :       return RgV_sub(x,y);
    1311     6534220 :     case t_COL:
    1312     6534220 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1313     6534220 :       return RgC_sub(x,y);
    1314      234822 :     case t_MAT: {
    1315      234822 :       long lx = lg(x);
    1316      234822 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1317      234821 :       if (lx == 1) return cgetg(1, t_MAT);
    1318      155735 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1319      155735 :       return RgM_sub(x,y);
    1320             :     }
    1321     2452440 :     case t_RFRAC: case t_SER: break;
    1322             : 
    1323           0 :     default: pari_err_TYPE2("+",x,y);
    1324             :   }
    1325    47891231 :   av = avma;
    1326    50343671 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1327             : }
    1328             : 
    1329             : /********************************************************************/
    1330             : /**                                                                **/
    1331             : /**                        MULTIPLICATION                          **/
    1332             : /**                                                                **/
    1333             : /********************************************************************/
    1334             : static GEN
    1335      302933 : mul_ser_scal(GEN y, GEN x)
    1336             : {
    1337             :   long l, i;
    1338             :   GEN z;
    1339      302933 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1340      299685 :   if (ser_isexactzero(y))
    1341             :   {
    1342         476 :     z = scalarser(lg(y) == 2? Rg_get_0(x): gmul(gel(y,2), x), varn(y), 1);
    1343         476 :     setvalp(z, valp(y)); return z;
    1344             :   }
    1345      299209 :   z = cgetg_copy(y, &l); z[1] = y[1];
    1346     4167984 :   for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
    1347      299209 :   return normalizeser(z);
    1348             : }
    1349             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1350             :  * [n/d a valid RFRAC]  */
    1351             : static GEN
    1352    10451122 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1353             : {
    1354    10451122 :   pari_sp av = avma;
    1355             :   GEN z;
    1356             : 
    1357    10451122 :   switch(typ(x))
    1358             :   {
    1359          21 :     case t_PADIC:
    1360          21 :       n = gmul(n, x);
    1361          21 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1362          21 :       return gerepileupto(av, gdiv(n,d));
    1363             : 
    1364         490 :     case t_INTMOD: case t_POLMOD:
    1365         490 :       n = gmul(n, x);
    1366         490 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1367         490 :       return gerepileupto(av, gdiv(n,d));
    1368             :   }
    1369    10450611 :   z = gred_rfrac2(x, d);
    1370    10450611 :   n = simplify_shallow(n);
    1371    10450611 :   if (typ(z) == t_RFRAC)
    1372             :   {
    1373     7933879 :     n = gmul(gel(z,1), n);
    1374     7933879 :     d = gel(z,2);
    1375     7933879 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1376           0 :       z = RgX_Rg_div(n, d);
    1377             :     else
    1378     7933879 :       z = gred_rfrac_simple(n, d);
    1379             :   }
    1380             :   else
    1381     2516732 :     z = gmul(z, n);
    1382    10450611 :   return gerepileupto(av, z);
    1383             : }
    1384             : static GEN
    1385   111560892 : mul_scal(GEN y, GEN x, long ty)
    1386             : {
    1387   111560892 :   switch(ty)
    1388             :   {
    1389   102679427 :     case t_POL:
    1390   102679427 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1391    99527959 :       return RgX_Rg_mul(y, x);
    1392      294848 :     case t_SER: return mul_ser_scal(y, x);
    1393     8586627 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1394          14 :     case t_QFB:
    1395          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1396             :   }
    1397          12 :   pari_err_TYPE2("*",x,y);
    1398             :   return NULL; /* LCOV_EXCL_LINE */
    1399             : }
    1400             : 
    1401             : static GEN
    1402      160865 : mul_gen_rfrac(GEN X, GEN Y)
    1403             : {
    1404      160865 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1405      160865 :   long vx = gvar(X), vy = varn(y2);
    1406      166626 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1407        5761 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1408             : }
    1409             : /* (x1/x2) * (y1/y2) */
    1410             : static GEN
    1411     7907093 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1412             : {
    1413             :   GEN z, X, Y;
    1414     7907093 :   pari_sp av = avma;
    1415             : 
    1416     7907093 :   X = gred_rfrac2(x1, y2);
    1417     7907093 :   Y = gred_rfrac2(y1, x2);
    1418     7907093 :   if (typ(X) == t_RFRAC)
    1419             :   {
    1420     6628032 :     if (typ(Y) == t_RFRAC) {
    1421     6562456 :       x1 = gel(X,1);
    1422     6562456 :       x2 = gel(X,2);
    1423     6562456 :       y1 = gel(Y,1);
    1424     6562456 :       y2 = gel(Y,2);
    1425     6562456 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1426             :     } else
    1427       65576 :       z = mul_gen_rfrac(Y, X);
    1428             :   }
    1429     1279061 :   else if (typ(Y) == t_RFRAC)
    1430       95289 :     z = mul_gen_rfrac(X, Y);
    1431             :   else
    1432     1183772 :     z = gmul(X, Y);
    1433     7907093 :   return gerepileupto(av, z);
    1434             : }
    1435             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1436             : static GEN
    1437      269687 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1438             : {
    1439      269687 :   pari_sp av = avma;
    1440      269687 :   GEN X = gred_rfrac2(x1, y2);
    1441      269687 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1442             :   {
    1443      263121 :     x2 = RgX_mul(gel(X,2), x2);
    1444      263121 :     x1 = gel(X,1);
    1445             :   }
    1446             :   else
    1447        6566 :     x1 = X;
    1448      269687 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1449             : }
    1450             : 
    1451             : /* Mod(y, Y) * x,  assuming x scalar */
    1452             : static GEN
    1453     2571802 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1454     2571802 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1455             : 
    1456             : /* cf mulqq */
    1457             : static GEN
    1458     5854148 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1459             : {
    1460     5854148 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1461     5854148 :   pari_sp tetpil, av = avma;
    1462     5854148 :   T[1] = x[1];
    1463     5854148 :   p2 = gmul(gel(x,2), gel(y,2));
    1464     5854148 :   p3 = gmul(gel(x,3), gel(y,3));
    1465     5854148 :   p1 = gmul(gneg_i(c),p3);
    1466             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1467     5854148 :   if (typ(b) == t_INT)
    1468             :   {
    1469     5854127 :     if (signe(b))
    1470             :     {
    1471     4284657 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1472     4284657 :       if (is_pm1(b))
    1473             :       {
    1474     4283999 :         if (signe(b) > 0) p3 = gneg(p3);
    1475             :       }
    1476             :       else
    1477         658 :         p3 = gmul(negi(b), p3);
    1478             :     }
    1479             :     else
    1480             :     {
    1481     1569470 :       p3 = gmul(gel(x,2),gel(y,3));
    1482     1569470 :       p4 = gmul(gel(x,3),gel(y,2));
    1483             :     }
    1484             :   }
    1485             :   else
    1486             :   {
    1487          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1488          21 :     p3 = gmul(gneg_i(b), p3);
    1489             :   }
    1490     5854148 :   tetpil = avma;
    1491     5854148 :   gel(T,2) = gadd(p2, p1);
    1492     5854148 :   gel(T,3) = gadd(p4, p3);
    1493     5854148 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1494     5854148 :   return normalizepol_lg(T,4);
    1495             : }
    1496             : /* Mod(x,T) * Mod(y,T) */
    1497             : static GEN
    1498     7782603 : mul_polmod_same(GEN T, GEN x, GEN y)
    1499             : {
    1500     7782603 :   GEN z = cgetg(3,t_POLMOD), a;
    1501     7782603 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1502     7782603 :   gel(z,1) = RgX_copy(T);
    1503             :   /* x * y mod T optimised */
    1504     7782603 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1505     7163285 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1506     1489363 :     a = gmul(x, y);
    1507             :   else
    1508             :   {
    1509     6293240 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1510     5849003 :       a = quad_polmod_mul(T, x, y);
    1511             :     else
    1512      444237 :       a = RgXQ_mul(x, y, gel(z,1));
    1513             :   }
    1514     7782603 :   gel(z,2) = a; return z;
    1515             : }
    1516             : static GEN
    1517       46880 : sqr_polmod(GEN T, GEN x)
    1518             : {
    1519       46880 :   GEN a, z = cgetg(3,t_POLMOD);
    1520       46880 :   gel(z,1) = RgX_copy(T);
    1521       46880 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1522       14226 :     a = gsqr(x);
    1523             :   else
    1524             :   {
    1525       32654 :     pari_sp av = avma;
    1526       32654 :     a = RgXQ_sqr(x, gel(z,1));
    1527       32654 :     a = gerepileupto(av, a);
    1528             :   }
    1529       46880 :   gel(z,2) = a; return z;
    1530             : }
    1531             : /* Mod(x,X) * Mod(y,Y) */
    1532             : static GEN
    1533     2678207 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1534             : {
    1535     2678207 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1536     2678207 :   long vx = varn(X), vy = varn(Y);
    1537     2678207 :   GEN z = cgetg(3,t_POLMOD);
    1538             : 
    1539     2678207 :   if (vx==vy) {
    1540             :     pari_sp av;
    1541          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1542          14 :     warn_coercion(X,Y,gel(z,1));
    1543          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1544          14 :     return z;
    1545             :   }
    1546     2678193 :   if (varncmp(vx, vy) < 0)
    1547      414316 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1548             :   else
    1549     2263877 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1550     2678193 :   gel(z,2) = gmul(x, y); return z;
    1551             : }
    1552             : 
    1553             : #if 0 /* used by 3M only */
    1554             : /* set z = x+y and return 1 if x,y have the same sign
    1555             :  * set z = x-y and return 0 otherwise */
    1556             : static int
    1557             : did_add(GEN x, GEN y, GEN *z)
    1558             : {
    1559             :   long tx = typ(x), ty = typ(y);
    1560             :   if (tx == ty) switch(tx)
    1561             :   {
    1562             :     case t_INT: *z = addii(x,y); return 1;
    1563             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1564             :     case t_REAL:
    1565             :       if (signe(x) == -signe(y))
    1566             :       { *z = subrr(x,y); return 0; }
    1567             :       else
    1568             :       { *z = addrr(x,y); return 1; }
    1569             :   }
    1570             :   if (tx == t_REAL) switch(ty)
    1571             :   {
    1572             :     case t_INT:
    1573             :       if (signe(x) == -signe(y))
    1574             :       { *z = subri(x,y); return 0; }
    1575             :       else
    1576             :       { *z = addri(x,y); return 1; }
    1577             :     case t_FRAC:
    1578             :       if (signe(x) == -signe(gel(y,1)))
    1579             :       { *z = gsub(x,y); return 0; }
    1580             :       else
    1581             :       { *z = gadd(x,y); return 1; }
    1582             :   }
    1583             :   else if (ty == t_REAL) switch(tx)
    1584             :   {
    1585             :     case t_INT:
    1586             :       if (signe(x) == -signe(y))
    1587             :       { *z = subir(x,y); return 0; }
    1588             :       else
    1589             :       { *z = addir(x,y); return 1; }
    1590             :     case t_FRAC:
    1591             :       if (signe(gel(x,1)) == -signe(y))
    1592             :       { *z = gsub(x,y); return 0; }
    1593             :       else
    1594             :       { *z = gadd(x,y); return 1; }
    1595             :   }
    1596             :   *z = gadd(x,y); return 1;
    1597             : }
    1598             : #endif
    1599             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1600             : static GEN
    1601    10995999 : mulcIR(GEN x, GEN y)
    1602             : {
    1603    10995999 :   GEN z = cgetg(3,t_COMPLEX);
    1604    10996196 :   pari_sp av = avma;
    1605    10996196 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1606    10996414 :   gel(z,2) = gmul(y, gel(x,1));
    1607    10996029 :   return z;
    1608             : 
    1609             : }
    1610             : /* x,y COMPLEX */
    1611             : static GEN
    1612   246951223 : mulcc(GEN x, GEN y)
    1613             : {
    1614   246951223 :   GEN xr = gel(x,1), xi = gel(x,2);
    1615   246951223 :   GEN yr = gel(y,1), yi = gel(y,2);
    1616             :   GEN p1, p2, p3, p4, z;
    1617             :   pari_sp tetpil, av;
    1618             : 
    1619   246951223 :   if (isintzero(xr))
    1620             :   {
    1621    15315067 :     if (isintzero(yr)) {
    1622     7160003 :       av = avma;
    1623     7160003 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1624             :     }
    1625     8154855 :     return mulcIR(y, xi);
    1626             :   }
    1627   231636687 :   if (isintzero(yr)) return mulcIR(x, yi);
    1628             : 
    1629   228805263 :   z = cgetg(3,t_COMPLEX); av = avma;
    1630             : #if 0
    1631             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1632             :    * e.g. xr + xi if exponents differ */
    1633             :   if (did_add(xr, xi, &p3))
    1634             :   {
    1635             :     if (did_add(yr, yi, &p4)) {
    1636             :     /* R = xr*yr - xi*yi
    1637             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1638             :       p1 = gmul(xr,yr);
    1639             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1640             :       p3 = gmul(p3, p4);
    1641             :       p4 = gsub(p2, p1);
    1642             :     } else {
    1643             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1644             :      * I = xr*yi + xi*yr */
    1645             :       p1 = gmul(p3,p4);
    1646             :       p3 = gmul(xr,yi);
    1647             :       p4 = gmul(xi,yr);
    1648             :       p2 = gsub(p3, p4);
    1649             :     }
    1650             :   } else {
    1651             :     if (did_add(yr, yi, &p4)) {
    1652             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1653             :       * I = xr*yi +xi*yr */
    1654             :       p1 = gmul(p3,p4);
    1655             :       p3 = gmul(xr,yi);
    1656             :       p4 = gmul(xi,yr);
    1657             :       p2 = gsub(p4, p3);
    1658             :     } else {
    1659             :     /* R = xr*yr - xi*yi
    1660             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1661             :       p3 = gneg( gmul(p3, p4) );
    1662             :       p1 = gmul(xr,yr);
    1663             :       p2 = gmul(xi,yi);
    1664             :       p4 = gadd(p1, p2);
    1665             : 
    1666             :       p2 = gneg(p2);
    1667             :     }
    1668             :   }
    1669             :   tetpil = avma;
    1670             :   gel(z,1) = gadd(p1,p2);
    1671             :   gel(z,2) = gadd(p3,p4);
    1672             : #else
    1673   228826886 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1674             :   { /* 3M formula */
    1675      557442 :     p3 = addii(xr,xi);
    1676      557442 :     p4 = addii(yr,yi);
    1677      557442 :     p1 = mulii(xr,yr);
    1678      557442 :     p2 = mulii(xi,yi);
    1679      557442 :     p3 = mulii(p3,p4);
    1680      557442 :     p4 = addii(p2,p1);
    1681      557442 :     tetpil = avma;
    1682      557442 :     gel(z,1) = subii(p1,p2);
    1683      557442 :     gel(z,2) = subii(p3,p4);
    1684      557442 :     if (!signe(gel(z,2)))
    1685      113057 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1686             :   }
    1687             :   else
    1688             :   { /* naive 4M formula: avoid all loss of accuracy */
    1689   228269444 :     p1 = gmul(xr,yr);
    1690   228196583 :     p2 = gmul(xi,yi);
    1691   228176803 :     p3 = gmul(xr,yi);
    1692   228180222 :     p4 = gmul(xi,yr);
    1693   228188600 :     tetpil = avma;
    1694   228188600 :     gel(z,1) = gsub(p1,p2);
    1695   228098981 :     gel(z,2) = gadd(p3,p4);
    1696   228104692 :     if (isintzero(gel(z,2)))
    1697             :     {
    1698       31387 :       cgiv(gel(z,2));
    1699       31387 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1700             :     }
    1701             :   }
    1702             : #endif
    1703             : 
    1704   228519035 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1705             : }
    1706             : /* x,y PADIC */
    1707             : static GEN
    1708    16846772 : mulpp(GEN x, GEN y) {
    1709    16846772 :   long l = valp(x) + valp(y);
    1710             :   pari_sp av;
    1711             :   GEN z, t;
    1712    16846772 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1713    16846730 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1714    16630093 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1715             : 
    1716    16293128 :   t = (precp(x) > precp(y))? y: x;
    1717    16293128 :   z = cgetp(t); setvalp(z,l); av = avma;
    1718    16293231 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1719    16292991 :   set_avma(av); return z;
    1720             : }
    1721             : /* x,y QUAD */
    1722             : static GEN
    1723        1106 : mulqq(GEN x, GEN y) {
    1724        1106 :   GEN z = cgetg(4,t_QUAD);
    1725        1106 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1726             :   pari_sp av, tetpil;
    1727        1106 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1728             : 
    1729        1106 :   gel(z,1) = ZX_copy(P); av = avma;
    1730        1106 :   p2 = gmul(gel(x,2),gel(y,2));
    1731        1106 :   p3 = gmul(gel(x,3),gel(y,3));
    1732        1106 :   p1 = gmul(gneg_i(c),p3);
    1733             : 
    1734        1106 :   if (signe(b))
    1735         987 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1736             :   else
    1737             :   {
    1738         119 :     p3 = gmul(gel(x,2),gel(y,3));
    1739         119 :     p4 = gmul(gel(x,3),gel(y,2));
    1740             :   }
    1741        1106 :   tetpil = avma;
    1742        1106 :   gel(z,2) = gadd(p2,p1);
    1743        1106 :   gel(z,3) = gadd(p4,p3);
    1744        1106 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1745             : }
    1746             : 
    1747             : GEN
    1748    10320816 : mulcxI(GEN x)
    1749             : {
    1750             :   GEN z;
    1751    10320816 :   switch(typ(x))
    1752             :   {
    1753     1561452 :     case t_INT: case t_REAL: case t_FRAC:
    1754     1561452 :       return mkcomplex(gen_0, x);
    1755     8759908 :     case t_COMPLEX:
    1756     8759908 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1757     8759278 :       z = cgetg(3,t_COMPLEX);
    1758     8759846 :       gel(z,1) = gneg(gel(x,2));
    1759     8760331 :       gel(z,2) = gel(x,1); return z;
    1760          48 :     default:
    1761          48 :       return gmul(gen_I(), x);
    1762             :   }
    1763             : }
    1764             : GEN
    1765      196014 : mulcxmI(GEN x)
    1766             : {
    1767             :   GEN z;
    1768      196014 :   switch(typ(x))
    1769             :   {
    1770       84922 :     case t_INT: case t_REAL: case t_FRAC:
    1771       84922 :       return mkcomplex(gen_0, gneg(x));
    1772      111029 :     case t_COMPLEX:
    1773      111029 :       if (isintzero(gel(x,1))) return gel(x,2);
    1774      108418 :       z = cgetg(3,t_COMPLEX);
    1775      108418 :       gel(z,1) = gel(x,2);
    1776      108418 :       gel(z,2) = gneg(gel(x,1)); return z;
    1777          63 :     default:
    1778          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1779             :   }
    1780             : }
    1781             : /* x * I^k */
    1782             : GEN
    1783      259031 : mulcxpowIs(GEN x, long k)
    1784             : {
    1785      259031 :   switch (k & 3)
    1786             :   {
    1787       49649 :     case 1: return mulcxI(x);
    1788       46337 :     case 2: return gneg(x);
    1789       49968 :     case 3: return mulcxmI(x);
    1790             :   }
    1791      113077 :   return x;
    1792             : }
    1793             : 
    1794             : static GEN
    1795     5394735 : init_ser(long l, long v, long e)
    1796             : {
    1797     5394735 :   GEN z = cgetg(l, t_SER);
    1798     5394735 :   z[1] = evalvalp(e) | evalvarn(v) | evalsigne(1); return z;
    1799             : }
    1800             : 
    1801             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1802             : static GEN
    1803     5379077 : fill_ser(GEN z, GEN y)
    1804             : {
    1805     5379077 :   long i, lx = lg(z), ly = lg(y);
    1806     5379077 :   if (ly >= lx) {
    1807    24020714 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1808             :   } else {
    1809      328331 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1810      237900 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1811             :   }
    1812     5379077 :   return normalizeser(z);
    1813             : }
    1814             : 
    1815             : GEN
    1816  7748018530 : gmul(GEN x, GEN y)
    1817             : {
    1818             :   long tx, ty, lx, ly, vx, vy, i, l;
    1819             :   pari_sp av, tetpil;
    1820             :   GEN z, p1;
    1821             : 
    1822  7748018530 :   if (x == y) return gsqr(x);
    1823  6875357121 :   tx = typ(x); ty = typ(y);
    1824  6875357121 :   if (tx == ty) switch(tx)
    1825             :   {
    1826  3418375397 :     case t_INT: return mulii(x,y);
    1827  1632440238 :     case t_REAL: return mulrr(x,y);
    1828     2048342 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1829     2048342 :       z = cgetg(3,t_INTMOD);
    1830     2048342 :       if (X==Y || equalii(X,Y))
    1831     2048306 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1832          35 :       gel(z,1) = gcdii(X,Y);
    1833          35 :       warn_coercion(X,Y,gel(z,1));
    1834          35 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1835          35 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1836             :     }
    1837    20622238 :     case t_FRAC:
    1838             :     {
    1839    20622238 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1840    20622238 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1841    20622238 :       z=cgetg(3,t_FRAC);
    1842    20622238 :       p1 = gcdii(x1, y2);
    1843    20622229 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1844    20622228 :       p1 = gcdii(x2, y1);
    1845    20622230 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1846    20622229 :       tetpil = avma;
    1847    20622229 :       gel(z,2) = mulii(x2,y2);
    1848    20622225 :       gel(z,1) = mulii(x1,y1);
    1849    20622225 :       fix_frac_if_int_GC(z,tetpil); return z;
    1850             :     }
    1851   239363270 :     case t_COMPLEX: return mulcc(x, y);
    1852    10344835 :     case t_PADIC: return mulpp(x, y);
    1853         882 :     case t_QUAD: return mulqq(x, y);
    1854    14523401 :     case t_FFELT: return FF_mul(x, y);
    1855    10213227 :     case t_POLMOD:
    1856    10213227 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1857     7535020 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1858     2678207 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1859    43948738 :     case t_POL:
    1860    43948738 :       vx = varn(x);
    1861    43948738 :       vy = varn(y);
    1862    43948738 :       if (vx != vy) {
    1863     4867710 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1864     2082600 :         else                     return RgX_Rg_mul(y, x);
    1865             :       }
    1866    39081028 :       return RgX_mul(x, y);
    1867             : 
    1868     4224360 :     case t_SER: {
    1869     4224360 :       vx = varn(x);
    1870     4224360 :       vy = varn(y);
    1871     4224360 :       if (vx != vy) {
    1872        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1873        3675 :         return mul_ser_scal(y, x);
    1874             :       }
    1875     4220685 :       lx = minss(lg(x), lg(y));
    1876     4220685 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1877     4220671 :       av = avma; z = init_ser(lx, vx, valp(x)+valp(y));
    1878     4220671 :       x = ser2pol_i(x, lx);
    1879     4220671 :       y = ser2pol_i(y, lx);
    1880     4220671 :       y = RgXn_mul(x, y, lx-2);
    1881     4220671 :       return gerepilecopy(av, fill_ser(z,y));
    1882             :     }
    1883       55815 :     case t_VEC: /* handle extended t_QFB */
    1884       55815 :     case t_QFB: return qfbcomp(x,y);
    1885     6720458 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1886     2749259 :     case t_MAT: return RgM_mul(x, y);
    1887             : 
    1888        1631 :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1889        1631 :       z = cgetg_copy(x, &l);
    1890        1631 :       if (l != lg(y)) break;
    1891       17325 :       for (i=1; i<l; i++)
    1892             :       {
    1893       15694 :         long yi = y[i];
    1894       15694 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1895       15694 :         z[i] = x[yi];
    1896             :       }
    1897        1631 :       return z;
    1898             : 
    1899           0 :     default:
    1900           0 :       pari_err_TYPE2("*",x,y);
    1901             :   }
    1902             :   /* tx != ty */
    1903  1471528078 :   if (is_const_t(ty) && is_const_t(tx))  {
    1904  1349497520 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1905  1349497520 :     switch(tx) {
    1906  1243608215 :     case t_INT:
    1907             :       switch(ty)
    1908             :       {
    1909   873849758 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1910     1169776 :         case t_INTMOD:
    1911     1169776 :           z = cgetg(3, t_INTMOD);
    1912     1169775 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1913    62194224 :         case t_FRAC:
    1914    62194224 :           if (!signe(x)) return gen_0;
    1915    44976705 :           z=cgetg(3,t_FRAC);
    1916    44976891 :           p1 = gcdii(x,gel(y,2));
    1917    44975585 :           if (equali1(p1))
    1918             :           {
    1919    27376779 :             set_avma((pari_sp)z);
    1920    27376809 :             gel(z,2) = icopy(gel(y,2));
    1921    27376951 :             gel(z,1) = mulii(gel(y,1), x);
    1922             :           }
    1923             :           else
    1924             :           {
    1925    17599156 :             x = diviiexact(x,p1); tetpil = avma;
    1926    17598424 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1927    17598294 :             gel(z,1) = mulii(gel(y,1), x);
    1928    17598714 :             fix_frac_if_int_GC(z,tetpil);
    1929             :           }
    1930    44976416 :           return z;
    1931   299284639 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1932     5671998 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1933        1701 :         case t_QUAD: return mulRq(x,y);
    1934     1622167 :         case t_FFELT: return FF_Z_mul(y,x);
    1935             :       }
    1936             : 
    1937             :     case t_REAL:
    1938    97616005 :       switch(ty)
    1939             :       {
    1940    29262858 :         case t_FRAC: return mulrfrac(x, y);
    1941    68353105 :         case t_COMPLEX: return mulRc(x, y);
    1942          21 :         case t_QUAD: return mulqf(y, x, realprec(x));
    1943          21 :         default: pari_err_TYPE2("*",x,y);
    1944             :       }
    1945             : 
    1946        8064 :     case t_INTMOD:
    1947             :       switch(ty)
    1948             :       {
    1949        7350 :         case t_FRAC: { GEN X = gel(x,1);
    1950        7350 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1951        7350 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1952             :         }
    1953          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1954         252 :         case t_PADIC: { GEN X = gel(x,1);
    1955         252 :           z = cgetg(3, t_INTMOD);
    1956         252 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1957             :         }
    1958          63 :         case t_QUAD: return mulRq(x, y);
    1959         350 :         case t_FFELT:
    1960         350 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1961           0 :             pari_err_OP("*",x,y);
    1962         350 :           return FF_Z_mul(y,gel(x,2));
    1963             :       }
    1964             : 
    1965             :     case t_FRAC:
    1966             :       switch(ty)
    1967             :       {
    1968     5103845 :         case t_COMPLEX: return mulRc(x, y);
    1969     3335979 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1970         343 :         case t_QUAD:  return mulRq(x, y);
    1971        1995 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1972             :       }
    1973             : 
    1974             :     case t_FFELT:
    1975          35 :       pari_err_TYPE2("*",x,y);
    1976             : 
    1977          21 :     case t_COMPLEX:
    1978             :       switch(ty)
    1979             :       {
    1980          14 :         case t_PADIC:
    1981          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1982           7 :         case t_QUAD:
    1983           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1984           7 :           return mulqf(y, x, lx);
    1985             :       }
    1986             : 
    1987             :     case t_PADIC: /* ty == t_QUAD */
    1988          28 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1989             :     }
    1990             :   }
    1991             : 
    1992   125740635 :   if (is_matvec_t(ty))
    1993             :   {
    1994     9001154 :     if (!is_matvec_t(tx))
    1995             :     {
    1996     8658933 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1997     8658931 :       z = cgetg_copy(y, &ly);
    1998  1407306797 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1999     8658471 :       return z;
    2000             :     }
    2001      342220 :     switch(tx)
    2002             :     {
    2003      120754 :       case t_VEC:
    2004             :         switch(ty) {
    2005       15456 :           case t_COL: return RgV_RgC_mul(x,y);
    2006      105298 :           case t_MAT: return RgV_RgM_mul(x,y);
    2007             :         }
    2008           0 :         break;
    2009        1687 :       case t_COL:
    2010             :         switch(ty) {
    2011        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    2012           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    2013             :         }
    2014           0 :         break;
    2015      219796 :       case t_MAT:
    2016             :         switch(ty) {
    2017           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    2018      219796 :           case t_COL: return RgM_RgC_mul(x,y);
    2019             :         }
    2020             :     }
    2021             :   }
    2022   117205253 :   if (is_matvec_t(tx))
    2023             :   {
    2024     1109328 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2025     1109542 :     z = cgetg_copy(x, &lx);
    2026     5290268 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2027     1109249 :     return z;
    2028             :   }
    2029   116095979 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2030             :   /* tx < ty, !ismatvec(x and y) */
    2031             : 
    2032   116095979 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2033     2524433 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2034   113571546 :   if (is_scalar_t(tx)) {
    2035   110281440 :     if (tx == t_POLMOD) {
    2036     3098711 :       vx = varn(gel(x,1));
    2037     3098711 :       vy = gvar(y);
    2038     3098711 :       if (vx != vy) {
    2039     2851688 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2040       47369 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2041             :       }
    2042             :       /* error if ty == t_SER */
    2043      247023 :       av = avma; y = gmod(y, gel(x,1));
    2044      247016 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2045             :     }
    2046   107182729 :     return mul_scal(y, x, ty);
    2047             :   }
    2048             : 
    2049             :   /* x and y are not scalars, nor matvec */
    2050     3289881 :   vx = gvar(x);
    2051     3289914 :   vy = gvar(y);
    2052     3289914 :   if (vx != vy) /* x or y is treated as a scalar */
    2053     2780047 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2054     2780047 :                                 : mul_scal(y, x, ty);
    2055             :   /* vx = vy */
    2056     1871138 :   switch(tx)
    2057             :   {
    2058     1871110 :     case t_POL:
    2059             :       switch (ty)
    2060             :       {
    2061        6685 :         case t_SER:
    2062             :         {
    2063             :           long v;
    2064        6685 :           av = avma; v = RgX_valrem(x, &x);
    2065        6685 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2066        6671 :           v += valp(y); ly = lg(y);
    2067        6671 :           if (ly == 2) { set_avma(av); return zeroser(vx, v); }
    2068        6440 :           if (degpol(x))
    2069             :           {
    2070        2030 :             z = init_ser(ly, vx, v);
    2071        2030 :             x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
    2072        2030 :             return gerepilecopy(av, fill_ser(z, x));
    2073             :           }
    2074             :           /* take advantage of x = c*t^v */
    2075        4410 :           set_avma(av); y = mul_ser_scal(y, gel(x,2));
    2076        4410 :           setvalp(y, v); return y;
    2077             :         }
    2078             : 
    2079     1864425 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2080             :       }
    2081           0 :       break;
    2082             : 
    2083          28 :     case t_SER:
    2084             :       switch (ty)
    2085             :       {
    2086          28 :         case t_RFRAC:
    2087          28 :           av = avma;
    2088          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2089             :       }
    2090           0 :       break;
    2091             :   }
    2092           0 :   pari_err_TYPE2("*",x,y);
    2093             :   return NULL; /* LCOV_EXCL_LINE */
    2094             : }
    2095             : 
    2096             : /* return a nonnormalized result */
    2097             : GEN
    2098      112829 : sqr_ser_part(GEN x, long l1, long l2)
    2099             : {
    2100             :   long i, j, l;
    2101             :   pari_sp av;
    2102             :   GEN Z, z, p1, p2;
    2103             :   long mi;
    2104      112829 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2105      112815 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2106      112815 :   Z = cgetg(l2-l1+3, t_SER);
    2107      112815 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2108      112815 :   z = Z + 2-l1;
    2109      112815 :   x += 2; mi = 0;
    2110      441896 :   for (i=0; i<l1; i++)
    2111             :   {
    2112      329081 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2113             :   }
    2114             : 
    2115      712344 :   for (i=l1; i<=l2; i++)
    2116             :   {
    2117      599529 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2118      599529 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2119     1248434 :     for (j=i-mi; j<=minss(l,mi); j++)
    2120      648905 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2121      599529 :     p1 = gshift(p1,1);
    2122      599529 :     if ((i&1) == 0 && p2[i>>1])
    2123       88356 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2124      599529 :     gel(z,i) = gerepileupto(av,p1);
    2125             :   }
    2126      112815 :   return Z;
    2127             : }
    2128             : 
    2129             : GEN
    2130  1055889643 : gsqr(GEN x)
    2131             : {
    2132             :   long i, lx;
    2133             :   pari_sp av, tetpil;
    2134             :   GEN z, p1, p2, p3, p4;
    2135             : 
    2136  1055889643 :   switch(typ(x))
    2137             :   {
    2138   882119462 :     case t_INT: return sqri(x);
    2139   160220509 :     case t_REAL: return sqrr(x);
    2140      142692 :     case t_INTMOD: { GEN X = gel(x,1);
    2141      142692 :       z = cgetg(3,t_INTMOD);
    2142      142693 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2143      142694 :       gel(z,1) = icopy(X); return z;
    2144             :     }
    2145     2976002 :     case t_FRAC: return sqrfrac(x);
    2146             : 
    2147     6592751 :     case t_COMPLEX:
    2148     6592751 :       if (isintzero(gel(x,1))) {
    2149      209263 :         av = avma;
    2150      209263 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2151             :       }
    2152     6383493 :       z = cgetg(3,t_COMPLEX); av = avma;
    2153     6383482 :       p1 = gadd(gel(x,1),gel(x,2));
    2154     6383321 :       p2 = gsub(gel(x,1), gel(x,2));
    2155     6383245 :       p3 = gmul(gel(x,1),gel(x,2));
    2156     6383347 :       tetpil = avma;
    2157     6383347 :       gel(z,1) = gmul(p1,p2);
    2158     6383370 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2159             : 
    2160        2583 :     case t_PADIC:
    2161        2583 :       z = cgetg(5,t_PADIC);
    2162        2583 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2163        2583 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2164        2583 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2165        2583 :       gel(z,2) = icopy(gel(x,2));
    2166        2583 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2167        2583 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2168        2583 :       return z;
    2169             : 
    2170          70 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2171          70 :       p1 = gel(x,1);
    2172          70 :       gel(z,1) = ZX_copy(p1); av = avma;
    2173          70 :       p2 = gsqr(gel(x,2));
    2174          70 :       p3 = gsqr(gel(x,3));
    2175          70 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2176             : 
    2177          70 :       if (gequal0(gel(p1,3)))
    2178             :       {
    2179           7 :         tetpil = avma;
    2180           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2181           7 :         av = avma;
    2182           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2183           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2184             :       }
    2185             : 
    2186          63 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2187          63 :       tetpil = avma;
    2188          63 :       gel(z,2) = gadd(p2,p4);
    2189          63 :       gel(z,3) = gadd(p1,p3);
    2190          63 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2191             : 
    2192       46880 :     case t_POLMOD:
    2193       46880 :       return sqr_polmod(gel(x,1), gel(x,2));
    2194             : 
    2195     2933018 :     case t_FFELT: return FF_sqr(x);
    2196             : 
    2197      756897 :     case t_POL: return RgX_sqr(x);
    2198             : 
    2199       27552 :     case t_SER:
    2200       27552 :       lx = lg(x);
    2201       27552 :       if (ser_isexactzero(x)) {
    2202          21 :         GEN z = gcopy(x);
    2203          21 :         setvalp(z, 2*valp(x));
    2204          21 :         return z;
    2205             :       }
    2206       27531 :       if (lx < 40)
    2207       27272 :         return normalizeser( sqr_ser_part(x, 0, lx-3) );
    2208             :       else
    2209             :       {
    2210         259 :         pari_sp av = avma;
    2211         259 :         GEN z = init_ser(lx, varn(x), 2*valp(x));
    2212         259 :         x = ser2pol_i(x, lx);
    2213         259 :         x = RgXn_sqr(x, lx-2);
    2214         259 :         return gerepilecopy(av, fill_ser(z,x));
    2215             :       }
    2216             : 
    2217          14 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2218          14 :       gel(z,1) = gsqr(gel(x,1));
    2219          14 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2220             : 
    2221         980 :     case t_MAT: return RgM_sqr(x);
    2222       95617 :     case t_VEC: /* handle extended t_QFB */
    2223       95617 :     case t_QFB: return qfbsqr(x);
    2224         658 :     case t_VECSMALL:
    2225         658 :       z = cgetg_copy(x, &lx);
    2226       16289 :       for (i=1; i<lx; i++)
    2227             :       {
    2228       15631 :         long xi = x[i];
    2229       15631 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2230       15631 :         z[i] = x[xi];
    2231             :       }
    2232         658 :       return z;
    2233             :   }
    2234           6 :   pari_err_TYPE2("*",x,x);
    2235             :   return NULL; /* LCOV_EXCL_LINE */
    2236             : }
    2237             : 
    2238             : /********************************************************************/
    2239             : /**                                                                **/
    2240             : /**                           DIVISION                             **/
    2241             : /**                                                                **/
    2242             : /********************************************************************/
    2243             : static GEN
    2244       30832 : div_rfrac_scal(GEN x, GEN y)
    2245             : {
    2246       30832 :   pari_sp av = avma;
    2247       30832 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2248       30832 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2249             : }
    2250             : static GEN
    2251       37459 : div_scal_rfrac(GEN x, GEN y)
    2252             : {
    2253       37459 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2254       37459 :   pari_sp av = avma;
    2255       37459 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2256             :   {
    2257           7 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2258           0 :     y1 = gel(y1,2);
    2259             :   }
    2260       37452 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2261             : }
    2262             : static GEN
    2263     1186635 : div_rfrac(GEN x, GEN y)
    2264     1186635 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2265             : 
    2266             : /* x != 0 */
    2267             : static GEN
    2268     1290013 : div_ser_scal(GEN y, GEN x)
    2269             : {
    2270             :   long i, l;
    2271             :   GEN z;
    2272     1290013 :   if (ser_isexactzero(y))
    2273             :   {
    2274          21 :     z = scalarser(lg(y) == 2? Rg_get_0(x): gdiv(gel(y,2), x), varn(y), 1);
    2275          21 :     setvalp(z, valp(y)); return z;
    2276             :   }
    2277     1289992 :   z = cgetg_copy(y, &l); z[1] = y[1];
    2278     5913412 :   for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
    2279     1289992 :   return normalizeser(z);
    2280             : }
    2281             : GEN
    2282        7553 : ser_normalize(GEN x)
    2283             : {
    2284        7553 :   long i, lx = lg(x);
    2285             :   GEN c, z;
    2286        7553 :   if (lx == 2) return x;
    2287        7553 :   c = gel(x,2); if (gequal1(c)) return x;
    2288        7476 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2289      107835 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2290        7476 :   return z;
    2291             : }
    2292             : 
    2293             : /* y != 0 */
    2294             : static GEN
    2295     4819083 : div_T_scal(GEN x, GEN y, long tx) {
    2296     4819083 :   switch(tx)
    2297             :   {
    2298     3502112 :     case t_POL: return RgX_Rg_div(x, y);
    2299     1290006 :     case t_SER: return div_ser_scal(x, y);
    2300       26968 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2301             :   }
    2302           0 :   pari_err_TYPE2("/",x,y);
    2303             :   return NULL; /* LCOV_EXCL_LINE */
    2304             : }
    2305             : 
    2306             : static GEN
    2307     9257631 : div_scal_pol(GEN x, GEN y) {
    2308     9257631 :   long ly = lg(y);
    2309             :   pari_sp av;
    2310     9257631 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2311     9179356 :   if (isrationalzero(x)) return zeropol(varn(y));
    2312     7130116 :   av = avma;
    2313     7130116 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2314             : }
    2315             : static GEN
    2316       18326 : div_scal_ser(GEN x, GEN y)
    2317       18326 : { pari_sp av = avma; return gerepileupto(av, gmul(x, ser_inv(y))); }
    2318             : static GEN
    2319     9266194 : div_scal_T(GEN x, GEN y, long ty) {
    2320     9266194 :   switch(ty)
    2321             :   {
    2322     9211172 :     case t_POL: return div_scal_pol(x, y);
    2323       18319 :     case t_SER: return div_scal_ser(x, y);
    2324       36703 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2325             :   }
    2326           0 :   pari_err_TYPE2("/",x,y);
    2327             :   return NULL; /* LCOV_EXCL_LINE */
    2328             : }
    2329             : 
    2330             : /* assume tx = ty = t_SER, same variable vx */
    2331             : static GEN
    2332      753667 : div_ser(GEN x, GEN y, long vx)
    2333             : {
    2334      753667 :   long e, v = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2335      753667 :   GEN y0 = y, z;
    2336      753667 :   pari_sp av = avma;
    2337             : 
    2338      753667 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2339      753660 :   if (ser_isexactzero(x))
    2340             :   {
    2341       59934 :     if (lx == 2) return zeroser(vx, v);
    2342         147 :     z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
    2343         147 :     setvalp(z, v); return z;
    2344             :   }
    2345      693726 :   if (lx < ly) ly = lx;
    2346      693726 :   y = ser2pol_approx(y, ly, &e);
    2347      693726 :   if (e) { v -= e; ly -= e; if (ly <= 2) pari_err_INV("div_ser", y0); }
    2348      693726 :   z = init_ser(ly, vx, v);
    2349      693726 :   if (ly == 3)
    2350             :   {
    2351       15658 :     gel(z,2) = gdiv(gel(x,2), gel(y,2));
    2352       15658 :     return gerepileupto(av, z);
    2353             :   }
    2354      678068 :   x = ser2pol_i(x, ly);
    2355      678068 :   y = RgXn_div_i(x, y, ly-2);
    2356      678068 :   return gerepilecopy(av, fill_ser(z,y));
    2357             : }
    2358             : /* x,y compatible PADIC */
    2359             : static GEN
    2360     2814678 : divpp(GEN x, GEN y) {
    2361             :   pari_sp av;
    2362             :   long a, b;
    2363             :   GEN z, M;
    2364             : 
    2365     2814678 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2366     2814680 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2367     2814365 :   a = precp(x);
    2368     2814365 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2369     2814365 :   z = cgetg(5, t_PADIC);
    2370     2814069 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2371     2813971 :   gel(z,2) = icopy(gel(x,2));
    2372     2813794 :   gel(z,3) = icopy(M); av = avma;
    2373     2813885 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2374     2814639 :   return z;
    2375             : }
    2376             : static GEN
    2377       27405 : div_polmod_same(GEN T, GEN x, GEN y)
    2378             : {
    2379       27405 :   long v = varn(T);
    2380       27405 :   GEN a, z = cgetg(3, t_POLMOD);
    2381       27405 :   gel(z,1) = RgX_copy(T);
    2382       27405 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2383        9772 :     a = gdiv(x, y);
    2384       17633 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2385         119 :   {
    2386         119 :     pari_sp av = avma;
    2387         119 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2388             :   }
    2389       17514 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2390        5145 :   {
    2391        5145 :     pari_sp av = avma;
    2392        5145 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2393        5145 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2394        5145 :     a = gerepileupto(av, a);
    2395             :   }
    2396             :   else
    2397             :   {
    2398       12369 :     pari_sp av = avma;
    2399       12369 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2400       12369 :     a = gerepileupto(av, a);
    2401             :   }
    2402       27405 :   gel(z,2) = a; return z;
    2403             : }
    2404             : GEN
    2405   374042008 : gdiv(GEN x, GEN y)
    2406             : {
    2407   374042008 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2408             :   pari_sp av, tetpil;
    2409             :   GEN z, p1, p2;
    2410             : 
    2411   374042008 :   if (tx == ty) switch(tx)
    2412             :   {
    2413    78044212 :     case t_INT:
    2414    78044212 :       return Qdivii(x,y);
    2415             : 
    2416    81952067 :     case t_REAL: return divrr(x,y);
    2417       21083 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2418       21083 :       z = cgetg(3,t_INTMOD);
    2419       21083 :       if (X==Y || equalii(X,Y))
    2420       21076 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2421           7 :       gel(z,1) = gcdii(X,Y);
    2422           7 :       warn_coercion(X,Y,gel(z,1));
    2423           7 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2424           7 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2425             :     }
    2426     3683729 :     case t_FRAC: {
    2427     3683729 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2428     3683729 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2429     3683729 :       z = cgetg(3, t_FRAC);
    2430     3683729 :       p1 = gcdii(x1, y1);
    2431     3683728 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2432     3683728 :       p1 = gcdii(x2, y2);
    2433     3683728 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2434     3683724 :       tetpil = avma;
    2435     3683724 :       gel(z,2) = mulii(x2,y1);
    2436     3683723 :       gel(z,1) = mulii(x1,y2);
    2437     3683726 :       normalize_frac(z);
    2438     3683727 :       fix_frac_if_int_GC(z,tetpil);
    2439     3683729 :       return z;
    2440             :     }
    2441     7697710 :     case t_COMPLEX:
    2442     7697710 :       if (isintzero(gel(y,1)))
    2443             :       {
    2444      110362 :         y = gel(y,2);
    2445      110362 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2446       93359 :         z = cgetg(3,t_COMPLEX);
    2447       93359 :         gel(z,1) = gdiv(gel(x,2), y);
    2448       93359 :         av = avma;
    2449       93359 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2450       93359 :         return z;
    2451             :       }
    2452     7587346 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2453     7587414 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2454             : 
    2455      767638 :     case t_PADIC:
    2456      767638 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2457      767638 :       return divpp(x, y);
    2458             : 
    2459         238 :     case t_QUAD:
    2460         238 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2461         224 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2462         224 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2463             : 
    2464      236484 :     case t_FFELT: return FF_div(x,y);
    2465             : 
    2466       31780 :     case t_POLMOD:
    2467       31780 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2468       27405 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2469             :       else {
    2470        4375 :         av = avma;
    2471        4375 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2472             :       }
    2473       31780 :       return z;
    2474             : 
    2475    18459397 :     case t_POL:
    2476    18459397 :       vx = varn(x);
    2477    18459397 :       vy = varn(y);
    2478    18459397 :       if (vx != vy) {
    2479      102557 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2480       46459 :                             else return div_scal_pol(x, y);
    2481             :       }
    2482    18356840 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2483    18356840 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2484    18233428 :       av = avma;
    2485    18233428 :       return gerepileupto(av, gred_rfrac2(x,y));
    2486             : 
    2487      753688 :     case t_SER:
    2488      753688 :       vx = varn(x);
    2489      753688 :       vy = varn(y);
    2490      753688 :       if (vx != vy) {
    2491          21 :         if (varncmp(vx, vy) < 0)
    2492             :         {
    2493          14 :           if (!signe(y)) pari_err_INV("gdiv",y);
    2494           7 :           return div_ser_scal(x, y);
    2495             :         }
    2496           7 :         return div_scal_ser(x, y);
    2497             :       }
    2498      753667 :       return div_ser(x, y, vx);
    2499     1191115 :     case t_RFRAC:
    2500     1191115 :       vx = varn(gel(x,2));
    2501     1191115 :       vy = varn(gel(y,2));
    2502     1191115 :       if (vx != vy) {
    2503        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2504         756 :                             else return div_scal_rfrac(x, y);
    2505             :       }
    2506     1186635 :       return div_rfrac(x,y);
    2507             : 
    2508          21 :     case t_VEC: /* handle extended t_QFB */
    2509          21 :     case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
    2510             : 
    2511          14 :     case t_MAT:
    2512          14 :       av = avma; p1 = RgM_inv(y);
    2513          14 :       if (!p1) pari_err_INV("gdiv",y);
    2514          14 :       return gerepileupto(av, RgM_mul(x, p1));
    2515             : 
    2516           0 :     default: pari_err_TYPE2("/",x,y);
    2517             :   }
    2518             : 
    2519   181208208 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2520             :   {
    2521     3845062 :     long s = signe(x);
    2522     3845062 :     if (!s) {
    2523      476184 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2524      476184 :       switch (ty)
    2525             :       {
    2526      473027 :       default: return gen_0;
    2527          21 :       case t_INTMOD:
    2528          21 :         z = cgetg(3,t_INTMOD);
    2529          21 :         gel(z,1) = icopy(gel(y,1));
    2530          21 :         gel(z,2) = gen_0; return z;
    2531        3136 :       case t_FFELT: return FF_zero(y);
    2532             :       }
    2533             :     }
    2534     3368878 :     if (is_pm1(x)) {
    2535     1095008 :       if (s > 0) return ginv(y);
    2536      189782 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2537             :     }
    2538     2273870 :     switch(ty)
    2539             :     {
    2540      798458 :       case t_REAL: return divir(x,y);
    2541          21 :       case t_INTMOD:
    2542          21 :         z = cgetg(3, t_INTMOD);
    2543          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2544      780681 :       case t_FRAC:
    2545      780681 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2546      780681 :         if (equali1(p1))
    2547             :         {
    2548      249857 :           set_avma((pari_sp)z);
    2549      249857 :           gel(z,2) = icopy(gel(y,1));
    2550      249857 :           gel(z,1) = mulii(gel(y,2), x);
    2551      249857 :           normalize_frac(z);
    2552      249857 :           fix_frac_if_int(z);
    2553             :         }
    2554             :         else
    2555             :         {
    2556      530824 :           x = diviiexact(x,p1); tetpil = avma;
    2557      530824 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2558      530824 :           gel(z,1) = mulii(gel(y,2), x);
    2559      530824 :           normalize_frac(z);
    2560      530824 :           fix_frac_if_int_GC(z,tetpil);
    2561             :         }
    2562      780681 :         return z;
    2563             : 
    2564         273 :       case t_FFELT: return Z_FF_div(x,y);
    2565      692713 :       case t_COMPLEX: return divRc(x,y);
    2566        1505 :       case t_PADIC: return divTp(x, y);
    2567         231 :       case t_QUAD:
    2568         231 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2569         231 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2570             :     }
    2571             :   }
    2572   177363133 :   if (gequal0(y))
    2573             :   {
    2574          49 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2575          28 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2576             :   }
    2577             : 
    2578   177373852 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2579             :   {
    2580   142413540 :     case t_REAL:
    2581   142413540 :       switch(ty)
    2582             :       {
    2583   140087005 :         case t_INT: return divri(x,y);
    2584      507168 :         case t_FRAC:
    2585      507168 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2586      507168 :           return gerepileuptoleaf(av, z);
    2587     1819430 :         case t_COMPLEX: return divRc(x, y);
    2588          42 :         case t_QUAD: return divfq(x, y, realprec(x));
    2589          18 :         default: pari_err_TYPE2("/",x,y);
    2590             :       }
    2591             : 
    2592         287 :     case t_INTMOD:
    2593             :       switch(ty)
    2594             :       {
    2595         203 :         case t_INT:
    2596         203 :           z = cgetg(3, t_INTMOD);
    2597         203 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2598          28 :         case t_FRAC: { GEN X = gel(x,1);
    2599          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2600          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2601             :         }
    2602           7 :         case t_FFELT:
    2603           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2604           0 :             pari_err_OP("/",x,y);
    2605           7 :           return Z_FF_div(gel(x,2),y);
    2606             : 
    2607           0 :         case t_COMPLEX:
    2608           0 :           av = avma;
    2609           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2610             : 
    2611          14 :         case t_QUAD:
    2612          14 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2613          14 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2614             : 
    2615           7 :         case t_PADIC: { GEN X = gel(x,1);
    2616           7 :           z = cgetg(3, t_INTMOD);
    2617           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2618             :         }
    2619          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2620             :       }
    2621             : 
    2622             :     case t_FRAC:
    2623             :       switch(ty)
    2624             :       {
    2625     1998009 :         case t_INT: z = cgetg(3, t_FRAC);
    2626     1998009 :         p1 = gcdii(y,gel(x,1));
    2627     1998009 :         if (equali1(p1))
    2628             :         {
    2629      755923 :           set_avma((pari_sp)z); tetpil = 0;
    2630      755923 :           gel(z,1) = icopy(gel(x,1));
    2631             :         }
    2632             :         else
    2633             :         {
    2634     1242086 :           y = diviiexact(y,p1); tetpil = avma;
    2635     1242086 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2636             :         }
    2637     1998009 :         gel(z,2) = mulii(gel(x,2),y);
    2638     1998009 :         normalize_frac(z);
    2639     1998009 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2640     1998009 :         return z;
    2641             : 
    2642       57745 :         case t_REAL:
    2643       57745 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2644       57745 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2645             : 
    2646           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2647           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2648           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2649             :         }
    2650             : 
    2651          28 :         case t_FFELT: av=avma;
    2652          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2653             : 
    2654         483 :         case t_COMPLEX: return divRc(x, y);
    2655             : 
    2656        2128 :         case t_PADIC:
    2657        2128 :           if (!signe(gel(x,1))) return gen_0;
    2658        2128 :           return divTp(x, y);
    2659             : 
    2660          42 :         case t_QUAD:
    2661          42 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2662          42 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2663             :       }
    2664             : 
    2665             :     case t_FFELT:
    2666         133 :       switch (ty)
    2667             :       {
    2668          49 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2669          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2670           7 :         case t_INTMOD:
    2671           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2672           0 :             pari_err_OP("/",x,y);
    2673           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2674          49 :         default:
    2675          49 :         pari_err_TYPE2("/",x,y);
    2676             :       }
    2677           0 :       break;
    2678             : 
    2679    11942807 :     case t_COMPLEX:
    2680             :       switch(ty)
    2681             :       {
    2682    11942788 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2683           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2684           0 :         case t_PADIC:
    2685           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2686           0 :         case t_QUAD:
    2687           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2688           0 :           return divfq(x, y, lx);
    2689             :       }
    2690             : 
    2691             :     case t_PADIC:
    2692             :       switch(ty)
    2693             :       {
    2694     1327291 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2695     1327200 :           return signe(gel(x,4))? divpT(x, y)
    2696     2654491 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2697             :         }
    2698           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2699           7 :           z = cgetg(3, t_INTMOD);
    2700           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2701             :         }
    2702          14 :         case t_COMPLEX: case t_QUAD:
    2703          14 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2704          14 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2705             : 
    2706          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2707             :       }
    2708             : 
    2709             :     case t_QUAD:
    2710             :       switch (ty)
    2711             :       {
    2712        1190 :         case t_INT: case t_INTMOD: case t_FRAC:
    2713        1190 :           z = cgetg(4,t_QUAD);
    2714        1190 :           gel(z,1) = ZX_copy(gel(x,1));
    2715        1190 :           gel(z,2) = gdiv(gel(x,2), y);
    2716        1190 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2717          63 :         case t_REAL: return divqf(x, y, realprec(y));
    2718          14 :         case t_PADIC: return divTp(x, y);
    2719           0 :         case t_COMPLEX:
    2720           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2721           0 :           return divqf(x, y, ly);
    2722             :       }
    2723             :   }
    2724    19629865 :   switch(ty) {
    2725      374309 :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2726      374309 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2727          42 :     case t_MAT:
    2728          42 :       av = avma; p1 = RgM_inv(y);
    2729          35 :       if (!p1) pari_err_INV("gdiv",y);
    2730          28 :       return gerepileupto(av, gmul(x, p1));
    2731           0 :     case t_VEC: case t_COL:
    2732             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2733           0 :       pari_err_TYPE2("/",x,y);
    2734             :   }
    2735    19255514 :   switch(tx) {
    2736     3169904 :     case t_VEC: case t_COL: case t_MAT:
    2737     3169904 :       z = cgetg_copy(x, &lx);
    2738    11081150 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2739     3169798 :       return z;
    2740           0 :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2741           0 :       pari_err_TYPE2("/",x,y);
    2742             :   }
    2743             : 
    2744    16085610 :   vy = gvar(y);
    2745    16086422 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2746      207456 :     vx = varn(X);
    2747      207456 :     if (vx != vy) {
    2748      206882 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2749      206553 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2750             :     }
    2751             :     /* y is POL, SER or RFRAC */
    2752         574 :     av = avma;
    2753         574 :     switch(ty)
    2754             :     {
    2755           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2756         574 :       default: y = ginvmod(gmod(y,X), X);
    2757             :     }
    2758         567 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2759             :   }
    2760             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2761             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2762             :    * variable NO_VARIABLE, then the operation is incorrect */
    2763    15878966 :   vx = gvar(x);
    2764    15878961 :   if (vx != vy) { /* includes cases where one is scalar */
    2765    14084948 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2766     9265863 :                         else return div_scal_T(x, y, ty);
    2767             :   }
    2768     1794013 :   switch(tx)
    2769             :   {
    2770     1046143 :     case t_POL:
    2771             :       switch(ty)
    2772             :       {
    2773          28 :         case t_SER:
    2774             :         {
    2775          28 :           GEN y0 = y;
    2776             :           long v;
    2777          28 :           ly = lg(y); /* > 2 */
    2778          28 :           av = avma; v = RgX_valrem(x, &x);
    2779          28 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2780          14 :           v -= valp(y);
    2781          14 :           y = ser2pol_approx(y, ly, &i); if (i) { ly -= i; v -= i; }
    2782          14 :           if (ly == 2) pari_err_INV("gdiv", y0);
    2783           7 :           z = init_ser(ly, vx, v);
    2784           7 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
    2785             :         }
    2786             : 
    2787     1046115 :         case t_RFRAC:
    2788             :         {
    2789     1046115 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2790     1046115 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2791          35 :             return mul_rfrac_scal(y2, y1, x);
    2792     1046080 :           av = avma;
    2793     1046080 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2794             :         }
    2795             :       }
    2796           0 :       break;
    2797             : 
    2798      478161 :     case t_SER:
    2799             :       switch(ty)
    2800             :       {
    2801      478147 :         case t_POL:
    2802             :         {
    2803      478147 :           long v = valp(x);
    2804      478147 :           lx = lg(x);
    2805      478147 :           if (lx == 2) return zeroser(vx, v - RgX_val(y));
    2806      478042 :           av = avma;
    2807      478042 :           x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
    2808      478042 :           z = init_ser(lx, vx, v);
    2809      478042 :           if (!signe(x)) setsigne(z,0);
    2810      478042 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
    2811             :         }
    2812          14 :         case t_RFRAC:
    2813          14 :           av = avma;
    2814          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2815             :       }
    2816           0 :       break;
    2817             : 
    2818      269687 :     case t_RFRAC:
    2819             :       switch(ty)
    2820             :       {
    2821      269687 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2822           0 :         case t_SER:
    2823           0 :           av = avma;
    2824           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
    2825             :       }
    2826           0 :       break;
    2827             :   }
    2828          22 :   pari_err_TYPE2("/",x,y);
    2829             :   return NULL; /* LCOV_EXCL_LINE */
    2830             : }
    2831             : 
    2832             : /********************************************************************/
    2833             : /**                                                                **/
    2834             : /**                     SIMPLE MULTIPLICATION                      **/
    2835             : /**                                                                **/
    2836             : /********************************************************************/
    2837             : GEN
    2838   246693993 : gmulsg(long s, GEN y)
    2839             : {
    2840             :   long ly, i;
    2841             :   pari_sp av;
    2842             :   GEN z;
    2843             : 
    2844   246693993 :   switch(typ(y))
    2845             :   {
    2846   156937086 :     case t_INT:  return mulsi(s,y);
    2847    69151131 :     case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
    2848      174905 :     case t_INTMOD: { GEN p = gel(y,1);
    2849      174905 :       z = cgetg(3,t_INTMOD);
    2850      174907 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2851      174907 :       gel(z,1) = icopy(p); return z;
    2852             :     }
    2853      547745 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2854     5292967 :     case t_FRAC:
    2855     5292967 :       if (!s) return gen_0;
    2856     5288053 :       z = cgetg(3,t_FRAC);
    2857     5288053 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2858     5288053 :       if (i == 1)
    2859             :       {
    2860     2021662 :         gel(z,2) = icopy(gel(y,2));
    2861     2021662 :         gel(z,1) = mulis(gel(y,1), s);
    2862             :       }
    2863             :       else
    2864             :       {
    2865     3266391 :         gel(z,2) = diviuexact(gel(y,2), (ulong)i);
    2866     3266391 :         gel(z,1) = mulis(gel(y,1), s/i);
    2867     3266391 :         fix_frac_if_int(z);
    2868             :       }
    2869     5288053 :       return z;
    2870             : 
    2871    13102396 :     case t_COMPLEX:
    2872    13102396 :       if (!s) return gen_0;
    2873    13028048 :       z = cgetg(3, t_COMPLEX);
    2874    13027861 :       gel(z,1) = gmulsg(s,gel(y,1));
    2875    13025930 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2876             : 
    2877        1435 :     case t_PADIC:
    2878        1435 :       if (!s) return gen_0;
    2879        1435 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2880             : 
    2881           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2882           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2883           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2884           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2885             : 
    2886       74802 :     case t_POLMOD:
    2887       74802 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2888             : 
    2889      585875 :     case t_POL:
    2890      585875 :       if (!signe(y)) return RgX_copy(y);
    2891      572295 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2892      570671 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2893     2342076 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2894      570671 :       return normalizepol_lg(z, ly);
    2895             : 
    2896         182 :     case t_SER:
    2897         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2898         182 :       if (!s) return Rg_get_0(y);
    2899         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2900        3864 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2901         182 :       return normalizeser(z);
    2902             : 
    2903          14 :     case t_RFRAC:
    2904          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2905          14 :       if (s == 1) return gcopy(y);
    2906          14 :       if (s == -1) return gneg(y);
    2907          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2908             : 
    2909      837166 :     case t_VEC: case t_COL: case t_MAT:
    2910      837166 :       z = cgetg_copy(y, &ly);
    2911     2480754 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2912      837154 :       return z;
    2913             :   }
    2914           0 :   pari_err_TYPE("gmulsg",y);
    2915             :   return NULL; /* LCOV_EXCL_LINE */
    2916             : }
    2917             : 
    2918             : GEN
    2919   132364991 : gmulug(ulong s, GEN y)
    2920             : {
    2921             :   long ly, i;
    2922             :   pari_sp av;
    2923             :   GEN z;
    2924             : 
    2925   132364991 :   switch(typ(y))
    2926             :   {
    2927   130456383 :     case t_INT:  return mului(s,y);
    2928     1046319 :     case t_REAL: return s? mulur(s,y): gen_0; /* gmul semantic */
    2929         364 :     case t_INTMOD: { GEN p = gel(y,1);
    2930         364 :       z = cgetg(3,t_INTMOD);
    2931         364 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(y,2)), p));
    2932         364 :       gel(z,1) = icopy(p); return z;
    2933             :     }
    2934         413 :     case t_FFELT: return FF_Z_mul(y,utoi(s));
    2935      802877 :     case t_FRAC:
    2936      802877 :       if (!s) return gen_0;
    2937      802863 :       z = cgetg(3,t_FRAC);
    2938      802863 :       i = ugcd(s, umodiu(gel(y,2), s));
    2939      802863 :       if (i == 1)
    2940             :       {
    2941      665972 :         gel(z,2) = icopy(gel(y,2));
    2942      665972 :         gel(z,1) = muliu(gel(y,1), s);
    2943             :       }
    2944             :       else
    2945             :       {
    2946      136891 :         gel(z,2) = diviuexact(gel(y,2), i);
    2947      136891 :         gel(z,1) = muliu(gel(y,1), s/i);
    2948      136891 :         fix_frac_if_int(z);
    2949             :       }
    2950      802863 :       return z;
    2951             : 
    2952       24750 :     case t_COMPLEX:
    2953       24750 :       if (!s) return gen_0;
    2954       24750 :       z = cgetg(3, t_COMPLEX);
    2955       24750 :       gel(z,1) = gmulug(s,gel(y,1));
    2956       24750 :       gel(z,2) = gmulug(s,gel(y,2)); return z;
    2957             : 
    2958        3898 :     case t_PADIC:
    2959        3898 :       if (!s) return gen_0;
    2960        3898 :       av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),y), y));
    2961             : 
    2962           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2963           0 :       gel(z,1) = ZX_copy(gel(y,1));
    2964           0 :       gel(z,2) = gmulug(s,gel(y,2));
    2965           0 :       gel(z,3) = gmulug(s,gel(y,3)); return z;
    2966             : 
    2967        3199 :     case t_POLMOD:
    2968        3199 :       retmkpolmod(gmulug(s,gel(y,2)), RgX_copy(gel(y,1)));
    2969             : 
    2970       12642 :     case t_POL:
    2971       12642 :       if (!signe(y)) return RgX_copy(y);
    2972       11865 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2973       11865 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2974       34923 :       for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    2975       11865 :       return normalizepol_lg(z, ly);
    2976             : 
    2977           0 :     case t_SER:
    2978           0 :       if (ser_isexactzero(y)) return gcopy(y);
    2979           0 :       if (!s) return Rg_get_0(y);
    2980           0 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2981           0 :       for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    2982           0 :       return normalizeser(z);
    2983             : 
    2984           0 :     case t_RFRAC:
    2985           0 :       if (!s) return zeropol(varn(gel(y,2)));
    2986           0 :       if (s == 1) return gcopy(y);
    2987           0 :       return mul_rfrac_scal(gel(y,1), gel(y,2), utoi(s));
    2988             : 
    2989       14147 :     case t_VEC: case t_COL: case t_MAT:
    2990       14147 :       z = cgetg_copy(y, &ly);
    2991       73922 :       for (i=1; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    2992       14147 :       return z;
    2993             :   }
    2994           0 :   pari_err_TYPE("gmulsg",y);
    2995             :   return NULL; /* LCOV_EXCL_LINE */
    2996             : }
    2997             : 
    2998             : /********************************************************************/
    2999             : /**                                                                **/
    3000             : /**                       SIMPLE DIVISION                          **/
    3001             : /**                                                                **/
    3002             : /********************************************************************/
    3003             : 
    3004             : GEN
    3005     7219500 : gdivgs(GEN x, long s)
    3006             : {
    3007     7219500 :   long tx = typ(x), lx, i;
    3008             :   pari_sp av;
    3009             :   GEN z;
    3010             : 
    3011     7219500 :   if (!s)
    3012             :   {
    3013           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3014           0 :     pari_err_INV("gdivgs",gen_0);
    3015             :   }
    3016     7219517 :   switch(tx)
    3017             :   {
    3018     1326121 :     case t_INT: return Qdivis(x, s);
    3019     4690896 :     case t_REAL: return divrs(x,s);
    3020             : 
    3021         357 :     case t_INTMOD:
    3022         357 :       z = cgetg(3, t_INTMOD);
    3023         357 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    3024             : 
    3025         735 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    3026             : 
    3027      383917 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3028      383917 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    3029      383917 :       if (i == 1)
    3030             :       {
    3031      294265 :         gel(z,2) = mulsi(s, gel(x,2));
    3032      294265 :         gel(z,1) = icopy(gel(x,1));
    3033             :       }
    3034             :       else
    3035             :       {
    3036       89652 :         gel(z,2) = mulsi(s/i, gel(x,2));
    3037       89652 :         gel(z,1) = divis(gel(x,1), i);
    3038             :       }
    3039      383917 :       normalize_frac(z);
    3040      383917 :       fix_frac_if_int(z); return z;
    3041             : 
    3042      756423 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3043      756423 :       gel(z,1) = gdivgs(gel(x,1),s);
    3044      756421 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3045             : 
    3046         133 :     case t_PADIC: /* divpT */
    3047             :     {
    3048         133 :       GEN p = gel(x,2);
    3049         133 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3050         133 :       av = avma;
    3051         133 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    3052             :     }
    3053             : 
    3054          28 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3055          28 :       gel(z,1) = ZX_copy(gel(x,1));
    3056          28 :       gel(z,2) = gdivgs(gel(x,2),s);
    3057          28 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    3058             : 
    3059       29659 :     case t_POLMOD:
    3060       29659 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    3061             : 
    3062          91 :     case t_RFRAC:
    3063          91 :       if (s == 1) return gcopy(x);
    3064          84 :       else if (s == -1) return gneg(x);
    3065          84 :       return div_rfrac_scal(x, stoi(s));
    3066             : 
    3067       29442 :     case t_POL: case t_SER:
    3068       29442 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3069      117191 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3070       29442 :       return z;
    3071        1715 :     case t_VEC: case t_COL: case t_MAT:
    3072        1715 :       z = cgetg_copy(x, &lx);
    3073       40971 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3074        1715 :       return z;
    3075             : 
    3076             :   }
    3077           0 :   pari_err_TYPE2("/",x, stoi(s));
    3078             :   return NULL; /* LCOV_EXCL_LINE */
    3079             : }
    3080             : 
    3081             : GEN
    3082    42784391 : gdivgu(GEN x, ulong s)
    3083             : {
    3084    42784391 :   long tx = typ(x), lx, i;
    3085             :   pari_sp av;
    3086             :   GEN z;
    3087             : 
    3088    42784391 :   if (!s)
    3089             :   {
    3090           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3091           0 :     pari_err_INV("gdivgu",gen_0);
    3092             :   }
    3093    42784327 :   switch(tx)
    3094             :   {
    3095    16943589 :     case t_INT: return Qdiviu(x, s);
    3096    12889670 :     case t_REAL: return divru(x,s);
    3097             : 
    3098      210315 :     case t_INTMOD:
    3099      210315 :       z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
    3100      210315 :       return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
    3101             : 
    3102         308 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
    3103             : 
    3104      280325 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3105      280325 :       i = ugcd(s, umodiu(gel(x,1), s));
    3106      280325 :       if (i == 1)
    3107             :       {
    3108      148028 :         gel(z,2) = mului(s, gel(x,2));
    3109      148028 :         gel(z,1) = icopy(gel(x,1));
    3110             :       }
    3111             :       else
    3112             :       {
    3113      132297 :         gel(z,2) = mului(s/i, gel(x,2));
    3114      132297 :         gel(z,1) = divis(gel(x,1), i);
    3115             :       }
    3116      280325 :       normalize_frac(z);
    3117      280325 :       fix_frac_if_int(z); return z;
    3118             : 
    3119    11556147 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3120    11556148 :       gel(z,1) = gdivgu(gel(x,1),s);
    3121    11556146 :       gel(z,2) = gdivgu(gel(x,2),s); return z;
    3122             : 
    3123      850395 :     case t_PADIC: /* divpT */
    3124             :     {
    3125      850395 :       GEN p = gel(x,2);
    3126      850395 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3127      716264 :       av = avma;
    3128      716264 :       return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
    3129             :     }
    3130             : 
    3131           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3132           0 :       gel(z,1) = ZX_copy(gel(x,1));
    3133           0 :       gel(z,2) = gdivgu(gel(x,2),s);
    3134           0 :       gel(z,3) = gdivgu(gel(x,3),s); return z;
    3135             : 
    3136        1456 :     case t_POLMOD:
    3137        1456 :       retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
    3138             : 
    3139          56 :     case t_RFRAC:
    3140          56 :       if (s == 1) return gcopy(x);
    3141          56 :       return div_rfrac_scal(x, utoi(s));
    3142             : 
    3143       51765 :     case t_POL: case t_SER:
    3144       51765 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3145      228326 :       for (i=2; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
    3146       51765 :       return z;
    3147         301 :     case t_VEC: case t_COL: case t_MAT:
    3148         301 :       z = cgetg_copy(x, &lx);
    3149        1092 :       for (i=1; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
    3150         301 :       return z;
    3151             :   }
    3152           0 :   pari_err_TYPE2("/",x, utoi(s));
    3153             :   return NULL; /* LCOV_EXCL_LINE */
    3154             : }
    3155             : 
    3156             : /* x / (i*(i+1)) */
    3157             : GEN
    3158   105279624 : divrunextu(GEN x, ulong i)
    3159             : {
    3160   105279624 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3161           0 :     return divri(x, muluu(i , i+1));
    3162             :   else
    3163   105279624 :     return divru(x, i*(i+1));
    3164             : }
    3165             : /* x / (i*(i+1)) */
    3166             : GEN
    3167      804801 : gdivgunextu(GEN x, ulong i)
    3168             : {
    3169      804801 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3170           0 :     return gdivgu(x, i*(i+1));
    3171             :   else
    3172      804801 :     return gdiv(x, muluu(i, i+1));
    3173             : }
    3174             : 
    3175             : /* True shift (exact multiplication by 2^n) */
    3176             : GEN
    3177    85971815 : gmul2n(GEN x, long n)
    3178             : {
    3179             :   long lx, i, k, l;
    3180             :   GEN z, a, b;
    3181             : 
    3182    85971815 :   switch(typ(x))
    3183             :   {
    3184    19111446 :     case t_INT:
    3185    19111446 :       if (n>=0) return shifti(x,n);
    3186     3383056 :       if (!signe(x)) return gen_0;
    3187     3291636 :       l = vali(x); n = -n;
    3188     3291721 :       if (n<=l) return shifti(x,-n);
    3189      375044 :       z = cgetg(3,t_FRAC);
    3190      375044 :       gel(z,1) = shifti(x,-l);
    3191      375044 :       gel(z,2) = int2n(n-l); return z;
    3192             : 
    3193    49073329 :     case t_REAL:
    3194    49073329 :       return shiftr(x,n);
    3195             : 
    3196      180800 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3197      180800 :       z = cgetg(3,t_INTMOD);
    3198      180800 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3199       82639 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3200       82642 :       gel(z,1) = icopy(b); return z;
    3201             : 
    3202      217234 :     case t_FFELT: return FF_mul2n(x,n);
    3203             : 
    3204      839637 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3205      839637 :       l = vali(a);
    3206      839637 :       k = vali(b);
    3207      839637 :       if (n+l >= k)
    3208             :       {
    3209      299999 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3210      236641 :         l = n-k; k = -k;
    3211             :       }
    3212             :       else
    3213             :       {
    3214      539638 :         k = -(l+n); l = -l;
    3215             :       }
    3216      776279 :       z = cgetg(3,t_FRAC);
    3217      776279 :       gel(z,1) = shifti(a,l);
    3218      776279 :       gel(z,2) = shifti(b,k); return z;
    3219             : 
    3220    12132834 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3221    12132833 :       gel(z,1) = gmul2n(gel(x,1),n);
    3222    12132836 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3223             : 
    3224         105 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3225         105 :       gel(z,1) = ZX_copy(gel(x,1));
    3226         105 :       gel(z,2) = gmul2n(gel(x,2),n);
    3227         105 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3228             : 
    3229       44030 :     case t_POLMOD:
    3230       44030 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3231             : 
    3232     1578729 :     case t_POL:
    3233     1578729 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3234     8586714 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3235     1578600 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3236      106158 :     case t_SER:
    3237      106158 :       if (ser_isexactzero(x)) return gcopy(x);
    3238      106130 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3239      642327 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3240      106130 :       return normalizeser(z); /* needed if char = 2 */
    3241     2684422 :     case t_VEC: case t_COL: case t_MAT:
    3242     2684422 :       z = cgetg_copy(x, &lx);
    3243    12155509 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3244     2684316 :       return z;
    3245             : 
    3246          21 :     case t_RFRAC: /* int2n wrong if n < 0 */
    3247          21 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3248             : 
    3249        4536 :     case t_PADIC: /* int2n wrong if n < 0 */
    3250        4536 :       return gmul(gmul2n(gen_1,n),x);
    3251             :   }
    3252           0 :   pari_err_TYPE("gmul2n",x);
    3253             :   return NULL; /* LCOV_EXCL_LINE */
    3254             : }
    3255             : 
    3256             : /*******************************************************************/
    3257             : /*                                                                 */
    3258             : /*                              INVERSE                            */
    3259             : /*                                                                 */
    3260             : /*******************************************************************/
    3261             : static GEN
    3262       79175 : inv_polmod(GEN T, GEN x)
    3263             : {
    3264       79175 :   GEN z = cgetg(3,t_POLMOD), a;
    3265       79173 :   gel(z,1) = RgX_copy(T);
    3266       79167 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3267       62251 :     a = ginv(x);
    3268             :   else
    3269             :   {
    3270       16916 :     if (lg(T) == 5) /* quadratic fields */
    3271       12705 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3272             :     else
    3273        4211 :       a = RgXQ_inv(x, T);
    3274             :   }
    3275       79167 :   gel(z,2) = a; return z;
    3276             : }
    3277             : 
    3278             : GEN
    3279    21599985 : ginv(GEN x)
    3280             : {
    3281             :   long s;
    3282             :   pari_sp av, tetpil;
    3283             :   GEN z, y, p1, p2;
    3284             : 
    3285    21599985 :   switch(typ(x))
    3286             :   {
    3287     3360888 :     case t_INT:
    3288     3360888 :       if (is_pm1(x)) return icopy(x);
    3289     2528513 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3290     2528499 :       z = cgetg(3,t_FRAC);
    3291     2528534 :       gel(z,1) = s<0? gen_m1: gen_1;
    3292     2528534 :       gel(z,2) = absi(x); return z;
    3293             : 
    3294     5265353 :     case t_REAL: return invr(x);
    3295             : 
    3296        4641 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3297        4641 :       gel(z,1) = icopy(gel(x,1));
    3298        4641 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3299             : 
    3300      386251 :     case t_FRAC: {
    3301      386251 :       GEN a = gel(x,1), b = gel(x,2);
    3302      386251 :       s = signe(a);
    3303      386251 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3304      205765 :       z = cgetg(3,t_FRAC);
    3305      205765 :       gel(z,1) = icopy(b);
    3306      205765 :       gel(z,2) = icopy(a);
    3307      205765 :       normalize_frac(z); return z;
    3308             :     }
    3309     6666599 :     case t_COMPLEX:
    3310     6666599 :       av=avma;
    3311     6666599 :       p1=cxnorm(x);
    3312     6665764 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3313     6666475 :       tetpil=avma;
    3314     6666475 :       return gerepile(av,tetpil,divcR(p2,p1));
    3315             : 
    3316         273 :     case t_QUAD:
    3317         273 :       av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
    3318         273 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3319             : 
    3320      346213 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3321      346213 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3322      346206 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3323      346206 :       gel(z,2) = icopy(gel(x,2));
    3324      346206 :       gel(z,3) = icopy(gel(x,3));
    3325      346206 :       gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
    3326             : 
    3327       79176 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3328       14269 :     case t_FFELT: return FF_inv(x);
    3329     5321021 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3330       18858 :     case t_SER: return ser_inv(x);
    3331        2926 :     case t_RFRAC:
    3332             :     {
    3333        2926 :       GEN n = gel(x,1), d = gel(x,2);
    3334        2926 :       pari_sp av = avma, ltop;
    3335        2926 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3336             : 
    3337        2926 :       n = simplify_shallow(n);
    3338        2926 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3339             :       {
    3340        2926 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3341         679 :         ltop = avma;
    3342         679 :         z = RgX_Rg_div(d,n);
    3343             :       } else {
    3344           0 :         ltop = avma;
    3345           0 :         z = cgetg(3,t_RFRAC);
    3346           0 :         gel(z,1) = RgX_copy(d);
    3347           0 :         gel(z,2) = RgX_copy(n);
    3348             :       }
    3349         679 :       stackdummy(av, ltop);
    3350         679 :       return z;
    3351             :     }
    3352             : 
    3353       99119 :     case t_VEC: /* handle extended t_QFB */
    3354             :     case t_QFB:
    3355       99119 :       return qfbpow(x, gen_m1);
    3356       34496 :     case t_MAT:
    3357       34496 :       y = RgM_inv(x);
    3358       34489 :       if (!y) pari_err_INV("ginv",x);
    3359       34419 :       return y;
    3360          28 :     case t_VECSMALL:
    3361             :     {
    3362          28 :       long i, lx = lg(x)-1;
    3363          28 :       y = zero_zv(lx);
    3364         112 :       for (i=1; i<=lx; i++)
    3365             :       {
    3366          84 :         long xi = x[i];
    3367          84 :         if (xi<1 || xi>lx || y[xi])
    3368           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3369          84 :         y[xi] = i;
    3370             :       }
    3371          28 :       return y;
    3372             :     }
    3373             :   }
    3374           0 :   pari_err_TYPE("inverse",x);
    3375             :   return NULL; /* LCOV_EXCL_LINE */
    3376             : }

Generated by: LCOV version 1.13