Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - gen1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23157-9c2a597eb) Lines: 1710 1812 94.4 %
Date: 2018-10-16 05:38:41 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          21 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
      43             : 
      44             : /* is -1 not a square in Zp, assume p prime */
      45             : INLINE int
      46          28 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      47             : 
      48             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      49             : static GEN addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN));
      50             : static GEN mulpp(GEN x, GEN y);
      51             : static GEN divpp(GEN x, GEN y);
      52             : /* Argument codes for inline routines
      53             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      54             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      55             :  * T: some type (to be converted to PADIC)
      56             :  */
      57             : static GEN
      58    74578534 : addRc(GEN x, GEN y) {
      59    74578534 :   GEN z = cgetg(3,t_COMPLEX);
      60    74578534 :   gel(z,1) = gadd(x,gel(y,1));
      61    74578534 :   gel(z,2) = gcopy(gel(y,2)); return z;
      62             : }
      63             : static GEN
      64    74805872 : mulRc(GEN x, GEN y) {
      65    74805872 :   GEN z = cgetg(3,t_COMPLEX);
      66    74805872 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      67    74805872 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      68             : }
      69             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      70             : static GEN
      71          49 : mulRc_direct(GEN x, GEN y) {
      72          49 :   GEN z = cgetg(3,t_COMPLEX);
      73          49 :   gel(z,1) = gmul(x,gel(y,1));
      74          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      75             : }
      76             : static GEN
      77      563814 : divRc(GEN x, GEN y) {
      78      563814 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      79      563814 :   GEN z = cgetg(3,t_COMPLEX);
      80      563814 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      81      563814 :   gel(z,2) = gmul(mt, gel(y,2));
      82      563814 :   return z;
      83             : }
      84             : static GEN
      85     6670776 : divcR(GEN x, GEN y) {
      86     6670776 :   GEN z = cgetg(3,t_COMPLEX);
      87     6670776 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      88     6670776 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      89             : }
      90             : static GEN
      91         763 : addRq(GEN x, GEN y) {
      92         763 :   GEN z = cgetg(4,t_QUAD);
      93         763 :   gel(z,1) = ZX_copy(gel(y,1));
      94         763 :   gel(z,2) = gadd(x, gel(y,2));
      95         763 :   gel(z,3) = gcopy(gel(y,3)); return z;
      96             : }
      97             : static GEN
      98        1736 : mulRq(GEN x, GEN y) {
      99        1736 :   GEN z = cgetg(4,t_QUAD);
     100        1736 :   gel(z,1) = ZX_copy(gel(y,1));
     101        1736 :   gel(z,2) = gmul(x,gel(y,2));
     102        1736 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     103             : }
     104             : static GEN
     105          28 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     106          28 :   long i = gexpo(x) - gexpo(y);
     107          28 :   if (i > 0) prec += nbits2extraprec( i );
     108          28 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     109             : }
     110             : static GEN
     111     9928780 : mulrfrac(GEN x, GEN y)
     112             : {
     113             :   pari_sp av;
     114     9928780 :   GEN z, a = gel(y,1), b = gel(y,2);
     115     9928780 :   if (is_pm1(a)) /* frequent special case */
     116             :   {
     117     2194595 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     118     2194595 :     return z;
     119             :   }
     120     7734185 :   av = avma;
     121     7734185 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     122             : }
     123             : static GEN
     124          14 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     125          14 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     126             : }
     127             : static GEN
     128          28 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     129          28 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     130             : }
     131             : static GEN
     132           7 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     133           7 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     134             : }
     135             : /* y PADIC, x + y by converting x to padic */
     136             : static GEN
     137           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     138             : 
     139           7 :   if (!valp(y)) z = cvtop2(x,y);
     140             :   else {
     141           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     142           7 :     z  = cvtop(x, gel(y,2), l);
     143             :   }
     144           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     145             : }
     146             : /* y PADIC, x * y by converting x to padic */
     147             : static GEN
     148      162718 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     149      162718 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     150             : }
     151             : /* y PADIC, non zero x / y by converting x to padic */
     152             : static GEN
     153       14119 : divTp(GEN x, GEN y) { pari_sp av = avma;
     154       14119 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     155             : }
     156             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     157             :  * converted to O(p^e) and division by 0 */
     158             : static GEN
     159        1820 : divpT(GEN x, GEN y) { pari_sp av = avma;
     160        1820 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     161             : }
     162             : 
     163             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     164             :  * clean memory from z on */
     165             : static GEN
     166     2236044 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     167     2236044 :   if (lgefint(X) == 3) {
     168     2183165 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     169     2183165 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     170             :   }
     171             :   else {
     172       52879 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     173       52873 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     174             :   }
     175     2236065 :   gel(z,1) = icopy(X); return z;
     176             : }
     177             : static GEN
     178     1102088 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     179     1102088 :   if (lgefint(X) == 3) {
     180      751075 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     181      751075 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     182             :   }
     183             :   else {
     184      351013 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     185      351003 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     186             :   }
     187     1102097 :   gel(z,1) = icopy(X); return z;
     188             : }
     189             : /* cf add_intmod_same */
     190             : static GEN
     191     2689532 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     192     2689532 :   if (lgefint(X) == 3) {
     193     2592547 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     194     2592547 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     195             :   }
     196             :   else
     197       96985 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     198     2689569 :   gel(z,1) = icopy(X); return z;
     199             : }
     200             : /* cf add_intmod_same */
     201             : static GEN
     202       59540 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     203             : {
     204       59540 :   if (lgefint(X) == 3) {
     205       51296 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     206       51289 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     207             :   }
     208             :   else
     209        8244 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     210       59538 :   gel(z,1) = icopy(X); return z;
     211             : }
     212             : 
     213             : /*******************************************************************/
     214             : /*                                                                 */
     215             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     216             : /*                                                                 */
     217             : /* (static routines are not memory clean, but OK for gerepileupto) */
     218             : /*******************************************************************/
     219             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     220             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     221             : static GEN
     222    10077139 : rfrac_denom_mul_scal(GEN d, GEN y)
     223             : {
     224    10077139 :   GEN D = RgX_Rg_mul(d, y);
     225    10077139 :   if (lg(D) != lg(d))
     226             :   { /* try to generate a meaningful diagnostic */
     227           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     228           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     229             :   }
     230    10077139 :   return D;
     231             : }
     232             : 
     233             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     234             : GEN
     235    58060793 : gred_rfrac_simple(GEN n, GEN d)
     236             : {
     237             :   GEN c, cn, cd, z;
     238    58060793 :   long dd = degpol(d);
     239             : 
     240    58060793 :   if (dd <= 0)
     241             :   {
     242       49427 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     243       49427 :     n = gdiv(n, gel(d,2));
     244       49427 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     245       49427 :     return n;
     246             :   }
     247             : 
     248    58011366 :   cd = content(d);
     249    58011366 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     250    58011366 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     251    58011366 :   if (!gequal1(cd)) {
     252     6576369 :     d = RgX_Rg_div(d,cd);
     253     6576369 :     if (!gequal1(cn))
     254             :     {
     255     1302605 :       if (gequal0(cn)) {
     256          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     257           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     258           0 :         c = gen_1;
     259             :       } else {
     260     1302556 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     261     1302556 :         c = gdiv(cn,cd);
     262             :       }
     263             :     }
     264             :     else
     265     5273764 :       c = ginv(cd);
     266             :   } else {
     267    51434997 :     if (!gequal1(cn))
     268             :     {
     269     3230189 :       if (gequal0(cn)) {
     270         413 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     271          21 :         c = gen_1;
     272             :       } else {
     273     3229776 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     274     3229776 :         c = cn;
     275             :       }
     276             :     } else {
     277    48204808 :       GEN y = cgetg(3,t_RFRAC);
     278    48204808 :       gel(y,1) = gcopy(n);
     279    48204808 :       gel(y,2) = RgX_copy(d); return y;
     280             :     }
     281             :   }
     282             : 
     283     9806117 :   if (typ(c) == t_POL)
     284             :   {
     285      905469 :     z = c;
     286      944907 :     do { z = content(z); } while (typ(z) == t_POL);
     287      905469 :     cd = denom_i(z);
     288      905469 :     cn = gmul(c, cd);
     289             :   }
     290             :   else
     291             :   {
     292     8900648 :     cn = numer_i(c);
     293     8900648 :     cd = denom_i(c);
     294             :   }
     295     9806117 :   z = cgetg(3,t_RFRAC);
     296     9806117 :   gel(z,1) = gmul(n, cn);
     297     9806117 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     298     9806117 :   return z;
     299             : }
     300             : 
     301             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     302             : static GEN
     303      249074 : fix_rfrac(GEN x, long d)
     304             : {
     305             :   GEN z, N, D;
     306      249074 :   if (!d || typ(x) == t_POL) return x;
     307      175630 :   z = cgetg(3, t_RFRAC);
     308      175630 :   N = gel(x,1);
     309      175630 :   D = gel(x,2);
     310      175630 :   if (d > 0) {
     311      344267 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     312      177975 :                                                    : monomialcopy(N,d,varn(D));
     313      166243 :     gel(z, 2) = RgX_copy(D);
     314             :   } else {
     315        9387 :     gel(z, 1) = gcopy(N);
     316        9387 :     gel(z, 2) = RgX_shift(D, -d);
     317             :   }
     318      175630 :   return z;
     319             : }
     320             : 
     321             : /* assume d != 0 */
     322             : static GEN
     323    44749348 : gred_rfrac2(GEN n, GEN d)
     324             : {
     325             :   GEN y, z;
     326             :   long v, vd, vn;
     327             : 
     328    44749348 :   n = simplify_shallow(n);
     329    44749348 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     330    37918783 :   d = simplify_shallow(d);
     331    37918783 :   if (typ(d) != t_POL) return gdiv(n,d);
     332    36739127 :   vd = varn(d);
     333    36739127 :   if (typ(n) != t_POL)
     334             :   {
     335    20613217 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     336    20611810 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     337           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     338             :   }
     339    16125910 :   vn = varn(n);
     340    16125910 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     341    15988292 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     342             : 
     343             :   /* now n and d are t_POLs in the same variable */
     344    15828540 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     345    15828540 :   if (!degpol(d))
     346             :   {
     347    12755842 :     n = RgX_Rg_div(n,gel(d,2));
     348    12755842 :     return v? RgX_mulXn(n,v): n;
     349             :   }
     350             : 
     351             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     352     3072698 :   if (!isinexact(n) && !isinexact(d))
     353             :   {
     354     3072488 :     y = RgX_divrem(n, d, &z);
     355     3072488 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     356      248864 :     z = RgX_gcd(d, z);
     357      248864 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     358             :   }
     359      249074 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     360             : }
     361             : 
     362             : /* x,y t_INT, return x/y in reduced form */
     363             : GEN
     364    70238565 : Qdivii(GEN x, GEN y)
     365             : {
     366    70238565 :   pari_sp av = avma;
     367             :   GEN r, q;
     368             : 
     369    70238565 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     370    32941719 :   if (equali1(x)) {
     371     5317726 :     long s = signe(y);
     372     5317726 :     if (!s) pari_err_INV("gdiv",y);
     373     5317621 :     if (signe(x) < 0) s = -s;
     374     5317621 :     q = cgetg(3, t_FRAC);
     375     5317621 :     gel(q,1) = s<0? gen_m1: gen_1;
     376     5317621 :     gel(q,2) = absi(y); return q;
     377             :   }
     378    27623993 :   q = dvmdii(x,y,&r);
     379    27623993 :   if (r == gen_0) return q; /* gen_0 intended */
     380    12034078 :   r = gcdii(y,r);
     381    12034078 :   if (lgefint(r) == 3)
     382             :   {
     383    11798804 :     ulong t = r[2];
     384    11798804 :     set_avma(av);
     385    11798804 :     if (t == 1) q = mkfraccopy(x,y);
     386             :     else
     387             :     {
     388     5162916 :       q = cgetg(3,t_FRAC);
     389     5162916 :       gel(q,1) = diviuexact(x,t);
     390     5162916 :       gel(q,2) = diviuexact(y,t);
     391             :     }
     392             :   }
     393             :   else
     394             :   { /* rare: r and q left on stack for efficiency */
     395      235274 :     q = cgetg(3,t_FRAC);
     396      235274 :     gel(q,1) = diviiexact(x,r);
     397      235274 :     gel(q,2) = diviiexact(y,r);
     398             :   }
     399    12034078 :   normalize_frac(q); return q;
     400             : }
     401             : 
     402             : /*******************************************************************/
     403             : /*                                                                 */
     404             : /*                          CONJUGATION                            */
     405             : /*                                                                 */
     406             : /*******************************************************************/
     407             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     408             : static GEN
     409       23877 : quad_polmod_conj(GEN x, GEN y)
     410             : {
     411             :   GEN z, u, v, a, b;
     412       23877 :   if (typ(x) != t_POL) return gcopy(x);
     413       23877 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     414       23877 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     415       23877 :   b = gel(y,3); v = gel(x,2);
     416       23877 :   z = cgetg(4, t_POL); z[1] = x[1];
     417       23877 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     418       23877 :   gel(z,3) = gneg(u); return z;
     419             : }
     420             : static GEN
     421       23877 : quad_polmod_norm(GEN x, GEN y)
     422             : {
     423             :   GEN z, u, v, a, b, c;
     424       23877 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     425       23877 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     426       23877 :   b = gel(y,3); v = gel(x,2);
     427       23877 :   c = gel(y,2);
     428       23877 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     429       23877 :   if (!gequal1(a)) z = gdiv(z, a);
     430       23877 :   return gadd(z, gsqr(v));
     431             : }
     432             : 
     433             : GEN
     434     8948576 : conj_i(GEN x)
     435             : {
     436             :   long lx, i;
     437             :   GEN y;
     438             : 
     439     8948576 :   switch(typ(x))
     440             :   {
     441             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     442     4567969 :       return x;
     443             : 
     444     4351442 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     445             : 
     446             :     case t_QUAD:
     447         595 :       y = cgetg(4,t_QUAD);
     448         595 :       gel(y,1) = gel(x,1);
     449        1190 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     450         595 :                                       : gadd(gel(x,2), gel(x,3));
     451         595 :       gel(y,3) = gneg(gel(x,3)); return y;
     452             : 
     453             :     case t_POL: case t_SER:
     454       10507 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     455       10507 :       for (i=2; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     456       10507 :       return y;
     457             : 
     458             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     459       18063 :       y = cgetg_copy(x, &lx);
     460       18063 :       for (i=1; i<lx; i++) gel(y,i) = conj_i(gel(x,i));
     461       18063 :       return y;
     462             : 
     463             :     case t_POLMOD:
     464             :     {
     465           0 :       GEN X = gel(x,1);
     466           0 :       long d = degpol(X);
     467           0 :       if (d < 2) return x;
     468           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     469             :     }
     470             :   }
     471           0 :   pari_err_TYPE("gconj",x);
     472             :   return NULL; /* LCOV_EXCL_LINE */
     473             : }
     474             : GEN
     475       30431 : gconj(GEN x)
     476       30431 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     477             : 
     478             : GEN
     479          84 : conjvec(GEN x,long prec)
     480             : {
     481             :   long lx, s, i;
     482             :   GEN z;
     483             : 
     484          84 :   switch(typ(x))
     485             :   {
     486             :     case t_INT: case t_INTMOD: case t_FRAC:
     487           0 :       return mkcolcopy(x);
     488             : 
     489             :     case t_COMPLEX: case t_QUAD:
     490           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     491             : 
     492             :     case t_FFELT:
     493          28 :       return FF_conjvec(x);
     494             : 
     495             :     case t_VEC: case t_COL:
     496           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     497           0 :       if (lx == 1) return z;
     498           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     499           0 :       s = lgcols(z);
     500           0 :       for (i=2; i<lx; i++)
     501             :       {
     502           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     503           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     504             :       }
     505           0 :       break;
     506             : 
     507             :     case t_POLMOD: {
     508          56 :       GEN T = gel(x,1), r;
     509             :       pari_sp av;
     510             : 
     511          56 :       lx = lg(T);
     512          56 :       if (lx <= 3) return cgetg(1,t_COL);
     513          56 :       x = gel(x,2);
     514         238 :       for (i=2; i<lx; i++)
     515             :       {
     516         189 :         GEN c = gel(T,i);
     517         189 :         switch(typ(c)) {
     518             :           case t_INTMOD: {
     519           7 :             GEN p = gel(c,1);
     520             :             pari_sp av;
     521           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     522           7 :             av = avma;
     523           7 :             T = RgX_to_FpX(T,p);
     524           7 :             x = RgX_to_FpX(x, p);
     525           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     526           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     527           7 :             return gerepileupto(av, z);
     528             :           }
     529             :           case t_INT:
     530         182 :           case t_FRAC: break;
     531           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     532             :         }
     533             :       }
     534          49 :       if (typ(x) != t_POL)
     535             :       {
     536           0 :         if (!is_rational_t(typ(x)))
     537           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     538           0 :         retconst_col(lx-3, gcopy(x));
     539             :       }
     540          49 :       RgX_check_QX(x,"conjvec");
     541          49 :       av = avma;
     542          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     543          49 :       r = cleanroots(T,prec);
     544          49 :       z = cgetg(lx-2,t_COL);
     545          49 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     546          49 :       return gerepileupto(av, z);
     547             :     }
     548             : 
     549             :     default:
     550           0 :       pari_err_TYPE("conjvec",x);
     551             :       return NULL; /* LCOV_EXCL_LINE */
     552             :   }
     553           0 :   return z;
     554             : }
     555             : 
     556             : 
     557             : /********************************************************************/
     558             : /**                                                                **/
     559             : /**                           ADDITION                             **/
     560             : /**                                                                **/
     561             : /********************************************************************/
     562             : /* x, y compatible PADIC, op = add or sub */
     563             : static GEN
     564     1178489 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     565             : {
     566     1178489 :   pari_sp av = avma;
     567             :   long d,e,r,rx,ry;
     568     1178489 :   GEN u, z, mod, p = gel(x,2);
     569             :   int swap;
     570             : 
     571     1178489 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     572     1178489 :   e = valp(x);
     573     1178489 :   r = valp(y); d = r-e;
     574     1178489 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     575     1178489 :   rx = precp(x);
     576     1178489 :   ry = precp(y);
     577     1178489 :   if (d) /* v(x) < v(y) */
     578             :   {
     579      548054 :     r = d+ry; z = powiu(p,d);
     580      548054 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     581      548054 :     z = mulii(z,gel(y,4));
     582      548054 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     583             :   }
     584             :   else
     585             :   {
     586             :     long c;
     587      630435 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     588      630435 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     589      630435 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     590             :     {
     591       59171 :       set_avma(av); return zeropadic(p, e+r);
     592             :     }
     593      571264 :     if (c)
     594             :     {
     595       95333 :       mod = diviiexact(mod, powiu(p,c));
     596       95333 :       r -= c;
     597       95333 :       e += c;
     598             :     }
     599             :   }
     600     1119318 :   u = modii(u, mod);
     601     1119318 :   set_avma(av); z = cgetg(5,t_PADIC);
     602     1119318 :   z[1] = evalprecp(r) | evalvalp(e);
     603     1119318 :   gel(z,2) = icopy(p);
     604     1119318 :   gel(z,3) = icopy(mod);
     605     1119318 :   gel(z,4) = icopy(u); return z;
     606             : }
     607             : /* Rg_to_Fp(t_FRAC) without GC */
     608             : static GEN
     609       56483 : Q_to_Fp(GEN x, GEN p)
     610       56483 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     611             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     612             : static GEN
     613      370414 : addQp(GEN x, GEN y)
     614             : {
     615      370414 :   pari_sp av = avma;
     616      370414 :   long d, r, e, vy = valp(y), py = precp(y);
     617      370414 :   GEN z, q, mod, u, p = gel(y,2);
     618             : 
     619      370414 :   e = Q_pvalrem(x, p, &x);
     620      370393 :   d = vy - e; r = d + py;
     621      370393 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     622      368811 :   mod = gel(y,3);
     623      368811 :   u   = gel(y,4);
     624      368811 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     625             : 
     626      368827 :   if (d > 0)
     627             :   {
     628      234628 :     q = powiu(p,d);
     629      234611 :     mod = mulii(mod, q);
     630      234599 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     631      234599 :     u = addii(x,  mulii(u, q));
     632             :   }
     633      134199 :   else if (d < 0)
     634             :   {
     635       25193 :     q = powiu(p,-d);
     636       25193 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     637       25193 :     u = addii(u, mulii(x, q));
     638       25194 :     r = py; e = vy;
     639             :   }
     640             :   else
     641             :   {
     642             :     long c;
     643      109006 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     644      109006 :     u = addii(u, x);
     645      109001 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     646             :     {
     647         966 :       set_avma(av); return zeropadic(p,e+r);
     648             :     }
     649      108036 :     if (c)
     650             :     {
     651       58813 :       mod = diviiexact(mod, powiu(p,c));
     652       58817 :       r -= c;
     653       58817 :       e += c;
     654             :     }
     655             :   }
     656      367829 :   u = modii(u, mod); set_avma(av);
     657      367846 :   z = cgetg(5,t_PADIC);
     658      367869 :   z[1] = evalprecp(r) | evalvalp(e);
     659      367868 :   gel(z,2) = icopy(p);
     660      367875 :   gel(z,3) = icopy(mod);
     661      367874 :   gel(z,4) = icopy(u); return z;
     662             : }
     663             : 
     664             : /* Mod(x,X) + Mod(y,X) */
     665             : #define addsub_polmod_same addsub_polmod_scal
     666             : /* Mod(x,X) +/- Mod(y,Y) */
     667             : static GEN
     668        1750 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     669             : {
     670        1750 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     671        1750 :   GEN z = cgetg(3,t_POLMOD);
     672        1750 :   long vx = varn(X), vy = varn(Y);
     673        1750 :   if (vx==vy) {
     674             :     pari_sp av;
     675          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     676          14 :     warn_coercion(X,Y,gel(z,1));
     677          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     678             :   }
     679        1736 :   if (varncmp(vx, vy) < 0)
     680        1736 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     681             :   else
     682           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     683        1736 :   gel(z,2) = op(x, y); return z;
     684             : }
     685             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     686             : static GEN
     687     2010051 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     688     2010051 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     689             : 
     690             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     691             : static GEN
     692      308833 : add_ser_scal(GEN y, GEN x)
     693             : {
     694             :   long i, l, ly, vy;
     695             :   GEN z;
     696             : 
     697      308833 :   if (isrationalzero(x)) return gcopy(y);
     698      282534 :   ly = lg(y);
     699      282534 :   l = valp(y);
     700      282534 :   if (l < 3-ly) return gcopy(y);
     701             :   /* l + ly >= 3 */
     702      282366 :   if (l < 0)
     703             :   {
     704         735 :     z = cgetg(ly,t_SER); z[1] = y[1];
     705         735 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     706         735 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     707         735 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     708         735 :     return z;
     709             :   }
     710      281631 :   vy = varn(y);
     711      281631 :   if (l > 0)
     712             :   {
     713       21315 :     if (ser_isexactzero(y))
     714       13580 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     715        7735 :     y -= l; ly += l;
     716        7735 :     z = cgetg(ly,t_SER);
     717        7735 :     x = gcopy(x);
     718        7735 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     719             :   }
     720             :   else
     721             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     722      260316 :     z = cgetg(ly,t_SER);
     723      260316 :     x = gadd(x, gel(y,2));
     724      260316 :     i = 3;
     725             :   }
     726      268051 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     727      268051 :   gel(z,2) = x;
     728      268051 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     729      268051 :   return gequal0(x)? normalize(z): z;
     730             : }
     731             : static long
     732     3539016 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     733             : /* x,y t_SER in the same variable: x+y */
     734             : static GEN
     735     1769900 : ser_add(GEN x, GEN y)
     736             : {
     737     1769900 :   long i, lx,ly, n = valp(y) - valp(x);
     738             :   GEN z;
     739     1769900 :   if (n < 0) { n = -n; swap(x,y); }
     740             :   /* valp(x) <= valp(y) */
     741     1769900 :   lx = _serprec(x);
     742     1769900 :   if (lx == 2) /* don't lose type information */
     743         784 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), valp(x));
     744     1769116 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     745     1769116 :   if (n)
     746             :   {
     747       90926 :     if (n+2 > lx) return gcopy(x);
     748       90373 :     z = cgetg(ly,t_SER);
     749       90373 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     750       90373 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     751             :   } else {
     752     1678190 :     z = cgetg(ly,t_SER);
     753     1678190 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     754             :   }
     755     1768563 :   z[1] = x[1]; return normalize(z);
     756             : }
     757             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     758             : static GEN
     759     8823497 : add_rfrac_scal(GEN y, GEN x)
     760             : {
     761             :   pari_sp av;
     762             :   GEN n;
     763             : 
     764     8823497 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     765     5149221 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     766     5149221 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     767             : }
     768             : 
     769             : /* x "scalar", ty != t_MAT and non-scalar */
     770             : static GEN
     771    19570903 : add_scal(GEN y, GEN x, long ty)
     772             : {
     773    19570903 :   switch(ty)
     774             :   {
     775    11933848 :     case t_POL: return RgX_Rg_add(y, x);
     776      308819 :     case t_SER: return add_ser_scal(y, x);
     777     4470248 :     case t_RFRAC: return add_rfrac_scal(y, x);
     778           0 :     case t_COL: return RgC_Rg_add(y, x);
     779             :     case t_VEC:
     780     2857981 :       if (isintzero(x)) return gcopy(y);
     781         168 :       break;
     782             :   }
     783         175 :   pari_err_TYPE2("+",x,y);
     784             :   return NULL; /* LCOV_EXCL_LINE */
     785             : }
     786             : 
     787             : static GEN
     788    17896614 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     789             : {
     790    17896614 :   pari_sp av = avma;
     791    17896614 :   GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
     792    17896614 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     793    17896614 :   int s = cmpii(x2, y2);
     794             : 
     795    17896614 :   if (!s)
     796             :   { /* common denominator: (x1 op y1) / x2 */
     797     5493228 :     n = op(x1, y1);
     798     5493228 :     if (!signe(n)) { set_avma(av); return gen_0; }
     799     5023766 :     d = x2;
     800     5023766 :     q = dvmdii(n, d, &r);
     801     5023766 :     if (r == gen_0) { set_avma(av); return icopy(q); }
     802     3459468 :     r = gcdii(d, r);
     803     3459468 :     if (!equali1(r)) { n = diviiexact(n, r); d = diviiexact(d, r); }
     804     3459468 :     gel(z,1) = icopy_avma(n, (pari_sp)z);
     805     3459468 :     gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     806     3459468 :     avma = (pari_sp)gel(z,2); return z;
     807             :   }
     808    12403386 :   if (s < 0)
     809             :   {
     810     6297251 :     GEN Q = dvmdii(y2, x2, &r);
     811     6297251 :     if (r == gen_0)
     812             :     { /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     813     4212563 :       pari_sp av = avma;
     814     4212563 :       n = op(mulii(Q,x1), y1);
     815     4212563 :       q = dvmdii(n, x2, &r);
     816     4212563 :       if (r == gen_0)
     817             :       {
     818      613828 :         gel(z,1) = gerepileuptoint(av, q);
     819      613828 :         gel(z,2) = Q; return z;
     820             :       }
     821     3598735 :       r = gcdii(x2, r);
     822     3598735 :       if (!equali1(r)) { n = diviiexact(n, r); x2 = diviiexact(x2, r); }
     823     3598735 :       d = mulii(x2,Q);
     824     3598735 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     825     3598735 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     826     3598735 :       avma = (pari_sp)gel(z,2); return z;
     827             :     }
     828     2084688 :     delta = gcdii(x2,r);
     829             :   }
     830             :   else
     831             :   {
     832     6106135 :     GEN Q = dvmdii(x2, y2, &r);
     833     6106135 :     if (r == gen_0)
     834             :     { /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     835     5003638 :       pari_sp av = avma;
     836     5003638 :       n = op(x1, mulii(Q,y1));
     837     5003638 :       q = dvmdii(n, y2, &r);
     838     5003638 :       if (r == gen_0)
     839             :       {
     840      131497 :         gel(z,1) = gerepileuptoint(av, q);
     841      131497 :         gel(z,2) = Q; return z;
     842             :       }
     843     4872141 :       r = gcdii(y2, r);
     844     4872141 :       if (!equali1(r)) { n = diviiexact(n, r); y2 = diviiexact(y2, r); }
     845     4872141 :       d = mulii(y2,Q);
     846     4872141 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     847     4872141 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     848     4872141 :       avma = (pari_sp)gel(z,2); return z;
     849             :     }
     850     1102497 :     delta = gcdii(y2,r);
     851             :   }
     852             :   /* delta = gcd(x2,y2) */
     853     3187185 :   if (equali1(delta))
     854             :   { /* numerator is non-zero */
     855     1244857 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     856     1244857 :     gel(z,2) = mulii(x2,y2); return z;
     857             :   }
     858     1942328 :   x2 = diviiexact(x2,delta);
     859     1942328 :   y2 = diviiexact(y2,delta);
     860     1942328 :   n = op(mulii(x1,y2), mulii(y1,x2));
     861     1942328 :   if (!signe(n)) { set_avma(av); return gen_0; }
     862     1942328 :   d = mulii(x2, y2);
     863     1942328 :   q = dvmdii(n, delta, &r);
     864     1942328 :   if (r == gen_0)
     865             :   {
     866      144334 :     if (equali1(d)) { set_avma(av); return icopy(q); }
     867      144334 :     avma = (pari_sp)z;
     868      144334 :     gel(z,2) = icopy(d);
     869      144334 :     gel(z,1) = icopy(q); return z;
     870             :   }
     871     1797994 :   r = gcdii(delta, r);
     872     1797994 :   if (!equali1(r))
     873             :   {
     874      402525 :     n     = diviiexact(n, r);
     875      402525 :     delta = diviiexact(delta, r);
     876             :   }
     877     1797994 :   d = mulii(d,delta); avma = (pari_sp)z;
     878     1797994 :   gel(z,1) = icopy(n);
     879     1797994 :   gel(z,2) = icopy(d); return z;
     880             : }
     881             : 
     882             : /* assume x2, y2 are t_POLs in the same variable */
     883             : static GEN
     884     3037311 : add_rfrac(GEN x, GEN y)
     885             : {
     886     3037311 :   pari_sp av = avma;
     887     3037311 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     888     3037311 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     889             : 
     890     3037311 :   delta = RgX_gcd(x2,y2);
     891     3037311 :   if (!degpol(delta))
     892             :   {
     893         658 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     894         658 :     d = RgX_mul(x2, y2);
     895         658 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     896             :   }
     897     3036653 :   x2 = RgX_div(x2,delta);
     898     3036653 :   y2 = RgX_div(y2,delta);
     899     3036653 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     900     3036653 :   if (!signe(n))
     901             :   {
     902      721054 :     n = simplify_shallow(n);
     903      721054 :     if (isrationalzero(n)) return gerepileupto(av, zeropol(varn(x2)));
     904           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     905             :   }
     906     2315599 :   if (degpol(n) == 0)
     907     1149721 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     908     1165878 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     909     1165878 :   if (isexactzero(r))
     910             :   {
     911             :     GEN z;
     912      228045 :     d = RgX_mul(x2, y2);
     913             :     /* "constant" denominator ? */
     914      228045 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     915      228045 :     return gerepileupto(av, z);
     916             :   }
     917      937833 :   r = RgX_gcd(delta, r);
     918      937833 :   if (degpol(r))
     919             :   {
     920      160705 :     n = RgX_div(n, r);
     921      160705 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     922             :   }
     923             :   else
     924      777128 :     d = RgX_mul(gel(x,2), y2);
     925      937833 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     926             : }
     927             : 
     928             : GEN
     929  2028472072 : gadd(GEN x, GEN y)
     930             : {
     931  2028472072 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     932             :   pari_sp av, tetpil;
     933             :   GEN z, p1;
     934             : 
     935  2028472072 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     936             :   {
     937  1114443578 :     case t_INT: return addii(x,y);
     938   537205853 :     case t_REAL: return addrr(x,y);
     939     1518258 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     940     1518258 :       z = cgetg(3,t_INTMOD);
     941     1518253 :       if (X==Y || equalii(X,Y))
     942     1518215 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     943          14 :       gel(z,1) = gcdii(X,Y);
     944          14 :       warn_coercion(X,Y,gel(z,1));
     945          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     946          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     947             :     }
     948    13772224 :     case t_FRAC: return addsub_frac(x,y,addii);
     949    77430514 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     950    77430514 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     951    77430514 :       if (isintzero(gel(z,2)))
     952             :       {
     953      328937 :         avma = (pari_sp)(z+3);
     954      328937 :         return gadd(gel(x,1),gel(y,1));
     955             :       }
     956    77101577 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     957    77101577 :       return z;
     958             :     case t_PADIC:
     959      949627 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     960      949627 :       return addsub_pp(x,y, addii);
     961         476 :     case t_QUAD: z = cgetg(4,t_QUAD);
     962         476 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     963         476 :       gel(z,1) = ZX_copy(gel(x,1));
     964         476 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     965         476 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     966             :     case t_POLMOD:
     967      760165 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     968      758450 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     969        1715 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     970     9194350 :     case t_FFELT: return FF_add(x,y);
     971             :     case t_POL:
     972    17556755 :       vx = varn(x);
     973    17556755 :       vy = varn(y);
     974    17556755 :       if (vx != vy) {
     975      806448 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     976       18240 :         else                     return RgX_Rg_add(y, x);
     977             :       }
     978    16750307 :       return RgX_add(x, y);
     979             :     case t_SER:
     980     1765091 :       vx = varn(x);
     981     1765091 :       vy = varn(y);
     982     1765091 :       if (vx != vy) {
     983          14 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
     984          14 :         else                     return add_ser_scal(y, x);
     985             :       }
     986     1765077 :       return ser_add(x, y);
     987             :     case t_RFRAC:
     988     4329198 :       vx = varn(gel(x,2));
     989     4329198 :       vy = varn(gel(y,2));
     990     4329198 :       if (vx != vy) {
     991     1291887 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
     992      538397 :         else                     return add_rfrac_scal(y, x);
     993             :       }
     994     3037311 :       return add_rfrac(x,y);
     995             :     case t_VEC:
     996      589572 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
     997      589572 :       return RgV_add(x,y);
     998             :     case t_COL:
     999      857777 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1000      857777 :       return RgC_add(x,y);
    1001             :     case t_MAT:
    1002      671902 :       lx = lg(x);
    1003      671902 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1004      671902 :       if (lx == 1) return cgetg(1, t_MAT);
    1005      671902 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1006      671895 :       return RgM_add(x,y);
    1007             : 
    1008           0 :     default: pari_err_TYPE2("+",x,y);
    1009             :   }
    1010             :   /* tx != ty */
    1011   247426747 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1012             : 
    1013   247426747 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1014             :   {
    1015             :     case t_INT:
    1016   191540308 :       switch(ty)
    1017             :       {
    1018   113466169 :         case t_REAL: return addir(x,y);
    1019             :         case t_INTMOD:
    1020      714261 :           z = cgetg(3, t_INTMOD);
    1021      714261 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1022    20617570 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1023    20617570 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1024    20617570 :           gel(z,2) = icopy(gel(y,2)); return z;
    1025    54753151 :         case t_COMPLEX: return addRc(x, y);
    1026             :         case t_PADIC:
    1027      330631 :           if (!signe(x)) return gcopy(y);
    1028      299660 :           return addQp(x,y);
    1029         721 :         case t_QUAD: return addRq(x, y);
    1030     1657808 :         case t_FFELT: return FF_Z_add(y,x);
    1031             :       }
    1032             : 
    1033             :     case t_REAL:
    1034    22522701 :       switch(ty)
    1035             :       {
    1036             :         case t_FRAC:
    1037    12610495 :           if (!signe(gel(y,1))) return rcopy(x);
    1038    12610495 :           if (!signe(x))
    1039             :           {
    1040        4825 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1041        4825 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1042             :           }
    1043    12605670 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1044    12605670 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1045     9912185 :         case t_COMPLEX: return addRc(x, y);
    1046          21 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1047             : 
    1048           0 :         default: pari_err_TYPE2("+",x,y);
    1049             :       }
    1050             : 
    1051             :     case t_INTMOD:
    1052        3633 :       switch(ty)
    1053             :       {
    1054        3500 :         case t_FRAC: { GEN X = gel(x,1);
    1055        3500 :           z = cgetg(3, t_INTMOD);
    1056        3500 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1057        3500 :           return add_intmod_same(z, X, p1, gel(x,2));
    1058             :         }
    1059             :         case t_FFELT:
    1060          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1061           0 :             pari_err_OP("+",x,y);
    1062          14 :           return FF_Z_add(y,gel(x,2));
    1063          91 :         case t_COMPLEX: return addRc(x, y);
    1064           0 :         case t_PADIC: { GEN X = gel(x,1);
    1065           0 :           z = cgetg(3, t_INTMOD);
    1066           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1067             :         }
    1068          28 :         case t_QUAD: return addRq(x, y);
    1069             :       }
    1070             : 
    1071             :     case t_FRAC:
    1072     9985214 :       switch (ty)
    1073             :       {
    1074     9913100 :         case t_COMPLEX: return addRc(x, y);
    1075             :         case t_PADIC:
    1076       70756 :           if (!signe(gel(x,1))) return gcopy(y);
    1077       70756 :           return addQp(x,y);
    1078          14 :         case t_QUAD: return addRq(x, y);
    1079        1344 :         case t_FFELT: return FF_Q_add(y, x);
    1080             :       }
    1081             : 
    1082             :     case t_FFELT:
    1083           0 :       pari_err_TYPE2("+",x,y);
    1084             : 
    1085             :     case t_COMPLEX:
    1086          21 :       switch(ty)
    1087             :       {
    1088             :         case t_PADIC:
    1089          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1090             :         case t_QUAD:
    1091           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1092           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1093             :       }
    1094             : 
    1095             :     case t_PADIC: /* ty == t_QUAD */
    1096           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1097             :   }
    1098             :   /* tx < ty, !is_const_t(y) */
    1099    23374870 :   switch(ty)
    1100             :   {
    1101             :     case t_MAT:
    1102       65548 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1103       65548 :       if (isrationalzero(x)) return gcopy(y);
    1104        5474 :       return RgM_Rg_add(y, x);
    1105             :     case t_COL:
    1106      331992 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1107      331992 :       return RgC_Rg_add(y, x);
    1108             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1109      295918 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1110             :   }
    1111    22681412 :   if (is_scalar_t(tx))  {
    1112    19592202 :     if (tx == t_POLMOD)
    1113             :     {
    1114       80612 :       vx = varn(gel(x,1));
    1115       80612 :       vy = gvar(y);
    1116       80612 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1117             :       else
    1118       47306 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1119       43890 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1120             :     }
    1121    19511590 :     return add_scal(y, x, ty);
    1122             :   }
    1123             :   /* x and y are not scalars, ty != t_MAT */
    1124     3089210 :   vx = gvar(x);
    1125     3089210 :   vy = gvar(y);
    1126     3089210 :   if (vx != vy) { /* x or y is treated as a scalar */
    1127       22598 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1128       22591 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1129       22591 :                                 : add_scal(y, x, ty);
    1130             :   }
    1131             :   /* vx = vy */
    1132     3066612 :   switch(tx)
    1133             :   {
    1134             :     case t_POL:
    1135     3066206 :       switch (ty)
    1136             :       {
    1137             :         case t_SER:
    1138        4844 :           if (lg(x) == 2) return gcopy(y);
    1139        4830 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1140        4830 :           i = lg(y) + valp(y) - i;
    1141        4830 :           if (i < 3) return gcopy(y);
    1142        4823 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1143        4823 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1144        4823 :           return y;
    1145             : 
    1146     3061362 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1147             :       }
    1148           0 :       break;
    1149             : 
    1150             :     case t_SER:
    1151         406 :       if (ty == t_RFRAC)
    1152             :       {
    1153             :         GEN n, d;
    1154             :         long vn, vd;
    1155         406 :         av = avma;
    1156         406 :         n = gel(y,1); vn = gval(n, vy);
    1157         406 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1158             : 
    1159         406 :         l = lg(x) + valp(x) - (vn - vd);
    1160         406 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1161             : 
    1162             :         /* take advantage of y = t^n ! */
    1163         406 :         if (degpol(d))
    1164          70 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1165             :         else {
    1166         336 :           y = gdiv(n, gel(d,2));
    1167         336 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1168             :         }
    1169         406 :         setvalp(y, valp(y) - vd);
    1170         406 :         return gerepileupto(av, gadd(y, x));
    1171             :       }
    1172           0 :       break;
    1173             :   }
    1174           0 :   pari_err_TYPE2("+",x,y);
    1175             :   return NULL; /* LCOV_EXCL_LINE */
    1176             : }
    1177             : 
    1178             : GEN
    1179    32080641 : gaddsg(long x, GEN y)
    1180             : {
    1181    32080641 :   long ty = typ(y);
    1182             :   GEN z;
    1183             : 
    1184    32080641 :   switch(ty)
    1185             :   {
    1186    10909769 :     case t_INT:  return addsi(x,y);
    1187     8188754 :     case t_REAL: return addsr(x,y);
    1188             :     case t_INTMOD:
    1189          14 :       z = cgetg(3, t_INTMOD);
    1190          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1191     4824397 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1192     4824397 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1193     4824397 :       gel(z,2) = icopy(gel(y,2)); return z;
    1194             :     case t_COMPLEX:
    1195     7569523 :       z = cgetg(3, t_COMPLEX);
    1196     7569523 :       gel(z,1) = gaddsg(x, gel(y,1));
    1197     7569523 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1198             : 
    1199      588184 :     default: return gadd(stoi(x), y);
    1200             :   }
    1201             : }
    1202             : 
    1203             : GEN
    1204     1048818 : gsubsg(long x, GEN y)
    1205             : {
    1206             :   GEN z, a, b;
    1207             :   pari_sp av;
    1208             : 
    1209     1048818 :   switch(typ(y))
    1210             :   {
    1211      100447 :     case t_INT:  return subsi(x,y);
    1212      560244 :     case t_REAL: return subsr(x,y);
    1213             :     case t_INTMOD:
    1214          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1215          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1216       41614 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1217       41614 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1218       41614 :       gel(z,2) = icopy(gel(y,2)); return z;
    1219             :     case t_COMPLEX:
    1220      318345 :       z = cgetg(3, t_COMPLEX);
    1221      318345 :       gel(z,1) = gsubsg(x, gel(y,1));
    1222      318345 :       gel(z,2) = gneg(gel(y,2)); return z;
    1223             :   }
    1224       28119 :   av = avma;
    1225       28119 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1226             : }
    1227             : 
    1228             : /********************************************************************/
    1229             : /**                                                                **/
    1230             : /**                          SUBTRACTION                           **/
    1231             : /**                                                                **/
    1232             : /********************************************************************/
    1233             : 
    1234             : GEN
    1235  1445951289 : gsub(GEN x, GEN y)
    1236             : {
    1237  1445951289 :   long tx = typ(x), ty = typ(y);
    1238             :   pari_sp av;
    1239             :   GEN z;
    1240  1445951289 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1241             :   {
    1242  1187518277 :     case t_INT: return subii(x,y);
    1243   196749030 :     case t_REAL: return subrr(x,y);
    1244     1102103 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1245     1102103 :       z = cgetg(3,t_INTMOD);
    1246     1102103 :       if (X==Y || equalii(X,Y))
    1247     1102083 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1248          14 :       gel(z,1) = gcdii(X,Y);
    1249          14 :       warn_coercion(X,Y,gel(z,1));
    1250          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1251          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1252             :     }
    1253     4124390 :     case t_FRAC: return addsub_frac(x,y, subii);
    1254    25810491 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1255    25810491 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1256    25810491 :       if (isintzero(gel(z,2)))
    1257             :       {
    1258       80779 :         avma = (pari_sp)(z+3);
    1259       80779 :         return gsub(gel(x,1),gel(y,1));
    1260             :       }
    1261    25729712 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1262    25729712 :       return z;
    1263             :     case t_PADIC:
    1264      228855 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1265      228855 :       return addsub_pp(x,y, subii);
    1266         595 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1267         595 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1268         595 :       gel(z,1) = ZX_copy(gel(x,1));
    1269         595 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1270         595 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1271             :     case t_POLMOD:
    1272      911828 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1273      911793 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1274          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1275      282781 :     case t_FFELT: return FF_sub(x,y);
    1276             :     case t_POL: {
    1277     7650624 :       long vx = varn(x);
    1278     7650624 :       long vy = varn(y);
    1279     7650624 :       if (vx != vy) {
    1280       12496 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1281        2072 :         else                     return Rg_RgX_sub(x, y);
    1282             :       }
    1283     7638128 :       return RgX_sub(x, y);
    1284             :     }
    1285             :     case t_VEC:
    1286      271041 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1287      271041 :       return RgV_sub(x,y);
    1288             :     case t_COL:
    1289     3245382 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1290     3245382 :       return RgC_sub(x,y);
    1291             :     case t_MAT: {
    1292       11992 :       long lx = lg(x);
    1293       11992 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1294       11992 :       if (lx == 1) return cgetg(1, t_MAT);
    1295       11992 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1296       11992 :       return RgM_sub(x,y);
    1297             :     }
    1298     2095577 :     case t_RFRAC: case t_SER: break;
    1299             : 
    1300           0 :     default: pari_err_TYPE2("+",x,y);
    1301             :   }
    1302    18043900 :   av = avma;
    1303    18043900 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1304             : }
    1305             : 
    1306             : /********************************************************************/
    1307             : /**                                                                **/
    1308             : /**                        MULTIPLICATION                          **/
    1309             : /**                                                                **/
    1310             : /********************************************************************/
    1311             : static GEN
    1312      515067 : mul_ser_scal(GEN y, GEN x) {
    1313             :   long ly, i;
    1314             :   GEN z;
    1315      515067 :   if (isexactzero(x)) return gmul(Rg_get_0(y), x);
    1316      511819 :   if (ser_isexactzero(y))
    1317             :   {
    1318         441 :     if (lg(y) == 2) return gcopy(y);
    1319           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1320             :   }
    1321      511378 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1322      511378 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1323      511378 :   return normalize(z);
    1324             : }
    1325             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1326             :  * [n/d a valid RFRAC]  */
    1327             : static GEN
    1328    10442582 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1329             : {
    1330    10442582 :   pari_sp av = avma;
    1331             :   GEN z;
    1332             : 
    1333    10442582 :   switch(typ(x))
    1334             :   {
    1335             :     case t_PADIC:
    1336          14 :       n = gmul(n, x);
    1337          14 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1338          14 :       return gerepileupto(av, gdiv(n,d));
    1339             : 
    1340             :     case t_INTMOD: case t_POLMOD:
    1341          49 :       n = gmul(n, x);
    1342          49 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1343          49 :       return gerepileupto(av, gdiv(n,d));
    1344             :   }
    1345    10442519 :   z = gred_rfrac2(x, d);
    1346    10442519 :   n = simplify_shallow(n);
    1347    10442519 :   if (typ(z) == t_RFRAC)
    1348             :   {
    1349     7927719 :     n = gmul(gel(z,1), n);
    1350     7927719 :     d = gel(z,2);
    1351     7927719 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1352           0 :       z = RgX_Rg_div(n, d);
    1353             :     else
    1354     7927719 :       z = gred_rfrac_simple(n, d);
    1355             :   }
    1356             :   else
    1357     2514800 :     z = gmul(z, n);
    1358    10442519 :   return gerepileupto(av, z);
    1359             : }
    1360             : static GEN
    1361    59427169 : mul_scal(GEN y, GEN x, long ty)
    1362             : {
    1363    59427169 :   switch(ty)
    1364             :   {
    1365             :     case t_POL:
    1366    50467943 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1367    47703246 :       return RgX_Rg_mul(y, x);
    1368      380341 :     case t_SER: return mul_ser_scal(y, x);
    1369     8578843 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1370             :     case t_QFI: case t_QFR:
    1371          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1372             :   }
    1373          28 :   pari_err_TYPE2("*",x,y);
    1374             :   return NULL; /* LCOV_EXCL_LINE */
    1375             : }
    1376             : 
    1377             : static GEN
    1378      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1379             : {
    1380      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1381      160830 :   long vx = gvar(X), vy = varn(y2);
    1382      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1383        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1384             : }
    1385             : /* (x1/x2) * (y1/y2) */
    1386             : static GEN
    1387     7903558 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1388             : {
    1389             :   GEN z, X, Y;
    1390     7903558 :   pari_sp av = avma;
    1391             : 
    1392     7903558 :   X = gred_rfrac2(x1, y2);
    1393     7903558 :   Y = gred_rfrac2(y1, x2);
    1394     7903558 :   if (typ(X) == t_RFRAC)
    1395             :   {
    1396     6624630 :     if (typ(Y) == t_RFRAC) {
    1397     6559054 :       x1 = gel(X,1);
    1398     6559054 :       x2 = gel(X,2);
    1399     6559054 :       y1 = gel(Y,1);
    1400     6559054 :       y2 = gel(Y,2);
    1401     6559054 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1402             :     } else
    1403       65576 :       z = mul_gen_rfrac(Y, X);
    1404             :   }
    1405     1278928 :   else if (typ(Y) == t_RFRAC)
    1406       95254 :     z = mul_gen_rfrac(X, Y);
    1407             :   else
    1408     1183674 :     z = gmul(X, Y);
    1409     7903558 :   return gerepileupto(av, z);
    1410             : }
    1411             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1412             : static GEN
    1413      269855 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1414             : {
    1415      269855 :   pari_sp av = avma;
    1416      269855 :   GEN X = gred_rfrac2(x1, y2);
    1417      269855 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1418             :   {
    1419      263303 :     x2 = RgX_mul(gel(X,2), x2);
    1420      263303 :     x1 = gel(X,1);
    1421             :   }
    1422             :   else
    1423        6552 :     x1 = X;
    1424      269855 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1425             : }
    1426             : 
    1427             : /* Mod(y, Y) * x,  assuming x scalar */
    1428             : static GEN
    1429      927530 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1430      927530 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1431             : 
    1432             : /* cf mulqq */
    1433             : static GEN
    1434      678734 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1435             : {
    1436      678734 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1437      678734 :   pari_sp tetpil, av = avma;
    1438      678734 :   T[1] = x[1];
    1439      678734 :   p2 = gmul(gel(x,2), gel(y,2));
    1440      678734 :   p3 = gmul(gel(x,3), gel(y,3));
    1441      678734 :   p1 = gmul(gneg_i(c),p3);
    1442             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1443      678734 :   if (typ(b) == t_INT)
    1444             :   {
    1445      678720 :     if (signe(b))
    1446             :     {
    1447      395864 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1448      395864 :       if (is_pm1(b))
    1449             :       {
    1450      395206 :         if (signe(b) > 0) p3 = gneg(p3);
    1451             :       }
    1452             :       else
    1453         658 :         p3 = gmul(negi(b), p3);
    1454             :     }
    1455             :     else
    1456             :     {
    1457      282856 :       p3 = gmul(gel(x,2),gel(y,3));
    1458      282856 :       p4 = gmul(gel(x,3),gel(y,2));
    1459             :     }
    1460             :   }
    1461             :   else
    1462             :   {
    1463          14 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1464          14 :     p3 = gmul(gneg_i(b), p3);
    1465             :   }
    1466      678734 :   tetpil = avma;
    1467      678734 :   gel(T,2) = gadd(p2, p1);
    1468      678734 :   gel(T,3) = gadd(p4, p3);
    1469      678734 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1470      678734 :   return normalizepol_lg(T,4);
    1471             : }
    1472             : /* Mod(x,T) * Mod(y,T) */
    1473             : static GEN
    1474     1802636 : mul_polmod_same(GEN T, GEN x, GEN y)
    1475             : {
    1476     1802636 :   GEN z = cgetg(3,t_POLMOD), a;
    1477     1802636 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1478     1802636 :   gel(z,1) = RgX_copy(T);
    1479             :   /* x * y mod T optimised */
    1480     1802636 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1481     1278893 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1482      713628 :     a = gmul(x, y);
    1483             :   else
    1484             :   {
    1485     1089008 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1486      672770 :       a = quad_polmod_mul(T, x, y);
    1487             :     else
    1488      416238 :       a = RgXQ_mul(x, y, gel(z,1));
    1489             :   }
    1490     1802636 :   gel(z,2) = a; return z;
    1491             : }
    1492             : static GEN
    1493       33873 : sqr_polmod(GEN T, GEN x)
    1494             : {
    1495       33873 :   GEN a, z = cgetg(3,t_POLMOD);
    1496       33873 :   gel(z,1) = RgX_copy(T);
    1497       33873 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1498        7525 :     a = gsqr(x);
    1499             :   else
    1500             :   {
    1501       26348 :     pari_sp av = avma;
    1502       26348 :     a = RgXQ_sqr(x, gel(z,1));
    1503       26348 :     a = gerepileupto(av, a);
    1504             :   }
    1505       33873 :   gel(z,2) = a; return z;
    1506             : }
    1507             : /* Mod(x,X) * Mod(y,Y) */
    1508             : static GEN
    1509      141148 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1510             : {
    1511      141148 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1512      141148 :   long vx = varn(X), vy = varn(Y);
    1513      141148 :   GEN z = cgetg(3,t_POLMOD);
    1514             : 
    1515      141148 :   if (vx==vy) {
    1516             :     pari_sp av;
    1517          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1518          14 :     warn_coercion(X,Y,gel(z,1));
    1519          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1520          14 :     return z;
    1521             :   }
    1522      141134 :   if (varncmp(vx, vy) < 0)
    1523      133196 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1524             :   else
    1525        7938 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1526      141134 :   gel(z,2) = gmul(x, y); return z;
    1527             : }
    1528             : 
    1529             : #if 0 /* used by 3M only */
    1530             : /* set z = x+y and return 1 if x,y have the same sign
    1531             :  * set z = x-y and return 0 otherwise */
    1532             : static int
    1533             : did_add(GEN x, GEN y, GEN *z)
    1534             : {
    1535             :   long tx = typ(x), ty = typ(y);
    1536             :   if (tx == ty) switch(tx)
    1537             :   {
    1538             :     case t_INT: *z = addii(x,y); return 1;
    1539             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1540             :     case t_REAL:
    1541             :       if (signe(x) == -signe(y))
    1542             :       { *z = subrr(x,y); return 0; }
    1543             :       else
    1544             :       { *z = addrr(x,y); return 1; }
    1545             :   }
    1546             :   if (tx == t_REAL) switch(ty)
    1547             :   {
    1548             :     case t_INT:
    1549             :       if (signe(x) == -signe(y))
    1550             :       { *z = subri(x,y); return 0; }
    1551             :       else
    1552             :       { *z = addri(x,y); return 1; }
    1553             :     case t_FRAC:
    1554             :       if (signe(x) == -signe(gel(y,1)))
    1555             :       { *z = gsub(x,y); return 0; }
    1556             :       else
    1557             :       { *z = gadd(x,y); return 1; }
    1558             :   }
    1559             :   else if (ty == t_REAL) switch(tx)
    1560             :   {
    1561             :     case t_INT:
    1562             :       if (signe(x) == -signe(y))
    1563             :       { *z = subir(x,y); return 0; }
    1564             :       else
    1565             :       { *z = addir(x,y); return 1; }
    1566             :     case t_FRAC:
    1567             :       if (signe(gel(x,1)) == -signe(y))
    1568             :       { *z = gsub(x,y); return 0; }
    1569             :       else
    1570             :       { *z = gadd(x,y); return 1; }
    1571             :   }
    1572             :   *z = gadd(x,y); return 1;
    1573             : }
    1574             : #endif
    1575             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1576             : static GEN
    1577     4253266 : mulcIR(GEN x, GEN y)
    1578             : {
    1579     4253266 :   GEN z = cgetg(3,t_COMPLEX);
    1580     4253266 :   pari_sp av = avma;
    1581     4253266 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1582     4253266 :   gel(z,2) = gmul(y, gel(x,1));
    1583     4253266 :   return z;
    1584             : 
    1585             : }
    1586             : /* x,y COMPLEX */
    1587             : static GEN
    1588   108475739 : mulcc(GEN x, GEN y)
    1589             : {
    1590   108475739 :   GEN xr = gel(x,1), xi = gel(x,2);
    1591   108475739 :   GEN yr = gel(y,1), yi = gel(y,2);
    1592             :   GEN p1, p2, p3, p4, z;
    1593             :   pari_sp tetpil, av;
    1594             : 
    1595   108475739 :   if (isintzero(xr))
    1596             :   {
    1597     4134241 :     if (isintzero(yr)) {
    1598      604559 :       av = avma;
    1599      604559 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1600             :     }
    1601     3529682 :     return mulcIR(y, xi);
    1602             :   }
    1603   104341498 :   if (isintzero(yr)) return mulcIR(x, yi);
    1604             : 
    1605   103617914 :   z = cgetg(3,t_COMPLEX); av = avma;
    1606             : #if 0
    1607             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1608             :    * e.g. xr + xi if exponents differ */
    1609             :   if (did_add(xr, xi, &p3))
    1610             :   {
    1611             :     if (did_add(yr, yi, &p4)) {
    1612             :     /* R = xr*yr - xi*yi
    1613             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1614             :       p1 = gmul(xr,yr);
    1615             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1616             :       p3 = gmul(p3, p4);
    1617             :       p4 = gsub(p2, p1);
    1618             :     } else {
    1619             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1620             :      * I = xr*yi + xi*yr */
    1621             :       p1 = gmul(p3,p4);
    1622             :       p3 = gmul(xr,yi);
    1623             :       p4 = gmul(xi,yr);
    1624             :       p2 = gsub(p3, p4);
    1625             :     }
    1626             :   } else {
    1627             :     if (did_add(yr, yi, &p4)) {
    1628             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1629             :       * I = xr*yi +xi*yr */
    1630             :       p1 = gmul(p3,p4);
    1631             :       p3 = gmul(xr,yi);
    1632             :       p4 = gmul(xi,yr);
    1633             :       p2 = gsub(p4, p3);
    1634             :     } else {
    1635             :     /* R = xr*yr - xi*yi
    1636             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1637             :       p3 = gneg( gmul(p3, p4) );
    1638             :       p1 = gmul(xr,yr);
    1639             :       p2 = gmul(xi,yi);
    1640             :       p4 = gadd(p1, p2);
    1641             : 
    1642             :       p2 = gneg(p2);
    1643             :     }
    1644             :   }
    1645             :   tetpil = avma;
    1646             :   gel(z,1) = gadd(p1,p2);
    1647             :   gel(z,2) = gadd(p3,p4);
    1648             : #else
    1649   103617914 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1650             :   { /* 3M formula */
    1651      615927 :     p3 = addii(xr,xi);
    1652      615927 :     p4 = addii(yr,yi);
    1653      615927 :     p1 = mulii(xr,yr);
    1654      615927 :     p2 = mulii(xi,yi);
    1655      615927 :     p3 = mulii(p3,p4);
    1656      615927 :     p4 = addii(p2,p1);
    1657      615927 :     tetpil = avma;
    1658      615927 :     gel(z,1) = subii(p1,p2);
    1659      615927 :     gel(z,2) = subii(p3,p4);
    1660     1118888 :     if (!signe(gel(z,2)))
    1661      112966 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1662             :   }
    1663             :   else
    1664             :   { /* naive 4M formula: avoid all loss of accuracy */
    1665   103001987 :     p1 = gmul(xr,yr);
    1666   103001987 :     p2 = gmul(xi,yi);
    1667   103001987 :     p3 = gmul(xr,yi);
    1668   103001987 :     p4 = gmul(xi,yr);
    1669   103001987 :     tetpil = avma;
    1670   103001987 :     gel(z,1) = gsub(p1,p2);
    1671   103001987 :     gel(z,2) = gadd(p3,p4);
    1672   103001987 :     if (isintzero(gel(z,2)))
    1673             :     {
    1674       13950 :       cgiv(gel(z,2));
    1675       13950 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1676             :     }
    1677             :   }
    1678             : #endif
    1679             : 
    1680   103490998 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1681             : }
    1682             : /* x,y PADIC */
    1683             : static GEN
    1684     1774712 : mulpp(GEN x, GEN y) {
    1685     1774712 :   long l = valp(x) + valp(y);
    1686             :   pari_sp av;
    1687             :   GEN z, t;
    1688     1774712 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1689     1774713 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1690     1614027 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1691             : 
    1692     1549162 :   t = (precp(x) > precp(y))? y: x;
    1693     1549162 :   z = cgetp(t); setvalp(z,l); av = avma;
    1694     1549168 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1695     1549156 :   set_avma(av); return z;
    1696             : }
    1697             : /* x,y QUAD */
    1698             : static GEN
    1699        1316 : mulqq(GEN x, GEN y) {
    1700        1316 :   GEN z = cgetg(4,t_QUAD);
    1701        1316 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1702             :   pari_sp av, tetpil;
    1703        1316 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1704             : 
    1705        1316 :   gel(z,1) = ZX_copy(P); av = avma;
    1706        1316 :   p2 = gmul(gel(x,2),gel(y,2));
    1707        1316 :   p3 = gmul(gel(x,3),gel(y,3));
    1708        1316 :   p1 = gmul(gneg_i(c),p3);
    1709             : 
    1710        1316 :   if (signe(b))
    1711        1204 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1712             :   else
    1713             :   {
    1714         112 :     p3 = gmul(gel(x,2),gel(y,3));
    1715         112 :     p4 = gmul(gel(x,3),gel(y,2));
    1716             :   }
    1717        1316 :   tetpil = avma;
    1718        1316 :   gel(z,2) = gadd(p2,p1);
    1719        1316 :   gel(z,3) = gadd(p4,p3);
    1720        1316 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1721             : }
    1722             : 
    1723             : GEN
    1724     2282146 : mulcxI(GEN x)
    1725             : {
    1726             :   GEN z;
    1727     2282146 :   switch(typ(x))
    1728             :   {
    1729             :     case t_INT: case t_REAL: case t_FRAC:
    1730      282549 :       return mkcomplex(gen_0, x);
    1731             :     case t_COMPLEX:
    1732     1999541 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1733     1997224 :       z = cgetg(3,t_COMPLEX);
    1734     1997224 :       gel(z,1) = gneg(gel(x,2));
    1735     1997224 :       gel(z,2) = gel(x,1); return z;
    1736             :     default:
    1737          56 :       return gmul(gen_I(), x);
    1738             :   }
    1739             : }
    1740             : GEN
    1741       49161 : mulcxmI(GEN x)
    1742             : {
    1743             :   GEN z;
    1744       49161 :   switch(typ(x))
    1745             :   {
    1746             :     case t_INT: case t_REAL: case t_FRAC:
    1747        1834 :       return mkcomplex(gen_0, gneg(x));
    1748             :     case t_COMPLEX:
    1749       47264 :       if (isintzero(gel(x,1))) return gel(x,2);
    1750       46165 :       z = cgetg(3,t_COMPLEX);
    1751       46165 :       gel(z,1) = gel(x,2);
    1752       46165 :       gel(z,2) = gneg(gel(x,1)); return z;
    1753             :     default:
    1754          63 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1755             :   }
    1756             : }
    1757             : /* x * I^k */
    1758             : GEN
    1759        1477 : mulcxpowIs(GEN x, long k)
    1760             : {
    1761        1477 :   switch (k & 3)
    1762             :   {
    1763         280 :     case 1: return mulcxI(x);
    1764         483 :     case 2: return gneg(x);
    1765         364 :     case 3: return mulcxmI(x);
    1766             :   }
    1767         350 :   return x;
    1768             : }
    1769             : 
    1770             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1771             : static GEN
    1772     2232673 : fill_ser(GEN z, GEN y)
    1773             : {
    1774     2232673 :   long i, lx = lg(z), ly = lg(y);
    1775     2232673 :   if (ly >= lx) {
    1776     2215943 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1777             :   } else {
    1778       16730 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1779       16730 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1780             :   }
    1781     2232673 :   return normalize(z);
    1782             : }
    1783             : 
    1784             : GEN
    1785  3499226775 : gmul(GEN x, GEN y)
    1786             : {
    1787             :   long tx, ty, lx, ly, vx, vy, i, l;
    1788             :   pari_sp av, tetpil;
    1789             :   GEN z, p1, p2;
    1790             : 
    1791  3499226775 :   if (x == y) return gsqr(x);
    1792  2565644876 :   tx = typ(x); ty = typ(y);
    1793  2565644876 :   if (tx == ty) switch(tx)
    1794             :   {
    1795  1182677226 :     case t_INT: return mulii(x,y);
    1796   721227654 :     case t_REAL: return mulrr(x,y);
    1797     1722462 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1798     1722462 :       z = cgetg(3,t_INTMOD);
    1799     1722452 :       if (X==Y || equalii(X,Y))
    1800     1722408 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1801          14 :       gel(z,1) = gcdii(X,Y);
    1802          14 :       warn_coercion(X,Y,gel(z,1));
    1803          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1804          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1805             :     }
    1806             :     case t_FRAC:
    1807             :     {
    1808     8393214 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1809     8393214 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1810     8393214 :       z=cgetg(3,t_FRAC);
    1811     8393214 :       p1 = gcdii(x1, y2);
    1812     8393214 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1813     8393214 :       p1 = gcdii(x2, y1);
    1814     8393214 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1815     8393214 :       tetpil = avma;
    1816     8393214 :       gel(z,2) = mulii(x2,y2);
    1817     8393214 :       gel(z,1) = mulii(x1,y1);
    1818     8393214 :       fix_frac_if_int_GC(z,tetpil); return z;
    1819             :     }
    1820   105896182 :     case t_COMPLEX: return mulcc(x, y);
    1821     1602792 :     case t_PADIC: return mulpp(x, y);
    1822         994 :     case t_QUAD: return mulqq(x, y);
    1823     9577414 :     case t_FFELT: return FF_mul(x, y);
    1824             :     case t_POLMOD:
    1825     1574723 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1826     1433575 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1827      141148 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1828             :     case t_POL:
    1829    40036051 :       vx = varn(x);
    1830    40036051 :       vy = varn(y);
    1831    40036051 :       if (vx != vy) {
    1832     4777971 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1833     2062217 :         else                     return RgX_Rg_mul(y, x);
    1834             :       }
    1835    35258080 :       return RgX_mul(x, y);
    1836             : 
    1837             :     case t_SER: {
    1838     2236215 :       vx = varn(x);
    1839     2236215 :       vy = varn(y);
    1840     2236215 :       if (vx != vy) {
    1841        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1842        3675 :         else                     return mul_ser_scal(y, x);
    1843             :       }
    1844     2232540 :       lx = lg(x);
    1845     2232540 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1846     2232540 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1847     2232484 :       av = avma; z = cgetg(lx,t_SER);
    1848     2232484 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1849     2232484 :       x = ser2pol_i(x, lx);
    1850     2232484 :       y = ser2pol_i(y, lx);
    1851     2232484 :       y = RgXn_mul(x, y, lx-2);
    1852     2232484 :       return gerepilecopy(av, fill_ser(z,y));
    1853             :     }
    1854    16035591 :     case t_QFI: return qficomp(x,y);
    1855          14 :     case t_QFR: return qfrcomp(x,y);
    1856     6717049 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1857     1773528 :     case t_MAT: return RgM_mul(x, y);
    1858             : 
    1859             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1860        1533 :       z = cgetg_copy(x, &l);
    1861        1533 :       if (l != lg(y)) break;
    1862       16933 :       for (i=1; i<l; i++)
    1863             :       {
    1864       15400 :         long yi = y[i];
    1865       15400 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1866       15400 :         z[i] = x[yi];
    1867             :       }
    1868        1533 :       return z;
    1869             : 
    1870             : 
    1871             :     default:
    1872           0 :       pari_err_TYPE2("*",x,y);
    1873             :   }
    1874             :   /* tx != ty */
    1875   466172234 :   if (is_const_t(ty) && is_const_t(tx))  {
    1876   397833234 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1877   397833234 :     switch(tx) {
    1878             :     case t_INT:
    1879   356079373 :       switch(ty)
    1880             :       {
    1881   231972927 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1882             :         case t_INTMOD:
    1883      966888 :           z = cgetg(3, t_INTMOD);
    1884      966879 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1885             :         case t_FRAC:
    1886    43616312 :           if (!signe(x)) return gen_0;
    1887    27169551 :           z=cgetg(3,t_FRAC);
    1888    27169551 :           p1 = gcdii(x,gel(y,2));
    1889    27169551 :           if (equali1(p1))
    1890             :           {
    1891    14828659 :             avma = (pari_sp)z;
    1892    14828659 :             gel(z,2) = icopy(gel(y,2));
    1893    14828659 :             gel(z,1) = mulii(gel(y,1), x);
    1894             :           }
    1895             :           else
    1896             :           {
    1897    12340892 :             x = diviiexact(x,p1); tetpil = avma;
    1898    12340892 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1899    12340892 :             gel(z,1) = mulii(gel(y,1), x);
    1900    12340892 :             fix_frac_if_int_GC(z,tetpil);
    1901             :           }
    1902    27169551 :           return z;
    1903    77503903 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1904      408939 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1905        1330 :         case t_QUAD: return mulRq(x,y);
    1906     1609075 :         case t_FFELT: return FF_Z_mul(y,x);
    1907             :       }
    1908             : 
    1909             :     case t_REAL:
    1910    39796193 :       switch(ty)
    1911             :       {
    1912     9928780 :         case t_FRAC: return mulrfrac(x, y);
    1913    29867364 :         case t_COMPLEX: return mulRc(x, y);
    1914           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1915          42 :         default: pari_err_TYPE2("*",x,y);
    1916             :       }
    1917             : 
    1918             :     case t_INTMOD:
    1919         826 :       switch(ty)
    1920             :       {
    1921         371 :         case t_FRAC: { GEN X = gel(x,1);
    1922         371 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1923         371 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1924             :         }
    1925          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1926         301 :         case t_PADIC: { GEN X = gel(x,1);
    1927         301 :           z = cgetg(3, t_INTMOD);
    1928         301 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1929             :         }
    1930          42 :         case t_QUAD: return mulRq(x, y);
    1931             :         case t_FFELT:
    1932          63 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1933           0 :             pari_err_OP("*",x,y);
    1934          63 :           return FF_Z_mul(y,gel(x,2));
    1935             :       }
    1936             : 
    1937             :     case t_FRAC:
    1938     1956778 :       switch(ty)
    1939             :       {
    1940     1924326 :         case t_COMPLEX: return mulRc(x, y);
    1941       30163 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1942         308 :         case t_QUAD:  return mulRq(x, y);
    1943        1981 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1944             :       }
    1945             : 
    1946             :     case t_FFELT:
    1947          35 :       pari_err_TYPE2("*",x,y);
    1948             : 
    1949             :     case t_COMPLEX:
    1950          21 :       switch(ty)
    1951             :       {
    1952             :         case t_PADIC:
    1953          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    1954             :         case t_QUAD:
    1955           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    1956           7 :           return mulqf(y, x, lx);
    1957             :       }
    1958             : 
    1959             :     case t_PADIC: /* ty == t_QUAD */
    1960          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    1961             :     }
    1962             :   }
    1963             : 
    1964    68339071 :   if (is_matvec_t(ty))
    1965             :   {
    1966     5406254 :     if (!is_matvec_t(tx))
    1967             :     {
    1968     5257232 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    1969     5257232 :       z = cgetg_copy(y, &ly);
    1970     5257232 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1971     5257218 :       return z;
    1972             :     }
    1973      149022 :     switch(tx)
    1974             :     {
    1975             :       case t_VEC:
    1976       41136 :         switch(ty) {
    1977       14861 :           case t_COL: return RgV_RgC_mul(x,y);
    1978       26275 :           case t_MAT: return RgV_RgM_mul(x,y);
    1979             :         }
    1980           0 :         break;
    1981             :       case t_COL:
    1982        1687 :         switch(ty) {
    1983        1687 :           case t_VEC: return RgC_RgV_mul(x,y);
    1984           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    1985             :         }
    1986           0 :         break;
    1987             :       case t_MAT:
    1988      106199 :         switch(ty) {
    1989           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    1990      106199 :           case t_COL: return RgM_RgC_mul(x,y);
    1991             :         }
    1992             :     }
    1993             :   }
    1994    62932844 :   if (is_matvec_t(tx))
    1995             :   {
    1996      367536 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    1997      367534 :     z = cgetg_copy(x, &lx);
    1998      367533 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    1999      367534 :     return z;
    2000             :   }
    2001    62565308 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2002             :   /* tx < ty, !ismatvec(x and y) */
    2003             : 
    2004    62565308 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2005      889737 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2006    61675571 :   if (is_scalar_t(tx)) {
    2007    58258708 :     if (tx == t_POLMOD) {
    2008      637222 :       vx = varn(gel(x,1));
    2009      637222 :       vy = gvar(y);
    2010      637222 :       if (vx != vy) {
    2011      268532 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2012       37793 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2013             :       }
    2014             :       /* error if ty == t_SER */
    2015      368690 :       av = avma; y = gmod(y, gel(x,1));
    2016      368683 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2017             :     }
    2018    57621486 :     return mul_scal(y, x, ty);
    2019             :   }
    2020             : 
    2021             :   /* x and y are not scalars, nor matvec */
    2022     3416863 :   vx = gvar(x);
    2023     3416863 :   vy = gvar(y);
    2024     3416863 :   if (vx != vy) /* x or y is treated as a scalar */
    2025     2781055 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2026     2781055 :                                 : mul_scal(y, x, ty);
    2027             :   /* vx = vy */
    2028     1997009 :   switch(tx)
    2029             :   {
    2030             :     case t_POL:
    2031     1996981 :       switch (ty)
    2032             :       {
    2033             :         case t_SER:
    2034             :         {
    2035             :           long vn;
    2036      133270 :           if (lg(x) == 2) return zeropol(vx);
    2037      133270 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2038      133039 :           av = avma;
    2039      133039 :           vn = RgX_valrem(x, &x);
    2040             :           /* take advantage of x = t^n ! */
    2041      133039 :           if (degpol(x)) {
    2042        1988 :             p1 = RgX_to_ser(x,lg(y));
    2043        1988 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2044        1988 :             p2 = gmul(p1,y);
    2045        1988 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2046             :           } else {
    2047      131051 :             set_avma(av);
    2048      131051 :             p2 = mul_ser_scal(y, gel(x,2));
    2049             :           }
    2050      133039 :           setvalp(p2, valp(p2) + vn);
    2051      133039 :           return p2;
    2052             :         }
    2053             : 
    2054     1863711 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2055             :       }
    2056           0 :       break;
    2057             : 
    2058             :     case t_SER:
    2059          28 :       switch (ty)
    2060             :       {
    2061             :         case t_RFRAC:
    2062          28 :           av = avma;
    2063          28 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2064             :       }
    2065           0 :       break;
    2066             :   }
    2067           0 :   pari_err_TYPE2("*",x,y);
    2068             :   return NULL; /* LCOV_EXCL_LINE */
    2069             : }
    2070             : 
    2071             : /* return a non-normalized result */
    2072             : GEN
    2073       88665 : sqr_ser_part(GEN x, long l1, long l2)
    2074             : {
    2075             :   long i, j, l;
    2076             :   pari_sp av;
    2077             :   GEN Z, z, p1, p2;
    2078             :   long mi;
    2079       88665 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2080       88658 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2081       88658 :   Z = cgetg(l2-l1+3, t_SER);
    2082       88658 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2083       88658 :   z = Z + 2-l1;
    2084       88658 :   x += 2; mi = 0;
    2085      340697 :   for (i=0; i<l1; i++)
    2086             :   {
    2087      252039 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2088             :   }
    2089             : 
    2090      608086 :   for (i=l1; i<=l2; i++)
    2091             :   {
    2092      519428 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2093      519428 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2094     1063858 :     for (j=i-mi; j<=minss(l,mi); j++)
    2095      544430 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2096      519428 :     p1 = gshift(p1,1);
    2097      519428 :     if ((i&1) == 0 && p2[i>>1])
    2098       50423 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2099      519428 :     gel(z,i) = gerepileupto(av,p1);
    2100             :   }
    2101       88658 :   return Z;
    2102             : }
    2103             : 
    2104             : GEN
    2105  1014935453 : gsqr(GEN x)
    2106             : {
    2107             :   long i, lx;
    2108             :   pari_sp av, tetpil;
    2109             :   GEN z, p1, p2, p3, p4;
    2110             : 
    2111  1014935453 :   switch(typ(x))
    2112             :   {
    2113   943450902 :     case t_INT: return sqri(x);
    2114    25614075 :     case t_REAL: return sqrr(x);
    2115      141999 :     case t_INTMOD: { GEN X = gel(x,1);
    2116      141999 :       z = cgetg(3,t_INTMOD);
    2117      141984 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2118      141987 :       gel(z,1) = icopy(X); return z;
    2119             :     }
    2120     1284280 :     case t_FRAC: return sqrfrac(x);
    2121             : 
    2122             :     case t_COMPLEX:
    2123     3097925 :       if (isintzero(gel(x,1))) {
    2124       30508 :         av = avma;
    2125       30508 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2126             :       }
    2127     3067417 :       z = cgetg(3,t_COMPLEX); av = avma;
    2128     3067417 :       p1 = gadd(gel(x,1),gel(x,2));
    2129     3067417 :       p2 = gsub(gel(x,1), gel(x,2));
    2130     3067417 :       p3 = gmul(gel(x,1),gel(x,2));
    2131     3067417 :       tetpil = avma;
    2132     3067417 :       gel(z,1) = gmul(p1,p2);
    2133     3067417 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2134             : 
    2135             :     case t_PADIC:
    2136       28063 :       z = cgetg(5,t_PADIC);
    2137       28063 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2138       28063 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2139       28063 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2140       28063 :       gel(z,2) = icopy(gel(x,2));
    2141       28063 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2142       28063 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2143       28063 :       return z;
    2144             : 
    2145          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2146          28 :       p1 = gel(x,1);
    2147          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2148          28 :       p2 = gsqr(gel(x,2));
    2149          28 :       p3 = gsqr(gel(x,3));
    2150          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2151             : 
    2152          28 :       if (gequal0(gel(p1,3)))
    2153             :       {
    2154           7 :         tetpil = avma;
    2155           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2156           7 :         av = avma;
    2157           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2158           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2159             :       }
    2160             : 
    2161          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2162          21 :       tetpil = avma;
    2163          21 :       gel(z,2) = gadd(p2,p4);
    2164          21 :       gel(z,3) = gadd(p1,p3);
    2165          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2166             : 
    2167             :     case t_POLMOD:
    2168       33873 :       return sqr_polmod(gel(x,1), gel(x,2));
    2169             : 
    2170     2892109 :     case t_FFELT: return FF_sqr(x);
    2171             : 
    2172      101682 :     case t_POL: return RgX_sqr(x);
    2173             : 
    2174             :     case t_SER:
    2175       24108 :       lx = lg(x);
    2176       24108 :       if (ser_isexactzero(x)) {
    2177           7 :         GEN z = gcopy(x);
    2178           7 :         setvalp(z, 2*valp(x));
    2179           7 :         return z;
    2180             :       }
    2181       24101 :       if (lx < 40)
    2182       23912 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2183             :       else
    2184             :       {
    2185         189 :         pari_sp av = avma;
    2186         189 :         GEN z = cgetg(lx,t_SER);
    2187         189 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2188         189 :         x = ser2pol_i(x,lx);
    2189         189 :         x = RgXn_sqr(x, lx-2);
    2190         189 :         return gerepilecopy(av, fill_ser(z,x));
    2191             :       }
    2192             : 
    2193           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2194           7 :       gel(z,1) = gsqr(gel(x,1));
    2195           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2196             : 
    2197         980 :     case t_MAT: return RgM_sqr(x);
    2198          14 :     case t_QFR: return qfrsqr(x);
    2199    38264792 :     case t_QFI: return qfisqr(x);
    2200             :     case t_VECSMALL:
    2201         616 :       z = cgetg_copy(x, &lx);
    2202        6720 :       for (i=1; i<lx; i++)
    2203             :       {
    2204        6104 :         long xi = x[i];
    2205        6104 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2206        6104 :         z[i] = x[xi];
    2207             :       }
    2208         616 :       return z;
    2209             :   }
    2210           0 :   pari_err_TYPE2("*",x,x);
    2211             :   return NULL; /* LCOV_EXCL_LINE */
    2212             : }
    2213             : 
    2214             : /********************************************************************/
    2215             : /**                                                                **/
    2216             : /**                           DIVISION                             **/
    2217             : /**                                                                **/
    2218             : /********************************************************************/
    2219             : static GEN
    2220      271022 : div_rfrac_scal(GEN x, GEN y)
    2221             : {
    2222      271022 :   pari_sp av = avma;
    2223      271022 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2224      271022 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2225             : }
    2226             : static GEN
    2227       37445 : div_scal_rfrac(GEN x, GEN y)
    2228             : {
    2229       37445 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2230       37445 :   pari_sp av = avma;
    2231       37445 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2232             :   {
    2233           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2234           0 :     y1 = gel(y1,2);
    2235             :   }
    2236       37445 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2237             : }
    2238             : static GEN
    2239     1186509 : div_rfrac(GEN x, GEN y)
    2240     1186509 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2241             : 
    2242             : static GEN
    2243      622137 : div_ser_scal(GEN x, GEN y) {
    2244             :   long i, lx;
    2245             :   GEN z;
    2246      622137 :   if (ser_isexactzero(x))
    2247             :   {
    2248           7 :     if (lg(x) == 2) return gcopy(x);
    2249           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2250             :   }
    2251      622130 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2252      622130 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2253      622130 :   return normalize(z);
    2254             : }
    2255             : GEN
    2256        7525 : ser_normalize(GEN x)
    2257             : {
    2258        7525 :   long i, lx = lg(x);
    2259             :   GEN c, z;
    2260        7525 :   if (lx == 2) return x;
    2261        7525 :   c = gel(x,2); if (gequal1(c)) return x;
    2262        7448 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2263        7448 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2264        7448 :   return z;
    2265             : }
    2266             : 
    2267             : static GEN
    2268     1482136 : div_T_scal(GEN x, GEN y, long tx) {
    2269     1482136 :   switch(tx)
    2270             :   {
    2271      592841 :     case t_POL: return RgX_Rg_div(x, y);
    2272      622130 :     case t_SER: return div_ser_scal(x, y);
    2273      267165 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2274             :   }
    2275           0 :   pari_err_TYPE2("/",x,y);
    2276             :   return NULL; /* LCOV_EXCL_LINE */
    2277             : }
    2278             : 
    2279             : static GEN
    2280     9216100 : div_scal_pol(GEN x, GEN y) {
    2281     9216100 :   long ly = lg(y);
    2282             :   pari_sp av;
    2283     9216100 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2284     9138133 :   if (isrationalzero(x)) return zeropol(varn(y));
    2285     7090930 :   av = avma;
    2286     7090930 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2287             : }
    2288             : static GEN
    2289       32592 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2290             :   GEN z;
    2291             :   long ly, i;
    2292       32592 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2293       32592 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2294       32592 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2295       32592 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2296       32592 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2297       32592 :   y = gdiv(z,y); pari_free(z); return y;
    2298             : }
    2299             : static GEN
    2300     9238950 : div_scal_T(GEN x, GEN y, long ty) {
    2301     9238950 :   switch(ty)
    2302             :   {
    2303     9169669 :     case t_POL: return div_scal_pol(x, y);
    2304       32592 :     case t_SER: return div_scal_ser(x, y);
    2305       36689 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2306             :   }
    2307           0 :   pari_err_TYPE2("/",x,y);
    2308             :   return NULL; /* LCOV_EXCL_LINE */
    2309             : }
    2310             : 
    2311             : /* assume tx = ty = t_SER, same variable vx */
    2312             : static GEN
    2313      876206 : div_ser(GEN x, GEN y, long vx)
    2314             : {
    2315      876206 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2316             :   GEN y_lead, p1, p2, z;
    2317             : 
    2318      876206 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2319      876199 :   if (ser_isexactzero(x))
    2320             :   {
    2321       59843 :     if (lx == 2) return zeroser(vx, l);
    2322         154 :     return scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), l);
    2323             :   }
    2324      816356 :   y_lead = gel(y,2);
    2325      816356 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2326             :   {
    2327          14 :     GEN y0 = y;
    2328          14 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2329          35 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2330             :     {
    2331          28 :       y_lead = gel(y,2);
    2332          28 :       if (!gequal0(y_lead)) break;
    2333             :     }
    2334          14 :     if (ly <= 2) pari_err_INV("div_ser", y0);
    2335             :   }
    2336      816349 :   if (ly < lx) lx = ly;
    2337      816349 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2338     2971431 :   for (i=3; i<lx; i++)
    2339             :   {
    2340     2155082 :     p1 = gel(y,i);
    2341     2155082 :     if (isrationalzero(p1)) p1 = NULL;
    2342     2155082 :     gel(p2,i) = p1;
    2343             :   }
    2344      816349 :   z = cgetg(lx,t_SER);
    2345      816349 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2346      816349 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2347     2971431 :   for (i=3; i<lx; i++)
    2348             :   {
    2349     2155082 :     pari_sp av = avma;
    2350     2155082 :     p1 = gel(x,i);
    2351    12823052 :     for (j=2, l=i; j<i; j++, l--)
    2352    10667970 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2353     2155082 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2354             :   }
    2355      816349 :   return normalize(z);
    2356             : }
    2357             : /* x,y compatible PADIC */
    2358             : static GEN
    2359      224562 : divpp(GEN x, GEN y) {
    2360             :   pari_sp av;
    2361             :   long a, b;
    2362             :   GEN z, M;
    2363             : 
    2364      224562 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2365      224556 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2366      223898 :   a = precp(x);
    2367      223898 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2368      223898 :   z = cgetg(5, t_PADIC);
    2369      223922 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2370      223921 :   gel(z,2) = icopy(gel(x,2));
    2371      223925 :   gel(z,3) = icopy(M); av = avma;
    2372      223928 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2373      223927 :   return z;
    2374             : }
    2375             : static GEN
    2376       29462 : div_polmod_same(GEN T, GEN x, GEN y)
    2377             : {
    2378       29462 :   long v = varn(T);
    2379       29462 :   GEN a, z = cgetg(3, t_POLMOD);
    2380       29462 :   gel(z,1) = RgX_copy(T);
    2381       29462 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2382       10864 :     a = gdiv(x, y);
    2383       18598 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2384         266 :   {
    2385         266 :     pari_sp av = avma;
    2386         266 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2387             :   }
    2388       18332 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2389        5964 :   {
    2390        5964 :     pari_sp av = avma;
    2391        5964 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2392        5964 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2393        5964 :     a = gerepileupto(av, a);
    2394             :   }
    2395             :   else
    2396             :   {
    2397       12368 :     pari_sp av = avma;
    2398       12368 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2399       12368 :     a = gerepileupto(av, a);
    2400             :   }
    2401       29462 :   gel(z,2) = a; return z;
    2402             : }
    2403             : GEN
    2404   171380734 : gdiv(GEN x, GEN y)
    2405             : {
    2406   171380734 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2407             :   pari_sp av, tetpil;
    2408             :   GEN z, p1, p2;
    2409             : 
    2410   171380734 :   if (tx == ty) switch(tx)
    2411             :   {
    2412             :     case t_INT:
    2413    65625335 :       return Qdivii(x,y);
    2414             : 
    2415    49762559 :     case t_REAL: return divrr(x,y);
    2416       37798 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2417       37798 :       z = cgetg(3,t_INTMOD);
    2418       37799 :       if (X==Y || equalii(X,Y))
    2419       37771 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2420          28 :       gel(z,1) = gcdii(X,Y);
    2421          28 :       warn_coercion(X,Y,gel(z,1));
    2422          28 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2423          28 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2424             :     }
    2425             :     case t_FRAC: {
    2426     1636296 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2427     1636296 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2428     1636296 :       z = cgetg(3, t_FRAC);
    2429     1636296 :       p1 = gcdii(x1, y1);
    2430     1636296 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2431     1636296 :       p1 = gcdii(x2, y2);
    2432     1636296 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2433     1636296 :       tetpil = avma;
    2434     1636296 :       gel(z,2) = mulii(x2,y1);
    2435     1636296 :       gel(z,1) = mulii(x1,y2);
    2436     1636296 :       normalize_frac(z);
    2437     1636296 :       fix_frac_if_int_GC(z,tetpil);
    2438     1636296 :       return z;
    2439             :     }
    2440             :     case t_COMPLEX:
    2441     2637524 :       if (isintzero(gel(y,1)))
    2442             :       {
    2443       57967 :         y = gel(y,2);
    2444       57967 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2445       40768 :         z = cgetg(3,t_COMPLEX);
    2446       40768 :         gel(z,1) = gdiv(gel(x,2), y);
    2447       40768 :         av = avma;
    2448       40768 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2449       40768 :         return z;
    2450             :       }
    2451     2579557 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2452     2579557 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2453             : 
    2454             :     case t_PADIC:
    2455      122137 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2456      122137 :       return divpp(x, y);
    2457             : 
    2458             :     case t_QUAD:
    2459         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2460         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2461         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2462             : 
    2463      234468 :     case t_FFELT: return FF_div(x,y);
    2464             : 
    2465             :     case t_POLMOD:
    2466       31093 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2467       29462 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2468             :       else {
    2469        1631 :         av = avma;
    2470        1631 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2471             :       }
    2472       31093 :       return z;
    2473             : 
    2474             :     case t_POL:
    2475    18514006 :       vx = varn(x);
    2476    18514006 :       vy = varn(y);
    2477    18514006 :       if (vx != vy) {
    2478       95088 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2479       46431 :                             else return div_scal_pol(x, y);
    2480             :       }
    2481    18418918 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2482    18418918 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2483    18229858 :       av = avma;
    2484    18229858 :       return gerepileupto(av, gred_rfrac2(x,y));
    2485             : 
    2486             :     case t_SER:
    2487      430474 :       vx = varn(x);
    2488      430474 :       vy = varn(y);
    2489      430474 :       if (vx != vy) {
    2490           7 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2491           0 :                             else return div_scal_ser(x, y);
    2492             :       }
    2493      430467 :       return div_ser(x, y, vx);
    2494             :     case t_RFRAC:
    2495     1190989 :       vx = varn(gel(x,2));
    2496     1190989 :       vy = varn(gel(y,2));
    2497     1190989 :       if (vx != vy) {
    2498        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2499         756 :                             else return div_scal_rfrac(x, y);
    2500             :       }
    2501     1186509 :       return div_rfrac(x,y);
    2502             : 
    2503           7 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2504           7 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2505             : 
    2506             :     case t_MAT:
    2507           7 :       av = avma; p1 = RgM_inv(y);
    2508           7 :       if (!p1) pari_err_INV("gdiv",y);
    2509           7 :       return gerepileupto(av, RgM_mul(x, p1));
    2510             : 
    2511           0 :     default: pari_err_TYPE2("/",x,y);
    2512             :   }
    2513             : 
    2514    31157712 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2515             :   {
    2516     2682265 :     long s = signe(x);
    2517     2682265 :     if (!s) {
    2518      688279 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2519      688279 :       switch (ty)
    2520             :       {
    2521      685094 :       default: return gen_0;
    2522             :       case t_INTMOD:
    2523          56 :         z = cgetg(3,t_INTMOD);
    2524          56 :         gel(z,1) = icopy(gel(y,1));
    2525          56 :         gel(z,2) = gen_0; return z;
    2526        3129 :       case t_FFELT: return FF_zero(y);
    2527             :       }
    2528             :     }
    2529     1993986 :     if (is_pm1(x)) {
    2530      540214 :       if (s > 0) return ginv(y);
    2531       35602 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2532             :     }
    2533     1453772 :     switch(ty)
    2534             :     {
    2535      149748 :       case t_REAL: return divir(x,y);
    2536             :       case t_INTMOD:
    2537          21 :         z = cgetg(3, t_INTMOD);
    2538          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2539             :       case t_FRAC:
    2540     1186481 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2541     1186481 :         if (equali1(p1))
    2542             :         {
    2543      809342 :           avma = (pari_sp)z;
    2544      809342 :           gel(z,2) = icopy(gel(y,1));
    2545      809342 :           gel(z,1) = mulii(gel(y,2), x);
    2546      809342 :           normalize_frac(z);
    2547      809342 :           fix_frac_if_int(z);
    2548             :         }
    2549             :         else
    2550             :         {
    2551      377139 :           x = diviiexact(x,p1); tetpil = avma;
    2552      377139 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2553      377139 :           gel(z,1) = mulii(gel(y,2), x);
    2554      377139 :           normalize_frac(z);
    2555      377139 :           fix_frac_if_int_GC(z,tetpil);
    2556             :         }
    2557     1186481 :         return z;
    2558             : 
    2559         224 :       case t_FFELT: return Z_FF_div(x,y);
    2560      105055 :       case t_COMPLEX: return divRc(x,y);
    2561       12201 :       case t_PADIC: return divTp(x, y);
    2562             :       case t_QUAD:
    2563          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2564          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2565             :     }
    2566             :   }
    2567    28475447 :   if (gequal0(y))
    2568             :   {
    2569          42 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2570          21 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2571             :   }
    2572             : 
    2573    28475419 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2574             :   {
    2575             :     case t_REAL:
    2576     9979513 :       switch(ty)
    2577             :       {
    2578     8729946 :         case t_INT: return divri(x,y);
    2579             :         case t_FRAC:
    2580      791207 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2581      791207 :           return gerepileuptoleaf(av, z);
    2582      458332 :         case t_COMPLEX: return divRc(x, y);
    2583           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2584          21 :         default: pari_err_TYPE2("/",x,y);
    2585             :       }
    2586             : 
    2587             :     case t_INTMOD:
    2588       13482 :       switch(ty)
    2589             :       {
    2590             :         case t_INT:
    2591       13405 :           z = cgetg(3, t_INTMOD);
    2592       13405 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2593          28 :         case t_FRAC: { GEN X = gel(x,1);
    2594          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2595          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2596             :         }
    2597             :         case t_FFELT:
    2598           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2599           0 :             pari_err_OP("/",x,y);
    2600           7 :           return Z_FF_div(gel(x,2),y);
    2601             : 
    2602             :         case t_COMPLEX:
    2603           0 :           av = avma;
    2604           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2605             : 
    2606             :         case t_QUAD:
    2607           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2608           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2609             : 
    2610           7 :         case t_PADIC: { GEN X = gel(x,1);
    2611           7 :           z = cgetg(3, t_INTMOD);
    2612           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2613             :         }
    2614          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2615             :       }
    2616             : 
    2617             :     case t_FRAC:
    2618      556752 :       switch(ty)
    2619             :       {
    2620      513988 :         case t_INT: z = cgetg(3, t_FRAC);
    2621      513988 :         p1 = gcdii(y,gel(x,1));
    2622      513988 :         if (equali1(p1))
    2623             :         {
    2624      338068 :           avma = (pari_sp)z; tetpil = 0;
    2625      338068 :           gel(z,1) = icopy(gel(x,1));
    2626             :         }
    2627             :         else
    2628             :         {
    2629      175920 :           y = diviiexact(y,p1); tetpil = avma;
    2630      175920 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2631             :         }
    2632      513988 :         gel(z,2) = mulii(gel(x,2),y);
    2633      513988 :         normalize_frac(z);
    2634      513988 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2635      513988 :         return z;
    2636             : 
    2637             :         case t_REAL:
    2638       40384 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2639       40384 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2640             : 
    2641           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2642           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2643           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2644             :         }
    2645             : 
    2646          28 :         case t_FFELT: av=avma;
    2647          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2648             : 
    2649         427 :         case t_COMPLEX: return divRc(x, y);
    2650             : 
    2651             :         case t_PADIC:
    2652        1911 :           if (!signe(gel(x,1))) return gen_0;
    2653        1911 :           return divTp(x, y);
    2654             : 
    2655             :         case t_QUAD:
    2656           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2657           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2658             :       }
    2659             : 
    2660             :     case t_FFELT:
    2661        1841 :       switch (ty)
    2662             :       {
    2663        1764 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2664          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2665             :         case t_INTMOD:
    2666           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2667           0 :             pari_err_OP("/",x,y);
    2668           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2669             :         default:
    2670          42 :         pari_err_TYPE2("/",x,y);
    2671             :       }
    2672           0 :       break;
    2673             : 
    2674             :     case t_COMPLEX:
    2675     4425557 :       switch(ty)
    2676             :       {
    2677     4425557 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2678           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2679             :         case t_PADIC:
    2680           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2681             :         case t_QUAD:
    2682           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2683           0 :           return divfq(x, y, lx);
    2684             :       }
    2685             : 
    2686             :     case t_PADIC:
    2687        1988 :       switch(ty)
    2688             :       {
    2689        1939 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2690        1939 :           return signe(gel(x,4))? divpT(x, y)
    2691        1939 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2692             :         }
    2693           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2694           7 :           z = cgetg(3, t_INTMOD);
    2695           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2696             :         }
    2697             :         case t_COMPLEX: case t_QUAD:
    2698           7 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2699           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2700             : 
    2701          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2702             :       }
    2703             : 
    2704             :     case t_QUAD:
    2705         896 :       switch (ty)
    2706             :       {
    2707             :         case t_INT: case t_INTMOD: case t_FRAC:
    2708         847 :           z = cgetg(4,t_QUAD);
    2709         847 :           gel(z,1) = ZX_copy(gel(x,1));
    2710         847 :           gel(z,2) = gdiv(gel(x,2), y);
    2711         847 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2712          28 :         case t_REAL: return divqf(x, y, lg(y));
    2713           7 :         case t_PADIC: return divTp(x, y);
    2714             :         case t_COMPLEX:
    2715           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2716           0 :           return divqf(x, y, ly);
    2717             :       }
    2718             :   }
    2719    13495411 :   switch(ty) {
    2720             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2721       31330 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2722             :     case t_MAT:
    2723          42 :       av = avma; p1 = RgM_inv(y);
    2724          35 :       if (!p1) pari_err_INV("gdiv",y);
    2725          28 :       return gerepileupto(av, gmul(x, p1));
    2726             :     case t_VEC: case t_COL:
    2727             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2728           0 :       pari_err_TYPE2("/",x,y);
    2729             :   }
    2730    13464039 :   switch(tx) {
    2731             :     case t_VEC: case t_COL: case t_MAT:
    2732      904587 :       z = cgetg_copy(x, &lx);
    2733      904587 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2734      904587 :       return z;
    2735             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2736           0 :       pari_err_TYPE2("/",x,y);
    2737             :   }
    2738             : 
    2739    12559452 :   vy = gvar(y);
    2740    12559452 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2741       76615 :     vx = varn(X);
    2742       76615 :     if (vx != vy) {
    2743       76230 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2744       76216 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2745             :     }
    2746             :     /* y is POL, SER or RFRAC */
    2747         385 :     av = avma;
    2748         385 :     switch(ty)
    2749             :     {
    2750           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2751         385 :       default: y = ginvmod(gmod(y,X), X);
    2752             :     }
    2753         378 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2754             :   }
    2755             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2756             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2757             :    * variable NO_VARIABLE, then the operation is incorrect */
    2758    12482837 :   vx = gvar(x);
    2759    12482837 :   if (vx != vy) { /* includes cases where one is scalar */
    2760    10721072 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2761     9238936 :                         else return div_scal_T(x, y, ty);
    2762             :   }
    2763     1761765 :   switch(tx)
    2764             :   {
    2765             :     case t_POL:
    2766     1046108 :       switch(ty)
    2767             :       {
    2768             :         case t_SER:
    2769          70 :           if (lg(y) == 2)
    2770           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2771          70 :           p1 = RgX_to_ser(x,lg(y));
    2772          70 :           p2 = div_ser(p1, y, vx);
    2773          63 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2774          63 :           return p2;
    2775             : 
    2776             :         case t_RFRAC:
    2777             :         {
    2778     1046038 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2779     1046038 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2780           0 :             return mul_rfrac_scal(y2, y1, x);
    2781     1046038 :           av = avma;
    2782     1046038 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2783             :         }
    2784             :       }
    2785           0 :       break;
    2786             : 
    2787             :     case t_SER:
    2788      445788 :       switch(ty)
    2789             :       {
    2790             :         case t_POL:
    2791      445774 :           if (lg(x) == 2)
    2792         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2793      445669 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2794      445669 :           p2 = div_ser(x, p1, vx);
    2795      445669 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2796      445669 :           return p2;
    2797             :         case t_RFRAC:
    2798          14 :           av = avma;
    2799          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2800             :       }
    2801           0 :       break;
    2802             : 
    2803             :     case t_RFRAC:
    2804      269855 :       switch(ty)
    2805             :       {
    2806      269855 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2807             :         case t_SER:
    2808           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2809           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2810             :       }
    2811           0 :       break;
    2812             :   }
    2813          14 :   pari_err_TYPE2("/",x,y);
    2814             :   return NULL; /* LCOV_EXCL_LINE */
    2815             : }
    2816             : 
    2817             : /********************************************************************/
    2818             : /**                                                                **/
    2819             : /**                     SIMPLE MULTIPLICATION                      **/
    2820             : /**                                                                **/
    2821             : /********************************************************************/
    2822             : GEN
    2823   104044253 : gmulsg(long s, GEN y)
    2824             : {
    2825             :   long ly, i;
    2826             :   pari_sp av;
    2827             :   GEN z;
    2828             : 
    2829   104044253 :   switch(typ(y))
    2830             :   {
    2831    68789545 :     case t_INT:  return mulsi(s,y);
    2832    26373703 :     case t_REAL: return mulsr(s,y);
    2833      172290 :     case t_INTMOD: { GEN p = gel(y,1);
    2834      172290 :       z = cgetg(3,t_INTMOD);
    2835      172287 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2836      172298 :       gel(z,1) = icopy(p); return z;
    2837             :     }
    2838      627124 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2839             :     case t_FRAC:
    2840     1016375 :       if (!s) return gen_0;
    2841     1011972 :       z = cgetg(3,t_FRAC);
    2842     1011972 :       i = labs(s); i = ugcd(i, umodiu(gel(y,2), i));
    2843     1011972 :       if (i == 1)
    2844             :       {
    2845      536805 :         gel(z,2) = icopy(gel(y,2));
    2846      536805 :         gel(z,1) = mulis(gel(y,1), s);
    2847             :       }
    2848             :       else
    2849             :       {
    2850      475167 :         gel(z,2) = divis(gel(y,2), i);
    2851      475167 :         gel(z,1) = mulis(gel(y,1), s/i);
    2852      475167 :         fix_frac_if_int(z);
    2853             :       }
    2854     1011972 :       return z;
    2855             : 
    2856     5922563 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2857     5922563 :       gel(z,1) = gmulsg(s,gel(y,1));
    2858     5922563 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2859             : 
    2860             :     case t_PADIC:
    2861        9212 :       if (!s) return gen_0;
    2862        9212 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2863             : 
    2864           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2865           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2866           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2867           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2868             : 
    2869             :     case t_POLMOD:
    2870       82922 :       retmkpolmod(gmulsg(s,gel(y,2)), RgX_copy(gel(y,1)));
    2871             : 
    2872             :     case t_POL:
    2873      665906 :       if (!signe(y)) return RgX_copy(y);
    2874      648091 :       if (!s) return scalarpol(Rg_get_0(y), varn(y));
    2875      585658 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2876      585668 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2877      585662 :       return normalizepol_lg(z, ly);
    2878             : 
    2879             :     case t_SER:
    2880         182 :       if (ser_isexactzero(y)) return gcopy(y);
    2881         182 :       if (!s) return Rg_get_0(y);
    2882         182 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2883         182 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2884         182 :       return normalize(z);
    2885             : 
    2886             :     case t_RFRAC:
    2887          14 :       if (!s) return zeropol(varn(gel(y,2)));
    2888          14 :       if (s == 1) return gcopy(y);
    2889          14 :       if (s == -1) return gneg(y);
    2890          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    2891             : 
    2892             :     case t_VEC: case t_COL: case t_MAT:
    2893      384410 :       z = cgetg_copy(y, &ly);
    2894      384410 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2895      384410 :       return z;
    2896             :   }
    2897           0 :   pari_err_TYPE("gmulsg",y);
    2898             :   return NULL; /* LCOV_EXCL_LINE */
    2899             : }
    2900             : 
    2901             : /********************************************************************/
    2902             : /**                                                                **/
    2903             : /**                       SIMPLE DIVISION                          **/
    2904             : /**                                                                **/
    2905             : /********************************************************************/
    2906             : 
    2907             : GEN
    2908    17919116 : gdivgs(GEN x, long s)
    2909             : {
    2910    17919116 :   long tx = typ(x), lx, i;
    2911             :   pari_sp av;
    2912             :   GEN z, y;
    2913             : 
    2914    17919116 :   if (!s)
    2915             :   {
    2916           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2917           0 :     pari_err_INV("gdivgs",gen_0);
    2918             :   }
    2919    17919117 :   switch(tx)
    2920             :   {
    2921             :     case t_INT:
    2922     6864173 :       av = avma; z = divis_rem(x,s,&i);
    2923     6864173 :       if (!i) return z;
    2924             : 
    2925     4808665 :       i = cgcd(s, i);
    2926     4808665 :       set_avma(av); z = cgetg(3,t_FRAC);
    2927     4808665 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    2928     4808665 :       gel(z,1) = y;
    2929     4808665 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    2930             : 
    2931             :     case t_REAL:
    2932     7710229 :       return divrs(x,s);
    2933             : 
    2934             :     case t_INTMOD:
    2935        7658 :       z = cgetg(3, t_INTMOD);
    2936        7658 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    2937             : 
    2938        1043 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    2939             : 
    2940      760753 :     case t_FRAC: z = cgetg(3, t_FRAC);
    2941      760753 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    2942      760753 :       if (i == 1)
    2943             :       {
    2944      389146 :         gel(z,2) = mulsi(s, gel(x,2));
    2945      389146 :         gel(z,1) = icopy(gel(x,1));
    2946             :       }
    2947             :       else
    2948             :       {
    2949      371607 :         gel(z,2) = mulsi(s/i, gel(x,2));
    2950      371607 :         gel(z,1) = divis(gel(x,1), i);
    2951             :       }
    2952      760753 :       normalize_frac(z);
    2953      760753 :       fix_frac_if_int(z); return z;
    2954             : 
    2955     2395954 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2956     2395954 :       gel(z,1) = gdivgs(gel(x,1),s);
    2957     2395954 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    2958             : 
    2959             :     case t_PADIC: /* divpT */
    2960             :     {
    2961       87243 :       GEN p = gel(x,2);
    2962       87243 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    2963       86522 :       av = avma;
    2964       86522 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    2965             :     }
    2966             : 
    2967           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2968           0 :       gel(z,1) = ZX_copy(gel(x,1));
    2969           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    2970           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    2971             : 
    2972             :     case t_POLMOD:
    2973       22568 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    2974             : 
    2975             :     case t_RFRAC:
    2976         140 :       if (s == 1) return gcopy(x);
    2977         133 :       else if (s == -1) return gneg(x);
    2978         133 :       return div_rfrac_scal(x, stoi(s));
    2979             : 
    2980             :     case t_POL: case t_SER:
    2981       69055 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    2982       69055 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2983       69055 :       return z;
    2984             :     case t_VEC: case t_COL: case t_MAT:
    2985         301 :       z = cgetg_copy(x, &lx);
    2986         301 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    2987         301 :       return z;
    2988             : 
    2989             :   }
    2990           0 :   pari_err_TYPE2("/",x, stoi(s));
    2991             :   return NULL; /* LCOV_EXCL_LINE */
    2992             : }
    2993             : 
    2994             : /* True shift (exact multiplication by 2^n) */
    2995             : GEN
    2996    49831108 : gmul2n(GEN x, long n)
    2997             : {
    2998             :   long lx, i, k, l;
    2999             :   GEN z, a, b;
    3000             : 
    3001    49831108 :   switch(typ(x))
    3002             :   {
    3003             :     case t_INT:
    3004    35496079 :       if (n>=0) return shifti(x,n);
    3005     1569680 :       if (!signe(x)) return gen_0;
    3006     1512003 :       l = vali(x); n = -n;
    3007     1512003 :       if (n<=l) return shifti(x,-n);
    3008      248978 :       z = cgetg(3,t_FRAC);
    3009      248978 :       gel(z,1) = shifti(x,-l);
    3010      248978 :       gel(z,2) = int2n(n-l); return z;
    3011             : 
    3012             :     case t_REAL:
    3013     7071654 :       return shiftr(x,n);
    3014             : 
    3015       75363 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3016       75363 :       z = cgetg(3,t_INTMOD);
    3017       75361 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3018       75095 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3019       75103 :       gel(z,1) = icopy(b); return z;
    3020             : 
    3021      216856 :     case t_FFELT: return FF_mul2n(x,n);
    3022             : 
    3023      477805 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3024      477805 :       l = vali(a);
    3025      477805 :       k = vali(b);
    3026      477805 :       if (n+l >= k)
    3027             :       {
    3028      112868 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3029       73220 :         l = n-k; k = -k;
    3030             :       }
    3031             :       else
    3032             :       {
    3033      364937 :         k = -(l+n); l = -l;
    3034             :       }
    3035      438157 :       z = cgetg(3,t_FRAC);
    3036      438157 :       gel(z,1) = shifti(a,l);
    3037      438157 :       gel(z,2) = shifti(b,k); return z;
    3038             : 
    3039     5013850 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3040     5013850 :       gel(z,1) = gmul2n(gel(x,1),n);
    3041     5013850 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3042             : 
    3043          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3044          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3045          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3046          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3047             : 
    3048             :     case t_POLMOD:
    3049       18368 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3050             : 
    3051             :     case t_POL:
    3052      429325 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3053      429325 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3054      429325 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3055             :     case t_SER:
    3056      121985 :       if (ser_isexactzero(x)) return gcopy(x);
    3057      121971 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3058      121971 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3059      121971 :       return normalize(z); /* needed if char = 2 */
    3060             :     case t_VEC: case t_COL: case t_MAT:
    3061      892796 :       z = cgetg_copy(x, &lx);
    3062      892796 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3063      892796 :       return z;
    3064             : 
    3065             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3066          14 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3067             : 
    3068             :     case t_PADIC: /* int2n wrong if n < 0 */
    3069       16999 :       return gmul(gmul2n(gen_1,n),x);
    3070             :   }
    3071           0 :   pari_err_TYPE("gmul2n",x);
    3072             :   return NULL; /* LCOV_EXCL_LINE */
    3073             : }
    3074             : 
    3075             : /*******************************************************************/
    3076             : /*                                                                 */
    3077             : /*                              INVERSE                            */
    3078             : /*                                                                 */
    3079             : /*******************************************************************/
    3080             : static GEN
    3081       85176 : inv_polmod(GEN T, GEN x)
    3082             : {
    3083       85176 :   GEN z = cgetg(3,t_POLMOD), a;
    3084       85176 :   gel(z,1) = RgX_copy(T);
    3085       85176 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3086       62055 :     a = ginv(x);
    3087             :   else
    3088             :   {
    3089       23121 :     if (lg(T) == 5) /* quadratic fields */
    3090       17913 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3091             :     else
    3092        5208 :       a = RgXQ_inv(x, T);
    3093             :   }
    3094       85176 :   gel(z,2) = a; return z;
    3095             : }
    3096             : 
    3097             : GEN
    3098    17804169 : ginv(GEN x)
    3099             : {
    3100             :   long s;
    3101             :   pari_sp av, tetpil;
    3102             :   GEN z, y, p1, p2;
    3103             : 
    3104    17804169 :   switch(typ(x))
    3105             :   {
    3106             :     case t_INT:
    3107     1910450 :       if (is_pm1(x)) return icopy(x);
    3108     1553485 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3109     1553485 :       z = cgetg(3,t_FRAC);
    3110     1553485 :       gel(z,1) = s<0? gen_m1: gen_1;
    3111     1553485 :       gel(z,2) = absi(x); return z;
    3112             : 
    3113     3177844 :     case t_REAL: return invr(x);
    3114             : 
    3115         679 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3116         679 :       gel(z,1) = icopy(gel(x,1));
    3117         679 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3118             : 
    3119             :     case t_FRAC: {
    3120      193784 :       GEN a = gel(x,1), b = gel(x,2);
    3121      193784 :       s = signe(a);
    3122      193784 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3123       56268 :       z = cgetg(3,t_FRAC);
    3124       56268 :       gel(z,1) = icopy(b);
    3125       56268 :       gel(z,2) = icopy(a);
    3126       56268 :       normalize_frac(z); return z;
    3127             :     }
    3128             :     case t_COMPLEX:
    3129     2245219 :       av=avma;
    3130     2245219 :       p1=cxnorm(x);
    3131     2245219 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3132     2245219 :       tetpil=avma;
    3133     2245219 :       return gerepile(av,tetpil,divcR(p2,p1));
    3134             : 
    3135             :     case t_QUAD:
    3136         203 :       av=avma; p1=gnorm(x); p2=conj_i(x); tetpil=avma;
    3137         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3138             : 
    3139        4417 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3140        4417 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3141        4410 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3142        4410 :       gel(z,2) = icopy(gel(x,2));
    3143        4410 :       gel(z,3) = icopy(gel(x,3));
    3144        4410 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3145             : 
    3146       85176 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3147        7668 :     case t_FFELT: return FF_inv(x);
    3148     5371932 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3149       15120 :     case t_SER: return gdiv(gen_1,x);
    3150             : 
    3151             :     case t_RFRAC:
    3152             :     {
    3153        2912 :       GEN n = gel(x,1), d = gel(x,2);
    3154        2912 :       pari_sp av = avma, ltop;
    3155        2912 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3156             : 
    3157        2912 :       n = simplify_shallow(n);
    3158        2912 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3159             :       {
    3160        2912 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3161         679 :         ltop = avma;
    3162         679 :         z = RgX_Rg_div(d,n);
    3163             :       } else {
    3164           0 :         ltop = avma;
    3165           0 :         z = cgetg(3,t_RFRAC);
    3166           0 :         gel(z,1) = RgX_copy(d);
    3167           0 :         gel(z,2) = RgX_copy(n);
    3168             :       }
    3169         679 :       stackdummy(av, ltop);
    3170         679 :       return z;
    3171             :     }
    3172             : 
    3173             :     case t_QFR:
    3174          14 :       av = avma; z = cgetg(5, t_QFR);
    3175          14 :       gel(z,1) = gel(x,1);
    3176          14 :       gel(z,2) = negi( gel(x,2) );
    3177          14 :       gel(z,3) = gel(x,3);
    3178          14 :       gel(z,4) = negr( gel(x,4) );
    3179          14 :       return gerepileupto(av, redreal(z));
    3180             : 
    3181             :     case t_QFI:
    3182     4784971 :       y = gcopy(x);
    3183     4784971 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3184     4670039 :         togglesign(gel(y,2));
    3185     4784971 :       return y;
    3186             :     case t_MAT:
    3187        3780 :       y = RgM_inv(x);
    3188        3773 :       if (!y) pari_err_INV("ginv",x);
    3189        3703 :       return y;
    3190             :     case t_VECSMALL:
    3191             :     {
    3192           0 :       long i, lx = lg(x)-1;
    3193           0 :       y = zero_zv(lx);
    3194           0 :       for (i=1; i<=lx; i++)
    3195             :       {
    3196           0 :         long xi = x[i];
    3197           0 :         if (xi<1 || xi>lx || y[xi])
    3198           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3199           0 :         y[xi] = i;
    3200             :       }
    3201           0 :       return y;
    3202             :     }
    3203             :   }
    3204           0 :   pari_err_TYPE("inverse",x);
    3205             :   return NULL; /* LCOV_EXCL_LINE */
    3206             : }

Generated by: LCOV version 1.13