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.16.2 lcov report (development 29464-65cd88daf0) Lines: 1936 2061 93.9 %
Date: 2024-07-23 09:03:50 Functions: 101 101 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         105 : warn_coercion(GEN x, GEN y, GEN z)
      38             : {
      39         105 :   if (DEBUGLEVEL)
      40          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      41         105 : }
      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   304454256 : addRc(GEN x, GEN y) {
      61   304454256 :   GEN z = cgetg(3,t_COMPLEX);
      62   304436556 :   gel(z,1) = gadd(x,gel(y,1));
      63   304418064 :   gel(z,2) = gcopy(gel(y,2)); return z;
      64             : }
      65             : static GEN
      66   351803418 : mulRc(GEN x, GEN y) {
      67   351803418 :   GEN z = cgetg(3,t_COMPLEX);
      68   351704564 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      69   352223026 :   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     2349732 : divRc(GEN x, GEN y) {
      80     2349732 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      81     2349728 :   GEN z = cgetg(3,t_COMPLEX);
      82     2349728 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      83     2349716 :   gel(z,2) = gmul(mt, gel(y,2));
      84     2349717 :   return z;
      85             : }
      86             : static GEN
      87    23100653 : divcR(GEN x, GEN y) {
      88    23100653 :   GEN z = cgetg(3,t_COMPLEX);
      89    23100588 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      90    23100252 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      91             : }
      92             : static GEN
      93        1274 : addRq(GEN x, GEN y) {
      94        1274 :   GEN z = cgetg(4,t_QUAD);
      95        1274 :   gel(z,1) = ZX_copy(gel(y,1));
      96        1274 :   gel(z,2) = gadd(x, gel(y,2));
      97        1274 :   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          77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     108          77 :   long i = gexpo(x) - gexpo(y);
     109          77 :   if (i > 0) prec += nbits2extraprec(i);
     110          77 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     111             : }
     112             : static GEN
     113    37482974 : mulrfrac(GEN x, GEN y)
     114             : {
     115             :   pari_sp av;
     116    37482974 :   GEN z, a = gel(y,1), b = gel(y,2);
     117    37482974 :   if (is_pm1(a)) /* frequent special case */
     118             :   {
     119    12687368 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     120    12688020 :     return z;
     121             :   }
     122    24795572 :   av = avma;
     123    24795572 :   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     4981668 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     151     4981668 :   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        3660 : divTp(GEN x, GEN y) { pari_sp av = avma;
     156        3660 :   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     1203979 : divpT(GEN x, GEN y) { pari_sp av = avma;
     162     1203979 :   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     3690618 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     169     3690618 :   if (lgefint(X) == 3) {
     170     3367474 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     171     3367474 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     172             :   }
     173             :   else {
     174      323144 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     175      323142 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     176             :   }
     177     3690616 :   gel(z,1) = icopy(X); return z;
     178             : }
     179             : static GEN
     180     1158553 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     181     1158553 :   if (lgefint(X) == 3) {
     182      784519 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     183      784519 :     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     1158553 :   gel(z,1) = icopy(X); return z;
     190             : }
     191             : /* cf add_intmod_same */
     192             : static GEN
     193     3103269 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     194     3103269 :   if (lgefint(X) == 3) {
     195     2961833 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     196     2961833 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     197             :   }
     198             :   else
     199      141436 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     200     3103277 :   gel(z,1) = icopy(X); return z;
     201             : }
     202             : /* cf add_intmod_same */
     203             : static GEN
     204      341803 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     205             : {
     206      341803 :   if (lgefint(X) == 3) {
     207      319286 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     208      319279 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     209             :   }
     210             :   else
     211       22517 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     212      341798 :   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     9857076 : rfrac_denom_mul_scal(GEN d, GEN y)
     225             : {
     226     9857076 :   GEN D = RgX_Rg_mul(d, y);
     227     9857076 :   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     9857076 :   return D;
     233             : }
     234             : 
     235             : static int
     236    57828812 : Leading_is_neg(GEN x)
     237             : {
     238   122277863 :   while (typ(x) == t_POL) x = leading_coeff(x);
     239    57828812 :   return is_real_t(typ(x))? gsigne(x) < 0: 0;
     240             : }
     241             : 
     242             : static int
     243   154479832 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
     244             : 
     245             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     246             : GEN
     247    57833139 : gred_rfrac_simple(GEN n, GEN d)
     248             : {
     249             :   GEN _1n, _1d, c, cn, cd, z;
     250    57833139 :   long dd = degpol(d);
     251             : 
     252    57833139 :   if (dd <= 0)
     253             :   {
     254        4327 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     255        4327 :     n = gdiv(n, gel(d,2));
     256        4327 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     257        4327 :     return n;
     258             :   }
     259    57828812 :   if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
     260    57828812 :   _1n = Rg_get_1(n);
     261    57828812 :   _1d = Rg_get_1(d);
     262    57828812 :   if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
     263    57828812 :   if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
     264    57828805 :   cd = content(d);
     265    59699245 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     266    57828805 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     267    57828805 :   if (!gequal1(cd)) {
     268     6568436 :     d = RgX_Rg_div(d,cd);
     269     6568436 :     if (!gequal1(cn))
     270             :     {
     271     1294462 :       if (gequal0(cn)) {
     272          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     273           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     274           0 :         c = gen_1;
     275             :       } else {
     276     1294413 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     277     1294413 :         c = gdiv(cn,cd);
     278             :       }
     279             :     }
     280             :     else
     281     5273974 :       c = ginv(cd);
     282             :   } else {
     283    51260369 :     if (!gequal1(cn))
     284             :     {
     285     3276983 :       if (gequal0(cn)) {
     286        1484 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     287          21 :         c = gen_1;
     288             :       } else {
     289     3275499 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     290     3275499 :         c = cn;
     291             :       }
     292             :     } else {
     293    47983386 :       GEN y = cgetg(3,t_RFRAC);
     294    47983386 :       gel(y,1) = gcopy(n);
     295    47983386 :       gel(y,2) = RgX_copy(d); return y;
     296             :     }
     297             :   }
     298             : 
     299     9843907 :   if (typ(c) == t_POL)
     300             :   {
     301      900394 :     z = c;
     302      938754 :     do { z = content(z); } while (typ(z) == t_POL);
     303      900394 :     cd = denom_i(z);
     304      900394 :     cn = gmul(c, cd);
     305             :   }
     306             :   else
     307             :   {
     308     8943513 :     cn = numer_i(c);
     309     8943513 :     cd = denom_i(c);
     310             :   }
     311     9843907 :   z = cgetg(3,t_RFRAC);
     312     9843907 :   gel(z,1) = gmul(n, cn);
     313     9843907 :   gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
     314             :   /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
     315     9843907 :   if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
     316     9843907 :   return z;
     317             : }
     318             : 
     319             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     320             : static GEN
     321      248325 : fix_rfrac(GEN x, long d)
     322             : {
     323             :   GEN z, N, D;
     324      248325 :   if (!d || typ(x) == t_POL) return x;
     325      171689 :   z = cgetg(3, t_RFRAC);
     326      171689 :   N = gel(x,1);
     327      171689 :   D = gel(x,2);
     328      171689 :   if (d > 0) {
     329        8701 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     330      165900 :                                                    : monomialcopy(N,d,varn(D));
     331      165851 :     gel(z, 2) = RgX_copy(D);
     332             :   } else {
     333        5838 :     gel(z, 1) = gcopy(N);
     334        5838 :     gel(z, 2) = RgX_shift(D, -d);
     335             :   }
     336      171689 :   return z;
     337             : }
     338             : 
     339             : /* assume d != 0 */
     340             : static GEN
     341    44766507 : gred_rfrac2(GEN n, GEN d)
     342             : {
     343             :   GEN y, z, _1n, _1d;
     344             :   long v, vd, vn;
     345             : 
     346    44766507 :   n = simplify_shallow(n);
     347    44766507 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     348    37937846 :   d = simplify_shallow(d);
     349    37937846 :   if (typ(d) != t_POL) return gdiv(n,d);
     350    36759172 :   vd = varn(d);
     351    36759172 :   if (typ(n) != t_POL)
     352             :   {
     353    20624450 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     354    20623043 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     355           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     356             :   }
     357    16134722 :   vn = varn(n);
     358    16134722 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     359    15997062 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     360    15837310 :   _1n = Rg_get_1(n);
     361    15837310 :   _1d = Rg_get_1(d);
     362    15837310 :   if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
     363    15837310 :   if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
     364             : 
     365             :   /* now n and d are t_POLs in the same variable */
     366    15837310 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     367    15837310 :   if (!degpol(d))
     368             :   {
     369    12757608 :     n = RgX_Rg_div(n,gel(d,2));
     370    12757608 :     return v? RgX_mulXn(n,v): n;
     371             :   }
     372             : 
     373             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     374     3079702 :   if (!isinexact(n) && !isinexact(d))
     375             :   {
     376     3079457 :     y = RgX_divrem(n, d, &z);
     377     3079457 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     378      248080 :     z = RgX_gcd(d, z);
     379      248080 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     380             :   }
     381      248325 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     382             : }
     383             : 
     384             : /* x,y t_INT, return x/y in reduced form */
     385             : GEN
     386   130639751 : Qdivii(GEN x, GEN y)
     387             : {
     388   130639751 :   pari_sp av = avma;
     389             :   GEN r, q;
     390             : 
     391   130639751 :   if (lgefint(y) == 3)
     392             :   {
     393   113817933 :     q = Qdiviu(x, y[2]);
     394   113815478 :     if (signe(y) > 0) return q;
     395    10841047 :     if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
     396    10841286 :     return q;
     397             :   }
     398    16821818 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     399    16823231 :   if (equali1(x))
     400             :   {
     401     5177897 :     if (!signe(y)) pari_err_INV("gdiv",y);
     402     5177785 :     retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
     403             :   }
     404    11645357 :   q = dvmdii(x,y,&r);
     405    11645273 :   if (r == gen_0) return q; /* gen_0 intended */
     406     8168631 :   r = gcdii(y, r);
     407     8168660 :   if (lgefint(r) == 3)
     408             :   {
     409     7301051 :     ulong t = r[2];
     410     7301051 :     set_avma(av);
     411     7301026 :     if (t == 1) q = mkfraccopy(x,y);
     412             :     else
     413             :     {
     414     2561366 :       q = cgetg(3,t_FRAC);
     415     2561485 :       gel(q,1) = diviuexact(x,t);
     416     2561403 :       gel(q,2) = diviuexact(y,t);
     417             :     }
     418             :   }
     419             :   else
     420             :   { /* rare: r and q left on stack for efficiency */
     421      867609 :     q = cgetg(3,t_FRAC);
     422      867615 :     gel(q,1) = diviiexact(x,r);
     423      867615 :     gel(q,2) = diviiexact(y,r);
     424             :   }
     425     8168585 :   normalize_frac(q); return q;
     426             : }
     427             : 
     428             : /* x t_INT, return x/y in reduced form */
     429             : GEN
     430   133125543 : Qdiviu(GEN x, ulong y)
     431             : {
     432   133125543 :   pari_sp av = avma;
     433             :   ulong r, t;
     434             :   GEN q;
     435             : 
     436   133125543 :   if (y == 1) return icopy(x);
     437   112839759 :   if (!y) pari_err_INV("gdiv",gen_0);
     438   112839759 :   if (equali1(x)) retmkfrac(gen_1, utoipos(y));
     439   105739898 :   q = absdiviu_rem(x,y,&r);
     440   105736489 :   if (!r)
     441             :   {
     442    60749668 :     if (signe(x) < 0) togglesign(q);
     443    60749629 :     return q;
     444             :   }
     445    44986821 :   t = ugcd(y, r); set_avma(av);
     446    44991613 :   if (t == 1) retmkfrac(icopy(x), utoipos(y));
     447    18132977 :   retmkfrac(diviuexact(x,t), utoipos(y / t));
     448             : }
     449             : 
     450             : /* x t_INT, return x/y in reduced form */
     451             : GEN
     452     1560284 : Qdivis(GEN x, long y)
     453             : {
     454     1560284 :   pari_sp av = avma;
     455             :   ulong r, t;
     456             :   long s;
     457             :   GEN q;
     458             : 
     459     1560284 :   if (y > 0) return Qdiviu(x, y);
     460      163110 :   if (!y) pari_err_INV("gdiv",gen_0);
     461      163110 :   s = signe(x);
     462      163110 :   if (!s) return gen_0;
     463      132118 :   if (y < 0) { y = -y; s = -s; }
     464      132118 :   if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
     465      131838 :   if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
     466      129696 :   q = absdiviu_rem(x,y,&r);
     467      129696 :   if (!r)
     468             :   {
     469       55328 :     if (s < 0) togglesign(q);
     470       55328 :     return q;
     471             :   }
     472       74368 :   t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
     473       74368 :   if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
     474       74368 :   gel(q,1) = x; setsigne(x, s);
     475       74368 :   gel(q,2) = utoipos(y); return q;
     476             : }
     477             : 
     478             : /*******************************************************************/
     479             : /*                                                                 */
     480             : /*                          CONJUGATION                            */
     481             : /*                                                                 */
     482             : /*******************************************************************/
     483             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     484             : static GEN
     485       18319 : quad_polmod_conj(GEN x, GEN y)
     486             : {
     487             :   GEN z, u, v, a, b;
     488       18319 :   if (typ(x) != t_POL) return gcopy(x);
     489       18319 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     490       18319 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     491       18319 :   b = gel(y,3); v = gel(x,2);
     492       18319 :   z = cgetg(4, t_POL); z[1] = x[1];
     493       18319 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     494       18319 :   gel(z,3) = gneg(u); return z;
     495             : }
     496             : static GEN
     497       18319 : quad_polmod_norm(GEN x, GEN y)
     498             : {
     499             :   GEN z, u, v, a, b, c;
     500       18319 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     501       18319 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     502       18319 :   b = gel(y,3); v = gel(x,2);
     503       18319 :   c = gel(y,2);
     504       18319 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     505       18319 :   if (!gequal1(a)) z = gdiv(z, a);
     506       18319 :   return gadd(z, gsqr(v));
     507             : }
     508             : 
     509             : GEN
     510    19340014 : conj_i(GEN x)
     511             : {
     512             :   long lx, i;
     513             :   GEN y;
     514             : 
     515    19340014 :   switch(typ(x))
     516             :   {
     517     6411731 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     518     6411731 :       return x;
     519             : 
     520    12765716 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     521             : 
     522         924 :     case t_QUAD:
     523         924 :       y = cgetg(4,t_QUAD);
     524         924 :       gel(y,1) = gel(x,1);
     525         924 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     526         924 :                                       : gadd(gel(x,2), gel(x,3));
     527         924 :       gel(y,3) = gneg(gel(x,3)); return y;
     528             : 
     529        9926 :     case t_POL: case t_SER:
     530        9926 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     531       31913 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     532        9926 :       return y;
     533             : 
     534      151732 :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     535      151732 :       y = cgetg_copy(x, &lx);
     536      546420 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     537      151732 :       return y;
     538             : 
     539           0 :     case t_POLMOD:
     540             :     {
     541           0 :       GEN X = gel(x,1);
     542           0 :       long d = degpol(X);
     543           0 :       if (d < 2) return x;
     544           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     545             :     }
     546             :   }
     547           0 :   pari_err_TYPE("gconj",x);
     548             :   return NULL; /* LCOV_EXCL_LINE */
     549             : }
     550             : GEN
     551      132038 : gconj(GEN x)
     552      132038 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     553             : 
     554             : GEN
     555          84 : conjvec(GEN x,long prec)
     556             : {
     557             :   long lx, s, i;
     558             :   GEN z;
     559             : 
     560          84 :   switch(typ(x))
     561             :   {
     562           0 :     case t_INT: case t_INTMOD: case t_FRAC:
     563           0 :       return mkcolcopy(x);
     564             : 
     565           0 :     case t_COMPLEX: case t_QUAD:
     566           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     567             : 
     568          28 :     case t_FFELT:
     569          28 :       return FF_conjvec(x);
     570             : 
     571           0 :     case t_VEC: case t_COL:
     572           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     573           0 :       if (lx == 1) return z;
     574           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     575           0 :       s = lgcols(z);
     576           0 :       for (i=2; i<lx; i++)
     577             :       {
     578           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     579           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     580             :       }
     581           0 :       break;
     582             : 
     583          56 :     case t_POLMOD: {
     584          56 :       GEN T = gel(x,1), r;
     585             :       pari_sp av;
     586             : 
     587          56 :       lx = lg(T);
     588          56 :       if (lx <= 3) return cgetg(1,t_COL);
     589          56 :       x = gel(x,2);
     590         238 :       for (i=2; i<lx; i++)
     591             :       {
     592         189 :         GEN c = gel(T,i);
     593         189 :         switch(typ(c)) {
     594           7 :           case t_INTMOD: {
     595           7 :             GEN p = gel(c,1);
     596             :             pari_sp av;
     597           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     598           7 :             av = avma;
     599           7 :             T = RgX_to_FpX(T,p);
     600           7 :             x = RgX_to_FpX(x, p);
     601           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     602           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     603           7 :             return gerepileupto(av, z);
     604             :           }
     605         182 :           case t_INT:
     606         182 :           case t_FRAC: break;
     607           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     608             :         }
     609             :       }
     610          49 :       if (typ(x) != t_POL)
     611             :       {
     612           0 :         if (!is_rational_t(typ(x)))
     613           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     614           0 :         retconst_col(lx-3, gcopy(x));
     615             :       }
     616          49 :       RgX_check_QX(x,"conjvec");
     617          49 :       av = avma;
     618          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     619          49 :       r = cleanroots(T,prec);
     620          49 :       z = cgetg(lx-2,t_COL);
     621         182 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     622          49 :       return gerepileupto(av, z);
     623             :     }
     624             : 
     625           0 :     default:
     626           0 :       pari_err_TYPE("conjvec",x);
     627             :       return NULL; /* LCOV_EXCL_LINE */
     628             :   }
     629           0 :   return z;
     630             : }
     631             : 
     632             : /********************************************************************/
     633             : /**                                                                **/
     634             : /**                           ADDITION                             **/
     635             : /**                                                                **/
     636             : /********************************************************************/
     637             : /* x, y compatible PADIC, op = add or sub */
     638             : static GEN
     639    19954762 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     640             : {
     641    19954762 :   pari_sp av = avma;
     642             :   long d,e,r,rx,ry;
     643    19954762 :   GEN u, z, mod, p = gel(x,2);
     644             :   int swap;
     645             : 
     646    19954762 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     647    19954317 :   e = valp(x);
     648    19954317 :   r = valp(y); d = r-e;
     649    19954317 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     650    19954317 :   rx = precp(x);
     651    19954317 :   ry = precp(y);
     652    19954317 :   if (d) /* v(x) < v(y) */
     653             :   {
     654    10689291 :     r = d+ry; z = powiu(p,d);
     655    10689376 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     656    10689404 :     z = mulii(z,gel(y,4));
     657    10689374 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     658             :   }
     659             :   else
     660             :   {
     661             :     long c;
     662     9265026 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     663     9265026 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     664     9265567 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     665             :     {
     666      138276 :       set_avma(av); return zeropadic(p, e+r);
     667             :     }
     668     9127261 :     if (c)
     669             :     {
     670     3351391 :       mod = diviiexact(mod, powiu(p,c));
     671     3351390 :       r -= c;
     672     3351390 :       e += c;
     673             :     }
     674             :   }
     675    19816579 :   u = modii(u, mod);
     676    19815547 :   set_avma(av); z = cgetg(5,t_PADIC);
     677    19814759 :   z[1] = evalprecp(r) | evalvalp(e);
     678    19814597 :   gel(z,2) = icopy(p);
     679    19814292 :   gel(z,3) = icopy(mod);
     680    19814068 :   gel(z,4) = icopy(u); return z;
     681             : }
     682             : /* Rg_to_Fp(t_FRAC) without GC */
     683             : static GEN
     684      226838 : Q_to_Fp(GEN x, GEN p)
     685      226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     686             : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
     687             : static GEN
     688     4760443 : addQp(GEN x, GEN y)
     689             : {
     690     4760443 :   pari_sp av = avma;
     691     4760443 :   long d, r, e, vy = valp(y), py = precp(y);
     692     4760443 :   GEN z, q, mod, u, p = gel(y,2);
     693             : 
     694     4760443 :   e = Q_pvalrem(x, p, &x);
     695     4760431 :   d = vy - e; r = d + py;
     696     4760431 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     697     4758471 :   mod = gel(y,3);
     698     4758471 :   u   = gel(y,4);
     699     4758471 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     700             : 
     701     4758460 :   if (d > 0)
     702             :   {
     703     1375881 :     q = powiu(p,d);
     704     1375964 :     mod = mulii(mod, q);
     705     1375963 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     706     1375963 :     u = addii(x,  mulii(u, q));
     707             :   }
     708     3382579 :   else if (d < 0)
     709             :   {
     710      405227 :     q = powiu(p,-d);
     711      405227 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     712      405227 :     u = addii(u, mulii(x, q));
     713      405227 :     r = py; e = vy;
     714             :   }
     715             :   else
     716             :   {
     717             :     long c;
     718     2977352 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     719     2977352 :     u = addii(u, x);
     720     2977360 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     721             :     {
     722        1205 :       set_avma(av); return zeropadic(p,e+r);
     723             :     }
     724     2976151 :     if (c)
     725             :     {
     726      970539 :       mod = diviiexact(mod, powiu(p,c));
     727      970538 :       r -= c;
     728      970538 :       e += c;
     729             :     }
     730             :   }
     731     4757347 :   u = modii(u, mod); set_avma(av);
     732     4757313 :   z = cgetg(5,t_PADIC);
     733     4757278 :   z[1] = evalprecp(r) | evalvalp(e);
     734     4757269 :   gel(z,2) = icopy(p);
     735     4757256 :   gel(z,3) = icopy(mod);
     736     4757242 :   gel(z,4) = icopy(u); return z;
     737             : }
     738             : 
     739             : /* Mod(x,X) + Mod(y,X) */
     740             : #define addsub_polmod_same addsub_polmod_scal
     741             : /* Mod(x,X) +/- Mod(y,Y) */
     742             : static GEN
     743        2142 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     744             : {
     745        2142 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     746        2142 :   GEN z = cgetg(3,t_POLMOD);
     747        2142 :   long vx = varn(X), vy = varn(Y);
     748        2142 :   if (vx==vy) {
     749             :     pari_sp av;
     750          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     751          14 :     warn_coercion(X,Y,gel(z,1));
     752          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     753             :   }
     754        2128 :   if (varncmp(vx, vy) < 0)
     755        2128 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     756             :   else
     757           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     758        2128 :   gel(z,2) = op(x, y); return z;
     759             : }
     760             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     761             : static GEN
     762    13516461 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     763    13516461 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     764             : 
     765             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     766             : static GEN
     767      408045 : add_ser_scal(GEN y, GEN x)
     768             : {
     769             :   long i, v, ly, vy;
     770             :   GEN z;
     771             : 
     772      408045 :   if (isrationalzero(x)) return gcopy(y);
     773      381802 :   ly = lg(y);
     774      381802 :   v = valser(y);
     775      381802 :   if (v < 3-ly) return gcopy(y);
     776             :   /* v + ly >= 3 */
     777      381550 :   if (v < 0)
     778             :   {
     779        1162 :     z = cgetg(ly,t_SER); z[1] = y[1];
     780        3276 :     for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
     781        1162 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     782        3157 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     783        1162 :     return normalizeser(z);
     784             :   }
     785      380388 :   vy = varn(y);
     786      380388 :   if (v > 0)
     787             :   {
     788       19278 :     if (ser_isexactzero(y))
     789        9310 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
     790        9968 :     y -= v; ly += v;
     791        9968 :     z = cgetg(ly,t_SER);
     792        9968 :     x = gcopy(x);
     793       20566 :     for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
     794             :   }
     795      361110 :   else if (ser_isexactzero(y)) return gcopy(y);
     796             :   else
     797             :   { /* v = 0, ly >= 3 */
     798      361103 :     z = cgetg(ly,t_SER);
     799      361103 :     x = gadd(x, gel(y,2));
     800      361103 :     i = 3;
     801             :   }
     802     1598979 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     803      371071 :   gel(z,2) = x;
     804      371071 :   z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
     805      371071 :   return gequal0(x)? normalizeser(z): z;
     806             : }
     807             : static long
     808     7161278 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     809             : /* x,y t_SER in the same variable: x+y */
     810             : static GEN
     811     3581038 : ser_add(GEN x, GEN y)
     812             : {
     813     3581038 :   long i, lx,ly, n = valser(y) - valser(x);
     814             :   GEN z;
     815     3581038 :   if (n < 0) { n = -n; swap(x,y); }
     816             :   /* valser(x) <= valser(y) */
     817     3581038 :   lx = _serprec(x);
     818     3581038 :   if (lx == 2) /* don't lose type information */
     819             :   {
     820         798 :     z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
     821         798 :     setvalser(z, valser(x)); return z;
     822             :   }
     823     3580240 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     824     3580240 :   if (n)
     825             :   {
     826      107432 :     if (n+2 > lx) return gcopy(x);
     827      106732 :     z = cgetg(ly,t_SER);
     828      817310 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     829      504972 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     830             :   } else {
     831     3472808 :     z = cgetg(ly,t_SER);
     832    17589354 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     833             :   }
     834     3579540 :   z[1] = x[1]; return normalizeser(z);
     835             : }
     836             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     837             : static GEN
     838     8830994 : add_rfrac_scal(GEN y, GEN x)
     839             : {
     840             :   pari_sp av;
     841             :   GEN n;
     842             : 
     843     8830994 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     844     5157187 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     845     5157187 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     846             : }
     847             : 
     848             : /* x "scalar", ty != t_MAT and nonscalar */
     849             : static GEN
     850    22904001 : add_scal(GEN y, GEN x, long ty)
     851             : {
     852    22904001 :   switch(ty)
     853             :   {
     854    17991981 :     case t_POL: return RgX_Rg_add(y, x);
     855      408017 :     case t_SER: return add_ser_scal(y, x);
     856     4468813 :     case t_RFRAC: return add_rfrac_scal(y, x);
     857           0 :     case t_COL: return RgC_Rg_add(y, x);
     858       35179 :     case t_VEC:
     859       35179 :       if (isintzero(x)) return gcopy(y);
     860         182 :       break;
     861             :   }
     862         193 :   pari_err_TYPE2("+",x,y);
     863             :   return NULL; /* LCOV_EXCL_LINE */
     864             : }
     865             : 
     866             : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
     867             : static GEN
     868    14819455 : setfrac(GEN z, GEN a, GEN b)
     869             : {
     870    14819455 :   gel(z,1) = icopy_avma(a, (pari_sp)z);
     871    14819442 :   gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
     872    14819521 :   set_avma((pari_sp)gel(z,2)); return z;
     873             : }
     874             : /* z <- a / (b*Q), (Q,a) = 1 */
     875             : static GEN
     876    14037434 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
     877             : {
     878    14037434 :   GEN q = Qdivii(a, b); /* != 0 */
     879    14037503 :   if (typ(q) == t_INT)
     880             :   {
     881     1975701 :     gel(z,1) = gerepileuptoint((pari_sp)Q, q);
     882     1975701 :     gel(z,2) = Q; return z;
     883             :   }
     884    12061802 :   return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
     885             : }
     886             : static GEN
     887    28492667 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     888             : {
     889    28492667 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     890    28492667 :   GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
     891    28492667 :   int s = cmpii(x2, y2);
     892             : 
     893             :   /* common denominator: (x1 op y1) / x2 */
     894    28492667 :   if (!s)
     895             :   {
     896    10142374 :     pari_sp av = avma;
     897    10142374 :     return gerepileupto(av, Qdivii(op(x1, y1), x2));
     898             :   }
     899    18350293 :   z = cgetg(3, t_FRAC);
     900    18350295 :   if (s < 0)
     901             :   {
     902    10414166 :     Q = dvmdii(y2, x2, &r);
     903             :     /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     904    10414127 :     if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
     905     2753005 :     delta = gcdii(x2,r);
     906             :   }
     907             :   else
     908             :   {
     909     7936129 :     Q = dvmdii(x2, y2, &r);
     910             :     /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     911     7936132 :     if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
     912     1559761 :     delta = gcdii(y2,r);
     913             :   }
     914             :   /* delta = gcd(x2,y2) */
     915     4312806 :   if (equali1(delta))
     916             :   { /* numerator is nonzero */
     917     1555079 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     918     1555079 :     gel(z,2) = mulii(x2,y2); return z;
     919             :   }
     920     2757722 :   x2 = diviiexact(x2,delta);
     921     2757722 :   y2 = diviiexact(y2,delta);
     922     2757724 :   n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
     923     2757726 :   q = dvmdii(n, delta, &r);
     924     2757718 :   if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
     925     2503196 :   r = gcdii(delta, r);
     926     2503202 :   if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
     927     2503202 :   return setfrac(z, n, mulii(mulii(x2, y2), delta));
     928             : }
     929             : 
     930             : /* assume x2, y2 are t_POLs in the same variable */
     931             : static GEN
     932     3039495 : add_rfrac(GEN x, GEN y)
     933             : {
     934     3039495 :   pari_sp av = avma;
     935     3039495 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     936     3039495 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     937             : 
     938     3039495 :   delta = RgX_gcd(x2,y2);
     939     3039495 :   if (!degpol(delta))
     940             :   {
     941         672 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     942         672 :     d = RgX_mul(x2, y2);
     943         672 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     944             :   }
     945     3038823 :   x2 = RgX_div(x2,delta);
     946     3038823 :   y2 = RgX_div(y2,delta);
     947     3038823 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     948     3038823 :   if (!signe(n))
     949             :   {
     950      721425 :     n = simplify_shallow(n);
     951      721425 :     if (isexactzero(n))
     952             :     {
     953      721418 :       if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
     954          14 :       return gerepileupto(av, scalarpol(n, varn(x2)));
     955             :     }
     956           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     957             :   }
     958     2317398 :   if (degpol(n) == 0)
     959     1150596 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     960     1166802 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     961     1166802 :   if (isexactzero(r))
     962             :   {
     963             :     GEN z;
     964      228087 :     d = RgX_mul(x2, y2);
     965             :     /* "constant" denominator ? */
     966      228087 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     967      228087 :     return gerepileupto(av, z);
     968             :   }
     969      938715 :   r = RgX_gcd(delta, r);
     970      938715 :   if (degpol(r))
     971             :   {
     972      160740 :     n = RgX_div(n, r);
     973      160740 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     974             :   }
     975             :   else
     976      777975 :     d = RgX_mul(gel(x,2), y2);
     977      938715 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     978             : }
     979             : 
     980             : GEN
     981  5813880945 : gadd(GEN x, GEN y)
     982             : {
     983  5813880945 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     984             :   pari_sp av, tetpil;
     985             :   GEN z, p1;
     986             : 
     987  5813880945 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     988             :   {
     989  2860036279 :     case t_INT: return addii(x,y);
     990  1627795857 :     case t_REAL: return addrr(x,y);
     991     1521939 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     992     1521939 :       z = cgetg(3,t_INTMOD);
     993     1521939 :       if (X==Y || equalii(X,Y))
     994     1521925 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     995          14 :       gel(z,1) = gcdii(X,Y);
     996          14 :       warn_coercion(X,Y,gel(z,1));
     997          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     998          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     999             :     }
    1000    23115606 :     case t_FRAC: return addsub_frac(x,y,addii);
    1001   354174256 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1002   353751806 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
    1003   354319712 :       if (isintzero(gel(z,2)))
    1004             :       {
    1005      521164 :         set_avma((pari_sp)(z+3));
    1006      521164 :         return gadd(gel(x,1),gel(y,1));
    1007             :       }
    1008   353633553 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
    1009   353738296 :       return z;
    1010    14517154 :     case t_PADIC:
    1011    14517154 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1012    14516971 :       return addsub_pp(x,y, addii);
    1013         672 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1014         672 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1015         672 :       gel(z,1) = ZX_copy(gel(x,1));
    1016         672 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
    1017         672 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
    1018    10215538 :     case t_POLMOD:
    1019    10215538 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1020    10213431 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
    1021        2107 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
    1022    14006261 :     case t_FFELT: return FF_add(x,y);
    1023    69254692 :     case t_POL:
    1024    69254692 :       vx = varn(x);
    1025    69254692 :       vy = varn(y);
    1026    69254692 :       if (vx != vy) {
    1027      872004 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
    1028       65574 :         else                     return RgX_Rg_add(y, x);
    1029             :       }
    1030    68382688 :       return RgX_add(x, y);
    1031     3575690 :     case t_SER:
    1032     3575690 :       vx = varn(x);
    1033     3575690 :       vy = varn(y);
    1034     3575690 :       if (vx != vy) {
    1035          28 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1036          21 :         else                     return add_ser_scal(y, x);
    1037             :       }
    1038     3575662 :       return ser_add(x, y);
    1039     4331396 :     case t_RFRAC:
    1040     4331396 :       vx = varn(gel(x,2));
    1041     4331396 :       vy = varn(gel(y,2));
    1042     4331396 :       if (vx != vy) {
    1043     1291901 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1044      538397 :         else                     return add_rfrac_scal(y, x);
    1045             :       }
    1046     3039495 :       return add_rfrac(x,y);
    1047     2055693 :     case t_VEC:
    1048     2055693 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1049     2055693 :       return RgV_add(x,y);
    1050     1100762 :     case t_COL:
    1051     1100762 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1052     1100762 :       return RgC_add(x,y);
    1053      671314 :     case t_MAT:
    1054      671314 :       lx = lg(x);
    1055      671314 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1056      671314 :       if (lx == 1) return cgetg(1, t_MAT);
    1057      671314 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1058      671307 :       return RgM_add(x,y);
    1059             : 
    1060           0 :     default: pari_err_TYPE2("+",x,y);
    1061             :   }
    1062             :   /* tx != ty */
    1063   831447526 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1064             : 
    1065   831447526 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1066             :   {
    1067   738895259 :     case t_INT:
    1068             :       switch(ty)
    1069             :       {
    1070   433797921 :         case t_REAL: return addir(x,y);
    1071     2139209 :         case t_INTMOD:
    1072     2139209 :           z = cgetg(3, t_INTMOD);
    1073     2139209 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1074    45195065 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1075    45195065 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1076    45195003 :           gel(z,2) = icopy(gel(y,2)); return z;
    1077   250852922 :         case t_COMPLEX: return addRc(x, y);
    1078     5594351 :         case t_PADIC:
    1079     5594351 :           if (!signe(x)) return gcopy(y);
    1080     4530703 :           return addQp(x,y);
    1081        1113 :         case t_QUAD: return addRq(x, y);
    1082     1356227 :         case t_FFELT: return FF_Z_add(y,x);
    1083             :       }
    1084             : 
    1085             :     case t_REAL:
    1086    56822001 :       switch(ty)
    1087             :       {
    1088    13500318 :         case t_FRAC:
    1089    13500318 :           if (!signe(gel(y,1))) return rcopy(x);
    1090    13500318 :           if (!signe(x))
    1091             :           {
    1092       11868 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1093       11868 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1094             :           }
    1095    13488450 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1096    13487140 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1097    43321613 :         case t_COMPLEX: return addRc(x, y);
    1098          70 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
    1099             : 
    1100           0 :         default: pari_err_TYPE2("+",x,y);
    1101             :       }
    1102             : 
    1103       17640 :     case t_INTMOD:
    1104             :       switch(ty)
    1105             :       {
    1106       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1107       17507 :           z = cgetg(3, t_INTMOD);
    1108       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1109       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1110             :         }
    1111          14 :         case t_FFELT:
    1112          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1113           0 :             pari_err_OP("+",x,y);
    1114          14 :           return FF_Z_add(y,gel(x,2));
    1115          91 :         case t_COMPLEX: return addRc(x, y);
    1116           0 :         case t_PADIC: { GEN X = gel(x,1);
    1117           0 :           z = cgetg(3, t_INTMOD);
    1118           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1119             :         }
    1120          28 :         case t_QUAD: return addRq(x, y);
    1121             :       }
    1122             : 
    1123             :     case t_FRAC:
    1124             :       switch (ty)
    1125             :       {
    1126    10294038 :         case t_COMPLEX: return addRc(x, y);
    1127      229736 :         case t_PADIC:
    1128      229736 :           if (!signe(gel(x,1))) return gcopy(y);
    1129      229736 :           return addQp(x,y);
    1130         133 :         case t_QUAD: return addRq(x, y);
    1131        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1132             :       }
    1133             : 
    1134             :     case t_FFELT:
    1135           0 :       pari_err_TYPE2("+",x,y);
    1136             : 
    1137          35 :     case t_COMPLEX:
    1138             :       switch(ty)
    1139             :       {
    1140          28 :         case t_PADIC:
    1141          28 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1142           7 :         case t_QUAD:
    1143           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1144           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1145             :       }
    1146             : 
    1147             :     case t_PADIC: /* ty == t_QUAD */
    1148           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1149             :   }
    1150             :   /* tx < ty, !is_const_t(y) */
    1151    28640354 :   switch(ty)
    1152             :   {
    1153        7925 :     case t_MAT:
    1154        7925 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1155        7925 :       if (isrationalzero(x)) return gcopy(y);
    1156        7848 :       return RgM_Rg_add(y, x);
    1157      206624 :     case t_COL:
    1158      206624 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1159      206624 :       return RgC_Rg_add(y, x);
    1160     2405026 :     case t_POLMOD: /* is_const_t(tx) in this case */
    1161     2405026 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1162             :   }
    1163    26020779 :   if (is_scalar_t(tx))  {
    1164    22921888 :     if (tx == t_POLMOD)
    1165             :     {
    1166      107447 :       vx = varn(gel(x,1));
    1167      107447 :       vy = gvar(y);
    1168      107447 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1169             :       else
    1170       81890 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1171       40600 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1172             :     }
    1173    22814441 :     return add_scal(y, x, ty);
    1174             :   }
    1175             :   /* x and y are not scalars, ty != t_MAT */
    1176     3098895 :   vx = gvar(x);
    1177     3098898 :   vy = gvar(y);
    1178     3098898 :   if (vx != vy) { /* x or y is treated as a scalar */
    1179       22724 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1180       32340 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1181       32340 :                                 : add_scal(y, x, ty);
    1182             :   }
    1183             :   /* vx = vy */
    1184     3076174 :   switch(tx)
    1185             :   {
    1186     3075677 :     case t_POL:
    1187             :       switch (ty)
    1188             :       {
    1189        5397 :         case t_SER:
    1190        5397 :           if (lg(x) == 2) return gcopy(y);
    1191        5383 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1192        5383 :           i = lg(y) + valser(y) - i;
    1193        5383 :           if (i < 3) return gcopy(y);
    1194        5376 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1195        5376 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1196        5376 :           return y;
    1197             : 
    1198     3070280 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1199             :       }
    1200           0 :       break;
    1201             : 
    1202         497 :     case t_SER:
    1203         497 :       if (ty == t_RFRAC)
    1204             :       {
    1205             :         long vn, vd;
    1206         497 :         av = avma;
    1207         497 :         vn = gval(gel(y,1), vy);
    1208         497 :         vd = RgX_valrem_inexact(gel(y,2), NULL);
    1209         497 :         if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
    1210             : 
    1211         497 :         l = lg(x) + valser(x) - (vn - vd);
    1212         497 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1213         497 :         return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
    1214             :       }
    1215           0 :       break;
    1216             :   }
    1217           0 :   pari_err_TYPE2("+",x,y);
    1218             :   return NULL; /* LCOV_EXCL_LINE */
    1219             : }
    1220             : 
    1221             : GEN
    1222   270575167 : gaddsg(long x, GEN y)
    1223             : {
    1224   270575167 :   long ty = typ(y);
    1225             :   GEN z;
    1226             : 
    1227   270575167 :   switch(ty)
    1228             :   {
    1229   124827075 :     case t_INT:  return addsi(x,y);
    1230   119199431 :     case t_REAL: return addsr(x,y);
    1231       11921 :     case t_INTMOD:
    1232       11921 :       z = cgetg(3, t_INTMOD);
    1233       11921 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1234    15245311 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1235    15245311 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1236    15245311 :       gel(z,2) = icopy(gel(y,2)); return z;
    1237     8163749 :     case t_COMPLEX:
    1238     8163749 :       z = cgetg(3, t_COMPLEX);
    1239     8163749 :       gel(z,1) = gaddsg(x, gel(y,1));
    1240     8163749 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1241             : 
    1242     3127680 :     default: return gadd(stoi(x), y);
    1243             :   }
    1244             : }
    1245             : 
    1246             : GEN
    1247     3184666 : gsubsg(long x, GEN y)
    1248             : {
    1249             :   GEN z, a, b;
    1250             :   pari_sp av;
    1251             : 
    1252     3184666 :   switch(typ(y))
    1253             :   {
    1254      276273 :     case t_INT:  return subsi(x,y);
    1255     1243095 :     case t_REAL: return subsr(x,y);
    1256          56 :     case t_INTMOD:
    1257          56 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1258          56 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1259      732304 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1260      732304 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1261      732304 :       gel(z,2) = icopy(gel(y,2)); return z;
    1262      895922 :     case t_COMPLEX:
    1263      895922 :       z = cgetg(3, t_COMPLEX);
    1264      895922 :       gel(z,1) = gsubsg(x, gel(y,1));
    1265      895922 :       gel(z,2) = gneg(gel(y,2)); return z;
    1266             :   }
    1267       37016 :   av = avma;
    1268       37016 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1269             : }
    1270             : 
    1271             : /********************************************************************/
    1272             : /**                                                                **/
    1273             : /**                          SUBTRACTION                           **/
    1274             : /**                                                                **/
    1275             : /********************************************************************/
    1276             : 
    1277             : GEN
    1278  2801344700 : gsub(GEN x, GEN y)
    1279             : {
    1280  2801344700 :   long tx = typ(x), ty = typ(y);
    1281             :   pari_sp av;
    1282             :   GEN z;
    1283  2801344700 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1284             :   {
    1285  2044649258 :     case t_INT: return subii(x,y);
    1286   569209632 :     case t_REAL: return subrr(x,y);
    1287     1158567 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1288     1158567 :       z = cgetg(3,t_INTMOD);
    1289     1158567 :       if (X==Y || equalii(X,Y))
    1290     1158553 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1291          14 :       gel(z,1) = gcdii(X,Y);
    1292          14 :       warn_coercion(X,Y,gel(z,1));
    1293          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1294          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1295             :     }
    1296     5377066 :     case t_FRAC: return addsub_frac(x,y, subii);
    1297   101659304 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1298   101692464 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1299   101607772 :       if (isintzero(gel(z,2)))
    1300             :       {
    1301       21224 :         set_avma((pari_sp)(z+3));
    1302       21224 :         return gsub(gel(x,1),gel(y,1));
    1303             :       }
    1304   101588780 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1305   101576493 :       return z;
    1306     5437800 :     case t_PADIC:
    1307     5437800 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1308     5437800 :       return addsub_pp(x,y, subii);
    1309         168 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1310         168 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1311         168 :       gel(z,1) = ZX_copy(gel(x,1));
    1312         168 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1313         168 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1314      857439 :     case t_POLMOD:
    1315      857439 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1316      857404 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1317          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1318      203219 :     case t_FFELT: return FF_sub(x,y);
    1319    10316965 :     case t_POL: {
    1320    10316965 :       long vx = varn(x);
    1321    10316965 :       long vy = varn(y);
    1322    10316965 :       if (vx != vy) {
    1323       29615 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1324        5600 :         else                     return Rg_RgX_sub(x, y);
    1325             :       }
    1326    10287350 :       return RgX_sub(x, y);
    1327             :     }
    1328      298793 :     case t_VEC:
    1329      298793 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1330      298793 :       return RgV_sub(x,y);
    1331     3365034 :     case t_COL:
    1332     3365034 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1333     3365034 :       return RgC_sub(x,y);
    1334      216440 :     case t_MAT: {
    1335      216440 :       long lx = lg(x);
    1336      216440 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1337      216440 :       if (lx == 1) return cgetg(1, t_MAT);
    1338      137048 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1339      137048 :       return RgM_sub(x,y);
    1340             :     }
    1341     2469471 :     case t_RFRAC: case t_SER: break;
    1342             : 
    1343           0 :     default: pari_err_TYPE2("+",x,y);
    1344             :   }
    1345    60093630 :   av = avma;
    1346    60093630 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1347             : }
    1348             : 
    1349             : /********************************************************************/
    1350             : /**                                                                **/
    1351             : /**                        MULTIPLICATION                          **/
    1352             : /**                                                                **/
    1353             : /********************************************************************/
    1354             : static GEN
    1355      316148 : mul_ser_scal(GEN y, GEN x)
    1356             : {
    1357             :   long l, i;
    1358             :   GEN z;
    1359      316148 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1360      312900 :   if (ser_isexactzero(y))
    1361             :   {
    1362         525 :     z = scalarser(lg(y) == 2? Rg_get_0(x): gmul(gel(y,2), x), varn(y), 1);
    1363         525 :     setvalser(z, valser(y)); return z;
    1364             :   }
    1365      312375 :   z = cgetg_copy(y, &l); z[1] = y[1];
    1366     4217857 :   for (i = 2; i < l; i++) gel(z,i) = gmul(gel(y,i), x);
    1367      312375 :   return normalizeser(z);
    1368             : }
    1369             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1370             :  * [n/d a valid RFRAC]  */
    1371             : static GEN
    1372    10453385 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1373             : {
    1374    10453385 :   pari_sp av = avma;
    1375             :   GEN z;
    1376             : 
    1377    10453385 :   switch(typ(x))
    1378             :   {
    1379          21 :     case t_PADIC:
    1380          21 :       n = gmul(n, x);
    1381          21 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1382          21 :       return gerepileupto(av, gdiv(n,d));
    1383             : 
    1384         518 :     case t_INTMOD: case t_POLMOD:
    1385         518 :       n = gmul(n, x);
    1386         518 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1387         518 :       return gerepileupto(av, gdiv(n,d));
    1388             :   }
    1389    10452846 :   z = gred_rfrac2(x, d);
    1390    10452846 :   n = simplify_shallow(n);
    1391    10452846 :   if (typ(z) == t_RFRAC)
    1392             :   {
    1393     7936833 :     n = gmul(gel(z,1), n);
    1394     7936833 :     d = gel(z,2);
    1395     7936833 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1396           0 :       z = RgX_Rg_div(n, d);
    1397             :     else
    1398     7936833 :       z = gred_rfrac_simple(n, d);
    1399             :   }
    1400             :   else
    1401     2516013 :     z = gmul(z, n);
    1402    10452846 :   return gerepileupto(av, z);
    1403             : }
    1404             : static GEN
    1405   117004123 : mul_scal(GEN y, GEN x, long ty)
    1406             : {
    1407   117004123 :   switch(ty)
    1408             :   {
    1409   108109308 :     case t_POL:
    1410   108109308 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1411   104957864 :       return RgX_Rg_mul(y, x);
    1412      308007 :     case t_SER: return mul_ser_scal(y, x);
    1413     8586788 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1414          14 :     case t_QFB:
    1415          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1416             :   }
    1417          12 :   pari_err_TYPE2("*",x,y);
    1418             :   return NULL; /* LCOV_EXCL_LINE */
    1419             : }
    1420             : 
    1421             : static GEN
    1422      160865 : mul_gen_rfrac(GEN X, GEN Y)
    1423             : {
    1424      160865 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1425      160865 :   long vx = gvar(X), vy = varn(y2);
    1426      166626 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1427        5761 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1428             : }
    1429             : /* (x1/x2) * (y1/y2) */
    1430             : static GEN
    1431     7907128 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1432             : {
    1433             :   GEN z, X, Y;
    1434     7907128 :   pari_sp av = avma;
    1435             : 
    1436     7907128 :   X = gred_rfrac2(x1, y2);
    1437     7907128 :   Y = gred_rfrac2(y1, x2);
    1438     7907128 :   if (typ(X) == t_RFRAC)
    1439             :   {
    1440     6628046 :     if (typ(Y) == t_RFRAC) {
    1441     6562470 :       x1 = gel(X,1);
    1442     6562470 :       x2 = gel(X,2);
    1443     6562470 :       y1 = gel(Y,1);
    1444     6562470 :       y2 = gel(Y,2);
    1445     6562470 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1446             :     } else
    1447       65576 :       z = mul_gen_rfrac(Y, X);
    1448             :   }
    1449     1279082 :   else if (typ(Y) == t_RFRAC)
    1450       95289 :     z = mul_gen_rfrac(X, Y);
    1451             :   else
    1452     1183793 :     z = gmul(X, Y);
    1453     7907128 :   return gerepileupto(av, z);
    1454             : }
    1455             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1456             : static GEN
    1457      271003 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1458             : {
    1459      271003 :   pari_sp av = avma;
    1460      271003 :   GEN X = gred_rfrac2(x1, y2);
    1461      271003 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1462             :   {
    1463      264423 :     x2 = RgX_mul(gel(X,2), x2);
    1464      264423 :     x1 = gel(X,1);
    1465             :   }
    1466             :   else
    1467        6580 :     x1 = X;
    1468      271003 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1469             : }
    1470             : 
    1471             : /* Mod(y, Y) * x,  assuming x scalar */
    1472             : static GEN
    1473     3440629 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1474     3440629 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1475             : 
    1476             : /* cf mulqq */
    1477             : static GEN
    1478     5908349 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1479             : {
    1480     5908349 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1481     5908349 :   pari_sp tetpil, av = avma;
    1482     5908349 :   T[1] = x[1];
    1483     5908349 :   p2 = gmul(gel(x,2), gel(y,2));
    1484     5908349 :   p3 = gmul(gel(x,3), gel(y,3));
    1485     5908349 :   p1 = gmul(gneg_i(c),p3);
    1486             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1487     5908349 :   if (typ(b) == t_INT)
    1488             :   {
    1489     5908328 :     if (signe(b))
    1490             :     {
    1491     4337028 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1492     4337028 :       if (is_pm1(b))
    1493             :       {
    1494     4336370 :         if (signe(b) > 0) p3 = gneg(p3);
    1495             :       }
    1496             :       else
    1497         658 :         p3 = gmul(negi(b), p3);
    1498             :     }
    1499             :     else
    1500             :     {
    1501     1571300 :       p3 = gmul(gel(x,2),gel(y,3));
    1502     1571300 :       p4 = gmul(gel(x,3),gel(y,2));
    1503             :     }
    1504             :   }
    1505             :   else
    1506             :   {
    1507          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1508          21 :     p3 = gmul(gneg_i(b), p3);
    1509             :   }
    1510     5908349 :   tetpil = avma;
    1511     5908349 :   gel(T,2) = gadd(p2, p1);
    1512     5908349 :   gel(T,3) = gadd(p4, p3);
    1513     5908349 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1514     5908349 :   return normalizepol_lg(T,4);
    1515             : }
    1516             : /* Mod(x,T) * Mod(y,T) */
    1517             : static GEN
    1518     8630564 : mul_polmod_same(GEN T, GEN x, GEN y)
    1519             : {
    1520     8630564 :   GEN z = cgetg(3,t_POLMOD), a;
    1521     8630564 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1522     8630564 :   gel(z,1) = RgX_copy(T);
    1523             :   /* x * y mod T optimised */
    1524     8630564 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1525     7937335 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1526     1610188 :     a = gmul(x, y);
    1527             :   else
    1528             :   {
    1529     7020376 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1530     5903078 :       a = quad_polmod_mul(T, x, y);
    1531             :     else
    1532     1117298 :       a = RgXQ_mul(x, y, gel(z,1));
    1533             :   }
    1534     8630564 :   gel(z,2) = a; return z;
    1535             : }
    1536             : static GEN
    1537      484317 : sqr_polmod(GEN T, GEN x)
    1538             : {
    1539      484317 :   GEN a, z = cgetg(3,t_POLMOD);
    1540      484317 :   gel(z,1) = RgX_copy(T);
    1541      484317 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1542       93668 :     a = gsqr(x);
    1543             :   else
    1544             :   {
    1545      390649 :     pari_sp av = avma;
    1546      390649 :     a = RgXQ_sqr(x, gel(z,1));
    1547      390649 :     a = gerepileupto(av, a);
    1548             :   }
    1549      484317 :   gel(z,2) = a; return z;
    1550             : }
    1551             : /* Mod(x,X) * Mod(y,Y) */
    1552             : static GEN
    1553     2678214 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1554             : {
    1555     2678214 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1556     2678214 :   long vx = varn(X), vy = varn(Y);
    1557     2678214 :   GEN z = cgetg(3,t_POLMOD);
    1558             : 
    1559     2678214 :   if (vx==vy) {
    1560             :     pari_sp av;
    1561          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1562          14 :     warn_coercion(X,Y,gel(z,1));
    1563          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1564          14 :     return z;
    1565             :   }
    1566     2678200 :   if (varncmp(vx, vy) < 0)
    1567      414316 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1568             :   else
    1569     2263884 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1570     2678200 :   gel(z,2) = gmul(x, y); return z;
    1571             : }
    1572             : 
    1573             : #if 0 /* used by 3M only */
    1574             : /* set z = x+y and return 1 if x,y have the same sign
    1575             :  * set z = x-y and return 0 otherwise */
    1576             : static int
    1577             : did_add(GEN x, GEN y, GEN *z)
    1578             : {
    1579             :   long tx = typ(x), ty = typ(y);
    1580             :   if (tx == ty) switch(tx)
    1581             :   {
    1582             :     case t_INT: *z = addii(x,y); return 1;
    1583             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1584             :     case t_REAL:
    1585             :       if (signe(x) == -signe(y))
    1586             :       { *z = subrr(x,y); return 0; }
    1587             :       else
    1588             :       { *z = addrr(x,y); return 1; }
    1589             :   }
    1590             :   if (tx == t_REAL) switch(ty)
    1591             :   {
    1592             :     case t_INT:
    1593             :       if (signe(x) == -signe(y))
    1594             :       { *z = subri(x,y); return 0; }
    1595             :       else
    1596             :       { *z = addri(x,y); return 1; }
    1597             :     case t_FRAC:
    1598             :       if (signe(x) == -signe(gel(y,1)))
    1599             :       { *z = gsub(x,y); return 0; }
    1600             :       else
    1601             :       { *z = gadd(x,y); return 1; }
    1602             :   }
    1603             :   else if (ty == t_REAL) switch(tx)
    1604             :   {
    1605             :     case t_INT:
    1606             :       if (signe(x) == -signe(y))
    1607             :       { *z = subir(x,y); return 0; }
    1608             :       else
    1609             :       { *z = addir(x,y); return 1; }
    1610             :     case t_FRAC:
    1611             :       if (signe(gel(x,1)) == -signe(y))
    1612             :       { *z = gsub(x,y); return 0; }
    1613             :       else
    1614             :       { *z = gadd(x,y); return 1; }
    1615             :   }
    1616             :   *z = gadd(x,y); return 1;
    1617             : }
    1618             : #endif
    1619             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1620             : static GEN
    1621    11632523 : mulcIR(GEN x, GEN y)
    1622             : {
    1623    11632523 :   GEN z = cgetg(3,t_COMPLEX);
    1624    11632566 :   pari_sp av = avma;
    1625    11632566 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1626    11632757 :   gel(z,2) = gmul(y, gel(x,1));
    1627    11632604 :   return z;
    1628             : 
    1629             : }
    1630             : /* x,y COMPLEX */
    1631             : static GEN
    1632   273493633 : mulcc(GEN x, GEN y)
    1633             : {
    1634   273493633 :   GEN xr = gel(x,1), xi = gel(x,2);
    1635   273493633 :   GEN yr = gel(y,1), yi = gel(y,2);
    1636             :   GEN p1, p2, p3, p4, z;
    1637             :   pari_sp tetpil, av;
    1638             : 
    1639   273493633 :   if (isintzero(xr))
    1640             :   {
    1641    15606648 :     if (isintzero(yr)) {
    1642     7116434 :       av = avma;
    1643     7116434 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1644             :     }
    1645     8490083 :     return mulcIR(y, xi);
    1646             :   }
    1647   257879846 :   if (isintzero(yr)) return mulcIR(x, yi);
    1648             : 
    1649   254727706 :   z = cgetg(3,t_COMPLEX); av = avma;
    1650             : #if 0
    1651             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1652             :    * e.g. xr + xi if exponents differ */
    1653             :   if (did_add(xr, xi, &p3))
    1654             :   {
    1655             :     if (did_add(yr, yi, &p4)) {
    1656             :     /* R = xr*yr - xi*yi
    1657             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1658             :       p1 = gmul(xr,yr);
    1659             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1660             :       p3 = gmul(p3, p4);
    1661             :       p4 = gsub(p2, p1);
    1662             :     } else {
    1663             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1664             :      * I = xr*yi + xi*yr */
    1665             :       p1 = gmul(p3,p4);
    1666             :       p3 = gmul(xr,yi);
    1667             :       p4 = gmul(xi,yr);
    1668             :       p2 = gsub(p3, p4);
    1669             :     }
    1670             :   } else {
    1671             :     if (did_add(yr, yi, &p4)) {
    1672             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1673             :       * I = xr*yi +xi*yr */
    1674             :       p1 = gmul(p3,p4);
    1675             :       p3 = gmul(xr,yi);
    1676             :       p4 = gmul(xi,yr);
    1677             :       p2 = gsub(p4, p3);
    1678             :     } else {
    1679             :     /* R = xr*yr - xi*yi
    1680             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1681             :       p3 = gneg( gmul(p3, p4) );
    1682             :       p1 = gmul(xr,yr);
    1683             :       p2 = gmul(xi,yi);
    1684             :       p4 = gadd(p1, p2);
    1685             : 
    1686             :       p2 = gneg(p2);
    1687             :     }
    1688             :   }
    1689             :   tetpil = avma;
    1690             :   gel(z,1) = gadd(p1,p2);
    1691             :   gel(z,2) = gadd(p3,p4);
    1692             : #else
    1693   254719398 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1694             :   { /* 3M formula */
    1695      559286 :     p3 = addii(xr,xi);
    1696      559286 :     p4 = addii(yr,yi);
    1697      559286 :     p1 = mulii(xr,yr);
    1698      559286 :     p2 = mulii(xi,yi);
    1699      559286 :     p3 = mulii(p3,p4);
    1700      559286 :     p4 = addii(p2,p1);
    1701      559286 :     tetpil = avma;
    1702      559286 :     gel(z,1) = subii(p1,p2);
    1703      559286 :     gel(z,2) = subii(p3,p4);
    1704      559286 :     if (!signe(gel(z,2)))
    1705      113225 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1706             :   }
    1707             :   else
    1708             :   { /* naive 4M formula: avoid all loss of accuracy */
    1709   254160112 :     p1 = gmul(xr,yr);
    1710   254128223 :     p2 = gmul(xi,yi);
    1711   254131587 :     p3 = gmul(xr,yi);
    1712   254134890 :     p4 = gmul(xi,yr);
    1713   254140948 :     tetpil = avma;
    1714   254140948 :     gel(z,1) = gsub(p1,p2);
    1715   254015562 :     gel(z,2) = gadd(p3,p4);
    1716   254024170 :     if (isintzero(gel(z,2)))
    1717             :     {
    1718       50575 :       cgiv(gel(z,2));
    1719       50575 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1720             :     }
    1721             :   }
    1722             : #endif
    1723             : 
    1724   254414119 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1725             : }
    1726             : /* x,y PADIC */
    1727             : static GEN
    1728    17470327 : mulpp(GEN x, GEN y) {
    1729    17470327 :   long l = valp(x) + valp(y);
    1730             :   pari_sp av;
    1731             :   GEN z, t;
    1732    17470327 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1733    17470322 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1734    17243668 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1735             : 
    1736    16897120 :   t = (precp(x) > precp(y))? y: x;
    1737    16897120 :   z = cgetp(t); setvalp(z,l); av = avma;
    1738    16897262 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1739    16897100 :   set_avma(av); return z;
    1740             : }
    1741             : /* x,y QUAD */
    1742             : static GEN
    1743        1106 : mulqq(GEN x, GEN y) {
    1744        1106 :   GEN z = cgetg(4,t_QUAD);
    1745        1106 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1746             :   pari_sp av, tetpil;
    1747        1106 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1748             : 
    1749        1106 :   gel(z,1) = ZX_copy(P); av = avma;
    1750        1106 :   p2 = gmul(gel(x,2),gel(y,2));
    1751        1106 :   p3 = gmul(gel(x,3),gel(y,3));
    1752        1106 :   p1 = gmul(gneg_i(c),p3);
    1753             : 
    1754        1106 :   if (signe(b))
    1755         987 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1756             :   else
    1757             :   {
    1758         119 :     p3 = gmul(gel(x,2),gel(y,3));
    1759         119 :     p4 = gmul(gel(x,3),gel(y,2));
    1760             :   }
    1761        1106 :   tetpil = avma;
    1762        1106 :   gel(z,2) = gadd(p2,p1);
    1763        1106 :   gel(z,3) = gadd(p4,p3);
    1764        1106 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1765             : }
    1766             : 
    1767             : GEN
    1768    15554387 : mulcxI(GEN x)
    1769             : {
    1770             :   GEN z;
    1771    15554387 :   switch(typ(x))
    1772             :   {
    1773     2618440 :     case t_INT: case t_REAL: case t_FRAC:
    1774     2618440 :       return mkcomplex(gen_0, x);
    1775    12936341 :     case t_COMPLEX:
    1776    12936341 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1777    12922869 :       z = cgetg(3,t_COMPLEX);
    1778    12923140 :       gel(z,1) = gneg(gel(x,2));
    1779    12923790 :       gel(z,2) = gel(x,1); return z;
    1780          48 :     default:
    1781          48 :       return gmul(gen_I(), x);
    1782             :   }
    1783             : }
    1784             : GEN
    1785     3142189 : mulcxmI(GEN x)
    1786             : {
    1787             :   GEN z;
    1788     3142189 :   switch(typ(x))
    1789             :   {
    1790      116143 :     case t_INT: case t_REAL: case t_FRAC:
    1791      116143 :       return mkcomplex(gen_0, gneg(x));
    1792     2830200 :     case t_COMPLEX:
    1793     2830200 :       if (isintzero(gel(x,1))) return gel(x,2);
    1794     2792012 :       z = cgetg(3,t_COMPLEX);
    1795     2791974 :       gel(z,1) = gel(x,2);
    1796     2791974 :       gel(z,2) = gneg(gel(x,1)); return z;
    1797      195846 :     default:
    1798      195846 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1799             :   }
    1800             : }
    1801             : /* x * I^k */
    1802             : GEN
    1803     5523672 : mulcxpowIs(GEN x, long k)
    1804             : {
    1805     5523672 :   switch (k & 3)
    1806             :   {
    1807     1368989 :     case 1: return mulcxI(x);
    1808     1361844 :     case 2: return gneg(x);
    1809     1363568 :     case 3: return mulcxmI(x);
    1810             :   }
    1811     1429271 :   return x;
    1812             : }
    1813             : 
    1814             : /* caller will assume l > 2 later */
    1815             : static GEN
    1816     5705829 : init_ser(long l, long v, long e)
    1817             : {
    1818     5705829 :   GEN z = cgetg(l, t_SER);
    1819     5705829 :   z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
    1820             : }
    1821             : 
    1822             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1823             : static GEN
    1824     5693636 : fill_ser(GEN z, GEN y)
    1825             : {
    1826     5693636 :   long i, lx = lg(z), ly = lg(y); /* lx > 2 */
    1827     5693636 :   if (ly >= lx) {
    1828    25569801 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1829             :   } else {
    1830      335275 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1831      238887 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1832             :   }
    1833             :   /* dangerous case (already normalized), don't use normalizeser */
    1834     5693636 :   if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
    1835     5693538 :   return normalizeser(z);
    1836             : }
    1837             : /* assume typ(x) = t_VEC */
    1838             : static int
    1839          98 : is_ext_qfr(GEN x)
    1840          91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
    1841         189 :                     && typ(gel(x,2)) == t_REAL; }
    1842             : 
    1843             : GEN
    1844  8420862479 : gmul(GEN x, GEN y)
    1845             : {
    1846             :   long tx, ty, lx, ly, vx, vy, i, l;
    1847             :   pari_sp av, tetpil;
    1848             :   GEN z, p1;
    1849             : 
    1850  8420862479 :   if (x == y) return gsqr(x);
    1851  7544069627 :   tx = typ(x); ty = typ(y);
    1852  7544069627 :   if (tx == ty) switch(tx)
    1853             :   {
    1854  3555846243 :     case t_INT: return mulii(x,y);
    1855  1850166104 :     case t_REAL: return mulrr(x,y);
    1856     1780325 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1857     1780325 :       z = cgetg(3,t_INTMOD);
    1858     1780325 :       if (X==Y || equalii(X,Y))
    1859     1780290 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1860          35 :       gel(z,1) = gcdii(X,Y);
    1861          35 :       warn_coercion(X,Y,gel(z,1));
    1862          35 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1863          35 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1864             :     }
    1865    23293183 :     case t_FRAC:
    1866             :     {
    1867    23293183 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1868    23293183 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1869    23293183 :       z=cgetg(3,t_FRAC);
    1870    23293182 :       p1 = gcdii(x1, y2);
    1871    23293181 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1872    23293179 :       p1 = gcdii(x2, y1);
    1873    23293181 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1874    23293181 :       tetpil = avma;
    1875    23293181 :       gel(z,2) = mulii(x2,y2);
    1876    23293181 :       gel(z,1) = mulii(x1,y1);
    1877    23293184 :       fix_frac_if_int_GC(z,tetpil); return z;
    1878             :     }
    1879   264527622 :     case t_COMPLEX: return mulcc(x, y);
    1880    12480101 :     case t_PADIC: return mulpp(x, y);
    1881         882 :     case t_QUAD: return mulqq(x, y);
    1882    14546069 :     case t_FFELT: return FF_mul(x, y);
    1883    11060950 :     case t_POLMOD:
    1884    11060950 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1885     8382736 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1886     2678214 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1887    45348230 :     case t_POL:
    1888    45348230 :       vx = varn(x);
    1889    45348230 :       vy = varn(y);
    1890    45348230 :       if (vx != vy) {
    1891     4866382 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1892     2080584 :         else                     return RgX_Rg_mul(y, x);
    1893             :       }
    1894    40481848 :       return RgX_mul(x, y);
    1895             : 
    1896     4411770 :     case t_SER: {
    1897     4411770 :       vx = varn(x);
    1898     4411770 :       vy = varn(y);
    1899     4411770 :       if (vx != vy) {
    1900        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1901        3675 :         return mul_ser_scal(y, x);
    1902             :       }
    1903     4408095 :       lx = minss(lg(x), lg(y));
    1904     4408095 :       if (lx == 2) return zeroser(vx, valser(x)+valser(y));
    1905     4408081 :       av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
    1906     4408081 :       x = ser2pol_i(x, lx);
    1907     4408081 :       y = ser2pol_i(y, lx);
    1908     4408081 :       y = RgXn_mul(x, y, lx-2);
    1909     4408081 :       return gerepilecopy(av, fill_ser(z,y));
    1910             :     }
    1911          42 :     case t_VEC:
    1912          42 :       if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
    1913             :     /* fall through, handle extended t_QFB */
    1914       47804 :     case t_QFB: return qfbcomp(x,y);
    1915     6720472 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1916     4063291 :     case t_MAT: return RgM_mul(x, y);
    1917             : 
    1918        1631 :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1919        1631 :       z = cgetg_copy(x, &l);
    1920        1631 :       if (l != lg(y)) break;
    1921       17325 :       for (i=1; i<l; i++)
    1922             :       {
    1923       15694 :         long yi = y[i];
    1924       15694 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1925       15694 :         z[i] = x[yi];
    1926             :       }
    1927        1631 :       return z;
    1928             : 
    1929           0 :     default:
    1930           0 :       pari_err_TYPE2("*",x,y);
    1931             :   }
    1932             :   /* tx != ty */
    1933  1752110825 :   if (is_const_t(ty) && is_const_t(tx))  {
    1934  1624725342 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1935  1624725342 :     switch(tx) {
    1936  1496193198 :     case t_INT:
    1937             :       switch(ty)
    1938             :       {
    1939   985704941 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1940     1322566 :         case t_INTMOD:
    1941     1322566 :           z = cgetg(3, t_INTMOD);
    1942     1322567 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1943    71717442 :         case t_FRAC:
    1944    71717442 :           if (!signe(x)) return gen_0;
    1945    54641989 :           z=cgetg(3,t_FRAC);
    1946    54642206 :           p1 = gcdii(x,gel(y,2));
    1947    54641467 :           if (equali1(p1))
    1948             :           {
    1949    28908536 :             set_avma((pari_sp)z);
    1950    28908525 :             gel(z,2) = icopy(gel(y,2));
    1951    28908542 :             gel(z,1) = mulii(gel(y,1), x);
    1952             :           }
    1953             :           else
    1954             :           {
    1955    25733208 :             x = diviiexact(x,p1); tetpil = avma;
    1956    25732773 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1957    25732622 :             gel(z,1) = mulii(gel(y,1), x);
    1958    25732991 :             fix_frac_if_int_GC(z,tetpil);
    1959             :           }
    1960    54641954 :           return z;
    1961   434329832 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1962     4196282 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1963        1701 :         case t_QUAD: return mulRq(x,y);
    1964     1625610 :         case t_FFELT: return FF_Z_mul(y,x);
    1965             :       }
    1966             : 
    1967             :     case t_REAL:
    1968   120928039 :       switch(ty)
    1969             :       {
    1970    37482941 :         case t_FRAC: return mulrfrac(x, y);
    1971    83445056 :         case t_COMPLEX: return mulRc(x, y);
    1972          21 :         case t_QUAD: return mulqf(y, x, realprec(x));
    1973          21 :         default: pari_err_TYPE2("*",x,y);
    1974             :       }
    1975             : 
    1976        8022 :     case t_INTMOD:
    1977             :       switch(ty)
    1978             :       {
    1979        7133 :         case t_FRAC: { GEN X = gel(x,1);
    1980        7133 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1981        7133 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1982             :         }
    1983          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1984         427 :         case t_PADIC: { GEN X = gel(x,1);
    1985         427 :           z = cgetg(3, t_INTMOD);
    1986         427 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1987             :         }
    1988          63 :         case t_QUAD: return mulRq(x, y);
    1989         350 :         case t_FFELT:
    1990         350 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1991           0 :             pari_err_OP("*",x,y);
    1992         350 :           return FF_Z_mul(y,gel(x,2));
    1993             :       }
    1994             : 
    1995             :     case t_FRAC:
    1996             :       switch(ty)
    1997             :       {
    1998     5207014 :         case t_COMPLEX: return mulRc(x, y);
    1999     2581571 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    2000         343 :         case t_QUAD:  return mulRq(x, y);
    2001        2051 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    2002             :       }
    2003             : 
    2004             :     case t_FFELT:
    2005          35 :       pari_err_TYPE2("*",x,y);
    2006             : 
    2007          21 :     case t_COMPLEX:
    2008             :       switch(ty)
    2009             :       {
    2010          14 :         case t_PADIC:
    2011          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    2012           7 :         case t_QUAD:
    2013           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    2014           7 :           return mulqf(y, x, lx);
    2015             :       }
    2016             : 
    2017             :     case t_PADIC: /* ty == t_QUAD */
    2018          28 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    2019             :     }
    2020             :   }
    2021             : 
    2022   130319878 :   if (is_matvec_t(ty))
    2023             :   {
    2024     8072398 :     if (!is_matvec_t(tx))
    2025             :     {
    2026     7123583 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    2027     7123586 :       z = cgetg_copy(y, &ly);
    2028   627575133 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2029     7123434 :       return z;
    2030             :     }
    2031      948816 :     switch(tx)
    2032             :     {
    2033      121973 :       case t_VEC:
    2034             :         switch(ty) {
    2035       15540 :           case t_COL: return RgV_RgC_mul(x,y);
    2036      106433 :           case t_MAT: return RgV_RgM_mul(x,y);
    2037             :         }
    2038           0 :         break;
    2039        1687 :       case t_COL:
    2040             :         switch(ty) {
    2041        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    2042           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    2043             :         }
    2044           0 :         break;
    2045      825166 :       case t_MAT:
    2046             :         switch(ty) {
    2047           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    2048      825166 :           case t_COL: return RgM_RgC_mul(x,y);
    2049             :         }
    2050             :     }
    2051             :   }
    2052   125348020 :   if (is_matvec_t(tx))
    2053             :   {
    2054     2939367 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2055     2939504 :     z = cgetg_copy(x, &lx);
    2056    10512858 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2057     2939134 :     return z;
    2058             :   }
    2059   122409083 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2060             :   /* tx < ty, !ismatvec(x and y) */
    2061             : 
    2062   122409083 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2063     3393246 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2064   119015837 :   if (is_scalar_t(tx)) {
    2065   115724665 :     if (tx == t_POLMOD) {
    2066     3142134 :       vx = varn(gel(x,1));
    2067     3142134 :       vy = gvar(y);
    2068     3142134 :       if (vx != vy) {
    2069     2894866 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2070       47383 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2071             :       }
    2072             :       /* error if ty == t_SER */
    2073      247268 :       av = avma; y = gmod(y, gel(x,1));
    2074      247261 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2075             :     }
    2076   112582531 :     return mul_scal(y, x, ty);
    2077             :   }
    2078             : 
    2079             :   /* x and y are not scalars, nor matvec */
    2080     3291057 :   vx = gvar(x);
    2081     3291090 :   vy = gvar(y);
    2082     3291090 :   if (vx != vy) /* x or y is treated as a scalar */
    2083     2780257 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2084     2780257 :                                 : mul_scal(y, x, ty);
    2085             :   /* vx = vy */
    2086     1872181 :   switch(tx)
    2087             :   {
    2088     1872146 :     case t_POL:
    2089             :       switch (ty)
    2090             :       {
    2091        6741 :         case t_SER:
    2092             :         {
    2093             :           long v;
    2094        6741 :           av = avma; v = RgX_valrem(x, &x);
    2095        6741 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2096        6727 :           v += valser(y); ly = lg(y);
    2097        6727 :           if (ly == 2) { set_avma(av); return zeroser(vx, v); }
    2098        6496 :           if (degpol(x))
    2099             :           {
    2100        2030 :             z = init_ser(ly, vx, v);
    2101        2030 :             x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
    2102        2030 :             return gerepilecopy(av, fill_ser(z, x));
    2103             :           }
    2104             :           /* take advantage of x = c*t^v */
    2105        4466 :           set_avma(av); y = mul_ser_scal(y, gel(x,2));
    2106        4466 :           setvalser(y, v); return y;
    2107             :         }
    2108             : 
    2109     1865405 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2110             :       }
    2111           0 :       break;
    2112             : 
    2113          35 :     case t_SER:
    2114             :       switch (ty)
    2115             :       {
    2116          35 :         case t_RFRAC:
    2117          35 :           av = avma;
    2118          35 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2119             :       }
    2120           0 :       break;
    2121             :   }
    2122           0 :   pari_err_TYPE2("*",x,y);
    2123             :   return NULL; /* LCOV_EXCL_LINE */
    2124             : }
    2125             : 
    2126             : /* return a nonnormalized result */
    2127             : GEN
    2128      112381 : sqr_ser_part(GEN x, long l1, long l2)
    2129             : {
    2130             :   long i, j, l;
    2131             :   pari_sp av;
    2132             :   GEN Z, z, p1, p2;
    2133             :   long mi;
    2134      112381 :   if (l2 < l1) return zeroser(varn(x), 2*valser(x));
    2135      112367 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2136      112367 :   Z = cgetg(l2-l1+3, t_SER);
    2137      112367 :   Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
    2138      112367 :   z = Z + 2-l1;
    2139      112367 :   x += 2; mi = 0;
    2140      423584 :   for (i=0; i<l1; i++)
    2141             :   {
    2142      311217 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2143             :   }
    2144             : 
    2145      718560 :   for (i=l1; i<=l2; i++)
    2146             :   {
    2147      606193 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2148      606193 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2149     1253138 :     for (j=i-mi; j<=minss(l,mi); j++)
    2150      646945 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2151      606193 :     p1 = gshift(p1,1);
    2152      606193 :     if ((i&1) == 0 && p2[i>>1])
    2153       92969 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2154      606193 :     gel(z,i) = gerepileupto(av,p1);
    2155             :   }
    2156      112367 :   return Z;
    2157             : }
    2158             : 
    2159             : GEN
    2160  1122445455 : gsqr(GEN x)
    2161             : {
    2162             :   long i, lx;
    2163             :   pari_sp av, tetpil;
    2164             :   GEN z, p1, p2, p3, p4;
    2165             : 
    2166  1122445455 :   switch(typ(x))
    2167             :   {
    2168   907900687 :     case t_INT: return sqri(x);
    2169   197556947 :     case t_REAL: return sqrr(x);
    2170      142395 :     case t_INTMOD: { GEN X = gel(x,1);
    2171      142395 :       z = cgetg(3,t_INTMOD);
    2172      142396 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2173      142398 :       gel(z,1) = icopy(X); return z;
    2174             :     }
    2175     4007902 :     case t_FRAC: return sqrfrac(x);
    2176             : 
    2177     8104168 :     case t_COMPLEX:
    2178     8104168 :       if (isintzero(gel(x,1))) {
    2179      237416 :         av = avma;
    2180      237416 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2181             :       }
    2182     7866750 :       z = cgetg(3,t_COMPLEX); av = avma;
    2183     7866705 :       p1 = gadd(gel(x,1),gel(x,2));
    2184     7866541 :       p2 = gsub(gel(x,1), gel(x,2));
    2185     7866543 :       p3 = gmul(gel(x,1),gel(x,2));
    2186     7866662 :       tetpil = avma;
    2187     7866662 :       gel(z,1) = gmul(p1,p2);
    2188     7866708 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2189             : 
    2190        4214 :     case t_PADIC:
    2191        4214 :       z = cgetg(5,t_PADIC);
    2192        4214 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2193        4214 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2194        4214 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2195        4214 :       gel(z,2) = icopy(gel(x,2));
    2196        4214 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2197        4214 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2198        4214 :       return z;
    2199             : 
    2200          70 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2201          70 :       p1 = gel(x,1);
    2202          70 :       gel(z,1) = ZX_copy(p1); av = avma;
    2203          70 :       p2 = gsqr(gel(x,2));
    2204          70 :       p3 = gsqr(gel(x,3));
    2205          70 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2206             : 
    2207          70 :       if (gequal0(gel(p1,3)))
    2208             :       {
    2209           7 :         tetpil = avma;
    2210           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2211           7 :         av = avma;
    2212           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2213           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2214             :       }
    2215             : 
    2216          63 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2217          63 :       tetpil = avma;
    2218          63 :       gel(z,2) = gadd(p2,p4);
    2219          63 :       gel(z,3) = gadd(p1,p3);
    2220          63 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2221             : 
    2222      484317 :     case t_POLMOD:
    2223      484317 :       return sqr_polmod(gel(x,1), gel(x,2));
    2224             : 
    2225     2957504 :     case t_FFELT: return FF_sqr(x);
    2226             : 
    2227     1186251 :     case t_POL: return RgX_sqr(x);
    2228             : 
    2229       35014 :     case t_SER:
    2230       35014 :       lx = lg(x);
    2231       35014 :       if (ser_isexactzero(x)) {
    2232          21 :         GEN z = gcopy(x);
    2233          21 :         setvalser(z, 2*valser(x));
    2234          21 :         return z;
    2235             :       }
    2236       34993 :       if (lx < 40)
    2237       34734 :         return normalizeser( sqr_ser_part(x, 0, lx-3) );
    2238             :       else
    2239             :       {
    2240         259 :         pari_sp av = avma;
    2241         259 :         GEN z = init_ser(lx, varn(x), 2*valser(x));
    2242         259 :         x = ser2pol_i(x, lx);
    2243         259 :         x = RgXn_sqr(x, lx-2);
    2244         259 :         return gerepilecopy(av, fill_ser(z,x));
    2245             :       }
    2246             : 
    2247          49 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2248          49 :       gel(z,1) = gsqr(gel(x,1));
    2249          49 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2250             : 
    2251        1099 :     case t_MAT: return RgM_sqr(x);
    2252          28 :     case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
    2253             :     /* fall through handle extended t_QFB */
    2254       81919 :     case t_QFB: return qfbsqr(x);
    2255         658 :     case t_VECSMALL:
    2256         658 :       z = cgetg_copy(x, &lx);
    2257       16289 :       for (i=1; i<lx; i++)
    2258             :       {
    2259       15631 :         long xi = x[i];
    2260       15631 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2261       15631 :         z[i] = x[xi];
    2262             :       }
    2263         658 :       return z;
    2264             :   }
    2265           7 :   pari_err_TYPE2("*",x,x);
    2266             :   return NULL; /* LCOV_EXCL_LINE */
    2267             : }
    2268             : 
    2269             : /********************************************************************/
    2270             : /**                                                                **/
    2271             : /**                           DIVISION                             **/
    2272             : /**                                                                **/
    2273             : /********************************************************************/
    2274             : static GEN
    2275       13169 : div_rfrac_scal(GEN x, GEN y)
    2276             : {
    2277       13169 :   pari_sp av = avma;
    2278       13169 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2279       13169 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2280             : }
    2281             : static GEN
    2282       37466 : div_scal_rfrac(GEN x, GEN y)
    2283             : {
    2284       37466 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2285       37466 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2286             :   {
    2287          14 :     if (degpol(y1))
    2288             :     {
    2289          14 :       pari_sp av = avma;
    2290          14 :       GEN _1 = Rg_get_1(x);
    2291          14 :       if (transtype(_1)) y1 = gmul(y1, _1);
    2292          14 :       return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2293             :     }
    2294           0 :     y1 = gel(y1,2);
    2295             :   }
    2296       37452 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2297             : }
    2298             : static GEN
    2299     1186656 : div_rfrac(GEN x, GEN y)
    2300     1186656 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2301             : 
    2302             : /* x != 0 */
    2303             : static GEN
    2304     1338593 : div_ser_scal(GEN y, GEN x)
    2305             : {
    2306             :   long i, l;
    2307             :   GEN z;
    2308     1338593 :   if (ser_isexactzero(y))
    2309             :   {
    2310          28 :     z = scalarser(lg(y) == 2? Rg_get_0(x): gdiv(gel(y,2), x), varn(y), 1);
    2311          28 :     setvalser(z, valser(y)); return z;
    2312             :   }
    2313     1338565 :   z = cgetg_copy(y, &l); z[1] = y[1];
    2314     6202113 :   for (i = 2; i < l; i++) gel(z,i) = gdiv(gel(y,i), x);
    2315     1338565 :   return normalizeser(z);
    2316             : }
    2317             : GEN
    2318        7560 : ser_normalize(GEN x)
    2319             : {
    2320        7560 :   long i, lx = lg(x);
    2321             :   GEN c, z;
    2322        7560 :   if (lx == 2) return x;
    2323        7560 :   c = gel(x,2); if (gequal1(c)) return x;
    2324        7483 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2325      107870 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2326        7483 :   return z;
    2327             : }
    2328             : 
    2329             : /* y != 0 */
    2330             : static GEN
    2331     5745964 : div_T_scal(GEN x, GEN y, long tx) {
    2332     5745964 :   switch(tx)
    2333             :   {
    2334     4396282 :     case t_POL: return RgX_Rg_div(x, y);
    2335     1338586 :     case t_SER: return div_ser_scal(x, y);
    2336       11097 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2337             :   }
    2338           0 :   pari_err_TYPE2("/",x,y);
    2339             :   return NULL; /* LCOV_EXCL_LINE */
    2340             : }
    2341             : 
    2342             : static GEN
    2343     9256385 : div_scal_pol(GEN x, GEN y) {
    2344     9256385 :   long ly = lg(y);
    2345             :   pari_sp av;
    2346             :   GEN _1;
    2347     9256385 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2348     9178817 :   if (isrationalzero(x)) return zeropol(varn(y));
    2349     7129024 :   av = avma;
    2350     7129024 :   _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
    2351     7129024 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2352             : }
    2353             : static GEN
    2354       18550 : div_scal_ser(GEN x, GEN y)
    2355             : {
    2356       18550 :   pari_sp av = avma;
    2357       18550 :   GEN _1 = Rg_get_1(x);
    2358       18550 :   if (transtype(_1)) y = gmul(y, _1);
    2359       18550 :   return gerepileupto(av, gmul(x, ser_inv(y)));
    2360             : }
    2361             : static GEN
    2362     9265697 : div_scal_T(GEN x, GEN y, long ty) {
    2363     9265697 :   switch(ty)
    2364             :   {
    2365     9210444 :     case t_POL: return div_scal_pol(x, y);
    2366       18543 :     case t_SER: return div_scal_ser(x, y);
    2367       36710 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2368             :   }
    2369           0 :   pari_err_TYPE2("/",x,y);
    2370             :   return NULL; /* LCOV_EXCL_LINE */
    2371             : }
    2372             : 
    2373             : /* assume tx = ty = t_SER, same variable vx */
    2374             : static GEN
    2375      771146 : div_ser(GEN x, GEN y, long vx)
    2376             : {
    2377      771146 :   long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
    2378      771146 :   GEN y0 = y, z;
    2379      771146 :   pari_sp av = avma;
    2380             : 
    2381      771146 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2382      771139 :   if (ser_isexactzero(x))
    2383             :   {
    2384       59892 :     if (lx == 2) return zeroser(vx, v);
    2385         147 :     z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
    2386         147 :     setvalser(z, v); return z;
    2387             :   }
    2388      711247 :   if (lx < ly) ly = lx;
    2389      711247 :   y = ser2pol_i_normalize(y, ly, &e);
    2390      711247 :   if (e)
    2391             :   {
    2392           7 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2393           7 :     ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
    2394             :   }
    2395      711247 :   z = init_ser(ly, vx, v);
    2396      711247 :   if (ly == 3)
    2397             :   {
    2398       12193 :     gel(z,2) = gdiv(gel(x,2), gel(y,2));
    2399       12193 :     if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
    2400       12193 :     return gerepileupto(av, z);
    2401             :   }
    2402      699054 :   x = ser2pol_i(x, ly);
    2403      699054 :   y = RgXn_div_i(x, y, ly-2);
    2404      699054 :   return gerepilecopy(av, fill_ser(z,y));
    2405             : }
    2406             : /* x,y compatible PADIC */
    2407             : static GEN
    2408    13211528 : divpp(GEN x, GEN y) {
    2409             :   pari_sp av;
    2410             :   long a, b;
    2411             :   GEN z, M;
    2412             : 
    2413    13211528 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2414    13211530 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2415    13211201 :   a = precp(x);
    2416    13211201 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2417    13211201 :   z = cgetg(5, t_PADIC);
    2418    13210736 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2419    13210648 :   gel(z,2) = icopy(gel(x,2));
    2420    13210417 :   gel(z,3) = icopy(M); av = avma;
    2421    13210137 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2422    13210905 :   return z;
    2423             : }
    2424             : static GEN
    2425       36211 : div_polmod_same(GEN T, GEN x, GEN y)
    2426             : {
    2427       36211 :   long v = varn(T);
    2428       36211 :   GEN a, z = cgetg(3, t_POLMOD);
    2429       36211 :   gel(z,1) = RgX_copy(T);
    2430       36211 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2431       10584 :     a = gdiv(x, y);
    2432       25627 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2433         119 :   {
    2434         119 :     pari_sp av = avma;
    2435         119 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2436             :   }
    2437       25508 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2438        5271 :   {
    2439        5271 :     pari_sp av = avma;
    2440        5271 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2441        5271 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2442        5271 :     a = gerepileupto(av, a);
    2443             :   }
    2444             :   else
    2445             :   {
    2446       20237 :     pari_sp av = avma;
    2447       20237 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2448       20237 :     a = gerepileupto(av, a);
    2449             :   }
    2450       36211 :   gel(z,2) = a; return z;
    2451             : }
    2452             : GEN
    2453   425638075 : gdiv(GEN x, GEN y)
    2454             : {
    2455   425638075 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2456             :   pari_sp av, tetpil;
    2457             :   GEN z, p1, p2;
    2458             : 
    2459   425638075 :   if (tx == ty) switch(tx)
    2460             :   {
    2461    85442926 :     case t_INT:
    2462    85442926 :       return Qdivii(x,y);
    2463             : 
    2464    91903229 :     case t_REAL: return divrr(x,y);
    2465       25585 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2466       25585 :       z = cgetg(3,t_INTMOD);
    2467       25585 :       if (X==Y || equalii(X,Y))
    2468       25571 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2469          14 :       gel(z,1) = gcdii(X,Y);
    2470          14 :       warn_coercion(X,Y,gel(z,1));
    2471          14 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2472          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2473             :     }
    2474     4025767 :     case t_FRAC: {
    2475     4025767 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2476     4025767 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2477     4025767 :       z = cgetg(3, t_FRAC);
    2478     4025767 :       p1 = gcdii(x1, y1);
    2479     4025766 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2480     4025766 :       p1 = gcdii(x2, y2);
    2481     4025766 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2482     4025765 :       tetpil = avma;
    2483     4025765 :       gel(z,2) = mulii(x2,y1);
    2484     4025768 :       gel(z,1) = mulii(x1,y2);
    2485     4025766 :       normalize_frac(z);
    2486     4025766 :       fix_frac_if_int_GC(z,tetpil);
    2487     4025768 :       return z;
    2488             :     }
    2489     9167155 :     case t_COMPLEX:
    2490     9167155 :       if (isintzero(gel(y,1)))
    2491             :       {
    2492      197967 :         y = gel(y,2);
    2493      197967 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2494      121079 :         z = cgetg(3,t_COMPLEX);
    2495      121079 :         gel(z,1) = gdiv(gel(x,2), y);
    2496      121079 :         av = avma;
    2497      121079 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2498      121079 :         return z;
    2499             :       }
    2500     8969183 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2501     8969206 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2502             : 
    2503      587154 :     case t_PADIC:
    2504      587154 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2505      587154 :       return divpp(x, y);
    2506             : 
    2507         238 :     case t_QUAD:
    2508         238 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2509         224 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2510         224 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2511             : 
    2512      243610 :     case t_FFELT: return FF_div(x,y);
    2513             : 
    2514       40586 :     case t_POLMOD:
    2515       40586 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2516       36211 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2517             :       else {
    2518        4375 :         av = avma;
    2519        4375 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2520             :       }
    2521       40586 :       return z;
    2522             : 
    2523    18452376 :     case t_POL:
    2524    18452376 :       vx = varn(x);
    2525    18452376 :       vy = varn(y);
    2526    18452376 :       if (vx != vy) {
    2527      101822 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2528       45941 :                             else return div_scal_pol(x, y);
    2529             :       }
    2530    18350554 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2531    18350554 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2532    18228402 :       av = avma;
    2533    18228402 :       return gerepileupto(av, gred_rfrac2(x,y));
    2534             : 
    2535      771167 :     case t_SER:
    2536      771167 :       vx = varn(x);
    2537      771167 :       vy = varn(y);
    2538      771167 :       if (vx != vy) {
    2539          21 :         if (varncmp(vx, vy) < 0)
    2540             :         {
    2541          14 :           if (!signe(y)) pari_err_INV("gdiv",y);
    2542           7 :           return div_ser_scal(x, y);
    2543             :         }
    2544           7 :         return div_scal_ser(x, y);
    2545             :       }
    2546      771146 :       return div_ser(x, y, vx);
    2547     1189344 :     case t_RFRAC:
    2548     1189344 :       vx = varn(gel(x,2));
    2549     1189344 :       vy = varn(gel(y,2));
    2550     1189344 :       if (vx != vy) {
    2551        2688 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2552         756 :                             else return div_scal_rfrac(x, y);
    2553             :       }
    2554     1186656 :       return div_rfrac(x,y);
    2555             : 
    2556          21 :     case t_VEC: /* handle extended t_QFB */
    2557          21 :     case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
    2558             : 
    2559          28 :     case t_MAT:
    2560          28 :       p1 = RgM_div(x,y);
    2561          28 :       if (!p1) pari_err_INV("gdiv",y);
    2562          21 :       return p1;
    2563             : 
    2564           0 :     default: pari_err_TYPE2("/",x,y);
    2565             :   }
    2566             : 
    2567   213793306 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2568             :   {
    2569     3553095 :     long s = signe(x);
    2570     3553095 :     if (!s) {
    2571      516087 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2572      516087 :       switch (ty)
    2573             :       {
    2574      512923 :       default: return gen_0;
    2575          28 :       case t_INTMOD:
    2576          28 :         z = cgetg(3,t_INTMOD);
    2577          28 :         gel(z,1) = icopy(gel(y,1));
    2578          28 :         gel(z,2) = gen_0; return z;
    2579        3136 :       case t_FFELT: return FF_zero(y);
    2580             :       }
    2581             :     }
    2582     3037008 :     if (is_pm1(x)) {
    2583     1086408 :       if (s > 0) return ginv(y);
    2584      193211 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2585             :     }
    2586     1950599 :     switch(ty)
    2587             :     {
    2588      639584 :       case t_REAL: return divir(x,y);
    2589          21 :       case t_INTMOD:
    2590          21 :         z = cgetg(3, t_INTMOD);
    2591          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2592      794481 :       case t_FRAC:
    2593      794481 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2594      794481 :         if (equali1(p1))
    2595             :         {
    2596      273702 :           set_avma((pari_sp)z);
    2597      273702 :           gel(z,2) = icopy(gel(y,1));
    2598      273702 :           gel(z,1) = mulii(gel(y,2), x);
    2599      273702 :           normalize_frac(z);
    2600      273702 :           fix_frac_if_int(z);
    2601             :         }
    2602             :         else
    2603             :         {
    2604      520779 :           x = diviiexact(x,p1); tetpil = avma;
    2605      520779 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2606      520779 :           gel(z,1) = mulii(gel(y,2), x);
    2607      520779 :           normalize_frac(z);
    2608      520779 :           fix_frac_if_int_GC(z,tetpil);
    2609             :         }
    2610      794481 :         return z;
    2611             : 
    2612         469 :       case t_FFELT: return Z_FF_div(x,y);
    2613      514317 :       case t_COMPLEX: return divRc(x,y);
    2614        1505 :       case t_PADIC: return divTp(x, y);
    2615         231 :       case t_QUAD:
    2616         231 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2617         231 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2618             :     }
    2619             :   }
    2620   210240201 :   if (gequal0(y))
    2621             :   {
    2622          49 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2623          28 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2624             :   }
    2625             : 
    2626   210248492 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2627             :   {
    2628   171884684 :     case t_REAL:
    2629   171884684 :       switch(ty)
    2630             :       {
    2631   169527654 :         case t_INT: return divri(x,y);
    2632      522503 :         case t_FRAC:
    2633      522503 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2634      522503 :           return gerepileuptoleaf(av, z);
    2635     1834571 :         case t_COMPLEX: return divRc(x, y);
    2636          42 :         case t_QUAD: return divfq(x, y, realprec(x));
    2637          18 :         default: pari_err_TYPE2("/",x,y);
    2638             :       }
    2639             : 
    2640         280 :     case t_INTMOD:
    2641             :       switch(ty)
    2642             :       {
    2643         196 :         case t_INT:
    2644         196 :           z = cgetg(3, t_INTMOD);
    2645         196 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2646          28 :         case t_FRAC: { GEN X = gel(x,1);
    2647          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2648          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2649             :         }
    2650           7 :         case t_FFELT:
    2651           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2652           0 :             pari_err_OP("/",x,y);
    2653           7 :           return Z_FF_div(gel(x,2),y);
    2654             : 
    2655           0 :         case t_COMPLEX:
    2656           0 :           av = avma;
    2657           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2658             : 
    2659          14 :         case t_QUAD:
    2660          14 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2661          14 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2662             : 
    2663           7 :         case t_PADIC: { GEN X = gel(x,1);
    2664           7 :           z = cgetg(3, t_INTMOD);
    2665           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2666             :         }
    2667          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2668             :       }
    2669             : 
    2670             :     case t_FRAC:
    2671             :       switch(ty)
    2672             :       {
    2673     2092246 :         case t_INT: z = cgetg(3, t_FRAC);
    2674     2092246 :         p1 = gcdii(y,gel(x,1));
    2675     2092246 :         if (equali1(p1))
    2676             :         {
    2677      823361 :           set_avma((pari_sp)z); tetpil = 0;
    2678      823361 :           gel(z,1) = icopy(gel(x,1));
    2679             :         }
    2680             :         else
    2681             :         {
    2682     1268885 :           y = diviiexact(y,p1); tetpil = avma;
    2683     1268885 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2684             :         }
    2685     2092246 :         gel(z,2) = mulii(gel(x,2),y);
    2686     2092246 :         normalize_frac(z);
    2687     2092246 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2688     2092246 :         return z;
    2689             : 
    2690       59502 :         case t_REAL:
    2691       59502 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2692       59502 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2693             : 
    2694           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2695           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2696           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2697             :         }
    2698             : 
    2699          28 :         case t_FFELT: av=avma;
    2700          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2701             : 
    2702         847 :         case t_COMPLEX: return divRc(x, y);
    2703             : 
    2704        2141 :         case t_PADIC:
    2705        2141 :           if (!signe(gel(x,1))) return gen_0;
    2706        2141 :           return divTp(x, y);
    2707             : 
    2708          42 :         case t_QUAD:
    2709          42 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2710          42 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2711             :       }
    2712             : 
    2713             :     case t_FFELT:
    2714         133 :       switch (ty)
    2715             :       {
    2716          49 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2717          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2718           7 :         case t_INTMOD:
    2719           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2720           0 :             pari_err_OP("/",x,y);
    2721           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2722          49 :         default:
    2723          49 :         pari_err_TYPE2("/",x,y);
    2724             :       }
    2725           0 :       break;
    2726             : 
    2727    13464246 :     case t_COMPLEX:
    2728             :       switch(ty)
    2729             :       {
    2730    13464244 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2731           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2732           0 :         case t_PADIC:
    2733           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2734           0 :         case t_QUAD:
    2735           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2736           0 :           return divfq(x, y, lx);
    2737             :       }
    2738             : 
    2739             :     case t_PADIC:
    2740             :       switch(ty)
    2741             :       {
    2742     1204070 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2743     1203979 :           return signe(gel(x,4))? divpT(x, y)
    2744     2408049 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2745             :         }
    2746           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2747           7 :           z = cgetg(3, t_INTMOD);
    2748           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2749             :         }
    2750          14 :         case t_COMPLEX: case t_QUAD:
    2751          14 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2752          14 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2753             : 
    2754          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2755             :       }
    2756             : 
    2757             :     case t_QUAD:
    2758             :       switch (ty)
    2759             :       {
    2760        1190 :         case t_INT: case t_INTMOD: case t_FRAC:
    2761        1190 :           z = cgetg(4,t_QUAD);
    2762        1190 :           gel(z,1) = ZX_copy(gel(x,1));
    2763        1190 :           gel(z,2) = gdiv(gel(x,2), y);
    2764        1190 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2765          63 :         case t_REAL: return divqf(x, y, realprec(y));
    2766          14 :         case t_PADIC: return divTp(x, y);
    2767           0 :         case t_COMPLEX:
    2768           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2769           0 :           return divqf(x, y, ly);
    2770             :       }
    2771             :   }
    2772    21538830 :   switch(ty) {
    2773      580865 :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2774      580865 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2775          42 :     case t_MAT:
    2776          42 :       av = avma; p1 = RgM_inv(y);
    2777          35 :       if (!p1) pari_err_INV("gdiv",y);
    2778          28 :       return gerepileupto(av, gmul(x, p1));
    2779           0 :     case t_VEC: case t_COL:
    2780             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2781           0 :       pari_err_TYPE2("/",x,y);
    2782             :   }
    2783    20957923 :   switch(tx) {
    2784     3837010 :     case t_VEC: case t_COL: case t_MAT:
    2785     3837010 :       z = cgetg_copy(x, &lx);
    2786    12936225 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2787     3836943 :       return z;
    2788           0 :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2789           0 :       pari_err_TYPE2("/",x,y);
    2790             :   }
    2791             : 
    2792    17120913 :   vy = gvar(y);
    2793    17121684 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2794      208835 :     vx = varn(X);
    2795      208835 :     if (vx != vy) {
    2796      208261 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2797      207932 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2798             :     }
    2799             :     /* y is POL, SER or RFRAC */
    2800         574 :     av = avma;
    2801         574 :     switch(ty)
    2802             :     {
    2803           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2804         574 :       default: y = ginvmod(gmod(y,X), X);
    2805             :     }
    2806         567 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2807             :   }
    2808             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2809             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2810             :    * variable NO_VARIABLE, then the operation is incorrect */
    2811    16912849 :   vx = gvar(x);
    2812    16912842 :   if (vx != vy) { /* includes cases where one is scalar */
    2813    15011330 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2814     9265367 :                         else return div_scal_T(x, y, ty);
    2815             :   }
    2816     1901512 :   switch(tx)
    2817             :   {
    2818     1046164 :     case t_POL:
    2819             :       switch(ty)
    2820             :       {
    2821          28 :         case t_SER:
    2822             :         {
    2823          28 :           GEN y0 = y;
    2824             :           long v;
    2825          28 :           av = avma; v = RgX_valrem(x, &x);
    2826          28 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2827          14 :           v -= valser(y); ly = lg(y); /* > 2 */
    2828          14 :           y = ser2pol_i_normalize(y, ly, &i);
    2829          14 :           if (i)
    2830             :           {
    2831           7 :             pari_warn(warner,"normalizing a series with 0 leading term");
    2832           7 :             ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
    2833             :           }
    2834           7 :           z = init_ser(ly, vx, v);
    2835           7 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
    2836             :         }
    2837             : 
    2838     1046136 :         case t_RFRAC:
    2839             :         {
    2840     1046136 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2841     1046136 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2842        1171 :             return mul_rfrac_scal(y2, y1, x);
    2843     1044965 :           av = avma;
    2844     1044965 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2845             :         }
    2846             :       }
    2847           0 :       break;
    2848             : 
    2849      584324 :     case t_SER:
    2850             :       switch(ty)
    2851             :       {
    2852      584310 :         case t_POL:
    2853             :         {
    2854      584310 :           long v = valser(x);
    2855      584310 :           lx = lg(x);
    2856      584310 :           if (lx == 2) return zeroser(vx, v - RgX_val(y));
    2857      584205 :           av = avma;
    2858      584205 :           x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
    2859      584205 :           z = init_ser(lx, vx, v);
    2860      584205 :           if (!signe(x)) setsigne(z,0);
    2861      584205 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
    2862             :         }
    2863          14 :         case t_RFRAC:
    2864          14 :           av = avma;
    2865          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2866             :       }
    2867           0 :       break;
    2868             : 
    2869      271003 :     case t_RFRAC:
    2870             :       switch(ty)
    2871             :       {
    2872      271003 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2873           0 :         case t_SER:
    2874           0 :           av = avma;
    2875           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
    2876             :       }
    2877           0 :       break;
    2878             :   }
    2879          21 :   pari_err_TYPE2("/",x,y);
    2880             :   return NULL; /* LCOV_EXCL_LINE */
    2881             : }
    2882             : 
    2883             : /********************************************************************/
    2884             : /**                                                                **/
    2885             : /**                     SIMPLE MULTIPLICATION                      **/
    2886             : /**                                                                **/
    2887             : /********************************************************************/
    2888             : GEN
    2889   227362307 : gmulsg(long s, GEN y)
    2890             : {
    2891             :   long ly, i;
    2892             :   pari_sp av;
    2893             :   GEN z;
    2894             : 
    2895   227362307 :   switch(typ(y))
    2896             :   {
    2897   135270605 :     case t_INT:  return mulsi(s,y);
    2898    68924966 :     case t_REAL: return s? mulsr(s,y): gen_0; /* gmul semantic */
    2899      174030 :     case t_INTMOD: { GEN p = gel(y,1);
    2900      174030 :       z = cgetg(3,t_INTMOD);
    2901      174031 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2902      174029 :       gel(z,1) = icopy(p); return z;
    2903             :     }
    2904      548214 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2905     6204523 :     case t_FRAC:
    2906     6204523 :       if (!s) return gen_0;
    2907     6189991 :       z = cgetg(3,t_FRAC);
    2908     6189991 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2909     6189991 :       if (i == 1)
    2910             :       {
    2911     2267501 :         gel(z,2) = icopy(gel(y,2));
    2912     2267501 :         gel(z,1) = mulis(gel(y,1), s);
    2913             :       }
    2914             :       else
    2915             :       {
    2916     3922490 :         gel(z,2) = diviuexact(gel(y,2), (ulong)i);
    2917     3922490 :         gel(z,1) = mulis(gel(y,1), s/i);
    2918     3922490 :         fix_frac_if_int(z);
    2919             :       }
    2920     6189991 :       return z;
    2921             : 
    2922    14060317 :     case t_COMPLEX:
    2923    14060317 :       if (!s) return gen_0;
    2924    12934499 :       z = cgetg(3, t_COMPLEX);
    2925    12934445 :       gel(z,1) = gmulsg(s,gel(y,1));
    2926    12932992 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2927             : 
    2928        1435 :     case t_PADIC:
    2929        1435 :       if (!s) return gen_0;
    2930        1435 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2931             : 
    2932           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2933           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2934           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2935           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2936             : 
    2937      205359 :     case t_POLMOD:
    2938      205359 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2939             : 
    2940      748697 :     case t_POL:
    2941      748697 :       if (!signe(y)) return RgX_copy(y);
    2942      732401 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2943      729188 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2944     3090773 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2945      729187 :       return normalizepol_lg(z, ly);
    2946             : 
    2947         182 :     case t_SER:
    2948         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2949         182 :       if (!s) return Rg_get_0(y);
    2950         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2951        3864 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2952         182 :       return normalizeser(z);
    2953             : 
    2954           0 :     case t_RFRAC:
    2955           0 :       if (!s) return zeropol(varn(gel(y,2)));
    2956           0 :       if (s == 1) return gcopy(y);
    2957           0 :       if (s == -1) return gneg(y);
    2958           0 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2959             : 
    2960     1233000 :     case t_VEC: case t_COL: case t_MAT:
    2961     1233000 :       z = cgetg_copy(y, &ly);
    2962     3880987 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2963     1233001 :       return z;
    2964             :   }
    2965           0 :   pari_err_TYPE("gmulsg",y);
    2966             :   return NULL; /* LCOV_EXCL_LINE */
    2967             : }
    2968             : 
    2969             : GEN
    2970   124315520 : gmulug(ulong s, GEN y)
    2971             : {
    2972             :   long ly, i;
    2973             :   pari_sp av;
    2974             :   GEN z;
    2975             : 
    2976   124315520 :   switch(typ(y))
    2977             :   {
    2978   122203851 :     case t_INT:  return mului(s,y);
    2979     1059371 :     case t_REAL: return s? mulur(s,y): gen_0; /* gmul semantic */
    2980         364 :     case t_INTMOD: { GEN p = gel(y,1);
    2981         364 :       z = cgetg(3,t_INTMOD);
    2982         364 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(y,2)), p));
    2983         364 :       gel(z,1) = icopy(p); return z;
    2984             :     }
    2985         413 :     case t_FFELT: return FF_Z_mul(y,utoi(s));
    2986      976139 :     case t_FRAC:
    2987      976139 :       if (!s) return gen_0;
    2988      976125 :       z = cgetg(3,t_FRAC);
    2989      976125 :       i = ugcd(s, umodiu(gel(y,2), s));
    2990      976125 :       if (i == 1)
    2991             :       {
    2992      819549 :         gel(z,2) = icopy(gel(y,2));
    2993      819549 :         gel(z,1) = muliu(gel(y,1), s);
    2994             :       }
    2995             :       else
    2996             :       {
    2997      156576 :         gel(z,2) = diviuexact(gel(y,2), i);
    2998      156576 :         gel(z,1) = muliu(gel(y,1), s/i);
    2999      156576 :         fix_frac_if_int(z);
    3000             :       }
    3001      976125 :       return z;
    3002             : 
    3003       30765 :     case t_COMPLEX:
    3004       30765 :       if (!s) return gen_0;
    3005       30765 :       z = cgetg(3, t_COMPLEX);
    3006       30765 :       gel(z,1) = gmulug(s,gel(y,1));
    3007       30765 :       gel(z,2) = gmulug(s,gel(y,2)); return z;
    3008             : 
    3009        7342 :     case t_PADIC:
    3010        7342 :       if (!s) return gen_0;
    3011        7342 :       av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),y), y));
    3012             : 
    3013           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3014           0 :       gel(z,1) = ZX_copy(gel(y,1));
    3015           0 :       gel(z,2) = gmulug(s,gel(y,2));
    3016           0 :       gel(z,3) = gmulug(s,gel(y,3)); return z;
    3017             : 
    3018        6783 :     case t_POLMOD:
    3019        6783 :       retmkpolmod(gmulug(s,gel(y,2)), RgX_copy(gel(y,1)));
    3020             : 
    3021       16226 :     case t_POL:
    3022       16226 :       if (!signe(y)) return RgX_copy(y);
    3023       15449 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    3024       15449 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    3025       52843 :       for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    3026       15449 :       return normalizepol_lg(z, ly);
    3027             : 
    3028           0 :     case t_SER:
    3029           0 :       if (ser_isexactzero(y)) return gcopy(y);
    3030           0 :       if (!s) return Rg_get_0(y);
    3031           0 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    3032           0 :       for (i=2; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    3033           0 :       return normalizeser(z);
    3034             : 
    3035           0 :     case t_RFRAC:
    3036           0 :       if (!s) return zeropol(varn(gel(y,2)));
    3037           0 :       if (s == 1) return gcopy(y);
    3038           0 :       return mul_rfrac_scal(gel(y,1), gel(y,2), utoi(s));
    3039             : 
    3040       14266 :     case t_VEC: case t_COL: case t_MAT:
    3041       14266 :       z = cgetg_copy(y, &ly);
    3042       74284 :       for (i=1; i<ly; i++) gel(z,i) = gmulug(s,gel(y,i));
    3043       14266 :       return z;
    3044             :   }
    3045           0 :   pari_err_TYPE("gmulsg",y);
    3046             :   return NULL; /* LCOV_EXCL_LINE */
    3047             : }
    3048             : 
    3049             : /********************************************************************/
    3050             : /**                                                                **/
    3051             : /**                       SIMPLE DIVISION                          **/
    3052             : /**                                                                **/
    3053             : /********************************************************************/
    3054             : 
    3055             : GEN
    3056    12911361 : gdivgs(GEN x, long s)
    3057             : {
    3058    12911361 :   long tx = typ(x), lx, i;
    3059             :   pari_sp av;
    3060             :   GEN z;
    3061             : 
    3062    12911361 :   if (!s)
    3063             :   {
    3064           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3065           0 :     pari_err_INV("gdivgs",gen_0);
    3066             :   }
    3067    12911366 :   switch(tx)
    3068             :   {
    3069     1560284 :     case t_INT: return Qdivis(x, s);
    3070     8489141 :     case t_REAL: return divrs(x,s);
    3071             : 
    3072         357 :     case t_INTMOD:
    3073         357 :       z = cgetg(3, t_INTMOD);
    3074         357 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    3075             : 
    3076         735 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    3077             : 
    3078      552008 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3079      552008 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    3080      552008 :       if (i == 1)
    3081             :       {
    3082      409201 :         gel(z,2) = mulsi(s, gel(x,2));
    3083      409201 :         gel(z,1) = icopy(gel(x,1));
    3084             :       }
    3085             :       else
    3086             :       {
    3087      142807 :         gel(z,2) = mulsi(s/i, gel(x,2));
    3088      142807 :         gel(z,1) = divis(gel(x,1), i);
    3089             :       }
    3090      552008 :       normalize_frac(z);
    3091      552008 :       fix_frac_if_int(z); return z;
    3092             : 
    3093     1839911 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3094     1839911 :       gel(z,1) = gdivgs(gel(x,1),s);
    3095     1839911 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3096             : 
    3097         133 :     case t_PADIC: /* divpT */
    3098             :     {
    3099         133 :       GEN p = gel(x,2);
    3100         133 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3101         133 :       av = avma;
    3102         133 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    3103             :     }
    3104             : 
    3105          28 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3106          28 :       gel(z,1) = ZX_copy(gel(x,1));
    3107          28 :       gel(z,2) = gdivgs(gel(x,2),s);
    3108          28 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    3109             : 
    3110       37702 :     case t_POLMOD:
    3111       37702 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    3112             : 
    3113          91 :     case t_RFRAC:
    3114          91 :       if (s == 1) return gcopy(x);
    3115          84 :       else if (s == -1) return gneg(x);
    3116          84 :       return div_rfrac_scal(x, stoi(s));
    3117             : 
    3118       37499 :     case t_POL: case t_SER:
    3119       37499 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3120      157174 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3121       37499 :       return z;
    3122      393477 :     case t_VEC: case t_COL: case t_MAT:
    3123      393477 :       z = cgetg_copy(x, &lx);
    3124      870049 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3125      393477 :       return z;
    3126             : 
    3127             :   }
    3128           0 :   pari_err_TYPE2("/",x, stoi(s));
    3129             :   return NULL; /* LCOV_EXCL_LINE */
    3130             : }
    3131             : 
    3132             : GEN
    3133    54702234 : gdivgu(GEN x, ulong s)
    3134             : {
    3135    54702234 :   long tx = typ(x), lx, i;
    3136             :   pari_sp av;
    3137             :   GEN z;
    3138             : 
    3139    54702234 :   if (!s)
    3140             :   {
    3141           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3142           0 :     pari_err_INV("gdivgu",gen_0);
    3143             :   }
    3144    54702409 :   switch(tx)
    3145             :   {
    3146    17909471 :     case t_INT: return Qdiviu(x, s);
    3147    12950984 :     case t_REAL: return divru(x,s);
    3148             : 
    3149      210315 :     case t_INTMOD:
    3150      210315 :       z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
    3151      210315 :       return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
    3152             : 
    3153         308 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
    3154             : 
    3155      641702 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3156      641702 :       i = ugcd(s, umodiu(gel(x,1), s));
    3157      641702 :       if (i == 1)
    3158             :       {
    3159      461980 :         gel(z,2) = mului(s, gel(x,2));
    3160      461980 :         gel(z,1) = icopy(gel(x,1));
    3161             :       }
    3162             :       else
    3163             :       {
    3164      179722 :         gel(z,2) = mului(s/i, gel(x,2));
    3165      179722 :         gel(z,1) = divis(gel(x,1), i);
    3166             :       }
    3167      641702 :       normalize_frac(z);
    3168      641702 :       fix_frac_if_int(z); return z;
    3169             : 
    3170    11383205 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3171    11383205 :       gel(z,1) = gdivgu(gel(x,1),s);
    3172    11383205 :       gel(z,2) = gdivgu(gel(x,2),s); return z;
    3173             : 
    3174    11551824 :     case t_PADIC: /* divpT */
    3175             :     {
    3176    11551824 :       GEN p = gel(x,2);
    3177    11551824 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3178    11416207 :       av = avma;
    3179    11416207 :       return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
    3180             :     }
    3181             : 
    3182           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3183           0 :       gel(z,1) = ZX_copy(gel(x,1));
    3184           0 :       gel(z,2) = gdivgu(gel(x,2),s);
    3185           0 :       gel(z,3) = gdivgu(gel(x,3),s); return z;
    3186             : 
    3187        1456 :     case t_POLMOD:
    3188        1456 :       retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
    3189             : 
    3190          56 :     case t_RFRAC:
    3191          56 :       if (s == 1) return gcopy(x);
    3192          56 :       return div_rfrac_scal(x, utoi(s));
    3193             : 
    3194       52759 :     case t_POL: case t_SER:
    3195       52759 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3196      221921 :       for (i=2; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
    3197       52759 :       return z;
    3198         329 :     case t_VEC: case t_COL: case t_MAT:
    3199         329 :       z = cgetg_copy(x, &lx);
    3200        1148 :       for (i=1; i<lx; i++) gel(z,i) = gdivgu(gel(x,i),s);
    3201         329 :       return z;
    3202             :   }
    3203           0 :   pari_err_TYPE2("/",x, utoi(s));
    3204             :   return NULL; /* LCOV_EXCL_LINE */
    3205             : }
    3206             : 
    3207             : /* x / (i*(i+1)) */
    3208             : GEN
    3209   223951658 : divrunextu(GEN x, ulong i)
    3210             : {
    3211   223951658 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3212           0 :     return divri(x, muluu(i , i+1));
    3213             :   else
    3214   223951658 :     return divru(x, i*(i+1));
    3215             : }
    3216             : /* x / (i*(i+1)) */
    3217             : GEN
    3218      813628 : gdivgunextu(GEN x, ulong i)
    3219             : {
    3220      813628 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3221           0 :     return gdivgu(x, i*(i+1));
    3222             :   else
    3223      813628 :     return gdiv(x, muluu(i, i+1));
    3224             : }
    3225             : 
    3226             : /* True shift (exact multiplication by 2^n) */
    3227             : GEN
    3228   124701207 : gmul2n(GEN x, long n)
    3229             : {
    3230             :   long lx, i, k, l;
    3231             :   GEN z, a, b;
    3232             : 
    3233   124701207 :   switch(typ(x))
    3234             :   {
    3235    31172189 :     case t_INT:
    3236    31172189 :       if (n>=0) return shifti(x,n);
    3237     5026979 :       if (!signe(x)) return gen_0;
    3238     3340309 :       l = vali(x); n = -n;
    3239     3340354 :       if (n<=l) return shifti(x,-n);
    3240      397660 :       z = cgetg(3,t_FRAC);
    3241      397660 :       gel(z,1) = shifti(x,-l);
    3242      397660 :       gel(z,2) = int2n(n-l); return z;
    3243             : 
    3244    62777824 :     case t_REAL:
    3245    62777824 :       return shiftr(x,n);
    3246             : 
    3247      180925 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3248      180925 :       z = cgetg(3,t_INTMOD);
    3249      180926 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3250       82765 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3251       82758 :       gel(z,1) = icopy(b); return z;
    3252             : 
    3253      217416 :     case t_FFELT: return FF_mul2n(x,n);
    3254             : 
    3255     1265583 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3256     1265583 :       l = vali(a);
    3257     1265583 :       k = vali(b);
    3258     1265583 :       if (n+l >= k)
    3259             :       {
    3260      393969 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3261      301459 :         l = n-k; k = -k;
    3262             :       }
    3263             :       else
    3264             :       {
    3265      871614 :         k = -(l+n); l = -l;
    3266             :       }
    3267     1173073 :       z = cgetg(3,t_FRAC);
    3268     1173073 :       gel(z,1) = shifti(a,l);
    3269     1173073 :       gel(z,2) = shifti(b,k); return z;
    3270             : 
    3271    10836088 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3272    10836085 :       gel(z,1) = gmul2n(gel(x,1),n);
    3273    10836085 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3274             : 
    3275         105 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3276         105 :       gel(z,1) = ZX_copy(gel(x,1));
    3277         105 :       gel(z,2) = gmul2n(gel(x,2),n);
    3278         105 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3279             : 
    3280      193774 :     case t_POLMOD:
    3281      193774 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3282             : 
    3283     1713106 :     case t_POL:
    3284     1713106 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3285     9089609 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3286     1713111 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3287      102483 :     case t_SER:
    3288      102483 :       if (ser_isexactzero(x)) return gcopy(x);
    3289      102455 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3290      632422 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3291      102455 :       return normalizeser(z); /* needed if char = 2 */
    3292    16237980 :     case t_VEC: case t_COL: case t_MAT:
    3293    16237980 :       z = cgetg_copy(x, &lx);
    3294    59817502 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3295    16237932 :       return z;
    3296             : 
    3297          21 :     case t_RFRAC: /* int2n wrong if n < 0 */
    3298          21 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3299             : 
    3300        5320 :     case t_PADIC: /* int2n wrong if n < 0 */
    3301        5320 :       return gmul(gmul2n(gen_1,n),x);
    3302             :   }
    3303           0 :   pari_err_TYPE("gmul2n",x);
    3304             :   return NULL; /* LCOV_EXCL_LINE */
    3305             : }
    3306             : 
    3307             : /*******************************************************************/
    3308             : /*                                                                 */
    3309             : /*                              INVERSE                            */
    3310             : /*                                                                 */
    3311             : /*******************************************************************/
    3312             : static GEN
    3313      208553 : inv_polmod(GEN T, GEN x)
    3314             : {
    3315      208553 :   GEN z = cgetg(3,t_POLMOD), a;
    3316      208554 :   gel(z,1) = RgX_copy(T);
    3317      208552 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3318       83933 :     a = ginv(x);
    3319             :   else
    3320             :   {
    3321      124619 :     if (lg(T) == 5) /* quadratic fields */
    3322       13048 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3323             :     else
    3324      111571 :       a = RgXQ_inv(x, T);
    3325             :   }
    3326      208550 :   gel(z,2) = a; return z;
    3327             : }
    3328             : 
    3329             : GEN
    3330    36124327 : ginv(GEN x)
    3331             : {
    3332             :   long s;
    3333             :   pari_sp av, tetpil;
    3334             :   GEN z, y, p1, p2;
    3335             : 
    3336    36124327 :   switch(typ(x))
    3337             :   {
    3338     9481049 :     case t_INT:
    3339     9481049 :       if (is_pm1(x)) return icopy(x);
    3340     7890382 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3341     7890368 :       z = cgetg(3,t_FRAC);
    3342     7895411 :       gel(z,1) = s<0? gen_m1: gen_1;
    3343     7895411 :       gel(z,2) = absi(x); return z;
    3344             : 
    3345    10374795 :     case t_REAL: return invr(x);
    3346             : 
    3347       12964 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3348       12964 :       gel(z,1) = icopy(gel(x,1));
    3349       12964 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3350             : 
    3351      550981 :     case t_FRAC: {
    3352      550981 :       GEN a = gel(x,1), b = gel(x,2);
    3353      550981 :       s = signe(a);
    3354      550981 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3355      248088 :       z = cgetg(3,t_FRAC);
    3356      248089 :       gel(z,1) = icopy(b);
    3357      248088 :       gel(z,2) = icopy(a);
    3358      248089 :       normalize_frac(z); return z;
    3359             :     }
    3360     9636605 :     case t_COMPLEX:
    3361     9636605 :       av=avma;
    3362     9636605 :       p1=cxnorm(x);
    3363     9635802 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3364     9636459 :       tetpil=avma;
    3365     9636459 :       return gerepile(av,tetpil,divcR(p2,p1));
    3366             : 
    3367         273 :     case t_QUAD:
    3368         273 :       av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
    3369         273 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3370             : 
    3371      358085 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3372      358085 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3373      358078 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3374      358078 :       gel(z,2) = icopy(gel(x,2));
    3375      358078 :       gel(z,3) = icopy(gel(x,3));
    3376      358078 :       gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
    3377             : 
    3378      208553 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3379       14220 :     case t_FFELT: return FF_inv(x);
    3380     5322786 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3381       26327 :     case t_SER: return ser_inv(x);
    3382        2926 :     case t_RFRAC:
    3383             :     {
    3384        2926 :       GEN n = gel(x,1), d = gel(x,2);
    3385        2926 :       pari_sp av = avma, ltop;
    3386        2926 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3387             : 
    3388        2926 :       n = simplify_shallow(n);
    3389        2926 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3390             :       {
    3391        2926 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3392         679 :         ltop = avma;
    3393         679 :         z = RgX_Rg_div(d,n);
    3394             :       } else {
    3395           0 :         ltop = avma;
    3396           0 :         z = cgetg(3,t_RFRAC);
    3397           0 :         gel(z,1) = RgX_copy(d);
    3398           0 :         gel(z,2) = RgX_copy(n);
    3399             :       }
    3400         679 :       stackdummy(av, ltop);
    3401         679 :       return z;
    3402             :     }
    3403             : 
    3404          21 :     case t_VEC: if (!is_ext_qfr(x)) break;
    3405             :     case t_QFB:
    3406       99470 :       return qfbpow(x, gen_m1);
    3407       36836 :     case t_MAT:
    3408       36836 :       y = RgM_inv(x);
    3409       36829 :       if (!y) pari_err_INV("ginv",x);
    3410       36759 :       return y;
    3411          28 :     case t_VECSMALL:
    3412             :     {
    3413          28 :       long i, lx = lg(x)-1;
    3414          28 :       y = zero_zv(lx);
    3415         112 :       for (i=1; i<=lx; i++)
    3416             :       {
    3417          84 :         long xi = x[i];
    3418          84 :         if (xi<1 || xi>lx || y[xi])
    3419           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3420          84 :         y[xi] = i;
    3421             :       }
    3422          28 :       return y;
    3423             :     }
    3424             :   }
    3425           6 :   pari_err_TYPE("inverse",x);
    3426             :   return NULL; /* LCOV_EXCL_LINE */
    3427             : }

Generated by: LCOV version 1.16