Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - gen1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19378-5335519) Lines: 1786 1934 92.3 %
Date: 2016-08-31 06:11:49 Functions: 89 89 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                      GENERIC OPERATIONS                        **/
      17             : /**                         (first part)                           **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /* assume z[1] was created last */
      24             : #define fix_frac_if_int(z) if (is_pm1(gel(z,2)))\
      25             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));
      26             : 
      27             : /* assume z[1] was created last */
      28             : #define fix_frac_if_int_GC(z,tetpil) { if (is_pm1(gel(z,2)))\
      29             :   z = gerepileupto((pari_sp)(z+3), gel(z,1));\
      30             : else\
      31             :   gerepilecoeffssp((pari_sp)z, tetpil, z+1, 2); }
      32             : 
      33             : static void
      34          77 : warn_coercion(GEN x, GEN y, GEN z)
      35             : {
      36          77 :   if (DEBUGLEVEL)
      37          56 :    pari_warn(warner,"coercing quotient rings; moduli %Ps and %Ps -> %Ps",x,y,z);
      38          77 : }
      39             : 
      40             : static long
      41          21 : kro_quad(GEN x, GEN y)
      42             : {
      43          21 :   pari_sp av=avma;
      44          21 :   long k = kronecker(quad_disc(x), y);
      45          21 :   avma = av; return k;
      46             : }
      47             : 
      48             : /* is -1 not a square in Zp, assume p prime */
      49             : INLINE int
      50          21 : Zp_nosquare_m1(GEN p) { return (mod4(p) & 2); /* 2 or 3 mod 4 */ }
      51             : 
      52             : static GEN addsub_pp(GEN x, GEN y, GEN(*op)(GEN,GEN));
      53             : static GEN addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN));
      54             : static GEN mulpp(GEN x, GEN y);
      55             : static GEN divpp(GEN x, GEN y);
      56             : /* Argument codes for inline routines
      57             :  * c: complex, p: padic, q: quad, f: floating point (REAL, some complex)
      58             :  * R: without imaginary part (INT, REAL, INTMOD, FRAC, PADIC if -1 not square)
      59             :  * T: some type (to be converted to PADIC)
      60             :  */
      61             : static GEN
      62    45859079 : addRc(GEN x, GEN y) {
      63    45859079 :   GEN z = cgetg(3,t_COMPLEX);
      64    45859079 :   gel(z,1) = gadd(x,gel(y,1));
      65    45859079 :   gel(z,2) = gcopy(gel(y,2)); return z;
      66             : }
      67             : static GEN
      68    47487296 : mulRc(GEN x, GEN y) {
      69    47487296 :   GEN z = cgetg(3,t_COMPLEX);
      70    47487296 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      71    47487296 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      72             : }
      73             : /* for INTMODs: can't simplify when Re(x) = gen_0 */
      74             : static GEN
      75          35 : mulRc_direct(GEN x, GEN y) {
      76          35 :   GEN z = cgetg(3,t_COMPLEX);
      77          35 :   gel(z,1) = gmul(x,gel(y,1));
      78          35 :   gel(z,2) = gmul(x,gel(y,2)); return z;
      79             : }
      80             : static GEN
      81      347196 : divRc(GEN x, GEN y) {
      82      347196 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      83      347196 :   GEN z = cgetg(3,t_COMPLEX);
      84      347196 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      85      347196 :   gel(z,2) = gmul(mt, gel(y,2));
      86      347196 :   return z;
      87             : }
      88             : static GEN
      89     3517089 : divcR(GEN x, GEN y) {
      90     3517089 :   GEN z = cgetg(3,t_COMPLEX);
      91     3517089 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      92     3517089 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      93             : }
      94             : static GEN
      95         679 : addRq(GEN x, GEN y) {
      96         679 :   GEN z = cgetg(4,t_QUAD);
      97         679 :   gel(z,1) = ZX_copy(gel(y,1));
      98         679 :   gel(z,2) = gadd(x, gel(y,2));
      99         679 :   gel(z,3) = gcopy(gel(y,3)); return z;
     100             : }
     101             : static GEN
     102        1834 : mulRq(GEN x, GEN y) {
     103        1834 :   GEN z = cgetg(4,t_QUAD);
     104        1834 :   gel(z,1) = ZX_copy(gel(y,1));
     105        1834 :   gel(z,2) = gmul(x,gel(y,2));
     106        1834 :   gel(z,3) = gmul(x,gel(y,3)); return z;
     107             : }
     108             : static GEN
     109          21 : addqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     110          21 :   long i = gexpo(x) - gexpo(y);
     111          21 :   if (i > 0) prec += nbits2extraprec( i );
     112          21 :   return gerepileupto(av, gadd(y, quadtofp(x, prec)));
     113             : }
     114             : static GEN
     115     1554570 : mulrfrac(GEN x, GEN y)
     116             : {
     117     1554570 :   pari_sp av = avma;
     118     1554570 :   GEN z, a = gel(y,1), b = gel(y,2);
     119     1554570 :   if (is_pm1(a)) /* frequent special case */
     120             :   {
     121      889359 :     z = divri(x, b);
     122      889359 :     if (signe(a) < 0) togglesign(z);
     123      889359 :     return z;
     124             :   }
     125      665211 :   return gerepileuptoleaf(av, divri(mulri(x,gel(y,1)), gel(y,2)));
     126             : }
     127             : static GEN
     128          14 : mulqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     129          14 :   return gerepileupto(av, gmul(y, quadtofp(x, prec)));
     130             : }
     131             : static GEN
     132          28 : divqf(GEN x, GEN y, long prec) { pari_sp av = avma;
     133          28 :   return gerepileupto(av, gdiv(quadtofp(x,prec), y));
     134             : }
     135             : static GEN
     136           7 : divfq(GEN x, GEN y, long prec) { pari_sp av = avma;
     137           7 :   return gerepileupto(av, gdiv(x, quadtofp(y,prec)));
     138             : }
     139             : /* y PADIC, x + y by converting x to padic */
     140             : static GEN
     141           7 : addTp(GEN x, GEN y) { pari_sp av = avma; GEN z;
     142             : 
     143           7 :   if (!valp(y)) z = cvtop2(x,y);
     144             :   else {
     145           7 :     long l = signe(gel(y,4))? valp(y) + precp(y): valp(y);
     146           7 :     z  = cvtop(x, gel(y,2), l);
     147             :   }
     148           7 :   return gerepileupto(av, addsub_pp(z, y, addii));
     149             : }
     150             : /* y PADIC, x * y by converting x to padic */
     151             : static GEN
     152      112567 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     153      112567 :   return gerepileupto(av, mulpp(cvtop2(x,y), y));
     154             : }
     155             : /* y PADIC, non zero x / y by converting x to padic */
     156             : static GEN
     157          91 : divTp(GEN x, GEN y) { pari_sp av = avma;
     158          91 :   return gerepileupto(av, divpp(cvtop2(x,y), y));
     159             : }
     160             : /* x PADIC, x / y by converting y to padic. Assume x != 0; otherwise y
     161             :  * converted to O(p^e) and division by 0 */
     162             : static GEN
     163        1134 : divpT(GEN x, GEN y) { pari_sp av = avma;
     164        1134 :   return gerepileupto(av, divpp(x, cvtop2(y,x)));
     165             : }
     166             : 
     167             : /* z := Mod(x,X) + Mod(y,X) [ t_INTMOD preallocated ], x,y,X INT, 0 <= x,y < X
     168             :  * clean memory from z on */
     169             : static GEN
     170     1813091 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     171     1813091 :   if (lgefint(X) == 3) {
     172     1777626 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     173     1777626 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     174             :   }
     175             :   else {
     176       35465 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     177       35476 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     178             :   }
     179     1813097 :   gel(z,1) = icopy(X); return z;
     180             : }
     181             : static GEN
     182      257695 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     183      257695 :   if (lgefint(X) == 3) {
     184      237151 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     185      237151 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     186             :   }
     187             :   else {
     188       20544 :     GEN u = subii(x,y); if (signe(u) < 0) u = addii(u, X);
     189       20547 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     190             :   }
     191      257697 :   gel(z,1) = icopy(X); return z;
     192             : }
     193             : /* cf add_intmod_same */
     194             : static GEN
     195     3008810 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     196     3008810 :   if (lgefint(X) == 3) {
     197     2893826 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     198     2893826 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     199             :   }
     200             :   else
     201      114984 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     202     3008834 :   gel(z,1) = icopy(X); return z;
     203             : }
     204             : /* cf add_intmod_same */
     205             : static GEN
     206       29769 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     207             : {
     208       29769 :   if (lgefint(X) == 3) {
     209       21840 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     210       21833 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     211             :   }
     212             :   else
     213        7929 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     214       29763 :   gel(z,1) = icopy(X); return z;
     215             : }
     216             : 
     217             : /*******************************************************************/
     218             : /*                                                                 */
     219             : /*        REDUCTION to IRREDUCIBLE TERMS (t_FRAC/t_RFRAC)          */
     220             : /*                                                                 */
     221             : /* (static routines are not memory clean, but OK for gerepileupto) */
     222             : /*******************************************************************/
     223             : /* Compute the denominator of (1/y) * (n/d) = n/yd, y a "scalar".
     224             :  * Sanity check : avoid (1/2) / (Mod(1,2)*x + 1) "=" 1 / (0 * x + 1) */
     225             : static GEN
     226     9406161 : rfrac_denom_mul_scal(GEN d, GEN y)
     227             : {
     228     9406161 :   GEN D = RgX_Rg_mul(d, y);
     229     9406161 :   if (lg(D) != lg(d))
     230             :   { /* try to generate a meaningful diagnostic */
     231           0 :     D = gdiv(leading_coeff(d), y); /* should fail */
     232           0 :     pari_err_INV("gred_rfrac", y); /* better than nothing */
     233             :   }
     234     9406161 :   return D;
     235             : }
     236             : 
     237             : /* d a t_POL, n a coprime t_POL of same var or "scalar". Not memory clean */
     238             : GEN
     239    54283772 : gred_rfrac_simple(GEN n, GEN d)
     240             : {
     241             :   GEN c, cn, cd, z;
     242    54283772 :   long dd = degpol(d);
     243             : 
     244    54283772 :   if (dd <= 0)
     245             :   {
     246        5796 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     247        5796 :     n = gdiv(n, gel(d,2));
     248        5796 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     249        5796 :     return n;
     250             :   }
     251             : 
     252    54277976 :   cd = content(d);
     253    54277976 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     254    54277976 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     255    54277976 :   if (!gequal1(cd)) {
     256     6548490 :     d = RgX_Rg_div(d,cd);
     257     6548490 :     if (!gequal1(cn))
     258             :     {
     259     1302285 :       if (gequal0(cn)) {
     260          42 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     261           0 :         n = (cn != n)? RgX_Rg_div(n,cd): gdiv(n, cd);
     262           0 :         c = gen_1;
     263             :       } else {
     264     1302243 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     265     1302243 :         c = gdiv(cn,cd);
     266             :       }
     267             :     }
     268             :     else
     269     5246205 :       c = ginv(cd);
     270             :   } else {
     271    47729486 :     if (!gequal1(cn))
     272             :     {
     273     2587347 :       if (gequal0(cn)) {
     274         301 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     275          21 :         c = gen_1;
     276             :       } else {
     277     2587046 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     278     2587046 :         c = cn;
     279             :       }
     280             :     } else {
     281    45142139 :       GEN y = cgetg(3,t_RFRAC);
     282    45142139 :       gel(y,1) = gcopy(n);
     283    45142139 :       gel(y,2) = RgX_copy(d); return y;
     284             :     }
     285             :   }
     286             : 
     287     9135515 :   if (typ(c) == t_POL)
     288             :   {
     289      905329 :     z = c;
     290      944746 :     do { z = content(z); } while (typ(z) == t_POL);
     291      905329 :     cd = denom(z);
     292      905329 :     cn = gmul(c, cd);
     293             :   }
     294             :   else
     295             :   {
     296     8230186 :     cn = numer(c);
     297     8230186 :     cd = denom(c);
     298             :   }
     299     9135515 :   z = cgetg(3,t_RFRAC);
     300     9135515 :   gel(z,1) = gmul(n, cn);
     301     9135515 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     302     9135515 :   return z;
     303             : }
     304             : 
     305             : /* in rare cases x may be a t_POL, after 0/x for instance -> pol_0() */
     306             : static GEN
     307      155400 : fix_rfrac(GEN x, long d)
     308             : {
     309             :   GEN z, N, D;
     310      155400 :   if (!d || typ(x) == t_POL) return x;
     311      153769 :   z = cgetg(3, t_RFRAC);
     312      153769 :   N = gel(x,1);
     313      153769 :   D = gel(x,2);
     314      153769 :   if (d > 0) {
     315      307601 :     gel(z, 1) = (typ(N)==t_POL && varn(N)==varn(D))? RgX_shift(N,d)
     316      153853 :                                                    : monomialcopy(N,d,varn(D));
     317      153699 :     gel(z, 2) = RgX_copy(D);
     318             :   } else {
     319          70 :     gel(z, 1) = gcopy(N);
     320          70 :     gel(z, 2) = RgX_shift(D, -d);
     321             :   }
     322      153769 :   return z;
     323             : }
     324             : 
     325             : /* assume d != 0 */
     326             : static GEN
     327    43370393 : gred_rfrac2_i(GEN n, GEN d)
     328             : {
     329             :   GEN y, z;
     330             :   long v, vd, vn;
     331             : 
     332    43370393 :   n = simplify_shallow(n);
     333    43370393 :   if (isrationalzero(n)) return gcopy(n);
     334    36576184 :   d = simplify_shallow(d);
     335    36576184 :   if (typ(d) != t_POL) return gdiv(n,d);
     336    35396687 :   vd = varn(d);
     337    35396687 :   if (typ(n) != t_POL)
     338             :   {
     339    19890718 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     340    19889311 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     341           0 :     pari_err_BUG("gred_rfrac2_i [incompatible variables]");
     342             :   }
     343    15505969 :   vn = varn(n);
     344    15505969 :   if (varncmp(vd, vn) < 0) return gred_rfrac_simple(n,d);
     345    15368393 :   if (varncmp(vd, vn) > 0) return RgX_Rg_div(n,d);
     346             : 
     347             :   /* now n and d are t_POLs in the same variable */
     348    15208641 :   v = RgX_valrem(n, &n) - RgX_valrem(d, &d);
     349    15208641 :   if (!degpol(d))
     350             :   {
     351    12229972 :     n = RgX_Rg_div(n,gel(d,2));
     352    12229972 :     return v? RgX_mulXn(n,v): n;
     353             :   }
     354             : 
     355             :   /* X does not divide gcd(n,d), deg(d) > 0 */
     356     2978669 :   if (!isinexact(n) && !isinexact(d))
     357             :   {
     358     2978543 :     y = RgX_divrem(n, d, &z);
     359     2978543 :     if (!signe(z)) { cgiv(z); return v? RgX_mulXn(y, v): y; }
     360      155274 :     z = RgX_gcd(d, z);
     361      155274 :     if (degpol(z)) { n = RgX_div(n,z); d = RgX_div(d,z); }
     362             :   }
     363      155400 :   return fix_rfrac(gred_rfrac_simple(n,d), v);
     364             : }
     365             : 
     366             : GEN
     367    18174684 : gred_rfrac2(GEN x1, GEN x2)
     368             : {
     369    18174684 :   pari_sp av = avma;
     370    18174684 :   return gerepileupto(av, gred_rfrac2_i(x1, x2));
     371             : }
     372             : 
     373             : /* x1,x2 t_INT, return x1/x2 in reduced form */
     374             : GEN
     375    18722616 : gred_frac2(GEN x1, GEN x2)
     376             : {
     377    18722616 :   GEN r, y = dvmdii(x1,x2,&r);
     378             :   pari_sp av;
     379             : 
     380    18722616 :   if (r == gen_0) return y; /* gen_0 intended */
     381     5240220 :   av = avma; r = gcdii(x2,r);
     382     5240220 :   if (lgefint(r) == 3)
     383             :   {
     384     5049129 :     ulong rr = r[2];
     385     5049129 :     avma = av;
     386     5049129 :     if (rr == 1) y = mkfraccopy(x1, x2);
     387             :     else
     388             :     {
     389     2626585 :       y = cgetg(3,t_FRAC);
     390     2626585 :       gel(y,1) = diviuexact(x1, rr);
     391     2626585 :       gel(y,2) = diviuexact(x2, rr);
     392             :     }
     393             :   }
     394             :   else
     395             :   { /* rare: r left on stack for efficiency */
     396      191091 :     y = cgetg(3,t_FRAC);
     397      191091 :     gel(y,1) = diviiexact(x1,r);
     398      191091 :     gel(y,2) = diviiexact(x2,r);
     399             :   }
     400     5240220 :   normalize_frac(y); return y;
     401             : }
     402             : 
     403             : /*******************************************************************/
     404             : /*                                                                 */
     405             : /*                          CONJUGATION                            */
     406             : /*                                                                 */
     407             : /*******************************************************************/
     408             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     409             : static GEN
     410       11844 : quad_polmod_conj(GEN x, GEN y)
     411             : {
     412             :   GEN z, u, v, a, b;
     413             :   pari_sp av;
     414       11844 :   if (typ(x) != t_POL) return gcopy(x);
     415       11844 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     416       11844 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     417       11844 :   b = gel(y,3); v = gel(x,2);
     418       11844 :   z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
     419       11844 :   gel(z,2) = gerepileupto(av, gsub(v, gdiv(gmul(u,b), a)));
     420       11844 :   gel(z,3) = gneg(u); return z;
     421             : }
     422             : static GEN
     423       11844 : quad_polmod_norm(GEN x, GEN y)
     424             : {
     425             :   GEN z, u, v, a, b, c;
     426             :   pari_sp av;
     427       11844 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0)
     428           0 :     return gsqr(x);
     429       11844 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     430       11844 :   b = gel(y,3); v = gel(x,2);
     431       11844 :   c = gel(y,2); av = avma;
     432       11844 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     433       11844 :   if (!gequal1(a)) z = gdiv(z, a);
     434       11844 :   return gerepileupto(av, gadd(z, gsqr(v)));
     435             : }
     436             : 
     437             : GEN
     438     1576736 : gconj(GEN x)
     439             : {
     440             :   long lx, i;
     441             :   GEN y;
     442             : 
     443     1576736 :   switch(typ(x))
     444             :   {
     445             :     case t_INT: case t_REAL:
     446       91345 :       return mpcopy(x);
     447             : 
     448             :     case t_INTMOD: case t_FRAC: case t_PADIC:
     449        1414 :       return gcopy(x);
     450             : 
     451             :     case t_COMPLEX:
     452     1474620 :       y = cgetg(3,t_COMPLEX);
     453     1474620 :       gel(y,1) = gcopy(gel(x,1));
     454     1474620 :       gel(y,2) = gneg(gel(x,2));
     455     1474620 :       break;
     456             : 
     457             :     case t_QUAD:
     458         595 :       y = cgetg(4,t_QUAD);
     459         595 :       gel(y,1) = ZX_copy(gel(x,1));
     460        1204 :       gel(y,2) = gequal0(gmael(x,1,3))? gcopy(gel(x,2))
     461         609 :                                     : gadd(gel(x,2), gel(x,3));
     462         595 :       gel(y,3) = gneg(gel(x,3));
     463         595 :       break;
     464             : 
     465             :     case t_POL: case t_SER:
     466        4704 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     467        4704 :       for (i=2; i<lx; i++) gel(y,i) = gconj(gel(x,i));
     468        4704 :       break;
     469             : 
     470             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     471        4058 :       y = cgetg_copy(x, &lx);
     472        4058 :       for (i=1; i<lx; i++) gel(y,i) = gconj(gel(x,i));
     473        4058 :       break;
     474             : 
     475             :     case t_POLMOD:
     476             :     {
     477           0 :       GEN y, X = gel(x,1);
     478           0 :       long d = degpol(X);
     479           0 :       if (d < 2) return RgX_copy(x);
     480           0 :       if (d == 2) {
     481           0 :         y = cgetg(3, t_POLMOD);
     482           0 :         gel(y,1) = RgX_copy(X);
     483           0 :         gel(y,2) = quad_polmod_conj(gel(x,2), X); return y;
     484             :       }
     485             :     }
     486             :     default:
     487           0 :       pari_err_TYPE("gconj",x);
     488           0 :       return NULL; /* not reached */
     489             :   }
     490     1483977 :   return y;
     491             : }
     492             : 
     493             : GEN
     494          42 : conjvec(GEN x,long prec)
     495             : {
     496             :   long lx, s, i;
     497             :   GEN z;
     498             : 
     499          42 :   switch(typ(x))
     500             :   {
     501             :     case t_INT: case t_INTMOD: case t_FRAC:
     502           0 :       return mkcolcopy(x);
     503             : 
     504             :     case t_COMPLEX: case t_QUAD:
     505           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     506             : 
     507             :     case t_FFELT:
     508          28 :       return FF_conjvec(x);
     509             : 
     510             :     case t_VEC: case t_COL:
     511           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     512           0 :       if (lx == 1) return z;
     513           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     514           0 :       s = lgcols(z);
     515           0 :       for (i=2; i<lx; i++)
     516             :       {
     517           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     518           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     519             :       }
     520           0 :       break;
     521             : 
     522             :     case t_POLMOD: {
     523          14 :       GEN T = gel(x,1), r;
     524             :       pari_sp av;
     525             : 
     526          14 :       lx = lg(T);
     527          14 :       if (lx <= 3) return cgetg(1,t_COL);
     528          14 :       x = gel(x,2);
     529          42 :       for (i=2; i<lx; i++)
     530             :       {
     531          35 :         GEN c = gel(T,i);
     532          35 :         switch(typ(c)) {
     533             :           case t_INTMOD: {
     534           7 :             GEN p = gel(c,1);
     535             :             pari_sp av;
     536           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     537           7 :             av = avma;
     538           7 :             T = RgX_to_FpX(T,p);
     539           7 :             x = RgX_to_FpX(x, p);
     540           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     541           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     542           7 :             return gerepileupto(av, z);
     543             :           }
     544             :           case t_INT:
     545          28 :           case t_FRAC: break;
     546           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     547             :         }
     548             :       }
     549           7 :       if (typ(x) != t_POL)
     550             :       {
     551           0 :         if (!is_rational_t(typ(x)))
     552           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     553           0 :         retconst_col(lx-3, gcopy(x));
     554             :       }
     555           7 :       RgX_check_QX(x,"conjvec");
     556           7 :       av = avma;
     557           7 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     558           7 :       r = cleanroots(T,prec);
     559           7 :       z = cgetg(lx-2,t_COL);
     560           7 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     561           7 :       return gerepileupto(av, z);
     562             :     }
     563             : 
     564             :     default:
     565           0 :       pari_err_TYPE("conjvec",x);
     566           0 :       return NULL; /* not reached */
     567             :   }
     568           0 :   return z;
     569             : }
     570             : 
     571             : 
     572             : /********************************************************************/
     573             : /**                                                                **/
     574             : /**                           ADDITION                             **/
     575             : /**                                                                **/
     576             : /********************************************************************/
     577             : /* x, y compatible PADIC, op = add or sub */
     578             : static GEN
     579      949522 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     580             : {
     581      949522 :   pari_sp av = avma;
     582             :   long d,e,r,rx,ry;
     583      949522 :   GEN u, z, mod, p = gel(x,2);
     584             :   int swap;
     585             : 
     586      949522 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     587      949522 :   e = valp(x);
     588      949522 :   r = valp(y); d = r-e;
     589      949522 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     590      949522 :   rx = precp(x);
     591      949522 :   ry = precp(y);
     592      949522 :   if (d) /* v(x) < v(y) */
     593             :   {
     594      404180 :     r = d+ry; z = powiu(p,d);
     595      404180 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     596      404180 :     z = mulii(z,gel(y,4));
     597      404180 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     598             :   }
     599             :   else
     600             :   {
     601             :     long c;
     602      545342 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     603      545342 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     604      545342 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     605             :     {
     606       56784 :       avma = av; return zeropadic(p, e+r);
     607             :     }
     608      488558 :     if (c)
     609             :     {
     610       45479 :       mod = diviiexact(mod, powiu(p,c));
     611       45479 :       r -= c;
     612       45479 :       e += c;
     613             :     }
     614             :   }
     615      892738 :   u = modii(u, mod);
     616      892738 :   avma = av; z = cgetg(5,t_PADIC);
     617      892738 :   z[1] = evalprecp(r) | evalvalp(e);
     618      892738 :   gel(z,2) = icopy(p);
     619      892738 :   gel(z,3) = icopy(mod);
     620      892738 :   gel(z,4) = icopy(u); return z;
     621             : }
     622             : /* Rg_to_Fp(t_FRAC) without GC */
     623             : static GEN
     624        1869 : Q_to_Fp(GEN x, GEN p)
     625        1869 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     626             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     627             : static GEN
     628      155361 : addQp(GEN x, GEN y)
     629             : {
     630      155361 :   pari_sp av = avma;
     631      155361 :   long d, r, e, vy = valp(y), py = precp(y);
     632      155361 :   GEN z, q, mod, u, p = gel(y,2);
     633             : 
     634      155361 :   e = Q_pvalrem(x, p, &x);
     635      155358 :   d = vy - e; r = d + py;
     636      155358 :   if (r <= 0) { avma = av; return gcopy(y); }
     637      153860 :   mod = gel(y,3);
     638      153860 :   u   = gel(y,4);
     639      153860 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     640             : 
     641      153869 :   if (d > 0)
     642             :   {
     643      126372 :     q = powiu(p,d);
     644      126367 :     mod = mulii(mod, q);
     645      126366 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     646      126366 :     u = addii(x,  mulii(u, q));
     647             :   }
     648       27497 :   else if (d < 0)
     649             :   {
     650        2024 :     q = powiu(p,-d);
     651        2024 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     652        2024 :     u = addii(u, mulii(x, q));
     653        2024 :     r = py; e = vy;
     654             :   }
     655             :   else
     656             :   {
     657             :     long c;
     658       25473 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     659       25473 :     u = addii(u, x);
     660       25472 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     661             :     {
     662         749 :       avma = av; return zeropadic(p,e+r);
     663             :     }
     664       24724 :     if (c)
     665             :     {
     666       13825 :       mod = diviiexact(mod, powiu(p,c));
     667       13825 :       r -= c;
     668       13825 :       e += c;
     669             :     }
     670             :   }
     671      153111 :   u = modii(u, mod); avma = av;
     672      153107 :   z = cgetg(5,t_PADIC);
     673      153113 :   z[1] = evalprecp(r) | evalvalp(e);
     674      153113 :   gel(z,2) = icopy(p);
     675      153117 :   gel(z,3) = icopy(mod);
     676      153116 :   gel(z,4) = icopy(u); return z;
     677             : }
     678             : 
     679             : /* Mod(x,X) + Mod(y,X) */
     680             : #define addsub_polmod_same addsub_polmod_scal
     681             : /* Mod(x,X) +/- Mod(y,Y) */
     682             : static GEN
     683         980 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     684             : {
     685         980 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     686         980 :   GEN z = cgetg(3,t_POLMOD);
     687         980 :   long vx = varn(X), vy = varn(Y);
     688         980 :   if (vx==vy) {
     689             :     pari_sp av;
     690          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     691          14 :     warn_coercion(X,Y,gel(z,1));
     692          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     693             :   }
     694         966 :   if (varncmp(vx, vy) < 0)
     695         966 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     696             :   else
     697           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     698         966 :   gel(z,2) = op(x, y); return z;
     699             : }
     700             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     701             : static GEN
     702     2569023 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     703             : {
     704     2569023 :   GEN z = cgetg(3,t_POLMOD);
     705     2569023 :   gel(z,1) = RgX_copy(Y);
     706     2569023 :   gel(z,2) = op(y, x); return z;
     707             : }
     708             : 
     709             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     710             : static GEN
     711      201072 : add_ser_scal(GEN y, GEN x)
     712             : {
     713             :   long i, l, ly, vy;
     714             :   GEN z;
     715             : 
     716      201072 :   if (isrationalzero(x)) return gcopy(y);
     717      196557 :   ly = lg(y);
     718      196557 :   l = valp(y);
     719      196557 :   if (l < 3-ly) return gcopy(y);
     720             :   /* l + ly >= 3 */
     721      196438 :   if (l < 0)
     722             :   {
     723         693 :     z = cgetg(ly,t_SER); z[1] = y[1];
     724         693 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     725         693 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     726         693 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     727         693 :     return z;
     728             :   }
     729      195745 :   vy = varn(y);
     730      195745 :   if (l > 0)
     731             :   {
     732        7938 :     if (ser_isexactzero(y))
     733        5509 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     734        2429 :     y -= l; ly += l;
     735        2429 :     z = cgetg(ly,t_SER);
     736        2429 :     x = gcopy(x);
     737        2429 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     738             :   }
     739             :   else
     740             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     741      187807 :     z = cgetg(ly,t_SER);
     742      187807 :     x = gadd(x, gel(y,2));
     743      187807 :     i = 3;
     744             :   }
     745      190236 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     746      190236 :   gel(z,2) = x;
     747      190236 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     748      190236 :   return gequal0(x)? normalize(z): z;
     749             : }
     750             : static long
     751     2312598 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     752             : /* x,y t_SER in the same variable: x+y */
     753             : static GEN
     754     1156670 : ser_add(GEN x, GEN y)
     755             : {
     756     1156670 :   long i, lx,ly, n = valp(y) - valp(x);
     757             :   GEN z;
     758     1156670 :   if (n < 0) { n = -n; swap(x,y); }
     759             :   /* valp(x) <= valp(y) */
     760     1156670 :   lx = _serprec(x);
     761     1156670 :   if (lx == 2) /* don't lose type information */
     762         742 :     return scalarser(gadd(RgX_get_0(x), RgX_get_0(y)), varn(x), valp(x));
     763     1155928 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     764     1155928 :   if (n)
     765             :   {
     766        8848 :     if (n+2 > lx) return gcopy(x);
     767        8148 :     z = cgetg(ly,t_SER);
     768        8148 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     769        8148 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     770             :   } else {
     771     1147080 :     z = cgetg(ly,t_SER);
     772     1147080 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     773             :   }
     774     1155228 :   z[1] = x[1]; return normalize(z);
     775             : }
     776             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     777             : static GEN
     778     8350089 : add_rfrac_scal(GEN y, GEN x)
     779             : {
     780             :   pari_sp av;
     781             :   GEN n;
     782             : 
     783     8350089 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     784     4625077 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     785     4625077 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     786             : }
     787             : 
     788             : /* x "scalar", ty != t_MAT and non-scalar */
     789             : static GEN
     790    14217356 : add_scal(GEN y, GEN x, long ty)
     791             : {
     792    14217356 :   switch(ty)
     793             :   {
     794     9565071 :     case t_POL: return RgX_Rg_add(y, x);
     795      201072 :     case t_SER: return add_ser_scal(y, x);
     796     4241003 :     case t_RFRAC: return add_rfrac_scal(y, x);
     797           0 :     case t_COL: return RgC_Rg_add(y, x);
     798             :     case t_VEC:
     799      210210 :       if (isintzero(x)) return gcopy(y);
     800         168 :       break;
     801             :   }
     802         168 :   pari_err_TYPE2("+",x,y);
     803           0 :   return NULL; /* not reached */
     804             : }
     805             : 
     806             : static GEN
     807    23489029 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     808             : {
     809    23489029 :   pari_sp av = avma;
     810    23489029 :   GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
     811    23489029 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     812    23489029 :   int s = cmpii(x2, y2);
     813             : 
     814    23489029 :   if (!s)
     815             :   { /* common denominator: (x1 op y1) / x2 */
     816    11040580 :     n = op(x1, y1);
     817    11040580 :     if (!signe(n)) { avma = av; return gen_0; }
     818    10504684 :     d = x2;
     819    10504684 :     q = dvmdii(n, d, &r);
     820    10504684 :     if (r == gen_0) { avma = av; return icopy(q); }
     821     9811722 :     r = gcdii(d, r);
     822     9811722 :     if (!is_pm1(r)) { n = diviiexact(n, r); d = diviiexact(d, r); }
     823     9811722 :     gel(z,1) = icopy_avma(n, (pari_sp)z);
     824     9811722 :     gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     825     9811722 :     avma = (pari_sp)gel(z,2); return z;
     826             :   }
     827    12448449 :   if (s < 0)
     828             :   {
     829     5426703 :     GEN Q = dvmdii(y2, x2, &r);
     830     5426703 :     if (r == gen_0)
     831             :     { /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     832     4388868 :       pari_sp av = avma;
     833     4388868 :       n = op(mulii(Q,x1), y1);
     834     4388868 :       q = dvmdii(n, x2, &r);
     835     4388868 :       if (r == gen_0)
     836             :       {
     837      700122 :         gel(z,1) = gerepileuptoint(av, q);
     838      700122 :         gel(z,2) = Q; return z;
     839             :       }
     840     3688746 :       r = gcdii(x2, r);
     841     3688746 :       if (!is_pm1(r)) { n = diviiexact(n, r); x2 = diviiexact(x2, r); }
     842     3688746 :       d = mulii(x2,Q);
     843     3688746 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     844     3688746 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     845     3688746 :       avma = (pari_sp)gel(z,2); return z;
     846             :     }
     847     1037835 :     delta = gcdii(x2,r);
     848             :   }
     849             :   else
     850             :   {
     851     7021746 :     GEN Q = dvmdii(x2, y2, &r);
     852     7021746 :     if (r == gen_0)
     853             :     { /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     854     5531071 :       pari_sp av = avma;
     855     5531071 :       n = op(x1, mulii(Q,y1));
     856     5531071 :       q = dvmdii(n, y2, &r);
     857     5531071 :       if (r == gen_0)
     858             :       {
     859      112670 :         gel(z,1) = gerepileuptoint(av, q);
     860      112670 :         gel(z,2) = Q; return z;
     861             :       }
     862     5418401 :       r = gcdii(y2, r);
     863     5418401 :       if (!is_pm1(r)) { n = diviiexact(n, r); y2 = diviiexact(y2, r); }
     864     5418401 :       d = mulii(y2,Q);
     865     5418401 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     866     5418401 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     867     5418401 :       avma = (pari_sp)gel(z,2); return z;
     868             :     }
     869     1490675 :     delta = gcdii(y2,r);
     870             :   }
     871             :   /* delta = gcd(x2,y2) */
     872     2528510 :   if (is_pm1(delta))
     873             :   { /* numerator is non-zero */
     874      191187 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     875      191187 :     gel(z,2) = mulii(x2,y2); return z;
     876             :   }
     877     2337323 :   x2 = diviiexact(x2,delta);
     878     2337323 :   y2 = diviiexact(y2,delta);
     879     2337323 :   n = op(mulii(x1,y2), mulii(y1,x2));
     880     2337323 :   if (!signe(n)) { avma = av; return gen_0; }
     881     2337323 :   d = mulii(x2, y2);
     882     2337323 :   q = dvmdii(n, delta, &r);
     883     2337323 :   if (r == gen_0)
     884             :   {
     885      142275 :     if (is_pm1(d)) { avma = av; return icopy(q); }
     886      142275 :     avma = (pari_sp)z;
     887      142275 :     gel(z,2) = icopy(d);
     888      142275 :     gel(z,1) = icopy(q); return z;
     889             :   }
     890     2195048 :   r = gcdii(delta, r);
     891     2195048 :   if (!is_pm1(r))
     892             :   {
     893      786081 :     n     = diviiexact(n, r);
     894      786081 :     delta = diviiexact(delta, r);
     895             :   }
     896     2195048 :   d = mulii(d,delta); avma = (pari_sp)z;
     897     2195048 :   gel(z,1) = icopy(n);
     898     2195048 :   gel(z,2) = icopy(d); return z;
     899             : }
     900             : 
     901             : /* assume x2, y2 are t_POLs in the same variable */
     902             : static GEN
     903     2781571 : add_rfrac(GEN x, GEN y)
     904             : {
     905     2781571 :   pari_sp av = avma;
     906     2781571 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     907     2781571 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     908             : 
     909     2781571 :   delta = RgX_gcd(x2,y2);
     910     2781571 :   if (!degpol(delta))
     911             :   {
     912         616 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     913         616 :     d = RgX_mul(x2, y2);
     914         616 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     915             :   }
     916     2780955 :   x2 = RgX_div(x2,delta);
     917     2780955 :   y2 = RgX_div(y2,delta);
     918     2780955 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     919     2780955 :   if (!signe(n))
     920             :   {
     921      680797 :     n = simplify_shallow(n);
     922      680797 :     if (isrationalzero(n)) return gerepileupto(av, n);
     923           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     924             :   }
     925     2100158 :   if (degpol(n) == 0)
     926     1138189 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     927      961969 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     928      961969 :   if (isexactzero(r))
     929             :   {
     930             :     GEN z;
     931      208991 :     d = RgX_mul(x2, y2);
     932             :     /* "constant" denominator ? */
     933      208991 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     934      208991 :     return gerepileupto(av, z);
     935             :   }
     936      752978 :   r = RgX_gcd(delta, r);
     937      752978 :   if (degpol(r))
     938             :   {
     939      148378 :     n = RgX_div(n, r);
     940      148378 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     941             :   }
     942             :   else
     943      604600 :     d = RgX_mul(gel(x,2), y2);
     944      752978 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     945             : }
     946             : 
     947             : GEN
     948  1680639840 : gadd(GEN x, GEN y)
     949             : {
     950  1680639840 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     951             :   pari_sp av, tetpil;
     952             :   GEN z, p1;
     953             : 
     954  1680639840 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     955             :   {
     956  1095896053 :     case t_INT: return addii(x,y);
     957   325291111 :     case t_REAL: return addrr(x,y);
     958     1346890 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     959     1346890 :       z = cgetg(3,t_INTMOD);
     960     1346878 :       if (X==Y || equalii(X,Y))
     961     1346869 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     962          14 :       gel(z,1) = gcdii(X,Y);
     963          14 :       warn_coercion(X,Y,gel(z,1));
     964          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     965          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     966             :     }
     967    16622488 :     case t_FRAC: return addsub_frac(x,y,addii);
     968    33078792 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     969    33078792 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     970    33078792 :       if (isintzero(gel(z,2)))
     971             :       {
     972         315 :         avma = (pari_sp)(z+3);
     973         315 :         return gadd(gel(x,1),gel(y,1));
     974             :       }
     975    33078477 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     976    33078477 :       return z;
     977             :     case t_PADIC:
     978      934367 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     979      934367 :       return addsub_pp(x,y, addii);
     980         301 :     case t_QUAD: z = cgetg(4,t_QUAD);
     981         301 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     982         301 :       gel(z,1) = ZX_copy(gel(x,1));
     983         301 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     984         301 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     985             :     case t_POLMOD:
     986     1160564 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     987     1159619 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     988         945 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     989    10102651 :     case t_FFELT: return FF_add(x,y);
     990             :     case t_POL:
     991    15411195 :       vx = varn(x);
     992    15411195 :       vy = varn(y);
     993    15411195 :       if (vx != vy) {
     994      806289 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     995       19772 :         else                     return RgX_Rg_add(y, x);
     996             :       }
     997    14604906 :       return RgX_add(x, y);
     998             :     case t_SER:
     999     1153744 :       vx = varn(x);
    1000     1153744 :       vy = varn(y);
    1001     1153744 :       if (vx != vy) {
    1002           0 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1003           0 :         else                     return add_ser_scal(y, x);
    1004             :       }
    1005     1153744 :       return ser_add(x, y);
    1006             :     case t_RFRAC:
    1007     4073474 :       vx = varn(gel(x,2));
    1008     4073474 :       vy = varn(gel(y,2));
    1009     4073474 :       if (vx != vy) {
    1010     1291903 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1011      536802 :         else                     return add_rfrac_scal(y, x);
    1012             :       }
    1013     2781571 :       return add_rfrac(x,y);
    1014             :     case t_VEC:
    1015      472688 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1016      472688 :       return RgV_add(x,y);
    1017             :     case t_COL:
    1018     5057054 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1019     5057054 :       return RgC_add(x,y);
    1020             :     case t_MAT:
    1021      495138 :       lx = lg(x);
    1022      495138 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1023      495138 :       if (lx == 1) return cgetg(1, t_MAT);
    1024      495138 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1025      495131 :       return RgM_add(x,y);
    1026             : 
    1027           0 :     default: pari_err_TYPE2("+",x,y);
    1028             :   }
    1029             :   /* tx != ty */
    1030   169545191 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1031             : 
    1032   169545191 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1033             :   {
    1034             :     case t_INT:
    1035   135431557 :       switch(ty)
    1036             :       {
    1037    81334610 :         case t_REAL: return addir(x,y);
    1038             :         case t_INTMOD:
    1039      462662 :           z = cgetg(3, t_INTMOD);
    1040      462662 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1041    21549671 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1042    21549671 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1043    21549671 :           gel(z,2) = icopy(gel(y,2)); return z;
    1044    30429371 :         case t_COMPLEX: return addRc(x, y);
    1045             :         case t_PADIC:
    1046      163399 :           if (!signe(x)) return gcopy(y);
    1047      152528 :           return addQp(x,y);
    1048         658 :         case t_QUAD: return addRq(x, y);
    1049     1491187 :         case t_FFELT: return FF_Z_add(y,x);
    1050             :       }
    1051             : 
    1052             :     case t_REAL:
    1053    16195931 :       switch(ty)
    1054             :       {
    1055             :         case t_FRAC:
    1056      996726 :           if (!signe(gel(y,1))) return rcopy(x);
    1057      996726 :           if (!signe(x))
    1058             :           {
    1059       11963 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1060       11963 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1061             :           }
    1062      984763 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1063      984763 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1064    15199191 :         case t_COMPLEX: return addRc(x, y);
    1065          14 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1066             : 
    1067           0 :         default: pari_err_TYPE2("+",x,y);
    1068             :       }
    1069             : 
    1070             :     case t_INTMOD:
    1071        4006 :       switch(ty)
    1072             :       {
    1073        3500 :         case t_FRAC: { GEN X = gel(x,1);
    1074        3500 :           z = cgetg(3, t_INTMOD);
    1075        3500 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1076        3500 :           return add_intmod_same(z, X, p1, gel(x,2));
    1077             :         }
    1078             :         case t_FFELT:
    1079         436 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1080           0 :             pari_err_OP("+",x,y);
    1081         436 :           return FF_Z_add(y,gel(x,2));
    1082          56 :         case t_COMPLEX: return addRc(x, y);
    1083           0 :         case t_PADIC: { GEN X = gel(x,1);
    1084           0 :           z = cgetg(3, t_INTMOD);
    1085           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1086             :         }
    1087          14 :         case t_QUAD: return addRq(x, y);
    1088             :       }
    1089             : 
    1090             :     case t_FRAC:
    1091      234738 :       switch (ty)
    1092             :       {
    1093      230454 :         case t_COMPLEX: return addRc(x, y);
    1094             :         case t_PADIC:
    1095        2835 :           if (!signe(gel(x,1))) return gcopy(y);
    1096        2835 :           return addQp(x,y);
    1097           7 :         case t_QUAD: return addRq(x, y);
    1098        1442 :         case t_FFELT: return FF_Q_add(y, x);
    1099             :       }
    1100             : 
    1101             :     case t_FFELT:
    1102           0 :       pari_err_TYPE2("+",x,y);
    1103             : 
    1104             :     case t_COMPLEX:
    1105          21 :       switch(ty)
    1106             :       {
    1107             :         case t_PADIC:
    1108          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1109             :         case t_QUAD:
    1110           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1111           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1112             :       }
    1113             : 
    1114             :     case t_PADIC: /* ty == t_QUAD */
    1115           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1116             :   }
    1117             :   /* tx < ty, !is_const_t(y) */
    1118    17678943 :   switch(ty)
    1119             :   {
    1120             :     case t_MAT:
    1121       17318 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1122       17318 :       if (isrationalzero(x)) return gcopy(y);
    1123       17192 :       return RgM_Rg_add(y, x);
    1124             :     case t_COL:
    1125       39472 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1126       39472 :       return RgC_Rg_add(y, x);
    1127             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1128      580131 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1129             :   }
    1130    17042022 :   if (is_scalar_t(tx))  {
    1131    14218241 :     if (tx == t_POLMOD)
    1132             :     {
    1133       23660 :       vx = varn(gel(x,1));
    1134       23660 :       vy = gvar(y);
    1135       23660 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1136             :       else
    1137       19747 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1138        3976 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1139             :     }
    1140    14194581 :     return add_scal(y, x, ty);
    1141             :   }
    1142             :   /* x and y are not scalars, ty != t_MAT */
    1143     2823781 :   vx = gvar(x);
    1144     2823781 :   vy = gvar(y);
    1145     2823781 :   if (vx != vy) { /* x or y is treated as a scalar */
    1146        3098 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1147        6182 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1148        3091 :                                 : add_scal(y, x, ty);
    1149             :   }
    1150             :   /* vx = vy */
    1151     2820683 :   switch(tx)
    1152             :   {
    1153             :     case t_POL:
    1154     2820130 :       switch (ty)
    1155             :       {
    1156             :         case t_SER:
    1157        2947 :           if (lg(x) == 2) return gcopy(y);
    1158        2933 :           i = lg(y) + valp(y) - RgX_val(x);
    1159        2933 :           if (i < 3) return gcopy(y);
    1160             : 
    1161        2926 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1162        2926 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1163        2926 :           return y;
    1164             : 
    1165     2817183 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1166             :       }
    1167           0 :       break;
    1168             : 
    1169             :     case t_SER:
    1170         553 :       if (ty == t_RFRAC)
    1171             :       {
    1172             :         GEN n, d;
    1173             :         long vn, vd;
    1174         553 :         av = avma;
    1175         553 :         n = gel(y,1); vn = gval(n, vy);
    1176         553 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1177             : 
    1178         553 :         l = lg(x) + valp(x) - (vn - vd);
    1179         553 :         if (l < 3) { avma = av; return gcopy(x); }
    1180             : 
    1181             :         /* take advantage of y = t^n ! */
    1182         553 :         if (degpol(d))
    1183          63 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1184             :         else {
    1185         490 :           y = gdiv(n, gel(d,2));
    1186         490 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1187             :         }
    1188         553 :         setvalp(y, valp(y) - vd);
    1189         553 :         return gerepileupto(av, gadd(y, x));
    1190             :       }
    1191           0 :       break;
    1192             :   }
    1193           0 :   pari_err_TYPE2("+",x,y);
    1194           0 :   return NULL; /* not reached */
    1195             : }
    1196             : 
    1197             : GEN
    1198    15430268 : gaddsg(long x, GEN y)
    1199             : {
    1200    15430268 :   long ty = typ(y);
    1201             :   GEN z;
    1202             : 
    1203    15430268 :   switch(ty)
    1204             :   {
    1205      778615 :     case t_INT:  return addsi(x,y);
    1206     7020387 :     case t_REAL: return addsr(x,y);
    1207             :     case t_INTMOD:
    1208          14 :       z = cgetg(3, t_INTMOD);
    1209          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1210      382984 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1211      382984 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1212      382984 :       gel(z,2) = icopy(gel(y,2)); return z;
    1213             :     case t_COMPLEX:
    1214     6826353 :       z = cgetg(3, t_COMPLEX);
    1215     6826353 :       gel(z,1) = gaddsg(x, gel(y,1));
    1216     6826353 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1217             : 
    1218      421915 :     default: return gadd(stoi(x), y);
    1219             :   }
    1220             : }
    1221             : 
    1222             : GEN
    1223      410617 : gsubsg(long x, GEN y)
    1224             : {
    1225             :   GEN z, a, b;
    1226             :   pari_sp av;
    1227             : 
    1228      410617 :   switch(typ(y))
    1229             :   {
    1230       37251 :     case t_INT:  return subsi(x,y);
    1231      292797 :     case t_REAL: return subsr(x,y);
    1232             :     case t_INTMOD:
    1233          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1234          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1235        7713 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1236        7713 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1237        7713 :       gel(z,2) = icopy(gel(y,2)); return z;
    1238             :     case t_COMPLEX:
    1239       46550 :       z = cgetg(3, t_COMPLEX);
    1240       46550 :       gel(z,1) = gsubsg(x, gel(y,1));
    1241       46550 :       gel(z,2) = gneg(gel(y,2)); return z;
    1242             :   }
    1243       26257 :   av = avma;
    1244       26257 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1245             : }
    1246             : 
    1247             : /********************************************************************/
    1248             : /**                                                                **/
    1249             : /**                          SUBTRACTION                           **/
    1250             : /**                                                                **/
    1251             : /********************************************************************/
    1252             : 
    1253             : GEN
    1254  1468627655 : gsub(GEN x, GEN y)
    1255             : {
    1256  1468627655 :   long tx = typ(x), ty = typ(y);
    1257             :   pari_sp av;
    1258             :   GEN z;
    1259  1468627655 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1260             :   {
    1261  1279985483 :     case t_INT: return subii(x,y);
    1262   140735295 :     case t_REAL: return subrr(x,y);
    1263      257713 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1264      257713 :       z = cgetg(3,t_INTMOD);
    1265      257710 :       if (X==Y || equalii(X,Y))
    1266      257694 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1267          14 :       gel(z,1) = gcdii(X,Y);
    1268          14 :       warn_coercion(X,Y,gel(z,1));
    1269          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1270          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1271             :     }
    1272     6866541 :     case t_FRAC: return addsub_frac(x,y, subii);
    1273    16366829 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1274    16366829 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1275    16366829 :       if (isintzero(gel(z,2)))
    1276             :       {
    1277        6627 :         avma = (pari_sp)(z+3);
    1278        6627 :         return gsub(gel(x,1),gel(y,1));
    1279             :       }
    1280    16360202 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1281    16360202 :       return z;
    1282             :     case t_PADIC:
    1283       15148 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1284       15148 :       return addsub_pp(x,y, subii);
    1285         581 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1286         581 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1287         581 :       gel(z,1) = ZX_copy(gel(x,1));
    1288         581 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1289         581 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1290             :     case t_POLMOD:
    1291      825332 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1292      825297 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1293          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1294      417800 :     case t_FFELT: return FF_sub(x,y);
    1295             :     case t_POL: {
    1296     6245457 :       long vx = varn(x);
    1297     6245457 :       long vy = varn(y);
    1298     6245457 :       if (vx != vy) {
    1299       13756 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1300        2611 :         else                     return Rg_RgX_sub(x, y);
    1301             :       }
    1302     6231701 :       return RgX_sub(x, y);
    1303             :     }
    1304             :     case t_VEC:
    1305       65588 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1306       65588 :       return RgV_sub(x,y);
    1307             :     case t_COL:
    1308     1220973 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1309     1220973 :       return RgC_sub(x,y);
    1310             :     case t_MAT: {
    1311        7570 :       long lx = lg(x);
    1312        7570 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1313        7570 :       if (lx == 1) return cgetg(1, t_MAT);
    1314        7570 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1315        7570 :       return RgM_sub(x,y);
    1316             :     }
    1317     1886431 :     case t_RFRAC: case t_SER: break;
    1318             : 
    1319           0 :     default: pari_err_TYPE2("+",x,y);
    1320             :   }
    1321    15617345 :   av = avma;
    1322    15617345 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1323             : }
    1324             : 
    1325             : /********************************************************************/
    1326             : /**                                                                **/
    1327             : /**                        MULTIPLICATION                          **/
    1328             : /**                                                                **/
    1329             : /********************************************************************/
    1330             : static GEN
    1331      229321 : mul_ser_scal(GEN y, GEN x) {
    1332             :   long ly, i;
    1333             :   GEN z;
    1334      229321 :   if (isexactzero(x)) return gmul(RgX_get_0(y), x);
    1335      228026 :   if (ser_isexactzero(y))
    1336             :   {
    1337        1393 :     if (lg(y) == 2) return gcopy(y);
    1338           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1339             :   }
    1340      226633 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1341      226633 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1342      226633 :   return normalize(z);
    1343             : }
    1344             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1345             :  * [n/d a valid RFRAC]  */
    1346             : static GEN
    1347     9454599 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1348             : {
    1349     9454599 :   pari_sp av = avma;
    1350             :   GEN z;
    1351             : 
    1352     9454599 :   switch(typ(x))
    1353             :   {
    1354             :     case t_PADIC:
    1355           7 :       n = gmul(n, x);
    1356           7 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1357           7 :       return gerepileupto(av, gdiv(n,d));
    1358             : 
    1359             :     case t_INTMOD: case t_POLMOD:
    1360          28 :       n = gmul(n, x);
    1361          28 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1362          28 :       return gerepileupto(av, gdiv(n,d));
    1363             :   }
    1364     9454564 :   z = gred_rfrac2_i(x, d);
    1365     9454564 :   n = simplify_shallow(n);
    1366     9454564 :   if (typ(z) == t_RFRAC)
    1367             :   {
    1368     7119872 :     n = gmul(gel(z,1), n);
    1369     7119872 :     d = gel(z,2);
    1370     7119872 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1371           0 :       z = RgX_Rg_div(n, d);
    1372             :     else
    1373     7119872 :       z = gred_rfrac_simple(n, d);
    1374             :   }
    1375             :   else
    1376     2334692 :     z = gmul(z, n);
    1377     9454564 :   return gerepileupto(av, z);
    1378             : }
    1379             : static GEN
    1380    57659343 : mul_scal(GEN y, GEN x, long ty)
    1381             : {
    1382    57659343 :   switch(ty)
    1383             :   {
    1384             :     case t_POL:
    1385    49394405 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1386    48580224 :       return RgX_Rg_mul(y, x);
    1387      153259 :     case t_SER: return mul_ser_scal(y, x);
    1388     8111679 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1389             :     case t_QFI: case t_QFR:
    1390           0 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1391             :   }
    1392           0 :   pari_err_TYPE2("*",x,y);
    1393           0 :   return NULL; /* not reached */
    1394             : }
    1395             : 
    1396             : static GEN
    1397      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1398             : {
    1399      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1400      160830 :   long vx = gvar(X), vy = varn(y2);
    1401      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1402        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1403             : }
    1404             : /* (x1/x2) * (y1/y2) */
    1405             : static GEN
    1406     7744009 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1407             : {
    1408             :   GEN z, X, Y;
    1409     7744009 :   pari_sp av = avma;
    1410             : 
    1411     7744009 :   X = gred_rfrac2_i(x1, y2);
    1412     7744009 :   Y = gred_rfrac2_i(y1, x2);
    1413     7744009 :   if (typ(X) == t_RFRAC)
    1414             :   {
    1415     6479589 :     if (typ(Y) == t_RFRAC) {
    1416     6399664 :       x1 = gel(X,1);
    1417     6399664 :       x2 = gel(X,2);
    1418     6399664 :       y1 = gel(Y,1);
    1419     6399664 :       y2 = gel(Y,2);
    1420     6399664 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1421             :     } else
    1422       79925 :       z = mul_gen_rfrac(Y, X);
    1423             :   }
    1424     1264420 :   else if (typ(Y) == t_RFRAC)
    1425       80905 :     z = mul_gen_rfrac(X, Y);
    1426             :   else
    1427     1183515 :     z = gmul(X, Y);
    1428     7744009 :   return gerepileupto(av, z);
    1429             : }
    1430             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1431             : static GEN
    1432      253127 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1433             : {
    1434      253127 :   pari_sp av = avma;
    1435      253127 :   GEN X = gred_rfrac2_i(x1, y2);
    1436      253127 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1437             :   {
    1438      246575 :     x2 = RgX_mul(gel(X,2), x2);
    1439      246575 :     x1 = gel(X,1);
    1440             :   }
    1441             :   else
    1442        6552 :     x1 = X;
    1443      253127 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1444             : }
    1445             : 
    1446             : /* Mod(y, Y) * x,  assuming x scalar */
    1447             : static GEN
    1448     1177885 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1449             : {
    1450     1177885 :   GEN z = cgetg(3,t_POLMOD);
    1451     1177885 :   gel(z,1) = RgX_copy(Y);
    1452     1177885 :   gel(z,2) = gmul(x,y); return z;
    1453             : }
    1454             : 
    1455             : /* cf mulqq */
    1456             : static GEN
    1457     1013349 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1458             : {
    1459     1013349 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1460     1013349 :   pari_sp tetpil, av = avma;
    1461     1013349 :   T[1] = x[1];
    1462     1013349 :   p2 = gmul(gel(x,2), gel(y,2));
    1463     1013349 :   p3 = gmul(gel(x,3), gel(y,3));
    1464     1013349 :   p1 = gmul(gneg_i(c),p3);
    1465             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1466     1013349 :   if (typ(b) == t_INT)
    1467             :   {
    1468     1013349 :     if (signe(b))
    1469             :     {
    1470      728671 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1471      728671 :       if (is_pm1(b))
    1472             :       {
    1473      727929 :         if (signe(b) > 0) p3 = gneg(p3);
    1474             :       }
    1475             :       else
    1476         742 :         p3 = gmul(negi(b), p3);
    1477             :     }
    1478             :     else
    1479             :     {
    1480      284678 :       p3 = gmul(gel(x,2),gel(y,3));
    1481      284678 :       p4 = gmul(gel(x,3),gel(y,2));
    1482             :     }
    1483             :   }
    1484             :   else
    1485             :   {
    1486           0 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1487           0 :     p3 = gmul(gneg_i(b), p3);
    1488             :   }
    1489     1013349 :   tetpil = avma;
    1490     1013349 :   gel(T,2) = gadd(p2, p1);
    1491     1013349 :   gel(T,3) = gadd(p4, p3);
    1492     1013349 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1493     1013349 :   return normalizepol_lg(T,4);
    1494             : }
    1495             : /* Mod(x,T) * Mod(y,T) */
    1496             : static GEN
    1497     1528675 : mul_polmod_same(GEN T, GEN x, GEN y)
    1498             : {
    1499     1528675 :   GEN z = cgetg(3,t_POLMOD), a;
    1500     1528675 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1501     1528675 :   gel(z,1) = RgX_copy(T);
    1502             :   /* x * y mod T optimised */
    1503     1528675 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1504     1313492 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1505      296501 :     a = gmul(x, y);
    1506             :   else
    1507             :   {
    1508     1232174 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1509     1011963 :       a = quad_polmod_mul(T, x, y);
    1510             :     else
    1511             :     {
    1512      220211 :       pari_sp av = avma;
    1513      220211 :       GEN p = NULL;
    1514      220211 :       if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1515             :       {
    1516       37121 :         T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    1517       37121 :         if (lgefint(p) == 3)
    1518             :         {
    1519       32284 :           ulong pp = p[2];
    1520       32284 :           a = Flxq_mul(ZX_to_Flx(x, pp), ZX_to_Flx(y, pp), ZX_to_Flx(T, pp), pp);
    1521       32284 :           a = Flx_to_ZX(a);
    1522             :         }
    1523             :         else
    1524        4837 :           a = FpXQ_mul(x, y, T, p);
    1525       37121 :         a = FpX_to_mod(a, p);
    1526             :       }
    1527             :       else
    1528      183090 :         a = RgXQ_mul(x, y, gel(z,1));
    1529      220211 :       a = gerepileupto(av, a);
    1530             :     }
    1531             :   }
    1532     1528675 :   gel(z,2) = a; return z;
    1533             : }
    1534             : static GEN
    1535       24351 : sqr_polmod(GEN T, GEN x)
    1536             : {
    1537       24351 :   GEN a, z = cgetg(3,t_POLMOD);
    1538       24351 :   gel(z,1) = RgX_copy(T);
    1539       24351 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1540        3341 :     a = gsqr(x);
    1541             :   else
    1542             :   {
    1543       21010 :     pari_sp av = avma;
    1544       21010 :     GEN p = NULL;
    1545       21010 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && p)
    1546             :     {
    1547          14 :       T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p);
    1548          14 :       if (lgefint(p) == 3)
    1549             :       {
    1550           7 :         ulong pp = p[2];
    1551           7 :         a = Flxq_sqr(ZX_to_Flx(x, pp), ZX_to_Flx(T, pp), pp);
    1552           7 :         a = Flx_to_ZX(a);
    1553             :       }
    1554             :       else
    1555           7 :         a = FpXQ_sqr(x, T, p);
    1556          14 :       a = FpX_to_mod(a, p);
    1557             :     }
    1558             :     else
    1559       20996 :       a = RgXQ_sqr(x, gel(z,1));
    1560       21010 :     a = gerepileupto(av, a);
    1561             :   }
    1562       24351 :   gel(z,2) = a; return z;
    1563             : }
    1564             : /* Mod(x,X) * Mod(y,Y) */
    1565             : static GEN
    1566        1386 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1567             : {
    1568        1386 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1569        1386 :   long vx = varn(X), vy = varn(Y);
    1570        1386 :   GEN z = cgetg(3,t_POLMOD);
    1571             : 
    1572        1386 :   if (vx==vy) {
    1573             :     pari_sp av;
    1574          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1575          14 :     warn_coercion(X,Y,gel(z,1));
    1576          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1577          14 :     return z;
    1578             :   }
    1579        1372 :   if (varncmp(vx, vy) < 0)
    1580          21 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1581             :   else
    1582        1351 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1583        1372 :   gel(z,2) = gmul(x, y); return z;
    1584             : }
    1585             : 
    1586             : #if 0 /* used by 3M only */
    1587             : /* set z = x+y and return 1 if x,y have the same sign
    1588             :  * set z = x-y and return 0 otherwise */
    1589             : static int
    1590             : did_add(GEN x, GEN y, GEN *z)
    1591             : {
    1592             :   long tx = typ(x), ty = typ(y);
    1593             :   if (tx == ty) switch(tx)
    1594             :   {
    1595             :     case t_INT: *z = addii(x,y); return 1;
    1596             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1597             :     case t_REAL:
    1598             :       if (signe(x) == -signe(y))
    1599             :       { *z = subrr(x,y); return 0; }
    1600             :       else
    1601             :       { *z = addrr(x,y); return 1; }
    1602             :   }
    1603             :   if (tx == t_REAL) switch(ty)
    1604             :   {
    1605             :     case t_INT:
    1606             :       if (signe(x) == -signe(y))
    1607             :       { *z = subri(x,y); return 0; }
    1608             :       else
    1609             :       { *z = addri(x,y); return 1; }
    1610             :     case t_FRAC:
    1611             :       if (signe(x) == -signe(gel(y,1)))
    1612             :       { *z = gsub(x,y); return 0; }
    1613             :       else
    1614             :       { *z = gadd(x,y); return 1; }
    1615             :   }
    1616             :   else if (ty == t_REAL) switch(tx)
    1617             :   {
    1618             :     case t_INT:
    1619             :       if (signe(x) == -signe(y))
    1620             :       { *z = subir(x,y); return 0; }
    1621             :       else
    1622             :       { *z = addir(x,y); return 1; }
    1623             :     case t_FRAC:
    1624             :       if (signe(gel(x,1)) == -signe(y))
    1625             :       { *z = gsub(x,y); return 0; }
    1626             :       else
    1627             :       { *z = gadd(x,y); return 1; }
    1628             :   }
    1629             :   *z = gadd(x,y); return 1;
    1630             : }
    1631             : #endif
    1632             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1633             : static GEN
    1634      897932 : mulcIR(GEN x, GEN y)
    1635             : {
    1636      897932 :   GEN z = cgetg(3,t_COMPLEX);
    1637      897932 :   pari_sp av = avma;
    1638      897932 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1639      897932 :   gel(z,2) = gmul(y, gel(x,1));
    1640      897932 :   return z;
    1641             : 
    1642             : }
    1643             : /* x,y COMPLEX */
    1644             : static GEN
    1645    70587275 : mulcc(GEN x, GEN y)
    1646             : {
    1647    70587275 :   GEN xr = gel(x,1), xi = gel(x,2);
    1648    70587275 :   GEN yr = gel(y,1), yi = gel(y,2);
    1649             :   GEN p1, p2, p3, p4, z;
    1650             :   pari_sp tetpil, av;
    1651             : 
    1652    70587275 :   if (isintzero(xr))
    1653             :   {
    1654      902127 :     if (isintzero(yr)) {
    1655      130779 :       av = avma;
    1656      130779 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1657             :     }
    1658      771348 :     return mulcIR(y, xi);
    1659             :   }
    1660    69685148 :   if (isintzero(yr)) return mulcIR(x, yi);
    1661             : 
    1662    69558564 :   z = cgetg(3,t_COMPLEX); av = avma;
    1663             : #if 0
    1664             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1665             :    * e.g. xr + xi if exponents differ */
    1666             :   if (did_add(xr, xi, &p3))
    1667             :   {
    1668             :     if (did_add(yr, yi, &p4)) {
    1669             :     /* R = xr*yr - xi*yi
    1670             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1671             :       p1 = gmul(xr,yr);
    1672             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1673             :       p3 = gmul(p3, p4);
    1674             :       p4 = gsub(p2, p1);
    1675             :     } else {
    1676             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1677             :      * I = xr*yi + xi*yr */
    1678             :       p1 = gmul(p3,p4);
    1679             :       p3 = gmul(xr,yi);
    1680             :       p4 = gmul(xi,yr);
    1681             :       p2 = gsub(p3, p4);
    1682             :     }
    1683             :   } else {
    1684             :     if (did_add(yr, yi, &p4)) {
    1685             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1686             :       * I = xr*yi +xi*yr */
    1687             :       p1 = gmul(p3,p4);
    1688             :       p3 = gmul(xr,yi);
    1689             :       p4 = gmul(xi,yr);
    1690             :       p2 = gsub(p4, p3);
    1691             :     } else {
    1692             :     /* R = xr*yr - xi*yi
    1693             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1694             :       p3 = gneg( gmul(p3, p4) );
    1695             :       p1 = gmul(xr,yr);
    1696             :       p2 = gmul(xi,yi);
    1697             :       p4 = gadd(p1, p2);
    1698             : 
    1699             :       p2 = gneg(p2);
    1700             :     }
    1701             :   }
    1702             :   tetpil = avma;
    1703             :   gel(z,1) = gadd(p1,p2);
    1704             :   gel(z,2) = gadd(p3,p4);
    1705             : #else
    1706    69558564 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1707             :   { /* 3M formula */
    1708         168 :     p3 = addii(xr,xi);
    1709         168 :     p4 = addii(yr,yi);
    1710         168 :     p1 = mulii(xr,yr);
    1711         168 :     p2 = mulii(xi,yi);
    1712         168 :     p3 = mulii(p3,p4);
    1713         168 :     p4 = addii(p2,p1);
    1714         168 :     tetpil = avma;
    1715         168 :     gel(z,1) = subii(p1,p2);
    1716         168 :     gel(z,2) = subii(p3,p4);
    1717         287 :     if (!signe(gel(z,2)))
    1718          49 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1719             :   }
    1720             :   else
    1721             :   { /* naive 4M formula: avoid all loss of accuracy */
    1722    69558396 :     p1 = gmul(xr,yr);
    1723    69558396 :     p2 = gmul(xi,yi);
    1724    69558396 :     p3 = gmul(xr,yi);
    1725    69558396 :     p4 = gmul(xi,yr);
    1726    69558396 :     tetpil = avma;
    1727    69558396 :     gel(z,1) = gsub(p1,p2);
    1728    69558396 :     gel(z,2) = gadd(p3,p4);
    1729    69558396 :     if (isintzero(gel(z,2)))
    1730             :     {
    1731        9132 :       cgiv(gel(z,2));
    1732        9132 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1733             :     }
    1734             :   }
    1735             : #endif
    1736             : 
    1737    69549383 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1738             : }
    1739             : /* x,y PADIC */
    1740             : static GEN
    1741     1236099 : mulpp(GEN x, GEN y) {
    1742     1236099 :   long l = valp(x) + valp(y);
    1743             :   pari_sp av;
    1744             :   GEN z, t;
    1745     1236099 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1746     1236099 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1747     1076688 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1748             : 
    1749     1016348 :   t = (precp(x) > precp(y))? y: x;
    1750     1016348 :   z = cgetp(t); setvalp(z,l); av = avma;
    1751     1016350 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1752     1016353 :   avma = av; return z;
    1753             : }
    1754             : /* x,y QUAD */
    1755             : static GEN
    1756        1344 : mulqq(GEN x, GEN y) {
    1757        1344 :   GEN z = cgetg(4,t_QUAD);
    1758        1344 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1759             :   pari_sp av, tetpil;
    1760        1344 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1761             : 
    1762        1344 :   gel(z,1) = ZX_copy(P); av = avma;
    1763        1344 :   p2 = gmul(gel(x,2),gel(y,2));
    1764        1344 :   p3 = gmul(gel(x,3),gel(y,3));
    1765        1344 :   p1 = gmul(gneg_i(c),p3);
    1766             : 
    1767        1344 :   if (signe(b))
    1768        1239 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1769             :   else
    1770             :   {
    1771         105 :     p3 = gmul(gel(x,2),gel(y,3));
    1772         105 :     p4 = gmul(gel(x,3),gel(y,2));
    1773             :   }
    1774        1344 :   tetpil = avma;
    1775        1344 :   gel(z,2) = gadd(p2,p1);
    1776        1344 :   gel(z,3) = gadd(p4,p3);
    1777        1344 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1778             : }
    1779             : 
    1780             : GEN
    1781     2119481 : mulcxI(GEN x)
    1782             : {
    1783             :   GEN z;
    1784     2119481 :   switch(typ(x))
    1785             :   {
    1786             :     case t_INT: case t_REAL: case t_FRAC:
    1787      213490 :       return mkcomplex(gen_0, x);
    1788             :     case t_COMPLEX:
    1789     1905935 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1790     1903814 :       z = cgetg(3,t_COMPLEX);
    1791     1903814 :       gel(z,1) = gneg(gel(x,2));
    1792     1903814 :       gel(z,2) = gel(x,1); return z;
    1793             :     default:
    1794          56 :       return gmul(gen_I(), x);
    1795             :   }
    1796             : }
    1797             : GEN
    1798       47845 : mulcxmI(GEN x)
    1799             : {
    1800             :   GEN z;
    1801       47845 :   switch(typ(x))
    1802             :   {
    1803             :     case t_INT: case t_REAL: case t_FRAC:
    1804        2051 :       return mkcomplex(gen_0, gneg(x));
    1805             :     case t_COMPLEX:
    1806       45738 :       if (isintzero(gel(x,1))) return gel(x,2);
    1807       44744 :       z = cgetg(3,t_COMPLEX);
    1808       44744 :       gel(z,1) = gel(x,2);
    1809       44744 :       gel(z,2) = gneg(gel(x,1)); return z;
    1810             :     default:
    1811          56 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1812             :   }
    1813             : }
    1814             : 
    1815             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1816             : static GEN
    1817     1645717 : fill_ser(GEN z, GEN y)
    1818             : {
    1819     1645717 :   long i, lx = lg(z), ly = lg(y);
    1820     1645717 :   if (ly >= lx) {
    1821     1640726 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1822             :   } else {
    1823        4991 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1824        4991 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1825             :   }
    1826     1645717 :   return normalize(z);
    1827             : }
    1828             : 
    1829             : GEN
    1830  3297792875 : gmul(GEN x, GEN y)
    1831             : {
    1832             :   long tx, ty, lx, ly, vx, vy, i, l;
    1833             :   pari_sp av, tetpil;
    1834             :   GEN z, p1, p2;
    1835             : 
    1836  3297792875 :   if (x == y) return gsqr(x);
    1837  2456308571 :   tx = typ(x); ty = typ(y);
    1838  2456308571 :   if (tx == ty) switch(tx)
    1839             :   {
    1840  1435531586 :     case t_INT: return mulii(x,y);
    1841   510916371 :     case t_REAL: return mulrr(x,y);
    1842     1424824 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1843     1424824 :       z = cgetg(3,t_INTMOD);
    1844     1424816 :       if (X==Y || equalii(X,Y))
    1845     1424801 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1846          14 :       gel(z,1) = gcdii(X,Y);
    1847          14 :       warn_coercion(X,Y,gel(z,1));
    1848          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1849          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1850             :     }
    1851             :     case t_FRAC:
    1852             :     {
    1853    11583196 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1854    11583196 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1855    11583196 :       z=cgetg(3,t_FRAC);
    1856    11583196 :       p1 = gcdii(x1, y2);
    1857    11583196 :       if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1858    11583196 :       p1 = gcdii(x2, y1);
    1859    11583196 :       if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1860    11583196 :       tetpil = avma;
    1861    11583196 :       gel(z,2) = mulii(x2,y2);
    1862    11583196 :       gel(z,1) = mulii(x1,y1);
    1863    11583196 :       fix_frac_if_int_GC(z,tetpil); return z;
    1864             :     }
    1865    69230710 :     case t_COMPLEX: return mulcc(x, y);
    1866     1122656 :     case t_PADIC: return mulpp(x, y);
    1867        1022 :     case t_QUAD: return mulqq(x, y);
    1868    10550899 :     case t_FFELT: return FF_mul(x, y);
    1869             :     case t_POLMOD:
    1870     1334306 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1871     1332920 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1872        1386 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1873             :     case t_POL:
    1874    33499974 :       vx = varn(x);
    1875    33499974 :       vy = varn(y);
    1876    33499974 :       if (vx != vy) {
    1877     4688564 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1878     1998808 :         else                     return RgX_Rg_mul(y, x);
    1879             :       }
    1880    28811410 :       if (RgX_is_ZX(x) && RgX_is_ZX(y)) return ZX_mul(x,y);
    1881     5178474 :       return RgX_mul(x, y);
    1882             : 
    1883             :     case t_SER: {
    1884     1645584 :       GEN p = NULL;
    1885     1645584 :       vx = varn(x);
    1886     1645584 :       vy = varn(y);
    1887     1645584 :       if (vx != vy) {
    1888           0 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1889           0 :         else                     return mul_ser_scal(y, x);
    1890             :       }
    1891     1645584 :       lx = lg(x);
    1892     1645584 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1893     1645584 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1894     1645570 :       z = cgetg(lx,t_SER);
    1895     1645570 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1896     1645570 :       x = ser2pol_i(x, lx);
    1897     1645570 :       y = ser2pol_i(y, lx);
    1898     1645570 :       if (RgX_is_FpX(x,&p) && RgX_is_FpX(y,&p))
    1899             :       {
    1900       29080 :         if (!p) y = ZX_mul(x,y);
    1901             :         else
    1902             :         {
    1903          21 :           x = RgX_to_FpX(x, p);
    1904          21 :           y = RgX_to_FpX(y, p);
    1905          21 :           y = FpX_to_mod(ZX_mul(x,y), p);
    1906             :         }
    1907             :       }
    1908             :       else
    1909     1631030 :         y = RgXn_mul(x, y, lx-2);
    1910     1645570 :       z = fill_ser(z, y);
    1911     1645570 :       return gerepilecopy((pari_sp)(z + lx), z);
    1912             :     }
    1913    13211787 :     case t_QFI: return qficomp(x,y);
    1914           7 :     case t_QFR: return qfrcomp(x,y);
    1915     6557680 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1916     1523000 :     case t_MAT: return RgM_mul(x, y);
    1917             : 
    1918             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1919      131621 :       z = cgetg_copy(x, &l);
    1920      131621 :       if (l != lg(y)) break;
    1921     2310854 :       for (i=1; i<l; i++)
    1922             :       {
    1923     2179233 :         long yi = y[i];
    1924     2179233 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1925     2179233 :         z[i] = x[yi];
    1926             :       }
    1927      131621 :       return z;
    1928             : 
    1929             : 
    1930             :     default:
    1931           0 :       pari_err_TYPE2("*",x,y);
    1932             :   }
    1933             :   /* tx != ty */
    1934   358043348 :   if (is_const_t(ty) && is_const_t(tx))  {
    1935   288916610 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1936   288916610 :     switch(tx) {
    1937             :     case t_INT:
    1938   261931742 :       switch(ty)
    1939             :       {
    1940   159626691 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1941             :         case t_INTMOD:
    1942     1584032 :           z = cgetg(3, t_INTMOD);
    1943     1584027 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1944             :         case t_FRAC:
    1945    55558260 :           if (!signe(x)) return gen_0;
    1946    24591343 :           z=cgetg(3,t_FRAC);
    1947    24591343 :           p1 = gcdii(x,gel(y,2));
    1948    24591343 :           if (is_pm1(p1))
    1949             :           {
    1950    14051417 :             avma = (pari_sp)z;
    1951    14051417 :             gel(z,2) = icopy(gel(y,2));
    1952    14051417 :             gel(z,1) = mulii(gel(y,1), x);
    1953             :           }
    1954             :           else
    1955             :           {
    1956    10539926 :             x = diviiexact(x,p1); tetpil = avma;
    1957    10539926 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1958    10539926 :             gel(z,1) = mulii(gel(y,1), x);
    1959    10539926 :             fix_frac_if_int_GC(z,tetpil);
    1960             :           }
    1961    24591343 :           return z;
    1962    43723081 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1963      104755 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1964        1281 :         case t_QUAD: return mulRq(x,y);
    1965     1333641 :         case t_FFELT: return FF_Z_mul(y,x);
    1966             :       }
    1967             : 
    1968             :     case t_REAL:
    1969    26894870 :       switch(ty)
    1970             :       {
    1971     1554570 :         case t_FRAC: return mulrfrac(x, y);
    1972    25340279 :         case t_COMPLEX: return mulRc(x, y);
    1973           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1974          14 :         default: pari_err_TYPE2("*",x,y);
    1975             :       }
    1976             : 
    1977             :     case t_INTMOD:
    1978        5359 :       switch(ty)
    1979             :       {
    1980         504 :         case t_FRAC: { GEN X = gel(x,1);
    1981         504 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1982         504 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1983             :         }
    1984          35 :         case t_COMPLEX: return mulRc_direct(x,y);
    1985           0 :         case t_PADIC: { GEN X = gel(x,1);
    1986           0 :           z = cgetg(3, t_INTMOD);
    1987           0 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1988             :         }
    1989          42 :         case t_QUAD: return mulRq(x, y);
    1990             :         case t_FFELT:
    1991        4778 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1992           0 :             pari_err_OP("*",x,y);
    1993        4778 :           return FF_Z_mul(y,gel(x,2));
    1994             :       }
    1995             : 
    1996             :     case t_FRAC:
    1997       84575 :       switch(ty)
    1998             :       {
    1999       69133 :         case t_COMPLEX: return mulRc(x, y);
    2000       13146 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    2001         455 :         case t_QUAD:  return mulRq(x, y);
    2002        1841 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    2003             :       }
    2004             : 
    2005             :     case t_FFELT:
    2006          35 :       pari_err_TYPE2("*",x,y);
    2007             : 
    2008             :     case t_COMPLEX:
    2009          14 :       switch(ty)
    2010             :       {
    2011             :         case t_PADIC:
    2012           7 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    2013             :         case t_QUAD:
    2014           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    2015           7 :           return mulqf(y, x, lx);
    2016             :       }
    2017             : 
    2018             :     case t_PADIC: /* ty == t_QUAD */
    2019          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    2020             :     }
    2021             :   }
    2022             : 
    2023    69194384 :   if (is_matvec_t(ty))
    2024             :   {
    2025     8754945 :     if (!is_matvec_t(tx))
    2026             :     {
    2027     8493182 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    2028     8493182 :       z = cgetg_copy(y, &ly);
    2029     8493182 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2030     8493182 :       return z;
    2031             :     }
    2032      261763 :     switch(tx)
    2033             :     {
    2034             :       case t_VEC:
    2035      170808 :         switch(ty) {
    2036      113355 :           case t_COL: return RgV_RgC_mul(x,y);
    2037       57453 :           case t_MAT: return RgV_RgM_mul(x,y);
    2038             :         }
    2039           0 :         break;
    2040             :       case t_COL:
    2041           7 :         switch(ty) {
    2042           7 :           case t_VEC: return RgC_RgV_mul(x,y);
    2043           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    2044             :         }
    2045           0 :         break;
    2046             :       case t_MAT:
    2047       90948 :         switch(ty) {
    2048           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    2049       90948 :           case t_COL: return RgM_RgC_mul(x,y);
    2050             :         }
    2051             :     }
    2052             :   }
    2053    60439444 :   if (is_matvec_t(tx))
    2054             :   {
    2055      140120 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2056      140120 :     z = cgetg_copy(x, &lx);
    2057      140118 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2058      140121 :     return z;
    2059             :   }
    2060    60299325 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2061             :   /* tx < ty, !ismatvec(x and y) */
    2062             : 
    2063    60299325 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2064     1176989 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2065    59122336 :   if (is_scalar_t(tx)) {
    2066    56286739 :     if (tx == t_POLMOD) {
    2067      230294 :       vx = varn(gel(x,1));
    2068      230294 :       vy = gvar(y);
    2069      230294 :       if (vx != vy) {
    2070       34672 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2071         896 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2072             :       }
    2073             :       /* error if ty == t_SER */
    2074      195622 :       av = avma; y = gmod(y, gel(x,1));
    2075      195615 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2076             :     }
    2077    56056445 :     return mul_scal(y, x, ty);
    2078             :   }
    2079             : 
    2080             :   /* x and y are not scalars, nor matvec */
    2081     2835597 :   vx = gvar(x);
    2082     2835597 :   vy = gvar(y);
    2083     2835597 :   if (vx != vy) /* x or y is treated as a scalar */
    2084     4184661 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2085     2770629 :                                 : mul_scal(y, x, ty);
    2086             :   /* vx = vy */
    2087     1421565 :   switch(tx)
    2088             :   {
    2089             :     case t_POL:
    2090     1421565 :       switch (ty)
    2091             :       {
    2092             :         case t_SER:
    2093             :         {
    2094             :           long vn;
    2095       78673 :           if (lg(x) == 2) return zeropol(vx);
    2096       78673 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2097       78463 :           av = avma;
    2098       78463 :           vn = RgX_valrem(x, &x);
    2099             :           /* take advantage of x = t^n ! */
    2100       78463 :           if (degpol(x)) {
    2101        2401 :             p1 = RgX_to_ser(x,lg(y));
    2102        2401 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2103        2401 :             p2 = gmul(p1,y);
    2104        2401 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2105             :           } else {
    2106       76062 :             avma = av;
    2107       76062 :             p2 = mul_ser_scal(y, gel(x,2));
    2108             :           }
    2109       78463 :           setvalp(p2, valp(p2) + vn);
    2110       78463 :           return p2;
    2111             :         }
    2112             : 
    2113     1342892 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2114             :       }
    2115           0 :       break;
    2116             : 
    2117             :     case t_SER:
    2118           0 :       switch (ty)
    2119             :       {
    2120             :         case t_RFRAC:
    2121           0 :           av = avma;
    2122           0 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2123             :       }
    2124           0 :       break;
    2125             :   }
    2126           0 :   pari_err_TYPE2("*",x,y);
    2127           0 :   return NULL; /* not reached */
    2128             : }
    2129             : 
    2130             : int
    2131       29834 : ff_poltype(GEN *x, GEN *p, GEN *pol)
    2132             : {
    2133       29834 :   GEN Q, P = *x;
    2134       29834 :   if (!signe(P)) return 0;
    2135       29687 :   if (!RgX_is_FpXQX(P,pol,p) || !*p || !*pol || typ(*pol)!=t_POL)
    2136       29498 :     return 0;
    2137         189 :   Q = RgX_to_FpX(*pol, *p);
    2138         189 :   P = RgX_to_FpXQX(P, Q, *p);
    2139         189 :   *x = ZXX_to_Kronecker(P, degpol(Q));
    2140         189 :   *pol = Q;
    2141         189 :   return 1;
    2142             : }
    2143             : 
    2144             : /* return a non-normalized result */
    2145             : GEN
    2146        9590 : sqr_ser_part(GEN x, long l1, long l2)
    2147             : {
    2148             :   long i, j, l;
    2149             :   pari_sp av;
    2150             :   GEN Z, z, p1, p2;
    2151             :   long mi;
    2152        9590 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2153        9583 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2154        9583 :   Z = cgetg(l2-l1+3, t_SER);
    2155        9583 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2156        9583 :   z = Z + 2-l1;
    2157        9583 :   x += 2; mi = 0;
    2158       31409 :   for (i=0; i<l1; i++)
    2159             :   {
    2160       21826 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2161             :   }
    2162             : 
    2163       57820 :   for (i=l1; i<=l2; i++)
    2164             :   {
    2165       48237 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2166       48237 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2167      406693 :     for (j=i-mi; j<=minss(l,mi); j++)
    2168      358456 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2169       48237 :     p1 = gshift(p1,1);
    2170       48237 :     if ((i&1) == 0 && p2[i>>1])
    2171       16583 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2172       48237 :     gel(z,i) = gerepileupto(av,p1);
    2173             :   }
    2174        9583 :   return Z;
    2175             : }
    2176             : 
    2177             : GEN
    2178   913160684 : gsqr(GEN x)
    2179             : {
    2180             :   long i, lx;
    2181             :   pari_sp av, tetpil;
    2182             :   GEN z, p1, p2, p3, p4;
    2183             : 
    2184   913160684 :   switch(typ(x))
    2185             :   {
    2186   855104482 :     case t_INT: return sqri(x);
    2187    17991832 :     case t_REAL: return sqrr(x);
    2188      144242 :     case t_INTMOD: { GEN X = gel(x,1);
    2189      144242 :       z = cgetg(3,t_INTMOD);
    2190      144236 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2191      144240 :       gel(z,1) = icopy(X); return z;
    2192             :     }
    2193     1006532 :     case t_FRAC: return sqrfrac(x);
    2194             : 
    2195             :     case t_COMPLEX:
    2196     2945646 :       if (isintzero(gel(x,1))) {
    2197       18593 :         av = avma;
    2198       18593 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2199             :       }
    2200     2927053 :       z = cgetg(3,t_COMPLEX); av = avma;
    2201     2927053 :       p1 = gadd(gel(x,1),gel(x,2));
    2202     2927053 :       p2 = gsub(gel(x,1), gel(x,2));
    2203     2927053 :       p3 = gmul(gel(x,1),gel(x,2));
    2204     2927053 :       tetpil = avma;
    2205     2927053 :       gel(z,1) = gmul(p1,p2);
    2206     2927053 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2207             : 
    2208             :     case t_PADIC:
    2209        8267 :       z = cgetg(5,t_PADIC);
    2210        8267 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2211        8267 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2212        8267 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2213        8267 :       gel(z,2) = icopy(gel(x,2));
    2214        8267 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2215        8267 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2216        8267 :       return z;
    2217             : 
    2218          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2219          28 :       p1 = gel(x,1);
    2220          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2221          28 :       p2 = gsqr(gel(x,2));
    2222          28 :       p3 = gsqr(gel(x,3));
    2223          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2224             : 
    2225          28 :       if (gequal0(gel(p1,3)))
    2226             :       {
    2227           7 :         tetpil = avma;
    2228           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2229           7 :         av = avma;
    2230           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2231           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2232             :       }
    2233             : 
    2234          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2235          21 :       tetpil = avma;
    2236          21 :       gel(z,2) = gadd(p2,p4);
    2237          21 :       gel(z,3) = gadd(p1,p3);
    2238          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2239             : 
    2240             :     case t_POLMOD:
    2241       24351 :       return sqr_polmod(gel(x,1), gel(x,2));
    2242             : 
    2243     2622959 :     case t_FFELT: return FF_sqr(x);
    2244             : 
    2245             :     case t_POL:
    2246             :     {
    2247     1729261 :       GEN a = x, p = NULL, pol = NULL;
    2248     1729261 :       av = avma;
    2249     1729261 :       if (RgX_is_ZX(x)) return ZX_sqr(x);
    2250       29834 :       if (ff_poltype(&x,&p,&pol))
    2251             :       {
    2252         189 :         z = ZX_sqr(x);
    2253         189 :         if (p) z = FpX_to_mod(z,p);
    2254         189 :         if (pol) z = Kronecker_to_mod(z,pol);
    2255         189 :         z = gerepileupto(av, z);
    2256             :       }
    2257       29645 :       else { avma = av; z = RgX_sqr(a); }
    2258       29834 :       return z;
    2259             :     }
    2260             : 
    2261             :     case t_SER:
    2262        6622 :       lx = lg(x);
    2263        6622 :       if (ser_isexactzero(x)) {
    2264           7 :         GEN z = gcopy(x);
    2265           7 :         setvalp(z, 2*valp(x));
    2266           7 :         return z;
    2267             :       }
    2268        6615 :       if (lx < 40)
    2269        6468 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2270             :       else
    2271             :       {
    2272         147 :         pari_sp av = avma;
    2273         147 :         GEN z = cgetg(lx,t_SER), p = NULL;
    2274         147 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2275         147 :         x = ser2pol_i(x,lx);
    2276         147 :         if (RgX_is_FpX(x,&p))
    2277             :         {
    2278          49 :           if (!p) x = ZX_sqr(x);
    2279             :           else
    2280             :           {
    2281           0 :             x = RgX_to_FpX(x, p);
    2282           0 :             x = FpX_to_mod(ZX_sqr(x), p);
    2283             :           }
    2284             :         }
    2285             :         else
    2286          98 :           x = RgXn_sqr(x, lx-2);
    2287         147 :         z = fill_ser(z, x);
    2288         147 :         return gerepilecopy(av, z);
    2289             :       }
    2290             : 
    2291           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2292           7 :       gel(z,1) = gsqr(gel(x,1));
    2293           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2294             : 
    2295         791 :     case t_MAT: return RgM_sqr(x);
    2296          14 :     case t_QFR: return qfrsqr(x);
    2297    31572192 :     case t_QFI: return qfisqr(x);
    2298             :     case t_VECSMALL:
    2299        3458 :       z = cgetg_copy(x, &lx);
    2300       56392 :       for (i=1; i<lx; i++)
    2301             :       {
    2302       52934 :         long xi = x[i];
    2303       52934 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2304       52934 :         z[i] = x[xi];
    2305             :       }
    2306        3458 :       return z;
    2307             :   }
    2308           0 :   pari_err_TYPE2("*",x,x);
    2309           0 :   return NULL; /* not reached */
    2310             : }
    2311             : 
    2312             : /********************************************************************/
    2313             : /**                                                                **/
    2314             : /**                           DIVISION                             **/
    2315             : /**                                                                **/
    2316             : /********************************************************************/
    2317             : static GEN
    2318      270646 : div_rfrac_scal(GEN x, GEN y)
    2319             : {
    2320      270646 :   pari_sp av = avma;
    2321      270646 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2322      270646 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2323             : }
    2324             : static GEN
    2325     1004057 : div_scal_rfrac(GEN x, GEN y)
    2326             : {
    2327     1004057 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2328     1004057 :   pari_sp av = avma;
    2329     1004057 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2330             :   {
    2331           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2332           0 :     y1 = gel(y1,2);
    2333             :   }
    2334     1004057 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2335             : }
    2336             : static GEN
    2337     1186329 : div_rfrac(GEN x, GEN y)
    2338     1186329 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2339             : 
    2340             : static GEN
    2341      451420 : div_ser_scal(GEN x, GEN y) {
    2342             :   long i, lx;
    2343             :   GEN z;
    2344      451420 :   if (ser_isexactzero(x))
    2345             :   {
    2346           7 :     if (lg(x) == 2) return gcopy(x);
    2347           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2348             :   }
    2349      451413 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2350      451413 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2351      451413 :   return normalize(z);
    2352             : }
    2353             : GEN
    2354         133 : ser_normalize(GEN x)
    2355             : {
    2356         133 :   long i, lx = lg(x);
    2357             :   GEN c, z;
    2358         133 :   if (lx == 2) return x;
    2359         133 :   c = gel(x,2); if (gequal1(c)) return x;
    2360         105 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2361         105 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2362         105 :   return z;
    2363             : }
    2364             : 
    2365             : static GEN
    2366     1060650 : div_T_scal(GEN x, GEN y, long tx) {
    2367     1060650 :   switch(tx)
    2368             :   {
    2369      342308 :     case t_POL: return RgX_Rg_div(x, y);
    2370      451420 :     case t_SER: return div_ser_scal(x, y);
    2371      266922 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2372             :   }
    2373           0 :   pari_err_TYPE2("/",x,y);
    2374           0 :   return NULL; /* not reached */
    2375             : }
    2376             : 
    2377             : static GEN
    2378     8183747 : div_scal_pol(GEN x, GEN y) {
    2379     8183747 :   long ly = lg(y);
    2380             :   pari_sp av;
    2381     8183747 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2382     8151881 :   if (isrationalzero(x)) return zeropol(varn(y));
    2383     6099930 :   av = avma;
    2384     6099930 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2385             : }
    2386             : static GEN
    2387        9485 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2388             :   GEN z;
    2389             :   long ly, i;
    2390        9485 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2391        9485 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2392        9485 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2393        9485 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2394        9485 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2395        9485 :   y = gdiv(z,y); pari_free(z); return y;
    2396             : }
    2397             : static GEN
    2398     9147631 : div_scal_T(GEN x, GEN y, long ty) {
    2399     9147631 :   switch(ty)
    2400             :   {
    2401     8134845 :     case t_POL: return div_scal_pol(x, y);
    2402        9485 :     case t_SER: return div_scal_ser(x, y);
    2403     1003301 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2404             :   }
    2405           0 :   pari_err_TYPE2("/",x,y);
    2406           0 :   return NULL; /* not reached */
    2407             : }
    2408             : 
    2409             : /* assume tx = ty = t_SER, same variable vx */
    2410             : static GEN
    2411      732049 : div_ser(GEN x, GEN y, long vx)
    2412             : {
    2413      732049 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2414             :   GEN y_lead, p1, p2, z;
    2415             : 
    2416      732049 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2417      732042 :   if (ser_isexactzero(x))
    2418             :   {
    2419         651 :     if (lx == 2) return zeroser(vx, l);
    2420           7 :     return scalarser(gmul(gel(x,2),RgX_get_0(y)), varn(x), l);
    2421             :   }
    2422      731391 :   y_lead = gel(y,2);
    2423      731391 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2424             :   {
    2425           0 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2426           0 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2427             :     {
    2428           0 :       y_lead = gel(y,2);
    2429           0 :       if (!gequal0(y_lead)) break;
    2430             :     }
    2431           0 :     if (ly <= 2) pari_err_INV("div_ser", y);
    2432             :   }
    2433      731391 :   if (ly < lx) lx = ly;
    2434      731391 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2435     2293736 :   for (i=3; i<lx; i++)
    2436             :   {
    2437     1562345 :     p1 = gel(y,i);
    2438     1562345 :     if (isrationalzero(p1)) p1 = NULL;
    2439     1562345 :     gel(p2,i) = p1;
    2440             :   }
    2441      731391 :   z = cgetg(lx,t_SER);
    2442      731391 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2443      731391 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2444     2293736 :   for (i=3; i<lx; i++)
    2445             :   {
    2446     1562345 :     pari_sp av = avma;
    2447     1562345 :     p1 = gel(x,i);
    2448    10838548 :     for (j=2, l=i; j<i; j++, l--)
    2449     9276203 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2450     1562345 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2451             :   }
    2452      731391 :   return normalize(z);
    2453             : }
    2454             : /* x,y compatible PADIC */
    2455             : static GEN
    2456       92500 : divpp(GEN x, GEN y) {
    2457             :   pari_sp av;
    2458             :   long a, b;
    2459             :   GEN z, M;
    2460             : 
    2461       92500 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2462       92496 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2463       92489 :   a = precp(x);
    2464       92489 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2465       92489 :   z = cgetg(5, t_PADIC);
    2466       92494 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2467       92494 :   gel(z,2) = icopy(gel(x,2));
    2468       92495 :   gel(z,3) = icopy(M); av = avma;
    2469       92499 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2470       92497 :   return z;
    2471             : }
    2472             : static GEN
    2473        6944 : div_polmod_same(GEN T, GEN x, GEN y)
    2474             : {
    2475        6944 :   long v = varn(T);
    2476        6944 :   GEN a, z = cgetg(3, t_POLMOD);
    2477        6944 :   gel(z,1) = RgX_copy(T);
    2478        6944 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2479        4907 :     a = gdiv(x, y);
    2480        2037 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2481         182 :   {
    2482         182 :     pari_sp av = avma;
    2483         182 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2484             :   }
    2485        1855 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2486        1386 :   {
    2487        1386 :     pari_sp av = avma;
    2488        1386 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2489        1386 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2490        1386 :     a = gerepileupto(av, a);
    2491             :   }
    2492             :   else
    2493             :   {
    2494         469 :     pari_sp av = avma;
    2495         469 :     GEN p = NULL;
    2496         469 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    2497             :     {
    2498          14 :       T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    2499          14 :       if (lgefint(p) == 3)
    2500             :       {
    2501           7 :         ulong pp = p[2];
    2502           7 :         x = ZX_to_Flx(x, pp);
    2503           7 :         y = ZX_to_Flx(y, pp);
    2504           7 :         T = ZX_to_Flx(T, pp);
    2505           7 :         a = Flxq_mul(x, Flxq_inv(y,T,pp),T,pp);
    2506           7 :         a = Flx_to_ZX(a);
    2507             :       }
    2508             :       else
    2509           7 :         a = FpXQ_div(x, y, T, p);
    2510          14 :       a = FpX_to_mod(a, p);
    2511             :     } else
    2512         455 :       a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2513         469 :     a = gerepileupto(av, a);
    2514             :   }
    2515        6944 :   gel(z,2) = a; return z;
    2516             : }
    2517             : GEN
    2518   130621003 : gdiv(GEN x, GEN y)
    2519             : {
    2520   130621003 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2521             :   pari_sp av, tetpil;
    2522             :   GEN z, p1, p2;
    2523             : 
    2524   130621003 :   if (tx == ty) switch(tx)
    2525             :   {
    2526             :     case t_INT:
    2527    57657922 :       if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
    2528    22784344 :       if (is_pm1(x)) {
    2529     4754693 :         long s = signe(y);
    2530     4754693 :         if (!s) pari_err_INV("gdiv",y);
    2531     4754616 :         if (signe(x) < 0) s = -s;
    2532     4754616 :         z = cgetg(3, t_FRAC);
    2533     4754616 :         gel(z,1) = s<0? gen_m1: gen_1;
    2534     4754616 :         gel(z,2) = absi(y); return z;
    2535             :       }
    2536    18029651 :       return gred_frac2(x,y);
    2537             : 
    2538    27009340 :     case t_REAL: return divrr(x,y);
    2539       21384 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2540       21384 :       z = cgetg(3,t_INTMOD);
    2541       21383 :       if (X==Y || equalii(X,Y))
    2542       21376 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2543           7 :       gel(z,1) = gcdii(X,Y);
    2544           7 :       warn_coercion(X,Y,gel(z,1));
    2545           7 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2546           7 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2547             :     }
    2548             :     case t_FRAC: {
    2549      866202 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2550      866202 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2551      866202 :       z = cgetg(3, t_FRAC);
    2552      866202 :       p1 = gcdii(x1, y1);
    2553      866202 :       if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2554      866202 :       p1 = gcdii(x2, y2);
    2555      866202 :       if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2556      866202 :       tetpil = avma;
    2557      866202 :       gel(z,2) = mulii(x2,y1);
    2558      866202 :       gel(z,1) = mulii(x1,y2);
    2559      866202 :       normalize_frac(z);
    2560      866202 :       fix_frac_if_int_GC(z,tetpil);
    2561      866202 :       return z;
    2562             :     }
    2563             :     case t_COMPLEX:
    2564     1364769 :       if (isintzero(gel(y,1)))
    2565             :       {
    2566        8204 :         y = gel(y,2);
    2567        8204 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2568        8141 :         z = cgetg(3,t_COMPLEX);
    2569        8141 :         gel(z,1) = gdiv(gel(x,2), y);
    2570        8141 :         av = avma;
    2571        8141 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2572        8141 :         return z;
    2573             :       }
    2574     1356565 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, gconj(y)); tetpil = avma;
    2575     1356565 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2576             : 
    2577             :     case t_PADIC:
    2578        9639 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2579        9639 :       return divpp(x, y);
    2580             : 
    2581             :     case t_QUAD:
    2582         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2583         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, gconj(y)); tetpil = avma;
    2584         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2585             : 
    2586       75096 :     case t_FFELT: return FF_div(x,y);
    2587             : 
    2588             :     case t_POLMOD:
    2589        6951 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2590        6944 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2591             :       else {
    2592           7 :         av = avma;
    2593           7 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2594             :       }
    2595        6951 :       return z;
    2596             : 
    2597             :     case t_POL:
    2598    18381245 :       vx = varn(x);
    2599    18381245 :       vy = varn(y);
    2600    18381245 :       if (vx != vy) {
    2601       97685 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2602       48902 :                             else return div_scal_pol(x, y);
    2603             :       }
    2604    18283560 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2605    18283560 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2606    18174684 :       return gred_rfrac2(x,y);
    2607             : 
    2608             :     case t_SER:
    2609      397091 :       vx = varn(x);
    2610      397091 :       vy = varn(y);
    2611      397091 :       if (vx != vy) {
    2612           0 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2613           0 :                             else return div_scal_ser(x, y);
    2614             :       }
    2615      397091 :       return div_ser(x, y, vx);
    2616             :     case t_RFRAC:
    2617     1190809 :       vx = varn(gel(x,2));
    2618     1190809 :       vy = varn(gel(y,2));
    2619     1190809 :       if (vx != vy) {
    2620        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2621         756 :                             else return div_scal_rfrac(x, y);
    2622             :       }
    2623     1186329 :       return div_rfrac(x,y);
    2624             : 
    2625           0 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2626           0 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2627             : 
    2628             :     case t_MAT:
    2629          14 :       av = avma; p1 = RgM_inv(y);
    2630          14 :       if (!p1) pari_err_INV("gdiv",y);
    2631          14 :       return gerepileupto(av, RgM_mul(x, p1));
    2632             : 
    2633           0 :     default: pari_err_TYPE2("/",x,y);
    2634             :   }
    2635             : 
    2636    23640218 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2637             :   {
    2638      606968 :     long s = signe(x);
    2639      606968 :     if (!s) {
    2640       34502 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2641       34502 :       switch (ty)
    2642             :       {
    2643       31359 :       default: return gen_0;
    2644             :       case t_INTMOD:
    2645          14 :         z = cgetg(3,t_INTMOD);
    2646          14 :         gel(z,1) = icopy(gel(y,1));
    2647          14 :         gel(z,2) = gen_0; return z;
    2648        3129 :       case t_FFELT: return FF_zero(y);
    2649             :       }
    2650             :     }
    2651      572466 :     if (is_pm1(x)) {
    2652      291190 :       if (s > 0) return ginv(y);
    2653       27094 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2654             :     }
    2655      281276 :     switch(ty)
    2656             :     {
    2657       32850 :       case t_REAL: return divir(x,y);
    2658             :       case t_INTMOD:
    2659          42 :         z = cgetg(3, t_INTMOD);
    2660          42 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2661             :       case t_FRAC:
    2662      226156 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2663      226156 :         if (is_pm1(p1))
    2664             :         {
    2665      221684 :           avma = (pari_sp)z;
    2666      221684 :           gel(z,2) = icopy(gel(y,1));
    2667      221684 :           gel(z,1) = mulii(gel(y,2), x);
    2668      221684 :           normalize_frac(z);
    2669      221684 :           fix_frac_if_int(z);
    2670             :         }
    2671             :         else
    2672             :         {
    2673        4472 :           x = diviiexact(x,p1); tetpil = avma;
    2674        4472 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2675        4472 :           gel(z,1) = mulii(gel(y,2), x);
    2676        4472 :           normalize_frac(z);
    2677        4472 :           fix_frac_if_int_GC(z,tetpil);
    2678             :         }
    2679      226156 :         return z;
    2680             : 
    2681         231 :       case t_FFELT: return Z_FF_div(x,y);
    2682       21913 :       case t_COMPLEX: return divRc(x,y);
    2683          42 :       case t_PADIC: return divTp(x, y);
    2684             :       case t_QUAD:
    2685          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, gconj(y)); tetpil = avma;
    2686          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2687             :     }
    2688             :   }
    2689    23033250 :   if (gequal0(y))
    2690             :   {
    2691          35 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2692          14 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2693             :   }
    2694             : 
    2695    23033222 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2696             :   {
    2697             :     case t_REAL:
    2698     9038018 :       switch(ty)
    2699             :       {
    2700     8045090 :         case t_INT: return divri(x,y);
    2701             :         case t_FRAC:
    2702      667848 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2703      667848 :           return gerepileuptoleaf(av, z);
    2704      325052 :         case t_COMPLEX: return divRc(x, y);
    2705           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2706          21 :         default: pari_err_TYPE2("/",x,y);
    2707             :       }
    2708             : 
    2709             :     case t_INTMOD:
    2710         574 :       switch(ty)
    2711             :       {
    2712             :         case t_INT:
    2713         441 :           z = cgetg(3, t_INTMOD);
    2714         441 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2715          28 :         case t_FRAC: { GEN X = gel(x,1);
    2716          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2717          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2718             :         }
    2719             :         case t_FFELT:
    2720          63 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2721           0 :             pari_err_OP("/",x,y);
    2722          63 :           return Z_FF_div(gel(x,2),y);
    2723             : 
    2724             :         case t_COMPLEX:
    2725           0 :           av = avma;
    2726           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), gconj(y)));
    2727             : 
    2728             :         case t_QUAD:
    2729           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,gconj(y)); tetpil = avma;
    2730           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2731             : 
    2732           7 :         case t_PADIC: { GEN X = gel(x,1);
    2733           7 :           z = cgetg(3, t_INTMOD);
    2734           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2735             :         }
    2736          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2737             :       }
    2738             : 
    2739             :     case t_FRAC:
    2740      311730 :       switch(ty)
    2741             :       {
    2742      279142 :         case t_INT: z = cgetg(3, t_FRAC);
    2743      279142 :         p1 = gcdii(y,gel(x,1));
    2744      279142 :         if (is_pm1(p1))
    2745             :         {
    2746      171670 :           avma = (pari_sp)z; tetpil = 0;
    2747      171670 :           gel(z,1) = icopy(gel(x,1));
    2748             :         }
    2749             :         else
    2750             :         {
    2751      107472 :           y = diviiexact(y,p1); tetpil = avma;
    2752      107472 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2753             :         }
    2754      279142 :         gel(z,2) = mulii(gel(x,2),y);
    2755      279142 :         normalize_frac(z);
    2756      279142 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2757      279142 :         return z;
    2758             : 
    2759             :         case t_REAL:
    2760       32273 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2761       32273 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2762             : 
    2763           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2764           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2765           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2766             :         }
    2767             : 
    2768          28 :         case t_FFELT: av=avma;
    2769          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2770             : 
    2771         231 :         case t_COMPLEX: return divRc(x, y);
    2772             : 
    2773             :         case t_PADIC:
    2774          42 :           if (!signe(gel(x,1))) return gen_0;
    2775          42 :           return divTp(x, y);
    2776             : 
    2777             :         case t_QUAD:
    2778           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
    2779           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2780             :       }
    2781             : 
    2782             :     case t_FFELT:
    2783        2527 :       switch (ty)
    2784             :       {
    2785        2408 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2786          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2787             :         case t_INTMOD:
    2788          49 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2789           0 :             pari_err_OP("/",x,y);
    2790          49 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2791             :         default:
    2792          42 :         pari_err_TYPE2("/",x,y);
    2793             :       }
    2794           0 :       break;
    2795             : 
    2796             :     case t_COMPLEX:
    2797     2351205 :       switch(ty)
    2798             :       {
    2799     2351205 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2800           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2801             :         case t_PADIC:
    2802           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2803             :         case t_QUAD:
    2804           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2805           0 :           return divfq(x, y, lx);
    2806             :       }
    2807             : 
    2808             :     case t_PADIC:
    2809        1267 :       switch(ty)
    2810             :       {
    2811        1218 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2812        2436 :           return signe(gel(x,4))? divpT(x, y)
    2813        1218 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2814             :         }
    2815           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2816           7 :           z = cgetg(3, t_INTMOD);
    2817           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2818             :         }
    2819             :         case t_COMPLEX: case t_QUAD:
    2820           7 :           av=avma; p1=gmul(x,gconj(y)); p2=gnorm(y); tetpil=avma;
    2821           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2822             : 
    2823          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2824             :       }
    2825             : 
    2826             :     case t_QUAD:
    2827         889 :       switch (ty)
    2828             :       {
    2829             :         case t_INT: case t_INTMOD: case t_FRAC:
    2830         840 :           z = cgetg(4,t_QUAD);
    2831         840 :           gel(z,1) = ZX_copy(gel(x,1));
    2832         840 :           gel(z,2) = gdiv(gel(x,2), y);
    2833         840 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2834          28 :         case t_REAL: return divqf(x, y, lg(y));
    2835           7 :         case t_PADIC: return divTp(x, y);
    2836             :         case t_COMPLEX:
    2837           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2838           0 :           return divqf(x, y, ly);
    2839             :       }
    2840             :   }
    2841    11327033 :   switch(ty) {
    2842             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2843       22441 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2844             :     case t_MAT:
    2845          28 :       av = avma; p1 = RgM_inv(y);
    2846          28 :       if (!p1) pari_err_INV("gdiv",y);
    2847          21 :       return gerepileupto(av, gmul(x, p1));
    2848             :     case t_VEC: case t_COL:
    2849             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2850           0 :       pari_err_TYPE2("/",x,y);
    2851             :   }
    2852    11304564 :   switch(tx) {
    2853             :     case t_VEC: case t_COL: case t_MAT:
    2854      417339 :       z = cgetg_copy(x, &lx);
    2855      417339 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2856      417339 :       return z;
    2857             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2858           0 :       pari_err_TYPE2("/",x,y);
    2859             :   }
    2860             : 
    2861    10887225 :   vy = gvar(y);
    2862    10887225 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2863       11445 :     vx = varn(X);
    2864       11445 :     if (vx != vy) {
    2865       11298 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2866       11067 :       z = cgetg(3,t_POLMOD);
    2867       11067 :       gel(z,1) = RgX_copy(X);
    2868       11067 :       gel(z,2) = gdiv(gel(x,2), y); return z;
    2869             :     }
    2870             :     /* y is POL, SER or RFRAC */
    2871         147 :     av = avma;
    2872         147 :     switch(ty)
    2873             :     {
    2874           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2875         147 :       default: y = ginvmod(gmod(y,X), X);
    2876             :     }
    2877         140 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2878             :   }
    2879             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2880             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2881             :    * variable NO_VARIABLE, then the operation is incorrect */
    2882    10875780 :   vx = gvar(x);
    2883    10875780 :   if (vx != vy) { /* includes cases where one is scalar */
    2884    10208050 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2885     9147400 :                         else return div_scal_T(x, y, ty);
    2886             :   }
    2887      667730 :   switch(tx)
    2888             :   {
    2889             :     case t_POL:
    2890       79554 :       switch(ty)
    2891             :       {
    2892             :         case t_SER:
    2893          28 :           if (lg(y) == 2)
    2894           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2895          28 :           p1 = RgX_to_ser(x,lg(y));
    2896          28 :           p2 = div_ser(p1, y, vx);
    2897          28 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2898          28 :           return p2;
    2899             : 
    2900             :         case t_RFRAC:
    2901             :         {
    2902       79526 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2903       79526 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2904           0 :             return mul_rfrac_scal(y2, y1, x);
    2905       79526 :           av = avma;
    2906       79526 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2907             :         }
    2908             :       }
    2909           0 :       break;
    2910             : 
    2911             :     case t_SER:
    2912      335035 :       switch(ty)
    2913             :       {
    2914             :         case t_POL:
    2915      335035 :           if (lg(x) == 2)
    2916         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2917      334930 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2918      334930 :           p2 = div_ser(x, p1, vx);
    2919      334930 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2920      334930 :           return p2;
    2921             :         case t_RFRAC:
    2922           0 :           av = avma;
    2923           0 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2924             :       }
    2925           0 :       break;
    2926             : 
    2927             :     case t_RFRAC:
    2928      253127 :       switch(ty)
    2929             :       {
    2930      253127 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2931             :         case t_SER:
    2932           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2933           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2934             :       }
    2935           0 :       break;
    2936             :   }
    2937          14 :   pari_err_TYPE2("/",x,y);
    2938           0 :   return NULL; /* not reached */
    2939             : }
    2940             : 
    2941             : /********************************************************************/
    2942             : /**                                                                **/
    2943             : /**                     SIMPLE MULTIPLICATION                      **/
    2944             : /**                                                                **/
    2945             : /********************************************************************/
    2946             : GEN
    2947    23374565 : gmulsg(long s, GEN y)
    2948             : {
    2949             :   long ly, i;
    2950             :   pari_sp av;
    2951             :   GEN z;
    2952             : 
    2953    23374565 :   switch(typ(y))
    2954             :   {
    2955    11824686 :     case t_INT:  return mulsi(s,y);
    2956     8496683 :     case t_REAL: return mulsr(s,y);
    2957      168505 :     case t_INTMOD: { GEN p = gel(y,1);
    2958      168505 :       z = cgetg(3,t_INTMOD);
    2959      168498 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2960      168509 :       gel(z,1) = icopy(p); return z;
    2961             :     }
    2962      391733 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2963             :     case t_FRAC:
    2964      911155 :       if (!s) return gen_0;
    2965      910483 :       z = cgetg(3,t_FRAC);
    2966      910483 :       i = cgcd(s, smodis(gel(y,2), s));
    2967      910483 :       if (i == 1)
    2968             :       {
    2969      332591 :         gel(z,2) = icopy(gel(y,2));
    2970      332591 :         gel(z,1) = mulis(gel(y,1), s);
    2971             :       }
    2972             :       else
    2973             :       {
    2974      577892 :         gel(z,2) = divis(gel(y,2), i);
    2975      577892 :         gel(z,1) = mulis(gel(y,1), s/i);
    2976      577892 :         fix_frac_if_int(z);
    2977             :       }
    2978      910483 :       return z;
    2979             : 
    2980     1085200 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2981     1085200 :       gel(z,1) = gmulsg(s,gel(y,1));
    2982     1085200 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2983             : 
    2984             :     case t_PADIC:
    2985         882 :       if (!s) return gen_0;
    2986         882 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2987             : 
    2988           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2989           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2990           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2991           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2992             : 
    2993        5908 :     case t_POLMOD: z = cgetg(3, t_POLMOD);
    2994        5908 :       gel(z,1) = RgX_copy(gel(y,1));
    2995        5908 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2996             : 
    2997             :     case t_POL:
    2998      322753 :       if (!signe(y)) return RgX_copy(y);
    2999      319085 :       if (!s) return scalarpol(RgX_get_0(y), varn(y));
    3000      315319 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    3001      315319 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    3002      315319 :       return normalizepol_lg(z, ly);
    3003             : 
    3004             :     case t_SER:
    3005         777 :       if (ser_isexactzero(y)) return gcopy(y);
    3006         777 :       if (!s) return RgX_get_0(y);
    3007         777 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    3008         777 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    3009         777 :       return normalize(z);
    3010             : 
    3011             :     case t_RFRAC:
    3012          14 :       if (!s) return zeropol(varn(gel(y,2)));
    3013          14 :       if (s == 1) return gcopy(y);
    3014          14 :       if (s == -1) return gneg(y);
    3015          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    3016             : 
    3017             :     case t_VEC: case t_COL: case t_MAT:
    3018      166262 :       z = cgetg_copy(y, &ly);
    3019      166262 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    3020      166262 :       return z;
    3021             :   }
    3022           0 :   pari_err_TYPE("gmulsg",y);
    3023           0 :   return NULL; /* not reached */
    3024             : }
    3025             : 
    3026             : /********************************************************************/
    3027             : /**                                                                **/
    3028             : /**                       SIMPLE DIVISION                          **/
    3029             : /**                                                                **/
    3030             : /********************************************************************/
    3031             : 
    3032             : GEN
    3033     6228340 : gdivgs(GEN x, long s)
    3034             : {
    3035     6228340 :   long tx = typ(x), lx, i;
    3036             :   pari_sp av;
    3037             :   GEN z, y, p1;
    3038             : 
    3039     6228340 :   if (!s)
    3040             :   {
    3041           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3042           0 :     pari_err_INV("gdivgs",gen_0);
    3043             :   }
    3044     6228340 :   switch(tx)
    3045             :   {
    3046             :     case t_INT:
    3047     2508681 :       av = avma; z = divis_rem(x,s,&i);
    3048     2508681 :       if (!i) return z;
    3049             : 
    3050     1374268 :       i = cgcd(s, i);
    3051     1374268 :       avma=av; z = cgetg(3,t_FRAC);
    3052     1374268 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    3053     1374268 :       gel(z,1) = y;
    3054     1374268 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    3055             : 
    3056             :     case t_REAL:
    3057     2501859 :       return divrs(x,s);
    3058             : 
    3059             :     case t_INTMOD:
    3060        7350 :       z = cgetg(3, t_INTMOD);
    3061        7350 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    3062             : 
    3063         399 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    3064             : 
    3065      215092 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3066      215092 :       i = cgcd(s, smodis(gel(x,1), s));
    3067      215092 :       if (i == 1)
    3068             :       {
    3069       39255 :         gel(z,2) = mulsi(s, gel(x,2));
    3070       39255 :         gel(z,1) = icopy(gel(x,1));
    3071             :       }
    3072             :       else
    3073             :       {
    3074      175837 :         gel(z,2) = mulsi(s/i, gel(x,2));
    3075      175837 :         gel(z,1) = divis(gel(x,1), i);
    3076             :       }
    3077      215092 :       normalize_frac(z);
    3078      215092 :       fix_frac_if_int(z); return z;
    3079             : 
    3080      893544 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3081      893544 :       gel(z,1) = gdivgs(gel(x,1),s);
    3082      893544 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3083             : 
    3084             :     case t_PADIC: /* divpT */
    3085             :     {
    3086       82375 :       GEN p = gel(x,2);
    3087       82375 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3088       81654 :       av = avma;
    3089       81654 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    3090             :     }
    3091             : 
    3092           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3093           0 :       gel(z,1) = ZX_copy(gel(x,1));
    3094           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    3095           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    3096             : 
    3097        1505 :     case t_POLMOD: z = cgetg(3, t_POLMOD);
    3098        1505 :       gel(z,1) = RgX_copy(gel(x,1));
    3099        1505 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3100             : 
    3101             :     case t_RFRAC:
    3102          91 :       av = avma;
    3103          91 :       p1 = ggcd(stoi(s),gel(x,1));
    3104          91 :       if (typ(p1) == t_INT)
    3105             :       {
    3106          42 :         avma = av;
    3107          42 :         z = cgetg(3, t_RFRAC);
    3108          42 :         i = p1[2];
    3109          42 :         if (i == 1)
    3110             :         {
    3111          42 :           gel(z,1) = gcopy(gel(x,1));
    3112          42 :           gel(z,2) = gmulsg(s,gel(x,2));
    3113             :         }
    3114             :         else
    3115             :         {
    3116           0 :           gel(z,1) = gdivgs(gel(x,1), i);
    3117           0 :           gel(z,2) = gmulgs(gel(x,2), s/i);
    3118             :         }
    3119             :       }
    3120             :       else /* t_FRAC */
    3121             :       {
    3122          49 :         z = cgetg(3, t_RFRAC);
    3123          49 :         gel(z,1) = gdiv(gel(x,1), p1);
    3124          49 :         gel(z,2) = RgX_Rg_mul(gel(x,2), gdivsg(s,p1));
    3125          49 :         z = gerepilecopy(av, z);
    3126             :       }
    3127          91 :       return z;
    3128             : 
    3129             :     case t_POL: case t_SER:
    3130       17437 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3131       17437 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3132       17437 :       return z;
    3133             :     case t_VEC: case t_COL: case t_MAT:
    3134           7 :       z = cgetg_copy(x, &lx);
    3135           7 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3136           7 :       return z;
    3137             : 
    3138             :   }
    3139           0 :   pari_err_TYPE2("/",x, stoi(s));
    3140           0 :   return NULL; /* not reached */
    3141             : }
    3142             : 
    3143             : /* True shift (exact multiplication by 2^n) */
    3144             : GEN
    3145    35589817 : gmul2n(GEN x, long n)
    3146             : {
    3147             :   long lx, i, k, l;
    3148             :   GEN z, a, b;
    3149             : 
    3150    35589817 :   switch(typ(x))
    3151             :   {
    3152             :     case t_INT:
    3153    27000746 :       if (n>=0) return shifti(x,n);
    3154      161406 :       if (!signe(x)) return gen_0;
    3155       65506 :       l = vali(x); n = -n;
    3156       65506 :       if (n<=l) return shifti(x,-n);
    3157       39417 :       z = cgetg(3,t_FRAC);
    3158       39417 :       gel(z,1) = shifti(x,-l);
    3159       39417 :       gel(z,2) = int2n(n-l); return z;
    3160             : 
    3161             :     case t_REAL:
    3162     4038202 :       return shiftr(x,n);
    3163             : 
    3164       78850 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3165       78850 :       z = cgetg(3,t_INTMOD);
    3166       78851 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3167       78844 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3168       78845 :       gel(z,1) = icopy(b); return z;
    3169             : 
    3170      153264 :     case t_FFELT: return FF_mul2n(x,n);
    3171             : 
    3172      401961 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3173      401961 :       l = vali(a);
    3174      401961 :       k = vali(b);
    3175      401961 :       if (n+l >= k)
    3176             :       {
    3177       59267 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3178       51133 :         l = n-k; k = -k;
    3179             :       }
    3180             :       else
    3181             :       {
    3182      342694 :         k = -(l+n); l = -l;
    3183             :       }
    3184      393827 :       z = cgetg(3,t_FRAC);
    3185      393827 :       gel(z,1) = shifti(a,l);
    3186      393827 :       gel(z,2) = shifti(b,k); return z;
    3187             : 
    3188     2755506 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3189     2755506 :       gel(z,1) = gmul2n(gel(x,1),n);
    3190     2755506 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3191             : 
    3192          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3193          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3194          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3195          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3196             : 
    3197       13896 :     case t_POLMOD: z = cgetg(3,t_POLMOD);
    3198       13896 :       gel(z,1) = RgX_copy(gel(x,1));
    3199       13896 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3200             : 
    3201             :     case t_POL:
    3202      266912 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3203      266912 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3204      266912 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3205             :     case t_SER:
    3206       57155 :       if (ser_isexactzero(x)) return gcopy(x);
    3207       57141 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3208       57141 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3209       57141 :       return normalize(z); /* needed if char = 2 */
    3210             :     case t_VEC: case t_COL: case t_MAT:
    3211      812727 :       z = cgetg_copy(x, &lx);
    3212      812727 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3213      812727 :       return z;
    3214             : 
    3215             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3216          14 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3217             : 
    3218             :     case t_PADIC: /* int2n wrong if n < 0 */
    3219       10570 :       return gmul(gmul2n(gen_1,n),x);
    3220             :   }
    3221           0 :   pari_err_TYPE("gmul2n",x);
    3222           0 :   return NULL; /* not reached */
    3223             : }
    3224             : 
    3225             : /*******************************************************************/
    3226             : /*                                                                 */
    3227             : /*                              INVERSE                            */
    3228             : /*                                                                 */
    3229             : /*******************************************************************/
    3230             : static GEN
    3231       20990 : inv_polmod(GEN T, GEN x)
    3232             : {
    3233       20990 :   GEN z = cgetg(3,t_POLMOD), a;
    3234       20990 :   gel(z,1) = RgX_copy(T);
    3235       20990 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3236        7140 :     a = ginv(x);
    3237             :   else
    3238             :   {
    3239       13850 :     pari_sp av = avma;
    3240       13850 :     if (lg(T) == 5) /* quadratic fields */
    3241       10458 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3242             :     else
    3243             :     {
    3244        3392 :       GEN p = NULL;
    3245        3392 :       if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && p)
    3246             :       {
    3247          14 :         T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p);
    3248          14 :         if (lgefint(p) == 3)
    3249             :         {
    3250           7 :           ulong pp = p[2];
    3251           7 :           a = Flxq_inv(ZX_to_Flx(x, pp), ZX_to_Flx(T, pp), pp);
    3252           7 :           a = Flx_to_ZX(a);
    3253             :         }
    3254             :         else
    3255           7 :           a = FpXQ_inv(x, T, p);
    3256          14 :         a = FpX_to_mod(a, p);
    3257          14 :         a = gerepileupto(av, a);
    3258             :       }
    3259        3378 :       else if (RgX_is_ZX(T) && RgX_is_QX(x))
    3260             :       {
    3261        2650 :         avma = av;
    3262        2650 :         a = QXQ_inv(x, T);
    3263             :       }
    3264             :       else {
    3265         728 :         avma = av;
    3266         728 :         a = RgXQ_inv(x, T);
    3267             :       }
    3268             :     }
    3269             :   }
    3270       20990 :   gel(z,2) = a; return z;
    3271             : }
    3272             : 
    3273             : GEN
    3274    12492888 : ginv(GEN x)
    3275             : {
    3276             :   long s;
    3277             :   pari_sp av, tetpil;
    3278             :   GEN z, y, p1, p2;
    3279             : 
    3280    12492888 :   switch(typ(x))
    3281             :   {
    3282             :     case t_INT:
    3283       97886 :       if (is_pm1(x)) return icopy(x);
    3284       61675 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3285       61675 :       z = cgetg(3,t_FRAC);
    3286       61675 :       gel(z,1) = s<0? gen_m1: gen_1;
    3287       61675 :       gel(z,2) = absi(x); return z;
    3288             : 
    3289     1700992 :     case t_REAL: return invr(x);
    3290             : 
    3291         798 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3292         798 :       gel(z,1) = icopy(gel(x,1));
    3293         798 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3294             : 
    3295             :     case t_FRAC: {
    3296      112230 :       GEN a = gel(x,1), b = gel(x,2);
    3297      112230 :       s = signe(a);
    3298      112230 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3299       48592 :       z = cgetg(3,t_FRAC);
    3300       48592 :       gel(z,1) = icopy(b);
    3301       48592 :       gel(z,2) = icopy(a);
    3302       48592 :       normalize_frac(z); return z;
    3303             :     }
    3304             :     case t_COMPLEX:
    3305     1165884 :       av=avma;
    3306     1165884 :       p1=cxnorm(x);
    3307     1165884 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3308     1165884 :       tetpil=avma;
    3309     1165884 :       return gerepile(av,tetpil,divcR(p2,p1));
    3310             : 
    3311             :     case t_QUAD:
    3312         203 :       av=avma; p1=gnorm(x); p2=gconj(x); tetpil=avma;
    3313         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3314             : 
    3315        1099 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3316        1099 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3317        1092 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3318        1092 :       gel(z,2) = icopy(gel(x,2));
    3319        1092 :       gel(z,3) = icopy(gel(x,3));
    3320        1092 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3321             : 
    3322       20990 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3323        4618 :     case t_FFELT: return FF_inv(x);
    3324     5436717 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3325        7700 :     case t_SER: return gdiv(gen_1,x);
    3326             : 
    3327             :     case t_RFRAC:
    3328             :     {
    3329        2905 :       GEN n = gel(x,1), d = gel(x,2);
    3330        2905 :       pari_sp av = avma, ltop;
    3331        2905 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3332             : 
    3333        2905 :       n = simplify_shallow(n);
    3334        2905 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3335             :       {
    3336        2905 :         if (gequal1(n)) { avma = av; return RgX_copy(d); }
    3337         679 :         ltop = avma;
    3338         679 :         z = RgX_Rg_div(d,n);
    3339             :       } else {
    3340           0 :         ltop = avma;
    3341           0 :         z = cgetg(3,t_RFRAC);
    3342           0 :         gel(z,1) = RgX_copy(d);
    3343           0 :         gel(z,2) = RgX_copy(n);
    3344             :       }
    3345         679 :       stackdummy(av, ltop);
    3346         679 :       return z;
    3347             :     }
    3348             : 
    3349             :     case t_QFR:
    3350           7 :       av = avma; z = cgetg(5, t_QFR);
    3351           7 :       gel(z,1) = gel(x,1);
    3352           7 :       gel(z,2) = negi( gel(x,2) );
    3353           7 :       gel(z,3) = gel(x,3);
    3354           7 :       gel(z,4) = negr( gel(x,4) );
    3355           7 :       return gerepileupto(av, redreal(z));
    3356             : 
    3357             :     case t_QFI:
    3358     3939928 :       y = gcopy(x);
    3359     3939928 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3360     3841632 :         togglesign(gel(y,2));
    3361     3939928 :       return y;
    3362             :     case t_MAT:
    3363         931 :       y = RgM_inv(x);
    3364         931 :       if (!y) pari_err_INV("ginv",x);
    3365         931 :       return y;
    3366             :     case t_VECSMALL:
    3367             :     {
    3368           0 :       long i, lx = lg(x)-1;
    3369           0 :       y = zero_zv(lx);
    3370           0 :       for (i=1; i<=lx; i++)
    3371             :       {
    3372           0 :         long xi = x[i];
    3373           0 :         if (xi<1 || xi>lx || y[xi])
    3374           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3375           0 :         y[xi] = i;
    3376             :       }
    3377           0 :       return y;
    3378             :     }
    3379             :   }
    3380           0 :   pari_err_TYPE("inverse",x);
    3381           0 :   return NULL; /* not reached */
    3382             : }

Generated by: LCOV version 1.11