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 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - gen1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 22875-74571cde7) Lines: 1713 1815 94.4 %
Date: 2018-07-21 05:35:54 Functions: 89 89 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                      GENERIC OPERATIONS                        **/
      17             : /**                         (first part)                           **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /* assume z[1] was created last */
      24             : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
      25             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));
      26             : 
      27             : /* assume z[1] was created last */
      28             : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
      29             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));\
      30             : else\
      31             :   gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
      32             : 
      33             : static void
      34          98 : warn_coercion(GEN x, GEN y, GEN z)
      35             : {
      36          98 :   if (DEBUGLEVEL)
      37          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      38          98 : }
      39             : 
      40             : static long
      41          21 : kro_quad(GEN x, GEN y)
      42             : {
      43          21 :   pari_sp av=avma;
      44          21 :   long k = kronecker(quad_disc(x), y);
      45          21 :   set_avma(av); return k;
      46             : }
      47             : 
      48             : /* is -1 not a square in Zp, assume p prime */
      49             : INLINE int
      50          28 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      51             : 
      52             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      53             : static GEN addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN));
      54             : static GEN mulpp(GEN x, GEN y);
      55             : static GEN divpp(GEN x, GEN y);
      56             : /* Argument codes for inline routines
      57             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      58             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      59             :  * T: some type (to be converted to PADIC)
      60             :  */
      61             : static GEN
      62    72296541 : addRc(GEN x, GEN y) {
      63    72296541 :   GEN z = cgetg(3,t_COMPLEX);
      64    72296541 :   gel(z,1) = gadd(x,gel(y,1));
      65    72296541 :   gel(z,2) = gcopy(gel(y,2)); return z;
      66             : }
      67             : static GEN
      68    69833840 : mulRc(GEN x, GEN y) {
      69    69833840 :   GEN z = cgetg(3,t_COMPLEX);
      70    69833840 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      71    69833840 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      72             : }
      73             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      74             : static GEN
      75          49 : mulRc_direct(GEN x, GEN y) {
      76          49 :   GEN z = cgetg(3,t_COMPLEX);
      77          49 :   gel(z,1) = gmul(x,gel(y,1));
      78          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      79             : }
      80             : static GEN
      81      488630 : divRc(GEN x, GEN y) {
      82      488630 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      83      488630 :   GEN z = cgetg(3,t_COMPLEX);
      84      488630 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      85      488630 :   gel(z,2) = gmul(mt, gel(y,2));
      86      488630 :   return z;
      87             : }
      88             : static GEN
      89     4667644 : divcR(GEN x, GEN y) {
      90     4667644 :   GEN z = cgetg(3,t_COMPLEX);
      91     4667644 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      92     4667644 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      93             : }
      94             : static GEN
      95         756 : addRq(GEN x, GEN y) {
      96         756 :   GEN z = cgetg(4,t_QUAD);
      97         756 :   gel(z,1) = ZX_copy(gel(y,1));
      98         756 :   gel(z,2) = gadd(x, gel(y,2));
      99         756 :   gel(z,3) = gcopy(gel(y,3)); return z;
     100             : }
     101             : static GEN
     102        1750 : mulRq(GEN x, GEN y) {
     103        1750 :   GEN z = cgetg(4,t_QUAD);
     104        1750 :   gel(z,1) = ZX_copy(gel(y,1));
     105        1750 :   gel(z,2) = gmul(x,gel(y,2));
     106        1750 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     107             : }
     108             : static GEN
     109          28 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     110          28 :   long i = gexpo(x) - gexpo(y);
     111          28 :   if (i > 0) prec += nbits2extraprec( i );
     112          28 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     113             : }
     114             : static GEN
     115     4461493 : mulrfrac(GEN x, GEN y)
     116             : {
     117     4461493 :   pari_sp av = avma;
     118     4461493 :   GEN z, a = gel(y,1), b = gel(y,2);
     119     4461493 :   if (is_pm1(a)) /* frequent special case */
     120             :   {
     121     1851566 :     z = divri(x, b);
     122     1851566 :     if (signe(a) < 0) togglesign(z);
     123     1851566 :     return z;
     124             :   }
     125     2609927 :   return gerepileuptoleaf(av, divri(mulri(x,gel(y,1)), gel(y,2)));
     126             : }
     127             : static GEN
     128          14 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     129          14 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     130             : }
     131             : static GEN
     132          28 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     133          28 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     134             : }
     135             : static GEN
     136           7 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     137           7 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     138             : }
     139             : /* y PADIC, x + y by converting x to padic */
     140             : static GEN
     141           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     142             : 
     143           7 :   if (!valp(y)) z = cvtop2(x,y);
     144             :   else {
     145           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     146           7 :     z  = cvtop(x, gel(y,2), l);
     147             :   }
     148           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     149             : }
     150             : /* y PADIC, x * y by converting x to padic */
     151             : static GEN
     152      162721 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     153      162721 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     154             : }
     155             : /* y PADIC, non zero x / y by converting x to padic */
     156             : static GEN
     157       14119 : divTp(GEN x, GEN y) { pari_sp av = avma;
     158       14119 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     159             : }
     160             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     161             :  * converted to O(p^e) and division by 0 */
     162             : static GEN
     163        1820 : divpT(GEN x, GEN y) { pari_sp av = avma;
     164        1820 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     165             : }
     166             : 
     167             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     168             :  * clean memory from z on */
     169             : static GEN
     170     2235945 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     171     2235945 :   if (lgefint(X) == 3) {
     172     2183200 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     173     2183200 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     174             :   }
     175             :   else {
     176       52745 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     177       52729 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     178             :   }
     179     2235948 :   gel(z,1) = icopy(X); return z;
     180             : }
     181             : static GEN
     182     1102001 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     183     1102001 :   if (lgefint(X) == 3) {
     184      751033 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     185      751033 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     186             :   }
     187             :   else {
     188      350968 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     189      350963 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     190             :   }
     191     1102005 :   gel(z,1) = icopy(X); return z;
     192             : }
     193             : /* cf add_intmod_same */
     194             : static GEN
     195     2689239 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     196     2689239 :   if (lgefint(X) == 3) {
     197     2592477 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     198     2592477 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     199             :   }
     200             :   else
     201       96762 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     202     2689253 :   gel(z,1) = icopy(X); return z;
     203             : }
     204             : /* cf add_intmod_same */
     205             : static GEN
     206       59522 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     207             : {
     208       59522 :   if (lgefint(X) == 3) {
     209       51296 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     210       51289 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     211             :   }
     212             :   else
     213        8226 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     214       59517 :   gel(z,1) = icopy(X); return z;
     215             : }
     216             : 
     217             : /*******************************************************************/
     218             : /*                                                                 */
     219             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     220             : /*                                                                 */
     221             : /* (static routines are not memory clean, but OK for gerepileupto) */
     222             : /*******************************************************************/
     223             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     224             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     225             : static GEN
     226    10074864 : rfrac_denom_mul_scal(GEN d, GEN y)
     227             : {
     228    10074864 :   GEN D = RgX_Rg_mul(d, y);
     229    10074864 :   if (lg(D) != lg(d))
     230             :   { /* try to generate a meaningful diagnostic */
     231           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     232           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     233             :   }
     234    10074864 :   return D;
     235             : }
     236             : 
     237             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     238             : GEN
     239    58058679 : gred_rfrac_simple(GEN n, GEN d)
     240             : {
     241             :   GEN c, cn, cd, z;
     242    58058679 :   long dd = degpol(d);
     243             : 
     244    58058679 :   if (dd <= 0)
     245             :   {
     246       49427 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     247       49427 :     n = gdiv(n, gel(d,2));
     248       49427 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     249       49427 :     return n;
     250             :   }
     251             : 
     252    58009252 :   cd = content(d);
     253    58009252 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     254    58009252 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     255    58009252 :   if (!gequal1(cd)) {
     256     6576369 :     d = RgX_Rg_div(d,cd);
     257     6576369 :     if (!gequal1(cn))
     258             :     {
     259     1302563 :       if (gequal0(cn)) {
     260          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     261           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     262           0 :         c = gen_1;
     263             :       } else {
     264     1302514 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     265     1302514 :         c = gdiv(cn,cd);
     266             :       }
     267             :     }
     268             :     else
     269     5273806 :       c = ginv(cd);
     270             :   } else {
     271    51432883 :     if (!gequal1(cn))
     272             :     {
     273     3227914 :       if (gequal0(cn)) {
     274         413 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     275          21 :         c = gen_1;
     276             :       } else {
     277     3227501 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     278     3227501 :         c = cn;
     279             :       }
     280             :     } else {
     281    48204969 :       GEN y = cgetg(3,t_RFRAC);
     282    48204969 :       gel(y,1) = gcopy(n);
     283    48204969 :       gel(y,2) = RgX_copy(d); return y;
     284             :     }
     285             :   }
     286             : 
     287     9803842 :   if (typ(c) == t_POL)
     288             :   {
     289      905469 :     z = c;
     290      944886 :     do { z = content(z); } while (typ(z) == t_POL);
     291      905469 :     cd = denom_i(z);
     292      905469 :     cn = gmul(c, cd);
     293             :   }
     294             :   else
     295             :   {
     296     8898373 :     cn = numer_i(c);
     297     8898373 :     cd = denom_i(c);
     298             :   }
     299     9803842 :   z = cgetg(3,t_RFRAC);
     300     9803842 :   gel(z,1) = gmul(n, cn);
     301     9803842 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     302     9803842 :   return z;
     303             : }
     304             : 
     305             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     306             : static GEN
     307      249067 : fix_rfrac(GEN x, long d)
     308             : {
     309             :   GEN z, N, D;
     310      249067 :   if (!d || typ(x) == t_POL) return x;
     311      175623 :   z = cgetg(3, t_RFRAC);
     312      175623 :   N = gel(x,1);
     313      175623 :   D = gel(x,2);
     314      175623 :   if (d > 0) {
     315      344253 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     316      177968 :                                                    : monomialcopy(N,d,varn(D));
     317      166236 :     gel(z, 2) = RgX_copy(D);
     318             :   } else {
     319        9387 :     gel(z, 1) = gcopy(N);
     320        9387 :     gel(z, 2) = RgX_shift(D, -d);
     321             :   }
     322      175623 :   return z;
     323             : }
     324             : 
     325             : /* assume d != 0 */
     326             : static GEN
     327    44746556 : gred_rfrac2(GEN n, GEN d)
     328             : {
     329             :   GEN y, z;
     330             :   long v, vd, vn;
     331             : 
     332    44746556 :   n = simplify_shallow(n);
     333    44746556 :   if (isrationalzero(n)) return gcopy(n);
     334    37918598 :   d = simplify_shallow(d);
     335    37918598 :   if (typ(d) != t_POL) return gdiv(n,d);
     336    36738942 :   vd = varn(d);
     337    36738942 :   if (typ(n) != t_POL)
     338             :   {
     339    20613081 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     340    20611674 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     341           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     342             :   }
     343    16125861 :   vn = varn(n);
     344    16125861 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     345    15988243 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     346             : 
     347             :   /* now n and d are t_POLs in the same variable */
     348    15828491 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     349    15828491 :   if (!degpol(d))
     350             :   {
     351    12755800 :     n = RgX_Rg_div(n,gel(d,2));
     352    12755800 :     return v? RgX_mulXn(n,v): n;
     353             :   }
     354             : 
     355             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     356     3072691 :   if (!isinexact(n) && !isinexact(d))
     357             :   {
     358     3072481 :     y = RgX_divrem(n, d, &z);
     359     3072481 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     360      248857 :     z = RgX_gcd(d, z);
     361      248857 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     362             :   }
     363      249067 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     364             : }
     365             : 
     366             : /* x,y t_INT, return x/y in reduced form */
     367             : GEN
     368    66443463 : Qdivii(GEN x, GEN y)
     369             : {
     370    66443463 :   pari_sp av = avma;
     371             :   GEN r, q;
     372             : 
     373    66443463 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     374    30155782 :   if (equali1(x)) {
     375     5308145 :     long s = signe(y);
     376     5308145 :     if (!s) pari_err_INV("gdiv",y);
     377     5308040 :     if (signe(x) < 0) s = -s;
     378     5308040 :     q = cgetg(3, t_FRAC);
     379     5308040 :     gel(q,1) = s<0? gen_m1: gen_1;
     380     5308040 :     gel(q,2) = absi(y); return q;
     381             :   }
     382    24847637 :   q = dvmdii(x,y,&r);
     383    24847637 :   if (r == gen_0) return q; /* gen_0 intended */
     384     9443031 :   r = gcdii(y,r);
     385     9443031 :   if (lgefint(r) == 3)
     386             :   {
     387     9219541 :     ulong t = r[2];
     388     9219541 :     set_avma(av);
     389     9219541 :     if (t == 1) q = mkfraccopy(x,y);
     390             :     else
     391             :     {
     392     4591725 :       q = cgetg(3,t_FRAC);
     393     4591725 :       gel(q,1) = diviuexact(x,t);
     394     4591725 :       gel(q,2) = diviuexact(y,t);
     395             :     }
     396             :   }
     397             :   else
     398             :   { /* rare: r and q left on stack for efficiency */
     399      223490 :     q = cgetg(3,t_FRAC);
     400      223490 :     gel(q,1) = diviiexact(x,r);
     401      223490 :     gel(q,2) = diviiexact(y,r);
     402             :   }
     403     9443031 :   normalize_frac(q); return q;
     404             : }
     405             : 
     406             : /*******************************************************************/
     407             : /*                                                                 */
     408             : /*                          CONJUGATION                            */
     409             : /*                                                                 */
     410             : /*******************************************************************/
     411             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     412             : static GEN
     413       23877 : quad_polmod_conj(GEN x, GEN y)
     414             : {
     415             :   GEN z, u, v, a, b;
     416       23877 :   if (typ(x) != t_POL) return gcopy(x);
     417       23877 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     418       23877 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     419       23877 :   b = gel(y,3); v = gel(x,2);
     420       23877 :   z = cgetg(4, t_POL); z[1] = x[1];
     421       23877 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     422       23877 :   gel(z,3) = gneg(u); return z;
     423             : }
     424             : static GEN
     425       23877 : quad_polmod_norm(GEN x, GEN y)
     426             : {
     427             :   GEN z, u, v, a, b, c;
     428       23877 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     429       23877 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     430       23877 :   b = gel(y,3); v = gel(x,2);
     431       23877 :   c = gel(y,2);
     432       23877 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     433       23877 :   if (!gequal1(a)) z = gdiv(z, a);
     434       23877 :   return gadd(z, gsqr(v));
     435             : }
     436             : 
     437             : GEN
     438     8220983 : conj_i(GEN x)
     439             : {
     440             :   long lx, i;
     441             :   GEN y;
     442             : 
     443     8220983 :   switch(typ(x))
     444             :   {
     445             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     446     4379238 :       return x;
     447             : 
     448     3812580 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     449             : 
     450             :     case t_QUAD:
     451         595 :       y = cgetg(4,t_QUAD);
     452         595 :       gel(y,1) = gel(x,1);
     453        1190 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     454         595 :                                       : gadd(gel(x,2), gel(x,3));
     455         595 :       gel(y,3) = gneg(gel(x,3)); return y;
     456             : 
     457             :     case t_POL: case t_SER:
     458       10507 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     459       10507 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     460       10507 :       return y;
     461             : 
     462             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     463       18063 :       y = cgetg_copy(x, &lx);
     464       18063 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     465       18063 :       return y;
     466             : 
     467             :     case t_POLMOD:
     468             :     {
     469           0 :       GEN X = gel(x,1);
     470           0 :       long d = degpol(X);
     471           0 :       if (d < 2) return x;
     472           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     473             :     }
     474             :   }
     475           0 :   pari_err_TYPE("gconj",x);
     476             :   return NULL; /* LCOV_EXCL_LINE */
     477             : }
     478             : GEN
     479       29521 : gconj(GEN x)
     480       29521 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     481             : 
     482             : GEN
     483          84 : conjvec(GEN x,long prec)
     484             : {
     485             :   long lx, s, i;
     486             :   GEN z;
     487             : 
     488          84 :   switch(typ(x))
     489             :   {
     490             :     case t_INT: case t_INTMOD: case t_FRAC:
     491           0 :       return mkcolcopy(x);
     492             : 
     493             :     case t_COMPLEX: case t_QUAD:
     494           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     495             : 
     496             :     case t_FFELT:
     497          28 :       return FF_conjvec(x);
     498             : 
     499             :     case t_VEC: case t_COL:
     500           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     501           0 :       if (lx == 1) return z;
     502           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     503           0 :       s = lgcols(z);
     504           0 :       for (i=2; i<lx; i++)
     505             :       {
     506           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     507           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     508             :       }
     509           0 :       break;
     510             : 
     511             :     case t_POLMOD: {
     512          56 :       GEN T = gel(x,1), r;
     513             :       pari_sp av;
     514             : 
     515          56 :       lx = lg(T);
     516          56 :       if (lx <= 3) return cgetg(1,t_COL);
     517          56 :       x = gel(x,2);
     518         238 :       for (i=2; i<lx; i++)
     519             :       {
     520         189 :         GEN c = gel(T,i);
     521         189 :         switch(typ(c)) {
     522             :           case t_INTMOD: {
     523           7 :             GEN p = gel(c,1);
     524             :             pari_sp av;
     525           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     526           7 :             av = avma;
     527           7 :             T = RgX_to_FpX(T,p);
     528           7 :             x = RgX_to_FpX(x, p);
     529           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     530           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     531           7 :             return gerepileupto(av, z);
     532             :           }
     533             :           case t_INT:
     534         182 :           case t_FRAC: break;
     535           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     536             :         }
     537             :       }
     538          49 :       if (typ(x) != t_POL)
     539             :       {
     540           0 :         if (!is_rational_t(typ(x)))
     541           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     542           0 :         retconst_col(lx-3, gcopy(x));
     543             :       }
     544          49 :       RgX_check_QX(x,"conjvec");
     545          49 :       av = avma;
     546          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     547          49 :       r = cleanroots(T,prec);
     548          49 :       z = cgetg(lx-2,t_COL);
     549          49 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     550          49 :       return gerepileupto(av, z);
     551             :     }
     552             : 
     553             :     default:
     554           0 :       pari_err_TYPE("conjvec",x);
     555             :       return NULL; /* LCOV_EXCL_LINE */
     556             :   }
     557           0 :   return z;
     558             : }
     559             : 
     560             : 
     561             : /********************************************************************/
     562             : /**                                                                **/
     563             : /**                           ADDITION                             **/
     564             : /**                                                                **/
     565             : /********************************************************************/
     566             : /* x, y compatible PADIC, op = add or sub */
     567             : static GEN
     568     1178489 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     569             : {
     570     1178489 :   pari_sp av = avma;
     571             :   long d,e,r,rx,ry;
     572     1178489 :   GEN u, z, mod, p = gel(x,2);
     573             :   int swap;
     574             : 
     575     1178489 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     576     1178489 :   e = valp(x);
     577     1178489 :   r = valp(y); d = r-e;
     578     1178489 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     579     1178489 :   rx = precp(x);
     580     1178489 :   ry = precp(y);
     581     1178489 :   if (d) /* v(x) < v(y) */
     582             :   {
     583      548054 :     r = d+ry; z = powiu(p,d);
     584      548054 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     585      548054 :     z = mulii(z,gel(y,4));
     586      548054 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     587             :   }
     588             :   else
     589             :   {
     590             :     long c;
     591      630435 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     592      630435 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     593      630435 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     594             :     {
     595       59171 :       set_avma(av); return zeropadic(p, e+r);
     596             :     }
     597      571264 :     if (c)
     598             :     {
     599       95333 :       mod = diviiexact(mod, powiu(p,c));
     600       95333 :       r -= c;
     601       95333 :       e += c;
     602             :     }
     603             :   }
     604     1119318 :   u = modii(u, mod);
     605     1119318 :   set_avma(av); z = cgetg(5,t_PADIC);
     606     1119318 :   z[1] = evalprecp(r) | evalvalp(e);
     607     1119318 :   gel(z,2) = icopy(p);
     608     1119318 :   gel(z,3) = icopy(mod);
     609     1119318 :   gel(z,4) = icopy(u); return z;
     610             : }
     611             : /* Rg_to_Fp(t_FRAC) without GC */
     612             : static GEN
     613       56483 : Q_to_Fp(GEN x, GEN p)
     614       56483 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     615             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     616             : static GEN
     617      370495 : addQp(GEN x, GEN y)
     618             : {
     619      370495 :   pari_sp av = avma;
     620      370495 :   long d, r, e, vy = valp(y), py = precp(y);
     621      370495 :   GEN z, q, mod, u, p = gel(y,2);
     622             : 
     623      370495 :   e = Q_pvalrem(x, p, &x);
     624      370463 :   d = vy - e; r = d + py;
     625      370463 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     626      368881 :   mod = gel(y,3);
     627      368881 :   u   = gel(y,4);
     628      368881 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     629             : 
     630      368913 :   if (d > 0)
     631             :   {
     632      234705 :     q = powiu(p,d);
     633      234660 :     mod = mulii(mod, q);
     634      234653 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     635      234653 :     u = addii(x,  mulii(u, q));
     636             :   }
     637      134208 :   else if (d < 0)
     638             :   {
     639       25198 :     q = powiu(p,-d);
     640       25198 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     641       25198 :     u = addii(u, mulii(x, q));
     642       25194 :     r = py; e = vy;
     643             :   }
     644             :   else
     645             :   {
     646             :     long c;
     647      109010 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     648      109010 :     u = addii(u, x);
     649      109010 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     650             :     {
     651         966 :       set_avma(av); return zeropadic(p,e+r);
     652             :     }
     653      108044 :     if (c)
     654             :     {
     655       58821 :       mod = diviiexact(mod, powiu(p,c));
     656       58821 :       r -= c;
     657       58821 :       e += c;
     658             :     }
     659             :   }
     660      367894 :   u = modii(u, mod); set_avma(av);
     661      367904 :   z = cgetg(5,t_PADIC);
     662      367935 :   z[1] = evalprecp(r) | evalvalp(e);
     663      367939 :   gel(z,2) = icopy(p);
     664      367943 :   gel(z,3) = icopy(mod);
     665      367945 :   gel(z,4) = icopy(u); return z;
     666             : }
     667             : 
     668             : /* Mod(x,X) + Mod(y,X) */
     669             : #define addsub_polmod_same addsub_polmod_scal
     670             : /* Mod(x,X) +/- Mod(y,Y) */
     671             : static GEN
     672        1750 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     673             : {
     674        1750 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     675        1750 :   GEN z = cgetg(3,t_POLMOD);
     676        1750 :   long vx = varn(X), vy = varn(Y);
     677        1750 :   if (vx==vy) {
     678             :     pari_sp av;
     679          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     680          14 :     warn_coercion(X,Y,gel(z,1));
     681          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     682             :   }
     683        1736 :   if (varncmp(vx, vy) < 0)
     684        1736 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     685             :   else
     686           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     687        1736 :   gel(z,2) = op(x, y); return z;
     688             : }
     689             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     690             : static GEN
     691     1979881 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     692     1979881 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     693             : 
     694             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     695             : static GEN
     696      308833 : add_ser_scal(GEN y, GEN x)
     697             : {
     698             :   long i, l, ly, vy;
     699             :   GEN z;
     700             : 
     701      308833 :   if (isrationalzero(x)) return gcopy(y);
     702      282534 :   ly = lg(y);
     703      282534 :   l = valp(y);
     704      282534 :   if (l < 3-ly) return gcopy(y);
     705             :   /* l + ly >= 3 */
     706      282366 :   if (l < 0)
     707             :   {
     708         735 :     z = cgetg(ly,t_SER); z[1] = y[1];
     709         735 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     710         735 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     711         735 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     712         735 :     return z;
     713             :   }
     714      281631 :   vy = varn(y);
     715      281631 :   if (l > 0)
     716             :   {
     717       21315 :     if (ser_isexactzero(y))
     718       13580 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     719        7735 :     y -= l; ly += l;
     720        7735 :     z = cgetg(ly,t_SER);
     721        7735 :     x = gcopy(x);
     722        7735 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     723             :   }
     724             :   else
     725             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     726      260316 :     z = cgetg(ly,t_SER);
     727      260316 :     x = gadd(x, gel(y,2));
     728      260316 :     i = 3;
     729             :   }
     730      268051 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     731      268051 :   gel(z,2) = x;
     732      268051 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     733      268051 :   return gequal0(x)? normalize(z): z;
     734             : }
     735             : static long
     736     3538960 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     737             : /* x,y t_SER in the same variable: x+y */
     738             : static GEN
     739     1769872 : ser_add(GEN x, GEN y)
     740             : {
     741     1769872 :   long i, lx,ly, n = valp(y) - valp(x);
     742             :   GEN z;
     743     1769872 :   if (n < 0) { n = -n; swap(x,y); }
     744             :   /* valp(x) <= valp(y) */
     745     1769872 :   lx = _serprec(x);
     746     1769872 :   if (lx == 2) /* don't lose type information */
     747         784 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
     748     1769088 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     749     1769088 :   if (n)
     750             :   {
     751       90905 :     if (n+2 > lx) return gcopy(x);
     752       90352 :     z = cgetg(ly,t_SER);
     753       90352 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     754       90352 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     755             :   } else {
     756     1678183 :     z = cgetg(ly,t_SER);
     757     1678183 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     758             :   }
     759     1768535 :   z[1] = x[1]; return normalize(z);
     760             : }
     761             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     762             : static GEN
     763     8824253 : add_rfrac_scal(GEN y, GEN x)
     764             : {
     765             :   pari_sp av;
     766             :   GEN n;
     767             : 
     768     8824253 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     769     5147142 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     770     5147142 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     771             : }
     772             : 
     773             : /* x "scalar", ty != t_MAT and non-scalar */
     774             : static GEN
     775    16082630 : add_scal(GEN y, GEN x, long ty)
     776             : {
     777    16082630 :   switch(ty)
     778             :   {
     779    11088586 :     case t_POL: return RgX_Rg_add(y, x);
     780      308819 :     case t_SER: return add_ser_scal(y, x);
     781     4474973 :     case t_RFRAC: return add_rfrac_scal(y, x);
     782           0 :     case t_COL: return RgC_Rg_add(y, x);
     783             :     case t_VEC:
     784      210245 :       if (isintzero(x)) return gcopy(y);
     785         168 :       break;
     786             :   }
     787         175 :   pari_err_TYPE2("+",x,y);
     788             :   return NULL; /* LCOV_EXCL_LINE */
     789             : }
     790             : 
     791             : static GEN
     792    15908981 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     793             : {
     794    15908981 :   pari_sp av = avma;
     795    15908981 :   GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
     796    15908981 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     797    15908981 :   int s = cmpii(x2, y2);
     798             : 
     799    15908981 :   if (!s)
     800             :   { /* common denominator: (x1 op y1) / x2 */
     801     5472935 :     n = op(x1, y1);
     802     5472935 :     if (!signe(n)) { set_avma(av); return gen_0; }
     803     5005422 :     d = x2;
     804     5005422 :     q = dvmdii(n, d, &r);
     805     5005422 :     if (r == gen_0) { set_avma(av); return icopy(q); }
     806     3448915 :     r = gcdii(d, r);
     807     3448915 :     if (!equali1(r)) { n = diviiexact(n, r); d = diviiexact(d, r); }
     808     3448915 :     gel(z,1) = icopy_avma(n, (pari_sp)z);
     809     3448915 :     gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     810     3448915 :     avma = (pari_sp)gel(z,2); return z;
     811             :   }
     812    10436046 :   if (s < 0)
     813             :   {
     814     4348381 :     GEN Q = dvmdii(y2, x2, &r);
     815     4348381 :     if (r == gen_0)
     816             :     { /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     817     3260206 :       pari_sp av = avma;
     818     3260206 :       n = op(mulii(Q,x1), y1);
     819     3260206 :       q = dvmdii(n, x2, &r);
     820     3260206 :       if (r == gen_0)
     821             :       {
     822      609054 :         gel(z,1) = gerepileuptoint(av, q);
     823      609054 :         gel(z,2) = Q; return z;
     824             :       }
     825     2651152 :       r = gcdii(x2, r);
     826     2651152 :       if (!equali1(r)) { n = diviiexact(n, r); x2 = diviiexact(x2, r); }
     827     2651152 :       d = mulii(x2,Q);
     828     2651152 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     829     2651152 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     830     2651152 :       avma = (pari_sp)gel(z,2); return z;
     831             :     }
     832     1088175 :     delta = gcdii(x2,r);
     833             :   }
     834             :   else
     835             :   {
     836     6087665 :     GEN Q = dvmdii(x2, y2, &r);
     837     6087665 :     if (r == gen_0)
     838             :     { /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     839     4998031 :       pari_sp av = avma;
     840     4998031 :       n = op(x1, mulii(Q,y1));
     841     4998031 :       q = dvmdii(n, y2, &r);
     842     4998031 :       if (r == gen_0)
     843             :       {
     844      131035 :         gel(z,1) = gerepileuptoint(av, q);
     845      131035 :         gel(z,2) = Q; return z;
     846             :       }
     847     4866996 :       r = gcdii(y2, r);
     848     4866996 :       if (!equali1(r)) { n = diviiexact(n, r); y2 = diviiexact(y2, r); }
     849     4866996 :       d = mulii(y2,Q);
     850     4866996 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     851     4866996 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     852     4866996 :       avma = (pari_sp)gel(z,2); return z;
     853             :     }
     854     1089634 :     delta = gcdii(y2,r);
     855             :   }
     856             :   /* delta = gcd(x2,y2) */
     857     2177809 :   if (equali1(delta))
     858             :   { /* numerator is non-zero */
     859      383243 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     860      383243 :     gel(z,2) = mulii(x2,y2); return z;
     861             :   }
     862     1794566 :   x2 = diviiexact(x2,delta);
     863     1794566 :   y2 = diviiexact(y2,delta);
     864     1794566 :   n = op(mulii(x1,y2), mulii(y1,x2));
     865     1794566 :   if (!signe(n)) { set_avma(av); return gen_0; }
     866     1794566 :   d = mulii(x2, y2);
     867     1794566 :   q = dvmdii(n, delta, &r);
     868     1794566 :   if (r == gen_0)
     869             :   {
     870      122149 :     if (equali1(d)) { set_avma(av); return icopy(q); }
     871      122149 :     avma = (pari_sp)z;
     872      122149 :     gel(z,2) = icopy(d);
     873      122149 :     gel(z,1) = icopy(q); return z;
     874             :   }
     875     1672417 :   r = gcdii(delta, r);
     876     1672417 :   if (!equali1(r))
     877             :   {
     878      401557 :     n     = diviiexact(n, r);
     879      401557 :     delta = diviiexact(delta, r);
     880             :   }
     881     1672417 :   d = mulii(d,delta); avma = (pari_sp)z;
     882     1672417 :   gel(z,1) = icopy(n);
     883     1672417 :   gel(z,2) = icopy(d); return z;
     884             : }
     885             : 
     886             : /* assume x2, y2 are t_POLs in the same variable */
     887             : static GEN
     888     3037311 : add_rfrac(GEN x, GEN y)
     889             : {
     890     3037311 :   pari_sp av = avma;
     891     3037311 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     892     3037311 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     893             : 
     894     3037311 :   delta = RgX_gcd(x2,y2);
     895     3037311 :   if (!degpol(delta))
     896             :   {
     897         658 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     898         658 :     d = RgX_mul(x2, y2);
     899         658 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     900             :   }
     901     3036653 :   x2 = RgX_div(x2,delta);
     902     3036653 :   y2 = RgX_div(y2,delta);
     903     3036653 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     904     3036653 :   if (!signe(n))
     905             :   {
     906      721054 :     n = simplify_shallow(n);
     907      721054 :     if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
     908           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     909             :   }
     910     2315599 :   if (degpol(n) == 0)
     911     1149721 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     912     1165878 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     913     1165878 :   if (isexactzero(r))
     914             :   {
     915             :     GEN z;
     916      228045 :     d = RgX_mul(x2, y2);
     917             :     /* "constant" denominator ? */
     918      228045 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     919      228045 :     return gerepileupto(av, z);
     920             :   }
     921      937833 :   r = RgX_gcd(delta, r);
     922      937833 :   if (degpol(r))
     923             :   {
     924      160705 :     n = RgX_div(n, r);
     925      160705 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     926             :   }
     927             :   else
     928      777128 :     d = RgX_mul(gel(x,2), y2);
     929      937833 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     930             : }
     931             : 
     932             : GEN
     933  1975501082 : gadd(GEN x, GEN y)
     934             : {
     935  1975501082 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     936             :   pari_sp av, tetpil;
     937             :   GEN z, p1;
     938             : 
     939  1975501082 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     940             :   {
     941  1110886899 :     case t_INT: return addii(x,y);
     942   516587132 :     case t_REAL: return addrr(x,y);
     943     1518154 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     944     1518154 :       z = cgetg(3,t_INTMOD);
     945     1518141 :       if (X==Y || equalii(X,Y))
     946     1518122 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     947          14 :       gel(z,1) = gcdii(X,Y);
     948          14 :       warn_coercion(X,Y,gel(z,1));
     949          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     950          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     951             :     }
     952    12157073 :     case t_FRAC: return addsub_frac(x,y,addii);
     953    69880313 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     954    69880313 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     955    69880313 :       if (isintzero(gel(z,2)))
     956             :       {
     957      324177 :         avma = (pari_sp)(z+3);
     958      324177 :         return gadd(gel(x,1),gel(y,1));
     959             :       }
     960    69556136 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     961    69556136 :       return z;
     962             :     case t_PADIC:
     963      949627 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     964      949627 :       return addsub_pp(x,y, addii);
     965         476 :     case t_QUAD: z = cgetg(4,t_QUAD);
     966         476 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     967         476 :       gel(z,1) = ZX_copy(gel(x,1));
     968         476 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     969         476 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     970             :     case t_POLMOD:
     971      730044 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     972      728329 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     973        1715 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     974     9194356 :     case t_FFELT: return FF_add(x,y);
     975             :     case t_POL:
     976    17120640 :       vx = varn(x);
     977    17120640 :       vy = varn(y);
     978    17120640 :       if (vx != vy) {
     979      811768 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     980       20851 :         else                     return RgX_Rg_add(y, x);
     981             :       }
     982    16308872 :       return RgX_add(x, y);
     983             :     case t_SER:
     984     1765084 :       vx = varn(x);
     985     1765084 :       vy = varn(y);
     986     1765084 :       if (vx != vy) {
     987          14 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
     988          14 :         else                     return add_ser_scal(y, x);
     989             :       }
     990     1765070 :       return ser_add(x, y);
     991             :     case t_RFRAC:
     992     4329198 :       vx = varn(gel(x,2));
     993     4329198 :       vy = varn(gel(y,2));
     994     4329198 :       if (vx != vy) {
     995     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
     996      538397 :         else                     return add_rfrac_scal(y, x);
     997             :       }
     998     3037311 :       return add_rfrac(x,y);
     999             :     case t_VEC:
    1000      589684 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1001      589684 :       return RgV_add(x,y);
    1002             :     case t_COL:
    1003      857714 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1004      857714 :       return RgC_add(x,y);
    1005             :     case t_MAT:
    1006      671902 :       lx = lg(x);
    1007      671902 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1008      671902 :       if (lx == 1) return cgetg(1, t_MAT);
    1009      671902 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1010      671895 :       return RgM_add(x,y);
    1011             : 
    1012           0 :     default: pari_err_TYPE2("+",x,y);
    1013             :   }
    1014             :   /* tx != ty */
    1015   228262794 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1016             : 
    1017   228262794 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1018             :   {
    1019             :     case t_INT:
    1020   189461810 :       switch(ty)
    1021             :       {
    1022   112259880 :         case t_REAL: return addir(x,y);
    1023             :         case t_INTMOD:
    1024      714261 :           z = cgetg(3, t_INTMOD);
    1025      714261 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1026    20328311 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1027    20328311 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1028    20328311 :           gel(z,2) = icopy(gel(y,2)); return z;
    1029    54170115 :         case t_COMPLEX: return addRc(x, y);
    1030             :         case t_PADIC:
    1031      330712 :           if (!signe(x)) return gcopy(y);
    1032      299741 :           return addQp(x,y);
    1033         714 :         case t_QUAD: return addRq(x, y);
    1034     1657820 :         case t_FFELT: return FF_Z_add(y,x);
    1035             :       }
    1036             : 
    1037             :     case t_REAL:
    1038    18604426 :       switch(ty)
    1039             :       {
    1040             :         case t_FRAC:
    1041      713774 :           if (!signe(gel(y,1))) return rcopy(x);
    1042      713774 :           if (!signe(x))
    1043             :           {
    1044        4753 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1045        4753 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1046             :           }
    1047      709021 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1048      709021 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1049    17890631 :         case t_COMPLEX: return addRc(x, y);
    1050          21 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1051             : 
    1052           0 :         default: pari_err_TYPE2("+",x,y);
    1053             :       }
    1054             : 
    1055             :     case t_INTMOD:
    1056        3633 :       switch(ty)
    1057             :       {
    1058        3500 :         case t_FRAC: { GEN X = gel(x,1);
    1059        3500 :           z = cgetg(3, t_INTMOD);
    1060        3500 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1061        3500 :           return add_intmod_same(z, X, p1, gel(x,2));
    1062             :         }
    1063             :         case t_FFELT:
    1064          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1065           0 :             pari_err_OP("+",x,y);
    1066          14 :           return FF_Z_add(y,gel(x,2));
    1067          91 :         case t_COMPLEX: return addRc(x, y);
    1068           0 :         case t_PADIC: { GEN X = gel(x,1);
    1069           0 :           z = cgetg(3, t_INTMOD);
    1070           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1071             :         }
    1072          28 :         case t_QUAD: return addRq(x, y);
    1073             :       }
    1074             : 
    1075             :     case t_FRAC:
    1076      307811 :       switch (ty)
    1077             :       {
    1078      235697 :         case t_COMPLEX: return addRc(x, y);
    1079             :         case t_PADIC:
    1080       70756 :           if (!signe(gel(x,1))) return gcopy(y);
    1081       70756 :           return addQp(x,y);
    1082          14 :         case t_QUAD: return addRq(x, y);
    1083        1344 :         case t_FFELT: return FF_Q_add(y, x);
    1084             :       }
    1085             : 
    1086             :     case t_FFELT:
    1087           0 :       pari_err_TYPE2("+",x,y);
    1088             : 
    1089             :     case t_COMPLEX:
    1090          21 :       switch(ty)
    1091             :       {
    1092             :         case t_PADIC:
    1093          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1094             :         case t_QUAD:
    1095           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1096           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1097             :       }
    1098             : 
    1099             :     case t_PADIC: /* ty == t_QUAD */
    1100           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1101             :   }
    1102             :   /* tx < ty, !is_const_t(y) */
    1103    19885092 :   switch(ty)
    1104             :   {
    1105             :     case t_MAT:
    1106       68075 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1107       68075 :       if (isrationalzero(x)) return gcopy(y);
    1108        5474 :       return RgM_Rg_add(y, x);
    1109             :     case t_COL:
    1110      331992 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1111      331992 :       return RgC_Rg_add(y, x);
    1112             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1113      295953 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1114             :   }
    1115    19189072 :   if (is_scalar_t(tx))  {
    1116    16101962 :     if (tx == t_POLMOD)
    1117             :     {
    1118       80479 :       vx = varn(gel(x,1));
    1119       80479 :       vy = gvar(y);
    1120       80479 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1121             :       else
    1122       47250 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1123       43813 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1124             :     }
    1125    16021483 :     return add_scal(y, x, ty);
    1126             :   }
    1127             :   /* x and y are not scalars, ty != t_MAT */
    1128     3087110 :   vx = gvar(x);
    1129     3087110 :   vy = gvar(y);
    1130     3087110 :   if (vx != vy) { /* x or y is treated as a scalar */
    1131       24488 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1132       24481 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1133       24481 :                                 : add_scal(y, x, ty);
    1134             :   }
    1135             :   /* vx = vy */
    1136     3062622 :   switch(tx)
    1137             :   {
    1138             :     case t_POL:
    1139     3062216 :       switch (ty)
    1140             :       {
    1141             :         case t_SER:
    1142        4823 :           if (lg(x) == 2) return gcopy(y);
    1143        4809 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1144        4809 :           i = lg(y) + valp(y) - i;
    1145        4809 :           if (i < 3) return gcopy(y);
    1146        4802 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1147        4802 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1148        4802 :           return y;
    1149             : 
    1150     3057393 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1151             :       }
    1152           0 :       break;
    1153             : 
    1154             :     case t_SER:
    1155         406 :       if (ty == t_RFRAC)
    1156             :       {
    1157             :         GEN n, d;
    1158             :         long vn, vd;
    1159         406 :         av = avma;
    1160         406 :         n = gel(y,1); vn = gval(n, vy);
    1161         406 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1162             : 
    1163         406 :         l = lg(x) + valp(x) - (vn - vd);
    1164         406 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1165             : 
    1166             :         /* take advantage of y = t^n ! */
    1167         406 :         if (degpol(d))
    1168          70 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1169             :         else {
    1170         336 :           y = gdiv(n, gel(d,2));
    1171         336 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1172             :         }
    1173         406 :         setvalp(y, valp(y) - vd);
    1174         406 :         return gerepileupto(av, gadd(y, x));
    1175             :       }
    1176           0 :       break;
    1177             :   }
    1178           0 :   pari_err_TYPE2("+",x,y);
    1179             :   return NULL; /* LCOV_EXCL_LINE */
    1180             : }
    1181             : 
    1182             : GEN
    1183    12825272 : gaddsg(long x, GEN y)
    1184             : {
    1185    12825272 :   long ty = typ(y);
    1186             :   GEN z;
    1187             : 
    1188    12825272 :   switch(ty)
    1189             :   {
    1190     3346000 :     case t_INT:  return addsi(x,y);
    1191     4244021 :     case t_REAL: return addsr(x,y);
    1192             :     case t_INTMOD:
    1193          14 :       z = cgetg(3, t_INTMOD);
    1194          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1195      745038 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1196      745038 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1197      745038 :       gel(z,2) = icopy(gel(y,2)); return z;
    1198             :     case t_COMPLEX:
    1199     3901854 :       z = cgetg(3, t_COMPLEX);
    1200     3901854 :       gel(z,1) = gaddsg(x, gel(y,1));
    1201     3901854 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1202             : 
    1203      588345 :     default: return gadd(stoi(x), y);
    1204             :   }
    1205             : }
    1206             : 
    1207             : GEN
    1208      382686 : gsubsg(long x, GEN y)
    1209             : {
    1210             :   GEN z, a, b;
    1211             :   pari_sp av;
    1212             : 
    1213      382686 :   switch(typ(y))
    1214             :   {
    1215       99600 :     case t_INT:  return subsi(x,y);
    1216      125555 :     case t_REAL: return subsr(x,y);
    1217             :     case t_INTMOD:
    1218          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1219          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1220       34502 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1221       34502 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1222       34502 :       gel(z,2) = icopy(gel(y,2)); return z;
    1223             :     case t_COMPLEX:
    1224       94861 :       z = cgetg(3, t_COMPLEX);
    1225       94861 :       gel(z,1) = gsubsg(x, gel(y,1));
    1226       94861 :       gel(z,2) = gneg(gel(y,2)); return z;
    1227             :   }
    1228       28119 :   av = avma;
    1229       28119 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1230             : }
    1231             : 
    1232             : /********************************************************************/
    1233             : /**                                                                **/
    1234             : /**                          SUBTRACTION                           **/
    1235             : /**                                                                **/
    1236             : /********************************************************************/
    1237             : 
    1238             : GEN
    1239  1434694516 : gsub(GEN x, GEN y)
    1240             : {
    1241  1434694516 :   long tx = typ(x), ty = typ(y);
    1242             :   pari_sp av;
    1243             :   GEN z;
    1244  1434694516 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1245             :   {
    1246  1187126313 :     case t_INT: return subii(x,y);
    1247   187707507 :     case t_REAL: return subrr(x,y);
    1248     1102018 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1249     1102018 :       z = cgetg(3,t_INTMOD);
    1250     1102018 :       if (X==Y || equalii(X,Y))
    1251     1102001 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1252          14 :       gel(z,1) = gcdii(X,Y);
    1253          14 :       warn_coercion(X,Y,gel(z,1));
    1254          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1255          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1256             :     }
    1257     3751908 :     case t_FRAC: return addsub_frac(x,y, subii);
    1258    25620883 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1259    25620883 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1260    25620883 :       if (isintzero(gel(z,2)))
    1261             :       {
    1262       80779 :         avma = (pari_sp)(z+3);
    1263       80779 :         return gsub(gel(x,1),gel(y,1));
    1264             :       }
    1265    25540104 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1266    25540104 :       return z;
    1267             :     case t_PADIC:
    1268      228855 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1269      228855 :       return addsub_pp(x,y, subii);
    1270         595 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1271         595 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1272         595 :       gel(z,1) = ZX_copy(gel(x,1));
    1273         595 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1274         595 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1275             :     case t_POLMOD:
    1276      911821 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1277      911786 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1278          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1279      282784 :     case t_FFELT: return FF_sub(x,y);
    1280             :     case t_POL: {
    1281     7310047 :       long vx = varn(x);
    1282     7310047 :       long vy = varn(y);
    1283     7310047 :       if (vx != vy) {
    1284       14953 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1285        1981 :         else                     return Rg_RgX_sub(x, y);
    1286             :       }
    1287     7295094 :       return RgX_sub(x, y);
    1288             :     }
    1289             :     case t_VEC:
    1290      271021 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1291      271021 :       return RgV_sub(x,y);
    1292             :     case t_COL:
    1293     3245340 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1294     3245340 :       return RgC_sub(x,y);
    1295             :     case t_MAT: {
    1296       11992 :       long lx = lg(x);
    1297       11992 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1298       11992 :       if (lx == 1) return cgetg(1, t_MAT);
    1299       11992 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1300       11992 :       return RgM_sub(x,y);
    1301             :     }
    1302     2095570 :     case t_RFRAC: case t_SER: break;
    1303             : 
    1304           0 :     default: pari_err_TYPE2("+",x,y);
    1305             :   }
    1306    17123432 :   av = avma;
    1307    17123432 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1308             : }
    1309             : 
    1310             : /********************************************************************/
    1311             : /**                                                                **/
    1312             : /**                        MULTIPLICATION                          **/
    1313             : /**                                                                **/
    1314             : /********************************************************************/
    1315             : static GEN
    1316      515039 : mul_ser_scal(GEN y, GEN x) {
    1317             :   long ly, i;
    1318             :   GEN z;
    1319      515039 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1320      511791 :   if (ser_isexactzero(y))
    1321             :   {
    1322         441 :     if (lg(y) == 2) return gcopy(y);
    1323           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1324             :   }
    1325      511350 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1326      511350 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1327      511350 :   return normalize(z);
    1328             : }
    1329             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1330             :  * [n/d a valid RFRAC]  */
    1331             : static GEN
    1332    10442554 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1333             : {
    1334    10442554 :   pari_sp av = avma;
    1335             :   GEN z;
    1336             : 
    1337    10442554 :   switch(typ(x))
    1338             :   {
    1339             :     case t_PADIC:
    1340          14 :       n = gmul(n, x);
    1341          14 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1342          14 :       return gerepileupto(av, gdiv(n,d));
    1343             : 
    1344             :     case t_INTMOD: case t_POLMOD:
    1345          49 :       n = gmul(n, x);
    1346          49 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1347          49 :       return gerepileupto(av, gdiv(n,d));
    1348             :   }
    1349    10442491 :   z = gred_rfrac2(x, d);
    1350    10442491 :   n = simplify_shallow(n);
    1351    10442491 :   if (typ(z) == t_RFRAC)
    1352             :   {
    1353     7927712 :     n = gmul(gel(z,1), n);
    1354     7927712 :     d = gel(z,2);
    1355     7927712 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1356           0 :       z = RgX_Rg_div(n, d);
    1357             :     else
    1358     7927712 :       z = gred_rfrac_simple(n, d);
    1359             :   }
    1360             :   else
    1361     2514779 :     z = gmul(z, n);
    1362    10442491 :   return gerepileupto(av, z);
    1363             : }
    1364             : static GEN
    1365    62223124 : mul_scal(GEN y, GEN x, long ty)
    1366             : {
    1367    62223124 :   switch(ty)
    1368             :   {
    1369             :     case t_POL:
    1370    53263783 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1371    52243699 :       return RgX_Rg_mul(y, x);
    1372      380313 :     case t_SER: return mul_ser_scal(y, x);
    1373     8578986 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1374             :     case t_QFI: case t_QFR:
    1375          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1376             :   }
    1377          28 :   pari_err_TYPE2("*",x,y);
    1378             :   return NULL; /* LCOV_EXCL_LINE */
    1379             : }
    1380             : 
    1381             : static GEN
    1382      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1383             : {
    1384      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1385      160830 :   long vx = gvar(X), vy = varn(y2);
    1386      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1387        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1388             : }
    1389             : /* (x1/x2) * (y1/y2) */
    1390             : static GEN
    1391     7903558 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1392             : {
    1393             :   GEN z, X, Y;
    1394     7903558 :   pari_sp av = avma;
    1395             : 
    1396     7903558 :   X = gred_rfrac2(x1, y2);
    1397     7903558 :   Y = gred_rfrac2(y1, x2);
    1398     7903558 :   if (typ(X) == t_RFRAC)
    1399             :   {
    1400     6624630 :     if (typ(Y) == t_RFRAC) {
    1401     6559054 :       x1 = gel(X,1);
    1402     6559054 :       x2 = gel(X,2);
    1403     6559054 :       y1 = gel(Y,1);
    1404     6559054 :       y2 = gel(Y,2);
    1405     6559054 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1406             :     } else
    1407       65576 :       z = mul_gen_rfrac(Y, X);
    1408             :   }
    1409     1278928 :   else if (typ(Y) == t_RFRAC)
    1410       95254 :     z = mul_gen_rfrac(X, Y);
    1411             :   else
    1412     1183674 :     z = gmul(X, Y);
    1413     7903558 :   return gerepileupto(av, z);
    1414             : }
    1415             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1416             : static GEN
    1417      269855 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1418             : {
    1419      269855 :   pari_sp av = avma;
    1420      269855 :   GEN X = gred_rfrac2(x1, y2);
    1421      269855 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1422             :   {
    1423      263303 :     x2 = RgX_mul(gel(X,2), x2);
    1424      263303 :     x1 = gel(X,1);
    1425             :   }
    1426             :   else
    1427        6552 :     x1 = X;
    1428      269855 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1429             : }
    1430             : 
    1431             : /* Mod(y, Y) * x,  assuming x scalar */
    1432             : static GEN
    1433      903177 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1434      903177 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1435             : 
    1436             : /* cf mulqq */
    1437             : static GEN
    1438      675612 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1439             : {
    1440      675612 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1441      675612 :   pari_sp tetpil, av = avma;
    1442      675612 :   T[1] = x[1];
    1443      675612 :   p2 = gmul(gel(x,2), gel(y,2));
    1444      675612 :   p3 = gmul(gel(x,3), gel(y,3));
    1445      675612 :   p1 = gmul(gneg_i(c),p3);
    1446             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1447      675612 :   if (typ(b) == t_INT)
    1448             :   {
    1449      675598 :     if (signe(b))
    1450             :     {
    1451      395864 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1452      395864 :       if (is_pm1(b))
    1453             :       {
    1454      395206 :         if (signe(b) > 0) p3 = gneg(p3);
    1455             :       }
    1456             :       else
    1457         658 :         p3 = gmul(negi(b), p3);
    1458             :     }
    1459             :     else
    1460             :     {
    1461      279734 :       p3 = gmul(gel(x,2),gel(y,3));
    1462      279734 :       p4 = gmul(gel(x,3),gel(y,2));
    1463             :     }
    1464             :   }
    1465             :   else
    1466             :   {
    1467          14 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1468          14 :     p3 = gmul(gneg_i(b), p3);
    1469             :   }
    1470      675612 :   tetpil = avma;
    1471      675612 :   gel(T,2) = gadd(p2, p1);
    1472      675612 :   gel(T,3) = gadd(p4, p3);
    1473      675612 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1474      675612 :   return normalizepol_lg(T,4);
    1475             : }
    1476             : /* Mod(x,T) * Mod(y,T) */
    1477             : static GEN
    1478     1796721 : mul_polmod_same(GEN T, GEN x, GEN y)
    1479             : {
    1480     1796721 :   GEN z = cgetg(3,t_POLMOD), a;
    1481     1796721 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1482     1796721 :   gel(z,1) = RgX_copy(T);
    1483             :   /* x * y mod T optimised */
    1484     1796721 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1485     1273510 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1486      711087 :     a = gmul(x, y);
    1487             :   else
    1488             :   {
    1489     1085634 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1490      669648 :       a = quad_polmod_mul(T, x, y);
    1491             :     else
    1492      415986 :       a = RgXQ_mul(x, y, gel(z,1));
    1493             :   }
    1494     1796721 :   gel(z,2) = a; return z;
    1495             : }
    1496             : static GEN
    1497       34048 : sqr_polmod(GEN T, GEN x)
    1498             : {
    1499       34048 :   GEN a, z = cgetg(3,t_POLMOD);
    1500       34048 :   gel(z,1) = RgX_copy(T);
    1501       34048 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1502        7707 :     a = gsqr(x);
    1503             :   else
    1504             :   {
    1505       26341 :     pari_sp av = avma;
    1506       26341 :     a = RgXQ_sqr(x, gel(z,1));
    1507       26341 :     a = gerepileupto(av, a);
    1508             :   }
    1509       34048 :   gel(z,2) = a; return z;
    1510             : }
    1511             : /* Mod(x,X) * Mod(y,Y) */
    1512             : static GEN
    1513      141148 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1514             : {
    1515      141148 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1516      141148 :   long vx = varn(X), vy = varn(Y);
    1517      141148 :   GEN z = cgetg(3,t_POLMOD);
    1518             : 
    1519      141148 :   if (vx==vy) {
    1520             :     pari_sp av;
    1521          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1522          14 :     warn_coercion(X,Y,gel(z,1));
    1523          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1524          14 :     return z;
    1525             :   }
    1526      141134 :   if (varncmp(vx, vy) < 0)
    1527      133196 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1528             :   else
    1529        7938 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1530      141134 :   gel(z,2) = gmul(x, y); return z;
    1531             : }
    1532             : 
    1533             : #if 0 /* used by 3M only */
    1534             : /* set z = x+y and return 1 if x,y have the same sign
    1535             :  * set z = x-y and return 0 otherwise */
    1536             : static int
    1537             : did_add(GEN x, GEN y, GEN *z)
    1538             : {
    1539             :   long tx = typ(x), ty = typ(y);
    1540             :   if (tx == ty) switch(tx)
    1541             :   {
    1542             :     case t_INT: *z = addii(x,y); return 1;
    1543             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1544             :     case t_REAL:
    1545             :       if (signe(x) == -signe(y))
    1546             :       { *z = subrr(x,y); return 0; }
    1547             :       else
    1548             :       { *z = addrr(x,y); return 1; }
    1549             :   }
    1550             :   if (tx == t_REAL) switch(ty)
    1551             :   {
    1552             :     case t_INT:
    1553             :       if (signe(x) == -signe(y))
    1554             :       { *z = subri(x,y); return 0; }
    1555             :       else
    1556             :       { *z = addri(x,y); return 1; }
    1557             :     case t_FRAC:
    1558             :       if (signe(x) == -signe(gel(y,1)))
    1559             :       { *z = gsub(x,y); return 0; }
    1560             :       else
    1561             :       { *z = gadd(x,y); return 1; }
    1562             :   }
    1563             :   else if (ty == t_REAL) switch(tx)
    1564             :   {
    1565             :     case t_INT:
    1566             :       if (signe(x) == -signe(y))
    1567             :       { *z = subir(x,y); return 0; }
    1568             :       else
    1569             :       { *z = addir(x,y); return 1; }
    1570             :     case t_FRAC:
    1571             :       if (signe(gel(x,1)) == -signe(y))
    1572             :       { *z = gsub(x,y); return 0; }
    1573             :       else
    1574             :       { *z = gadd(x,y); return 1; }
    1575             :   }
    1576             :   *z = gadd(x,y); return 1;
    1577             : }
    1578             : #endif
    1579             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1580             : static GEN
    1581     4238064 : mulcIR(GEN x, GEN y)
    1582             : {
    1583     4238064 :   GEN z = cgetg(3,t_COMPLEX);
    1584     4238064 :   pari_sp av = avma;
    1585     4238064 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1586     4238064 :   gel(z,2) = gmul(y, gel(x,1));
    1587     4238064 :   return z;
    1588             : 
    1589             : }
    1590             : /* x,y COMPLEX */
    1591             : static GEN
    1592    99675528 : mulcc(GEN x, GEN y)
    1593             : {
    1594    99675528 :   GEN xr = gel(x,1), xi = gel(x,2);
    1595    99675528 :   GEN yr = gel(y,1), yi = gel(y,2);
    1596             :   GEN p1, p2, p3, p4, z;
    1597             :   pari_sp tetpil, av;
    1598             : 
    1599    99675528 :   if (isintzero(xr))
    1600             :   {
    1601     4112823 :     if (isintzero(yr)) {
    1602      590372 :       av = avma;
    1603      590372 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1604             :     }
    1605     3522451 :     return mulcIR(y, xi);
    1606             :   }
    1607    95562705 :   if (isintzero(yr)) return mulcIR(x, yi);
    1608             : 
    1609    94847092 :   z = cgetg(3,t_COMPLEX); av = avma;
    1610             : #if 0
    1611             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1612             :    * e.g. xr + xi if exponents differ */
    1613             :   if (did_add(xr, xi, &p3))
    1614             :   {
    1615             :     if (did_add(yr, yi, &p4)) {
    1616             :     /* R = xr*yr - xi*yi
    1617             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1618             :       p1 = gmul(xr,yr);
    1619             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1620             :       p3 = gmul(p3, p4);
    1621             :       p4 = gsub(p2, p1);
    1622             :     } else {
    1623             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1624             :      * I = xr*yi + xi*yr */
    1625             :       p1 = gmul(p3,p4);
    1626             :       p3 = gmul(xr,yi);
    1627             :       p4 = gmul(xi,yr);
    1628             :       p2 = gsub(p3, p4);
    1629             :     }
    1630             :   } else {
    1631             :     if (did_add(yr, yi, &p4)) {
    1632             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1633             :       * I = xr*yi +xi*yr */
    1634             :       p1 = gmul(p3,p4);
    1635             :       p3 = gmul(xr,yi);
    1636             :       p4 = gmul(xi,yr);
    1637             :       p2 = gsub(p4, p3);
    1638             :     } else {
    1639             :     /* R = xr*yr - xi*yi
    1640             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1641             :       p3 = gneg( gmul(p3, p4) );
    1642             :       p1 = gmul(xr,yr);
    1643             :       p2 = gmul(xi,yi);
    1644             :       p4 = gadd(p1, p2);
    1645             : 
    1646             :       p2 = gneg(p2);
    1647             :     }
    1648             :   }
    1649             :   tetpil = avma;
    1650             :   gel(z,1) = gadd(p1,p2);
    1651             :   gel(z,2) = gadd(p3,p4);
    1652             : #else
    1653    94847092 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1654             :   { /* 3M formula */
    1655      198961 :     p3 = addii(xr,xi);
    1656      198961 :     p4 = addii(yr,yi);
    1657      198961 :     p1 = mulii(xr,yr);
    1658      198961 :     p2 = mulii(xi,yi);
    1659      198961 :     p3 = mulii(p3,p4);
    1660      198961 :     p4 = addii(p2,p1);
    1661      198961 :     tetpil = avma;
    1662      198961 :     gel(z,1) = subii(p1,p2);
    1663      198961 :     gel(z,2) = subii(p3,p4);
    1664      284956 :     if (!signe(gel(z,2)))
    1665      112966 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1666             :   }
    1667             :   else
    1668             :   { /* naive 4M formula: avoid all loss of accuracy */
    1669    94648131 :     p1 = gmul(xr,yr);
    1670    94648131 :     p2 = gmul(xi,yi);
    1671    94648131 :     p3 = gmul(xr,yi);
    1672    94648131 :     p4 = gmul(xi,yr);
    1673    94648131 :     tetpil = avma;
    1674    94648131 :     gel(z,1) = gsub(p1,p2);
    1675    94648131 :     gel(z,2) = gadd(p3,p4);
    1676    94648131 :     if (isintzero(gel(z,2)))
    1677             :     {
    1678       13950 :       cgiv(gel(z,2));
    1679       13950 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1680             :     }
    1681             :   }
    1682             : #endif
    1683             : 
    1684    94720176 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1685             : }
    1686             : /* x,y PADIC */
    1687             : static GEN
    1688     1774789 : mulpp(GEN x, GEN y) {
    1689     1774789 :   long l = valp(x) + valp(y);
    1690             :   pari_sp av;
    1691             :   GEN z, t;
    1692     1774789 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1693     1774798 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1694     1614112 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1695             : 
    1696     1549247 :   t = (precp(x) > precp(y))? y: x;
    1697     1549247 :   z = cgetp(t); setvalp(z,l); av = avma;
    1698     1549238 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1699     1549230 :   set_avma(av); return z;
    1700             : }
    1701             : /* x,y QUAD */
    1702             : static GEN
    1703        1316 : mulqq(GEN x, GEN y) {
    1704        1316 :   GEN z = cgetg(4,t_QUAD);
    1705        1316 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1706             :   pari_sp av, tetpil;
    1707        1316 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1708             : 
    1709        1316 :   gel(z,1) = ZX_copy(P); av = avma;
    1710        1316 :   p2 = gmul(gel(x,2),gel(y,2));
    1711        1316 :   p3 = gmul(gel(x,3),gel(y,3));
    1712        1316 :   p1 = gmul(gneg_i(c),p3);
    1713             : 
    1714        1316 :   if (signe(b))
    1715        1204 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1716             :   else
    1717             :   {
    1718         112 :     p3 = gmul(gel(x,2),gel(y,3));
    1719         112 :     p4 = gmul(gel(x,3),gel(y,2));
    1720             :   }
    1721        1316 :   tetpil = avma;
    1722        1316 :   gel(z,2) = gadd(p2,p1);
    1723        1316 :   gel(z,3) = gadd(p4,p3);
    1724        1316 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1725             : }
    1726             : 
    1727             : GEN
    1728     2277488 : mulcxI(GEN x)
    1729             : {
    1730             :   GEN z;
    1731     2277488 :   switch(typ(x))
    1732             :   {
    1733             :     case t_INT: case t_REAL: case t_FRAC:
    1734      281875 :       return mkcomplex(gen_0, x);
    1735             :     case t_COMPLEX:
    1736     1995557 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1737     1993240 :       z = cgetg(3,t_COMPLEX);
    1738     1993240 :       gel(z,1) = gneg(gel(x,2));
    1739     1993240 :       gel(z,2) = gel(x,1); return z;
    1740             :     default:
    1741          56 :       return gmul(gen_I(), x);
    1742             :   }
    1743             : }
    1744             : GEN
    1745       49133 : mulcxmI(GEN x)
    1746             : {
    1747             :   GEN z;
    1748       49133 :   switch(typ(x))
    1749             :   {
    1750             :     case t_INT: case t_REAL: case t_FRAC:
    1751        1806 :       return mkcomplex(gen_0, gneg(x));
    1752             :     case t_COMPLEX:
    1753       47264 :       if (isintzero(gel(x,1))) return gel(x,2);
    1754       46165 :       z = cgetg(3,t_COMPLEX);
    1755       46165 :       gel(z,1) = gel(x,2);
    1756       46165 :       gel(z,2) = gneg(gel(x,1)); return z;
    1757             :     default:
    1758          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1759             :   }
    1760             : }
    1761             : /* x * I^k */
    1762             : GEN
    1763        1260 : mulcxpowIs(GEN x, long k)
    1764             : {
    1765        1260 :   switch (k & 3)
    1766             :   {
    1767         266 :     case 1: return mulcxI(x);
    1768         392 :     case 2: return gneg(x);
    1769         336 :     case 3: return mulcxmI(x);
    1770             :   }
    1771         266 :   return x;
    1772             : }
    1773             : 
    1774             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1775             : static GEN
    1776     2232603 : fill_ser(GEN z, GEN y)
    1777             : {
    1778     2232603 :   long i, lx = lg(z), ly = lg(y);
    1779     2232603 :   if (ly >= lx) {
    1780     2215901 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1781             :   } else {
    1782       16702 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1783       16702 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1784             :   }
    1785     2232603 :   return normalize(z);
    1786             : }
    1787             : 
    1788             : GEN
    1789  3431374960 : gmul(GEN x, GEN y)
    1790             : {
    1791             :   long tx, ty, lx, ly, vx, vy, i, l;
    1792             :   pari_sp av, tetpil;
    1793             :   GEN z, p1, p2;
    1794             : 
    1795  3431374960 :   if (x == y) return gsqr(x);
    1796  2492981446 :   tx = typ(x); ty = typ(y);
    1797  2492981446 :   if (tx == ty) switch(tx)
    1798             :   {
    1799  1182187689 :     case t_INT: return mulii(x,y);
    1800   683814974 :     case t_REAL: return mulrr(x,y);
    1801     1722274 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1802     1722274 :       z = cgetg(3,t_INTMOD);
    1803     1722259 :       if (X==Y || equalii(X,Y))
    1804     1722230 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1805          14 :       gel(z,1) = gcdii(X,Y);
    1806          14 :       warn_coercion(X,Y,gel(z,1));
    1807          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1808          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1809             :     }
    1810             :     case t_FRAC:
    1811             :     {
    1812     7033188 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1813     7033188 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1814     7033188 :       z=cgetg(3,t_FRAC);
    1815     7033188 :       p1 = gcdii(x1, y2);
    1816     7033188 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1817     7033188 :       p1 = gcdii(x2, y1);
    1818     7033188 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1819     7033188 :       tetpil = avma;
    1820     7033188 :       gel(z,2) = mulii(x2,y2);
    1821     7033188 :       gel(z,1) = mulii(x1,y1);
    1822     7033188 :       fix_frac_if_int_GC(z,tetpil); return z;
    1823             :     }
    1824    97607484 :     case t_COMPLEX: return mulcc(x, y);
    1825     1602864 :     case t_PADIC: return mulpp(x, y);
    1826         994 :     case t_QUAD: return mulqq(x, y);
    1827     9577425 :     case t_FFELT: return FF_mul(x, y);
    1828             :     case t_POLMOD:
    1829     1570033 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1830     1428885 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1831      141148 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1832             :     case t_POL:
    1833    35784278 :       vx = varn(x);
    1834    35784278 :       vy = varn(y);
    1835    35784278 :       if (vx != vy) {
    1836     4761409 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1837     2049484 :         else                     return RgX_Rg_mul(y, x);
    1838             :       }
    1839    31022869 :       return RgX_mul(x, y);
    1840             : 
    1841             :     case t_SER: {
    1842     2236145 :       vx = varn(x);
    1843     2236145 :       vy = varn(y);
    1844     2236145 :       if (vx != vy) {
    1845        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1846        3675 :         else                     return mul_ser_scal(y, x);
    1847             :       }
    1848     2232470 :       lx = lg(x);
    1849     2232470 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1850     2232470 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1851     2232414 :       av = avma; z = cgetg(lx,t_SER);
    1852     2232414 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1853     2232414 :       x = ser2pol_i(x, lx);
    1854     2232414 :       y = ser2pol_i(y, lx);
    1855     2232414 :       y = RgXn_mul(x, y, lx-2);
    1856     2232414 :       return gerepilecopy(av, fill_ser(z,y));
    1857             :     }
    1858    16035591 :     case t_QFI: return qficomp(x,y);
    1859          14 :     case t_QFR: return qfrcomp(x,y);
    1860     6717049 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1861     1773528 :     case t_MAT: return RgM_mul(x, y);
    1862             : 
    1863             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1864        1533 :       z = cgetg_copy(x, &l);
    1865        1533 :       if (l != lg(y)) break;
    1866       16933 :       for (i=1; i<l; i++)
    1867             :       {
    1868       15400 :         long yi = y[i];
    1869       15400 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1870       15400 :         z[i] = x[yi];
    1871             :       }
    1872        1533 :       return z;
    1873             : 
    1874             : 
    1875             :     default:
    1876           0 :       pari_err_TYPE2("*",x,y);
    1877             :   }
    1878             :   /* tx != ty */
    1879   445316383 :   if (is_const_t(ty) && is_const_t(tx))  {
    1880   374195094 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1881   374195094 :     switch(tx) {
    1882             :     case t_INT:
    1883   340557994 :       switch(ty)
    1884             :       {
    1885   223873520 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1886             :         case t_INTMOD:
    1887      966759 :           z = cgetg(3, t_INTMOD);
    1888      966753 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1889             :         case t_FRAC:
    1890    38555253 :           if (!signe(x)) return gen_0;
    1891    22109696 :           z=cgetg(3,t_FRAC);
    1892    22109696 :           p1 = gcdii(x,gel(y,2));
    1893    22109696 :           if (equali1(p1))
    1894             :           {
    1895    10754350 :             avma = (pari_sp)z;
    1896    10754350 :             gel(z,2) = icopy(gel(y,2));
    1897    10754350 :             gel(z,1) = mulii(gel(y,1), x);
    1898             :           }
    1899             :           else
    1900             :           {
    1901    11355346 :             x = diviiexact(x,p1); tetpil = avma;
    1902    11355346 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1903    11355346 :             gel(z,1) = mulii(gel(y,1), x);
    1904    11355346 :             fix_frac_if_int_GC(z,tetpil);
    1905             :           }
    1906    22109696 :           return z;
    1907    75143059 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1908      408942 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1909        1344 :         case t_QUAD: return mulRq(x,y);
    1910     1609113 :         case t_FFELT: return FF_Z_mul(y,x);
    1911             :       }
    1912             : 
    1913             :     case t_REAL:
    1914    32765278 :       switch(ty)
    1915             :       {
    1916     4461493 :         case t_FRAC: return mulrfrac(x, y);
    1917    28303736 :         case t_COMPLEX: return mulRc(x, y);
    1918           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1919          42 :         default: pari_err_TYPE2("*",x,y);
    1920             :       }
    1921             : 
    1922             :     case t_INTMOD:
    1923         826 :       switch(ty)
    1924             :       {
    1925         371 :         case t_FRAC: { GEN X = gel(x,1);
    1926         371 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1927         371 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1928             :         }
    1929          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1930         301 :         case t_PADIC: { GEN X = gel(x,1);
    1931         301 :           z = cgetg(3, t_INTMOD);
    1932         301 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1933             :         }
    1934          42 :         case t_QUAD: return mulRq(x, y);
    1935             :         case t_FFELT:
    1936          63 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1937           0 :             pari_err_OP("*",x,y);
    1938          63 :           return FF_Z_mul(y,gel(x,2));
    1939             :       }
    1940             : 
    1941             :     case t_FRAC:
    1942      870928 :       switch(ty)
    1943             :       {
    1944      838476 :         case t_COMPLEX: return mulRc(x, y);
    1945       30163 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1946         308 :         case t_QUAD:  return mulRq(x, y);
    1947        1981 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1948             :       }
    1949             : 
    1950             :     case t_FFELT:
    1951          35 :       pari_err_TYPE2("*",x,y);
    1952             : 
    1953             :     case t_COMPLEX:
    1954          21 :       switch(ty)
    1955             :       {
    1956             :         case t_PADIC:
    1957          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1958             :         case t_QUAD:
    1959           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1960           7 :           return mulqf(y, x, lx);
    1961             :       }
    1962             : 
    1963             :     case t_PADIC: /* ty == t_QUAD */
    1964          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1965             :     }
    1966             :   }
    1967             : 
    1968    71123308 :   if (is_matvec_t(ty))
    1969             :   {
    1970     5420197 :     if (!is_matvec_t(tx))
    1971             :     {
    1972     5257309 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1973     5257309 :       z = cgetg_copy(y, &ly);
    1974     5257309 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1975     5257295 :       return z;
    1976             :     }
    1977      162887 :     switch(tx)
    1978             :     {
    1979             :       case t_VEC:
    1980       41123 :         switch(ty) {
    1981       14861 :           case t_COL: return RgV_RgC_mul(x,y);
    1982       26262 :           case t_MAT: return RgV_RgM_mul(x,y);
    1983             :         }
    1984           0 :         break;
    1985             :       case t_COL:
    1986        1687 :         switch(ty) {
    1987        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    1988           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    1989             :         }
    1990           0 :         break;
    1991             :       case t_MAT:
    1992      120077 :         switch(ty) {
    1993           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    1994      120077 :           case t_COL: return RgM_RgC_mul(x,y);
    1995             :         }
    1996             :     }
    1997             :   }
    1998    65703129 :   if (is_matvec_t(tx))
    1999             :   {
    2000      367613 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2001      367613 :     z = cgetg_copy(x, &lx);
    2002      367619 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2003      367612 :     return z;
    2004             :   }
    2005    65335514 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2006             :   /* tx < ty, !ismatvec(x and y) */
    2007             : 
    2008    65335514 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2009      865384 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2010    64470130 :   if (is_scalar_t(tx)) {
    2011    61058016 :     if (tx == t_POLMOD) {
    2012      635990 :       vx = varn(gel(x,1));
    2013      635990 :       vy = gvar(y);
    2014      635990 :       if (vx != vy) {
    2015      268525 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2016       37793 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2017             :       }
    2018             :       /* error if ty == t_SER */
    2019      367465 :       av = avma; y = gmod(y, gel(x,1));
    2020      367458 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2021             :     }
    2022    60422026 :     return mul_scal(y, x, ty);
    2023             :   }
    2024             : 
    2025             :   /* x and y are not scalars, nor matvec */
    2026     3412114 :   vx = gvar(x);
    2027     3412114 :   vy = gvar(y);
    2028     3412114 :   if (vx != vy) /* x or y is treated as a scalar */
    2029     2771899 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2030     2771899 :                                 : mul_scal(y, x, ty);
    2031             :   /* vx = vy */
    2032     1996838 :   switch(tx)
    2033             :   {
    2034             :     case t_POL:
    2035     1996810 :       switch (ty)
    2036             :       {
    2037             :         case t_SER:
    2038             :         {
    2039             :           long vn;
    2040      133270 :           if (lg(x) == 2) return zeropol(vx);
    2041      133270 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2042      133039 :           av = avma;
    2043      133039 :           vn = RgX_valrem(x, &x);
    2044             :           /* take advantage of x = t^n ! */
    2045      133039 :           if (degpol(x)) {
    2046        1988 :             p1 = RgX_to_ser(x,lg(y));
    2047        1988 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2048        1988 :             p2 = gmul(p1,y);
    2049        1988 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2050             :           } else {
    2051      131051 :             set_avma(av);
    2052      131051 :             p2 = mul_ser_scal(y, gel(x,2));
    2053             :           }
    2054      133039 :           setvalp(p2, valp(p2) + vn);
    2055      133039 :           return p2;
    2056             :         }
    2057             : 
    2058     1863540 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2059             :       }
    2060           0 :       break;
    2061             : 
    2062             :     case t_SER:
    2063          28 :       switch (ty)
    2064             :       {
    2065             :         case t_RFRAC:
    2066          28 :           av = avma;
    2067          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2068             :       }
    2069           0 :       break;
    2070             :   }
    2071           0 :   pari_err_TYPE2("*",x,y);
    2072             :   return NULL; /* LCOV_EXCL_LINE */
    2073             : }
    2074             : 
    2075             : /* return a non-normalized result */
    2076             : GEN
    2077       88658 : sqr_ser_part(GEN x, long l1, long l2)
    2078             : {
    2079             :   long i, j, l;
    2080             :   pari_sp av;
    2081             :   GEN Z, z, p1, p2;
    2082             :   long mi;
    2083       88658 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2084       88651 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2085       88651 :   Z = cgetg(l2-l1+3, t_SER);
    2086       88651 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2087       88651 :   z = Z + 2-l1;
    2088       88651 :   x += 2; mi = 0;
    2089      340690 :   for (i=0; i<l1; i++)
    2090             :   {
    2091      252039 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2092             :   }
    2093             : 
    2094      608058 :   for (i=l1; i<=l2; i++)
    2095             :   {
    2096      519407 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2097      519407 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2098     1063830 :     for (j=i-mi; j<=minss(l,mi); j++)
    2099      544423 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2100      519407 :     p1 = gshift(p1,1);
    2101      519407 :     if ((i&1) == 0 && p2[i>>1])
    2102       50416 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2103      519407 :     gel(z,i) = gerepileupto(av,p1);
    2104             :   }
    2105       88651 :   return Z;
    2106             : }
    2107             : 
    2108             : GEN
    2109  1024456925 : gsqr(GEN x)
    2110             : {
    2111             :   long i, lx;
    2112             :   pari_sp av, tetpil;
    2113             :   GEN z, p1, p2, p3, p4;
    2114             : 
    2115  1024456925 :   switch(typ(x))
    2116             :   {
    2117   952829006 :     case t_INT: return sqri(x);
    2118    24290050 :     case t_REAL: return sqrr(x);
    2119      141891 :     case t_INTMOD: { GEN X = gel(x,1);
    2120      141891 :       z = cgetg(3,t_INTMOD);
    2121      141891 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2122      141894 :       gel(z,1) = icopy(X); return z;
    2123             :     }
    2124     1251289 :     case t_FRAC: return sqrfrac(x);
    2125             : 
    2126             :     case t_COMPLEX:
    2127     2966626 :       if (isintzero(gel(x,1))) {
    2128       30389 :         av = avma;
    2129       30389 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2130             :       }
    2131     2936237 :       z = cgetg(3,t_COMPLEX); av = avma;
    2132     2936237 :       p1 = gadd(gel(x,1),gel(x,2));
    2133     2936237 :       p2 = gsub(gel(x,1), gel(x,2));
    2134     2936237 :       p3 = gmul(gel(x,1),gel(x,2));
    2135     2936237 :       tetpil = avma;
    2136     2936237 :       gel(z,1) = gmul(p1,p2);
    2137     2936237 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2138             : 
    2139             :     case t_PADIC:
    2140       28063 :       z = cgetg(5,t_PADIC);
    2141       28063 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2142       28063 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2143       28063 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2144       28063 :       gel(z,2) = icopy(gel(x,2));
    2145       28063 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2146       28063 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2147       28063 :       return z;
    2148             : 
    2149          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2150          28 :       p1 = gel(x,1);
    2151          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2152          28 :       p2 = gsqr(gel(x,2));
    2153          28 :       p3 = gsqr(gel(x,3));
    2154          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2155             : 
    2156          28 :       if (gequal0(gel(p1,3)))
    2157             :       {
    2158           7 :         tetpil = avma;
    2159           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2160           7 :         av = avma;
    2161           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2162           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2163             :       }
    2164             : 
    2165          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2166          21 :       tetpil = avma;
    2167          21 :       gel(z,2) = gadd(p2,p4);
    2168          21 :       gel(z,3) = gadd(p1,p3);
    2169          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2170             : 
    2171             :     case t_POLMOD:
    2172       34048 :       return sqr_polmod(gel(x,1), gel(x,2));
    2173             : 
    2174     2892139 :     case t_FFELT: return FF_sqr(x);
    2175             : 
    2176     1733288 :     case t_POL: return RgX_sqr(x);
    2177             : 
    2178             :     case t_SER:
    2179       24101 :       lx = lg(x);
    2180       24101 :       if (ser_isexactzero(x)) {
    2181           7 :         GEN z = gcopy(x);
    2182           7 :         setvalp(z, 2*valp(x));
    2183           7 :         return z;
    2184             :       }
    2185       24094 :       if (lx < 40)
    2186       23905 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2187             :       else
    2188             :       {
    2189         189 :         pari_sp av = avma;
    2190         189 :         GEN z = cgetg(lx,t_SER);
    2191         189 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2192         189 :         x = ser2pol_i(x,lx);
    2193         189 :         x = RgXn_sqr(x, lx-2);
    2194         189 :         return gerepilecopy(av, fill_ser(z,x));
    2195             :       }
    2196             : 
    2197           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2198           7 :       gel(z,1) = gsqr(gel(x,1));
    2199           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2200             : 
    2201         980 :     case t_MAT: return RgM_sqr(x);
    2202          14 :     case t_QFR: return qfrsqr(x);
    2203    38264779 :     case t_QFI: return qfisqr(x);
    2204             :     case t_VECSMALL:
    2205         616 :       z = cgetg_copy(x, &lx);
    2206        6720 :       for (i=1; i<lx; i++)
    2207             :       {
    2208        6104 :         long xi = x[i];
    2209        6104 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2210        6104 :         z[i] = x[xi];
    2211             :       }
    2212         616 :       return z;
    2213             :   }
    2214           0 :   pari_err_TYPE2("*",x,x);
    2215             :   return NULL; /* LCOV_EXCL_LINE */
    2216             : }
    2217             : 
    2218             : /********************************************************************/
    2219             : /**                                                                **/
    2220             : /**                           DIVISION                             **/
    2221             : /**                                                                **/
    2222             : /********************************************************************/
    2223             : static GEN
    2224      271022 : div_rfrac_scal(GEN x, GEN y)
    2225             : {
    2226      271022 :   pari_sp av = avma;
    2227      271022 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2228      271022 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2229             : }
    2230             : static GEN
    2231     1004062 : div_scal_rfrac(GEN x, GEN y)
    2232             : {
    2233     1004062 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2234     1004062 :   pari_sp av = avma;
    2235     1004062 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2236             :   {
    2237           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2238           0 :     y1 = gel(y1,2);
    2239             :   }
    2240     1004062 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2241             : }
    2242             : static GEN
    2243     1186509 : div_rfrac(GEN x, GEN y)
    2244     1186509 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2245             : 
    2246             : static GEN
    2247      622137 : div_ser_scal(GEN x, GEN y) {
    2248             :   long i, lx;
    2249             :   GEN z;
    2250      622137 :   if (ser_isexactzero(x))
    2251             :   {
    2252           7 :     if (lg(x) == 2) return gcopy(x);
    2253           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2254             :   }
    2255      622130 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2256      622130 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2257      622130 :   return normalize(z);
    2258             : }
    2259             : GEN
    2260        7525 : ser_normalize(GEN x)
    2261             : {
    2262        7525 :   long i, lx = lg(x);
    2263             :   GEN c, z;
    2264        7525 :   if (lx == 2) return x;
    2265        7525 :   c = gel(x,2); if (gequal1(c)) return x;
    2266        7448 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2267        7448 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2268        7448 :   return z;
    2269             : }
    2270             : 
    2271             : static GEN
    2272     1480981 : div_T_scal(GEN x, GEN y, long tx) {
    2273     1480981 :   switch(tx)
    2274             :   {
    2275      591686 :     case t_POL: return RgX_Rg_div(x, y);
    2276      622130 :     case t_SER: return div_ser_scal(x, y);
    2277      267165 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2278             :   }
    2279           0 :   pari_err_TYPE2("/",x,y);
    2280             :   return NULL; /* LCOV_EXCL_LINE */
    2281             : }
    2282             : 
    2283             : static GEN
    2284     9222307 : div_scal_pol(GEN x, GEN y) {
    2285     9222307 :   long ly = lg(y);
    2286             :   pari_sp av;
    2287     9222307 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2288     9142395 :   if (isrationalzero(x)) return zeropol(varn(y));
    2289     7091003 :   av = avma;
    2290     7091003 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2291             : }
    2292             : static GEN
    2293       32571 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2294             :   GEN z;
    2295             :   long ly, i;
    2296       32571 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2297       32571 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2298       32571 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2299       32571 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2300       32571 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2301       32571 :   y = gdiv(z,y); pari_free(z); return y;
    2302             : }
    2303             : static GEN
    2304    10209296 : div_scal_T(GEN x, GEN y, long ty) {
    2305    10209296 :   switch(ty)
    2306             :   {
    2307     9173419 :     case t_POL: return div_scal_pol(x, y);
    2308       32571 :     case t_SER: return div_scal_ser(x, y);
    2309     1003306 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2310             :   }
    2311           0 :   pari_err_TYPE2("/",x,y);
    2312             :   return NULL; /* LCOV_EXCL_LINE */
    2313             : }
    2314             : 
    2315             : /* assume tx = ty = t_SER, same variable vx */
    2316             : static GEN
    2317      876010 : div_ser(GEN x, GEN y, long vx)
    2318             : {
    2319      876010 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2320             :   GEN y_lead, p1, p2, z;
    2321             : 
    2322      876010 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2323      876003 :   if (ser_isexactzero(x))
    2324             :   {
    2325       59843 :     if (lx == 2) return zeroser(vx, l);
    2326         154 :     return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), l);
    2327             :   }
    2328      816160 :   y_lead = gel(y,2);
    2329      816160 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2330             :   {
    2331          14 :     GEN y0 = y;
    2332          14 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2333          35 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2334             :     {
    2335          28 :       y_lead = gel(y,2);
    2336          28 :       if (!gequal0(y_lead)) break;
    2337             :     }
    2338          14 :     if (ly <= 2) pari_err_INV("div_ser", y0);
    2339             :   }
    2340      816153 :   if (ly < lx) lx = ly;
    2341      816153 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2342     2970010 :   for (i=3; i<lx; i++)
    2343             :   {
    2344     2153857 :     p1 = gel(y,i);
    2345     2153857 :     if (isrationalzero(p1)) p1 = NULL;
    2346     2153857 :     gel(p2,i) = p1;
    2347             :   }
    2348      816153 :   z = cgetg(lx,t_SER);
    2349      816153 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2350      816153 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2351     2970010 :   for (i=3; i<lx; i++)
    2352             :   {
    2353     2153857 :     pari_sp av = avma;
    2354     2153857 :     p1 = gel(x,i);
    2355    12797502 :     for (j=2, l=i; j<i; j++, l--)
    2356    10643645 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2357     2153857 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2358             :   }
    2359      816153 :   return normalize(z);
    2360             : }
    2361             : /* x,y compatible PADIC */
    2362             : static GEN
    2363      224634 : divpp(GEN x, GEN y) {
    2364             :   pari_sp av;
    2365             :   long a, b;
    2366             :   GEN z, M;
    2367             : 
    2368      224634 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2369      224626 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2370      223968 :   a = precp(x);
    2371      223968 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2372      223968 :   z = cgetg(5, t_PADIC);
    2373      224000 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2374      223999 :   gel(z,2) = icopy(gel(x,2));
    2375      224002 :   gel(z,3) = icopy(M); av = avma;
    2376      224003 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2377      224008 :   return z;
    2378             : }
    2379             : static GEN
    2380       29469 : div_polmod_same(GEN T, GEN x, GEN y)
    2381             : {
    2382       29469 :   long v = varn(T);
    2383       29469 :   GEN a, z = cgetg(3, t_POLMOD);
    2384       29469 :   gel(z,1) = RgX_copy(T);
    2385       29469 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2386       10857 :     a = gdiv(x, y);
    2387       18612 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2388         280 :   {
    2389         280 :     pari_sp av = avma;
    2390         280 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2391             :   }
    2392       18332 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2393        5964 :   {
    2394        5964 :     pari_sp av = avma;
    2395        5964 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2396        5964 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2397        5964 :     a = gerepileupto(av, a);
    2398             :   }
    2399             :   else
    2400             :   {
    2401       12368 :     pari_sp av = avma;
    2402       12368 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2403       12368 :     a = gerepileupto(av, a);
    2404             :   }
    2405       29469 :   gel(z,2) = a; return z;
    2406             : }
    2407             : GEN
    2408   160673404 : gdiv(GEN x, GEN y)
    2409             : {
    2410   160673404 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2411             :   pari_sp av, tetpil;
    2412             :   GEN z, p1, p2;
    2413             : 
    2414   160673404 :   if (tx == ty) switch(tx)
    2415             :   {
    2416             :     case t_INT:
    2417    61830233 :       return Qdivii(x,y);
    2418             : 
    2419    47436125 :     case t_REAL: return divrr(x,y);
    2420       37780 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2421       37780 :       z = cgetg(3,t_INTMOD);
    2422       37781 :       if (X==Y || equalii(X,Y))
    2423       37752 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2424          28 :       gel(z,1) = gcdii(X,Y);
    2425          28 :       warn_coercion(X,Y,gel(z,1));
    2426          28 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2427          28 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2428             :     }
    2429             :     case t_FRAC: {
    2430      379857 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2431      379857 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2432      379857 :       z = cgetg(3, t_FRAC);
    2433      379857 :       p1 = gcdii(x1, y1);
    2434      379857 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2435      379857 :       p1 = gcdii(x2, y2);
    2436      379857 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2437      379857 :       tetpil = avma;
    2438      379857 :       gel(z,2) = mulii(x2,y1);
    2439      379857 :       gel(z,1) = mulii(x1,y2);
    2440      379857 :       normalize_frac(z);
    2441      379857 :       fix_frac_if_int_GC(z,tetpil);
    2442      379857 :       return z;
    2443             :     }
    2444             :     case t_COMPLEX:
    2445     2105837 :       if (isintzero(gel(y,1)))
    2446             :       {
    2447       37793 :         y = gel(y,2);
    2448       37793 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2449       20594 :         z = cgetg(3,t_COMPLEX);
    2450       20594 :         gel(z,1) = gdiv(gel(x,2), y);
    2451       20594 :         av = avma;
    2452       20594 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2453       20594 :         return z;
    2454             :       }
    2455     2068044 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2456     2068044 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2457             : 
    2458             :     case t_PADIC:
    2459      122137 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2460      122137 :       return divpp(x, y);
    2461             : 
    2462             :     case t_QUAD:
    2463         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2464         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2465         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2466             : 
    2467      234471 :     case t_FFELT: return FF_div(x,y);
    2468             : 
    2469             :     case t_POLMOD:
    2470       31100 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2471       29469 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2472             :       else {
    2473        1631 :         av = avma;
    2474        1631 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2475             :       }
    2476       31100 :       return z;
    2477             : 
    2478             :     case t_POL:
    2479    18512006 :       vx = varn(x);
    2480    18512006 :       vy = varn(y);
    2481    18512006 :       if (vx != vy) {
    2482       97545 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2483       48888 :                             else return div_scal_pol(x, y);
    2484             :       }
    2485    18414461 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2486    18414461 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2487    18227094 :       av = avma;
    2488    18227094 :       return gerepileupto(av, gred_rfrac2(x,y));
    2489             : 
    2490             :     case t_SER:
    2491      430453 :       vx = varn(x);
    2492      430453 :       vy = varn(y);
    2493      430453 :       if (vx != vy) {
    2494           7 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2495           0 :                             else return div_scal_ser(x, y);
    2496             :       }
    2497      430446 :       return div_ser(x, y, vx);
    2498             :     case t_RFRAC:
    2499     1190989 :       vx = varn(gel(x,2));
    2500     1190989 :       vy = varn(gel(y,2));
    2501     1190989 :       if (vx != vy) {
    2502        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2503         756 :                             else return div_scal_rfrac(x, y);
    2504             :       }
    2505     1186509 :       return div_rfrac(x,y);
    2506             : 
    2507           7 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2508           7 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2509             : 
    2510             :     case t_MAT:
    2511           7 :       av = avma; p1 = RgM_inv(y);
    2512           7 :       if (!p1) pari_err_INV("gdiv",y);
    2513           7 :       return gerepileupto(av, RgM_mul(x, p1));
    2514             : 
    2515           0 :     default: pari_err_TYPE2("/",x,y);
    2516             :   }
    2517             : 
    2518    28362073 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2519             :   {
    2520     1422536 :     long s = signe(x);
    2521     1422536 :     if (!s) {
    2522      679389 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2523      679389 :       switch (ty)
    2524             :       {
    2525      676204 :       default: return gen_0;
    2526             :       case t_INTMOD:
    2527          56 :         z = cgetg(3,t_INTMOD);
    2528          56 :         gel(z,1) = icopy(gel(y,1));
    2529          56 :         gel(z,2) = gen_0; return z;
    2530        3129 :       case t_FFELT: return FF_zero(y);
    2531             :       }
    2532             :     }
    2533      743147 :     if (is_pm1(x)) {
    2534      535245 :       if (s > 0) return ginv(y);
    2535       33236 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2536             :     }
    2537      207902 :     switch(ty)
    2538             :     {
    2539       57923 :       case t_REAL: return divir(x,y);
    2540             :       case t_INTMOD:
    2541          21 :         z = cgetg(3, t_INTMOD);
    2542          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2543             :       case t_FRAC:
    2544       90987 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2545       90987 :         if (equali1(p1))
    2546             :         {
    2547       88077 :           avma = (pari_sp)z;
    2548       88077 :           gel(z,2) = icopy(gel(y,1));
    2549       88077 :           gel(z,1) = mulii(gel(y,2), x);
    2550       88077 :           normalize_frac(z);
    2551       88077 :           fix_frac_if_int(z);
    2552             :         }
    2553             :         else
    2554             :         {
    2555        2910 :           x = diviiexact(x,p1); tetpil = avma;
    2556        2910 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2557        2910 :           gel(z,1) = mulii(gel(y,2), x);
    2558        2910 :           normalize_frac(z);
    2559        2910 :           fix_frac_if_int_GC(z,tetpil);
    2560             :         }
    2561       90987 :         return z;
    2562             : 
    2563         224 :       case t_FFELT: return Z_FF_div(x,y);
    2564       46504 :       case t_COMPLEX: return divRc(x,y);
    2565       12201 :       case t_PADIC: return divTp(x, y);
    2566             :       case t_QUAD:
    2567          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2568          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2569             :     }
    2570             :   }
    2571    26939537 :   if (gequal0(y))
    2572             :   {
    2573          42 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2574          21 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2575             :   }
    2576             : 
    2577    26939509 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2578             :   {
    2579             :     case t_REAL:
    2580     9830349 :       switch(ty)
    2581             :       {
    2582     8669935 :         case t_INT: return divri(x,y);
    2583             :         case t_FRAC:
    2584      718687 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2585      718687 :           return gerepileuptoleaf(av, z);
    2586      441699 :         case t_COMPLEX: return divRc(x, y);
    2587           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2588          21 :         default: pari_err_TYPE2("/",x,y);
    2589             :       }
    2590             : 
    2591             :     case t_INTMOD:
    2592       13482 :       switch(ty)
    2593             :       {
    2594             :         case t_INT:
    2595       13405 :           z = cgetg(3, t_INTMOD);
    2596       13405 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2597          28 :         case t_FRAC: { GEN X = gel(x,1);
    2598          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2599          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2600             :         }
    2601             :         case t_FFELT:
    2602           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2603           0 :             pari_err_OP("/",x,y);
    2604           7 :           return Z_FF_div(gel(x,2),y);
    2605             : 
    2606             :         case t_COMPLEX:
    2607           0 :           av = avma;
    2608           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2609             : 
    2610             :         case t_QUAD:
    2611           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2612           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2613             : 
    2614           7 :         case t_PADIC: { GEN X = gel(x,1);
    2615           7 :           z = cgetg(3, t_INTMOD);
    2616           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2617             :         }
    2618          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2619             :       }
    2620             : 
    2621             :     case t_FRAC:
    2622      384271 :       switch(ty)
    2623             :       {
    2624      341822 :         case t_INT: z = cgetg(3, t_FRAC);
    2625      341822 :         p1 = gcdii(y,gel(x,1));
    2626      341822 :         if (equali1(p1))
    2627             :         {
    2628      214514 :           avma = (pari_sp)z; tetpil = 0;
    2629      214514 :           gel(z,1) = icopy(gel(x,1));
    2630             :         }
    2631             :         else
    2632             :         {
    2633      127308 :           y = diviiexact(y,p1); tetpil = avma;
    2634      127308 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2635             :         }
    2636      341822 :         gel(z,2) = mulii(gel(x,2),y);
    2637      341822 :         normalize_frac(z);
    2638      341822 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2639      341822 :         return z;
    2640             : 
    2641             :         case t_REAL:
    2642       40069 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2643       40069 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2644             : 
    2645           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2646           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2647           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2648             :         }
    2649             : 
    2650          28 :         case t_FFELT: av=avma;
    2651          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2652             : 
    2653         427 :         case t_COMPLEX: return divRc(x, y);
    2654             : 
    2655             :         case t_PADIC:
    2656        1911 :           if (!signe(gel(x,1))) return gen_0;
    2657        1911 :           return divTp(x, y);
    2658             : 
    2659             :         case t_QUAD:
    2660           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2661           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2662             :       }
    2663             : 
    2664             :     case t_FFELT:
    2665        1841 :       switch (ty)
    2666             :       {
    2667        1764 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2668          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2669             :         case t_INTMOD:
    2670           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2671           0 :             pari_err_OP("/",x,y);
    2672           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2673             :         default:
    2674          42 :         pari_err_TYPE2("/",x,y);
    2675             :       }
    2676           0 :       break;
    2677             : 
    2678             :     case t_COMPLEX:
    2679     3209858 :       switch(ty)
    2680             :       {
    2681     3209858 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2682           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2683             :         case t_PADIC:
    2684           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2685             :         case t_QUAD:
    2686           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2687           0 :           return divfq(x, y, lx);
    2688             :       }
    2689             : 
    2690             :     case t_PADIC:
    2691        1988 :       switch(ty)
    2692             :       {
    2693        1939 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2694        1939 :           return signe(gel(x,4))? divpT(x, y)
    2695        1939 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2696             :         }
    2697           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2698           7 :           z = cgetg(3, t_INTMOD);
    2699           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2700             :         }
    2701             :         case t_COMPLEX: case t_QUAD:
    2702           7 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2703           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2704             : 
    2705          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2706             :       }
    2707             : 
    2708             :     case t_QUAD:
    2709         896 :       switch (ty)
    2710             :       {
    2711             :         case t_INT: case t_INTMOD: case t_FRAC:
    2712         847 :           z = cgetg(4,t_QUAD);
    2713         847 :           gel(z,1) = ZX_copy(gel(x,1));
    2714         847 :           gel(z,2) = gdiv(gel(x,2), y);
    2715         847 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2716          28 :         case t_REAL: return divqf(x, y, lg(y));
    2717           7 :         case t_PADIC: return divTp(x, y);
    2718             :         case t_COMPLEX:
    2719           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2720           0 :           return divqf(x, y, ly);
    2721             :       }
    2722             :   }
    2723    13496845 :   switch(ty) {
    2724             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2725       31275 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2726             :     case t_MAT:
    2727          42 :       av = avma; p1 = RgM_inv(y);
    2728          35 :       if (!p1) pari_err_INV("gdiv",y);
    2729          28 :       return gerepileupto(av, gmul(x, p1));
    2730             :     case t_VEC: case t_COL:
    2731             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2732           0 :       pari_err_TYPE2("/",x,y);
    2733             :   }
    2734    13465528 :   switch(tx) {
    2735             :     case t_VEC: case t_COL: case t_MAT:
    2736      904601 :       z = cgetg_copy(x, &lx);
    2737      904601 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2738      904601 :       return z;
    2739             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2740           0 :       pari_err_TYPE2("/",x,y);
    2741             :   }
    2742             : 
    2743    12560927 :   vy = gvar(y);
    2744    12560927 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2745       75411 :     vx = varn(X);
    2746       75411 :     if (vx != vy) {
    2747       75026 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2748       75012 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2749             :     }
    2750             :     /* y is POL, SER or RFRAC */
    2751         385 :     av = avma;
    2752         385 :     switch(ty)
    2753             :     {
    2754           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2755         385 :       default: y = ginvmod(gmod(y,X), X);
    2756             :     }
    2757         378 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2758             :   }
    2759             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2760             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2761             :    * variable NO_VARIABLE, then the operation is incorrect */
    2762    12485516 :   vx = gvar(x);
    2763    12485516 :   if (vx != vy) { /* includes cases where one is scalar */
    2764    11690263 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2765    10209282 :                         else return div_scal_T(x, y, ty);
    2766             :   }
    2767      795253 :   switch(tx)
    2768             :   {
    2769             :     case t_POL:
    2770       79764 :       switch(ty)
    2771             :       {
    2772             :         case t_SER:
    2773          63 :           if (lg(y) == 2)
    2774           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2775          63 :           p1 = RgX_to_ser(x,lg(y));
    2776          63 :           p2 = div_ser(p1, y, vx);
    2777          56 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2778          56 :           return p2;
    2779             : 
    2780             :         case t_RFRAC:
    2781             :         {
    2782       79701 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2783       79701 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2784           0 :             return mul_rfrac_scal(y2, y1, x);
    2785       79701 :           av = avma;
    2786       79701 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2787             :         }
    2788             :       }
    2789           0 :       break;
    2790             : 
    2791             :     case t_SER:
    2792      445620 :       switch(ty)
    2793             :       {
    2794             :         case t_POL:
    2795      445606 :           if (lg(x) == 2)
    2796         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2797      445501 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2798      445501 :           p2 = div_ser(x, p1, vx);
    2799      445501 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2800      445501 :           return p2;
    2801             :         case t_RFRAC:
    2802          14 :           av = avma;
    2803          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2804             :       }
    2805           0 :       break;
    2806             : 
    2807             :     case t_RFRAC:
    2808      269855 :       switch(ty)
    2809             :       {
    2810      269855 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2811             :         case t_SER:
    2812           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2813           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2814             :       }
    2815           0 :       break;
    2816             :   }
    2817          14 :   pari_err_TYPE2("/",x,y);
    2818             :   return NULL; /* LCOV_EXCL_LINE */
    2819             : }
    2820             : 
    2821             : /********************************************************************/
    2822             : /**                                                                **/
    2823             : /**                     SIMPLE MULTIPLICATION                      **/
    2824             : /**                                                                **/
    2825             : /********************************************************************/
    2826             : GEN
    2827    98787402 : gmulsg(long s, GEN y)
    2828             : {
    2829             :   long ly, i;
    2830             :   pari_sp av;
    2831             :   GEN z;
    2832             : 
    2833    98787402 :   switch(typ(y))
    2834             :   {
    2835    68297403 :     case t_INT:  return mulsi(s,y);
    2836    23013086 :     case t_REAL: return mulsr(s,y);
    2837      172207 :     case t_INTMOD: { GEN p = gel(y,1);
    2838      172207 :       z = cgetg(3,t_INTMOD);
    2839      172201 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2840      172193 :       gel(z,1) = icopy(p); return z;
    2841             :     }
    2842      627134 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2843             :     case t_FRAC:
    2844      922596 :       if (!s) return gen_0;
    2845      918193 :       z = cgetg(3,t_FRAC);
    2846      918193 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2847      918193 :       if (i == 1)
    2848             :       {
    2849      474400 :         gel(z,2) = icopy(gel(y,2));
    2850      474400 :         gel(z,1) = mulis(gel(y,1), s);
    2851             :       }
    2852             :       else
    2853             :       {
    2854      443793 :         gel(z,2) = divis(gel(y,2), i);
    2855      443793 :         gel(z,1) = mulis(gel(y,1), s/i);
    2856      443793 :         fix_frac_if_int(z);
    2857             :       }
    2858      918193 :       return z;
    2859             : 
    2860     4665986 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2861     4665986 :       gel(z,1) = gmulsg(s,gel(y,1));
    2862     4665986 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2863             : 
    2864             :     case t_PADIC:
    2865        9212 :       if (!s) return gen_0;
    2866        9212 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2867             : 
    2868           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2869           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2870           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2871           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2872             : 
    2873             :     case t_POLMOD:
    2874       75474 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2875             : 
    2876             :     case t_POL:
    2877      619400 :       if (!signe(y)) return RgX_copy(y);
    2878      603307 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2879      540874 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2880      540872 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2881      540869 :       return normalizepol_lg(z, ly);
    2882             : 
    2883             :     case t_SER:
    2884         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2885         182 :       if (!s) return Rg_get_0(y);
    2886         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2887         182 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2888         182 :       return normalize(z);
    2889             : 
    2890             :     case t_RFRAC:
    2891          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2892          14 :       if (s == 1) return gcopy(y);
    2893          14 :       if (s == -1) return gneg(y);
    2894          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2895             : 
    2896             :     case t_VEC: case t_COL: case t_MAT:
    2897      384701 :       z = cgetg_copy(y, &ly);
    2898      384701 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2899      384701 :       return z;
    2900             :   }
    2901           0 :   pari_err_TYPE("gmulsg",y);
    2902             :   return NULL; /* LCOV_EXCL_LINE */
    2903             : }
    2904             : 
    2905             : /********************************************************************/
    2906             : /**                                                                **/
    2907             : /**                       SIMPLE DIVISION                          **/
    2908             : /**                                                                **/
    2909             : /********************************************************************/
    2910             : 
    2911             : GEN
    2912    10623793 : gdivgs(GEN x, long s)
    2913             : {
    2914    10623793 :   long tx = typ(x), lx, i;
    2915             :   pari_sp av;
    2916             :   GEN z, y;
    2917             : 
    2918    10623793 :   if (!s)
    2919             :   {
    2920           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2921           0 :     pari_err_INV("gdivgs",gen_0);
    2922             :   }
    2923    10623793 :   switch(tx)
    2924             :   {
    2925             :     case t_INT:
    2926     6059981 :       av = avma; z = divis_rem(x,s,&i);
    2927     6059981 :       if (!i) return z;
    2928             : 
    2929     4051719 :       i = cgcd(s, i);
    2930     4051719 :       avma=av; z = cgetg(3,t_FRAC);
    2931     4051719 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    2932     4051719 :       gel(z,1) = y;
    2933     4051719 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    2934             : 
    2935             :     case t_REAL:
    2936     3148424 :       return divrs(x,s);
    2937             : 
    2938             :     case t_INTMOD:
    2939        7658 :       z = cgetg(3, t_INTMOD);
    2940        7658 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    2941             : 
    2942        1043 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    2943             : 
    2944      305387 :     case t_FRAC: z = cgetg(3, t_FRAC);
    2945      305387 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    2946      305387 :       if (i == 1)
    2947             :       {
    2948       64856 :         gel(z,2) = mulsi(s, gel(x,2));
    2949       64856 :         gel(z,1) = icopy(gel(x,1));
    2950             :       }
    2951             :       else
    2952             :       {
    2953      240531 :         gel(z,2) = mulsi(s/i, gel(x,2));
    2954      240531 :         gel(z,1) = divis(gel(x,1), i);
    2955             :       }
    2956      305387 :       normalize_frac(z);
    2957      305387 :       fix_frac_if_int(z); return z;
    2958             : 
    2959      922861 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2960      922861 :       gel(z,1) = gdivgs(gel(x,1),s);
    2961      922861 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    2962             : 
    2963             :     case t_PADIC: /* divpT */
    2964             :     {
    2965       87320 :       GEN p = gel(x,2);
    2966       87320 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    2967       86599 :       av = avma;
    2968       86599 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    2969             :     }
    2970             : 
    2971           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2972           0 :       gel(z,1) = ZX_copy(gel(x,1));
    2973           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    2974           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    2975             : 
    2976             :     case t_POLMOD:
    2977       22568 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    2978             : 
    2979             :     case t_RFRAC:
    2980         140 :       if (s == 1) return gcopy(x);
    2981         133 :       else if (s == -1) return gneg(x);
    2982         133 :       return div_rfrac_scal(x, stoi(s));
    2983             : 
    2984             :     case t_POL: case t_SER:
    2985       68110 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    2986       68110 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2987       68110 :       return z;
    2988             :     case t_VEC: case t_COL: case t_MAT:
    2989         301 :       z = cgetg_copy(x, &lx);
    2990         301 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2991         301 :       return z;
    2992             : 
    2993             :   }
    2994           0 :   pari_err_TYPE2("/",x, stoi(s));
    2995             :   return NULL; /* LCOV_EXCL_LINE */
    2996             : }
    2997             : 
    2998             : /* True shift (exact multiplication by 2^n) */
    2999             : GEN
    3000    48982419 : gmul2n(GEN x, long n)
    3001             : {
    3002             :   long lx, i, k, l;
    3003             :   GEN z, a, b;
    3004             : 
    3005    48982419 :   switch(typ(x))
    3006             :   {
    3007             :     case t_INT:
    3008    35311142 :       if (n>=0) return shifti(x,n);
    3009     1512039 :       if (!signe(x)) return gen_0;
    3010     1454404 :       l = vali(x); n = -n;
    3011     1454404 :       if (n<=l) return shifti(x,-n);
    3012      231518 :       z = cgetg(3,t_FRAC);
    3013      231518 :       gel(z,1) = shifti(x,-l);
    3014      231518 :       gel(z,2) = int2n(n-l); return z;
    3015             : 
    3016             :     case t_REAL:
    3017     6526118 :       return shiftr(x,n);
    3018             : 
    3019       75320 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3020       75320 :       z = cgetg(3,t_INTMOD);
    3021       75320 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3022       75054 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3023       75056 :       gel(z,1) = icopy(b); return z;
    3024             : 
    3025      216859 :     case t_FFELT: return FF_mul2n(x,n);
    3026             : 
    3027      449988 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3028      449988 :       l = vali(a);
    3029      449988 :       k = vali(b);
    3030      449988 :       if (n+l >= k)
    3031             :       {
    3032       85862 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3033       72527 :         l = n-k; k = -k;
    3034             :       }
    3035             :       else
    3036             :       {
    3037      364126 :         k = -(l+n); l = -l;
    3038             :       }
    3039      436653 :       z = cgetg(3,t_FRAC);
    3040      436653 :       gel(z,1) = shifti(a,l);
    3041      436653 :       gel(z,2) = shifti(b,k); return z;
    3042             : 
    3043     4997026 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3044     4997026 :       gel(z,1) = gmul2n(gel(x,1),n);
    3045     4997026 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3046             : 
    3047          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3048          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3049          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3050          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3051             : 
    3052             :     case t_POLMOD:
    3053       18361 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3054             : 
    3055             :     case t_POL:
    3056      392855 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3057      392855 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3058      392855 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3059             :     case t_SER:
    3060      121978 :       if (ser_isexactzero(x)) return gcopy(x);
    3061      121964 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3062      121964 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3063      121964 :       return normalize(z); /* needed if char = 2 */
    3064             :     case t_VEC: case t_COL: case t_MAT:
    3065      855745 :       z = cgetg_copy(x, &lx);
    3066      855745 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3067      855745 :       return z;
    3068             : 
    3069             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3070          14 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3071             : 
    3072             :     case t_PADIC: /* int2n wrong if n < 0 */
    3073       16999 :       return gmul(gmul2n(gen_1,n),x);
    3074             :   }
    3075           0 :   pari_err_TYPE("gmul2n",x);
    3076             :   return NULL; /* LCOV_EXCL_LINE */
    3077             : }
    3078             : 
    3079             : /*******************************************************************/
    3080             : /*                                                                 */
    3081             : /*                              INVERSE                            */
    3082             : /*                                                                 */
    3083             : /*******************************************************************/
    3084             : static GEN
    3085       85155 : inv_polmod(GEN T, GEN x)
    3086             : {
    3087       85155 :   GEN z = cgetg(3,t_POLMOD), a;
    3088       85155 :   gel(z,1) = RgX_copy(T);
    3089       85155 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3090       62048 :     a = ginv(x);
    3091             :   else
    3092             :   {
    3093       23107 :     if (lg(T) == 5) /* quadratic fields */
    3094       17913 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3095             :     else
    3096        5194 :       a = RgXQ_inv(x, T);
    3097             :   }
    3098       85155 :   gel(z,2) = a; return z;
    3099             : }
    3100             : 
    3101             : GEN
    3102    15490633 : ginv(GEN x)
    3103             : {
    3104             :   long s;
    3105             :   pari_sp av, tetpil;
    3106             :   GEN z, y, p1, p2;
    3107             : 
    3108    15490633 :   switch(typ(x))
    3109             :   {
    3110             :     case t_INT:
    3111      548950 :       if (is_pm1(x)) return icopy(x);
    3112      200511 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3113      200511 :       z = cgetg(3,t_FRAC);
    3114      200511 :       gel(z,1) = s<0? gen_m1: gen_1;
    3115      200511 :       gel(z,2) = absi(x); return z;
    3116             : 
    3117     2971796 :     case t_REAL: return invr(x);
    3118             : 
    3119         679 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3120         679 :       gel(z,1) = icopy(gel(x,1));
    3121         679 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3122             : 
    3123             :     case t_FRAC: {
    3124      235253 :       GEN a = gel(x,1), b = gel(x,2);
    3125      235253 :       s = signe(a);
    3126      235253 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3127       99956 :       z = cgetg(3,t_FRAC);
    3128       99956 :       gel(z,1) = icopy(b);
    3129       99956 :       gel(z,2) = icopy(a);
    3130       99956 :       normalize_frac(z); return z;
    3131             :     }
    3132             :     case t_COMPLEX:
    3133     1457786 :       av=avma;
    3134     1457786 :       p1=cxnorm(x);
    3135     1457786 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3136     1457786 :       tetpil=avma;
    3137     1457786 :       return gerepile(av,tetpil,divcR(p2,p1));
    3138             : 
    3139             :     case t_QUAD:
    3140         203 :       av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
    3141         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3142             : 
    3143        4417 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3144        4417 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3145        4410 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3146        4410 :       gel(z,2) = icopy(gel(x,2));
    3147        4410 :       gel(z,3) = icopy(gel(x,3));
    3148        4410 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3149             : 
    3150       85155 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3151        7671 :     case t_FFELT: return FF_inv(x);
    3152     5371974 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3153       15106 :     case t_SER: return gdiv(gen_1,x);
    3154             : 
    3155             :     case t_RFRAC:
    3156             :     {
    3157        2912 :       GEN n = gel(x,1), d = gel(x,2);
    3158        2912 :       pari_sp av = avma, ltop;
    3159        2912 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3160             : 
    3161        2912 :       n = simplify_shallow(n);
    3162        2912 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3163             :       {
    3164        2912 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3165         679 :         ltop = avma;
    3166         679 :         z = RgX_Rg_div(d,n);
    3167             :       } else {
    3168           0 :         ltop = avma;
    3169           0 :         z = cgetg(3,t_RFRAC);
    3170           0 :         gel(z,1) = RgX_copy(d);
    3171           0 :         gel(z,2) = RgX_copy(n);
    3172             :       }
    3173         679 :       stackdummy(av, ltop);
    3174         679 :       return z;
    3175             :     }
    3176             : 
    3177             :     case t_QFR:
    3178          14 :       av = avma; z = cgetg(5, t_QFR);
    3179          14 :       gel(z,1) = gel(x,1);
    3180          14 :       gel(z,2) = negi( gel(x,2) );
    3181          14 :       gel(z,3) = gel(x,3);
    3182          14 :       gel(z,4) = negr( gel(x,4) );
    3183          14 :       return gerepileupto(av, redreal(z));
    3184             : 
    3185             :     case t_QFI:
    3186     4784951 :       y = gcopy(x);
    3187     4784951 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3188     4670019 :         togglesign(gel(y,2));
    3189     4784951 :       return y;
    3190             :     case t_MAT:
    3191        3766 :       y = RgM_inv(x);
    3192        3759 :       if (!y) pari_err_INV("ginv",x);
    3193        3689 :       return y;
    3194             :     case t_VECSMALL:
    3195             :     {
    3196           0 :       long i, lx = lg(x)-1;
    3197           0 :       y = zero_zv(lx);
    3198           0 :       for (i=1; i<=lx; i++)
    3199             :       {
    3200           0 :         long xi = x[i];
    3201           0 :         if (xi<1 || xi>lx || y[xi])
    3202           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3203           0 :         y[xi] = i;
    3204             :       }
    3205           0 :       return y;
    3206             :     }
    3207             :   }
    3208           0 :   pari_err_TYPE("inverse",x);
    3209             :   return NULL; /* LCOV_EXCL_LINE */
    3210             : }

Generated by: LCOV version 1.13