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.10.0 lcov report (development 19825-b77c7f8) Lines: 1785 1933 92.3 %
Date: 2016-12-05 05:49:04 Functions: 88 88 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    51969365 : addRc(GEN x, GEN y) {
      63    51969365 :   GEN z = cgetg(3,t_COMPLEX);
      64    51969365 :   gel(z,1) = gadd(x,gel(y,1));
      65    51969365 :   gel(z,2) = gcopy(gel(y,2)); return z;
      66             : }
      67             : static GEN
      68    52268059 : mulRc(GEN x, GEN y) {
      69    52268059 :   GEN z = cgetg(3,t_COMPLEX);
      70    52268059 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(x,gel(y,1));
      71    52268059 :   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      361365 : divRc(GEN x, GEN y) {
      82      361365 :   GEN t = gdiv(x, cxnorm(y)), mt = gneg(t); /* left on stack for efficiency */
      83      361365 :   GEN z = cgetg(3,t_COMPLEX);
      84      361365 :   gel(z,1) = isintzero(gel(y,1))? gen_0: gmul(t, gel(y,1));
      85      361365 :   gel(z,2) = gmul(mt, gel(y,2));
      86      361365 :   return z;
      87             : }
      88             : static GEN
      89     3526704 : divcR(GEN x, GEN y) {
      90     3526704 :   GEN z = cgetg(3,t_COMPLEX);
      91     3526704 :   gel(z,1) = isintzero(gel(x,1))? gen_0: gdiv(gel(x,1), y);
      92     3526704 :   gel(z,2) = gdiv(gel(x,2), y); return z;
      93             : }
      94             : static GEN
      95         490 : addRq(GEN x, GEN y) {
      96         490 :   GEN z = cgetg(4,t_QUAD);
      97         490 :   gel(z,1) = ZX_copy(gel(y,1));
      98         490 :   gel(z,2) = gadd(x, gel(y,2));
      99         490 :   gel(z,3) = gcopy(gel(y,3)); return z;
     100             : }
     101             : static GEN
     102        1624 : mulRq(GEN x, GEN y) {
     103        1624 :   GEN z = cgetg(4,t_QUAD);
     104        1624 :   gel(z,1) = ZX_copy(gel(y,1));
     105        1624 :   gel(z,2) = gmul(x,gel(y,2));
     106        1624 :   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     1583329 : mulrfrac(GEN x, GEN y)
     116             : {
     117     1583329 :   pari_sp av = avma;
     118     1583329 :   GEN z, a = gel(y,1), b = gel(y,2);
     119     1583329 :   if (is_pm1(a)) /* frequent special case */
     120             :   {
     121      917180 :     z = divri(x, b);
     122      917180 :     if (signe(a) < 0) togglesign(z);
     123      917180 :     return z;
     124             :   }
     125      666149 :   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      156661 : mulTp(GEN x, GEN y) { pari_sp av = avma;
     153      156661 :   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       14000 : divTp(GEN x, GEN y) { pari_sp av = avma;
     158       14000 :   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        1218 : divpT(GEN x, GEN y) { pari_sp av = avma;
     164        1218 :   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     1869748 : add_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     171     1869748 :   if (lgefint(X) == 3) {
     172     1834300 :     ulong u = Fl_add(itou(x),itou(y), X[2]);
     173     1834300 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     174             :   }
     175             :   else {
     176       35448 :     GEN u = addii(x,y); if (cmpii(u, X) >= 0) u = subii(u, X);
     177       35444 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     178             :   }
     179     1869759 :   gel(z,1) = icopy(X); return z;
     180             : }
     181             : static GEN
     182      262748 : sub_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     183      262748 :   if (lgefint(X) == 3) {
     184      242204 :     ulong u = Fl_sub(itou(x),itou(y), X[2]);
     185      242204 :     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       20539 :     gel(z,2) = gerepileuptoint((pari_sp)z, u);
     190             :   }
     191      262746 :   gel(z,1) = icopy(X); return z;
     192             : }
     193             : /* cf add_intmod_same */
     194             : static GEN
     195     3062379 : mul_intmod_same(GEN z, GEN X, GEN x, GEN y) {
     196     3062379 :   if (lgefint(X) == 3) {
     197     2947379 :     ulong u = Fl_mul(itou(x),itou(y), X[2]);
     198     2947379 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     199             :   }
     200             :   else
     201      115000 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x,y), X) );
     202     3062382 :   gel(z,1) = icopy(X); return z;
     203             : }
     204             : /* cf add_intmod_same */
     205             : static GEN
     206       29683 : div_intmod_same(GEN z, GEN X, GEN x, GEN y)
     207             : {
     208       29683 :   if (lgefint(X) == 3) {
     209       21756 :     ulong m = uel(X,2), u = Fl_div(itou(x), itou(y), m);
     210       21749 :     avma = (pari_sp)z; gel(z,2) = utoi(u);
     211             :   }
     212             :   else
     213        7927 :     gel(z,2) = gerepileuptoint((pari_sp)z, remii(mulii(x, Fp_inv(y,X)), X) );
     214       29677 :   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     9406168 : rfrac_denom_mul_scal(GEN d, GEN y)
     227             : {
     228     9406168 :   GEN D = RgX_Rg_mul(d, y);
     229     9406168 :   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     9406168 :   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    54283897 : gred_rfrac_simple(GEN n, GEN d)
     240             : {
     241             :   GEN c, cn, cd, z;
     242    54283897 :   long dd = degpol(d);
     243             : 
     244    54283897 :   if (dd <= 0)
     245             :   {
     246        5907 :     if (dd < 0) pari_err_INV("gred_rfrac_simple", d);
     247        5907 :     n = gdiv(n, gel(d,2));
     248        5907 :     if (typ(n) != t_POL || varn(n) != varn(d)) n = scalarpol(n, varn(d));
     249        5907 :     return n;
     250             :   }
     251             : 
     252    54277990 :   cd = content(d);
     253    54277990 :   while (typ(n) == t_POL && !degpol(n)) n = gel(n,2);
     254    54277990 :   cn = (typ(n) == t_POL && varn(n) == varn(d))? content(n): n;
     255    54277990 :   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    47729500 :     if (!gequal1(cn))
     272             :     {
     273     2587354 :       if (gequal0(cn)) {
     274         301 :         if (isexactzero(cn)) return scalarpol(cn, varn(d));
     275          21 :         c = gen_1;
     276             :       } else {
     277     2587053 :         n = (cn != n)? RgX_Rg_div(n,cn): gen_1;
     278     2587053 :         c = cn;
     279             :       }
     280             :     } else {
     281    45142146 :       GEN y = cgetg(3,t_RFRAC);
     282    45142146 :       gel(y,1) = gcopy(n);
     283    45142146 :       gel(y,2) = RgX_copy(d); return y;
     284             :     }
     285             :   }
     286             : 
     287     9135522 :   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     8230193 :     cn = numer(c);
     297     8230193 :     cd = denom(c);
     298             :   }
     299     9135522 :   z = cgetg(3,t_RFRAC);
     300     9135522 :   gel(z,1) = gmul(n, cn);
     301     9135522 :   gel(z,2) = rfrac_denom_mul_scal(d, cd);
     302     9135522 :   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    43370400 : gred_rfrac2(GEN n, GEN d)
     328             : {
     329             :   GEN y, z;
     330             :   long v, vd, vn;
     331             : 
     332    43370400 :   n = simplify_shallow(n);
     333    43370400 :   if (isrationalzero(n)) return gcopy(n);
     334    36576191 :   d = simplify_shallow(d);
     335    36576191 :   if (typ(d) != t_POL) return gdiv(n,d);
     336    35396694 :   vd = varn(d);
     337    35396694 :   if (typ(n) != t_POL)
     338             :   {
     339    19890725 :     if (varncmp(vd, gvar(n)) >= 0) return gdiv(n,d);
     340    19889318 :     if (varncmp(vd, gvar2(n)) < 0) return gred_rfrac_simple(n,d);
     341           0 :     pari_err_BUG("gred_rfrac2 [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             : /* x,y t_INT, return x/y in reduced form */
     367             : GEN
     368    56758193 : Qdivii(GEN x, GEN y)
     369             : {
     370    56758193 :   pari_sp av = avma;
     371             :   GEN r, q;
     372             : 
     373    56758193 :   if (is_pm1(y)) return (signe(y) < 0)? negi(x): icopy(x);
     374    22342706 :   if (is_pm1(x)) {
     375     4690825 :     long s = signe(y);
     376     4690825 :     if (!s) pari_err_INV("gdiv",y);
     377     4690741 :     if (signe(x) < 0) s = -s;
     378     4690741 :     q = cgetg(3, t_FRAC);
     379     4690741 :     gel(q,1) = s<0? gen_m1: gen_1;
     380     4690741 :     gel(q,2) = absi(y); return q;
     381             :   }
     382    17651881 :   q = dvmdii(x,y,&r);
     383    17651881 :   if (r == gen_0) return q; /* gen_0 intended */
     384     5245902 :   r = gcdii(y,r);
     385     5245902 :   if (lgefint(r) == 3)
     386             :   {
     387     5057260 :     ulong t = r[2];
     388     5057260 :     avma = av;
     389     5057260 :     if (t == 1) q = mkfraccopy(x,y);
     390             :     else
     391             :     {
     392     2618767 :       q = cgetg(3,t_FRAC);
     393     2618767 :       gel(q,1) = diviuexact(x,t);
     394     2618767 :       gel(q,2) = diviuexact(y,t);
     395             :     }
     396             :   }
     397             :   else
     398             :   { /* rare: r and q left on stack for efficiency */
     399      188642 :     q = cgetg(3,t_FRAC);
     400      188642 :     gel(q,1) = diviiexact(x,r);
     401      188642 :     gel(q,2) = diviiexact(y,r);
     402             :   }
     403     5245902 :   normalize_frac(q); return q;
     404             : }
     405             : 
     406             : /*******************************************************************/
     407             : /*                                                                 */
     408             : /*                          CONJUGATION                            */
     409             : /*                                                                 */
     410             : /*******************************************************************/
     411             : /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
     412             : static GEN
     413       11963 : quad_polmod_conj(GEN x, GEN y)
     414             : {
     415             :   GEN z, u, v, a, b;
     416             :   pari_sp av;
     417       11963 :   if (typ(x) != t_POL) return gcopy(x);
     418       11963 :   if (varn(x) != varn(y) || degpol(x) <= 0) return RgX_copy(x);
     419       11963 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     420       11963 :   b = gel(y,3); v = gel(x,2);
     421       11963 :   z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
     422       11963 :   gel(z,2) = gerepileupto(av, gsub(v, gdiv(gmul(u,b), a)));
     423       11963 :   gel(z,3) = gneg(u); return z;
     424             : }
     425             : static GEN
     426       11963 : quad_polmod_norm(GEN x, GEN y)
     427             : {
     428             :   GEN z, u, v, a, b, c;
     429             :   pari_sp av;
     430       11963 :   if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0)
     431           0 :     return gsqr(x);
     432       11963 :   a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
     433       11963 :   b = gel(y,3); v = gel(x,2);
     434       11963 :   c = gel(y,2); av = avma;
     435       11963 :   z = gmul(u, gsub(gmul(c,u), gmul(b,v)));
     436       11963 :   if (!gequal1(a)) z = gdiv(z, a);
     437       11963 :   return gerepileupto(av, gadd(z, gsqr(v)));
     438             : }
     439             : 
     440             : GEN
     441     1579670 : gconj(GEN x)
     442             : {
     443             :   long lx, i;
     444             :   GEN y;
     445             : 
     446     1579670 :   switch(typ(x))
     447             :   {
     448             :     case t_INT: case t_REAL:
     449       92287 :       return mpcopy(x);
     450             : 
     451             :     case t_INTMOD: case t_FRAC: case t_PADIC:
     452        1421 :       return gcopy(x);
     453             : 
     454             :     case t_COMPLEX:
     455     1476080 :       y = cgetg(3,t_COMPLEX);
     456     1476080 :       gel(y,1) = gcopy(gel(x,1));
     457     1476080 :       gel(y,2) = gneg(gel(x,2));
     458     1476080 :       break;
     459             : 
     460             :     case t_QUAD:
     461         595 :       y = cgetg(4,t_QUAD);
     462         595 :       gel(y,1) = ZX_copy(gel(x,1));
     463        1204 :       gel(y,2) = gequal0(gmael(x,1,3))? gcopy(gel(x,2))
     464         609 :                                     : gadd(gel(x,2), gel(x,3));
     465         595 :       gel(y,3) = gneg(gel(x,3));
     466         595 :       break;
     467             : 
     468             :     case t_POL: case t_SER:
     469        4732 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     470        4732 :       for (i=2; i<lx; i++) gel(y,i) = gconj(gel(x,i));
     471        4732 :       break;
     472             : 
     473             :     case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     474        4555 :       y = cgetg_copy(x, &lx);
     475        4555 :       for (i=1; i<lx; i++) gel(y,i) = gconj(gel(x,i));
     476        4555 :       break;
     477             : 
     478             :     case t_POLMOD:
     479             :     {
     480           0 :       GEN y, X = gel(x,1);
     481           0 :       long d = degpol(X);
     482           0 :       if (d < 2) return RgX_copy(x);
     483           0 :       if (d == 2) {
     484           0 :         y = cgetg(3, t_POLMOD);
     485           0 :         gel(y,1) = RgX_copy(X);
     486           0 :         gel(y,2) = quad_polmod_conj(gel(x,2), X); return y;
     487             :       }
     488             :     }
     489             :     default:
     490           0 :       pari_err_TYPE("gconj",x);
     491           0 :       return NULL; /* not reached */
     492             :   }
     493     1485962 :   return y;
     494             : }
     495             : 
     496             : GEN
     497          42 : conjvec(GEN x,long prec)
     498             : {
     499             :   long lx, s, i;
     500             :   GEN z;
     501             : 
     502          42 :   switch(typ(x))
     503             :   {
     504             :     case t_INT: case t_INTMOD: case t_FRAC:
     505           0 :       return mkcolcopy(x);
     506             : 
     507             :     case t_COMPLEX: case t_QUAD:
     508           0 :       z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
     509             : 
     510             :     case t_FFELT:
     511          28 :       return FF_conjvec(x);
     512             : 
     513             :     case t_VEC: case t_COL:
     514           0 :       lx = lg(x); z = cgetg(lx,t_MAT);
     515           0 :       if (lx == 1) return z;
     516           0 :       gel(z,1) = conjvec(gel(x,1),prec);
     517           0 :       s = lgcols(z);
     518           0 :       for (i=2; i<lx; i++)
     519             :       {
     520           0 :         gel(z,i) = conjvec(gel(x,i),prec);
     521           0 :         if (lg(gel(z,i)) != s) pari_err_OP("conjvec", gel(z,1), gel(z,i));
     522             :       }
     523           0 :       break;
     524             : 
     525             :     case t_POLMOD: {
     526          14 :       GEN T = gel(x,1), r;
     527             :       pari_sp av;
     528             : 
     529          14 :       lx = lg(T);
     530          14 :       if (lx <= 3) return cgetg(1,t_COL);
     531          14 :       x = gel(x,2);
     532          42 :       for (i=2; i<lx; i++)
     533             :       {
     534          35 :         GEN c = gel(T,i);
     535          35 :         switch(typ(c)) {
     536             :           case t_INTMOD: {
     537           7 :             GEN p = gel(c,1);
     538             :             pari_sp av;
     539           7 :             if (typ(x) != t_POL) retconst_col(lx-3, Rg_to_Fp(x, p));
     540           7 :             av = avma;
     541           7 :             T = RgX_to_FpX(T,p);
     542           7 :             x = RgX_to_FpX(x, p);
     543           7 :             if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     544           7 :             z = FpXQC_to_mod(FpXQ_conjvec(x, T , p), T, p);
     545           7 :             return gerepileupto(av, z);
     546             :           }
     547             :           case t_INT:
     548          28 :           case t_FRAC: break;
     549           0 :           default: pari_err_TYPE("conjvec [not a rational t_POL]",T);
     550             :         }
     551             :       }
     552           7 :       if (typ(x) != t_POL)
     553             :       {
     554           0 :         if (!is_rational_t(typ(x)))
     555           0 :           pari_err_TYPE("conjvec [not a rational t_POL]",x);
     556           0 :         retconst_col(lx-3, gcopy(x));
     557             :       }
     558           7 :       RgX_check_QX(x,"conjvec");
     559           7 :       av = avma;
     560           7 :       if (varn(x) != varn(T)) pari_err_VAR("conjvec",x,T);
     561           7 :       r = cleanroots(T,prec);
     562           7 :       z = cgetg(lx-2,t_COL);
     563           7 :       for (i=1; i<=lx-3; i++) gel(z,i) = poleval(x, gel(r,i));
     564           7 :       return gerepileupto(av, z);
     565             :     }
     566             : 
     567             :     default:
     568           0 :       pari_err_TYPE("conjvec",x);
     569           0 :       return NULL; /* not reached */
     570             :   }
     571           0 :   return z;
     572             : }
     573             : 
     574             : 
     575             : /********************************************************************/
     576             : /**                                                                **/
     577             : /**                           ADDITION                             **/
     578             : /**                                                                **/
     579             : /********************************************************************/
     580             : /* x, y compatible PADIC, op = add or sub */
     581             : static GEN
     582     1173655 : addsub_pp(GEN x, GEN y, GEN (*op)(GEN,GEN))
     583             : {
     584     1173655 :   pari_sp av = avma;
     585             :   long d,e,r,rx,ry;
     586     1173655 :   GEN u, z, mod, p = gel(x,2);
     587             :   int swap;
     588             : 
     589     1173655 :   (void)new_chunk(5 + lgefint(gel(x,3)) + lgefint(gel(y,3)));
     590     1173655 :   e = valp(x);
     591     1173655 :   r = valp(y); d = r-e;
     592     1173655 :   if (d < 0) { swap = 1; swap(x,y); e = r; d = -d; } else swap = 0;
     593     1173655 :   rx = precp(x);
     594     1173655 :   ry = precp(y);
     595     1173655 :   if (d) /* v(x) < v(y) */
     596             :   {
     597      548232 :     r = d+ry; z = powiu(p,d);
     598      548232 :     if (r < rx) mod = mulii(z,gel(y,3)); else { r = rx; mod = gel(x,3); }
     599      548232 :     z = mulii(z,gel(y,4));
     600      548232 :     u = swap? op(z, gel(x,4)): op(gel(x,4), z);
     601             :   }
     602             :   else
     603             :   {
     604             :     long c;
     605      625423 :     if (ry < rx) { r=ry; mod = gel(y,3); } else { r=rx; mod = gel(x,3); }
     606      625423 :     u = swap? op(gel(y,4), gel(x,4)): op(gel(x,4), gel(y,4));
     607      625423 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     608             :     {
     609       59199 :       avma = av; return zeropadic(p, e+r);
     610             :     }
     611      566224 :     if (c)
     612             :     {
     613       94886 :       mod = diviiexact(mod, powiu(p,c));
     614       94886 :       r -= c;
     615       94886 :       e += c;
     616             :     }
     617             :   }
     618     1114456 :   u = modii(u, mod);
     619     1114456 :   avma = av; z = cgetg(5,t_PADIC);
     620     1114456 :   z[1] = evalprecp(r) | evalvalp(e);
     621     1114456 :   gel(z,2) = icopy(p);
     622     1114456 :   gel(z,3) = icopy(mod);
     623     1114456 :   gel(z,4) = icopy(u); return z;
     624             : }
     625             : /* Rg_to_Fp(t_FRAC) without GC */
     626             : static GEN
     627       56510 : Q_to_Fp(GEN x, GEN p)
     628       56510 : { return mulii(gel(x,1), Fp_inv(gel(x,2),p)); }
     629             : /* return x + y, where y t_PADIC and x is a non-zero t_INT or t_FRAC */
     630             : static GEN
     631      359374 : addQp(GEN x, GEN y)
     632             : {
     633      359374 :   pari_sp av = avma;
     634      359374 :   long d, r, e, vy = valp(y), py = precp(y);
     635      359374 :   GEN z, q, mod, u, p = gel(y,2);
     636             : 
     637      359374 :   e = Q_pvalrem(x, p, &x);
     638      359370 :   d = vy - e; r = d + py;
     639      359370 :   if (r <= 0) { avma = av; return gcopy(y); }
     640      357844 :   mod = gel(y,3);
     641      357844 :   u   = gel(y,4);
     642      357844 :   (void)new_chunk(5 + ((lgefint(mod) + lgefint(p)*labs(d)) << 1));
     643             : 
     644      357847 :   if (d > 0)
     645             :   {
     646      233567 :     q = powiu(p,d);
     647      233568 :     mod = mulii(mod, q);
     648      233565 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     649      233565 :     u = addii(x,  mulii(u, q));
     650             :   }
     651      124280 :   else if (d < 0)
     652             :   {
     653       23138 :     q = powiu(p,-d);
     654       23138 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     655       23138 :     u = addii(u, mulii(x, q));
     656       23138 :     r = py; e = vy;
     657             :   }
     658             :   else
     659             :   {
     660             :     long c;
     661      101142 :     if (typ(x) != t_INT) x = Q_to_Fp(x, mod);
     662      101142 :     u = addii(u, x);
     663      101141 :     if (!signe(u) || (c = Z_pvalrem(u,p,&u)) >= r)
     664             :     {
     665         959 :       avma = av; return zeropadic(p,e+r);
     666             :     }
     667      100183 :     if (c)
     668             :     {
     669       55789 :       mod = diviiexact(mod, powiu(p,c));
     670       55789 :       r -= c;
     671       55789 :       e += c;
     672             :     }
     673             :   }
     674      356882 :   u = modii(u, mod); avma = av;
     675      356879 :   z = cgetg(5,t_PADIC);
     676      356887 :   z[1] = evalprecp(r) | evalvalp(e);
     677      356885 :   gel(z,2) = icopy(p);
     678      356887 :   gel(z,3) = icopy(mod);
     679      356886 :   gel(z,4) = icopy(u); return z;
     680             : }
     681             : 
     682             : /* Mod(x,X) + Mod(y,X) */
     683             : #define addsub_polmod_same addsub_polmod_scal
     684             : /* Mod(x,X) +/- Mod(y,Y) */
     685             : static GEN
     686         980 : addsub_polmod(GEN X, GEN Y, GEN x, GEN y, GEN(*op)(GEN,GEN))
     687             : {
     688         980 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
     689         980 :   GEN z = cgetg(3,t_POLMOD);
     690         980 :   long vx = varn(X), vy = varn(Y);
     691         980 :   if (vx==vy) {
     692             :     pari_sp av;
     693          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
     694          14 :     warn_coercion(X,Y,gel(z,1));
     695          14 :     gel(z,2) = gerepileupto(av, gmod(op(x, y), gel(z,1))); return z;
     696             :   }
     697         966 :   if (varncmp(vx, vy) < 0)
     698         966 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
     699             :   else
     700           0 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
     701         966 :   gel(z,2) = op(x, y); return z;
     702             : }
     703             : /* Mod(y, Y) +/- x,  x scalar or polynomial in same var and reduced degree */
     704             : static GEN
     705     2593484 : addsub_polmod_scal(GEN Y, GEN y, GEN x, GEN(*op)(GEN,GEN))
     706             : {
     707     2593484 :   GEN z = cgetg(3,t_POLMOD);
     708     2593484 :   gel(z,1) = RgX_copy(Y);
     709     2593484 :   gel(z,2) = op(y, x); return z;
     710             : }
     711             : 
     712             : /* typ(y) == t_SER, x "scalar" [e.g object in lower variable] */
     713             : static GEN
     714      201597 : add_ser_scal(GEN y, GEN x)
     715             : {
     716             :   long i, l, ly, vy;
     717             :   GEN z;
     718             : 
     719      201597 :   if (isrationalzero(x)) return gcopy(y);
     720      197040 :   ly = lg(y);
     721      197040 :   l = valp(y);
     722      197040 :   if (l < 3-ly) return gcopy(y);
     723             :   /* l + ly >= 3 */
     724      196921 :   if (l < 0)
     725             :   {
     726         693 :     z = cgetg(ly,t_SER); z[1] = y[1];
     727         693 :     for (i = 2; i <= 1-l; i++) gel(z,i) = gcopy(gel(y,i));
     728         693 :     gel(z,i) = gadd(x,gel(y,i)); i++;
     729         693 :     for (     ; i < ly; i++)   gel(z,i) = gcopy(gel(y,i));
     730         693 :     return z;
     731             :   }
     732      196228 :   vy = varn(y);
     733      196228 :   if (l > 0)
     734             :   {
     735        8267 :     if (ser_isexactzero(y))
     736        5838 :       return scalarser(ly == 2? x: gadd(x,gel(y,2)), vy, l);
     737        2429 :     y -= l; ly += l;
     738        2429 :     z = cgetg(ly,t_SER);
     739        2429 :     x = gcopy(x);
     740        2429 :     for (i=3; i<=l+1; i++) gel(z,i) = gen_0;
     741             :   }
     742             :   else
     743             :   { /* l = 0, ly >= 3. Also OK if ser_isexactzero(y) */
     744      187961 :     z = cgetg(ly,t_SER);
     745      187961 :     x = gadd(x, gel(y,2));
     746      187961 :     i = 3;
     747             :   }
     748      190390 :   for (; i<ly; i++) gel(z,i) = gcopy(gel(y,i));
     749      190390 :   gel(z,2) = x;
     750      190390 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(vy);
     751      190390 :   return gequal0(x)? normalize(z): z;
     752             : }
     753             : static long
     754     2339594 : _serprec(GEN x) { return ser_isexactzero(x)? 2: lg(x); }
     755             : /* x,y t_SER in the same variable: x+y */
     756             : static GEN
     757     1170168 : ser_add(GEN x, GEN y)
     758             : {
     759     1170168 :   long i, lx,ly, n = valp(y) - valp(x);
     760             :   GEN z;
     761     1170168 :   if (n < 0) { n = -n; swap(x,y); }
     762             :   /* valp(x) <= valp(y) */
     763     1170168 :   lx = _serprec(x);
     764     1170168 :   if (lx == 2) /* don't lose type information */
     765         742 :     return scalarser(gadd(RgX_get_0(x), RgX_get_0(y)), varn(x), valp(x));
     766     1169426 :   ly = _serprec(y) + n; if (lx < ly) ly = lx;
     767     1169426 :   if (n)
     768             :   {
     769        9387 :     if (n+2 > lx) return gcopy(x);
     770        8687 :     z = cgetg(ly,t_SER);
     771        8687 :     for (i=2; i<=n+1; i++) gel(z,i) = gcopy(gel(x,i));
     772        8687 :     for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-n));
     773             :   } else {
     774     1160039 :     z = cgetg(ly,t_SER);
     775     1160039 :     for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
     776             :   }
     777     1168726 :   z[1] = x[1]; return normalize(z);
     778             : }
     779             : /* typ(y) == RFRAC, x polynomial in same variable or "scalar" */
     780             : static GEN
     781     8350089 : add_rfrac_scal(GEN y, GEN x)
     782             : {
     783             :   pari_sp av;
     784             :   GEN n;
     785             : 
     786     8350089 :   if (isintzero(x)) return gcopy(y); /* frequent special case */
     787     4625077 :   av = avma; n = gadd(gmul(x, gel(y,2)), gel(y,1));
     788     4625077 :   return gerepileupto(av, gred_rfrac_simple(n, gel(y,2)));
     789             : }
     790             : 
     791             : /* x "scalar", ty != t_MAT and non-scalar */
     792             : static GEN
     793    14195522 : add_scal(GEN y, GEN x, long ty)
     794             : {
     795    14195522 :   switch(ty)
     796             :   {
     797     9542712 :     case t_POL: return RgX_Rg_add(y, x);
     798      201597 :     case t_SER: return add_ser_scal(y, x);
     799     4241003 :     case t_RFRAC: return add_rfrac_scal(y, x);
     800           0 :     case t_COL: return RgC_Rg_add(y, x);
     801             :     case t_VEC:
     802      210210 :       if (isintzero(x)) return gcopy(y);
     803         168 :       break;
     804             :   }
     805         168 :   pari_err_TYPE2("+",x,y);
     806           0 :   return NULL; /* not reached */
     807             : }
     808             : 
     809             : static GEN
     810    25152108 : addsub_frac(GEN x, GEN y, GEN (*op)(GEN,GEN))
     811             : {
     812    25152108 :   pari_sp av = avma;
     813    25152108 :   GEN x1 = gel(x,1), x2 = gel(x,2), z = cgetg(3,t_FRAC);
     814    25152108 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     815    25152108 :   int s = cmpii(x2, y2);
     816             : 
     817    25152108 :   if (!s)
     818             :   { /* common denominator: (x1 op y1) / x2 */
     819    11588180 :     n = op(x1, y1);
     820    11588180 :     if (!signe(n)) { avma = av; return gen_0; }
     821    11033111 :     d = x2;
     822    11033111 :     q = dvmdii(n, d, &r);
     823    11033111 :     if (r == gen_0) { avma = av; return icopy(q); }
     824    10052991 :     r = gcdii(d, r);
     825    10052991 :     if (!is_pm1(r)) { n = diviiexact(n, r); d = diviiexact(d, r); }
     826    10052991 :     gel(z,1) = icopy_avma(n, (pari_sp)z);
     827    10052991 :     gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     828    10052991 :     avma = (pari_sp)gel(z,2); return z;
     829             :   }
     830    13563928 :   if (s < 0)
     831             :   {
     832     5867840 :     GEN Q = dvmdii(y2, x2, &r);
     833     5867840 :     if (r == gen_0)
     834             :     { /* y2 = Q x2: 1/x2 . (Q x1 op y1)/Q, where latter is in coprime form */
     835     4740646 :       pari_sp av = avma;
     836     4740646 :       n = op(mulii(Q,x1), y1);
     837     4740646 :       q = dvmdii(n, x2, &r);
     838     4740646 :       if (r == gen_0)
     839             :       {
     840      708589 :         gel(z,1) = gerepileuptoint(av, q);
     841      708589 :         gel(z,2) = Q; return z;
     842             :       }
     843     4032057 :       r = gcdii(x2, r);
     844     4032057 :       if (!is_pm1(r)) { n = diviiexact(n, r); x2 = diviiexact(x2, r); }
     845     4032057 :       d = mulii(x2,Q);
     846     4032057 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     847     4032057 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     848     4032057 :       avma = (pari_sp)gel(z,2); return z;
     849             :     }
     850     1127194 :     delta = gcdii(x2,r);
     851             :   }
     852             :   else
     853             :   {
     854     7696088 :     GEN Q = dvmdii(x2, y2, &r);
     855     7696088 :     if (r == gen_0)
     856             :     { /* x2 = Q y2: 1/y2 . (x1 op Q y1)/Q, where latter is in coprime form */
     857     6102024 :       pari_sp av = avma;
     858     6102024 :       n = op(x1, mulii(Q,y1));
     859     6102024 :       q = dvmdii(n, y2, &r);
     860     6102024 :       if (r == gen_0)
     861             :       {
     862      121749 :         gel(z,1) = gerepileuptoint(av, q);
     863      121749 :         gel(z,2) = Q; return z;
     864             :       }
     865     5980275 :       r = gcdii(y2, r);
     866     5980275 :       if (!is_pm1(r)) { n = diviiexact(n, r); y2 = diviiexact(y2, r); }
     867     5980275 :       d = mulii(y2,Q);
     868     5980275 :       gel(z,1) = icopy_avma(n, (pari_sp)z);
     869     5980275 :       gel(z,2) = icopy_avma(d, (pari_sp)gel(z,1));
     870     5980275 :       avma = (pari_sp)gel(z,2); return z;
     871             :     }
     872     1594064 :     delta = gcdii(y2,r);
     873             :   }
     874             :   /* delta = gcd(x2,y2) */
     875     2721258 :   if (is_pm1(delta))
     876             :   { /* numerator is non-zero */
     877      206306 :     gel(z,1) = gerepileuptoint((pari_sp)z, op(mulii(x1,y2), mulii(y1,x2)));
     878      206306 :     gel(z,2) = mulii(x2,y2); return z;
     879             :   }
     880     2514952 :   x2 = diviiexact(x2,delta);
     881     2514952 :   y2 = diviiexact(y2,delta);
     882     2514952 :   n = op(mulii(x1,y2), mulii(y1,x2));
     883     2514952 :   if (!signe(n)) { avma = av; return gen_0; }
     884     2514952 :   d = mulii(x2, y2);
     885     2514952 :   q = dvmdii(n, delta, &r);
     886     2514952 :   if (r == gen_0)
     887             :   {
     888      142498 :     if (is_pm1(d)) { avma = av; return icopy(q); }
     889      142498 :     avma = (pari_sp)z;
     890      142498 :     gel(z,2) = icopy(d);
     891      142498 :     gel(z,1) = icopy(q); return z;
     892             :   }
     893     2372454 :   r = gcdii(delta, r);
     894     2372454 :   if (!is_pm1(r))
     895             :   {
     896      802771 :     n     = diviiexact(n, r);
     897      802771 :     delta = diviiexact(delta, r);
     898             :   }
     899     2372454 :   d = mulii(d,delta); avma = (pari_sp)z;
     900     2372454 :   gel(z,1) = icopy(n);
     901     2372454 :   gel(z,2) = icopy(d); return z;
     902             : }
     903             : 
     904             : /* assume x2, y2 are t_POLs in the same variable */
     905             : static GEN
     906     2781571 : add_rfrac(GEN x, GEN y)
     907             : {
     908     2781571 :   pari_sp av = avma;
     909     2781571 :   GEN x1 = gel(x,1), x2 = gel(x,2);
     910     2781571 :   GEN y1 = gel(y,1), y2 = gel(y,2), q, r, n, d, delta;
     911             : 
     912     2781571 :   delta = RgX_gcd(x2,y2);
     913     2781571 :   if (!degpol(delta))
     914             :   {
     915         616 :     n = simplify_shallow( gadd(gmul(x1,y2), gmul(y1,x2)) );
     916         616 :     d = RgX_mul(x2, y2);
     917         616 :     return gerepileupto(av, gred_rfrac_simple(n, d));
     918             :   }
     919     2780955 :   x2 = RgX_div(x2,delta);
     920     2780955 :   y2 = RgX_div(y2,delta);
     921     2780955 :   n = gadd(gmul(x1,y2), gmul(y1,x2));
     922     2780955 :   if (!signe(n))
     923             :   {
     924      680797 :     n = simplify_shallow(n);
     925      680797 :     if (isrationalzero(n)) return gerepileupto(av, n);
     926           7 :     return gerepilecopy(av, mkrfrac(n, RgX_mul(gel(x,2),y2)));
     927             :   }
     928     2100158 :   if (degpol(n) == 0)
     929     1138189 :     return gerepileupto(av, gred_rfrac_simple(gel(n,2), RgX_mul(gel(x,2),y2)));
     930      961969 :   q = RgX_divrem(n, delta, &r); /* we want gcd(n,delta) */
     931      961969 :   if (isexactzero(r))
     932             :   {
     933             :     GEN z;
     934      208991 :     d = RgX_mul(x2, y2);
     935             :     /* "constant" denominator ? */
     936      208991 :     z = lg(d) == 3? RgX_Rg_div(q, gel(d,2)): gred_rfrac_simple(q, d);
     937      208991 :     return gerepileupto(av, z);
     938             :   }
     939      752978 :   r = RgX_gcd(delta, r);
     940      752978 :   if (degpol(r))
     941             :   {
     942      148378 :     n = RgX_div(n, r);
     943      148378 :     d = RgX_mul(RgX_mul(x2,y2), RgX_div(delta, r));
     944             :   }
     945             :   else
     946      604600 :     d = RgX_mul(gel(x,2), y2);
     947      752978 :   return gerepileupto(av, gred_rfrac_simple(n, d));
     948             : }
     949             : 
     950             : GEN
     951  1705895431 : gadd(GEN x, GEN y)
     952             : {
     953  1705895431 :   long tx = typ(x), ty = typ(y), vx, vy, lx, i, l;
     954             :   pari_sp av, tetpil;
     955             :   GEN z, p1;
     956             : 
     957  1705895431 :   if (tx == ty) switch(tx) /* shortcut to generic case */
     958             :   {
     959  1090812037 :     case t_INT: return addii(x,y);
     960   336229772 :     case t_REAL: return addrr(x,y);
     961     1402911 :     case t_INTMOD:  { GEN X = gel(x,1), Y = gel(y,1);
     962     1402911 :       z = cgetg(3,t_INTMOD);
     963     1402900 :       if (X==Y || equalii(X,Y))
     964     1402889 :         return add_intmod_same(z, X, gel(x,2), gel(y,2));
     965          14 :       gel(z,1) = gcdii(X,Y);
     966          14 :       warn_coercion(X,Y,gel(z,1));
     967          14 :       av = avma; p1 = addii(gel(x,2),gel(y,2));
     968          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
     969             :     }
     970    17746462 :     case t_FRAC: return addsub_frac(x,y,addii);
     971    37056294 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
     972    37056294 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     973    37056294 :       if (isintzero(gel(z,2)))
     974             :       {
     975         672 :         avma = (pari_sp)(z+3);
     976         672 :         return gadd(gel(x,1),gel(y,1));
     977             :       }
     978    37055622 :       gel(z,1) = gadd(gel(x,1),gel(y,1));
     979    37055622 :       return z;
     980             :     case t_PADIC:
     981      944832 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
     982      944832 :       return addsub_pp(x,y, addii);
     983         476 :     case t_QUAD: z = cgetg(4,t_QUAD);
     984         476 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
     985         476 :       gel(z,1) = ZX_copy(gel(x,1));
     986         476 :       gel(z,2) = gadd(gel(x,2),gel(y,2));
     987         476 :       gel(z,3) = gadd(gel(x,3),gel(y,3)); return z;
     988             :     case t_POLMOD:
     989     1164921 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
     990     1163976 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gadd);
     991         945 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gadd);
     992    10332878 :     case t_FFELT: return FF_add(x,y);
     993             :     case t_POL:
     994    14828859 :       vx = varn(x);
     995    14828859 :       vy = varn(y);
     996    14828859 :       if (vx != vy) {
     997      802832 :         if (varncmp(vx, vy) < 0) return RgX_Rg_add(x, y);
     998       20143 :         else                     return RgX_Rg_add(y, x);
     999             :       }
    1000    14026027 :       return RgX_add(x, y);
    1001             :     case t_SER:
    1002     1167193 :       vx = varn(x);
    1003     1167193 :       vy = varn(y);
    1004     1167193 :       if (vx != vy) {
    1005           0 :         if (varncmp(vx, vy) < 0) return add_ser_scal(x, y);
    1006           0 :         else                     return add_ser_scal(y, x);
    1007             :       }
    1008     1167193 :       return ser_add(x, y);
    1009             :     case t_RFRAC:
    1010     4073474 :       vx = varn(gel(x,2));
    1011     4073474 :       vy = varn(gel(y,2));
    1012     4073474 :       if (vx != vy) {
    1013     1291903 :         if (varncmp(vx, vy) < 0) return add_rfrac_scal(x, y);
    1014      536802 :         else                     return add_rfrac_scal(y, x);
    1015             :       }
    1016     2781571 :       return add_rfrac(x,y);
    1017             :     case t_VEC:
    1018      478005 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1019      478005 :       return RgV_add(x,y);
    1020             :     case t_COL:
    1021     4429083 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1022     4429083 :       return RgC_add(x,y);
    1023             :     case t_MAT:
    1024      495145 :       lx = lg(x);
    1025      495145 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1026      495145 :       if (lx == 1) return cgetg(1, t_MAT);
    1027      495145 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1028      495138 :       return RgM_add(x,y);
    1029             : 
    1030           0 :     default: pari_err_TYPE2("+",x,y);
    1031             :   }
    1032             :   /* tx != ty */
    1033   184734234 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1034             : 
    1035   184734234 :   if (is_const_t(ty)) switch(tx) /* tx < ty, is_const_t(tx) && is_const_t(ty) */
    1036             :   {
    1037             :     case t_INT:
    1038   150578788 :       switch(ty)
    1039             :       {
    1040    89178473 :         case t_REAL: return addir(x,y);
    1041             :         case t_INTMOD:
    1042      463295 :           z = cgetg(3, t_INTMOD);
    1043      463295 :           return add_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1044    22488677 :         case t_FRAC: z = cgetg(3,t_FRAC);
    1045    22488677 :           gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulii(gel(y,2),x)));
    1046    22488677 :           gel(z,2) = icopy(gel(y,2)); return z;
    1047    36559612 :         case t_COMPLEX: return addRc(x, y);
    1048             :         case t_PADIC:
    1049      319188 :           if (!signe(x)) return gcopy(y);
    1050      288388 :           return addQp(x,y);
    1051         469 :         case t_QUAD: return addRq(x, y);
    1052     1569075 :         case t_FFELT: return FF_Z_add(y,x);
    1053             :       }
    1054             : 
    1055             :     case t_REAL:
    1056    16176284 :       switch(ty)
    1057             :       {
    1058             :         case t_FRAC:
    1059      994164 :           if (!signe(gel(y,1))) return rcopy(x);
    1060      994164 :           if (!signe(x))
    1061             :           {
    1062       11963 :             lx = expi(gel(y,1)) - expi(gel(y,2)) - expo(x);
    1063       11963 :             return lx <= 0? rcopy(x): fractor(y, nbits2prec(lx));
    1064             :           }
    1065      982201 :           av=avma; z=addir(gel(y,1),mulir(gel(y,2),x)); tetpil=avma;
    1066      982201 :           return gerepile(av,tetpil,divri(z,gel(y,2)));
    1067    15182106 :         case t_COMPLEX: return addRc(x, y);
    1068          14 :         case t_QUAD: return gequal0(y)? rcopy(x): addqf(y, x, lg(x));
    1069             : 
    1070           0 :         default: pari_err_TYPE2("+",x,y);
    1071             :       }
    1072             : 
    1073             :     case t_INTMOD:
    1074        4006 :       switch(ty)
    1075             :       {
    1076        3500 :         case t_FRAC: { GEN X = gel(x,1);
    1077        3500 :           z = cgetg(3, t_INTMOD);
    1078        3500 :           p1 = Fp_div(gel(y,1), gel(y,2), X);
    1079        3500 :           return add_intmod_same(z, X, p1, gel(x,2));
    1080             :         }
    1081             :         case t_FFELT:
    1082         436 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1083           0 :             pari_err_OP("+",x,y);
    1084         436 :           return FF_Z_add(y,gel(x,2));
    1085          56 :         case t_COMPLEX: return addRc(x, y);
    1086           0 :         case t_PADIC: { GEN X = gel(x,1);
    1087           0 :           z = cgetg(3, t_INTMOD);
    1088           0 :           return add_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1089             :         }
    1090          14 :         case t_QUAD: return addRq(x, y);
    1091             :       }
    1092             : 
    1093             :     case t_FRAC:
    1094      300019 :       switch (ty)
    1095             :       {
    1096      227584 :         case t_COMPLEX: return addRc(x, y);
    1097             :         case t_PADIC:
    1098       70986 :           if (!signe(gel(x,1))) return gcopy(y);
    1099       70986 :           return addQp(x,y);
    1100           7 :         case t_QUAD: return addRq(x, y);
    1101        1442 :         case t_FFELT: return FF_Q_add(y, x);
    1102             :       }
    1103             : 
    1104             :     case t_FFELT:
    1105           0 :       pari_err_TYPE2("+",x,y);
    1106             : 
    1107             :     case t_COMPLEX:
    1108          21 :       switch(ty)
    1109             :       {
    1110             :         case t_PADIC:
    1111          14 :           return Zp_nosquare_m1(gel(y,2))? addRc(y, x): addTp(x, y);
    1112             :         case t_QUAD:
    1113           7 :           lx = precision(x); if (!lx) pari_err_OP("+",x,y);
    1114           7 :           return gequal0(y)? gcopy(x): addqf(y, x, lx);
    1115             :       }
    1116             : 
    1117             :     case t_PADIC: /* ty == t_QUAD */
    1118           0 :       return (kro_quad(y,gel(x,2)) == -1)? addRq(x, y): addTp(y, x);
    1119             :   }
    1120             :   /* tx < ty, !is_const_t(y) */
    1121    17675117 :   switch(ty)
    1122             :   {
    1123             :     case t_MAT:
    1124       16233 :       if (is_matvec_t(tx)) pari_err_TYPE2("+",x,y);
    1125       16233 :       if (isrationalzero(x)) return gcopy(y);
    1126       16107 :       return RgM_Rg_add(y, x);
    1127             :     case t_COL:
    1128       39027 :       if (tx == t_VEC) pari_err_TYPE2("+",x,y);
    1129       39027 :       return RgC_Rg_add(y, x);
    1130             :     case t_POLMOD: /* is_const_t(tx) in this case */
    1131      592711 :       return addsub_polmod_scal(gel(y,1), gel(y,2), x, &gadd);
    1132             :   }
    1133    17027146 :   if (is_scalar_t(tx))  {
    1134    14203316 :     if (tx == t_POLMOD)
    1135             :     {
    1136       30667 :       vx = varn(gel(x,1));
    1137       30667 :       vy = gvar(y);
    1138       30667 :       if (vx == vy) y = gmod(y, gel(x,1)); /* error if ty == t_SER */
    1139             :       else
    1140       26754 :         if (varncmp(vx,vy) > 0) return add_scal(y, x, ty);
    1141       10885 :       return addsub_polmod_scal(gel(x,1), gel(x,2), y, &gadd);
    1142             :     }
    1143    14172649 :     return add_scal(y, x, ty);
    1144             :   }
    1145             :   /* x and y are not scalars, ty != t_MAT */
    1146     2823830 :   vx = gvar(x);
    1147     2823830 :   vy = gvar(y);
    1148     2823830 :   if (vx != vy) { /* x or y is treated as a scalar */
    1149        3098 :     if (is_vec_t(tx) || is_vec_t(ty)) pari_err_TYPE2("+",x,y);
    1150        6182 :     return (varncmp(vx, vy) < 0)? add_scal(x, y, tx)
    1151        3091 :                                 : add_scal(y, x, ty);
    1152             :   }
    1153             :   /* vx = vy */
    1154     2820732 :   switch(tx)
    1155             :   {
    1156             :     case t_POL:
    1157     2820179 :       switch (ty)
    1158             :       {
    1159             :         case t_SER:
    1160        2996 :           if (lg(x) == 2) return gcopy(y);
    1161        2982 :           i = lg(y) + valp(y) - RgX_val(x);
    1162        2982 :           if (i < 3) return gcopy(y);
    1163             : 
    1164        2975 :           p1 = RgX_to_ser(x,i); y = ser_add(p1,y);
    1165        2975 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    1166        2975 :           return y;
    1167             : 
    1168     2817183 :         case t_RFRAC: return add_rfrac_scal(y, x);
    1169             :       }
    1170           0 :       break;
    1171             : 
    1172             :     case t_SER:
    1173         553 :       if (ty == t_RFRAC)
    1174             :       {
    1175             :         GEN n, d;
    1176             :         long vn, vd;
    1177         553 :         av = avma;
    1178         553 :         n = gel(y,1); vn = gval(n, vy);
    1179         553 :         d = gel(y,2); vd = RgX_valrem(d, &d);
    1180             : 
    1181         553 :         l = lg(x) + valp(x) - (vn - vd);
    1182         553 :         if (l < 3) { avma = av; return gcopy(x); }
    1183             : 
    1184             :         /* take advantage of y = t^n ! */
    1185         553 :         if (degpol(d))
    1186          63 :           y = gdiv(n, RgX_to_ser_inexact(d,l));
    1187             :         else {
    1188         490 :           y = gdiv(n, gel(d,2));
    1189         490 :           if (gvar(y) == vy) y = RgX_to_ser(y,l); else y = scalarser(y, vy, l);
    1190             :         }
    1191         553 :         setvalp(y, valp(y) - vd);
    1192         553 :         return gerepileupto(av, gadd(y, x));
    1193             :       }
    1194           0 :       break;
    1195             :   }
    1196           0 :   pari_err_TYPE2("+",x,y);
    1197           0 :   return NULL; /* not reached */
    1198             : }
    1199             : 
    1200             : GEN
    1201    15462969 : gaddsg(long x, GEN y)
    1202             : {
    1203    15462969 :   long ty = typ(y);
    1204             :   GEN z;
    1205             : 
    1206    15462969 :   switch(ty)
    1207             :   {
    1208      779357 :     case t_INT:  return addsi(x,y);
    1209     7020338 :     case t_REAL: return addsr(x,y);
    1210             :     case t_INTMOD:
    1211          14 :       z = cgetg(3, t_INTMOD);
    1212          14 :       return add_intmod_same(z, gel(y,1), gel(y,2), modsi(x, gel(y,1)));
    1213      383271 :     case t_FRAC: z = cgetg(3,t_FRAC);
    1214      383271 :       gel(z,1) = gerepileuptoint((pari_sp)z, addii(gel(y,1), mulis(gel(y,2),x)));
    1215      383271 :       gel(z,2) = icopy(gel(y,2)); return z;
    1216             :     case t_COMPLEX:
    1217     6826353 :       z = cgetg(3, t_COMPLEX);
    1218     6826353 :       gel(z,1) = gaddsg(x, gel(y,1));
    1219     6826353 :       gel(z,2) = gcopy(gel(y,2)); return z;
    1220             : 
    1221      453636 :     default: return gadd(stoi(x), y);
    1222             :   }
    1223             : }
    1224             : 
    1225             : GEN
    1226      411436 : gsubsg(long x, GEN y)
    1227             : {
    1228             :   GEN z, a, b;
    1229             :   pari_sp av;
    1230             : 
    1231      411436 :   switch(typ(y))
    1232             :   {
    1233       37573 :     case t_INT:  return subsi(x,y);
    1234      292797 :     case t_REAL: return subsr(x,y);
    1235             :     case t_INTMOD:
    1236          49 :       z = cgetg(3, t_INTMOD); a = gel(y,1); b = gel(y,2);
    1237          49 :       return add_intmod_same(z, a, Fp_neg(b,a), modsi(x, a));
    1238        7706 :     case t_FRAC: z = cgetg(3,t_FRAC); a = gel(y,1); b = gel(y,2);
    1239        7706 :       gel(z,1) = gerepileuptoint((pari_sp)z, subii(mulis(b,x), a));
    1240        7706 :       gel(z,2) = icopy(gel(y,2)); return z;
    1241             :     case t_COMPLEX:
    1242       46550 :       z = cgetg(3, t_COMPLEX);
    1243       46550 :       gel(z,1) = gsubsg(x, gel(y,1));
    1244       46550 :       gel(z,2) = gneg(gel(y,2)); return z;
    1245             :   }
    1246       26761 :   av = avma;
    1247       26761 :   return gerepileupto(av, gadd(stoi(x), gneg_i(y)));
    1248             : }
    1249             : 
    1250             : /********************************************************************/
    1251             : /**                                                                **/
    1252             : /**                          SUBTRACTION                           **/
    1253             : /**                                                                **/
    1254             : /********************************************************************/
    1255             : 
    1256             : GEN
    1257  1456557629 : gsub(GEN x, GEN y)
    1258             : {
    1259  1456557629 :   long tx = typ(x), ty = typ(y);
    1260             :   pari_sp av;
    1261             :   GEN z;
    1262  1456557629 :   if (tx == ty) switch(tx) /* shortcut to generic case */
    1263             :   {
    1264  1261534748 :     case t_INT: return subii(x,y);
    1265   145063554 :     case t_REAL: return subrr(x,y);
    1266      262762 :     case t_INTMOD:  { GEN p1, X = gel(x,1), Y = gel(y,1);
    1267      262762 :       z = cgetg(3,t_INTMOD);
    1268      262762 :       if (X==Y || equalii(X,Y))
    1269      262748 :         return sub_intmod_same(z, X, gel(x,2), gel(y,2));
    1270          14 :       gel(z,1) = gcdii(X,Y);
    1271          14 :       warn_coercion(X,Y,gel(z,1));
    1272          14 :       av = avma; p1 = subii(gel(x,2),gel(y,2));
    1273          14 :       gel(z,2) = gerepileuptoint(av, modii(p1, gel(z,1))); return z;
    1274             :     }
    1275     7405646 :     case t_FRAC: return addsub_frac(x,y, subii);
    1276    17352150 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    1277    17352150 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1278    17352150 :       if (isintzero(gel(z,2)))
    1279             :       {
    1280        7378 :         avma = (pari_sp)(z+3);
    1281        7378 :         return gsub(gel(x,1),gel(y,1));
    1282             :       }
    1283    17344772 :       gel(z,1) = gsub(gel(x,1),gel(y,1));
    1284    17344772 :       return z;
    1285             :     case t_PADIC:
    1286      228816 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("+",x,y);
    1287      228816 :       return addsub_pp(x,y, subii);
    1288         581 :     case t_QUAD: z = cgetg(4,t_QUAD);
    1289         581 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("+",x,y);
    1290         581 :       gel(z,1) = ZX_copy(gel(x,1));
    1291         581 :       gel(z,2) = gsub(gel(x,2),gel(y,2));
    1292         581 :       gel(z,3) = gsub(gel(x,3),gel(y,3)); return z;
    1293             :     case t_POLMOD:
    1294      825947 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1295      825912 :         return addsub_polmod_same(gel(x,1), gel(x,2), gel(y,2), &gsub);
    1296          35 :       return addsub_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2), &gsub);
    1297      473591 :     case t_FFELT: return FF_sub(x,y);
    1298             :     case t_POL: {
    1299     6062731 :       long vx = varn(x);
    1300     6062731 :       long vy = varn(y);
    1301     6062731 :       if (vx != vy) {
    1302       14127 :         if (varncmp(vx, vy) < 0) return RgX_Rg_sub(x, y);
    1303        1939 :         else                     return Rg_RgX_sub(x, y);
    1304             :       }
    1305     6048604 :       return RgX_sub(x, y);
    1306             :     }
    1307             :     case t_VEC:
    1308       65637 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1309       65637 :       return RgV_sub(x,y);
    1310             :     case t_COL:
    1311     1270486 :       if (lg(y) != lg(x)) pari_err_OP("+",x,y);
    1312     1270486 :       return RgC_sub(x,y);
    1313             :     case t_MAT: {
    1314        8437 :       long lx = lg(x);
    1315        8437 :       if (lg(y) != lx) pari_err_OP("+",x,y);
    1316        8437 :       if (lx == 1) return cgetg(1, t_MAT);
    1317        8437 :       if (lgcols(y) != lgcols(x)) pari_err_OP("+",x,y);
    1318        8437 :       return RgM_sub(x,y);
    1319             :     }
    1320     1889391 :     case t_RFRAC: case t_SER: break;
    1321             : 
    1322           0 :     default: pari_err_TYPE2("+",x,y);
    1323             :   }
    1324    16002543 :   av = avma;
    1325    16002543 :   return gerepileupto(av, gadd(x,gneg_i(y)));
    1326             : }
    1327             : 
    1328             : /********************************************************************/
    1329             : /**                                                                **/
    1330             : /**                        MULTIPLICATION                          **/
    1331             : /**                                                                **/
    1332             : /********************************************************************/
    1333             : static GEN
    1334      231029 : mul_ser_scal(GEN y, GEN x) {
    1335             :   long ly, i;
    1336             :   GEN z;
    1337      231029 :   if (isexactzero(x)) return gmul(RgX_get_0(y), x);
    1338      229734 :   if (ser_isexactzero(y))
    1339             :   {
    1340        1519 :     if (lg(y) == 2) return gcopy(y);
    1341           0 :     return scalarser(gmul(x,gel(y,2)), varn(y), valp(y));
    1342             :   }
    1343      228215 :   z = cgetg_copy(y, &ly); z[1] = y[1];
    1344      228215 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    1345      228215 :   return normalize(z);
    1346             : }
    1347             : /* (n/d) * x, x "scalar" or polynomial in the same variable as d
    1348             :  * [n/d a valid RFRAC]  */
    1349             : static GEN
    1350     9454606 : mul_rfrac_scal(GEN n, GEN d, GEN x)
    1351             : {
    1352     9454606 :   pari_sp av = avma;
    1353             :   GEN z;
    1354             : 
    1355     9454606 :   switch(typ(x))
    1356             :   {
    1357             :     case t_PADIC:
    1358           7 :       n = gmul(n, x);
    1359           7 :       d = gcvtop(d, gel(x,2), signe(gel(x,4))? precp(x): 1);
    1360           7 :       return gerepileupto(av, gdiv(n,d));
    1361             : 
    1362             :     case t_INTMOD: case t_POLMOD:
    1363          28 :       n = gmul(n, x);
    1364          28 :       d = gmul(d, gmodulo(gen_1, gel(x,1)));
    1365          28 :       return gerepileupto(av, gdiv(n,d));
    1366             :   }
    1367     9454571 :   z = gred_rfrac2(x, d);
    1368     9454571 :   n = simplify_shallow(n);
    1369     9454571 :   if (typ(z) == t_RFRAC)
    1370             :   {
    1371     7119879 :     n = gmul(gel(z,1), n);
    1372     7119879 :     d = gel(z,2);
    1373     7119879 :     if (typ(n) == t_POL && varncmp(varn(n), varn(d)) < 0)
    1374           0 :       z = RgX_Rg_div(n, d);
    1375             :     else
    1376     7119879 :       z = gred_rfrac_simple(n, d);
    1377             :   }
    1378             :   else
    1379     2334692 :     z = gmul(z, n);
    1380     9454571 :   return gerepileupto(av, z);
    1381             : }
    1382             : static GEN
    1383    56715290 : mul_scal(GEN y, GEN x, long ty)
    1384             : {
    1385    56715290 :   switch(ty)
    1386             :   {
    1387             :     case t_POL:
    1388    48448637 :       if (lg(y) == 2) return scalarpol(gmul(gen_0,x), varn(y));
    1389    47659002 :       return RgX_Rg_mul(y, x);
    1390      154967 :     case t_SER: return mul_ser_scal(y, x);
    1391     8111686 :     case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    1392             :     case t_QFI: case t_QFR:
    1393           0 :       if (typ(x) == t_INT && gequal1(x)) return gcopy(y); /* fall through */
    1394             :   }
    1395           0 :   pari_err_TYPE2("*",x,y);
    1396           0 :   return NULL; /* not reached */
    1397             : }
    1398             : 
    1399             : static GEN
    1400      160830 : mul_gen_rfrac(GEN X, GEN Y)
    1401             : {
    1402      160830 :   GEN y1 = gel(Y,1), y2 = gel(Y,2);
    1403      160830 :   long vx = gvar(X), vy = varn(y2);
    1404      166570 :   return (varncmp(vx, vy) <= 0)? mul_scal(Y, X, typ(Y)):
    1405        5740 :                                  gred_rfrac_simple(gmul(y1,X), y2);
    1406             : }
    1407             : /* (x1/x2) * (y1/y2) */
    1408             : static GEN
    1409     7744009 : mul_rfrac(GEN x1, GEN x2, GEN y1, GEN y2)
    1410             : {
    1411             :   GEN z, X, Y;
    1412     7744009 :   pari_sp av = avma;
    1413             : 
    1414     7744009 :   X = gred_rfrac2(x1, y2);
    1415     7744009 :   Y = gred_rfrac2(y1, x2);
    1416     7744009 :   if (typ(X) == t_RFRAC)
    1417             :   {
    1418     6479589 :     if (typ(Y) == t_RFRAC) {
    1419     6399664 :       x1 = gel(X,1);
    1420     6399664 :       x2 = gel(X,2);
    1421     6399664 :       y1 = gel(Y,1);
    1422     6399664 :       y2 = gel(Y,2);
    1423     6399664 :       z = gred_rfrac_simple(gmul(x1,y1), gmul(x2,y2));
    1424             :     } else
    1425       79925 :       z = mul_gen_rfrac(Y, X);
    1426             :   }
    1427     1264420 :   else if (typ(Y) == t_RFRAC)
    1428       80905 :     z = mul_gen_rfrac(X, Y);
    1429             :   else
    1430     1183515 :     z = gmul(X, Y);
    1431     7744009 :   return gerepileupto(av, z);
    1432             : }
    1433             : /* (x1/x2) /y2, x2 and y2 are t_POL in the same variable */
    1434             : static GEN
    1435      253127 : div_rfrac_pol(GEN x1, GEN x2, GEN y2)
    1436             : {
    1437      253127 :   pari_sp av = avma;
    1438      253127 :   GEN X = gred_rfrac2(x1, y2);
    1439      253127 :   if (typ(X) == t_RFRAC && varn(gel(X,2)) == varn(x2))
    1440             :   {
    1441      246575 :     x2 = RgX_mul(gel(X,2), x2);
    1442      246575 :     x1 = gel(X,1);
    1443             :   }
    1444             :   else
    1445        6552 :     x1 = X;
    1446      253127 :   return gerepileupto(av, gred_rfrac_simple(x1, x2));
    1447             : }
    1448             : 
    1449             : /* Mod(y, Y) * x,  assuming x scalar */
    1450             : static GEN
    1451     1197036 : mul_polmod_scal(GEN Y, GEN y, GEN x)
    1452             : {
    1453     1197036 :   GEN z = cgetg(3,t_POLMOD);
    1454     1197036 :   gel(z,1) = RgX_copy(Y);
    1455     1197036 :   gel(z,2) = gmul(x,y); return z;
    1456             : }
    1457             : 
    1458             : /* cf mulqq */
    1459             : static GEN
    1460     1013687 : quad_polmod_mul(GEN P, GEN x, GEN y)
    1461             : {
    1462     1013687 :   GEN T = cgetg(4, t_POL), b = gel(P,3), c = gel(P,2), p1, p2, p3, p4;
    1463     1013687 :   pari_sp tetpil, av = avma;
    1464     1013687 :   T[1] = x[1];
    1465     1013687 :   p2 = gmul(gel(x,2), gel(y,2));
    1466     1013687 :   p3 = gmul(gel(x,3), gel(y,3));
    1467     1013687 :   p1 = gmul(gneg_i(c),p3);
    1468             :   /* operands are usually small: gadd ~ gmul and Karatsuba is a waste */
    1469     1013687 :   if (typ(b) == t_INT)
    1470             :   {
    1471     1013687 :     if (signe(b))
    1472             :     {
    1473      728308 :       p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1474      728308 :       if (is_pm1(b))
    1475             :       {
    1476      727552 :         if (signe(b) > 0) p3 = gneg(p3);
    1477             :       }
    1478             :       else
    1479         756 :         p3 = gmul(negi(b), p3);
    1480             :     }
    1481             :     else
    1482             :     {
    1483      285379 :       p3 = gmul(gel(x,2),gel(y,3));
    1484      285379 :       p4 = gmul(gel(x,3),gel(y,2));
    1485             :     }
    1486             :   }
    1487             :   else
    1488             :   {
    1489           0 :     p4 = gadd(gmul(gel(x,2), gel(y,3)), gmul(gel(x,3), gel(y,2)));
    1490           0 :     p3 = gmul(gneg_i(b), p3);
    1491             :   }
    1492     1013687 :   tetpil = avma;
    1493     1013687 :   gel(T,2) = gadd(p2, p1);
    1494     1013687 :   gel(T,3) = gadd(p4, p3);
    1495     1013687 :   gerepilecoeffssp(av,tetpil,T+2,2);
    1496     1013687 :   return normalizepol_lg(T,4);
    1497             : }
    1498             : /* Mod(x,T) * Mod(y,T) */
    1499             : static GEN
    1500     1536367 : mul_polmod_same(GEN T, GEN x, GEN y)
    1501             : {
    1502     1536367 :   GEN z = cgetg(3,t_POLMOD), a;
    1503     1536367 :   long v = varn(T), lx = lg(x), ly = lg(y);
    1504     1536367 :   gel(z,1) = RgX_copy(T);
    1505             :   /* x * y mod T optimised */
    1506     1536367 :   if (typ(x) != t_POL || varn(x) != v || lx <= 3
    1507     1307777 :    || typ(y) != t_POL || varn(y) != v || ly <= 3)
    1508      310217 :     a = gmul(x, y);
    1509             :   else
    1510             :   {
    1511     1226150 :     if (lg(T) == 5 && isint1(gel(T,4))) /* quadratic fields */
    1512     1012322 :       a = quad_polmod_mul(T, x, y);
    1513             :     else
    1514             :     {
    1515      213828 :       pari_sp av = avma;
    1516      213828 :       GEN p = NULL;
    1517      213828 :       if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    1518             :       {
    1519       37121 :         T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    1520       37121 :         if (lgefint(p) == 3)
    1521             :         {
    1522       32284 :           ulong pp = p[2];
    1523       32284 :           a = Flxq_mul(ZX_to_Flx(x, pp), ZX_to_Flx(y, pp), ZX_to_Flx(T, pp), pp);
    1524       32284 :           a = Flx_to_ZX(a);
    1525             :         }
    1526             :         else
    1527        4837 :           a = FpXQ_mul(x, y, T, p);
    1528       37121 :         a = FpX_to_mod(a, p);
    1529             :       }
    1530             :       else
    1531      176707 :         a = RgXQ_mul(x, y, gel(z,1));
    1532      213828 :       a = gerepileupto(av, a);
    1533             :     }
    1534             :   }
    1535     1536367 :   gel(z,2) = a; return z;
    1536             : }
    1537             : static GEN
    1538       24107 : sqr_polmod(GEN T, GEN x)
    1539             : {
    1540       24107 :   GEN a, z = cgetg(3,t_POLMOD);
    1541       24107 :   gel(z,1) = RgX_copy(T);
    1542       24107 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    1543        3710 :     a = gsqr(x);
    1544             :   else
    1545             :   {
    1546       20397 :     pari_sp av = avma;
    1547       20397 :     GEN p = NULL;
    1548       20397 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && p)
    1549             :     {
    1550          14 :       T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p);
    1551          14 :       if (lgefint(p) == 3)
    1552             :       {
    1553           7 :         ulong pp = p[2];
    1554           7 :         a = Flxq_sqr(ZX_to_Flx(x, pp), ZX_to_Flx(T, pp), pp);
    1555           7 :         a = Flx_to_ZX(a);
    1556             :       }
    1557             :       else
    1558           7 :         a = FpXQ_sqr(x, T, p);
    1559          14 :       a = FpX_to_mod(a, p);
    1560             :     }
    1561             :     else
    1562       20383 :       a = RgXQ_sqr(x, gel(z,1));
    1563       20397 :     a = gerepileupto(av, a);
    1564             :   }
    1565       24107 :   gel(z,2) = a; return z;
    1566             : }
    1567             : /* Mod(x,X) * Mod(y,Y) */
    1568             : static GEN
    1569        1386 : mul_polmod(GEN X, GEN Y, GEN x, GEN y)
    1570             : {
    1571        1386 :   long T[3] = { evaltyp(t_POLMOD) | _evallg(3),0,0 };
    1572        1386 :   long vx = varn(X), vy = varn(Y);
    1573        1386 :   GEN z = cgetg(3,t_POLMOD);
    1574             : 
    1575        1386 :   if (vx==vy) {
    1576             :     pari_sp av;
    1577          14 :     gel(z,1) = RgX_gcd(X,Y); av = avma;
    1578          14 :     warn_coercion(X,Y,gel(z,1));
    1579          14 :     gel(z,2) = gerepileupto(av, gmod(gmul(x, y), gel(z,1)));
    1580          14 :     return z;
    1581             :   }
    1582        1372 :   if (varncmp(vx, vy) < 0)
    1583          21 :   { gel(z,1) = RgX_copy(X); gel(T,1) = Y; gel(T,2) = y; y = T; }
    1584             :   else
    1585        1351 :   { gel(z,1) = RgX_copy(Y); gel(T,1) = X; gel(T,2) = x; x = T; }
    1586        1372 :   gel(z,2) = gmul(x, y); return z;
    1587             : }
    1588             : 
    1589             : #if 0 /* used by 3M only */
    1590             : /* set z = x+y and return 1 if x,y have the same sign
    1591             :  * set z = x-y and return 0 otherwise */
    1592             : static int
    1593             : did_add(GEN x, GEN y, GEN *z)
    1594             : {
    1595             :   long tx = typ(x), ty = typ(y);
    1596             :   if (tx == ty) switch(tx)
    1597             :   {
    1598             :     case t_INT: *z = addii(x,y); return 1;
    1599             :     case t_FRAC: *z = addsub_frac(x,y,addii); return 1;
    1600             :     case t_REAL:
    1601             :       if (signe(x) == -signe(y))
    1602             :       { *z = subrr(x,y); return 0; }
    1603             :       else
    1604             :       { *z = addrr(x,y); return 1; }
    1605             :   }
    1606             :   if (tx == t_REAL) switch(ty)
    1607             :   {
    1608             :     case t_INT:
    1609             :       if (signe(x) == -signe(y))
    1610             :       { *z = subri(x,y); return 0; }
    1611             :       else
    1612             :       { *z = addri(x,y); return 1; }
    1613             :     case t_FRAC:
    1614             :       if (signe(x) == -signe(gel(y,1)))
    1615             :       { *z = gsub(x,y); return 0; }
    1616             :       else
    1617             :       { *z = gadd(x,y); return 1; }
    1618             :   }
    1619             :   else if (ty == t_REAL) switch(tx)
    1620             :   {
    1621             :     case t_INT:
    1622             :       if (signe(x) == -signe(y))
    1623             :       { *z = subir(x,y); return 0; }
    1624             :       else
    1625             :       { *z = addir(x,y); return 1; }
    1626             :     case t_FRAC:
    1627             :       if (signe(gel(x,1)) == -signe(y))
    1628             :       { *z = gsub(x,y); return 0; }
    1629             :       else
    1630             :       { *z = gadd(x,y); return 1; }
    1631             :   }
    1632             :   *z = gadd(x,y); return 1;
    1633             : }
    1634             : #endif
    1635             : /* x * I * y, x t_COMPLEX with non-intzero real part, y non-intzero "scalar" */
    1636             : static GEN
    1637      944471 : mulcIR(GEN x, GEN y)
    1638             : {
    1639      944471 :   GEN z = cgetg(3,t_COMPLEX);
    1640      944471 :   pari_sp av = avma;
    1641      944471 :   gel(z,1) = gerepileupto(av, gneg(gmul(y,gel(x,2))));
    1642      944471 :   gel(z,2) = gmul(y, gel(x,1));
    1643      944471 :   return z;
    1644             : 
    1645             : }
    1646             : /* x,y COMPLEX */
    1647             : static GEN
    1648    73666156 : mulcc(GEN x, GEN y)
    1649             : {
    1650    73666156 :   GEN xr = gel(x,1), xi = gel(x,2);
    1651    73666156 :   GEN yr = gel(y,1), yi = gel(y,2);
    1652             :   GEN p1, p2, p3, p4, z;
    1653             :   pari_sp tetpil, av;
    1654             : 
    1655    73666156 :   if (isintzero(xr))
    1656             :   {
    1657      921055 :     if (isintzero(yr)) {
    1658      133187 :       av = avma;
    1659      133187 :       return gerepileupto(av, gneg(gmul(xi,yi)));
    1660             :     }
    1661      787868 :     return mulcIR(y, xi);
    1662             :   }
    1663    72745101 :   if (isintzero(yr)) return mulcIR(x, yi);
    1664             : 
    1665    72588498 :   z = cgetg(3,t_COMPLEX); av = avma;
    1666             : #if 0
    1667             :   /* 3M method avoiding catastrophic cancellation, BUT loses accuracy due to
    1668             :    * e.g. xr + xi if exponents differ */
    1669             :   if (did_add(xr, xi, &p3))
    1670             :   {
    1671             :     if (did_add(yr, yi, &p4)) {
    1672             :     /* R = xr*yr - xi*yi
    1673             :      * I = (xr+xi)(yr+yi) - xr*yr - xi*yi */
    1674             :       p1 = gmul(xr,yr);
    1675             :       p2 = gmul(xi,yi); p2 = gneg(p2);
    1676             :       p3 = gmul(p3, p4);
    1677             :       p4 = gsub(p2, p1);
    1678             :     } else {
    1679             :     /* R = (xr + xi) * (yr - yi) + (xr * yi - xi * yr)
    1680             :      * I = xr*yi + xi*yr */
    1681             :       p1 = gmul(p3,p4);
    1682             :       p3 = gmul(xr,yi);
    1683             :       p4 = gmul(xi,yr);
    1684             :       p2 = gsub(p3, p4);
    1685             :     }
    1686             :   } else {
    1687             :     if (did_add(yr, yi, &p4)) {
    1688             :      /* R = (xr - xi) * (yr + yi) + (xi * yr - xr * yi)
    1689             :       * I = xr*yi +xi*yr */
    1690             :       p1 = gmul(p3,p4);
    1691             :       p3 = gmul(xr,yi);
    1692             :       p4 = gmul(xi,yr);
    1693             :       p2 = gsub(p4, p3);
    1694             :     } else {
    1695             :     /* R = xr*yr - xi*yi
    1696             :      * I = -(xr-xi)(yr-yi) + xr*yr + xi*yi */
    1697             :       p3 = gneg( gmul(p3, p4) );
    1698             :       p1 = gmul(xr,yr);
    1699             :       p2 = gmul(xi,yi);
    1700             :       p4 = gadd(p1, p2);
    1701             : 
    1702             :       p2 = gneg(p2);
    1703             :     }
    1704             :   }
    1705             :   tetpil = avma;
    1706             :   gel(z,1) = gadd(p1,p2);
    1707             :   gel(z,2) = gadd(p3,p4);
    1708             : #else
    1709    72588498 :   if (typ(xr)==t_INT && typ(yr)==t_INT && typ(xi)==t_INT && typ(yi)==t_INT)
    1710             :   { /* 3M formula */
    1711         168 :     p3 = addii(xr,xi);
    1712         168 :     p4 = addii(yr,yi);
    1713         168 :     p1 = mulii(xr,yr);
    1714         168 :     p2 = mulii(xi,yi);
    1715         168 :     p3 = mulii(p3,p4);
    1716         168 :     p4 = addii(p2,p1);
    1717         168 :     tetpil = avma;
    1718         168 :     gel(z,1) = subii(p1,p2);
    1719         168 :     gel(z,2) = subii(p3,p4);
    1720         287 :     if (!signe(gel(z,2)))
    1721          49 :       return gerepileuptoint((pari_sp)(z+3), gel(z,1));
    1722             :   }
    1723             :   else
    1724             :   { /* naive 4M formula: avoid all loss of accuracy */
    1725    72588330 :     p1 = gmul(xr,yr);
    1726    72588330 :     p2 = gmul(xi,yi);
    1727    72588330 :     p3 = gmul(xr,yi);
    1728    72588330 :     p4 = gmul(xi,yr);
    1729    72588330 :     tetpil = avma;
    1730    72588330 :     gel(z,1) = gsub(p1,p2);
    1731    72588330 :     gel(z,2) = gadd(p3,p4);
    1732    72588330 :     if (isintzero(gel(z,2)))
    1733             :     {
    1734        9389 :       cgiv(gel(z,2));
    1735        9389 :       return gerepileupto((pari_sp)(z+3), gel(z,1));
    1736             :     }
    1737             :   }
    1738             : #endif
    1739             : 
    1740    72579060 :   gerepilecoeffssp(av,tetpil, z+1,2); return z;
    1741             : }
    1742             : /* x,y PADIC */
    1743             : static GEN
    1744     1756606 : mulpp(GEN x, GEN y) {
    1745     1756606 :   long l = valp(x) + valp(y);
    1746             :   pari_sp av;
    1747             :   GEN z, t;
    1748     1756606 :   if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("*",x,y);
    1749     1756606 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), l);
    1750     1596425 :   if (!signe(gel(y,4))) return zeropadic(gel(x,2), l);
    1751             : 
    1752     1532396 :   t = (precp(x) > precp(y))? y: x;
    1753     1532396 :   z = cgetp(t); setvalp(z,l); av = avma;
    1754     1532395 :   affii(remii(mulii(gel(x,4),gel(y,4)), gel(t,3)), gel(z,4));
    1755     1532397 :   avma = av; return z;
    1756             : }
    1757             : /* x,y QUAD */
    1758             : static GEN
    1759        1295 : mulqq(GEN x, GEN y) {
    1760        1295 :   GEN z = cgetg(4,t_QUAD);
    1761        1295 :   GEN p1, p2, p3, p4, P = gel(x,1), b = gel(P,3), c = gel(P,2);
    1762             :   pari_sp av, tetpil;
    1763        1295 :   if (!ZX_equal(P, gel(y,1))) pari_err_OP("*",x,y);
    1764             : 
    1765        1295 :   gel(z,1) = ZX_copy(P); av = avma;
    1766        1295 :   p2 = gmul(gel(x,2),gel(y,2));
    1767        1295 :   p3 = gmul(gel(x,3),gel(y,3));
    1768        1295 :   p1 = gmul(gneg_i(c),p3);
    1769             : 
    1770        1295 :   if (signe(b))
    1771        1183 :     p4 = gadd(gmul(gel(x,2),gel(y,3)), gmul(gel(x,3),gel(y,2)));
    1772             :   else
    1773             :   {
    1774         112 :     p3 = gmul(gel(x,2),gel(y,3));
    1775         112 :     p4 = gmul(gel(x,3),gel(y,2));
    1776             :   }
    1777        1295 :   tetpil = avma;
    1778        1295 :   gel(z,2) = gadd(p2,p1);
    1779        1295 :   gel(z,3) = gadd(p4,p3);
    1780        1295 :   gerepilecoeffssp(av,tetpil,z+2,2); return z;
    1781             : }
    1782             : 
    1783             : GEN
    1784     2147568 : mulcxI(GEN x)
    1785             : {
    1786             :   GEN z;
    1787     2147568 :   switch(typ(x))
    1788             :   {
    1789             :     case t_INT: case t_REAL: case t_FRAC:
    1790      216818 :       return mkcomplex(gen_0, x);
    1791             :     case t_COMPLEX:
    1792     1930694 :       if (isintzero(gel(x,1))) return gneg(gel(x,2));
    1793     1928573 :       z = cgetg(3,t_COMPLEX);
    1794     1928573 :       gel(z,1) = gneg(gel(x,2));
    1795     1928573 :       gel(z,2) = gel(x,1); return z;
    1796             :     default:
    1797          56 :       return gmul(gen_I(), x);
    1798             :   }
    1799             : }
    1800             : GEN
    1801       47845 : mulcxmI(GEN x)
    1802             : {
    1803             :   GEN z;
    1804       47845 :   switch(typ(x))
    1805             :   {
    1806             :     case t_INT: case t_REAL: case t_FRAC:
    1807        2051 :       return mkcomplex(gen_0, gneg(x));
    1808             :     case t_COMPLEX:
    1809       45738 :       if (isintzero(gel(x,1))) return gel(x,2);
    1810       44744 :       z = cgetg(3,t_COMPLEX);
    1811       44744 :       gel(z,1) = gel(x,2);
    1812       44744 :       gel(z,2) = gneg(gel(x,1)); return z;
    1813             :     default:
    1814          56 :       return gmul(mkcomplex(gen_0, gen_m1), x);
    1815             :   }
    1816             : }
    1817             : 
    1818             : /* fill in coefficients of t_SER z from coeffs of t_POL y */
    1819             : static GEN
    1820     1661419 : fill_ser(GEN z, GEN y)
    1821             : {
    1822     1661419 :   long i, lx = lg(z), ly = lg(y);
    1823     1661419 :   if (ly >= lx) {
    1824     1656435 :     for (i = 2; i < lx; i++) gel(z,i) = gel(y,i);
    1825             :   } else {
    1826        4984 :     for (i = 2; i < ly; i++) gel(z,i) = gel(y,i);
    1827        4984 :     for (     ; i < lx; i++) gel(z,i) = gen_0;
    1828             :   }
    1829     1661419 :   return normalize(z);
    1830             : }
    1831             : 
    1832             : GEN
    1833  3309151266 : gmul(GEN x, GEN y)
    1834             : {
    1835             :   long tx, ty, lx, ly, vx, vy, i, l;
    1836             :   pari_sp av, tetpil;
    1837             :   GEN z, p1, p2;
    1838             : 
    1839  3309151266 :   if (x == y) return gsqr(x);
    1840  2473307517 :   tx = typ(x); ty = typ(y);
    1841  2473307517 :   if (tx == ty) switch(tx)
    1842             :   {
    1843  1409407149 :     case t_INT: return mulii(x,y);
    1844   526041924 :     case t_REAL: return mulrr(x,y);
    1845     1478292 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    1846     1478292 :       z = cgetg(3,t_INTMOD);
    1847     1478289 :       if (X==Y || equalii(X,Y))
    1848     1478276 :         return mul_intmod_same(z, X, gel(x,2), gel(y,2));
    1849          14 :       gel(z,1) = gcdii(X,Y);
    1850          14 :       warn_coercion(X,Y,gel(z,1));
    1851          14 :       av = avma; p1 = mulii(gel(x,2),gel(y,2));
    1852          14 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    1853             :     }
    1854             :     case t_FRAC:
    1855             :     {
    1856    12580430 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    1857    12580430 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    1858    12580430 :       z=cgetg(3,t_FRAC);
    1859    12580430 :       p1 = gcdii(x1, y2);
    1860    12580430 :       if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y2 = diviiexact(y2,p1); }
    1861    12580430 :       p1 = gcdii(x2, y1);
    1862    12580430 :       if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y1 = diviiexact(y1,p1); }
    1863    12580430 :       tetpil = avma;
    1864    12580430 :       gel(z,2) = mulii(x2,y2);
    1865    12580430 :       gel(z,1) = mulii(x1,y1);
    1866    12580430 :       fix_frac_if_int_GC(z,tetpil); return z;
    1867             :     }
    1868    72319276 :     case t_COMPLEX: return mulcc(x, y);
    1869     1590859 :     case t_PADIC: return mulpp(x, y);
    1870         973 :     case t_QUAD: return mulqq(x, y);
    1871    10652401 :     case t_FFELT: return FF_mul(x, y);
    1872             :     case t_POLMOD:
    1873     1343209 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    1874     1341823 :         return mul_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    1875        1386 :       return mul_polmod(gel(x,1), gel(y,1), gel(x,2), gel(y,2));
    1876             :     case t_POL:
    1877    33282203 :       vx = varn(x);
    1878    33282203 :       vy = varn(y);
    1879    33282203 :       if (vx != vy) {
    1880     4696467 :         if (varncmp(vx, vy) < 0) return RgX_Rg_mul(x, y);
    1881     2005647 :         else                     return RgX_Rg_mul(y, x);
    1882             :       }
    1883    28585736 :       if (RgX_is_ZX(x) && RgX_is_ZX(y)) return ZX_mul(x,y);
    1884     5182672 :       return RgX_mul(x, y);
    1885             : 
    1886             :     case t_SER: {
    1887     1661286 :       GEN p = NULL;
    1888     1661286 :       vx = varn(x);
    1889     1661286 :       vy = varn(y);
    1890     1661286 :       if (vx != vy) {
    1891           0 :         if (varncmp(vx, vy) < 0) return mul_ser_scal(x, y);
    1892           0 :         else                     return mul_ser_scal(y, x);
    1893             :       }
    1894     1661286 :       lx = lg(x);
    1895     1661286 :       ly = lg(y); if (lx > ly) { lx = ly; swap(x, y); }
    1896     1661286 :       if (lx == 2) return zeroser(vx, valp(x)+valp(y));
    1897     1661272 :       z = cgetg(lx,t_SER);
    1898     1661272 :       z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
    1899     1661272 :       x = ser2pol_i(x, lx);
    1900     1661272 :       y = ser2pol_i(y, lx);
    1901     1661272 :       if (RgX_is_FpX(x,&p) && RgX_is_FpX(y,&p))
    1902             :       {
    1903       30760 :         if (!p) y = ZX_mul(x,y);
    1904             :         else
    1905             :         {
    1906          21 :           x = RgX_to_FpX(x, p);
    1907          21 :           y = RgX_to_FpX(y, p);
    1908          21 :           y = FpX_to_mod(ZX_mul(x,y), p);
    1909             :         }
    1910             :       }
    1911             :       else
    1912     1645892 :         y = RgXn_mul(x, y, lx-2);
    1913     1661272 :       z = fill_ser(z, y);
    1914     1661272 :       return gerepilecopy((pari_sp)(z + lx), z);
    1915             :     }
    1916    13130685 :     case t_QFI: return qficomp(x,y);
    1917           7 :     case t_QFR: return qfrcomp(x,y);
    1918     6557680 :     case t_RFRAC: return mul_rfrac(gel(x,1),gel(x,2), gel(y,1),gel(y,2));
    1919     1522788 :     case t_MAT: return RgM_mul(x, y);
    1920             : 
    1921             :     case t_VECSMALL: /* multiply as permutation. cf perm_mul */
    1922      131621 :       z = cgetg_copy(x, &l);
    1923      131621 :       if (l != lg(y)) break;
    1924     2310854 :       for (i=1; i<l; i++)
    1925             :       {
    1926     2179233 :         long yi = y[i];
    1927     2179233 :         if (yi < 1 || yi >= l) pari_err_TYPE2("*",x,y);
    1928     2179233 :         z[i] = x[yi];
    1929             :       }
    1930      131621 :       return z;
    1931             : 
    1932             : 
    1933             :     default:
    1934           0 :       pari_err_TYPE2("*",x,y);
    1935             :   }
    1936             :   /* tx != ty */
    1937   381606734 :   if (is_const_t(ty) && is_const_t(tx))  {
    1938   314492169 :     if (tx > ty) { swap(x,y); lswap(tx,ty); }
    1939   314492169 :     switch(tx) {
    1940             :     case t_INT:
    1941   286328860 :       switch(ty)
    1942             :       {
    1943   169135134 :         case t_REAL: return signe(x)? mulir(x,y): gen_0;
    1944             :         case t_INTMOD:
    1945     1584118 :           z = cgetg(3, t_INTMOD);
    1946     1584109 :           return mul_intmod_same(z, gel(y,1), gel(y,2), modii(x, gel(y,1)));
    1947             :         case t_FRAC:
    1948    58279332 :           if (!signe(x)) return gen_0;
    1949    28477617 :           z=cgetg(3,t_FRAC);
    1950    28477617 :           p1 = gcdii(x,gel(y,2));
    1951    28477617 :           if (is_pm1(p1))
    1952             :           {
    1953    15261251 :             avma = (pari_sp)z;
    1954    15261251 :             gel(z,2) = icopy(gel(y,2));
    1955    15261251 :             gel(z,1) = mulii(gel(y,1), x);
    1956             :           }
    1957             :           else
    1958             :           {
    1959    13216366 :             x = diviiexact(x,p1); tetpil = avma;
    1960    13216366 :             gel(z,2) = diviiexact(gel(y,2), p1);
    1961    13216366 :             gel(z,1) = mulii(gel(y,1), x);
    1962    13216366 :             fix_frac_if_int_GC(z,tetpil);
    1963             :           }
    1964    28477617 :           return z;
    1965    55362892 :         case t_COMPLEX: return signe(x)? mulRc(x, y): gen_0;
    1966      397909 :         case t_PADIC: return signe(x)? mulTp(x, y): gen_0;
    1967        1218 :         case t_QUAD: return mulRq(x,y);
    1968     1568260 :         case t_FFELT: return FF_Z_mul(y,x);
    1969             :       }
    1970             : 
    1971             :     case t_REAL:
    1972    28057338 :       switch(ty)
    1973             :       {
    1974     1583329 :         case t_FRAC: return mulrfrac(x, y);
    1975    26473988 :         case t_COMPLEX: return mulRc(x, y);
    1976           7 :         case t_QUAD: return mulqf(y, x, lg(x));
    1977          14 :         default: pari_err_TYPE2("*",x,y);
    1978             :       }
    1979             : 
    1980             :     case t_INTMOD:
    1981        5443 :       switch(ty)
    1982             :       {
    1983         588 :         case t_FRAC: { GEN X = gel(x,1);
    1984         588 :           z = cgetg(3, t_INTMOD); p1 = Fp_mul(gel(y,1), gel(x,2), X);
    1985         588 :           return div_intmod_same(z, X, p1, remii(gel(y,2), X));
    1986             :         }
    1987          35 :         case t_COMPLEX: return mulRc_direct(x,y);
    1988           0 :         case t_PADIC: { GEN X = gel(x,1);
    1989           0 :           z = cgetg(3, t_INTMOD);
    1990           0 :           return mul_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    1991             :         }
    1992          42 :         case t_QUAD: return mulRq(x, y);
    1993             :         case t_FFELT:
    1994        4778 :           if (!equalii(gel(x,1),FF_p_i(y)))
    1995           0 :             pari_err_OP("*",x,y);
    1996        4778 :           return FF_Z_mul(y,gel(x,2));
    1997             :       }
    1998             : 
    1999             :     case t_FRAC:
    2000      100465 :       switch(ty)
    2001             :       {
    2002       69182 :         case t_COMPLEX: return mulRc(x, y);
    2003       28546 :         case t_PADIC: return signe(gel(x,1))? mulTp(x, y): gen_0;
    2004         308 :         case t_QUAD:  return mulRq(x, y);
    2005        2429 :         case t_FFELT: return FF_Z_Z_muldiv(y, gel(x,1),gel(x,2));
    2006             :       }
    2007             : 
    2008             :     case t_FFELT:
    2009          35 :       pari_err_TYPE2("*",x,y);
    2010             : 
    2011             :     case t_COMPLEX:
    2012          14 :       switch(ty)
    2013             :       {
    2014             :         case t_PADIC:
    2015           7 :           return Zp_nosquare_m1(gel(y,2))? mulRc(y, x): mulTp(x, y);
    2016             :         case t_QUAD:
    2017           7 :           lx = precision(x); if (!lx) pari_err_OP("*",x,y);
    2018           7 :           return mulqf(y, x, lx);
    2019             :       }
    2020             : 
    2021             :     case t_PADIC: /* ty == t_QUAD */
    2022          21 :       return (kro_quad(y,gel(x,2))== -1)? mulRq(x, y): mulTp(y, x);
    2023             :     }
    2024             :   }
    2025             : 
    2026    67187122 :   if (is_matvec_t(ty))
    2027             :   {
    2028     7689284 :     if (!is_matvec_t(tx))
    2029             :     {
    2030     7424466 :       if (is_noncalc_t(tx)) pari_err_TYPE2( "*",x,y); /* necessary if ly = 1 */
    2031     7424466 :       z = cgetg_copy(y, &ly);
    2032     7424466 :       for (i=1; i<ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2033     7424466 :       return z;
    2034             :     }
    2035      264818 :     switch(tx)
    2036             :     {
    2037             :       case t_VEC:
    2038      176107 :         switch(ty) {
    2039      123022 :           case t_COL: return RgV_RgC_mul(x,y);
    2040       53085 :           case t_MAT: return RgV_RgM_mul(x,y);
    2041             :         }
    2042           0 :         break;
    2043             :       case t_COL:
    2044           7 :         switch(ty) {
    2045           7 :           case t_VEC: return RgC_RgV_mul(x,y);
    2046           0 :           case t_MAT: return RgC_RgM_mul(x,y);
    2047             :         }
    2048           0 :         break;
    2049             :       case t_MAT:
    2050       88704 :         switch(ty) {
    2051           0 :           case t_VEC: return RgM_RgV_mul(x,y);
    2052       88704 :           case t_COL: return RgM_RgC_mul(x,y);
    2053             :         }
    2054             :     }
    2055             :   }
    2056    59497852 :   if (is_matvec_t(tx))
    2057             :   {
    2058      124640 :     if (is_noncalc_t(ty)) pari_err_TYPE2( "*",x,y); /* necessary if lx = 1 */
    2059      124641 :     z = cgetg_copy(x, &lx);
    2060      124640 :     for (i=1; i<lx; i++) gel(z,i) = gmul(y,gel(x,i));
    2061      124640 :     return z;
    2062             :   }
    2063    59373212 :   if (tx > ty) { swap(x,y); lswap(tx,ty); }
    2064             :   /* tx < ty, !ismatvec(x and y) */
    2065             : 
    2066    59373212 :   if (ty == t_POLMOD) /* is_const_t(tx) in this case */
    2067     1189798 :     return mul_polmod_scal(gel(y,1), gel(y,2), x);
    2068    58183414 :   if (is_scalar_t(tx)) {
    2069    55347810 :     if (tx == t_POLMOD) {
    2070      235445 :       vx = varn(gel(x,1));
    2071      235445 :       vy = gvar(y);
    2072      235445 :       if (vx != vy) {
    2073       41034 :         if (varncmp(vx,vy) > 0) return mul_scal(y, x, ty);
    2074        7238 :         return mul_polmod_scal(gel(x,1), gel(x,2), y);
    2075             :       }
    2076             :       /* error if ty == t_SER */
    2077      194411 :       av = avma; y = gmod(y, gel(x,1));
    2078      194404 :       return gerepileupto(av, mul_polmod_same(gel(x,1), gel(x,2), y));
    2079             :     }
    2080    55112365 :     return mul_scal(y, x, ty);
    2081             :   }
    2082             : 
    2083             :   /* x and y are not scalars, nor matvec */
    2084     2835604 :   vx = gvar(x);
    2085     2835604 :   vy = gvar(y);
    2086     2835604 :   if (vx != vy) /* x or y is treated as a scalar */
    2087     4184675 :     return (varncmp(vx, vy) < 0)? mul_scal(x, y, tx)
    2088     2770636 :                                 : mul_scal(y, x, ty);
    2089             :   /* vx = vy */
    2090     1421565 :   switch(tx)
    2091             :   {
    2092             :     case t_POL:
    2093     1421565 :       switch (ty)
    2094             :       {
    2095             :         case t_SER:
    2096             :         {
    2097             :           long vn;
    2098       78673 :           if (lg(x) == 2) return zeropol(vx);
    2099       78673 :           if (lg(y) == 2) return zeroser(vx, valp(y)+RgX_val(x));
    2100       78463 :           av = avma;
    2101       78463 :           vn = RgX_valrem(x, &x);
    2102             :           /* take advantage of x = t^n ! */
    2103       78463 :           if (degpol(x)) {
    2104        2401 :             p1 = RgX_to_ser(x,lg(y));
    2105        2401 :             if (vn) settyp(x, t_VECSMALL); /* *new* x left on stack */
    2106        2401 :             p2 = gmul(p1,y);
    2107        2401 :             settyp(p1, t_VECSMALL); /* p1 left on stack */
    2108             :           } else {
    2109       76062 :             avma = av;
    2110       76062 :             p2 = mul_ser_scal(y, gel(x,2));
    2111             :           }
    2112       78463 :           setvalp(p2, valp(p2) + vn);
    2113       78463 :           return p2;
    2114             :         }
    2115             : 
    2116     1342892 :         case t_RFRAC: return mul_rfrac_scal(gel(y,1),gel(y,2), x);
    2117             :       }
    2118           0 :       break;
    2119             : 
    2120             :     case t_SER:
    2121           0 :       switch (ty)
    2122             :       {
    2123             :         case t_RFRAC:
    2124           0 :           av = avma;
    2125           0 :           return gerepileupto(av, gdiv(gmul(gel(y,1),x), gel(y,2)));
    2126             :       }
    2127           0 :       break;
    2128             :   }
    2129           0 :   pari_err_TYPE2("*",x,y);
    2130           0 :   return NULL; /* not reached */
    2131             : }
    2132             : 
    2133             : int
    2134       30359 : ff_poltype(GEN *x, GEN *p, GEN *pol)
    2135             : {
    2136       30359 :   GEN Q, P = *x;
    2137       30359 :   if (!signe(P)) return 0;
    2138       30212 :   if (!RgX_is_FpXQX(P,pol,p) || !*p || !*pol || typ(*pol)!=t_POL)
    2139       30023 :     return 0;
    2140         189 :   Q = RgX_to_FpX(*pol, *p);
    2141         189 :   P = RgX_to_FpXQX(P, Q, *p);
    2142         189 :   *x = ZXX_to_Kronecker(P, degpol(Q));
    2143         189 :   *pol = Q;
    2144         189 :   return 1;
    2145             : }
    2146             : 
    2147             : /* return a non-normalized result */
    2148             : GEN
    2149       10122 : sqr_ser_part(GEN x, long l1, long l2)
    2150             : {
    2151             :   long i, j, l;
    2152             :   pari_sp av;
    2153             :   GEN Z, z, p1, p2;
    2154             :   long mi;
    2155       10122 :   if (l2 < l1) return zeroser(varn(x), 2*valp(x));
    2156       10115 :   p2 = cgetg(l2+2, t_VECSMALL)+1; /* left on stack on exit */
    2157       10115 :   Z = cgetg(l2-l1+3, t_SER);
    2158       10115 :   Z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x));
    2159       10115 :   z = Z + 2-l1;
    2160       10115 :   x += 2; mi = 0;
    2161       43337 :   for (i=0; i<l1; i++)
    2162             :   {
    2163       33222 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2164             :   }
    2165             : 
    2166       69398 :   for (i=l1; i<=l2; i++)
    2167             :   {
    2168       59283 :     p2[i] = !isrationalzero(gel(x,i)); if (p2[i]) mi = i;
    2169       59283 :     p1=gen_0; av=avma; l=((i+1)>>1) - 1;
    2170      599284 :     for (j=i-mi; j<=minss(l,mi); j++)
    2171      540001 :       if (p2[j] && p2[i-j]) p1 = gadd(p1, gmul(gel(x,j),gel(x,i-j)));
    2172       59283 :     p1 = gshift(p1,1);
    2173       59283 :     if ((i&1) == 0 && p2[i>>1])
    2174       22197 :       p1 = gadd(p1, gsqr(gel(x,i>>1)));
    2175       59283 :     gel(z,i) = gerepileupto(av,p1);
    2176             :   }
    2177       10115 :   return Z;
    2178             : }
    2179             : 
    2180             : GEN
    2181   908165772 : gsqr(GEN x)
    2182             : {
    2183             :   long i, lx;
    2184             :   pari_sp av, tetpil;
    2185             :   GEN z, p1, p2, p3, p4;
    2186             : 
    2187   908165772 :   switch(typ(x))
    2188             :   {
    2189   849199337 :     case t_INT: return sqri(x);
    2190    19129548 :     case t_REAL: return sqrr(x);
    2191      144232 :     case t_INTMOD: { GEN X = gel(x,1);
    2192      144232 :       z = cgetg(3,t_INTMOD);
    2193      144225 :       gel(z,2) = gerepileuptoint((pari_sp)z, remii(sqri(gel(x,2)), X));
    2194      144223 :       gel(z,1) = icopy(X); return z;
    2195             :     }
    2196     1417951 :     case t_FRAC: return sqrfrac(x);
    2197             : 
    2198             :     case t_COMPLEX:
    2199     2322568 :       if (isintzero(gel(x,1))) {
    2200       19432 :         av = avma;
    2201       19432 :         return gerepileupto(av, gneg(gsqr(gel(x,2))));
    2202             :       }
    2203     2303136 :       z = cgetg(3,t_COMPLEX); av = avma;
    2204     2303136 :       p1 = gadd(gel(x,1),gel(x,2));
    2205     2303136 :       p2 = gsub(gel(x,1), gel(x,2));
    2206     2303136 :       p3 = gmul(gel(x,1),gel(x,2));
    2207     2303136 :       tetpil = avma;
    2208     2303136 :       gel(z,1) = gmul(p1,p2);
    2209     2303136 :       gel(z,2) = gshift(p3,1); gerepilecoeffssp(av,tetpil,z+1,2); return z;
    2210             : 
    2211             :     case t_PADIC:
    2212       26810 :       z = cgetg(5,t_PADIC);
    2213       26810 :       i = (absequaliu(gel(x,2), 2) && signe(gel(x,4)))? 1: 0;
    2214       26810 :       if (i && precp(x) == 1) i = 2; /* (1 + O(2))^2 = 1 + O(2^3) */
    2215       26810 :       z[1] = evalprecp(precp(x)+i) | evalvalp(valp(x)*2);
    2216       26810 :       gel(z,2) = icopy(gel(x,2));
    2217       26810 :       gel(z,3) = shifti(gel(x,3), i); av = avma;
    2218       26810 :       gel(z,4) = gerepileuptoint(av, remii(sqri(gel(x,4)), gel(z,3)));
    2219       26810 :       return z;
    2220             : 
    2221          28 :     case t_QUAD: z = cgetg(4,t_QUAD);
    2222          28 :       p1 = gel(x,1);
    2223          28 :       gel(z,1) = ZX_copy(p1); av = avma;
    2224          28 :       p2 = gsqr(gel(x,2));
    2225          28 :       p3 = gsqr(gel(x,3));
    2226          28 :       p4 = gmul(gneg_i(gel(p1,2)),p3);
    2227             : 
    2228          28 :       if (gequal0(gel(p1,3)))
    2229             :       {
    2230           7 :         tetpil = avma;
    2231           7 :         gel(z,2) = gerepile(av,tetpil,gadd(p4,p2));
    2232           7 :         av = avma;
    2233           7 :         p2 = gmul(gel(x,2),gel(x,3)); tetpil = avma;
    2234           7 :         gel(z,3) = gerepile(av,tetpil,gmul2n(p2,1)); return z;
    2235             :       }
    2236             : 
    2237          21 :       p1 = gmul2n(gmul(gel(x,2),gel(x,3)), 1);
    2238          21 :       tetpil = avma;
    2239          21 :       gel(z,2) = gadd(p2,p4);
    2240          21 :       gel(z,3) = gadd(p1,p3);
    2241          21 :       gerepilecoeffssp(av,tetpil,z+2,2); return z;
    2242             : 
    2243             :     case t_POLMOD:
    2244       24107 :       return sqr_polmod(gel(x,1), gel(x,2));
    2245             : 
    2246     2778292 :     case t_FFELT: return FF_sqr(x);
    2247             : 
    2248             :     case t_POL:
    2249             :     {
    2250     1728409 :       GEN a = x, p = NULL, pol = NULL;
    2251     1728409 :       av = avma;
    2252     1728409 :       if (RgX_is_ZX(x)) return ZX_sqr(x);
    2253       30359 :       if (ff_poltype(&x,&p,&pol))
    2254             :       {
    2255         189 :         z = ZX_sqr(x);
    2256         189 :         if (p) z = FpX_to_mod(z,p);
    2257         189 :         if (pol) z = Kronecker_to_mod(z,pol);
    2258         189 :         z = gerepileupto(av, z);
    2259             :       }
    2260       30170 :       else { avma = av; z = RgX_sqr(a); }
    2261       30359 :       return z;
    2262             :     }
    2263             : 
    2264             :     case t_SER:
    2265        6664 :       lx = lg(x);
    2266        6664 :       if (ser_isexactzero(x)) {
    2267           7 :         GEN z = gcopy(x);
    2268           7 :         setvalp(z, 2*valp(x));
    2269           7 :         return z;
    2270             :       }
    2271        6657 :       if (lx < 40)
    2272        6510 :         return normalize( sqr_ser_part(x, 0, lx-3) );
    2273             :       else
    2274             :       {
    2275         147 :         pari_sp av = avma;
    2276         147 :         GEN z = cgetg(lx,t_SER), p = NULL;
    2277         147 :         z[1] = evalvalp(2*valp(x)) | evalvarn(varn(x)) | evalsigne(1);
    2278         147 :         x = ser2pol_i(x,lx);
    2279         147 :         if (RgX_is_FpX(x,&p))
    2280             :         {
    2281          49 :           if (!p) x = ZX_sqr(x);
    2282             :           else
    2283             :           {
    2284           0 :             x = RgX_to_FpX(x, p);
    2285           0 :             x = FpX_to_mod(ZX_sqr(x), p);
    2286             :           }
    2287             :         }
    2288             :         else
    2289          98 :           x = RgXn_sqr(x, lx-2);
    2290         147 :         z = fill_ser(z, x);
    2291         147 :         return gerepilecopy(av, z);
    2292             :       }
    2293             : 
    2294           7 :     case t_RFRAC: z = cgetg(3,t_RFRAC);
    2295           7 :       gel(z,1) = gsqr(gel(x,1));
    2296           7 :       gel(z,2) = gsqr(gel(x,2)); return z;
    2297             : 
    2298         791 :     case t_MAT: return RgM_sqr(x);
    2299          14 :     case t_QFR: return qfrsqr(x);
    2300    31383556 :     case t_QFI: return qfisqr(x);
    2301             :     case t_VECSMALL:
    2302        3458 :       z = cgetg_copy(x, &lx);
    2303       56392 :       for (i=1; i<lx; i++)
    2304             :       {
    2305       52934 :         long xi = x[i];
    2306       52934 :         if (xi < 1 || xi >= lx) pari_err_TYPE2("*",x,x);
    2307       52934 :         z[i] = x[xi];
    2308             :       }
    2309        3458 :       return z;
    2310             :   }
    2311           0 :   pari_err_TYPE2("*",x,x);
    2312           0 :   return NULL; /* not reached */
    2313             : }
    2314             : 
    2315             : /********************************************************************/
    2316             : /**                                                                **/
    2317             : /**                           DIVISION                             **/
    2318             : /**                                                                **/
    2319             : /********************************************************************/
    2320             : static GEN
    2321      270646 : div_rfrac_scal(GEN x, GEN y)
    2322             : {
    2323      270646 :   pari_sp av = avma;
    2324      270646 :   GEN d = rfrac_denom_mul_scal(gel(x,2), y);
    2325      270646 :   return gerepileupto(av, gred_rfrac_simple(gel(x,1), d));
    2326             : }
    2327             : static GEN
    2328     1004057 : div_scal_rfrac(GEN x, GEN y)
    2329             : {
    2330     1004057 :   GEN y1 = gel(y,1), y2 = gel(y,2);
    2331     1004057 :   pari_sp av = avma;
    2332     1004057 :   if (typ(y1) == t_POL && varn(y2) == varn(y1))
    2333             :   {
    2334           0 :     if (degpol(y1)) return gerepileupto(av, gred_rfrac_simple(gmul(x, y2), y1));
    2335           0 :     y1 = gel(y1,2);
    2336             :   }
    2337     1004057 :   return RgX_Rg_mul(y2, gdiv(x,y1));
    2338             : }
    2339             : static GEN
    2340     1186329 : div_rfrac(GEN x, GEN y)
    2341     1186329 : { return mul_rfrac(gel(x,1),gel(x,2), gel(y,2),gel(y,1)); }
    2342             : 
    2343             : static GEN
    2344      456227 : div_ser_scal(GEN x, GEN y) {
    2345             :   long i, lx;
    2346             :   GEN z;
    2347      456227 :   if (ser_isexactzero(x))
    2348             :   {
    2349           7 :     if (lg(x) == 2) return gcopy(x);
    2350           0 :     return scalarser(gdiv(gel(x,2), y), varn(x), valp(x));
    2351             :   }
    2352      456220 :   z = cgetg_copy(x, &lx); z[1] = x[1];
    2353      456220 :   for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2354      456220 :   return normalize(z);
    2355             : }
    2356             : GEN
    2357         133 : ser_normalize(GEN x)
    2358             : {
    2359         133 :   long i, lx = lg(x);
    2360             :   GEN c, z;
    2361         133 :   if (lx == 2) return x;
    2362         133 :   c = gel(x,2); if (gequal1(c)) return x;
    2363         105 :   z = cgetg(lx, t_SER); z[1] = x[1]; gel(z,2) = gen_1;
    2364         105 :   for (i=3; i<lx; i++) gel(z,i) = gdiv(gel(x,i),c);
    2365         105 :   return z;
    2366             : }
    2367             : 
    2368             : static GEN
    2369     1047007 : div_T_scal(GEN x, GEN y, long tx) {
    2370     1047007 :   switch(tx)
    2371             :   {
    2372      323858 :     case t_POL: return RgX_Rg_div(x, y);
    2373      456227 :     case t_SER: return div_ser_scal(x, y);
    2374      266922 :     case t_RFRAC: return div_rfrac_scal(x,y);
    2375             :   }
    2376           0 :   pari_err_TYPE2("/",x,y);
    2377           0 :   return NULL; /* not reached */
    2378             : }
    2379             : 
    2380             : static GEN
    2381     8183754 : div_scal_pol(GEN x, GEN y) {
    2382     8183754 :   long ly = lg(y);
    2383             :   pari_sp av;
    2384     8183754 :   if (ly == 3) return scalarpol(gdiv(x,gel(y,2)), varn(y));
    2385     8151881 :   if (isrationalzero(x)) return zeropol(varn(y));
    2386     6099930 :   av = avma;
    2387     6099930 :   return gerepileupto(av, gred_rfrac_simple(x,y));
    2388             : }
    2389             : static GEN
    2390        9569 : div_scal_ser(GEN x, GEN y) { /* TODO: improve */
    2391             :   GEN z;
    2392             :   long ly, i;
    2393        9569 :   if (gequal0(x)) { pari_sp av=avma; return gerepileupto(av, gmul(x, ginv(y))); }
    2394        9569 :   ly = lg(y); z = (GEN)pari_malloc(ly*sizeof(long));
    2395        9569 :   z[0] = evaltyp(t_SER) | evallg(ly);
    2396        9569 :   z[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    2397        9569 :   gel(z,2) = x; for (i=3; i<ly; i++) gel(z,i) = gen_0;
    2398        9569 :   y = gdiv(z,y); pari_free(z); return y;
    2399             : }
    2400             : static GEN
    2401     9147722 : div_scal_T(GEN x, GEN y, long ty) {
    2402     9147722 :   switch(ty)
    2403             :   {
    2404     8134852 :     case t_POL: return div_scal_pol(x, y);
    2405        9569 :     case t_SER: return div_scal_ser(x, y);
    2406     1003301 :     case t_RFRAC: return div_scal_rfrac(x, y);
    2407             :   }
    2408           0 :   pari_err_TYPE2("/",x,y);
    2409           0 :   return NULL; /* not reached */
    2410             : }
    2411             : 
    2412             : /* assume tx = ty = t_SER, same variable vx */
    2413             : static GEN
    2414      735625 : div_ser(GEN x, GEN y, long vx)
    2415             : {
    2416      735625 :   long i, j, l = valp(x) - valp(y), lx = lg(x), ly = lg(y);
    2417             :   GEN y_lead, p1, p2, z;
    2418             : 
    2419      735625 :   if (!signe(y)) pari_err_INV("div_ser", y);
    2420      735618 :   if (ser_isexactzero(x))
    2421             :   {
    2422         651 :     if (lx == 2) return zeroser(vx, l);
    2423           7 :     return scalarser(gmul(gel(x,2),RgX_get_0(y)), varn(x), l);
    2424             :   }
    2425      734967 :   y_lead = gel(y,2);
    2426      734967 :   if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
    2427             :   {
    2428           0 :     pari_warn(warner,"normalizing a series with 0 leading term");
    2429           0 :     for (l--, ly--,y++; ly > 2; l--, ly--, y++)
    2430             :     {
    2431           0 :       y_lead = gel(y,2);
    2432           0 :       if (!gequal0(y_lead)) break;
    2433             :     }
    2434           0 :     if (ly <= 2) pari_err_INV("div_ser", y);
    2435             :   }
    2436      734967 :   if (ly < lx) lx = ly;
    2437      734967 :   p2 = cgetg(lx, t_VECSMALL); /* left on stack for efficiency */
    2438     2317229 :   for (i=3; i<lx; i++)
    2439             :   {
    2440     1582262 :     p1 = gel(y,i);
    2441     1582262 :     if (isrationalzero(p1)) p1 = NULL;
    2442     1582262 :     gel(p2,i) = p1;
    2443             :   }
    2444      734967 :   z = cgetg(lx,t_SER);
    2445      734967 :   z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
    2446      734967 :   gel(z,2) = gdiv(gel(x,2), y_lead);
    2447     2317229 :   for (i=3; i<lx; i++)
    2448             :   {
    2449     1582262 :     pari_sp av = avma;
    2450     1582262 :     p1 = gel(x,i);
    2451    11377748 :     for (j=2, l=i; j<i; j++, l--)
    2452     9795486 :       if (p2[l]) p1 = gsub(p1, gmul(gel(z,j), gel(p2,l)));
    2453     1582262 :     gel(z,i) = gerepileupto(av, gdiv(p1, y_lead));
    2454             :   }
    2455      734967 :   return normalize(z);
    2456             : }
    2457             : /* x,y compatible PADIC */
    2458             : static GEN
    2459      216064 : divpp(GEN x, GEN y) {
    2460             :   pari_sp av;
    2461             :   long a, b;
    2462             :   GEN z, M;
    2463             : 
    2464      216064 :   if (!signe(gel(y,4))) pari_err_INV("divpp",y);
    2465      216057 :   if (!signe(gel(x,4))) return zeropadic(gel(x,2), valp(x)-valp(y));
    2466      215357 :   a = precp(x);
    2467      215357 :   b = precp(y); if (a > b) { M = gel(y,3); } else { M = gel(x,3); b = a; }
    2468      215357 :   z = cgetg(5, t_PADIC);
    2469      215362 :   z[1] = _evalprecp(b) | evalvalp(valp(x) - valp(y));
    2470      215360 :   gel(z,2) = icopy(gel(x,2));
    2471      215362 :   gel(z,3) = icopy(M); av = avma;
    2472      215363 :   gel(z,4) = gerepileuptoint(av, remii(mulii(gel(x,4), Fp_inv(gel(y,4), M)), M) );
    2473      215363 :   return z;
    2474             : }
    2475             : static GEN
    2476        6789 : div_polmod_same(GEN T, GEN x, GEN y)
    2477             : {
    2478        6789 :   long v = varn(T);
    2479        6789 :   GEN a, z = cgetg(3, t_POLMOD);
    2480        6789 :   gel(z,1) = RgX_copy(T);
    2481        6789 :   if (typ(y) != t_POL || varn(y) != v || lg(y) <= 3)
    2482        4767 :     a = gdiv(x, y);
    2483        2022 :   else if (typ(x) != t_POL || varn(x) != v || lg(x) <= 3)
    2484         189 :   {
    2485         189 :     pari_sp av = avma;
    2486         189 :     a = gerepileupto(av, gmul(x, RgXQ_inv(y, T)));
    2487             :   }
    2488        1833 :   else if (degpol(T) == 2 && isint1(gel(T,4))) /* quadratic fields */
    2489        1365 :   {
    2490        1365 :     pari_sp av = avma;
    2491        1365 :     a = quad_polmod_mul(T, x, quad_polmod_conj(y, T));
    2492        1365 :     a = RgX_Rg_div(a, quad_polmod_norm(y, T));
    2493        1365 :     a = gerepileupto(av, a);
    2494             :   }
    2495             :   else
    2496             :   {
    2497         468 :     pari_sp av = avma;
    2498         468 :     GEN p = NULL;
    2499         468 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
    2500             :     {
    2501          14 :       T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
    2502          14 :       if (lgefint(p) == 3)
    2503             :       {
    2504           7 :         ulong pp = p[2];
    2505           7 :         x = ZX_to_Flx(x, pp);
    2506           7 :         y = ZX_to_Flx(y, pp);
    2507           7 :         T = ZX_to_Flx(T, pp);
    2508           7 :         a = Flxq_mul(x, Flxq_inv(y,T,pp),T,pp);
    2509           7 :         a = Flx_to_ZX(a);
    2510             :       }
    2511             :       else
    2512           7 :         a = FpXQ_div(x, y, T, p);
    2513          14 :       a = FpX_to_mod(a, p);
    2514             :     } else
    2515         454 :       a = RgXQ_mul(x, ginvmod(y, gel(z,1)), gel(z,1));
    2516         468 :     a = gerepileupto(av, a);
    2517             :   }
    2518        6789 :   gel(z,2) = a; return z;
    2519             : }
    2520             : GEN
    2521   128951273 : gdiv(GEN x, GEN y)
    2522             : {
    2523   128951273 :   long tx = typ(x), ty = typ(y), lx, ly, vx, vy, i;
    2524             :   pari_sp av, tetpil;
    2525             :   GEN z, p1, p2;
    2526             : 
    2527   128951273 :   if (tx == ty) switch(tx)
    2528             :   {
    2529             :     case t_INT:
    2530    56065186 :       return Qdivii(x,y);
    2531             : 
    2532    27334306 :     case t_REAL: return divrr(x,y);
    2533       21318 :     case t_INTMOD: { GEN X = gel(x,1), Y = gel(y,1);
    2534       21318 :       z = cgetg(3,t_INTMOD);
    2535       21318 :       if (X==Y || equalii(X,Y))
    2536       21312 :         return div_intmod_same(z, X, gel(x,2), gel(y,2));
    2537           7 :       gel(z,1) = gcdii(X,Y);
    2538           7 :       warn_coercion(X,Y,gel(z,1));
    2539           7 :       av = avma; p1 = mulii(gel(x,2), Fp_inv(gel(y,2), gel(z,1)));
    2540           7 :       gel(z,2) = gerepileuptoint(av, remii(p1, gel(z,1))); return z;
    2541             :     }
    2542             :     case t_FRAC: {
    2543      224481 :       GEN x1 = gel(x,1), x2 = gel(x,2);
    2544      224481 :       GEN y1 = gel(y,1), y2 = gel(y,2);
    2545      224481 :       z = cgetg(3, t_FRAC);
    2546      224481 :       p1 = gcdii(x1, y1);
    2547      224481 :       if (!is_pm1(p1)) { x1 = diviiexact(x1,p1); y1 = diviiexact(y1,p1); }
    2548      224481 :       p1 = gcdii(x2, y2);
    2549      224481 :       if (!is_pm1(p1)) { x2 = diviiexact(x2,p1); y2 = diviiexact(y2,p1); }
    2550      224481 :       tetpil = avma;
    2551      224481 :       gel(z,2) = mulii(x2,y1);
    2552      224481 :       gel(z,1) = mulii(x1,y2);
    2553      224481 :       normalize_frac(z);
    2554      224481 :       fix_frac_if_int_GC(z,tetpil);
    2555      224481 :       return z;
    2556             :     }
    2557             :     case t_COMPLEX:
    2558     1355784 :       if (isintzero(gel(y,1)))
    2559             :       {
    2560        8904 :         y = gel(y,2);
    2561        8904 :         if (isintzero(gel(x,1))) return gdiv(gel(x,2), y);
    2562        8141 :         z = cgetg(3,t_COMPLEX);
    2563        8141 :         gel(z,1) = gdiv(gel(x,2), y);
    2564        8141 :         av = avma;
    2565        8141 :         gel(z,2) = gerepileupto(av, gneg(gdiv(gel(x,1), y)));
    2566        8141 :         return z;
    2567             :       }
    2568     1346880 :       av = avma; p1 = cxnorm(y); p2 = mulcc(x, gconj(y)); tetpil = avma;
    2569     1346880 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2570             : 
    2571             :     case t_PADIC:
    2572      118664 :       if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("/",x,y);
    2573      118664 :       return divpp(x, y);
    2574             : 
    2575             :     case t_QUAD:
    2576         322 :       if (!ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("/",x,y);
    2577         322 :       av = avma; p1 = quadnorm(y); p2 = mulqq(x, gconj(y)); tetpil = avma;
    2578         322 :       return gerepile(av, tetpil, gdiv(p2,p1));
    2579             : 
    2580      111860 :     case t_FFELT: return FF_div(x,y);
    2581             : 
    2582             :     case t_POLMOD:
    2583        6796 :       if (RgX_equal_var(gel(x,1), gel(y,1)))
    2584        6789 :         z = div_polmod_same(gel(x,1), gel(x,2), gel(y,2));
    2585             :       else {
    2586           7 :         av = avma;
    2587           7 :         z = gerepileupto(av, gmul(x, ginv(y)));
    2588             :       }
    2589        6796 :       return z;
    2590             : 
    2591             :     case t_POL:
    2592    18381420 :       vx = varn(x);
    2593    18381420 :       vy = varn(y);
    2594    18381420 :       if (vx != vy) {
    2595       97846 :         if (varncmp(vx, vy) < 0) return RgX_Rg_div(x, y);
    2596       48902 :                             else return div_scal_pol(x, y);
    2597             :       }
    2598    18283574 :       if (!signe(y)) pari_err_INV("gdiv",y);
    2599    18283574 :       if (lg(y) == 3) return RgX_Rg_div(x,gel(y,2));
    2600    18174684 :       av = avma;
    2601    18174684 :       return gerepileupto(av, gred_rfrac2(x,y));
    2602             : 
    2603             :     case t_SER:
    2604      400247 :       vx = varn(x);
    2605      400247 :       vy = varn(y);
    2606      400247 :       if (vx != vy) {
    2607           0 :         if (varncmp(vx, vy) < 0) return div_ser_scal(x, y);
    2608           0 :                             else return div_scal_ser(x, y);
    2609             :       }
    2610      400247 :       return div_ser(x, y, vx);
    2611             :     case t_RFRAC:
    2612     1190809 :       vx = varn(gel(x,2));
    2613     1190809 :       vy = varn(gel(y,2));
    2614     1190809 :       if (vx != vy) {
    2615        4480 :         if (varncmp(vx, vy) < 0) return div_rfrac_scal(x, y);
    2616         756 :                             else return div_scal_rfrac(x, y);
    2617             :       }
    2618     1186329 :       return div_rfrac(x,y);
    2619             : 
    2620           0 :     case t_QFI: av = avma; return gerepileupto(av, qficomp(x, ginv(y)));
    2621           0 :     case t_QFR: av = avma; return gerepileupto(av, qfrcomp(x, ginv(y)));
    2622             : 
    2623             :     case t_MAT:
    2624          14 :       av = avma; p1 = RgM_inv(y);
    2625          14 :       if (!p1) pari_err_INV("gdiv",y);
    2626          14 :       return gerepileupto(av, RgM_mul(x, p1));
    2627             : 
    2628           0 :     default: pari_err_TYPE2("/",x,y);
    2629             :   }
    2630             : 
    2631    23740066 :   if (tx==t_INT && is_const_t(ty)) /* optimized for speed */
    2632             :   {
    2633      419402 :     long s = signe(x);
    2634      419402 :     if (!s) {
    2635       31380 :       if (gequal0(y)) pari_err_INV("gdiv",y);
    2636       31380 :       switch (ty)
    2637             :       {
    2638       28237 :       default: return gen_0;
    2639             :       case t_INTMOD:
    2640          14 :         z = cgetg(3,t_INTMOD);
    2641          14 :         gel(z,1) = icopy(gel(y,1));
    2642          14 :         gel(z,2) = gen_0; return z;
    2643        3129 :       case t_FFELT: return FF_zero(y);
    2644             :       }
    2645             :     }
    2646      388022 :     if (is_pm1(x)) {
    2647      296797 :       if (s > 0) return ginv(y);
    2648       30760 :       av = avma; return gerepileupto(av, ginv(gneg(y)));
    2649             :     }
    2650       91225 :     switch(ty)
    2651             :     {
    2652       39873 :       case t_REAL: return divir(x,y);
    2653             :       case t_INTMOD:
    2654          21 :         z = cgetg(3, t_INTMOD);
    2655          21 :         return div_intmod_same(z, gel(y,1), modii(x, gel(y,1)), gel(y,2));
    2656             :       case t_FRAC:
    2657       10215 :         z = cgetg(3,t_FRAC); p1 = gcdii(x,gel(y,1));
    2658       10215 :         if (is_pm1(p1))
    2659             :         {
    2660        6142 :           avma = (pari_sp)z;
    2661        6142 :           gel(z,2) = icopy(gel(y,1));
    2662        6142 :           gel(z,1) = mulii(gel(y,2), x);
    2663        6142 :           normalize_frac(z);
    2664        6142 :           fix_frac_if_int(z);
    2665             :         }
    2666             :         else
    2667             :         {
    2668        4073 :           x = diviiexact(x,p1); tetpil = avma;
    2669        4073 :           gel(z,2) = diviiexact(gel(y,1), p1);
    2670        4073 :           gel(z,1) = mulii(gel(y,2), x);
    2671        4073 :           normalize_frac(z);
    2672        4073 :           fix_frac_if_int_GC(z,tetpil);
    2673             :         }
    2674       10215 :         return z;
    2675             : 
    2676         231 :       case t_FFELT: return Z_FF_div(x,y);
    2677       28901 :       case t_COMPLEX: return divRc(x,y);
    2678       11942 :       case t_PADIC: return divTp(x, y);
    2679             :       case t_QUAD:
    2680          42 :         av = avma; p1 = quadnorm(y); p2 = mulRq(x, gconj(y)); tetpil = avma;
    2681          42 :         return gerepile(av, tetpil, gdiv(p2,p1));
    2682             :     }
    2683             :   }
    2684    23320664 :   if (gequal0(y))
    2685             :   {
    2686          35 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    2687          14 :     if (ty != t_MAT) pari_err_INV("gdiv",y);
    2688             :   }
    2689             : 
    2690    23320636 :   if (is_const_t(tx) && is_const_t(ty)) switch(tx)
    2691             :   {
    2692             :     case t_REAL:
    2693     9069948 :       switch(ty)
    2694             :       {
    2695     8069111 :         case t_INT: return divri(x,y);
    2696             :         case t_FRAC:
    2697      668576 :           av = avma; z = divri(mulri(x,gel(y,2)), gel(y,1));
    2698      668576 :           return gerepileuptoleaf(av, z);
    2699      332233 :         case t_COMPLEX: return divRc(x, y);
    2700           7 :         case t_QUAD: return divfq(x, y, lg(x));
    2701          21 :         default: pari_err_TYPE2("/",x,y);
    2702             :       }
    2703             : 
    2704             :     case t_INTMOD:
    2705         490 :       switch(ty)
    2706             :       {
    2707             :         case t_INT:
    2708         357 :           z = cgetg(3, t_INTMOD);
    2709         357 :           return div_intmod_same(z, gel(x,1), gel(x,2), modii(y, gel(x,1)));
    2710          28 :         case t_FRAC: { GEN X = gel(x,1);
    2711          28 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2), gel(x,2)), X);
    2712          28 :           return div_intmod_same(z, X, p1, modii(gel(y,1), X));
    2713             :         }
    2714             :         case t_FFELT:
    2715          63 :           if (!equalii(gel(x,1),FF_p_i(y)))
    2716           0 :             pari_err_OP("/",x,y);
    2717          63 :           return Z_FF_div(gel(x,2),y);
    2718             : 
    2719             :         case t_COMPLEX:
    2720           0 :           av = avma;
    2721           0 :           return gerepileupto(av, mulRc_direct(gdiv(x,cxnorm(y)), gconj(y)));
    2722             : 
    2723             :         case t_QUAD:
    2724           7 :           av = avma; p1 = quadnorm(y); p2 = gmul(x,gconj(y)); tetpil = avma;
    2725           7 :           return gerepile(av,tetpil, gdiv(p2,p1));
    2726             : 
    2727           7 :         case t_PADIC: { GEN X = gel(x,1);
    2728           7 :           z = cgetg(3, t_INTMOD);
    2729           7 :           return div_intmod_same(z, X, gel(x,2), padic_to_Fp(y, X));
    2730             :         }
    2731          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2732             :       }
    2733             : 
    2734             :     case t_FRAC:
    2735      310879 :       switch(ty)
    2736             :       {
    2737      276282 :         case t_INT: z = cgetg(3, t_FRAC);
    2738      276282 :         p1 = gcdii(y,gel(x,1));
    2739      276282 :         if (is_pm1(p1))
    2740             :         {
    2741      166996 :           avma = (pari_sp)z; tetpil = 0;
    2742      166996 :           gel(z,1) = icopy(gel(x,1));
    2743             :         }
    2744             :         else
    2745             :         {
    2746      109286 :           y = diviiexact(y,p1); tetpil = avma;
    2747      109286 :           gel(z,1) = diviiexact(gel(x,1), p1);
    2748             :         }
    2749      276282 :         gel(z,2) = mulii(gel(x,2),y);
    2750      276282 :         normalize_frac(z);
    2751      276282 :         if (tetpil) fix_frac_if_int_GC(z,tetpil);
    2752      276282 :         return z;
    2753             : 
    2754             :         case t_REAL:
    2755       32273 :           av=avma; p1=mulri(y,gel(x,2)); tetpil=avma;
    2756       32273 :           return gerepile(av, tetpil, divir(gel(x,1), p1));
    2757             : 
    2758           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2759           7 :           z = cgetg(3,t_INTMOD); p1 = remii(mulii(gel(y,2),gel(x,2)), Y);
    2760           7 :           return div_intmod_same(z, Y, modii(gel(x,1), Y), p1);
    2761             :         }
    2762             : 
    2763          28 :         case t_FFELT: av=avma;
    2764          28 :           return gerepileupto(av,Z_FF_div(gel(x,1),FF_Z_mul(y,gel(x,2))));
    2765             : 
    2766         231 :         case t_COMPLEX: return divRc(x, y);
    2767             : 
    2768             :         case t_PADIC:
    2769        2051 :           if (!signe(gel(x,1))) return gen_0;
    2770        2051 :           return divTp(x, y);
    2771             : 
    2772             :         case t_QUAD:
    2773           7 :           av=avma; p1=quadnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
    2774           7 :           return gerepile(av,tetpil,gdiv(p2,p1));
    2775             :       }
    2776             : 
    2777             :     case t_FFELT:
    2778        1939 :       switch (ty)
    2779             :       {
    2780        1820 :         case t_INT: return FF_Z_Z_muldiv(x,gen_1,y);
    2781          28 :         case t_FRAC: return FF_Z_Z_muldiv(x,gel(y,2),gel(y,1));
    2782             :         case t_INTMOD:
    2783          49 :           if (!equalii(gel(y,1),FF_p_i(x)))
    2784           0 :             pari_err_OP("/",x,y);
    2785          49 :           return FF_Z_Z_muldiv(x,gen_1,gel(y,2));
    2786             :         default:
    2787          42 :         pari_err_TYPE2("/",x,y);
    2788             :       }
    2789           0 :       break;
    2790             : 
    2791             :     case t_COMPLEX:
    2792     2346525 :       switch(ty)
    2793             :       {
    2794     2346525 :         case t_INT: case t_REAL: case t_FRAC: return divcR(x,y);
    2795           0 :         case t_INTMOD: return mulRc_direct(ginv(y), x);
    2796             :         case t_PADIC:
    2797           0 :           return Zp_nosquare_m1(gel(y,2))? divcR(x,y): divTp(x, y);
    2798             :         case t_QUAD:
    2799           0 :           lx = precision(x); if (!lx) pari_err_OP("/",x,y);
    2800           0 :           return divfq(x, y, lx);
    2801             :       }
    2802             : 
    2803             :     case t_PADIC:
    2804        1351 :       switch(ty)
    2805             :       {
    2806        1302 :         case t_INT: case t_FRAC: { GEN p = gel(x,2);
    2807        2604 :           return signe(gel(x,4))? divpT(x, y)
    2808        1302 :                             : zeropadic(p, valp(x) - Q_pval(y,p));
    2809             :         }
    2810           7 :         case t_INTMOD: { GEN Y = gel(y,1);
    2811           7 :           z = cgetg(3, t_INTMOD);
    2812           7 :           return div_intmod_same(z, Y, padic_to_Fp(x, Y), gel(y,2));
    2813             :         }
    2814             :         case t_COMPLEX: case t_QUAD:
    2815           7 :           av=avma; p1=gmul(x,gconj(y)); p2=gnorm(y); tetpil=avma;
    2816           7 :           return gerepile(av,tetpil,gdiv(p1,p2));
    2817             : 
    2818          28 :         case t_REAL: pari_err_TYPE2("/",x,y);
    2819             :       }
    2820             : 
    2821             :     case t_QUAD:
    2822         889 :       switch (ty)
    2823             :       {
    2824             :         case t_INT: case t_INTMOD: case t_FRAC:
    2825         840 :           z = cgetg(4,t_QUAD);
    2826         840 :           gel(z,1) = ZX_copy(gel(x,1));
    2827         840 :           gel(z,2) = gdiv(gel(x,2), y);
    2828         840 :           gel(z,3) = gdiv(gel(x,3), y); return z;
    2829          28 :         case t_REAL: return divqf(x, y, lg(y));
    2830           7 :         case t_PADIC: return divTp(x, y);
    2831             :         case t_COMPLEX:
    2832           0 :           ly = precision(y); if (!ly) pari_err_OP("/",x,y);
    2833           0 :           return divqf(x, y, ly);
    2834             :       }
    2835             :   }
    2836    11588636 :   switch(ty) {
    2837             :     case t_REAL: case t_INTMOD: case t_PADIC: case t_POLMOD:
    2838       22806 :       return gmul(x, ginv(y)); /* missing gerepile, for speed */
    2839             :     case t_MAT:
    2840          28 :       av = avma; p1 = RgM_inv(y);
    2841          28 :       if (!p1) pari_err_INV("gdiv",y);
    2842          21 :       return gerepileupto(av, gmul(x, p1));
    2843             :     case t_VEC: case t_COL:
    2844             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2845           0 :       pari_err_TYPE2("/",x,y);
    2846             :   }
    2847    11565802 :   switch(tx) {
    2848             :     case t_VEC: case t_COL: case t_MAT:
    2849      692108 :       z = cgetg_copy(x, &lx);
    2850      692108 :       for (i=1; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
    2851      692108 :       return z;
    2852             :     case t_LIST: case t_STR: case t_VECSMALL: case t_CLOSURE:
    2853           0 :       pari_err_TYPE2("/",x,y);
    2854             :   }
    2855             : 
    2856    10873694 :   vy = gvar(y);
    2857    10873694 :   if (tx == t_POLMOD) { GEN X = gel(x,1);
    2858       11046 :     vx = varn(X);
    2859       11046 :     if (vx != vy) {
    2860       10899 :       if (varncmp(vx, vy) > 0) return div_scal_T(x, y, ty);
    2861       10668 :       z = cgetg(3,t_POLMOD);
    2862       10668 :       gel(z,1) = RgX_copy(X);
    2863       10668 :       gel(z,2) = gdiv(gel(x,2), y); return z;
    2864             :     }
    2865             :     /* y is POL, SER or RFRAC */
    2866         147 :     av = avma;
    2867         147 :     switch(ty)
    2868             :     {
    2869           0 :       case t_RFRAC: y = gmod(ginv(y), X); break;
    2870         147 :       default: y = ginvmod(gmod(y,X), X);
    2871             :     }
    2872         140 :     return gerepileupto(av, mul_polmod_same(X, gel(x,2), y));
    2873             :   }
    2874             :   /* x and y are not both is_scalar_t. If one of them is scalar, it's not a
    2875             :    * POLMOD (done already), hence its variable is NO_VARIABLE. If the other has
    2876             :    * variable NO_VARIABLE, then the operation is incorrect */
    2877    10862648 :   vx = gvar(x);
    2878    10862648 :   if (vx != vy) { /* includes cases where one is scalar */
    2879    10194498 :     if (varncmp(vx, vy) < 0) return div_T_scal(x, y, tx);
    2880     9147491 :                         else return div_scal_T(x, y, ty);
    2881             :   }
    2882      668150 :   switch(tx)
    2883             :   {
    2884             :     case t_POL:
    2885       79554 :       switch(ty)
    2886             :       {
    2887             :         case t_SER:
    2888          28 :           if (lg(y) == 2)
    2889           0 :             return zeroser(vx, RgX_val(x) - valp(y));
    2890          28 :           p1 = RgX_to_ser(x,lg(y));
    2891          28 :           p2 = div_ser(p1, y, vx);
    2892          28 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2893          28 :           return p2;
    2894             : 
    2895             :         case t_RFRAC:
    2896             :         {
    2897       79526 :           GEN y1 = gel(y,1), y2 = gel(y,2);
    2898       79526 :           if (typ(y1) == t_POL && varn(y1) == vx)
    2899           0 :             return mul_rfrac_scal(y2, y1, x);
    2900       79526 :           av = avma;
    2901       79526 :           return gerepileupto(av, RgX_Rg_div(RgX_mul(y2, x), y1));
    2902             :         }
    2903             :       }
    2904           0 :       break;
    2905             : 
    2906             :     case t_SER:
    2907      335455 :       switch(ty)
    2908             :       {
    2909             :         case t_POL:
    2910      335455 :           if (lg(x) == 2)
    2911         105 :             return zeroser(vx, valp(x) - RgX_val(y));
    2912      335350 :           p1 = RgX_to_ser_inexact(y,lg(x));
    2913      335350 :           p2 = div_ser(x, p1, vx);
    2914      335350 :           settyp(p1, t_VECSMALL); /* p1 left on stack */
    2915      335350 :           return p2;
    2916             :         case t_RFRAC:
    2917           0 :           av = avma;
    2918           0 :           return gerepileupto(av, gdiv(gmul(x,gel(y,2)), gel(y,1)));
    2919             :       }
    2920           0 :       break;
    2921             : 
    2922             :     case t_RFRAC:
    2923      253127 :       switch(ty)
    2924             :       {
    2925      253127 :         case t_POL: return div_rfrac_pol(gel(x,1),gel(x,2), y);
    2926             :         case t_SER:
    2927           0 :           av = avma; z = RgX_to_ser_inexact(gel(x,2), lg(y));
    2928           0 :           return gerepileupto(av, gdiv(gel(x,1), gmul(z,y)));
    2929             :       }
    2930           0 :       break;
    2931             :   }
    2932          14 :   pari_err_TYPE2("/",x,y);
    2933           0 :   return NULL; /* not reached */
    2934             : }
    2935             : 
    2936             : /********************************************************************/
    2937             : /**                                                                **/
    2938             : /**                     SIMPLE MULTIPLICATION                      **/
    2939             : /**                                                                **/
    2940             : /********************************************************************/
    2941             : GEN
    2942    23801596 : gmulsg(long s, GEN y)
    2943             : {
    2944             :   long ly, i;
    2945             :   pari_sp av;
    2946             :   GEN z;
    2947             : 
    2948    23801596 :   switch(typ(y))
    2949             :   {
    2950    11808377 :     case t_INT:  return mulsi(s,y);
    2951     8634609 :     case t_REAL: return mulsr(s,y);
    2952      168494 :     case t_INTMOD: { GEN p = gel(y,1);
    2953      168494 :       z = cgetg(3,t_INTMOD);
    2954      168495 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(mulsi(s,gel(y,2)), p));
    2955      168495 :       gel(z,1) = icopy(p); return z;
    2956             :     }
    2957      626307 :     case t_FFELT: return FF_Z_mul(y,stoi(s));
    2958             :     case t_FRAC:
    2959      928767 :       if (!s) return gen_0;
    2960      927864 :       z = cgetg(3,t_FRAC);
    2961      927864 :       i = cgcd(s, smodis(gel(y,2), s));
    2962      927864 :       if (i == 1)
    2963             :       {
    2964      334502 :         gel(z,2) = icopy(gel(y,2));
    2965      334502 :         gel(z,1) = mulis(gel(y,1), s);
    2966             :       }
    2967             :       else
    2968             :       {
    2969      593362 :         gel(z,2) = divis(gel(y,2), i);
    2970      593362 :         gel(z,1) = mulis(gel(y,1), s/i);
    2971      593362 :         fix_frac_if_int(z);
    2972             :       }
    2973      927864 :       return z;
    2974             : 
    2975     1128545 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    2976     1128545 :       gel(z,1) = gmulsg(s,gel(y,1));
    2977     1128545 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2978             : 
    2979             :     case t_PADIC:
    2980        9093 :       if (!s) return gen_0;
    2981        9093 :       av = avma; return gerepileupto(av, mulpp(cvtop2(stoi(s),y), y));
    2982             : 
    2983           7 :     case t_QUAD: z = cgetg(4, t_QUAD);
    2984           7 :       gel(z,1) = ZX_copy(gel(y,1));
    2985           7 :       gel(z,2) = gmulsg(s,gel(y,2));
    2986           7 :       gel(z,3) = gmulsg(s,gel(y,3)); return z;
    2987             : 
    2988        6223 :     case t_POLMOD: z = cgetg(3, t_POLMOD);
    2989        6223 :       gel(z,1) = RgX_copy(gel(y,1));
    2990        6223 :       gel(z,2) = gmulsg(s,gel(y,2)); return z;
    2991             : 
    2992             :     case t_POL:
    2993      322228 :       if (!signe(y)) return RgX_copy(y);
    2994      318987 :       if (!s) return scalarpol(RgX_get_0(y), varn(y));
    2995      315186 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    2996      315186 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    2997      315186 :       return normalizepol_lg(z, ly);
    2998             : 
    2999             :     case t_SER:
    3000         819 :       if (ser_isexactzero(y)) return gcopy(y);
    3001         819 :       if (!s) return RgX_get_0(y);
    3002         819 :       z = cgetg_copy(y, &ly); z[1]=y[1];
    3003         819 :       for (i=2; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    3004         819 :       return normalize(z);
    3005             : 
    3006             :     case t_RFRAC:
    3007          14 :       if (!s) return zeropol(varn(gel(y,2)));
    3008          14 :       if (s == 1) return gcopy(y);
    3009          14 :       if (s == -1) return gneg(y);
    3010          14 :       return mul_rfrac_scal(gel(y,1), gel(y,2), stoi(s));
    3011             : 
    3012             :     case t_VEC: case t_COL: case t_MAT:
    3013      168113 :       z = cgetg_copy(y, &ly);
    3014      168113 :       for (i=1; i<ly; i++) gel(z,i) = gmulsg(s,gel(y,i));
    3015      168113 :       return z;
    3016             :   }
    3017           0 :   pari_err_TYPE("gmulsg",y);
    3018           0 :   return NULL; /* not reached */
    3019             : }
    3020             : 
    3021             : /********************************************************************/
    3022             : /**                                                                **/
    3023             : /**                       SIMPLE DIVISION                          **/
    3024             : /**                                                                **/
    3025             : /********************************************************************/
    3026             : 
    3027             : GEN
    3028     6323741 : gdivgs(GEN x, long s)
    3029             : {
    3030     6323741 :   long tx = typ(x), lx, i;
    3031             :   pari_sp av;
    3032             :   GEN z, y, p1;
    3033             : 
    3034     6323741 :   if (!s)
    3035             :   {
    3036           0 :     if (is_matvec_t(tx) && lg(x) == 1) return gcopy(x);
    3037           0 :     pari_err_INV("gdivgs",gen_0);
    3038             :   }
    3039     6323741 :   switch(tx)
    3040             :   {
    3041             :     case t_INT:
    3042     2600926 :       av = avma; z = divis_rem(x,s,&i);
    3043     2600926 :       if (!i) return z;
    3044             : 
    3045     1459548 :       i = cgcd(s, i);
    3046     1459548 :       avma=av; z = cgetg(3,t_FRAC);
    3047     1459548 :       if (i == 1) y = icopy(x); else { s /= i; y = diviuexact(x, i); }
    3048     1459548 :       gel(z,1) = y;
    3049     1459548 :       gel(z,2) = stoi(s); normalize_frac(z); return z;
    3050             : 
    3051             :     case t_REAL:
    3052     2505177 :       return divrs(x,s);
    3053             : 
    3054             :     case t_INTMOD:
    3055        7350 :       z = cgetg(3, t_INTMOD);
    3056        7350 :       return div_intmod_same(z, gel(x,1), gel(x,2), modsi(s, gel(x,1)));
    3057             : 
    3058         399 :     case t_FFELT: return FF_Z_Z_muldiv(x,gen_1,stoi(s));
    3059             : 
    3060      218389 :     case t_FRAC: z = cgetg(3, t_FRAC);
    3061      218389 :       i = cgcd(s, smodis(gel(x,1), s));
    3062      218389 :       if (i == 1)
    3063             :       {
    3064       41411 :         gel(z,2) = mulsi(s, gel(x,2));
    3065       41411 :         gel(z,1) = icopy(gel(x,1));
    3066             :       }
    3067             :       else
    3068             :       {
    3069      176978 :         gel(z,2) = mulsi(s/i, gel(x,2));
    3070      176978 :         gel(z,1) = divis(gel(x,1), i);
    3071             :       }
    3072      218389 :       normalize_frac(z);
    3073      218389 :       fix_frac_if_int(z); return z;
    3074             : 
    3075      889300 :     case t_COMPLEX: z = cgetg(3, t_COMPLEX);
    3076      889300 :       gel(z,1) = gdivgs(gel(x,1),s);
    3077      889300 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3078             : 
    3079             :     case t_PADIC: /* divpT */
    3080             :     {
    3081       82915 :       GEN p = gel(x,2);
    3082       82915 :       if (!signe(gel(x,4))) return zeropadic(p, valp(x) - u_pval(s,p));
    3083       82194 :       av = avma;
    3084       82194 :       return gerepileupto(av, divpp(x, cvtop2(stoi(s),x)));
    3085             :     }
    3086             : 
    3087           0 :     case t_QUAD: z = cgetg(4, t_QUAD);
    3088           0 :       gel(z,1) = ZX_copy(gel(x,1));
    3089           0 :       gel(z,2) = gdivgs(gel(x,2),s);
    3090           0 :       gel(z,3) = gdivgs(gel(x,3),s); return z;
    3091             : 
    3092        1505 :     case t_POLMOD: z = cgetg(3, t_POLMOD);
    3093        1505 :       gel(z,1) = RgX_copy(gel(x,1));
    3094        1505 :       gel(z,2) = gdivgs(gel(x,2),s); return z;
    3095             : 
    3096             :     case t_RFRAC:
    3097          91 :       av = avma;
    3098          91 :       p1 = ggcd(stoi(s),gel(x,1));
    3099          91 :       if (typ(p1) == t_INT)
    3100             :       {
    3101          42 :         avma = av;
    3102          42 :         z = cgetg(3, t_RFRAC);
    3103          42 :         i = p1[2];
    3104          42 :         if (i == 1)
    3105             :         {
    3106          42 :           gel(z,1) = gcopy(gel(x,1));
    3107          42 :           gel(z,2) = gmulsg(s,gel(x,2));
    3108             :         }
    3109             :         else
    3110             :         {
    3111           0 :           gel(z,1) = gdivgs(gel(x,1), i);
    3112           0 :           gel(z,2) = gmulgs(gel(x,2), s/i);
    3113             :         }
    3114             :       }
    3115             :       else /* t_FRAC */
    3116             :       {
    3117          49 :         z = cgetg(3, t_RFRAC);
    3118          49 :         gel(z,1) = gdiv(gel(x,1), p1);
    3119          49 :         gel(z,2) = RgX_Rg_mul(gel(x,2), gdivsg(s,p1));
    3120          49 :         z = gerepilecopy(av, z);
    3121             :       }
    3122          91 :       return z;
    3123             : 
    3124             :     case t_POL: case t_SER:
    3125       17675 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3126       17675 :       for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3127       17675 :       return z;
    3128             :     case t_VEC: case t_COL: case t_MAT:
    3129          14 :       z = cgetg_copy(x, &lx);
    3130          14 :       for (i=1; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),s);
    3131          14 :       return z;
    3132             : 
    3133             :   }
    3134           0 :   pari_err_TYPE2("/",x, stoi(s));
    3135           0 :   return NULL; /* not reached */
    3136             : }
    3137             : 
    3138             : /* True shift (exact multiplication by 2^n) */
    3139             : GEN
    3140    36633498 : gmul2n(GEN x, long n)
    3141             : {
    3142             :   long lx, i, k, l;
    3143             :   GEN z, a, b;
    3144             : 
    3145    36633498 :   switch(typ(x))
    3146             :   {
    3147             :     case t_INT:
    3148    27099597 :       if (n>=0) return shifti(x,n);
    3149      163793 :       if (!signe(x)) return gen_0;
    3150       66507 :       l = vali(x); n = -n;
    3151       66507 :       if (n<=l) return shifti(x,-n);
    3152       40019 :       z = cgetg(3,t_FRAC);
    3153       40019 :       gel(z,1) = shifti(x,-l);
    3154       40019 :       gel(z,2) = int2n(n-l); return z;
    3155             : 
    3156             :     case t_REAL:
    3157     4594902 :       return shiftr(x,n);
    3158             : 
    3159       78846 :     case t_INTMOD: b = gel(x,1); a = gel(x,2);
    3160       78846 :       z = cgetg(3,t_INTMOD);
    3161       78845 :       if (n <= 0) return div_intmod_same(z, b, a, modii(int2n(-n), b));
    3162       78838 :       gel(z,2) = gerepileuptoint((pari_sp)z, modii(shifti(a,n), b));
    3163       78842 :       gel(z,1) = icopy(b); return z;
    3164             : 
    3165      231463 :     case t_FFELT: return FF_mul2n(x,n);
    3166             : 
    3167      437076 :     case t_FRAC: a = gel(x,1); b = gel(x,2);
    3168      437076 :       l = vali(a);
    3169      437076 :       k = vali(b);
    3170      437076 :       if (n+l >= k)
    3171             :       {
    3172       68363 :         if (expi(b) == k) return shifti(a,n-k); /* b power of 2 */
    3173       49904 :         l = n-k; k = -k;
    3174             :       }
    3175             :       else
    3176             :       {
    3177      368713 :         k = -(l+n); l = -l;
    3178             :       }
    3179      418617 :       z = cgetg(3,t_FRAC);
    3180      418617 :       gel(z,1) = shifti(a,l);
    3181      418617 :       gel(z,2) = shifti(b,k); return z;
    3182             : 
    3183     3010165 :     case t_COMPLEX: z = cgetg(3,t_COMPLEX);
    3184     3010165 :       gel(z,1) = gmul2n(gel(x,1),n);
    3185     3010165 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3186             : 
    3187          14 :     case t_QUAD: z = cgetg(4,t_QUAD);
    3188          14 :       gel(z,1) = ZX_copy(gel(x,1));
    3189          14 :       gel(z,2) = gmul2n(gel(x,2),n);
    3190          14 :       gel(z,3) = gmul2n(gel(x,3),n); return z;
    3191             : 
    3192       14202 :     case t_POLMOD: z = cgetg(3,t_POLMOD);
    3193       14202 :       gel(z,1) = RgX_copy(gel(x,1));
    3194       14202 :       gel(z,2) = gmul2n(gel(x,2),n); return z;
    3195             : 
    3196             :     case t_POL:
    3197      264809 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3198      264809 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3199      264809 :       return normalizepol_lg(z, lx); /* needed if char = 2 */
    3200             :     case t_SER:
    3201       57799 :       if (ser_isexactzero(x)) return gcopy(x);
    3202       57785 :       z = cgetg_copy(x, &lx); z[1] = x[1];
    3203       57785 :       for (i=2; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3204       57785 :       return normalize(z); /* needed if char = 2 */
    3205             :     case t_VEC: case t_COL: case t_MAT:
    3206      829582 :       z = cgetg_copy(x, &lx);
    3207      829582 :       for (i=1; i<lx; i++) gel(z,i) = gmul2n(gel(x,i),n);
    3208      829582 :       return z;
    3209             : 
    3210             :     case t_RFRAC: /* int2n wrong if n < 0 */
    3211          14 :       return mul_rfrac_scal(gel(x,1),gel(x,2), gmul2n(gen_1,n));
    3212             : 
    3213             :     case t_PADIC: /* int2n wrong if n < 0 */
    3214       15029 :       return gmul(gmul2n(gen_1,n),x);
    3215             :   }
    3216           0 :   pari_err_TYPE("gmul2n",x);
    3217           0 :   return NULL; /* not reached */
    3218             : }
    3219             : 
    3220             : /*******************************************************************/
    3221             : /*                                                                 */
    3222             : /*                              INVERSE                            */
    3223             : /*                                                                 */
    3224             : /*******************************************************************/
    3225             : static GEN
    3226       21343 : inv_polmod(GEN T, GEN x)
    3227             : {
    3228       21343 :   GEN z = cgetg(3,t_POLMOD), a;
    3229       21343 :   gel(z,1) = RgX_copy(T);
    3230       21343 :   if (typ(x) != t_POL || varn(x) != varn(T) || lg(x) <= 3)
    3231        7335 :     a = ginv(x);
    3232             :   else
    3233             :   {
    3234       14008 :     pari_sp av = avma;
    3235       14008 :     if (lg(T) == 5) /* quadratic fields */
    3236       10598 :       a = RgX_Rg_div(quad_polmod_conj(x,T), quad_polmod_norm(x,T));
    3237             :     else
    3238             :     {
    3239        3410 :       GEN p = NULL;
    3240        3410 :       if (RgX_is_FpX(T, &p) && RgX_is_FpX(x, &p) && p)
    3241             :       {
    3242          14 :         T = RgX_to_FpX(T, p); x = RgX_to_FpX(x, p);
    3243          14 :         if (lgefint(p) == 3)
    3244             :         {
    3245           7 :           ulong pp = p[2];
    3246           7 :           a = Flxq_inv(ZX_to_Flx(x, pp), ZX_to_Flx(T, pp), pp);
    3247           7 :           a = Flx_to_ZX(a);
    3248             :         }
    3249             :         else
    3250           7 :           a = FpXQ_inv(x, T, p);
    3251          14 :         a = FpX_to_mod(a, p);
    3252          14 :         a = gerepileupto(av, a);
    3253             :       }
    3254        3396 :       else if (RgX_is_ZX(T) && RgX_is_QX(x))
    3255             :       {
    3256        2670 :         avma = av;
    3257        2670 :         a = QXQ_inv(x, T);
    3258             :       }
    3259             :       else {
    3260         726 :         avma = av;
    3261         726 :         a = RgXQ_inv(x, T);
    3262             :       }
    3263             :     }
    3264             :   }
    3265       21343 :   gel(z,2) = a; return z;
    3266             : }
    3267             : 
    3268             : GEN
    3269    12712716 : ginv(GEN x)
    3270             : {
    3271             :   long s;
    3272             :   pari_sp av, tetpil;
    3273             :   GEN z, y, p1, p2;
    3274             : 
    3275    12712716 :   switch(typ(x))
    3276             :   {
    3277             :     case t_INT:
    3278      242770 :       if (is_pm1(x)) return icopy(x);
    3279      208848 :       s = signe(x); if (!s) pari_err_INV("ginv",gen_0);
    3280      208848 :       z = cgetg(3,t_FRAC);
    3281      208848 :       gel(z,1) = s<0? gen_m1: gen_1;
    3282      208848 :       gel(z,2) = absi(x); return z;
    3283             : 
    3284     1721533 :     case t_REAL: return invr(x);
    3285             : 
    3286         805 :     case t_INTMOD: z=cgetg(3,t_INTMOD);
    3287         805 :       gel(z,1) = icopy(gel(x,1));
    3288         805 :       gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
    3289             : 
    3290             :     case t_FRAC: {
    3291      184113 :       GEN a = gel(x,1), b = gel(x,2);
    3292      184113 :       s = signe(a);
    3293      184113 :       if (is_pm1(a)) return s > 0? icopy(b): negi(b);
    3294       48564 :       z = cgetg(3,t_FRAC);
    3295       48564 :       gel(z,1) = icopy(b);
    3296       48564 :       gel(z,2) = icopy(a);
    3297       48564 :       normalize_frac(z); return z;
    3298             :     }
    3299             :     case t_COMPLEX:
    3300     1180179 :       av=avma;
    3301     1180179 :       p1=cxnorm(x);
    3302     1180179 :       p2=mkcomplex(gel(x,1), gneg(gel(x,2)));
    3303     1180179 :       tetpil=avma;
    3304     1180179 :       return gerepile(av,tetpil,divcR(p2,p1));
    3305             : 
    3306             :     case t_QUAD:
    3307         203 :       av=avma; p1=gnorm(x); p2=gconj(x); tetpil=avma;
    3308         203 :       return gerepile(av,tetpil,gdiv(p2,p1));
    3309             : 
    3310        3094 :     case t_PADIC: z = cgetg(5,t_PADIC);
    3311        3094 :       if (!signe(gel(x,4))) pari_err_INV("ginv",x);
    3312        3087 :       z[1] = _evalprecp(precp(x)) | evalvalp(-valp(x));
    3313        3087 :       gel(z,2) = icopy(gel(x,2));
    3314        3087 :       gel(z,3) = icopy(gel(x,3));
    3315        3087 :       gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
    3316             : 
    3317       21343 :     case t_POLMOD: return inv_polmod(gel(x,1), gel(x,2));
    3318        4705 :     case t_FFELT: return FF_inv(x);
    3319     5436828 :     case t_POL: return gred_rfrac_simple(gen_1,x);
    3320        7777 :     case t_SER: return gdiv(gen_1,x);
    3321             : 
    3322             :     case t_RFRAC:
    3323             :     {
    3324        2905 :       GEN n = gel(x,1), d = gel(x,2);
    3325        2905 :       pari_sp av = avma, ltop;
    3326        2905 :       if (gequal0(n)) pari_err_INV("ginv",x);
    3327             : 
    3328        2905 :       n = simplify_shallow(n);
    3329        2905 :       if (typ(n) != t_POL || varn(n) != varn(d))
    3330             :       {
    3331        2905 :         if (gequal1(n)) { avma = av; return RgX_copy(d); }
    3332         679 :         ltop = avma;
    3333         679 :         z = RgX_Rg_div(d,n);
    3334             :       } else {
    3335           0 :         ltop = avma;
    3336           0 :         z = cgetg(3,t_RFRAC);
    3337           0 :         gel(z,1) = RgX_copy(d);
    3338           0 :         gel(z,2) = RgX_copy(n);
    3339             :       }
    3340         679 :       stackdummy(av, ltop);
    3341         679 :       return z;
    3342             :     }
    3343             : 
    3344             :     case t_QFR:
    3345           7 :       av = avma; z = cgetg(5, t_QFR);
    3346           7 :       gel(z,1) = gel(x,1);
    3347           7 :       gel(z,2) = negi( gel(x,2) );
    3348           7 :       gel(z,3) = gel(x,3);
    3349           7 :       gel(z,4) = negr( gel(x,4) );
    3350           7 :       return gerepileupto(av, redreal(z));
    3351             : 
    3352             :     case t_QFI:
    3353     3905523 :       y = gcopy(x);
    3354     3905523 :       if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
    3355     3810776 :         togglesign(gel(y,2));
    3356     3905523 :       return y;
    3357             :     case t_MAT:
    3358         931 :       y = RgM_inv(x);
    3359         931 :       if (!y) pari_err_INV("ginv",x);
    3360         931 :       return y;
    3361             :     case t_VECSMALL:
    3362             :     {
    3363           0 :       long i, lx = lg(x)-1;
    3364           0 :       y = zero_zv(lx);
    3365           0 :       for (i=1; i<=lx; i++)
    3366             :       {
    3367           0 :         long xi = x[i];
    3368           0 :         if (xi<1 || xi>lx || y[xi])
    3369           0 :           pari_err_TYPE("ginv [not a permutation]", x);
    3370           0 :         y[xi] = i;
    3371             :       }
    3372           0 :       return y;
    3373             :     }
    3374             :   }
    3375           0 :   pari_err_TYPE("inverse",x);
    3376           0 :   return NULL; /* not reached */
    3377             : }

Generated by: LCOV version 1.11