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.18.0 lcov report (development 29792-e7296d7d08) Lines: 1898 2022 93.9 %
Date: 2024-12-05 09:09:40 Functions: 102 102 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                      GENERIC OPERATIONS                        **/
      18             : /**                         (first part)                           **/
      19             : /**                                                                **/
      20             : /********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_mod
      25             : 
      26             : /* assume z[1] was created last */
      27             : #define fix_frac_if_int(z) if (equali1(gel(z,2)))\
      28             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));
      29             : 
      30             : /* assume z[1] was created last */
      31             : #define fix_frac_if_int_GC(z,tetpil) { if (equali1(gel(z,2)))\
      32             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));\
      33             : else\
      34             :   gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
      35             : 
      36             : static void
      37         105 : warn_coercion(GEN x, GEN y, GEN z)
      38             : {
      39         105 :   if (DEBUGLEVEL)
      40          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      41         105 : }
      42             : 
      43             : static long
      44          28 : kro_quad(GEN x, GEN y)
      45          28 : { pari_sp av=avma; return gc_long(av, kronecker(quad_disc(x), y)); }
      46             : 
      47             : /* is -1 not a square in Zp, assume p prime */
      48             : INLINE int
      49          42 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      50             : 
      51             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      52             : static GEN mulpp(GEN x, GEN y);
      53             : static GEN divpp(GEN x, GEN y);
      54             : /* Argument codes for inline routines
      55             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      56             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      57             :  * T: some type (to be converted to PADIC)
      58             :  */
      59             : static GEN
      60   303890888 : addRc(GEN x, GEN y) {
      61   303890888 :   GEN z = cgetg(3,t_COMPLEX);
      62   303873578 :   gel(z,1) = gadd(x,gel(y,1));
      63   303855214 :   gel(z,2) = gcopy(gel(y,2)); return z;
      64             : }
      65             : static GEN
      66   351214560 : mulRc(GEN x, GEN y) {
      67   351214560 :   GEN z = cgetg(3,t_COMPLEX);
      68   351190394 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      69   351633161 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      70             : }
      71             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      72             : static GEN
      73          49 : mulRc_direct(GEN x, GEN y) {
      74          49 :   GEN z = cgetg(3,t_COMPLEX);
      75          49 :   gel(z,1) = gmul(x,gel(y,1));
      76          49 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      77             : }
      78             : static GEN
      79     2352670 : divRc(GEN x, GEN y) {
      80     2352670 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      81     2352662 :   GEN z = cgetg(3,t_COMPLEX);
      82     2352665 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      83     2352654 :   gel(z,2) = gmul(mt, gel(y,2));
      84     2352655 :   return z;
      85             : }
      86             : static GEN
      87    23109710 : divcR(GEN x, GEN y) {
      88    23109710 :   GEN z = cgetg(3,t_COMPLEX);
      89    23109671 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      90    23109146 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      91             : }
      92             : static GEN
      93        1274 : addRq(GEN x, GEN y) {
      94        1274 :   GEN z = cgetg(4,t_QUAD);
      95        1274 :   gel(z,1) = ZX_copy(gel(y,1));
      96        1274 :   gel(z,2) = gadd(x, gel(y,2));
      97        1274 :   gel(z,3) = gcopy(gel(y,3)); return z;
      98             : }
      99             : static GEN
     100        2352 : mulRq(GEN x, GEN y) {
     101        2352 :   GEN z = cgetg(4,t_QUAD);
     102        2352 :   gel(z,1) = ZX_copy(gel(y,1));
     103        2352 :   gel(z,2) = gmul(x,gel(y,2));
     104        2352 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     105             : }
     106             : static GEN
     107          77 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     108          77 :   long i = gexpo(x) - gexpo(y);
     109          77 :   if (i > 0) prec += nbits2extraprec(i);
     110          77 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     111             : }
     112             : static GEN
     113    37510904 : mulrfrac(GEN x, GEN y)
     114             : {
     115             :   pari_sp av;
     116    37510904 :   GEN z, a = gel(y,1), b = gel(y,2);
     117    37510904 :   if (is_pm1(a)) /* frequent special case */
     118             :   {
     119    12654789 :     z = divri(x, b); if (signe(a) < 0) togglesign(z);
     120    12653993 :     return z;
     121             :   }
     122    24855981 :   av = avma;
     123    24855981 :   return gerepileuptoleaf(av, divri(mulri(x,a), b));
     124             : }
     125             : static GEN
     126          28 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     127          28 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     128             : }
     129             : static GEN
     130          63 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     131          63 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     132             : }
     133             : static GEN
     134          42 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     135          42 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     136             : }
     137             : /* y PADIC, x + y by converting x to padic */
     138             : static GEN
     139           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     140             : 
     141           7 :   if (!valp(y)) z = cvtop2(x,y);
     142             :   else {
     143           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     144           7 :     z  = cvtop(x, gel(y,2), l);
     145             :   }
     146           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     147             : }
     148             : /* y PADIC, x * y by converting x to padic */
     149             : static GEN
     150     4982894 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     151     4982894 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     152             : }
     153             : /* y PADIC, non zero x / y by converting x to padic */
     154             : static GEN
     155        3660 : divTp(GEN x, GEN y) { pari_sp av = avma;
     156        3660 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     157             : }
     158             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     159             :  * converted to O(p^e) and division by 0 */
     160             : static GEN
     161     1204910 : divpT(GEN x, GEN y) { pari_sp av = avma;
     162     1204910 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     163             : }
     164             : 
     165             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     166             :  * clean memory from z on */
     167             : static GEN
     168     3690733 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     169     3690733 :   if (lgefint(X) == 3) {
     170     3367566 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     171     3367566 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     172             :   }
     173             :   else {
     174      323167 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     175      323166 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     176             :   }
     177     3690732 :   gel(z,1) = icopy(X); return z;
     178             : }
     179             : static GEN
     180     1158553 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     181     1158553 :   if (lgefint(X) == 3) {
     182      784519 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     183      784519 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     184             :   }
     185             :   else {
     186      374034 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     187      374034 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     188             :   }
     189     1158553 :   gel(z,1) = icopy(X); return z;
     190             : }
     191             : /* cf add_intmod_same */
     192             : static GEN
     193     3005203 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     194     3005203 :   if (lgefint(X) == 3) {
     195     2863550 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     196     2863550 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     197             :   }
     198             :   else
     199      141653 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     200     3005206 :   gel(z,1) = icopy(X); return z;
     201             : }
     202             : /* cf add_intmod_same */
     203             : static GEN
     204      341964 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     205             : {
     206      341964 :   if (lgefint(X) == 3) {
     207      319422 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     208      319359 :     set_avma((pari_sp)z); gel(z,2) = utoi(u);
     209             :   }
     210             :   else
     211       22542 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     212      341901 :   gel(z,1) = icopy(X); return z;
     213             : }
     214             : 
     215             : /*******************************************************************/
     216             : /*                                                                 */
     217             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     218             : /*                                                                 */
     219             : /* (static routines are not memory clean, but OK for gerepileupto) */
     220             : /*******************************************************************/
     221             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     222             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     223             : static GEN
     224     9857237 : rfrac_denom_mul_scal(GEN d, GEN y)
     225             : {
     226     9857237 :   GEN D = RgX_Rg_mul(d, y);
     227     9857237 :   if (lg(D) != lg(d))
     228             :   { /* try to generate a meaningful diagnostic */
     229           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     230           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     231             :   }
     232     9857237 :   return D;
     233             : }
     234             : 
     235             : static int
     236    57829141 : Leading_is_neg(GEN x)
     237             : {
     238   122278619 :   while (typ(x) == t_POL) x = leading_coeff(x);
     239    57829141 :   return is_real_t(typ(x))? gsigne(x) < 0: 0;
     240             : }
     241             : 
     242             : static int
     243   154482632 : transtype(GEN x) { return x != gen_1 && typ(x) != t_PADIC; }
     244             : 
     245             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     246             : GEN
     247    57833391 : gred_rfrac_simple(GEN n, GEN d)
     248             : {
     249             :   GEN _1n, _1d, c, cn, cd, z;
     250    57833391 :   long dd = degpol(d);
     251             : 
     252    57833391 :   if (dd <= 0)
     253             :   {
     254        4250 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     255        4250 :     n = gdiv(n, gel(d,2));
     256        4250 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     257        4250 :     return n;
     258             :   }
     259    57829141 :   if (Leading_is_neg(d)) { d = gneg(d); n = gneg(n); }
     260    57829141 :   _1n = Rg_get_1(n);
     261    57829141 :   _1d = Rg_get_1(d);
     262    57829141 :   if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
     263    57829141 :   if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
     264    57829134 :   cd = content(d);
     265    59699609 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     266    57829134 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     267    57829134 :   if (!gequal1(cd)) {
     268     6568093 :     d = RgX_Rg_div(d,cd);
     269     6568093 :     if (!gequal1(cn))
     270             :     {
     271     1294238 :       if (gequal0(cn)) {
     272          49 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     273           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     274           0 :         c = gen_1;
     275             :       } else {
     276     1294189 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     277     1294189 :         c = gdiv(cn,cd);
     278             :       }
     279             :     }
     280             :     else
     281     5273855 :       c = ginv(cd);
     282             :   } else {
     283    51261041 :     if (!gequal1(cn))
     284             :     {
     285     3277487 :       if (gequal0(cn)) {
     286        1484 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     287          21 :         c = gen_1;
     288             :       } else {
     289     3276003 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     290     3276003 :         c = cn;
     291             :       }
     292             :     } else {
     293    47983554 :       GEN y = cgetg(3,t_RFRAC);
     294    47983554 :       gel(y,1) = gcopy(n);
     295    47983554 :       gel(y,2) = RgX_copy(d); return y;
     296             :     }
     297             :   }
     298             : 
     299     9844068 :   if (typ(c) == t_POL)
     300             :   {
     301      900639 :     z = c;
     302      939069 :     do { z = content(z); } while (typ(z) == t_POL);
     303      900639 :     cd = denom_i(z);
     304      900639 :     cn = gmul(c, cd);
     305             :   }
     306             :   else
     307             :   {
     308     8943429 :     cn = numer_i(c);
     309     8943429 :     cd = denom_i(c);
     310             :   }
     311     9844068 :   z = cgetg(3,t_RFRAC);
     312     9844068 :   gel(z,1) = gmul(n, cn);
     313     9844068 :   gel(z,2) = d = rfrac_denom_mul_scal(d, cd);
     314             :   /* can happen: Pol(O(17^-1)) / Pol([Mod(9,23), O(23^-3)]) */
     315     9844068 :   if (!signe(d)) pari_err_INV("gred_rfrac_simple", d);
     316     9844068 :   return z;
     317             : }
     318             : 
     319             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     320             : static GEN
     321      248374 : fix_rfrac(GEN x, long d)
     322             : {
     323             :   GEN z, N, D;
     324      248374 :   if (!d || typ(x) == t_POL) return x;
     325      171710 :   z = cgetg(3, t_RFRAC);
     326      171710 :   N = gel(x,1);
     327      171710 :   D = gel(x,2);
     328      171710 :   if (d > 0) {
     329        8722 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     330      165935 :                                                    : monomialcopy(N,d,varn(D));
     331      165872 :     gel(z, 2) = RgX_copy(D);
     332             :   } else {
     333        5838 :     gel(z, 1) = gcopy(N);
     334        5838 :     gel(z, 2) = RgX_shift(D, -d);
     335             :   }
     336      171710 :   return z;
     337             : }
     338             : 
     339             : /* assume d != 0 */
     340             : static GEN
     341    44766283 : gred_rfrac2(GEN n, GEN d)
     342             : {
     343             :   GEN y, z, _1n, _1d;
     344             :   long v, vd, vn;
     345             : 
     346    44766283 :   n = simplify_shallow(n);
     347    44766283 :   if (isintzero(n)) return scalarpol(Rg_get_0(d), varn(d));
     348    37937776 :   d = simplify_shallow(d);
     349    37937776 :   if (typ(d) != t_POL) return gdiv(n,d);
     350    36758990 :   vd = varn(d);
     351    36758990 :   if (typ(n) != t_POL)
     352             :   {
     353    20623729 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     354    20622322 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     355           0 :     pari_err_BUG("gred_rfrac2 [incompatible variables]");
     356             :   }
     357    16135261 :   vn = varn(n);
     358    16135261 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     359    15997552 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     360    15837779 :   _1n = Rg_get_1(n);
     361    15837779 :   _1d = Rg_get_1(d);
     362    15837779 :   if (transtype(_1n) && !gidentical(_1n, _1d)) d = gmul(d, _1n);
     363    15837779 :   if (transtype(_1d) && !gidentical(_1n, _1d)) n = gmul(n, _1d);
     364             : 
     365             :   /* now n and d are t_POLs in the same variable */
     366    15837779 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     367    15837779 :   if (!degpol(d))
     368             :   {
     369    12757958 :     n = RgX_Rg_div(n,gel(d,2));
     370    12757958 :     return v? RgX_mulXn(n,v): n;
     371             :   }
     372             : 
     373             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     374     3079821 :   if (!isinexact(n) && !isinexact(d))
     375             :   {
     376     3079576 :     y = RgX_divrem(n, d, &z);
     377     3079576 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     378      248129 :     z = RgX_gcd(d, z);
     379      248129 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     380             :   }
     381      248374 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     382             : }
     383             : 
     384             : /* x,y t_INT, return x/y in reduced form */
     385             : GEN
     386   134835266 : Qdivii(GEN x, GEN y)
     387             : {
     388   134835266 :   pari_sp av = avma;
     389             :   GEN r, q;
     390             : 
     391   134835266 :   if (lgefint(y) == 3)
     392             :   {
     393   116989720 :     q = Qdiviu(x, y[2]);
     394   116987955 :     if (signe(y) > 0) return q;
     395    11012005 :     if (typ(q) == t_INT) togglesign(q); else togglesign_safe(&gel(q,1));
     396    11012283 :     return q;
     397             :   }
     398    17845546 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     399    17846959 :   if (equali1(x))
     400             :   {
     401     5178187 :     if (!signe(y)) pari_err_INV("gdiv",y);
     402     5178082 :     retmkfrac(signe(y) < 0? gen_m1: gen_1, absi(y));
     403             :   }
     404    12668804 :   q = dvmdii(x,y,&r);
     405    12668662 :   if (r == gen_0) return q; /* gen_0 intended */
     406     8316077 :   r = gcdii(y, r);
     407     8316118 :   if (lgefint(r) == 3)
     408             :   {
     409     7443086 :     ulong t = r[2];
     410     7443086 :     set_avma(av);
     411     7443049 :     if (t == 1) q = mkfraccopy(x,y);
     412             :     else
     413             :     {
     414     2648142 :       q = cgetg(3,t_FRAC);
     415     2648322 :       gel(q,1) = diviuexact(x,t);
     416     2648226 :       gel(q,2) = diviuexact(y,t);
     417             :     }
     418             :   }
     419             :   else
     420             :   { /* rare: r and q left on stack for efficiency */
     421      873032 :     q = cgetg(3,t_FRAC);
     422      873054 :     gel(q,1) = diviiexact(x,r);
     423      873055 :     gel(q,2) = diviiexact(y,r);
     424             :   }
     425     8316046 :   normalize_frac(q); return q;
     426             : }
     427             : 
     428             : /* x t_INT, return x/y in reduced form */
     429             : GEN
     430   136305028 : Qdiviu(GEN x, ulong y)
     431             : {
     432   136305028 :   pari_sp av = avma;
     433             :   ulong r, t;
     434             :   GEN q;
     435             : 
     436   136305028 :   if (y == 1) return icopy(x);
     437   115518561 :   if (!y) pari_err_INV("gdiv",gen_0);
     438   115518561 :   if (equali1(x)) retmkfrac(gen_1, utoipos(y));
     439   108355210 :   q = absdiviu_rem(x,y,&r);
     440   108350476 :   if (!r)
     441             :   {
     442    62515815 :     if (signe(x) < 0) togglesign(q);
     443    62515813 :     return q;
     444             :   }
     445    45834661 :   t = ugcd(y, r); set_avma(av);
     446    45841625 :   if (t == 1) retmkfrac(icopy(x), utoipos(y));
     447    18372029 :   retmkfrac(diviuexact(x,t), utoipos(y / t));
     448             : }
     449             : 
     450             : /* x t_INT, return x/y in reduced form */
     451             : GEN
     452     1560428 : Qdivis(GEN x, long y)
     453             : {
     454     1560428 :   pari_sp av = avma;
     455             :   ulong r, t;
     456             :   long s;
     457             :   GEN q;
     458             : 
     459     1560428 :   if (y > 0) return Qdiviu(x, y);
     460      162429 :   if (!y) pari_err_INV("gdiv",gen_0);
     461      162429 :   s = signe(x);
     462      162429 :   if (!s) return gen_0;
     463      132202 :   if (y < 0) { y = -y; s = -s; }
     464      132202 :   if (y == 1) { x = icopy(x); setsigne(x,s); return x; }
     465      131922 :   if (equali1(x)) retmkfrac(s > 0? gen_1: gen_m1, utoipos(y));
     466      129759 :   q = absdiviu_rem(x,y,&r);
     467      129759 :   if (!r)
     468             :   {
     469       55356 :     if (s < 0) togglesign(q);
     470       55356 :     return q;
     471             :   }
     472       74403 :   t = ugcd(y, r); set_avma(av); q = cgetg(3, t_FRAC);
     473       74403 :   if (t != 1) { x = diviuexact(x,t); y /= t; } else x = icopy(x);
     474       74403 :   gel(q,1) = x; setsigne(x, s);
     475       74403 :   gel(q,2) = utoipos(y); return q;
     476             : }
     477             : 
     478             : /*******************************************************************/
     479             : /*                                                                 */
     480             : /*                          CONJUGATION                            */
     481             : /*                                                                 */
     482             : /*******************************************************************/
     483             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     484             : static GEN
     485       18326 : quad_polmod_conj(GEN x, GEN y)
     486             : {
     487             :   GEN z, u, v, a, b;
     488       18326 :   if (typ(x) != t_POL) return gcopy(x);
     489       18326 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     490       18326 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     491       18326 :   b = gel(y,3); v = gel(x,2);
     492       18326 :   z = cgetg(4, t_POL); z[1] = x[1];
     493       18326 :   gel(z,2) = gsub(v, gdiv(gmul(u,b), a));
     494       18326 :   gel(z,3) = gneg(u); return z;
     495             : }
     496             : static GEN
     497       18326 : quad_polmod_norm(GEN x, GEN y)
     498             : {
     499             :   GEN z, u, v, a, b, c;
     500       18326 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0) return gsqr(x);
     501       18326 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     502       18326 :   b = gel(y,3); v = gel(x,2);
     503       18326 :   c = gel(y,2);
     504       18326 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     505       18326 :   if (!gequal1(a)) z = gdiv(z, a);
     506       18326 :   return gadd(z, gsqr(v));
     507             : }
     508             : 
     509             : GEN
     510    19572822 : conj_i(GEN x)
     511             : {
     512    19572822 :   switch(typ(x))
     513             :   {
     514     6536622 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
     515     6536622 :       return x;
     516             : 
     517    12872584 :     case t_COMPLEX: return mkcomplex(gel(x,1), gneg(gel(x,2)));
     518             : 
     519         924 :     case t_QUAD:
     520             :     {
     521         924 :       GEN y = cgetg(4,t_QUAD);
     522         924 :       gel(y,1) = gel(x,1);
     523         924 :       gel(y,2) = gequal0(gmael(x,1,3))? gel(x,2)
     524         924 :                                       : gadd(gel(x,2), gel(x,3));
     525         924 :       gel(y,3) = gneg(gel(x,3)); return y;
     526             :     }
     527         350 :     case t_POL: pari_APPLY_pol_normalized(conj_i(gel(x,i)));
     528       31591 :     case t_SER: pari_APPLY_ser_normalized(conj_i(gel(x,i)));
     529             : 
     530      152780 :     case t_RFRAC:
     531             :     case t_VEC:
     532             :     case t_COL:
     533      550461 :     case t_MAT: pari_APPLY_same(conj_i(gel(x,i)));
     534             : 
     535           0 :     case t_POLMOD:
     536             :     {
     537           0 :       GEN X = gel(x,1);
     538           0 :       long d = degpol(X);
     539           0 :       if (d < 2) return x;
     540           0 :       if (d == 2) return mkpolmod(quad_polmod_conj(gel(x,2), X), X);
     541             :     }
     542             :   }
     543           0 :   pari_err_TYPE("gconj",x);
     544             :   return NULL; /* LCOV_EXCL_LINE */
     545             : }
     546             : GEN
     547      133181 : gconj(GEN x)
     548      133181 : { pari_sp av = avma; return gerepilecopy(av, conj_i(x)); }
     549             : 
     550             : GEN
     551          84 : conjvec(GEN x,long prec)
     552             : {
     553             :   long lx, s, i;
     554             :   GEN z;
     555             : 
     556          84 :   switch(typ(x))
     557             :   {
     558           0 :     case t_INT: case t_INTMOD: case t_FRAC:
     559           0 :       return mkcolcopy(x);
     560             : 
     561           0 :     case t_COMPLEX: case t_QUAD:
     562           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     563             : 
     564          28 :     case t_FFELT:
     565          28 :       return FF_conjvec(x);
     566             : 
     567           0 :     case t_VEC: case t_COL:
     568           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     569           0 :       if (lx == 1) return z;
     570           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     571           0 :       s = lgcols(z);
     572           0 :       for (i=2; i<lx; i++)
     573             :       {
     574           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     575           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     576             :       }
     577           0 :       break;
     578             : 
     579          56 :     case t_POLMOD: {
     580          56 :       GEN T = gel(x,1), r;
     581             :       pari_sp av;
     582             : 
     583          56 :       lx = lg(T);
     584          56 :       if (lx <= 3) return cgetg(1,t_COL);
     585          56 :       x = gel(x,2);
     586         238 :       for (i=2; i<lx; i++)
     587             :       {
     588         189 :         GEN c = gel(T,i);
     589         189 :         switch(typ(c)) {
     590           7 :           case t_INTMOD: {
     591           7 :             GEN p = gel(c,1);
     592             :             pari_sp av;
     593           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     594           7 :             av = avma;
     595           7 :             T = RgX_to_FpX(T,p);
     596           7 :             x = RgX_to_FpX(x, p);
     597           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     598           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     599           7 :             return gerepileupto(av, z);
     600             :           }
     601         182 :           case t_INT:
     602         182 :           case t_FRAC: break;
     603           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     604             :         }
     605             :       }
     606          49 :       if (typ(x) != t_POL)
     607             :       {
     608           0 :         if (!is_rational_t(typ(x)))
     609           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     610           0 :         retconst_col(lx-3, gcopy(x));
     611             :       }
     612          49 :       RgX_check_QX(x,"conjvec");
     613          49 :       av = avma;
     614          49 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     615          49 :       r = cleanroots(T,prec);
     616          49 :       z = cgetg(lx-2,t_COL);
     617         182 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     618          49 :       return gerepileupto(av, z);
     619             :     }
     620             : 
     621           0 :     default:
     622           0 :       pari_err_TYPE("conjvec",x);
     623             :       return NULL; /* LCOV_EXCL_LINE */
     624             :   }
     625           0 :   return z;
     626             : }
     627             : 
     628             : /********************************************************************/
     629             : /**                                                                **/
     630             : /**                           ADDITION                             **/
     631             : /**                                                                **/
     632             : /********************************************************************/
     633             : /* x, y compatible PADIC, op = add or sub */
     634             : static GEN
     635    19968884 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     636             : {
     637    19968884 :   pari_sp av = avma;
     638             :   long d,e,r,rx,ry;
     639    19968884 :   GEN u, z, mod, p = gel(x,2);
     640             :   int swap;
     641             : 
     642    19968884 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     643    19968564 :   e = valp(x);
     644    19968564 :   r = valp(y); d = r-e;
     645    19968564 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     646    19968564 :   rx = precp(x);
     647    19968564 :   ry = precp(y);
     648    19968564 :   if (d) /* v(x) < v(y) */
     649             :   {
     650    10703647 :     r = d+ry; z = powiu(p,d);
     651    10703733 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     652    10703752 :     z = mulii(z,gel(y,4));
     653    10703706 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     654             :   }
     655             :   else
     656             :   {
     657             :     long c;
     658     9264917 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     659     9264917 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     660     9265594 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     661             :     {
     662      138237 :       set_avma(av); return zeropadic(p, e+r);
     663             :     }
     664     9126756 :     if (c)
     665             :     {
     666     3351429 :       mod = diviiexact(mod, powiu(p,c));
     667     3351432 :       r -= c;
     668     3351432 :       e += c;
     669             :     }
     670             :   }
     671    19830414 :   u = modii(u, mod);
     672    19829531 :   set_avma(av); z = cgetg(5,t_PADIC);
     673    19829040 :   z[1] = evalprecp(r) | evalvalp(e);
     674    19828930 :   gel(z,2) = icopy(p);
     675    19828683 :   gel(z,3) = icopy(mod);
     676    19828361 :   gel(z,4) = icopy(u); return z;
     677             : }
     678             : /* Rg_to_Fp(t_FRAC) without GC */
     679             : static GEN
     680      226838 : Q_to_Fp(GEN x, GEN p)
     681      226838 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     682             : /* return x + y, where y t_PADIC and x is a nonzero t_INT or t_FRAC */
     683             : static GEN
     684     4762872 : addQp(GEN x, GEN y)
     685             : {
     686     4762872 :   pari_sp av = avma;
     687     4762872 :   long d, r, e, vy = valp(y), py = precp(y);
     688     4762872 :   GEN z, q, mod, u, p = gel(y,2);
     689             : 
     690     4762872 :   e = Q_pvalrem(x, p, &x);
     691     4762843 :   d = vy - e; r = d + py;
     692     4762843 :   if (r <= 0) { set_avma(av); return gcopy(y); }
     693     4760883 :   mod = gel(y,3);
     694     4760883 :   u   = gel(y,4);
     695     4760883 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     696             : 
     697     4760881 :   if (d > 0)
     698             :   {
     699     1376442 :     q = powiu(p,d);
     700     1376538 :     mod = mulii(mod, q);
     701     1376531 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     702     1376531 :     u = addii(x,  mulii(u, q));
     703             :   }
     704     3384439 :   else if (d < 0)
     705             :   {
     706      405332 :     q = powiu(p,-d);
     707      405332 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     708      405332 :     u = addii(u, mulii(x, q));
     709      405332 :     r = py; e = vy;
     710             :   }
     711             :   else
     712             :   {
     713             :     long c;
     714     2979107 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     715     2979107 :     u = addii(u, x);
     716     2979108 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     717             :     {
     718        1205 :       set_avma(av); return zeropadic(p,e+r);
     719             :     }
     720     2977887 :     if (c)
     721             :     {
     722      970645 :       mod = diviiexact(mod, powiu(p,c));
     723      970645 :       r -= c;
     724      970645 :       e += c;
     725             :     }
     726             :   }
     727     4759753 :   u = modii(u, mod); set_avma(av);
     728     4759721 :   z = cgetg(5,t_PADIC);
     729     4759688 :   z[1] = evalprecp(r) | evalvalp(e);
     730     4759675 :   gel(z,2) = icopy(p);
     731     4759665 :   gel(z,3) = icopy(mod);
     732     4759660 :   gel(z,4) = icopy(u); return z;
     733             : }
     734             : 
     735             : /* Mod(x,X) + Mod(y,X) */
     736             : #define addsub_polmod_same addsub_polmod_scal
     737             : /* Mod(x,X) +/- Mod(y,Y) */
     738             : static GEN
     739        7203 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     740             : {
     741        7203 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     742        7203 :   GEN z = cgetg(3,t_POLMOD);
     743        7203 :   long vx = varn(X), vy = varn(Y);
     744        7203 :   if (vx==vy) {
     745             :     pari_sp av;
     746          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     747          14 :     warn_coercion(X,Y,gel(z,1));
     748          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     749             :   }
     750        7189 :   if (varncmp(vx, vy) < 0)
     751        7189 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     752             :   else
     753           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     754        7189 :   gel(z,2) = op(x, y); return z;
     755             : }
     756             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     757             : static GEN
     758    13425153 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     759    13425153 : { retmkpolmod(degpol(Y)? op(y, x): gen_0, RgX_copy(Y)); }
     760             : 
     761             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     762             : static GEN
     763      408094 : add_ser_scal(GEN y, GEN x)
     764             : {
     765             :   long i, v, ly, vy;
     766             :   GEN z;
     767             : 
     768      408094 :   if (isrationalzero(x)) return gcopy(y);
     769      381837 :   ly = lg(y);
     770      381837 :   v = valser(y);
     771      381837 :   if (v < 3-ly) return gcopy(y);
     772             :   /* v + ly >= 3 */
     773      381585 :   if (v < 0)
     774             :   {
     775        1162 :     z = cgetg(ly,t_SER); z[1] = y[1];
     776        3276 :     for (i = 2; i <= 1-v; i++) gel(z,i) = gcopy(gel(y,i));
     777        1162 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     778        3157 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     779        1162 :     return normalizeser(z);
     780             :   }
     781      380423 :   vy = varn(y);
     782      380423 :   if (v > 0)
     783             :   {
     784       19299 :     if (ser_isexactzero(y))
     785        9338 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, v);
     786        9961 :     y -= v; ly += v;
     787        9961 :     z = cgetg(ly,t_SER);
     788        9961 :     x = gcopy(x);
     789       20559 :     for (i=3; i<=v+1; i++) gel(z,i) = gen_0;
     790             :   }
     791      361124 :   else if (ser_isexactzero(y)) return gcopy(y);
     792             :   else
     793             :   { /* v = 0, ly >= 3 */
     794      361117 :     z = cgetg(ly,t_SER);
     795      361117 :     x = gadd(x, gel(y,2));
     796      361117 :     i = 3;
     797             :   }
     798     1598979 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     799      371078 :   gel(z,2) = x;
     800      371078 :   z[1] = evalsigne(1) | _evalvalser(0) | evalvarn(vy);
     801      371078 :   return gequal0(x)? normalizeser(z): z;
     802             : }
     803             : static long
     804     7161978 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     805             : /* x,y t_SER in the same variable: x+y */
     806             : static GEN
     807     3581388 : ser_add(GEN x, GEN y)
     808             : {
     809     3581388 :   long i, lx,ly, n = valser(y) - valser(x);
     810             :   GEN z;
     811     3581388 :   if (n < 0) { n = -n; swap(x,y); }
     812             :   /* valser(x) <= valser(y) */
     813     3581388 :   lx = _serprec(x);
     814     3581388 :   if (lx == 2) /* don't lose type information */
     815             :   {
     816         798 :     z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
     817         798 :     setvalser(z, valser(x)); return z;
     818             :   }
     819     3580590 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     820     3580590 :   if (n)
     821             :   {
     822      107782 :     if (n+2 > lx) return gcopy(x);
     823      107082 :     z = cgetg(ly,t_SER);
     824      819116 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     825      506624 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     826             :   } else {
     827     3472808 :     z = cgetg(ly,t_SER);
     828    17589354 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     829             :   }
     830     3579890 :   z[1] = x[1]; return normalizeser(z);
     831             : }
     832             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     833             : static GEN
     834     8830966 : add_rfrac_scal(GEN y, GEN x)
     835             : {
     836             :   pari_sp av;
     837             :   GEN n;
     838             : 
     839     8830966 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     840     5157194 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     841     5157194 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     842             : }
     843             : 
     844             : /* x "scalar", ty != t_MAT and nonscalar */
     845             : static GEN
     846    26105032 : add_scal(GEN y, GEN x, long ty)
     847             : {
     848    26105032 :   switch(ty)
     849             :   {
     850    21192808 :     case t_POL: return RgX_Rg_add(y, x);
     851      408066 :     case t_SER: return add_ser_scal(y, x);
     852     4468785 :     case t_RFRAC: return add_rfrac_scal(y, x);
     853           0 :     case t_COL: return RgC_Rg_add(y, x);
     854       35364 :     case t_VEC:
     855       35364 :       if (isintzero(x)) return gcopy(y);
     856         175 :       break;
     857             :   }
     858         184 :   pari_err_TYPE2("+",x,y);
     859             :   return NULL; /* LCOV_EXCL_LINE */
     860             : }
     861             : 
     862             : /* assumes z = cget(3, t_FRAC) comes first in stack, then a, then b */
     863             : static GEN
     864    14921658 : setfrac(GEN z, GEN a, GEN b)
     865             : {
     866    14921658 :   gel(z,1) = icopy_avma(a, (pari_sp)z);
     867    14921649 :   gel(z,2) = icopy_avma(b, (pari_sp)gel(z,1));
     868    14921712 :   set_avma((pari_sp)gel(z,2)); return z;
     869             : }
     870             : /* z <- a / (b*Q), (Q,a) = 1 */
     871             : static GEN
     872    14115685 : addsub_frac_i(GEN z, GEN Q, GEN a, GEN b)
     873             : {
     874    14115685 :   GEN q = Qdivii(a, b); /* != 0 */
     875    14115744 :   if (typ(q) == t_INT)
     876             :   {
     877     1976920 :     gel(z,1) = gerepileuptoint((pari_sp)Q, q);
     878     1976920 :     gel(z,2) = Q; return z;
     879             :   }
     880    12138824 :   return setfrac(z, gel(q,1), mulii(gel(q,2), Q));
     881             : }
     882             : static GEN
     883    28640424 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     884             : {
     885    28640424 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     886    28640424 :   GEN y1 = gel(y,1), y2 = gel(y,2), z, Q, q, r, n, delta;
     887    28640424 :   int s = cmpii(x2, y2);
     888             : 
     889             :   /* common denominator: (x1 op y1) / x2 */
     890    28640428 :   if (!s)
     891             :   {
     892    10185264 :     pari_sp av = avma;
     893    10185264 :     return gerepileupto(av, Qdivii(op(x1, y1), x2));
     894             :   }
     895    18455164 :   z = cgetg(3, t_FRAC);
     896    18455166 :   if (s < 0)
     897             :   {
     898    10460795 :     Q = dvmdii(y2, x2, &r);
     899             :     /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     900    10460772 :     if (r == gen_0) return addsub_frac_i(z, Q, op(mulii(Q,x1), y1), x2);
     901     2762606 :     delta = gcdii(x2,r);
     902             :   }
     903             :   else
     904             :   {
     905     7994371 :     Q = dvmdii(x2, y2, &r);
     906             :     /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     907     7994386 :     if (r == gen_0) return addsub_frac_i(z, Q, op(x1, mulii(Q,y1)), y2);
     908     1576804 :     delta = gcdii(y2,r);
     909             :   }
     910             :   /* delta = gcd(x2,y2) */
     911     4339429 :   if (equali1(delta))
     912             :   { /* numerator is nonzero */
     913     1556526 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     914     1556526 :     gel(z,2) = mulii(x2,y2); return z;
     915             :   }
     916     2782903 :   x2 = diviiexact(x2,delta);
     917     2782901 :   y2 = diviiexact(y2,delta);
     918     2782902 :   n = op(mulii(x1,y2), mulii(y1,x2)); /* != 0 */
     919     2782901 :   q = dvmdii(n, delta, &r);
     920     2782900 :   if (r == gen_0) return setfrac(z, q, mulii(x2, y2));
     921     2528142 :   r = gcdii(delta, r);
     922     2528144 :   if (!equali1(r)) { n = diviiexact(n, r); delta = diviiexact(delta, r); }
     923     2528144 :   return setfrac(z, n, mulii(mulii(x2, y2), delta));
     924             : }
     925             : 
     926             : /* assume x2, y2 are t_POLs in the same variable */
     927             : static GEN
     928     3039747 : add_rfrac(GEN x, GEN y)
     929             : {
     930     3039747 :   pari_sp av = avma;
     931     3039747 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     932     3039747 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     933             : 
     934     3039747 :   delta = RgX_gcd(x2,y2);
     935     3039747 :   if (!degpol(delta))
     936             :   {
     937         672 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     938         672 :     d = RgX_mul(x2, y2);
     939         672 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     940             :   }
     941     3039075 :   x2 = RgX_div(x2,delta);
     942     3039075 :   y2 = RgX_div(y2,delta);
     943     3039075 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     944     3039075 :   if (!signe(n))
     945             :   {
     946      721425 :     n = simplify_shallow(n);
     947      721425 :     if (isexactzero(n))
     948             :     {
     949      721418 :       if (isrationalzero(n)) { set_avma(av); return zeropol(varn(x2)); }
     950          14 :       return gerepileupto(av, scalarpol(n, varn(x2)));
     951             :     }
     952           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     953             :   }
     954     2317650 :   if (degpol(n) == 0)
     955     1150603 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     956     1167047 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     957     1167047 :   if (isexactzero(r))
     958             :   {
     959             :     GEN z;
     960      228087 :     d = RgX_mul(x2, y2);
     961             :     /* "constant" denominator ? */
     962      228087 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     963      228087 :     return gerepileupto(av, z);
     964             :   }
     965      938960 :   r = RgX_gcd(delta, r);
     966      938960 :   if (degpol(r))
     967             :   {
     968      160740 :     n = RgX_div(n, r);
     969      160740 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     970             :   }
     971             :   else
     972      778220 :     d = RgX_mul(gel(x,2), y2);
     973      938960 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     974             : }
     975             : 
     976             : GEN
     977  5884005335 : gadd(GEN x, GEN y)
     978             : {
     979  5884005335 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     980             :   pari_sp av, tetpil;
     981             :   GEN z, p1;
     982             : 
     983  5884005335 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     984             :   {
     985  2910823644 :     case t_INT: return addii(x,y);
     986  1637221367 :     case t_REAL: return addrr(x,y);
     987     1521962 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     988     1521962 :       z = cgetg(3,t_INTMOD);
     989     1521962 :       if (X==Y || equalii(X,Y))
     990     1521948 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     991          14 :       gel(z,1) = gcdii(X,Y);
     992          14 :       warn_coercion(X,Y,gel(z,1));
     993          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     994          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     995             :     }
     996    23196376 :     case t_FRAC: return addsub_frac(x,y,addii);
     997   354654368 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     998   354164921 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     999   354822320 :       if (isintzero(gel(z,2)))
    1000             :       {
    1001      521178 :         set_avma((pari_sp)(z+3));
    1002      521178 :         return gadd(gel(x,1),gel(y,1));
    1003             :       }
    1004   354181042 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
    1005   354202709 :       return z;
    1006    14531380 :     case t_PADIC:
    1007    14531380 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1008    14531193 :       return addsub_pp(x,y, addii);
    1009         672 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1010         672 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1011         672 :       gel(z,1) = ZX_copy(gel(x,1));
    1012         672 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
    1013         672 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
    1014    10161232 :     case t_POLMOD:
    1015    10161232 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1016    10154064 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
    1017        7168 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
    1018    14006793 :     case t_FFELT: return FF_add(x,y);
    1019    72966668 :     case t_POL:
    1020    72966668 :       vx = varn(x);
    1021    72966668 :       vy = varn(y);
    1022    72966668 :       if (vx != vy) {
    1023      855113 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
    1024       66302 :         else                     return RgX_Rg_add(y, x);
    1025             :       }
    1026    72111555 :       return RgX_add(x, y);
    1027     3576033 :     case t_SER:
    1028     3576033 :       vx = varn(x);
    1029     3576033 :       vy = varn(y);
    1030     3576033 :       if (vx != vy) {
    1031          28 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1032          21 :         else                     return add_ser_scal(y, x);
    1033             :       }
    1034     3576005 :       return ser_add(x, y);
    1035     4331648 :     case t_RFRAC:
    1036     4331648 :       vx = varn(gel(x,2));
    1037     4331648 :       vy = varn(gel(y,2));
    1038     4331648 :       if (vx != vy) {
    1039     1291901 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1040      538397 :         else                     return add_rfrac_scal(y, x);
    1041             :       }
    1042     3039747 :       return add_rfrac(x,y);
    1043     2060709 :     case t_VEC:
    1044     2060709 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1045     2060709 :       return RgV_add(x,y);
    1046     1108574 :     case t_COL:
    1047     1108574 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1048     1108574 :       return RgC_add(x,y);
    1049      671356 :     case t_MAT:
    1050      671356 :       lx = lg(x);
    1051      671356 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1052      671356 :       if (lx == 1) return cgetg(1, t_MAT);
    1053      671356 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1054      671349 :       return RgM_add(x,y);
    1055             : 
    1056           0 :     default: pari_err_TYPE2("+",x,y);
    1057             :   }
    1058             :   /* tx != ty */
    1059   837467418 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1060             : 
    1061   837467418 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1062             :   {
    1063   740730611 :     case t_INT:
    1064             :       switch(ty)
    1065             :       {
    1066   436536458 :         case t_REAL: return addir(x,y);
    1067     2139209 :         case t_INTMOD:
    1068     2139209 :           z = cgetg(3, t_INTMOD);
    1069     2139209 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1070    45115176 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1071    45115177 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1072    45115135 :           gel(z,2) = icopy(gel(y,2)); return z;
    1073   250029984 :         case t_COMPLEX: return addRc(x, y);
    1074     5596860 :         case t_PADIC:
    1075     5596860 :           if (!signe(x)) return gcopy(y);
    1076     4533135 :           return addQp(x,y);
    1077        1113 :         case t_QUAD: return addRq(x, y);
    1078     1357772 :         case t_FFELT: return FF_Z_add(y,x);
    1079             :       }
    1080             : 
    1081             :     case t_REAL:
    1082    57094310 :       switch(ty)
    1083             :       {
    1084    13516821 :         case t_FRAC:
    1085    13516821 :           if (!signe(gel(y,1))) return rcopy(x);
    1086    13516821 :           if (!signe(x))
    1087             :           {
    1088       12089 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1089       12089 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1090             :           }
    1091    13504732 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1092    13503139 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1093    43577419 :         case t_COMPLEX: return addRc(x, y);
    1094          70 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, realprec(x));
    1095             : 
    1096           0 :         default: pari_err_TYPE2("+",x,y);
    1097             :       }
    1098             : 
    1099       17640 :     case t_INTMOD:
    1100             :       switch(ty)
    1101             :       {
    1102       17507 :         case t_FRAC: { GEN X = gel(x,1);
    1103       17507 :           z = cgetg(3, t_INTMOD);
    1104       17507 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1105       17507 :           return add_intmod_same(z, X, p1, gel(x,2));
    1106             :         }
    1107          14 :         case t_FFELT:
    1108          14 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1109           0 :             pari_err_OP("+",x,y);
    1110          14 :           return FF_Z_add(y,gel(x,2));
    1111          91 :         case t_COMPLEX: return addRc(x, y);
    1112           0 :         case t_PADIC: { GEN X = gel(x,1);
    1113           0 :           z = cgetg(3, t_INTMOD);
    1114           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1115             :         }
    1116          28 :         case t_QUAD: return addRq(x, y);
    1117             :       }
    1118             : 
    1119             :     case t_FRAC:
    1120             :       switch (ty)
    1121             :       {
    1122    10299316 :         case t_COMPLEX: return addRc(x, y);
    1123      229736 :         case t_PADIC:
    1124      229736 :           if (!signe(gel(x,1))) return gcopy(y);
    1125      229736 :           return addQp(x,y);
    1126         133 :         case t_QUAD: return addRq(x, y);
    1127        1337 :         case t_FFELT: return FF_Q_add(y, x);
    1128             :       }
    1129             : 
    1130             :     case t_FFELT:
    1131           0 :       pari_err_TYPE2("+",x,y);
    1132             : 
    1133          35 :     case t_COMPLEX:
    1134             :       switch(ty)
    1135             :       {
    1136          28 :         case t_PADIC:
    1137          28 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1138           7 :         case t_QUAD:
    1139           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1140           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1141             :       }
    1142             : 
    1143             :     case t_PADIC: /* ty == t_QUAD */
    1144           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1145             :   }
    1146             :   /* tx < ty, !is_const_t(y) */
    1147    31816203 :   switch(ty)
    1148             :   {
    1149        7951 :     case t_MAT:
    1150        7951 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1151        7951 :       if (isrationalzero(x)) return gcopy(y);
    1152        7874 :       return RgM_Rg_add(y, x);
    1153      206994 :     case t_COL:
    1154      206994 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1155      206994 :       return RgC_Rg_add(y, x);
    1156     2392279 :     case t_POLMOD: /* is_const_t(tx) in this case */
    1157     2392279 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1158             :   }
    1159    29208979 :   if (is_scalar_t(tx))  {
    1160    26110042 :     if (tx == t_POLMOD)
    1161             :     {
    1162      112774 :       vx = varn(gel(x,1));
    1163      112774 :       vy = gvar(y);
    1164      112774 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1165             :       else
    1166       85544 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1167       27755 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1168             :     }
    1169    25997268 :     return add_scal(y, x, ty);
    1170             :   }
    1171             :   /* x and y are not scalars, ty != t_MAT */
    1172     3098942 :   vx = gvar(x);
    1173     3098947 :   vy = gvar(y);
    1174     3098947 :   if (vx != vy) { /* x or y is treated as a scalar */
    1175       22759 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1176       32417 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1177       32417 :                                 : add_scal(y, x, ty);
    1178             :   }
    1179             :   /* vx = vy */
    1180     3076188 :   switch(tx)
    1181             :   {
    1182     3075691 :     case t_POL:
    1183             :       switch (ty)
    1184             :       {
    1185        5411 :         case t_SER:
    1186        5411 :           if (lg(x) == 2) return gcopy(y);
    1187        5390 :           i = RgX_val(x); if (i == LONG_MAX) i = 0; /* e.g. x = Mod(0,3)*x^0 */
    1188        5390 :           i = lg(y) + valser(y) - i;
    1189        5390 :           if (i < 3) return gcopy(y);
    1190        5383 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1191        5383 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1192        5383 :           return y;
    1193             : 
    1194     3070280 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1195             :       }
    1196           0 :       break;
    1197             : 
    1198         497 :     case t_SER:
    1199         497 :       if (ty == t_RFRAC)
    1200             :       {
    1201             :         long vn, vd;
    1202         497 :         av = avma;
    1203         497 :         vn = gval(gel(y,1), vy);
    1204         497 :         vd = RgX_valrem_inexact(gel(y,2), NULL);
    1205         497 :         if (vd == LONG_MAX) pari_err_INV("gadd", gel(y,2));
    1206             : 
    1207         497 :         l = lg(x) + valser(x) - (vn - vd);
    1208         497 :         if (l < 3) { set_avma(av); return gcopy(x); }
    1209         497 :         return gerepileupto(av, gadd(x, rfrac_to_ser_i(y, l)));
    1210             :       }
    1211           0 :       break;
    1212             :   }
    1213           0 :   pari_err_TYPE2("+",x,y);
    1214             :   return NULL; /* LCOV_EXCL_LINE */
    1215             : }
    1216             : 
    1217             : GEN
    1218   270570293 : gaddsg(long x, GEN y)
    1219             : {
    1220   270570293 :   long ty = typ(y);
    1221             :   GEN z;
    1222             : 
    1223   270570293 :   switch(ty)
    1224             :   {
    1225   124840942 :     case t_INT:  return addsi(x,y);
    1226   119175277 :     case t_REAL: return addsr(x,y);
    1227       12013 :     case t_INTMOD:
    1228       12013 :       z = cgetg(3, t_INTMOD);
    1229       12013 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1230    15247663 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1231    15247663 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1232    15247663 :       gel(z,2) = icopy(gel(y,2)); return z;
    1233     8164231 :     case t_COMPLEX:
    1234     8164231 :       z = cgetg(3, t_COMPLEX);
    1235     8164231 :       gel(z,1) = gaddsg(x, gel(y,1));
    1236     8164231 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1237             : 
    1238     3130167 :     default: return gadd(stoi(x), y);
    1239             :   }
    1240             : }
    1241             : 
    1242             : GEN
    1243     3167486 : gsubsg(long x, GEN y)
    1244             : {
    1245             :   GEN z, a, b;
    1246             :   pari_sp av;
    1247             : 
    1248     3167486 :   switch(typ(y))
    1249             :   {
    1250      276647 :     case t_INT:  return subsi(x,y);
    1251     1224995 :     case t_REAL: return subsr(x,y);
    1252          56 :     case t_INTMOD:
    1253          56 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1254          56 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1255      732333 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1256      732333 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1257      732333 :       gel(z,2) = icopy(gel(y,2)); return z;
    1258      896341 :     case t_COMPLEX:
    1259      896341 :       z = cgetg(3, t_COMPLEX);
    1260      896341 :       gel(z,1) = gsubsg(x, gel(y,1));
    1261      896341 :       gel(z,2) = gneg(gel(y,2)); return z;
    1262             :   }
    1263       37114 :   av = avma;
    1264       37114 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1265             : }
    1266             : 
    1267             : /********************************************************************/
    1268             : /**                                                                **/
    1269             : /**                          SUBTRACTION                           **/
    1270             : /**                                                                **/
    1271             : /********************************************************************/
    1272             : 
    1273             : GEN
    1274  2828345364 : gsub(GEN x, GEN y)
    1275             : {
    1276  2828345364 :   long tx = typ(x), ty = typ(y);
    1277             :   pari_sp av;
    1278             :   GEN z;
    1279  2828345364 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1280             :   {
    1281  2059566534 :     case t_INT: return subii(x,y);
    1282   577246802 :     case t_REAL: return subrr(x,y);
    1283     1158567 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1284     1158567 :       z = cgetg(3,t_INTMOD);
    1285     1158567 :       if (X==Y || equalii(X,Y))
    1286     1158553 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1287          14 :       gel(z,1) = gcdii(X,Y);
    1288          14 :       warn_coercion(X,Y,gel(z,1));
    1289          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1290          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1291             :     }
    1292     5444059 :     case t_FRAC: return addsub_frac(x,y, subii);
    1293   101975263 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1294   102019009 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1295   101915911 :       if (isintzero(gel(z,2)))
    1296             :       {
    1297       21301 :         set_avma((pari_sp)(z+3));
    1298       21301 :         return gsub(gel(x,1),gel(y,1));
    1299             :       }
    1300   101901027 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1301   101879555 :       return z;
    1302     5437807 :     case t_PADIC:
    1303     5437807 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1304     5437807 :       return addsub_pp(x,y, subii);
    1305         168 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1306         168 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1307         168 :       gel(z,1) = ZX_copy(gel(x,1));
    1308         168 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1309         168 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1310      851090 :     case t_POLMOD:
    1311      851090 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1312      851055 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1313          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1314      203212 :     case t_FFELT: return FF_sub(x,y);
    1315    10083912 :     case t_POL: {
    1316    10083912 :       long vx = varn(x);
    1317    10083912 :       long vy = varn(y);
    1318    10083912 :       if (vx != vy) {
    1319       23627 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1320        2261 :         else                     return Rg_RgX_sub(x, y);
    1321             :       }
    1322    10060285 :       return RgX_sub(x, y);
    1323             :     }
    1324      299522 :     case t_VEC:
    1325      299522 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1326      299522 :       return RgV_sub(x,y);
    1327     3567194 :     case t_COL:
    1328     3567194 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1329     3567194 :       return RgC_sub(x,y);
    1330      251754 :     case t_MAT: {
    1331      251754 :       long lx = lg(x);
    1332      251754 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1333      251751 :       if (lx == 1) return cgetg(1, t_MAT);
    1334      172177 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1335      172182 :       return RgM_sub(x,y);
    1336             :     }
    1337     2469828 :     case t_RFRAC: case t_SER: break;
    1338             : 
    1339           0 :     default: pari_err_TYPE2("+",x,y);
    1340             :   }
    1341    63962836 :   av = avma;
    1342    63962836 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1343             : }
    1344             : 
    1345             : /********************************************************************/
    1346             : /**                                                                **/
    1347             : /**                        MULTIPLICATION                          **/
    1348             : /**                                                                **/
    1349             : /********************************************************************/
    1350             : static GEN
    1351      316176 : mul_ser_scal(GEN x, GEN t)
    1352             : {
    1353      316176 :   if (isexactzero(t)) return gmul(Rg_get_0(x), t);
    1354      312928 :   if (isint1(t)) return gcopy(x);
    1355      262360 :   if (ser_isexactzero(x))
    1356             :   {
    1357         378 :     GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gmul(gel(x,2), t), varn(x), 1);
    1358         378 :     setvalser(z, valser(x)); return z;
    1359             :   }
    1360     3442166 :   pari_APPLY_ser(gmul(gel(x,i), t));
    1361             : }
    1362             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1363             :  * [n/d a valid RFRAC]  */
    1364             : static GEN
    1365    10452363 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1366             : {
    1367    10452363 :   pari_sp av = avma;
    1368             :   GEN z;
    1369             : 
    1370    10452363 :   switch(typ(x))
    1371             :   {
    1372          21 :     case t_PADIC:
    1373          21 :       n = gmul(n, x);
    1374          21 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1375          21 :       return gerepileupto(av, gdiv(n,d));
    1376             : 
    1377         518 :     case t_INTMOD: case t_POLMOD:
    1378         518 :       n = gmul(n, x);
    1379         518 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1380         518 :       return gerepileupto(av, gdiv(n,d));
    1381             :   }
    1382    10451824 :   z = gred_rfrac2(x, d);
    1383    10451824 :   n = simplify_shallow(n);
    1384    10451824 :   if (typ(z) == t_RFRAC)
    1385             :   {
    1386     7936098 :     n = gmul(gel(z,1), n);
    1387     7936098 :     d = gel(z,2);
    1388     7936098 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1389           0 :       z = RgX_Rg_div(n, d);
    1390             :     else
    1391     7936098 :       z = gred_rfrac_simple(n, d);
    1392             :   }
    1393             :   else
    1394     2515726 :     z = gmul(z, n);
    1395    10451824 :   return gerepileupto(av, z);
    1396             : }
    1397             : static GEN
    1398   122126642 : mul_scal(GEN y, GEN x, long ty)
    1399             : {
    1400   122126642 :   switch(ty)
    1401             :   {
    1402   113232924 :     case t_POL:
    1403   113232924 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1404   110079444 :       return RgX_Rg_mul(y, x);
    1405      308035 :     case t_SER: return mul_ser_scal(y, x);
    1406     8585675 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1407          14 :     case t_QFB:
    1408          14 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1409             :   }
    1410          12 :   pari_err_TYPE2("*",x,y);
    1411             :   return NULL; /* LCOV_EXCL_LINE */
    1412             : }
    1413             : static GEN
    1414    10542359 : mul_self_scal(GEN x, GEN y)
    1415   643176000 : { pari_APPLY_same(gmul(y,gel(x,i))); }
    1416             : 
    1417             : static GEN
    1418      160886 : mul_gen_rfrac(GEN X, GEN Y)
    1419             : {
    1420      160886 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1421      160886 :   long vx = gvar(X), vy = varn(y2);
    1422      166647 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1423        5761 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1424             : }
    1425             : /* (x1/x2) * (y1/y2) */
    1426             : static GEN
    1427     7907310 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1428             : {
    1429             :   GEN z, X, Y;
    1430     7907310 :   pari_sp av = avma;
    1431             : 
    1432     7907310 :   X = gred_rfrac2(x1, y2);
    1433     7907310 :   Y = gred_rfrac2(y1, x2);
    1434     7907310 :   if (typ(X) == t_RFRAC)
    1435             :   {
    1436     6628116 :     if (typ(Y) == t_RFRAC) {
    1437     6562519 :       x1 = gel(X,1);
    1438     6562519 :       x2 = gel(X,2);
    1439     6562519 :       y1 = gel(Y,1);
    1440     6562519 :       y2 = gel(Y,2);
    1441     6562519 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1442             :     } else
    1443       65597 :       z = mul_gen_rfrac(Y, X);
    1444             :   }
    1445     1279194 :   else if (typ(Y) == t_RFRAC)
    1446       95289 :     z = mul_gen_rfrac(X, Y);
    1447             :   else
    1448     1183905 :     z = gmul(X, Y);
    1449     7907310 :   return gerepileupto(av, z);
    1450             : }
    1451             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1452             : static GEN
    1453      271003 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1454             : {
    1455      271003 :   pari_sp av = avma;
    1456      271003 :   GEN X = gred_rfrac2(x1, y2);
    1457      271003 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1458             :   {
    1459      264423 :     x2 = RgX_mul(gel(X,2), x2);
    1460      264423 :     x1 = gel(X,1);
    1461             :   }
    1462             :   else
    1463        6580 :     x1 = X;
    1464      271003 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1465             : }
    1466             : 
    1467             : /* Mod(y, Y) * x,  assuming x scalar */
    1468             : static GEN
    1469     3427238 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1470     3427238 : { retmkpolmod(gmul(x,y), RgX_copy(Y)); }
    1471             : 
    1472             : /* cf mulqq */
    1473             : static GEN
    1474     5858691 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1475             : {
    1476     5858691 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1477     5858691 :   pari_sp tetpil, av = avma;
    1478     5858691 :   T[1] = x[1];
    1479     5858691 :   p2 = gmul(gel(x,2), gel(y,2));
    1480     5858691 :   p3 = gmul(gel(x,3), gel(y,3));
    1481     5858691 :   p1 = gmul(gneg_i(c),p3);
    1482             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1483     5858691 :   if (typ(b) == t_INT)
    1484             :   {
    1485     5858670 :     if (signe(b))
    1486             :     {
    1487     4284101 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1488     4284101 :       if (is_pm1(b))
    1489             :       {
    1490     4283443 :         if (signe(b) > 0) p3 = gneg(p3);
    1491             :       }
    1492             :       else
    1493         658 :         p3 = gmul(negi(b), p3);
    1494             :     }
    1495             :     else
    1496             :     {
    1497     1574569 :       p3 = gmul(gel(x,2),gel(y,3));
    1498     1574569 :       p4 = gmul(gel(x,3),gel(y,2));
    1499             :     }
    1500             :   }
    1501             :   else
    1502             :   {
    1503          21 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1504          21 :     p3 = gmul(gneg_i(b), p3);
    1505             :   }
    1506     5858691 :   tetpil = avma;
    1507     5858691 :   gel(T,2) = gadd(p2, p1);
    1508     5858691 :   gel(T,3) = gadd(p4, p3);
    1509     5858691 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1510     5858691 :   return normalizepol_lg(T,4);
    1511             : }
    1512             : /* Mod(x,T) * Mod(y,T) */
    1513             : static GEN
    1514     8562573 : mul_polmod_same(GEN T, GEN x, GEN y)
    1515             : {
    1516     8562573 :   GEN z = cgetg(3,t_POLMOD), a;
    1517     8562573 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1518     8562573 :   gel(z,1) = RgX_copy(T);
    1519             :   /* x * y mod T optimised */
    1520     8562573 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1521     7891737 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1522     1587375 :     a = gmul(x, y);
    1523             :   else
    1524             :   {
    1525     6975198 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1526     5853420 :       a = quad_polmod_mul(T, x, y);
    1527             :     else
    1528     1121778 :       a = RgXQ_mul(x, y, gel(z,1));
    1529             :   }
    1530     8562573 :   gel(z,2) = a; return z;
    1531             : }
    1532             : static GEN
    1533      484415 : sqr_polmod(GEN T, GEN x)
    1534             : {
    1535      484415 :   GEN a, z = cgetg(3,t_POLMOD);
    1536      484415 :   gel(z,1) = RgX_copy(T);
    1537      484415 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1538       93745 :     a = gsqr(x);
    1539             :   else
    1540             :   {
    1541      390670 :     pari_sp av = avma;
    1542      390670 :     a = RgXQ_sqr(x, gel(z,1));
    1543      390670 :     a = gerepileupto(av, a);
    1544             :   }
    1545      484415 :   gel(z,2) = a; return z;
    1546             : }
    1547             : /* Mod(x,X) * Mod(y,Y) */
    1548             : static GEN
    1549     2675106 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1550             : {
    1551     2675106 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1552     2675106 :   long vx = varn(X), vy = varn(Y);
    1553     2675106 :   GEN z = cgetg(3,t_POLMOD);
    1554             : 
    1555     2675106 :   if (vx==vy) {
    1556             :     pari_sp av;
    1557          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1558          14 :     warn_coercion(X,Y,gel(z,1));
    1559          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1560          14 :     return z;
    1561             :   }
    1562     2675092 :   if (varncmp(vx, vy) < 0)
    1563      410662 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1564             :   else
    1565     2264430 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1566     2675092 :   gel(z,2) = gmul(x, y); return z;
    1567             : }
    1568             : 
    1569             : #if 0 /* used by 3M only */
    1570             : /* set z = x+y and return 1 if x,y have the same sign
    1571             :  * set z = x-y and return 0 otherwise */
    1572             : static int
    1573             : did_add(GEN x, GEN y, GEN *z)
    1574             : {
    1575             :   long tx = typ(x), ty = typ(y);
    1576             :   if (tx == ty) switch(tx)
    1577             :   {
    1578             :     case t_INT: *z = addii(x,y); return 1;
    1579             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1580             :     case t_REAL:
    1581             :       if (signe(x) == -signe(y))
    1582             :       { *z = subrr(x,y); return 0; }
    1583             :       else
    1584             :       { *z = addrr(x,y); return 1; }
    1585             :   }
    1586             :   if (tx == t_REAL) switch(ty)
    1587             :   {
    1588             :     case t_INT:
    1589             :       if (signe(x) == -signe(y))
    1590             :       { *z = subri(x,y); return 0; }
    1591             :       else
    1592             :       { *z = addri(x,y); return 1; }
    1593             :     case t_FRAC:
    1594             :       if (signe(x) == -signe(gel(y,1)))
    1595             :       { *z = gsub(x,y); return 0; }
    1596             :       else
    1597             :       { *z = gadd(x,y); return 1; }
    1598             :   }
    1599             :   else if (ty == t_REAL) switch(tx)
    1600             :   {
    1601             :     case t_INT:
    1602             :       if (signe(x) == -signe(y))
    1603             :       { *z = subir(x,y); return 0; }
    1604             :       else
    1605             :       { *z = addir(x,y); return 1; }
    1606             :     case t_FRAC:
    1607             :       if (signe(gel(x,1)) == -signe(y))
    1608             :       { *z = gsub(x,y); return 0; }
    1609             :       else
    1610             :       { *z = gadd(x,y); return 1; }
    1611             :   }
    1612             :   *z = gadd(x,y); return 1;
    1613             : }
    1614             : #endif
    1615             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1616             : static GEN
    1617    11682047 : mulcIR(GEN x, GEN y)
    1618             : {
    1619    11682047 :   GEN z = cgetg(3,t_COMPLEX);
    1620    11682284 :   pari_sp av = avma;
    1621    11682284 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1622    11682249 :   gel(z,2) = gmul(y, gel(x,1));
    1623    11682189 :   return z;
    1624             : 
    1625             : }
    1626             : /* x,y COMPLEX */
    1627             : static GEN
    1628   273773748 : mulcc(GEN x, GEN y)
    1629             : {
    1630   273773748 :   GEN xr = gel(x,1), xi = gel(x,2);
    1631   273773748 :   GEN yr = gel(y,1), yi = gel(y,2);
    1632             :   GEN p1, p2, p3, p4, z;
    1633             :   pari_sp tetpil, av;
    1634             : 
    1635   273773748 :   if (isintzero(xr))
    1636             :   {
    1637    15649471 :     if (isintzero(yr)) {
    1638     7122559 :       av = avma;
    1639     7122559 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1640             :     }
    1641     8526755 :     return mulcIR(y, xi);
    1642             :   }
    1643   258115683 :   if (isintzero(yr)) return mulcIR(x, yi);
    1644             : 
    1645   254954831 :   z = cgetg(3,t_COMPLEX); av = avma;
    1646             : #if 0
    1647             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1648             :    * e.g. xr + xi if exponents differ */
    1649             :   if (did_add(xr, xi, &p3))
    1650             :   {
    1651             :     if (did_add(yr, yi, &p4)) {
    1652             :     /* R = xr*yr - xi*yi
    1653             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1654             :       p1 = gmul(xr,yr);
    1655             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1656             :       p3 = gmul(p3, p4);
    1657             :       p4 = gsub(p2, p1);
    1658             :     } else {
    1659             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1660             :      * I = xr*yi + xi*yr */
    1661             :       p1 = gmul(p3,p4);
    1662             :       p3 = gmul(xr,yi);
    1663             :       p4 = gmul(xi,yr);
    1664             :       p2 = gsub(p3, p4);
    1665             :     }
    1666             :   } else {
    1667             :     if (did_add(yr, yi, &p4)) {
    1668             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1669             :       * I = xr*yi +xi*yr */
    1670             :       p1 = gmul(p3,p4);
    1671             :       p3 = gmul(xr,yi);
    1672             :       p4 = gmul(xi,yr);
    1673             :       p2 = gsub(p4, p3);
    1674             :     } else {
    1675             :     /* R = xr*yr - xi*yi
    1676             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1677             :       p3 = gneg( gmul(p3, p4) );
    1678             :       p1 = gmul(xr,yr);
    1679             :       p2 = gmul(xi,yi);
    1680             :       p4 = gadd(p1, p2);
    1681             : 
    1682             :       p2 = gneg(p2);
    1683             :     }
    1684             :   }
    1685             :   tetpil = avma;
    1686             :   gel(z,1) = gadd(p1,p2);
    1687             :   gel(z,2) = gadd(p3,p4);
    1688             : #else
    1689   254964652 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1690             :   { /* 3M formula */
    1691      559286 :     p3 = addii(xr,xi);
    1692      559286 :     p4 = addii(yr,yi);
    1693      559286 :     p1 = mulii(xr,yr);
    1694      559286 :     p2 = mulii(xi,yi);
    1695      559286 :     p3 = mulii(p3,p4);
    1696      559286 :     p4 = addii(p2,p1);
    1697      559286 :     tetpil = avma;
    1698      559286 :     gel(z,1) = subii(p1,p2);
    1699      559286 :     gel(z,2) = subii(p3,p4);
    1700      559286 :     if (!signe(gel(z,2)))
    1701      113225 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1702             :   }
    1703             :   else
    1704             :   { /* naive 4M formula: avoid all loss of accuracy */
    1705   254405366 :     p1 = gmul(xr,yr);
    1706   254337682 :     p2 = gmul(xi,yi);
    1707   254332640 :     p3 = gmul(xr,yi);
    1708   254336605 :     p4 = gmul(xi,yr);
    1709   254348074 :     tetpil = avma;
    1710   254348074 :     gel(z,1) = gsub(p1,p2);
    1711   254189554 :     gel(z,2) = gadd(p3,p4);
    1712   254201759 :     if (isintzero(gel(z,2)))
    1713             :     {
    1714       50575 :       cgiv(gel(z,2));
    1715       50575 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1716             :     }
    1717             :   }
    1718             : #endif
    1719             : 
    1720   254592561 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1721             : }
    1722             : /* x,y PADIC */
    1723             : static GEN
    1724    17502794 : mulpp(GEN x, GEN y) {
    1725    17502794 :   long l = valp(x) + valp(y);
    1726             :   pari_sp av;
    1727             :   GEN z, t;
    1728    17502794 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1729    17502753 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1730    17276085 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1731             : 
    1732    16929495 :   t = (precp(x) > precp(y))? y: x;
    1733    16929495 :   z = cgetp(t); setvalp(z,l); av = avma;
    1734    16929350 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1735    16929297 :   set_avma(av); return z;
    1736             : }
    1737             : /* x,y QUAD */
    1738             : static GEN
    1739        1106 : mulqq(GEN x, GEN y) {
    1740        1106 :   GEN z = cgetg(4,t_QUAD);
    1741        1106 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1742             :   pari_sp av, tetpil;
    1743        1106 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1744             : 
    1745        1106 :   gel(z,1) = ZX_copy(P); av = avma;
    1746        1106 :   p2 = gmul(gel(x,2),gel(y,2));
    1747        1106 :   p3 = gmul(gel(x,3),gel(y,3));
    1748        1106 :   p1 = gmul(gneg_i(c),p3);
    1749             : 
    1750        1106 :   if (signe(b))
    1751         987 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1752             :   else
    1753             :   {
    1754         119 :     p3 = gmul(gel(x,2),gel(y,3));
    1755         119 :     p4 = gmul(gel(x,3),gel(y,2));
    1756             :   }
    1757        1106 :   tetpil = avma;
    1758        1106 :   gel(z,2) = gadd(p2,p1);
    1759        1106 :   gel(z,3) = gadd(p4,p3);
    1760        1106 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1761             : }
    1762             : 
    1763             : GEN
    1764    15627401 : mulcxI(GEN x)
    1765             : {
    1766             :   GEN z;
    1767    15627401 :   switch(typ(x))
    1768             :   {
    1769     2627992 :     case t_INT: case t_REAL: case t_FRAC:
    1770     2627992 :       return mkcomplex(gen_0, x);
    1771    12999899 :     case t_COMPLEX:
    1772    12999899 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1773    12986449 :       z = cgetg(3,t_COMPLEX);
    1774    12986883 :       gel(z,1) = gneg(gel(x,2));
    1775    12987306 :       gel(z,2) = gel(x,1); return z;
    1776          48 :     default:
    1777          48 :       return gmul(gen_I(), x);
    1778             :   }
    1779             : }
    1780             : GEN
    1781     3079362 : mulcxmI(GEN x)
    1782             : {
    1783             :   GEN z;
    1784     3079362 :   switch(typ(x))
    1785             :   {
    1786      116080 :     case t_INT: case t_REAL: case t_FRAC:
    1787      116080 :       return mkcomplex(gen_0, gneg(x));
    1788     2760378 :     case t_COMPLEX:
    1789     2760378 :       if (isintzero(gel(x,1))) return gel(x,2);
    1790     2722192 :       z = cgetg(3,t_COMPLEX);
    1791     2722166 :       gel(z,1) = gel(x,2);
    1792     2722166 :       gel(z,2) = gneg(gel(x,1)); return z;
    1793      202904 :     default:
    1794      202904 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1795             :   }
    1796             : }
    1797             : /* x * I^k */
    1798             : GEN
    1799     5524084 : mulcxpowIs(GEN x, long k)
    1800             : {
    1801     5524084 :   switch (k & 3)
    1802             :   {
    1803     1369183 :     case 1: return mulcxI(x);
    1804     1362136 :     case 2: return gneg(x);
    1805     1363717 :     case 3: return mulcxmI(x);
    1806             :   }
    1807     1429048 :   return x;
    1808             : }
    1809             : 
    1810             : /* caller will assume l > 2 later */
    1811             : static GEN
    1812     5706235 : init_ser(long l, long v, long e)
    1813             : {
    1814     5706235 :   GEN z = cgetg(l, t_SER);
    1815     5706235 :   z[1] = evalvalser(e) | evalvarn(v) | evalsigne(1); return z;
    1816             : }
    1817             : 
    1818             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1819             : static GEN
    1820     5693888 : fill_ser(GEN z, GEN y)
    1821             : {
    1822     5693888 :   long i, lx = lg(z), ly = lg(y); /* lx > 2 */
    1823     5693888 :   if (ly >= lx) {
    1824    25570501 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1825             :   } else {
    1826      335415 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1827      239573 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1828             :   }
    1829             :   /* dangerous case (already normalized), don't use normalizeser */
    1830     5693888 :   if (ly == 3 && !signe(y)) { setsigne(z, 0); return z; }
    1831     5693790 :   return normalizeser(z);
    1832             : }
    1833             : /* assume typ(x) = t_VEC */
    1834             : static int
    1835          98 : is_ext_qfr(GEN x)
    1836          91 : { return lg(x) == 3 && typ(gel(x,1)) == t_QFB && !qfb_is_qfi(gel(x,1))
    1837         189 :                     && typ(gel(x,2)) == t_REAL; }
    1838             : 
    1839             : GEN
    1840  8517538209 : gmul(GEN x, GEN y)
    1841             : {
    1842             :   long tx, ty, lx, ly, vx, vy, i, l;
    1843             :   pari_sp av, tetpil;
    1844             :   GEN z, p1;
    1845             : 
    1846  8517538209 :   if (x == y) return gsqr(x);
    1847  7639324969 :   tx = typ(x); ty = typ(y);
    1848  7639324969 :   if (tx == ty) switch(tx)
    1849             :   {
    1850  3621952084 :     case t_INT: return mulii(x,y);
    1851  1867769735 :     case t_REAL: return mulrr(x,y);
    1852     1780474 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1853     1780474 :       z = cgetg(3,t_INTMOD);
    1854     1780473 :       if (X==Y || equalii(X,Y))
    1855     1780438 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1856          35 :       gel(z,1) = gcdii(X,Y);
    1857          35 :       warn_coercion(X,Y,gel(z,1));
    1858          35 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1859          35 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1860             :     }
    1861    23411066 :     case t_FRAC:
    1862             :     {
    1863    23411066 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1864    23411066 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1865    23411066 :       z=cgetg(3,t_FRAC);
    1866    23411066 :       p1 = gcdii(x1, y2);
    1867    23411065 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1868    23411065 :       p1 = gcdii(x2, y1);
    1869    23411063 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1870    23411063 :       tetpil = avma;
    1871    23411063 :       gel(z,2) = mulii(x2,y2);
    1872    23411065 :       gel(z,1) = mulii(x1,y1);
    1873    23411065 :       fix_frac_if_int_GC(z,tetpil); return z;
    1874             :     }
    1875   264764018 :     case t_COMPLEX: return mulcc(x, y);
    1876    12511323 :     case t_PADIC: return mulpp(x, y);
    1877         882 :     case t_QUAD: return mulqq(x, y);
    1878    14545949 :     case t_FFELT: return FF_mul(x, y);
    1879    11002829 :     case t_POLMOD:
    1880    11002829 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1881     8327723 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1882     2675106 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1883    48071430 :     case t_POL:
    1884    48071430 :       vx = varn(x);
    1885    48071430 :       vy = varn(y);
    1886    48071430 :       if (vx != vy) {
    1887     4854818 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1888     2071778 :         else                     return RgX_Rg_mul(y, x);
    1889             :       }
    1890    43216612 :       return RgX_mul(x, y);
    1891             : 
    1892     4411812 :     case t_SER: {
    1893     4411812 :       vx = varn(x);
    1894     4411812 :       vy = varn(y);
    1895     4411812 :       if (vx != vy) {
    1896        3675 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1897        3675 :         return mul_ser_scal(y, x);
    1898             :       }
    1899     4408137 :       lx = minss(lg(x), lg(y));
    1900     4408137 :       if (lx == 2) return zeroser(vx, valser(x)+valser(y));
    1901     4408123 :       av = avma; z = init_ser(lx, vx, valser(x)+valser(y));
    1902     4408123 :       x = ser2pol_i(x, lx);
    1903     4408123 :       y = ser2pol_i(y, lx);
    1904     4408123 :       y = RgXn_mul(x, y, lx-2);
    1905     4408123 :       return gerepilecopy(av, fill_ser(z,y));
    1906             :     }
    1907          42 :     case t_VEC:
    1908          42 :       if (!is_ext_qfr(x) || !is_ext_qfr(y)) pari_err_TYPE2("*",x,y);
    1909             :     /* fall through, handle extended t_QFB */
    1910       47804 :     case t_QFB: return qfbcomp(x,y);
    1911     6720535 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1912     4114519 :     case t_MAT: return RgM_mul(x, y);
    1913             : 
    1914        1631 :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1915        1631 :       z = cgetg_copy(x, &l);
    1916        1631 :       if (l != lg(y)) break;
    1917       17325 :       for (i=1; i<l; i++)
    1918             :       {
    1919       15694 :         long yi = y[i];
    1920       15694 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1921       15694 :         z[i] = x[yi];
    1922             :       }
    1923        1631 :       return z;
    1924             : 
    1925           0 :     default:
    1926           0 :       pari_err_TYPE2("*",x,y);
    1927             :   }
    1928             :   /* tx != ty */
    1929  1760503430 :   if (is_const_t(ty) && is_const_t(tx))  {
    1930  1626782679 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1931  1626782679 :     switch(tx) {
    1932  1498143681 :     case t_INT:
    1933             :       switch(ty)
    1934             :       {
    1935   987755231 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1936     1224347 :         case t_INTMOD:
    1937     1224347 :           z = cgetg(3, t_INTMOD);
    1938     1224347 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1939    72414469 :         case t_FRAC:
    1940    72414469 :           if (!signe(x)) return gen_0;
    1941    55236262 :           z=cgetg(3,t_FRAC);
    1942    55236503 :           p1 = gcdii(x,gel(y,2));
    1943    55235628 :           if (equali1(p1))
    1944             :           {
    1945    28962939 :             set_avma((pari_sp)z);
    1946    28962921 :             gel(z,2) = icopy(gel(y,2));
    1947    28963020 :             gel(z,1) = mulii(gel(y,1), x);
    1948             :           }
    1949             :           else
    1950             :           {
    1951    26273019 :             x = diviiexact(x,p1); tetpil = avma;
    1952    26272473 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1953    26272251 :             gel(z,1) = mulii(gel(y,1), x);
    1954    26272660 :             fix_frac_if_int_GC(z,tetpil);
    1955             :           }
    1956    55236167 :           return z;
    1957   433864184 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1958     4196326 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1959        1701 :         case t_QUAD: return mulRq(x,y);
    1960     1627534 :         case t_FFELT: return FF_Z_mul(y,x);
    1961             :       }
    1962             : 
    1963             :     case t_REAL:
    1964   121135119 :       switch(ty)
    1965             :       {
    1966    37510900 :         case t_FRAC: return mulrfrac(x, y);
    1967    83624177 :         case t_COMPLEX: return mulRc(x, y);
    1968          21 :         case t_QUAD: return mulqf(y, x, realprec(x));
    1969          21 :         default: pari_err_TYPE2("*",x,y);
    1970             :       }
    1971             : 
    1972        8022 :     case t_INTMOD:
    1973             :       switch(ty)
    1974             :       {
    1975        7133 :         case t_FRAC: { GEN X = gel(x,1);
    1976        7133 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1977        7133 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1978             :         }
    1979          49 :         case t_COMPLEX: return mulRc_direct(x,y);
    1980         427 :         case t_PADIC: { GEN X = gel(x,1);
    1981         427 :           z = cgetg(3, t_INTMOD);
    1982         427 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1983             :         }
    1984          63 :         case t_QUAD: return mulRq(x, y);
    1985         350 :         case t_FFELT:
    1986         350 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1987           0 :             pari_err_OP("*",x,y);
    1988         350 :           return FF_Z_mul(y,gel(x,2));
    1989             :       }
    1990             : 
    1991             :     case t_FRAC:
    1992             :       switch(ty)
    1993             :       {
    1994     5235240 :         case t_COMPLEX: return mulRc(x, y);
    1995     2582750 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    1996         343 :         case t_QUAD:  return mulRq(x, y);
    1997        2051 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    1998             :       }
    1999             : 
    2000             :     case t_FFELT:
    2001          35 :       pari_err_TYPE2("*",x,y);
    2002             : 
    2003          21 :     case t_COMPLEX:
    2004             :       switch(ty)
    2005             :       {
    2006          14 :         case t_PADIC:
    2007          14 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    2008           7 :         case t_QUAD:
    2009           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    2010           7 :           return mulqf(y, x, lx);
    2011             :       }
    2012             : 
    2013             :     case t_PADIC: /* ty == t_QUAD */
    2014          28 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    2015             :     }
    2016             :   }
    2017             : 
    2018   135547999 :   if (is_matvec_t(ty)) switch(tx)
    2019             :   {
    2020      122585 :     case t_VEC:
    2021             :       switch(ty) {
    2022       15540 :         case t_COL: return RgV_RgC_mul(x,y);
    2023      107045 :         case t_MAT: return RgV_RgM_mul(x,y);
    2024             :       }
    2025           0 :       break;
    2026        1687 :     case t_COL:
    2027             :       switch(ty) {
    2028        1687 :         case t_VEC: return RgC_RgV_mul(x,y);
    2029           0 :         case t_MAT: return RgC_RgM_mul(x,y);
    2030             :       }
    2031           0 :       break;
    2032      825143 :     case t_MAT:
    2033             :       switch(ty) {
    2034           0 :         case t_VEC: return RgM_RgV_mul(x,y);
    2035      825143 :         case t_COL: return RgM_RgC_mul(x,y);
    2036             :       }
    2037             :     default:
    2038     7457743 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    2039     7457774 :       return mul_self_scal(y, x);
    2040             :   }
    2041   130588711 :   if (is_matvec_t(tx))
    2042             :   {
    2043     3084134 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2044     3084590 :     return mul_self_scal(x, y);
    2045             :   }
    2046   127505147 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2047             :   /* tx < ty, !ismatvec(x and y) */
    2048             : 
    2049   127505147 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2050     3403914 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2051   124101233 :   if (is_scalar_t(tx)) {
    2052   120809910 :     if (tx == t_POLMOD) {
    2053     3106910 :       vx = varn(gel(x,1));
    2054     3106910 :       vy = gvar(y);
    2055     3106910 :       if (vx != vy) {
    2056     2872620 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2057       23324 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2058             :       }
    2059             :       /* error if ty == t_SER */
    2060      234290 :       av = avma; y = gmod(y, gel(x,1));
    2061      234283 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2062             :     }
    2063   117703000 :     return mul_scal(y, x, ty);
    2064             :   }
    2065             : 
    2066             :   /* x and y are not scalars, nor matvec */
    2067     3291267 :   vx = gvar(x);
    2068     3291286 :   vy = gvar(y);
    2069     3291286 :   if (vx != vy) /* x or y is treated as a scalar */
    2070     2780439 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2071     2780439 :                                 : mul_scal(y, x, ty);
    2072             :   /* vx = vy */
    2073     1872272 :   switch(tx)
    2074             :   {
    2075     1872237 :     case t_POL:
    2076             :       switch (ty)
    2077             :       {
    2078        6741 :         case t_SER:
    2079             :         {
    2080             :           long v;
    2081        6741 :           av = avma; v = RgX_valrem(x, &x);
    2082        6741 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2083        6727 :           v += valser(y); ly = lg(y);
    2084        6727 :           if (ly == 2) { set_avma(av); return zeroser(vx, v); }
    2085        6496 :           if (degpol(x))
    2086             :           {
    2087        2030 :             z = init_ser(ly, vx, v);
    2088        2030 :             x = RgXn_mul(x, ser2pol_i(y, ly), ly-2);
    2089        2030 :             return gerepilecopy(av, fill_ser(z, x));
    2090             :           }
    2091             :           /* take advantage of x = c*t^v */
    2092        4466 :           set_avma(av); y = mul_ser_scal(y, gel(x,2));
    2093        4466 :           setvalser(y, v); return y;
    2094             :         }
    2095             : 
    2096     1865496 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2097             :       }
    2098           0 :       break;
    2099             : 
    2100          35 :     case t_SER:
    2101             :       switch (ty)
    2102             :       {
    2103          35 :         case t_RFRAC:
    2104          35 :           av = avma;
    2105          35 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2106             :       }
    2107           0 :       break;
    2108             :   }
    2109           0 :   pari_err_TYPE2("*",x,y);
    2110             :   return NULL; /* LCOV_EXCL_LINE */
    2111             : }
    2112             : 
    2113             : /* return a nonnormalized result */
    2114             : GEN
    2115      112731 : sqr_ser_part(GEN x, long l1, long l2)
    2116             : {
    2117             :   long i, j, l;
    2118             :   pari_sp av;
    2119             :   GEN Z, z, p1, p2;
    2120             :   long mi;
    2121      112731 :   if (l2 < l1) return zeroser(varn(x), 2*valser(x));
    2122      112717 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2123      112717 :   Z = cgetg(l2-l1+3, t_SER);
    2124      112717 :   Z[1] = evalvalser(2*valser(x)) | evalvarn(varn(x));
    2125      112717 :   z = Z + 2-l1;
    2126      112717 :   x += 2; mi = 0;
    2127      425390 :   for (i=0; i<l1; i++)
    2128             :   {
    2129      312673 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2130             :   }
    2131             : 
    2132      720212 :   for (i=l1; i<=l2; i++)
    2133             :   {
    2134      607495 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2135      607495 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2136     1256134 :     for (j=i-mi; j<=minss(l,mi); j++)
    2137      648639 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2138      607495 :     p1 = gshift(p1,1);
    2139      607495 :     if ((i&1) == 0 && p2[i>>1])
    2140       93585 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2141      607495 :     gel(z,i) = gerepileupto(av,p1);
    2142             :   }
    2143      112717 :   return Z;
    2144             : }
    2145             : 
    2146             : GEN
    2147  1122310403 : gsqr(GEN x)
    2148             : {
    2149             :   long i, lx;
    2150             :   pari_sp av, tetpil;
    2151             :   GEN z, p1, p2, p3, p4;
    2152             : 
    2153  1122310403 :   switch(typ(x))
    2154             :   {
    2155   909804649 :     case t_INT: return sqri(x);
    2156   195239758 :     case t_REAL: return sqrr(x);
    2157      142489 :     case t_INTMOD: { GEN X = gel(x,1);
    2158      142489 :       z = cgetg(3,t_INTMOD);
    2159      142489 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2160      142491 :       gel(z,1) = icopy(X); return z;
    2161             :     }
    2162     4140919 :     case t_FRAC: return sqrfrac(x);
    2163             : 
    2164     8109179 :     case t_COMPLEX:
    2165     8109179 :       if (isintzero(gel(x,1))) {
    2166      238814 :         av = avma;
    2167      238814 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2168             :       }
    2169     7870360 :       z = cgetg(3,t_COMPLEX); av = avma;
    2170     7870351 :       p1 = gadd(gel(x,1),gel(x,2));
    2171     7870078 :       p2 = gsub(gel(x,1), gel(x,2));
    2172     7870040 :       p3 = gmul(gel(x,1),gel(x,2));
    2173     7870195 :       tetpil = avma;
    2174     7870195 :       gel(z,1) = gmul(p1,p2);
    2175     7870300 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2176             : 
    2177        4746 :     case t_PADIC:
    2178        4746 :       z = cgetg(5,t_PADIC);
    2179        4746 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2180        4746 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2181        4746 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2182        4746 :       gel(z,2) = icopy(gel(x,2));
    2183        4746 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2184        4746 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2185        4746 :       return z;
    2186             : 
    2187          70 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2188          70 :       p1 = gel(x,1);
    2189          70 :       gel(z,1) = ZX_copy(p1); av = avma;
    2190          70 :       p2 = gsqr(gel(x,2));
    2191          70 :       p3 = gsqr(gel(x,3));
    2192          70 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2193             : 
    2194          70 :       if (gequal0(gel(p1,3)))
    2195             :       {
    2196           7 :         tetpil = avma;
    2197           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2198           7 :         av = avma;
    2199           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2200           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2201             :       }
    2202             : 
    2203          63 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2204          63 :       tetpil = avma;
    2205          63 :       gel(z,2) = gadd(p2,p4);
    2206          63 :       gel(z,3) = gadd(p1,p3);
    2207          63 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2208             : 
    2209      484415 :     case t_POLMOD:
    2210      484415 :       return sqr_polmod(gel(x,1), gel(x,2));
    2211             : 
    2212     2956533 :     case t_FFELT: return FF_sqr(x);
    2213             : 
    2214     1327674 :     case t_POL: return RgX_sqr(x);
    2215             : 
    2216       35014 :     case t_SER:
    2217       35014 :       lx = lg(x);
    2218       35014 :       if (ser_isexactzero(x)) {
    2219          21 :         GEN z = gcopy(x);
    2220          21 :         setvalser(z, 2*valser(x));
    2221          21 :         return z;
    2222             :       }
    2223       34993 :       if (lx < 40)
    2224       34734 :         return normalizeser( sqr_ser_part(x, 0, lx-3) );
    2225             :       else
    2226             :       {
    2227         259 :         pari_sp av = avma;
    2228         259 :         GEN z = init_ser(lx, varn(x), 2*valser(x));
    2229         259 :         x = ser2pol_i(x, lx);
    2230         259 :         x = RgXn_sqr(x, lx-2);
    2231         259 :         return gerepilecopy(av, fill_ser(z,x));
    2232             :       }
    2233             : 
    2234          70 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2235          70 :       gel(z,1) = gsqr(gel(x,1));
    2236          70 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2237             : 
    2238        1113 :     case t_MAT: return RgM_sqr(x);
    2239          28 :     case t_VEC: if (!is_ext_qfr(x)) pari_err_TYPE2("*",x,x);
    2240             :     /* fall through handle extended t_QFB */
    2241       81919 :     case t_QFB: return qfbsqr(x);
    2242         658 :     case t_VECSMALL:
    2243         658 :       z = cgetg_copy(x, &lx);
    2244       16289 :       for (i=1; i<lx; i++)
    2245             :       {
    2246       15631 :         long xi = x[i];
    2247       15631 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2248       15631 :         z[i] = x[xi];
    2249             :       }
    2250         658 :       return z;
    2251             :   }
    2252           7 :   pari_err_TYPE2("*",x,x);
    2253             :   return NULL; /* LCOV_EXCL_LINE */
    2254             : }
    2255             : 
    2256             : /********************************************************************/
    2257             : /**                                                                **/
    2258             : /**                           DIVISION                             **/
    2259             : /**                                                                **/
    2260             : /********************************************************************/
    2261             : static GEN
    2262       13169 : div_rfrac_scal(GEN x, GEN y)
    2263             : {
    2264       13169 :   pari_sp av = avma;
    2265       13169 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2266       13169 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2267             : }
    2268             : static GEN
    2269       37515 : div_scal_rfrac(GEN x, GEN y)
    2270             : {
    2271       37515 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2272       37515 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2273             :   {
    2274          14 :     if (degpol(y1))
    2275             :     {
    2276          14 :       pari_sp av = avma;
    2277          14 :       GEN _1 = Rg_get_1(x);
    2278          14 :       if (transtype(_1)) y1 = gmul(y1, _1);
    2279          14 :       return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2280             :     }
    2281           0 :     y1 = gel(y1,2);
    2282             :   }
    2283       37501 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2284             : }
    2285             : static GEN
    2286     1186775 : div_rfrac(GEN x, GEN y)
    2287     1186775 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2288             : 
    2289             : /* x != 0 */
    2290             : static GEN
    2291     1338593 : div_ser_scal(GEN x, GEN t)
    2292             : {
    2293     1338593 :   if (ser_isexactzero(x))
    2294             :   {
    2295          28 :     GEN z = scalarser(lg(x) == 2? Rg_get_0(t): gdiv(gel(x,2), t), varn(x), 1);
    2296          28 :     setvalser(z, valser(x)); return z;
    2297             :   }
    2298     6202113 :   pari_APPLY_ser(gdiv(gel(x,i), t));
    2299             : }
    2300             : GEN
    2301        7630 : ser_normalize(GEN x)
    2302             : {
    2303        7630 :   long i, lx = lg(x);
    2304             :   GEN c, z;
    2305        7630 :   if (lx == 2) return x;
    2306        7630 :   c = gel(x,2); if (gequal1(c)) return x;
    2307        7553 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2308      108528 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2309        7553 :   return z;
    2310             : }
    2311             : 
    2312             : /* y != 0 */
    2313             : static GEN
    2314     6684995 : div_T_scal(GEN x, GEN y, long tx) {
    2315     6684995 :   switch(tx)
    2316             :   {
    2317     5335315 :     case t_POL: return RgX_Rg_div(x, y);
    2318     1338586 :     case t_SER: return div_ser_scal(x, y);
    2319       11097 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2320             :   }
    2321           0 :   pari_err_TYPE2("/",x,y);
    2322             :   return NULL; /* LCOV_EXCL_LINE */
    2323             : }
    2324             : 
    2325             : static GEN
    2326     9257890 : div_scal_pol(GEN x, GEN y) {
    2327     9257890 :   long ly = lg(y);
    2328             :   pari_sp av;
    2329             :   GEN _1;
    2330     9257890 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2331     9180196 :   if (isrationalzero(x)) return zeropol(varn(y));
    2332     7130228 :   av = avma;
    2333     7130228 :   _1 = Rg_get_1(x); if (transtype(_1)) y = gmul(y, _1);
    2334     7130228 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2335             : }
    2336             : static GEN
    2337       18550 : div_scal_ser(GEN x, GEN y)
    2338             : {
    2339       18550 :   pari_sp av = avma;
    2340       18550 :   GEN _1 = Rg_get_1(x);
    2341       18550 :   if (transtype(_1)) y = gmul(y, _1);
    2342       18550 :   return gerepileupto(av, gmul(x, ser_inv(y)));
    2343             : }
    2344             : static GEN
    2345     9267125 : div_scal_T(GEN x, GEN y, long ty) {
    2346     9267125 :   switch(ty)
    2347             :   {
    2348     9211823 :     case t_POL: return div_scal_pol(x, y);
    2349       18543 :     case t_SER: return div_scal_ser(x, y);
    2350       36759 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2351             :   }
    2352           0 :   pari_err_TYPE2("/",x,y);
    2353             :   return NULL; /* LCOV_EXCL_LINE */
    2354             : }
    2355             : 
    2356             : /* assume tx = ty = t_SER, same variable vx */
    2357             : static GEN
    2358      771510 : div_ser(GEN x, GEN y, long vx)
    2359             : {
    2360      771510 :   long e, v = valser(x) - valser(y), lx = lg(x), ly = lg(y);
    2361      771510 :   GEN y0 = y, z;
    2362      771510 :   pari_sp av = avma;
    2363             : 
    2364      771510 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2365      771503 :   if (ser_isexactzero(x))
    2366             :   {
    2367       59892 :     if (lx == 2) return zeroser(vx, v);
    2368         147 :     z = scalarser(gmul(gel(x,2),Rg_get_0(y)), varn(x), 1);
    2369         147 :     setvalser(z, v); return z;
    2370             :   }
    2371      711611 :   if (lx < ly) ly = lx;
    2372      711611 :   y = ser2pol_i_normalize(y, ly, &e);
    2373      711611 :   if (e)
    2374             :   {
    2375           7 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2376           7 :     ly -= e; v -= e; if (ly <= 2) pari_err_INV("div_ser", y0);
    2377             :   }
    2378      711611 :   z = init_ser(ly, vx, v);
    2379      711611 :   if (ly == 3)
    2380             :   {
    2381       12347 :     gel(z,2) = gdiv(gel(x,2), gel(y,2));
    2382       12347 :     if (gequal0(gel(z,2))) setsigne(z, 0); /* can happen: Mod(1,2)/Mod(1,3) */
    2383       12347 :     return gerepileupto(av, z);
    2384             :   }
    2385      699264 :   x = ser2pol_i(x, ly);
    2386      699264 :   y = RgXn_div_i(x, y, ly-2);
    2387      699264 :   return gerepilecopy(av, fill_ser(z,y));
    2388             : }
    2389             : /* x,y compatible PADIC */
    2390             : static GEN
    2391    13212137 : divpp(GEN x, GEN y) {
    2392             :   pari_sp av;
    2393             :   long a, b;
    2394             :   GEN z, M;
    2395             : 
    2396    13212137 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2397    13212130 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2398    13211801 :   a = precp(x);
    2399    13211801 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2400    13211801 :   z = cgetg(5, t_PADIC);
    2401    13211718 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2402    13211614 :   gel(z,2) = icopy(gel(x,2));
    2403    13211547 :   gel(z,3) = icopy(M); av = avma;
    2404    13211697 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2405    13212189 :   return z;
    2406             : }
    2407             : static GEN
    2408       36764 : div_polmod_same(GEN T, GEN x, GEN y)
    2409             : {
    2410       36764 :   long v = varn(T);
    2411       36764 :   GEN a, z = cgetg(3, t_POLMOD);
    2412       36764 :   gel(z,1) = RgX_copy(T);
    2413       36764 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2414       11130 :     a = gdiv(x, y);
    2415       25634 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2416         119 :   {
    2417         119 :     pari_sp av = avma;
    2418         119 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2419             :   }
    2420       25515 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2421        5271 :   {
    2422        5271 :     pari_sp av = avma;
    2423        5271 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2424        5271 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2425        5271 :     a = gerepileupto(av, a);
    2426             :   }
    2427             :   else
    2428             :   {
    2429       20244 :     pari_sp av = avma;
    2430       20244 :     a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2431       20244 :     a = gerepileupto(av, a);
    2432             :   }
    2433       36764 :   gel(z,2) = a; return z;
    2434             : }
    2435             : GEN
    2436   432521117 : gdiv(GEN x, GEN y)
    2437             : {
    2438   432521117 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2439             :   pari_sp av, tetpil;
    2440             :   GEN z, p1, p2;
    2441             : 
    2442   432521117 :   if (tx == ty) switch(tx)
    2443             :   {
    2444    88658019 :     case t_INT:
    2445    88658019 :       return Qdivii(x,y);
    2446             : 
    2447    93618880 :     case t_REAL: return divrr(x,y);
    2448       25746 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2449       25746 :       z = cgetg(3,t_INTMOD);
    2450       25746 :       if (X==Y || equalii(X,Y))
    2451       25732 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2452          14 :       gel(z,1) = gcdii(X,Y);
    2453          14 :       warn_coercion(X,Y,gel(z,1));
    2454          14 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2455          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2456             :     }
    2457     4040955 :     case t_FRAC: {
    2458     4040955 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2459     4040955 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2460     4040955 :       z = cgetg(3, t_FRAC);
    2461     4040955 :       p1 = gcdii(x1, y1);
    2462     4040954 :       if (!equali1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2463     4040954 :       p1 = gcdii(x2, y2);
    2464     4040954 :       if (!equali1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2465     4040954 :       tetpil = avma;
    2466     4040954 :       gel(z,2) = mulii(x2,y1);
    2467     4040953 :       gel(z,1) = mulii(x1,y2);
    2468     4040954 :       normalize_frac(z);
    2469     4040954 :       fix_frac_if_int_GC(z,tetpil);
    2470     4040954 :       return z;
    2471             :     }
    2472     9200919 :     case t_COMPLEX:
    2473     9200919 :       if (isintzero(gel(y,1)))
    2474             :       {
    2475      198044 :         y = gel(y,2);
    2476      198044 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2477      121135 :         z = cgetg(3,t_COMPLEX);
    2478      121135 :         gel(z,1) = gdiv(gel(x,2), y);
    2479      121135 :         av = avma;
    2480      121135 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2481      121135 :         return z;
    2482             :       }
    2483     9002871 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, conj_i(y)); tetpil = avma;
    2484     9002886 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2485             : 
    2486      587420 :     case t_PADIC:
    2487      587420 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2488      587420 :       return divpp(x, y);
    2489             : 
    2490         238 :     case t_QUAD:
    2491         238 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2492         224 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, conj_i(y)); tetpil = avma;
    2493         224 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2494             : 
    2495      243512 :     case t_FFELT: return FF_div(x,y);
    2496             : 
    2497       41139 :     case t_POLMOD:
    2498       41139 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2499       36764 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2500             :       else {
    2501        4375 :         av = avma;
    2502        4375 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2503             :       }
    2504       41139 :       return z;
    2505             : 
    2506    18453300 :     case t_POL:
    2507    18453300 :       vx = varn(x);
    2508    18453300 :       vy = varn(y);
    2509    18453300 :       if (vx != vy) {
    2510      102172 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2511       46067 :                             else return div_scal_pol(x, y);
    2512             :       }
    2513    18351128 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2514    18351128 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2515    18228836 :       av = avma;
    2516    18228836 :       return gerepileupto(av, gred_rfrac2(x,y));
    2517             : 
    2518      771531 :     case t_SER:
    2519      771531 :       vx = varn(x);
    2520      771531 :       vy = varn(y);
    2521      771531 :       if (vx != vy) {
    2522          21 :         if (varncmp(vx, vy) < 0)
    2523             :         {
    2524          14 :           if (!signe(y)) pari_err_INV("gdiv",y);
    2525           7 :           return div_ser_scal(x, y);
    2526             :         }
    2527           7 :         return div_scal_ser(x, y);
    2528             :       }
    2529      771510 :       return div_ser(x, y, vx);
    2530     1189463 :     case t_RFRAC:
    2531     1189463 :       vx = varn(gel(x,2));
    2532     1189463 :       vy = varn(gel(y,2));
    2533     1189463 :       if (vx != vy) {
    2534        2688 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2535         756 :                             else return div_scal_rfrac(x, y);
    2536             :       }
    2537     1186775 :       return div_rfrac(x,y);
    2538             : 
    2539          21 :     case t_VEC: /* handle extended t_QFB */
    2540          21 :     case t_QFB: av = avma; return gerepileupto(av, qfbcomp(x, ginv(y)));
    2541             : 
    2542          28 :     case t_MAT:
    2543          28 :       p1 = RgM_div(x,y);
    2544          28 :       if (!p1) pari_err_INV("gdiv",y);
    2545          21 :       return p1;
    2546             : 
    2547           0 :     default: pari_err_TYPE2("/",x,y);
    2548             :   }
    2549             : 
    2550   215695124 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2551             :   {
    2552     3777684 :     long s = signe(x);
    2553     3777684 :     if (!s) {
    2554      520743 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2555      520743 :       switch (ty)
    2556             :       {
    2557      517579 :       default: return gen_0;
    2558          28 :       case t_INTMOD:
    2559          28 :         z = cgetg(3,t_INTMOD);
    2560          28 :         gel(z,1) = icopy(gel(y,1));
    2561          28 :         gel(z,2) = gen_0; return z;
    2562        3136 :       case t_FFELT: return FF_zero(y);
    2563             :       }
    2564             :     }
    2565     3256941 :     if (is_pm1(x)) {
    2566     1307621 :       if (s > 0) return ginv(y);
    2567      221022 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2568             :     }
    2569     1949320 :     switch(ty)
    2570             :     {
    2571      639053 :       case t_REAL: return divir(x,y);
    2572          21 :       case t_INTMOD:
    2573          21 :         z = cgetg(3, t_INTMOD);
    2574          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2575      794523 :       case t_FRAC:
    2576      794523 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2577      794523 :         if (equali1(p1))
    2578             :         {
    2579      273744 :           set_avma((pari_sp)z);
    2580      273744 :           gel(z,2) = icopy(gel(y,1));
    2581      273744 :           gel(z,1) = mulii(gel(y,2), x);
    2582      273744 :           normalize_frac(z);
    2583      273744 :           fix_frac_if_int(z);
    2584             :         }
    2585             :         else
    2586             :         {
    2587      520779 :           x = diviiexact(x,p1); tetpil = avma;
    2588      520779 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2589      520779 :           gel(z,1) = mulii(gel(y,2), x);
    2590      520779 :           normalize_frac(z);
    2591      520779 :           fix_frac_if_int_GC(z,tetpil);
    2592             :         }
    2593      794523 :         return z;
    2594             : 
    2595         462 :       case t_FFELT: return Z_FF_div(x,y);
    2596      513529 :       case t_COMPLEX: return divRc(x,y);
    2597        1505 :       case t_PADIC: return divTp(x, y);
    2598         231 :       case t_QUAD:
    2599         231 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, conj_i(y)); tetpil = avma;
    2600         231 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2601             :     }
    2602             :   }
    2603   211917435 :   if (gequal0(y))
    2604             :   {
    2605          49 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2606          28 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2607             :   }
    2608             : 
    2609   211927313 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2610             :   {
    2611   172564656 :     case t_REAL:
    2612   172564656 :       switch(ty)
    2613             :       {
    2614   170203926 :         case t_INT: return divri(x,y);
    2615      522503 :         case t_FRAC:
    2616      522503 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2617      522503 :           return gerepileuptoleaf(av, z);
    2618     1838281 :         case t_COMPLEX: return divRc(x, y);
    2619          42 :         case t_QUAD: return divfq(x, y, realprec(x));
    2620          18 :         default: pari_err_TYPE2("/",x,y);
    2621             :       }
    2622             : 
    2623         280 :     case t_INTMOD:
    2624             :       switch(ty)
    2625             :       {
    2626         196 :         case t_INT:
    2627         196 :           z = cgetg(3, t_INTMOD);
    2628         196 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2629          28 :         case t_FRAC: { GEN X = gel(x,1);
    2630          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2631          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2632             :         }
    2633           7 :         case t_FFELT:
    2634           7 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2635           0 :             pari_err_OP("/",x,y);
    2636           7 :           return Z_FF_div(gel(x,2),y);
    2637             : 
    2638           0 :         case t_COMPLEX:
    2639           0 :           av = avma;
    2640           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), conj_i(y)));
    2641             : 
    2642          14 :         case t_QUAD:
    2643          14 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,conj_i(y)); tetpil = avma;
    2644          14 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2645             : 
    2646           7 :         case t_PADIC: { GEN X = gel(x,1);
    2647           7 :           z = cgetg(3, t_INTMOD);
    2648           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2649             :         }
    2650          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2651             :       }
    2652             : 
    2653             :     case t_FRAC:
    2654             :       switch(ty)
    2655             :       {
    2656     2097610 :         case t_INT: z = cgetg(3, t_FRAC);
    2657     2097610 :         p1 = gcdii(y,gel(x,1));
    2658     2097610 :         if (equali1(p1))
    2659             :         {
    2660      828157 :           set_avma((pari_sp)z); tetpil = 0;
    2661      828157 :           gel(z,1) = icopy(gel(x,1));
    2662             :         }
    2663             :         else
    2664             :         {
    2665     1269453 :           y = diviiexact(y,p1); tetpil = avma;
    2666     1269453 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2667             :         }
    2668     2097610 :         gel(z,2) = mulii(gel(x,2),y);
    2669     2097610 :         normalize_frac(z);
    2670     2097610 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2671     2097610 :         return z;
    2672             : 
    2673       59516 :         case t_REAL:
    2674       59516 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2675       59516 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2676             : 
    2677           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2678           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2679           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2680             :         }
    2681             : 
    2682          28 :         case t_FFELT: av=avma;
    2683          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2684             : 
    2685         861 :         case t_COMPLEX: return divRc(x, y);
    2686             : 
    2687        2141 :         case t_PADIC:
    2688        2141 :           if (!signe(gel(x,1))) return gen_0;
    2689        2141 :           return divTp(x, y);
    2690             : 
    2691          42 :         case t_QUAD:
    2692          42 :           av=avma; p1=quadnorm(y); p2=gmul(x,conj_i(y)); tetpil=avma;
    2693          42 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2694             :       }
    2695             : 
    2696             :     case t_FFELT:
    2697         161 :       switch (ty)
    2698             :       {
    2699          77 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2700          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2701           7 :         case t_INTMOD:
    2702           7 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2703           0 :             pari_err_OP("/",x,y);
    2704           7 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2705          49 :         default:
    2706          49 :         pari_err_TYPE2("/",x,y);
    2707             :       }
    2708           0 :       break;
    2709             : 
    2710    13507845 :     case t_COMPLEX:
    2711             :       switch(ty)
    2712             :       {
    2713    13507851 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2714           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2715           0 :         case t_PADIC:
    2716           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2717           0 :         case t_QUAD:
    2718           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2719           0 :           return divfq(x, y, lx);
    2720             :       }
    2721             : 
    2722             :     case t_PADIC:
    2723             :       switch(ty)
    2724             :       {
    2725     1205001 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2726     1204910 :           return signe(gel(x,4))? divpT(x, y)
    2727     2409911 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2728             :         }
    2729           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2730           7 :           z = cgetg(3, t_INTMOD);
    2731           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2732             :         }
    2733          14 :         case t_COMPLEX: case t_QUAD:
    2734          14 :           av=avma; p1=gmul(x,conj_i(y)); p2=gnorm(y); tetpil=avma;
    2735          14 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2736             : 
    2737          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2738             :       }
    2739             : 
    2740             :     case t_QUAD:
    2741             :       switch (ty)
    2742             :       {
    2743        1190 :         case t_INT: case t_INTMOD: case t_FRAC:
    2744        1190 :           z = cgetg(4,t_QUAD);
    2745        1190 :           gel(z,1) = ZX_copy(gel(x,1));
    2746        1190 :           gel(z,2) = gdiv(gel(x,2), y);
    2747        1190 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2748          63 :         case t_REAL: return divqf(x, y, realprec(y));
    2749          14 :         case t_PADIC: return divTp(x, y);
    2750           0 :         case t_COMPLEX:
    2751           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2752           0 :           return divqf(x, y, ly);
    2753             :       }
    2754             :   }
    2755    22487661 :   switch(ty) {
    2756      581830 :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2757      581830 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2758          42 :     case t_MAT:
    2759          42 :       av = avma; p1 = RgM_inv(y);
    2760          35 :       if (!p1) pari_err_INV("gdiv",y);
    2761          28 :       return gerepileupto(av, gmul(x, p1));
    2762           0 :     case t_VEC: case t_COL:
    2763             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2764           0 :       pari_err_TYPE2("/",x,y);
    2765             :   }
    2766    21905789 :   switch(tx) {
    2767     3843706 :     case t_VEC: case t_COL: case t_MAT:
    2768    12997008 :       pari_APPLY_same(gdiv(gel(x,i),y));
    2769           0 :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2770           0 :       pari_err_TYPE2("/",x,y);
    2771             :   }
    2772             : 
    2773    18062083 :   vy = gvar(y);
    2774    18063031 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2775      209668 :     vx = varn(X);
    2776      209668 :     if (vx != vy) {
    2777      209094 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2778      208765 :       retmkpolmod(gdiv(gel(x,2), y), RgX_copy(X));
    2779             :     }
    2780             :     /* y is POL, SER or RFRAC */
    2781         574 :     av = avma;
    2782         574 :     switch(ty)
    2783             :     {
    2784           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2785         574 :       default: y = ginvmod(gmod(y,X), X);
    2786             :     }
    2787         567 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2788             :   }
    2789             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2790             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2791             :    * variable NO_VARIABLE, then the operation is incorrect */
    2792    17853363 :   vx = gvar(x);
    2793    17853361 :   if (vx != vy) { /* includes cases where one is scalar */
    2794    15951800 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2795     9266800 :                         else return div_scal_T(x, y, ty);
    2796             :   }
    2797     1901561 :   switch(tx)
    2798             :   {
    2799     1046213 :     case t_POL:
    2800             :       switch(ty)
    2801             :       {
    2802          28 :         case t_SER:
    2803             :         {
    2804          28 :           GEN y0 = y;
    2805             :           long v;
    2806          28 :           av = avma; v = RgX_valrem(x, &x);
    2807          28 :           if (v == LONG_MAX) return gerepileupto(av, Rg_get_0(x));
    2808          14 :           v -= valser(y); ly = lg(y); /* > 2 */
    2809          14 :           y = ser2pol_i_normalize(y, ly, &i);
    2810          14 :           if (i)
    2811             :           {
    2812           7 :             pari_warn(warner,"normalizing a series with 0 leading term");
    2813           7 :             ly -= i; v -= i; if (ly <= 2) pari_err_INV("gdiv", y0);
    2814             :           }
    2815           7 :           z = init_ser(ly, vx, v);
    2816           7 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, ly-2)));
    2817             :         }
    2818             : 
    2819     1046185 :         case t_RFRAC:
    2820             :         {
    2821     1046185 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2822     1046185 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2823        1171 :             return mul_rfrac_scal(y2, y1, x);
    2824     1045014 :           av = avma;
    2825     1045014 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2826             :         }
    2827             :       }
    2828           0 :       break;
    2829             : 
    2830      584324 :     case t_SER:
    2831             :       switch(ty)
    2832             :       {
    2833      584310 :         case t_POL:
    2834             :         {
    2835      584310 :           long v = valser(x);
    2836      584310 :           lx = lg(x);
    2837      584310 :           if (lx == 2) return zeroser(vx, v - RgX_val(y));
    2838      584205 :           av = avma;
    2839      584205 :           x = ser2pol_i(x, lx); v -= RgX_valrem_inexact(y, &y);
    2840      584205 :           z = init_ser(lx, vx, v);
    2841      584205 :           if (!signe(x)) setsigne(z,0);
    2842      584205 :           return gerepilecopy(av, fill_ser(z, RgXn_div(x, y, lx - 2)));
    2843             :         }
    2844          14 :         case t_RFRAC:
    2845          14 :           av = avma;
    2846          14 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2847             :       }
    2848           0 :       break;
    2849             : 
    2850      271003 :     case t_RFRAC:
    2851             :       switch(ty)
    2852             :       {
    2853      271003 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2854           0 :         case t_SER:
    2855           0 :           av = avma;
    2856           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(gel(x,2), y)));
    2857             :       }
    2858           0 :       break;
    2859             :   }
    2860          21 :   pari_err_TYPE2("/",x,y);
    2861             :   return NULL; /* LCOV_EXCL_LINE */
    2862             : }
    2863             : 
    2864             : /********************************************************************/
    2865             : /**                                                                **/
    2866             : /**                     SIMPLE MULTIPLICATION                      **/
    2867             : /**                                                                **/
    2868             : /********************************************************************/
    2869             : GEN
    2870   229511809 : gmulsg(long s, GEN x)
    2871             : {
    2872             :   pari_sp av;
    2873             :   long i;
    2874             :   GEN z;
    2875             : 
    2876   229511809 :   switch(typ(x))
    2877             :   {
    2878   136834323 :     case t_INT:  return mulsi(s,x);
    2879    69368296 :     case t_REAL: return s? mulsr(s,x): gen_0; /* gmul semantic */
    2880      174123 :     case t_INTMOD: { GEN p = gel(x,1);
    2881      174123 :       z = cgetg(3,t_INTMOD);
    2882      174123 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(x,2)), p));
    2883      174122 :       gel(z,1) = icopy(p); return z;
    2884             :     }
    2885      548192 :     case t_FFELT: return FF_Z_mul(x,stoi(s));
    2886     6207351 :     case t_FRAC:
    2887     6207351 :       if (!s) return gen_0;
    2888     6192707 :       z = cgetg(3,t_FRAC);
    2889     6192707 :       i = labs(s); i = ugcd(i, umodiu(gel(x,2), i));
    2890     6192707 :       if (i == 1)
    2891             :       {
    2892     2267753 :         gel(z,2) = icopy(gel(x,2));
    2893     2267753 :         gel(z,1) = mulis(gel(x,1), s);
    2894             :       }
    2895             :       else
    2896             :       {
    2897     3924954 :         gel(z,2) = diviuexact(gel(x,2), (ulong)i);
    2898     3924954 :         gel(z,1) = mulis(gel(x,1), s/i);
    2899     3924954 :         fix_frac_if_int(z);
    2900             :       }
    2901     6192707 :       return z;
    2902             : 
    2903    14201143 :     case t_COMPLEX:
    2904    14201143 :       if (!s) return gen_0;
    2905    13075339 :       z = cgetg(3, t_COMPLEX);
    2906    13075249 :       gel(z,1) = gmulsg(s,gel(x,1));
    2907    13073666 :       gel(z,2) = gmulsg(s,gel(x,2)); return z;
    2908             : 
    2909        1435 :     case t_PADIC:
    2910        1435 :       if (!s) return gen_0;
    2911        1435 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),x), x));
    2912             : 
    2913           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2914           7 :       gel(z,1) = ZX_copy(gel(x,1));
    2915           7 :       gel(z,2) = gmulsg(s,gel(x,2));
    2916           7 :       gel(z,3) = gmulsg(s,gel(x,3)); return z;
    2917             : 
    2918      205408 :     case t_POLMOD:
    2919      205408 :       retmkpolmod(gmulsg(s,gel(x,2)), RgX_copy(gel(x,1)));
    2920             : 
    2921      746910 :     case t_POL:
    2922      746910 :       if (!signe(x)) return RgX_copy(x);
    2923      730579 :       if (!s) return scalarpol(Rg_get_0(x), varn(x));
    2924     3085451 :       pari_APPLY_pol(gmulsg(s,gel(x,i)));
    2925             : 
    2926         182 :     case t_SER:
    2927         182 :       if (ser_isexactzero(x)) return gcopy(x);
    2928         182 :       if (!s) return Rg_get_0(x);
    2929        3864 :       pari_APPLY_ser(gmulsg(s,gel(x,i)));
    2930             : 
    2931           0 :     case t_RFRAC:
    2932           0 :       if (!s) return zeropol(varn(gel(x,2)));
    2933           0 :       if (s == 1) return gcopy(x);
    2934           0 :       if (s == -1) return gneg(x);
    2935           0 :       return mul_rfrac_scal(gel(x,1), gel(x,2), stoi(s));
    2936             : 
    2937     1233506 :     case t_VEC: case t_COL: case t_MAT:
    2938     3892943 :       pari_APPLY_same(gmulsg(s,gel(x,i)));
    2939             :   }
    2940           0 :   pari_err_TYPE("gmulsg",x);
    2941             :   return NULL; /* LCOV_EXCL_LINE */
    2942             : }
    2943             : 
    2944             : GEN
    2945   124338900 : gmulug(ulong s, GEN x)
    2946             : {
    2947             :   pari_sp av;
    2948             :   long i;
    2949             :   GEN z;
    2950             : 
    2951   124338900 :   switch(typ(x))
    2952             :   {
    2953   122223416 :     case t_INT:  return mului(s,x);
    2954     1059371 :     case t_REAL: return s? mulur(s,x): gen_0; /* gmul semantic */
    2955         364 :     case t_INTMOD: { GEN p = gel(x,1);
    2956         364 :       z = cgetg(3,t_INTMOD);
    2957         364 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mului(s,gel(x,2)), p));
    2958         364 :       gel(z,1) = icopy(p); return z;
    2959             :     }
    2960         413 :     case t_FFELT: return FF_Z_mul(x,utoi(s));
    2961      979933 :     case t_FRAC:
    2962      979933 :       if (!s) return gen_0;
    2963      979919 :       z = cgetg(3,t_FRAC);
    2964      979919 :       i = ugcd(s, umodiu(gel(x,2), s));
    2965      979919 :       if (i == 1)
    2966             :       {
    2967      821635 :         gel(z,2) = icopy(gel(x,2));
    2968      821635 :         gel(z,1) = muliu(gel(x,1), s);
    2969             :       }
    2970             :       else
    2971             :       {
    2972      158284 :         gel(z,2) = diviuexact(gel(x,2), i);
    2973      158284 :         gel(z,1) = muliu(gel(x,1), s/i);
    2974      158284 :         fix_frac_if_int(z);
    2975             :       }
    2976      979919 :       return z;
    2977             : 
    2978       30765 :     case t_COMPLEX:
    2979       30765 :       if (!s) return gen_0;
    2980       30765 :       z = cgetg(3, t_COMPLEX);
    2981       30765 :       gel(z,1) = gmulug(s,gel(x,1));
    2982       30765 :       gel(z,2) = gmulug(s,gel(x,2)); return z;
    2983             : 
    2984        7342 :     case t_PADIC:
    2985        7342 :       if (!s) return gen_0;
    2986        7342 :       av = avma; return gerepileupto(av, mulpp(cvtop2(utoi(s),x), x));
    2987             : 
    2988           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2989           0 :       gel(z,1) = ZX_copy(gel(x,1));
    2990           0 :       gel(z,2) = gmulug(s,gel(x,2));
    2991           0 :       gel(z,3) = gmulug(s,gel(x,3)); return z;
    2992             : 
    2993        6783 :     case t_POLMOD:
    2994        6783 :       retmkpolmod(gmulug(s,gel(x,2)), RgX_copy(gel(x,1)));
    2995             : 
    2996       16226 :     case t_POL:
    2997       16226 :       if (!signe(x)) return RgX_copy(x);
    2998       15449 :       if (!s) return scalarpol(Rg_get_0(x), varn(x));
    2999       52843 :       pari_APPLY_pol(gmulug(s,gel(x,i)));
    3000             : 
    3001           0 :     case t_SER:
    3002           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3003           0 :       if (!s) return Rg_get_0(x);
    3004           0 :       pari_APPLY_ser(gmulug(s,gel(x,i)));
    3005             : 
    3006           0 :     case t_RFRAC:
    3007           0 :       if (!s) return zeropol(varn(gel(x,2)));
    3008           0 :       if (s == 1) return gcopy(x);
    3009           0 :       return mul_rfrac_scal(gel(x,1), gel(x,2), utoi(s));
    3010             : 
    3011       14287 :     case t_VEC: case t_COL: case t_MAT:
    3012       74473 :       pari_APPLY_same(gmulug(s,gel(x,i)));
    3013             :   }
    3014           0 :   pari_err_TYPE("gmulsg",x);
    3015             :   return NULL; /* LCOV_EXCL_LINE */
    3016             : }
    3017             : 
    3018             : /********************************************************************/
    3019             : /**                                                                **/
    3020             : /**                       SIMPLE DIVISION                          **/
    3021             : /**                                                                **/
    3022             : /********************************************************************/
    3023             : 
    3024             : GEN
    3025    12932867 : gdivgs(GEN x, long s)
    3026             : {
    3027    12932867 :   long tx = typ(x), i;
    3028             :   pari_sp av;
    3029             :   GEN z;
    3030             : 
    3031    12932867 :   if (!s)
    3032             :   {
    3033           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3034           0 :     pari_err_INV("gdivgs",gen_0);
    3035             :   }
    3036    12932872 :   switch(tx)
    3037             :   {
    3038     1560428 :     case t_INT: return Qdivis(x, s);
    3039     8495331 :     case t_REAL: return divrs(x,s);
    3040             : 
    3041         357 :     case t_INTMOD:
    3042         357 :       z = cgetg(3, t_INTMOD);
    3043         357 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    3044             : 
    3045         735 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    3046             : 
    3047      552123 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3048      552123 :       i = labs(s); i = ugcd(i, umodiu(gel(x,1), i));
    3049      552123 :       if (i == 1)
    3050             :       {
    3051      409302 :         gel(z,2) = mulsi(s, gel(x,2));
    3052      409302 :         gel(z,1) = icopy(gel(x,1));
    3053             :       }
    3054             :       else
    3055             :       {
    3056      142821 :         gel(z,2) = mulsi(s/i, gel(x,2));
    3057      142821 :         gel(z,1) = divis(gel(x,1), i);
    3058             :       }
    3059      552123 :       normalize_frac(z);
    3060      552123 :       fix_frac_if_int(z); return z;
    3061             : 
    3062     1840600 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3063     1840600 :       gel(z,1) = gdivgs(gel(x,1),s);
    3064     1840597 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3065             : 
    3066         133 :     case t_PADIC: /* divpT */
    3067             :     {
    3068         133 :       GEN p = gel(x,2);
    3069         133 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3070         133 :       av = avma;
    3071         133 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    3072             :     }
    3073             : 
    3074          28 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3075          28 :       gel(z,1) = ZX_copy(gel(x,1));
    3076          28 :       gel(z,2) = gdivgs(gel(x,2),s);
    3077          28 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    3078             : 
    3079       37828 :     case t_POLMOD:
    3080       37828 :       retmkpolmod(gdivgs(gel(x,2),s), RgX_copy(gel(x,1)));
    3081             : 
    3082          91 :     case t_RFRAC:
    3083          91 :       if (s == 1) return gcopy(x);
    3084          84 :       else if (s == -1) return gneg(x);
    3085          84 :       return div_rfrac_scal(x, stoi(s));
    3086             : 
    3087      157454 :     case t_POL: pari_APPLY_pol_normalized(gdivgs(gel(x,i),s));
    3088           0 :     case t_SER: pari_APPLY_ser_normalized(gdivgs(gel(x,i),s));
    3089      407593 :     case t_VEC:
    3090             :     case t_COL:
    3091      894977 :     case t_MAT: pari_APPLY_same(gdivgs(gel(x,i),s));
    3092             :   }
    3093           0 :   pari_err_TYPE2("/",x, stoi(s));
    3094             :   return NULL; /* LCOV_EXCL_LINE */
    3095             : }
    3096             : 
    3097             : GEN
    3098    54713126 : gdivgu(GEN x, ulong s)
    3099             : {
    3100    54713126 :   long tx = typ(x), i;
    3101             :   pari_sp av;
    3102             :   GEN z;
    3103             : 
    3104    54713126 :   if (!s)
    3105             :   {
    3106           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3107           0 :     pari_err_INV("gdivgu",gen_0);
    3108             :   }
    3109    54713031 :   switch(tx)
    3110             :   {
    3111    17915716 :     case t_INT: return Qdiviu(x, s);
    3112    12954304 :     case t_REAL: return divru(x,s);
    3113             : 
    3114      210315 :     case t_INTMOD:
    3115      210315 :       z = cgetg(3, t_INTMOD); s = umodui(s, gel(x,1));
    3116      210315 :       return div_intmod_same(z, gel(x,1), gel(x,2), utoi(s));
    3117             : 
    3118         308 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,utoi(s));
    3119             : 
    3120      641968 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3121      641968 :       i = ugcd(s, umodiu(gel(x,1), s));
    3122      641968 :       if (i == 1)
    3123             :       {
    3124      462120 :         gel(z,2) = mului(s, gel(x,2));
    3125      462120 :         gel(z,1) = icopy(gel(x,1));
    3126             :       }
    3127             :       else
    3128             :       {
    3129      179848 :         gel(z,2) = mului(s/i, gel(x,2));
    3130      179848 :         gel(z,1) = divis(gel(x,1), i);
    3131             :       }
    3132      641968 :       normalize_frac(z);
    3133      641968 :       fix_frac_if_int(z); return z;
    3134             : 
    3135    11383806 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3136    11383806 :       gel(z,1) = gdivgu(gel(x,1),s);
    3137    11383806 :       gel(z,2) = gdivgu(gel(x,2),s); return z;
    3138             : 
    3139    11551909 :     case t_PADIC: /* divpT */
    3140             :     {
    3141    11551909 :       GEN p = gel(x,2);
    3142    11551909 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3143    11416250 :       av = avma;
    3144    11416250 :       return gerepileupto(av, divpp(x, cvtop2(utoi(s),x)));
    3145             :     }
    3146             : 
    3147           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3148           0 :       gel(z,1) = ZX_copy(gel(x,1));
    3149           0 :       gel(z,2) = gdivgu(gel(x,2),s);
    3150           0 :       gel(z,3) = gdivgu(gel(x,3),s); return z;
    3151             : 
    3152        1456 :     case t_POLMOD:
    3153        1456 :       retmkpolmod(gdivgu(gel(x,2),s), RgX_copy(gel(x,1)));
    3154             : 
    3155          56 :     case t_RFRAC:
    3156          56 :       if (s == 1) return gcopy(x);
    3157          56 :       return div_rfrac_scal(x, utoi(s));
    3158             : 
    3159      205660 :     case t_POL: pari_APPLY_pol_normalized(gdivgu(gel(x,i),s));
    3160       16730 :     case t_SER: pari_APPLY_ser_normalized(gdivgu(gel(x,i),s));
    3161         329 :     case t_VEC:
    3162             :     case t_COL:
    3163        1148 :     case t_MAT: pari_APPLY_same(gdivgu(gel(x,i),s));
    3164             :   }
    3165           0 :   pari_err_TYPE2("/",x, utoi(s));
    3166             :   return NULL; /* LCOV_EXCL_LINE */
    3167             : }
    3168             : 
    3169             : /* x / (i*(i+1)) */
    3170             : GEN
    3171   224328053 : divrunextu(GEN x, ulong i)
    3172             : {
    3173   224328053 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3174           0 :     return divri(x, muluu(i , i+1));
    3175             :   else
    3176   224328053 :     return divru(x, i*(i+1));
    3177             : }
    3178             : /* x / (i*(i+1)) */
    3179             : GEN
    3180      814552 : gdivgunextu(GEN x, ulong i)
    3181             : {
    3182      814552 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
    3183           0 :     return gdivgu(x, i*(i+1));
    3184             :   else
    3185      814552 :     return gdiv(x, muluu(i, i+1));
    3186             : }
    3187             : 
    3188             : /* True shift (exact multiplication by 2^n) */
    3189             : GEN
    3190   124662497 : gmul2n(GEN x, long n)
    3191             : {
    3192             :   GEN z, a, b;
    3193             :   long k, l;
    3194             : 
    3195   124662497 :   switch(typ(x))
    3196             :   {
    3197    31200579 :     case t_INT:
    3198    31200579 :       if (n>=0) return shifti(x,n);
    3199     5035047 :       if (!signe(x)) return gen_0;
    3200     3341649 :       l = vali(x); n = -n;
    3201     3341710 :       if (n<=l) return shifti(x,-n);
    3202      398367 :       z = cgetg(3,t_FRAC);
    3203      398367 :       gel(z,1) = shifti(x,-l);
    3204      398367 :       gel(z,2) = int2n(n-l); return z;
    3205             : 
    3206    62786365 :     case t_REAL:
    3207    62786365 :       return shiftr(x,n);
    3208             : 
    3209      180971 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3210      180971 :       z = cgetg(3,t_INTMOD);
    3211      180971 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3212       82810 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3213       82808 :       gel(z,1) = icopy(b); return z;
    3214             : 
    3215      217409 :     case t_FFELT: return FF_mul2n(x,n);
    3216             : 
    3217     1268243 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3218     1268243 :       l = vali(a);
    3219     1268243 :       k = vali(b);
    3220     1268243 :       if (n+l >= k)
    3221             :       {
    3222      395186 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3223      302676 :         l = n-k; k = -k;
    3224             :       }
    3225             :       else
    3226             :       {
    3227      873057 :         k = -(l+n); l = -l;
    3228             :       }
    3229     1175733 :       z = cgetg(3,t_FRAC);
    3230     1175733 :       gel(z,1) = shifti(a,l);
    3231     1175733 :       gel(z,2) = shifti(b,k); return z;
    3232             : 
    3233    10726654 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3234    10726650 :       gel(z,1) = gmul2n(gel(x,1),n);
    3235    10726646 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3236             : 
    3237         105 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3238         105 :       gel(z,1) = ZX_copy(gel(x,1));
    3239         105 :       gel(z,2) = gmul2n(gel(x,2),n);
    3240         105 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3241             : 
    3242      193970 :     case t_POLMOD:
    3243      193970 :       retmkpolmod(gmul2n(gel(x,2),n), RgX_copy(gel(x,1)));
    3244             : 
    3245     1717801 :     case t_POL:
    3246     9105395 :       pari_APPLY_pol(gmul2n(gel(x,i),n));
    3247      102847 :     case t_SER:
    3248      102847 :       if (ser_isexactzero(x)) return gcopy(x);
    3249      634256 :       pari_APPLY_ser(gmul2n(gel(x,i),n));
    3250    16263790 :     case t_VEC: case t_COL: case t_MAT:
    3251    59916568 :       pari_APPLY_same(gmul2n(gel(x,i),n));
    3252             : 
    3253          21 :     case t_RFRAC: /* int2n wrong if n < 0 */
    3254          21 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3255             : 
    3256        5579 :     case t_PADIC: /* int2n wrong if n < 0 */
    3257        5579 :       return gmul(gmul2n(gen_1,n),x);
    3258             :   }
    3259           0 :   pari_err_TYPE("gmul2n",x);
    3260             :   return NULL; /* LCOV_EXCL_LINE */
    3261             : }
    3262             : 
    3263             : /*******************************************************************/
    3264             : /*                                                                 */
    3265             : /*                              INVERSE                            */
    3266             : /*                                                                 */
    3267             : /*******************************************************************/
    3268             : static GEN
    3269      209317 : inv_polmod(GEN T, GEN x)
    3270             : {
    3271      209317 :   GEN z = cgetg(3,t_POLMOD), a;
    3272      209317 :   gel(z,1) = RgX_copy(T);
    3273      209316 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3274       84683 :     a = ginv(x);
    3275             :   else
    3276             :   {
    3277      124633 :     if (lg(T) == 5) /* quadratic fields */
    3278       13055 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3279             :     else
    3280      111578 :       a = RgXQ_inv(x, T);
    3281             :   }
    3282      209316 :   gel(z,2) = a; return z;
    3283             : }
    3284             : 
    3285             : GEN
    3286    36432205 : ginv(GEN x)
    3287             : {
    3288             :   long s;
    3289             :   pari_sp av, tetpil;
    3290             :   GEN z, y, p1, p2;
    3291             : 
    3292    36432205 :   switch(typ(x))
    3293             :   {
    3294     9452538 :     case t_INT:
    3295     9452538 :       if (is_pm1(x)) return icopy(x);
    3296     7897626 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3297     7897612 :       z = cgetg(3,t_FRAC);
    3298     7904062 :       gel(z,1) = s<0? gen_m1: gen_1;
    3299     7904062 :       gel(z,2) = absi(x); return z;
    3300             : 
    3301    10741412 :     case t_REAL: return invr(x);
    3302             : 
    3303       12964 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3304       12964 :       gel(z,1) = icopy(gel(x,1));
    3305       12964 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3306             : 
    3307      552048 :     case t_FRAC: {
    3308      552048 :       GEN a = gel(x,1), b = gel(x,2);
    3309      552048 :       s = signe(a);
    3310      552048 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3311      248266 :       z = cgetg(3,t_FRAC);
    3312      248266 :       gel(z,1) = icopy(b);
    3313      248266 :       gel(z,2) = icopy(a);
    3314      248265 :       normalize_frac(z); return z;
    3315             :     }
    3316     9602086 :     case t_COMPLEX:
    3317     9602086 :       av=avma;
    3318     9602086 :       p1=cxnorm(x);
    3319     9601102 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3320     9601928 :       tetpil=avma;
    3321     9601928 :       return gerepile(av,tetpil,divcR(p2,p1));
    3322             : 
    3323         273 :     case t_QUAD:
    3324         273 :       av=avma; p1=quadnorm(x); p2=conj_i(x); tetpil=avma;
    3325         273 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3326             : 
    3327      358358 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3328      358358 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3329      358351 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3330      358351 :       gel(z,2) = icopy(gel(x,2));
    3331      358351 :       gel(z,3) = icopy(gel(x,3));
    3332      358351 :       gel(z,4) = Zp_inv(gel(x,4),gel(z,2),precp(x)); return z;
    3333             : 
    3334      209317 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3335       14213 :     case t_FFELT: return FF_inv(x);
    3336     5322842 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3337       26341 :     case t_SER: return ser_inv(x);
    3338        2933 :     case t_RFRAC:
    3339             :     {
    3340        2933 :       GEN n = gel(x,1), d = gel(x,2);
    3341        2933 :       pari_sp av = avma, ltop;
    3342        2933 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3343             : 
    3344        2933 :       n = simplify_shallow(n);
    3345        2933 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3346             :       {
    3347        2933 :         if (gequal1(n)) { set_avma(av); return RgX_copy(d); }
    3348         679 :         ltop = avma;
    3349         679 :         z = RgX_Rg_div(d,n);
    3350             :       } else {
    3351           0 :         ltop = avma;
    3352           0 :         z = cgetg(3,t_RFRAC);
    3353           0 :         gel(z,1) = RgX_copy(d);
    3354           0 :         gel(z,2) = RgX_copy(n);
    3355             :       }
    3356         679 :       stackdummy(av, ltop);
    3357         679 :       return z;
    3358             :     }
    3359             : 
    3360          21 :     case t_VEC: if (!is_ext_qfr(x)) break;
    3361             :     case t_QFB:
    3362       99551 :       return qfbpow(x, gen_m1);
    3363       38820 :     case t_MAT:
    3364       38820 :       y = RgM_inv(x);
    3365       38813 :       if (!y) pari_err_INV("ginv",x);
    3366       38743 :       return y;
    3367          28 :     case t_VECSMALL:
    3368             :     {
    3369          28 :       long i, lx = lg(x)-1;
    3370          28 :       y = zero_zv(lx);
    3371         112 :       for (i=1; i<=lx; i++)
    3372             :       {
    3373          84 :         long xi = x[i];
    3374          84 :         if (xi<1 || xi>lx || y[xi])
    3375           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3376          84 :         y[xi] = i;
    3377             :       }
    3378          28 :       return y;
    3379             :     }
    3380             :   }
    3381           6 :   pari_err_TYPE("inverse",x);
    3382             :   return NULL; /* LCOV_EXCL_LINE */
    3383             : }

Generated by: LCOV version 1.16