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-bordeaux1.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 - polarit2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16746-c2cb716) Lines: 1504 1729 87.0 %
Date: 2014-08-31 Functions: 123 133 92.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1153 1601 72.0 %

           Branch data     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                 :            : /**               ARITHMETIC OPERATIONS ON POLYNOMIALS                **/
      17                 :            : /**                         (second part)                             **/
      18                 :            : /**                                                                   **/
      19                 :            : /***********************************************************************/
      20                 :            : #include "pari.h"
      21                 :            : #include "paripriv.h"
      22                 :            : 
      23                 :            : #define addshift(x,y) addshiftpol((x),(y),1)
      24                 :            : 
      25                 :            : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
      26                 :            :  * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
      27                 :            :  * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
      28                 :            :  * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
      29                 :            :  * Not memory clean in the latter case */
      30                 :            : GEN
      31                 :       6795 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      32                 :            : {
      33                 :       6795 :   long dP=degpol(P), i, k, m;
      34                 :            :   pari_sp av1, av2;
      35                 :            :   GEN s,y,P_lead;
      36                 :            : 
      37         [ -  + ]:       6795 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      38         [ -  + ]:       6795 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      39         [ -  + ]:       6795 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      40                 :       6795 :   y = cgetg(n+2,t_COL);
      41         [ +  + ]:       6795 :   if (y0)
      42                 :            :   {
      43         [ -  + ]:       3000 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      44                 :       3000 :     m = lg(y0)-1;
      45         [ +  + ]:      13275 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      46                 :            :   }
      47                 :            :   else
      48                 :            :   {
      49                 :       3795 :     m = 1;
      50                 :       3795 :     gel(y,1) = stoi(dP);
      51                 :            :   }
      52                 :       6795 :   P += 2; /* strip codewords */
      53                 :            : 
      54         [ +  + ]:       6795 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      55         [ +  + ]:       6795 :   if (P_lead)
      56                 :            :   {
      57         [ -  + ]:          5 :     if (N) P_lead = Fq_inv(P_lead,T,N);
      58         [ -  + ]:          5 :     else if (T) P_lead = QXQ_inv(P_lead,T);
      59                 :            :   }
      60         [ +  + ]:      20630 :   for (k=m; k<=n; k++)
      61                 :            :   {
      62         [ +  + ]:      13835 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      63 [ +  + ][ +  + ]:      83825 :     for (i=1; i<k && i<=dP; i++)
      64                 :      69990 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      65         [ +  + ]:      13835 :     if (N)
      66                 :            :     {
      67                 :       4025 :       s = Fq_red(s, T, N);
      68         [ -  + ]:       4025 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      69                 :            :     }
      70         [ +  + ]:       9810 :     else if (T)
      71                 :            :     {
      72                 :        410 :       s = grem(s, T);
      73         [ -  + ]:        410 :       if (P_lead) s = grem(gmul(s, P_lead), T);
      74                 :            :     }
      75                 :            :     else
      76         [ +  + ]:       9400 :       if (P_lead) s = gdiv(s, P_lead);
      77                 :      13835 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      78                 :            :   }
      79                 :       6795 :   return y;
      80                 :            : }
      81                 :            : 
      82                 :            : GEN
      83                 :       2645 : polsym(GEN x, long n)
      84                 :            : {
      85                 :       2645 :   return polsym_gen(x, NULL, n, NULL,NULL);
      86                 :            : }
      87                 :            : 
      88                 :            : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
      89                 :            : GEN
      90                 :   66743711 : centermodii(GEN x, GEN p, GEN po2)
      91                 :            : {
      92                 :   66743711 :   GEN y = remii(x, p);
      93   [ +  +  +  - ]:   66743711 :   switch(signe(y))
      94                 :            :   {
      95                 :   17593074 :     case 0: break;
      96 [ +  + ][ +  + ]:   31654510 :     case 1: if (po2 && absi_cmp(y,po2) > 0) y = subii(y, p);
      97                 :   31654510 :       break;
      98 [ +  + ][ +  + ]:   17496127 :     case -1: if (!po2 || absi_cmp(y,po2) > 0) y = addii(y, p);
      99                 :   17496127 :       break;
     100                 :            :   }
     101                 :   66743711 :   return y;
     102                 :            : }
     103                 :            : 
     104                 :            : static long
     105                 :          0 : s_centermod(long x, ulong pp, ulong pps2)
     106                 :            : {
     107                 :          0 :   long y = x % (long)pp;
     108         [ #  # ]:          0 :   if (y < 0) y += pp;
     109                 :          0 :   return Fl_center(y, pp,pps2);
     110                 :            : }
     111                 :            : 
     112                 :            : /* for internal use */
     113                 :            : GEN
     114                 :    4490989 : centermod_i(GEN x, GEN p, GEN ps2)
     115                 :            : {
     116                 :            :   long i, lx;
     117                 :            :   pari_sp av;
     118                 :            :   GEN y;
     119                 :            : 
     120         [ +  + ]:    4490989 :   if (!ps2) ps2 = shifti(p,-1);
     121   [ +  +  +  +  :    4490989 :   switch(typ(x))
                   -  - ]
     122                 :            :   {
     123                 :    2590873 :     case t_INT: return centermodii(x,p,ps2);
     124                 :            : 
     125                 :    1373656 :     case t_POL: lx = lg(x);
     126                 :    1373656 :       y = cgetg(lx,t_POL); y[1] = x[1];
     127         [ +  + ]:    9768704 :       for (i=2; i<lx; i++)
     128                 :            :       {
     129                 :    8395048 :         av = avma;
     130                 :    8395048 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     131                 :            :       }
     132                 :    1373656 :       return normalizepol_lg(y, lx);
     133                 :            : 
     134                 :     526300 :     case t_COL: lx = lg(x);
     135                 :     526300 :       y = cgetg(lx,t_COL);
     136         [ +  + ]:    2486782 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     137                 :     526300 :       return y;
     138                 :            : 
     139                 :        160 :     case t_MAT: lx = lg(x);
     140                 :        160 :       y = cgetg(lx,t_MAT);
     141         [ +  + ]:       3150 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     142                 :        160 :       return y;
     143                 :            : 
     144                 :          0 :     case t_VECSMALL: lx = lg(x);
     145                 :            :     {
     146                 :          0 :       ulong pp = itou(p), pps2 = itou(ps2);
     147                 :          0 :       y = cgetg(lx,t_VECSMALL);
     148         [ #  # ]:          0 :       for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
     149                 :          0 :       return y;
     150                 :            :     }
     151                 :            :   }
     152                 :    4490989 :   return x;
     153                 :            : }
     154                 :            : 
     155                 :            : GEN
     156                 :    3336744 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
     157                 :            : 
     158                 :            : /***********************************************************************/
     159                 :            : /**                                                                   **/
     160                 :            : /**                          FACTORIZATION                            **/
     161                 :            : /**                                                                   **/
     162                 :            : /***********************************************************************/
     163                 :            : #define assign_or_fail(x,y) { GEN __x = x;\
     164                 :            :   if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
     165                 :            : }
     166                 :            : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
     167                 :            : 
     168                 :            : static const long tsh = 6;
     169                 :            : static long
     170                 :         60 : code(long t1, long t2) { return (t1 << tsh) | t2; }
     171                 :            : void
     172                 :         65 : RgX_type_decode(long x, long *t1, long *t2)
     173                 :            : {
     174                 :         65 :   *t1 = x >> tsh;
     175                 :         65 :   *t2 = (x & ((1L<<tsh)-1));
     176                 :         65 : }
     177                 :            : int
     178                 :     112585 : RgX_type_is_composite(long t) { return t >= tsh; }
     179                 :            : 
     180                 :            : long
     181                 :     112900 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     182                 :            : {
     183                 :     112900 :   long t[] = {0,0,0,0,0,0,0,0,0,0};
     184                 :     112900 :   long tx = typ(x), lx, i, j, t2 = 0;
     185                 :     112900 :   GEN ff = NULL;
     186                 :     112900 :   *p = *pol = NULL; *pa = LONG_MAX;
     187         [ -  + ]:     112900 :   if (is_scalar_t(tx))
     188                 :            :   {
     189         [ #  # ]:          0 :     if (tx == t_POLMOD) return 0;
     190                 :          0 :     x = scalarpol(x,0);
     191                 :            :   }
     192                 :            :   /* t[0..1] unused. Other values, if set, indicate a coefficient of type
     193                 :            :    * t[2] : t_REAL
     194                 :            :    * t[3] : t_INTMOD
     195                 :            :    * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
     196                 :            :    * t[5] : t_FFELT
     197                 :            :    * t[6] : t_COMPLEX of t_REAL
     198                 :            :    * t[7] : t_PADIC
     199                 :            :    * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
     200                 :            :    * t[9]: t_POLMOD of rationals (t_INT/t_FRAC) */
     201                 :            :   /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
     202                 :            :    * given by t) */
     203                 :     112900 :   lx = lg(x);
     204         [ +  + ]:     414260 :   for (i=2; i<lx; i++)
     205                 :            :   {
     206                 :     301405 :     GEN c = gel(x,i);
     207   [ +  +  +  +  :     301405 :     switch(typ(c))
             +  +  +  +  
                      + ]
     208                 :            :     {
     209                 :            :       case t_INT: case t_FRAC:
     210                 :     300045 :         break;
     211                 :            :       case t_REAL:
     212         [ +  + ]:        750 :         update_prec(precision(c), pa);
     213                 :        750 :         t[2]=1; break;
     214                 :            :       case t_INTMOD:
     215 [ +  + ][ -  + ]:        165 :         assign_or_fail(gel(c,1),p);
     216                 :        165 :         t[3]=1; break;
     217                 :            :       case t_FFELT:
     218 [ +  + ][ -  + ]:        305 :         if (!ff) ff=c; else if (!FF_samefield(c,ff)) return 0;
     219 [ +  + ][ -  + ]:        305 :         assign_or_fail(FF_p_i(c),p);
     220                 :        305 :         t[5]=1; break;
     221                 :            :       case t_COMPLEX:
     222         [ +  + ]:         65 :         for (j=1; j<=2; j++)
     223                 :            :         {
     224                 :         45 :           GEN d = gel(c,j);
     225   [ +  +  +  +  :         45 :           switch(typ(d))
                      - ]
     226                 :            :           {
     227                 :            :             case t_INT: case t_FRAC:
     228                 :         25 :               t[4]=1; break;
     229                 :            :             case t_REAL:
     230         [ +  - ]:          5 :               update_prec(precision(d), pa);
     231                 :          5 :               t[6]=1; break;
     232                 :            :             case t_INTMOD:
     233 [ +  - ][ #  # ]:         10 :               assign_or_fail(gel(d,1),p);
     234 [ +  - ][ +  + ]:         10 :               if (!signe(*p) || mod4(*p) != 3) return 0;
     235         [ +  - ]:          5 :               if (!t2) t2 = t_COMPLEX;
     236                 :          5 :               t[3]=1; break;
     237                 :            :             case t_PADIC:
     238         [ +  - ]:          5 :               update_prec(precp(d)+valp(d), pa);
     239 [ +  - ][ #  # ]:          5 :               assign_or_fail(gel(d,2),p);
     240         [ +  - ]:          5 :               if (!t2) t2 = t_COMPLEX;
     241                 :          5 :               t[7]=1; break;
     242                 :          0 :             default: return 0;
     243                 :            :           }
     244                 :            :         }
     245 [ +  + ][ +  - ]:         20 :         if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
                 [ #  # ]
     246                 :         20 :         break;
     247                 :            :       case t_PADIC:
     248         [ +  - ]:         20 :         update_prec(precp(c)+valp(c), pa);
     249 [ +  - ][ #  # ]:         20 :         assign_or_fail(gel(c,2),p);
     250                 :         20 :         t[7]=1; break;
     251                 :            :       case t_QUAD:
     252 [ +  - ][ #  # ]:         20 :         assign_or_fail(gel(c,1),pol);
     253         [ +  + ]:         60 :         for (j=2; j<=3; j++)
     254                 :            :         {
     255                 :         40 :           GEN d = gel(c,j);
     256   [ +  +  +  - ]:         40 :           switch(typ(d))
     257                 :            :           {
     258                 :            :             case t_INT: case t_FRAC:
     259                 :         15 :               t[8]=1; break;
     260                 :            :             case t_INTMOD:
     261 [ +  + ][ -  + ]:         20 :               assign_or_fail(gel(d,1),p);
     262         [ +  - ]:         20 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     263                 :         20 :               t[3]=1; break;
     264                 :            :             case t_PADIC:
     265         [ +  - ]:          5 :               update_prec(precp(d)+valp(d), pa);
     266 [ +  - ][ #  # ]:          5 :               assign_or_fail(gel(d,2),p);
     267         [ +  - ]:          5 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     268                 :          5 :               t[7]=1; break;
     269                 :          0 :             default: return 0;
     270                 :            :           }
     271                 :            :         }
     272                 :         20 :         break;
     273                 :            :       case t_POLMOD:
     274 [ +  + ][ -  + ]:         35 :         assign_or_fail(gel(c,1),pol);
     275         [ +  + ]:        105 :         for (j=1; j<=2; j++)
     276                 :            :         {
     277                 :         70 :           GEN pbis = NULL, polbis = NULL;
     278                 :            :           long pabis;
     279   [ +  +  +  - ]:         70 :           switch(RgX_type(gel(c,j),&pbis,&polbis,&pabis))
     280                 :            :           {
     281                 :         55 :             case t_INT: t[9]=1; break;
     282                 :         10 :             case t_INTMOD: t[3]=1; t2 = t_POLMOD; break;
     283         [ +  - ]:          5 :             case t_PADIC: t[7]=1; t2 = t_POLMOD; update_prec(pabis,pa); break;
     284                 :          0 :             default: return 0;
     285                 :            :           }
     286 [ +  + ][ +  - ]:         70 :           if (pbis) assign_or_fail(pbis,p);
                 [ #  # ]
     287 [ -  + ][ #  # ]:         70 :           if (polbis) assign_or_fail(polbis,pol);
                 [ #  # ]
     288                 :            :         }
     289                 :         35 :         break;
     290                 :         40 :       default: return 0;
     291                 :            :     }
     292                 :            :   }
     293         [ +  + ]:     112855 :   if (t[5]) /* ffelt */
     294                 :            :   {
     295 [ +  - ][ +  - ]:         75 :     if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
         [ +  - ][ +  - ]
         [ +  - ][ -  + ]
     296                 :         75 :     *pol=ff; return t_FFELT;
     297                 :            :   }
     298         [ +  + ]:     112780 :   if (t[6]) /* inexact, complex */
     299                 :            :   {
     300 [ +  - ][ +  - ]:          5 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
         [ +  - ][ -  + ]
     301                 :          5 :     return t_COMPLEX;
     302                 :            :   }
     303         [ +  + ]:     112775 :   if (t[2]) /* inexact, real */
     304                 :            :   {
     305 [ +  - ][ +  - ]:        100 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
         [ +  - ][ -  + ]
     306         [ -  + ]:        100 :     return t[4]?t_COMPLEX:t_REAL;
     307                 :            :   }
     308         [ +  + ]:     112675 :   if (t2) /* polmod/quad/complex of intmod/padic */
     309                 :            :   {
     310         [ +  + ]:         40 :     if (t[3]) return code(t2,t_INTMOD);
     311         [ +  - ]:         15 :     if (t[7]) return code(t2,t_PADIC);
     312                 :            :   }
     313         [ +  + ]:     112635 :   if (t[9]) return code(t_POLMOD,t_INT);
     314         [ +  + ]:     112625 :   if (t[8]) return code(t_QUAD,t_INT);
     315         [ +  + ]:     112620 :   if (t[4]) return code(t_COMPLEX,t_INT);
     316         [ +  + ]:     112615 :   if (t[3]) return t_INTMOD;
     317         [ +  + ]:     112555 :   if (t[7]) return t_PADIC;
     318                 :     112900 :   return t_INT;
     319                 :            : }
     320                 :            : 
     321                 :            : GEN
     322                 :          0 : factor0(GEN x,long flag)
     323                 :            : {
     324         [ #  # ]:          0 :   return (flag<0)? factor(x): boundfact(x,flag);
     325                 :            : }
     326                 :            : 
     327                 :            : /* only present for interface with GP */
     328                 :            : GEN
     329                 :      25845 : gp_factor0(GEN x, GEN flag)
     330                 :            : {
     331                 :            :   ulong B;
     332         [ +  + ]:      25845 :   if (!flag) return factor(x);
     333 [ +  - ][ -  + ]:         25 :   if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
     334      [ +  +  - ]:         25 :   switch(lgefint(flag))
     335                 :            :   {
     336                 :          5 :     case 2: B = 0; break;
     337                 :         20 :     case 3: B = flag[2]; break;
     338                 :          0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     339                 :          0 :              return NULL; /*not reached*/
     340                 :            :   }
     341                 :      25815 :   return boundfact(x, B);
     342                 :            : }
     343                 :            : 
     344                 :            : GEN
     345                 :      27420 : deg1_from_roots(GEN L, long v)
     346                 :            : {
     347                 :      27420 :   long i, l = lg(L);
     348                 :      27420 :   GEN z = cgetg(l,t_COL);
     349         [ +  + ]:      55225 :   for (i=1; i<l; i++)
     350                 :      27805 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     351                 :      27420 :   return z;
     352                 :            : }
     353                 :            : GEN
     354                 :        435 : roots_from_deg1(GEN x)
     355                 :            : {
     356                 :        435 :   long i,l = lg(x);
     357                 :        435 :   GEN r = cgetg(l,t_VEC);
     358         [ +  + ]:       6295 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     359                 :        435 :   return r;
     360                 :            : }
     361                 :            : 
     362                 :            : static GEN
     363                 :         25 : gauss_factor_p(GEN p)
     364                 :            : {
     365                 :         25 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     366                 :         25 :   return mkcomplex(a, b);
     367                 :            : }
     368                 :            : 
     369                 :            : static GEN
     370                 :         30 : gauss_primpart(GEN x, GEN *c)
     371                 :            : {
     372                 :         30 :   GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
     373         [ -  + ]:         30 :   *c = n; if (n == gen_1) return x;
     374                 :         30 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     375                 :            : }
     376                 :            : 
     377                 :            : static GEN
     378                 :         50 : gauss_primpart_try(GEN x, GEN c)
     379                 :            : {
     380                 :            :   GEN r, y;
     381         [ +  + ]:         50 :   if (typ(x) == t_INT)
     382                 :            :   {
     383         [ -  + ]:         30 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     384                 :            :   }
     385                 :            :   else
     386                 :            :   {
     387                 :         20 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     388         [ +  + ]:         20 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     389         [ -  + ]:         10 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     390                 :            :   }
     391                 :         50 :   return y;
     392                 :            : }
     393                 :            : 
     394                 :            : static int
     395                 :         55 : gauss_cmp(GEN x, GEN y)
     396                 :            : {
     397                 :            :   int v;
     398         [ -  + ]:         55 :   if (typ(x) != t_COMPLEX)
     399         [ #  # ]:          0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     400         [ +  + ]:         55 :   if (typ(y) != t_COMPLEX) return 1;
     401                 :         35 :   v = cmpii(gel(x,2), gel(y,2));
     402         [ +  + ]:         35 :   if (v) return v;
     403                 :         55 :   return gcmp(gel(x,1), gel(y,1));
     404                 :            : }
     405                 :            : 
     406                 :            : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
     407                 :            : static GEN
     408                 :         80 : gauss_normal(GEN x)
     409                 :            : {
     410 [ +  + ][ +  - ]:         80 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     411         [ -  + ]:         75 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     412         [ +  + ]:         75 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     413                 :         80 :   return x;
     414                 :            : }
     415                 :            : 
     416                 :            : static GEN
     417                 :         30 : gauss_factor(GEN x)
     418                 :            : {
     419                 :         30 :   pari_sp av = avma;
     420                 :         30 :   GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
     421                 :         30 :   long t1 = typ(a);
     422                 :         30 :   long t2 = typ(b), i, j, l, exp = 0;
     423         [ +  + ]:         30 :   if (t1 == t_FRAC) d = gel(a,2);
     424         [ +  + ]:         30 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     425         [ +  + ]:         30 :   if (d == gen_1) y = x;
     426                 :            :   else
     427                 :            :   {
     428                 :         10 :     y = gmul(x, d);
     429                 :         10 :     a = gel(y,1); t1 = typ(a);
     430                 :         10 :     b = gel(y,2); t2 = typ(b);
     431                 :            :   }
     432 [ +  - ][ -  + ]:         30 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     433                 :         30 :   y = gauss_primpart(y, &n);
     434                 :         30 :   fa = factor(cxnorm(y));
     435                 :         30 :   P = gel(fa,1);
     436                 :         30 :   E = gel(fa,2); l = lg(P);
     437                 :         30 :   P2 = cgetg(l, t_COL);
     438                 :         30 :   E2 = cgetg(l, t_COL);
     439         [ +  + ]:         70 :   for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
     440                 :            :   { /* either p = 2 (ramified) or those factors split in Q(i) */
     441                 :         40 :     GEN p = gel(P,i), w, w2, t, we, pe;
     442                 :         40 :     long v, e = itos(gel(E,i));
     443                 :         40 :     int is2 = equaliu(p, 2);
     444         [ +  + ]:         40 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     445                 :         40 :     w2 = gauss_normal( gconj(w) );
     446                 :            :     /* w * w2 * I^3 = p, w2 = gconj(w) * I */
     447                 :         40 :     pe = powiu(p, e);
     448                 :         40 :     we = gpowgs(w, e);
     449                 :         40 :     t = gauss_primpart_try( gmul(y, gconj(we)), pe );
     450         [ +  + ]:         40 :     if (t) y = t; /* y /= w^e */
     451                 :            :     else {
     452                 :            :       /* y /= conj(w)^e, should be y /= w2^e */
     453                 :         10 :       y = gauss_primpart_try( gmul(y, we), pe );
     454                 :         10 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     455                 :            :     }
     456                 :         40 :     gel(P,i) = w;
     457                 :         40 :     v = Z_pvalrem(n, p, &n);
     458         [ +  + ]:         40 :     if (v) {
     459                 :          5 :       exp -= v; /* += 3*v mod 4 */
     460         [ +  - ]:          5 :       if (is2) v <<= 1; /* 2 = w^2 I^3 */
     461                 :            :       else {
     462                 :          0 :         gel(P2,j) = w2;
     463                 :          0 :         gel(E2,j) = utoipos(v); j++;
     464                 :            :       }
     465                 :          5 :       gel(E,i) = stoi(e + v);
     466                 :            :     }
     467                 :         40 :     v = Z_pvalrem(d, p, &d);
     468         [ +  + ]:         40 :     if (v) {
     469                 :          5 :       exp += v; /* -= 3*v mod 4 */
     470         [ -  + ]:          5 :       if (is2) v <<= 1; /* 2 is ramified */
     471                 :            :       else {
     472                 :          5 :         gel(P2,j) = w2;
     473                 :          5 :         gel(E2,j) = utoineg(v); j++;
     474                 :            :       }
     475                 :          5 :       gel(E,i) = stoi(e - v);
     476                 :            :     }
     477                 :         40 :     exp &= 3;
     478                 :            :   }
     479         [ +  + ]:         30 :   if (j > 1) {
     480                 :          5 :     long k = 1;
     481                 :          5 :     GEN P1 = cgetg(l, t_COL);
     482                 :          5 :     GEN E1 = cgetg(l, t_COL);
     483                 :            :     /* remove factors with exponent 0 */
     484         [ +  + ]:         10 :     for (i = 1; i < l; i++)
     485         [ -  + ]:          5 :       if (signe(gel(E,i)))
     486                 :            :       {
     487                 :          0 :         gel(P1,k) = gel(P,i);
     488                 :          0 :         gel(E1,k) = gel(E,i);
     489                 :          0 :         k++;
     490                 :            :       }
     491                 :          5 :     setlg(P1, k); setlg(E1, k);
     492                 :          5 :     setlg(P2, j); setlg(E2, j);
     493                 :          5 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     494                 :            :   }
     495 [ +  + ][ +  + ]:         30 :   if (!is_pm1(n) || !is_pm1(d))
     496                 :            :   {
     497                 :         15 :     GEN Fa = factor(gdiv(n, d));
     498                 :         15 :     P = gel(Fa,1); l = lg(P);
     499                 :         15 :     E = gel(Fa,2);
     500         [ +  + ]:         35 :     for (i = 1; i < l; i++)
     501                 :            :     {
     502                 :         20 :       GEN w, p = gel(P,i);
     503                 :            :       long e;
     504                 :            :       int is2;
     505      [ +  +  + ]:         20 :       switch(mod4(p))
     506                 :            :       {
     507                 :         10 :         case 3: continue;
     508                 :          5 :         case 2: is2 = 1; break;
     509                 :          5 :         default:is2 = 0; break;
     510                 :            :       }
     511                 :         10 :       e = itos(gel(E,i));
     512         [ +  + ]:         10 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     513                 :         10 :       gel(P,i) = w;
     514         [ +  + ]:         10 :       if (is2)
     515                 :          5 :         gel(E,i) = stoi(e << 1);
     516                 :            :       else
     517                 :            :       {
     518                 :          5 :         P = shallowconcat(P, gauss_normal( gconj(w) ));
     519                 :          5 :         E = shallowconcat(E, gel(E,i));
     520                 :            :       }
     521                 :         10 :       exp -= e; /* += 3*e mod 4 */
     522                 :         10 :       exp &= 3;
     523                 :            :     }
     524                 :         15 :     gel(Fa,1) = P;
     525                 :         15 :     gel(Fa,2) = E;
     526                 :         15 :     fa = famat_mul_shallow(fa, Fa);
     527                 :            :   }
     528                 :         30 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     529                 :            : 
     530                 :         30 :   y = gmul(y, powIs(exp));
     531         [ +  + ]:         30 :   if (!gequal1(y)) {
     532                 :         25 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     533                 :         25 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     534                 :            :   }
     535                 :         30 :   return gerepilecopy(av, fa);
     536                 :            : }
     537                 :            : 
     538                 :            : GEN
     539                 :      26850 : factor(GEN x)
     540                 :            : {
     541                 :      26850 :   long tx=typ(x), lx, i, pa, v, r1;
     542                 :            :   pari_sp av, tetpil;
     543                 :            :   GEN  y, p, p1, p2, pol;
     544                 :            : 
     545         [ +  + ]:      26850 :   if (gequal0(x))
     546         [ +  - ]:         15 :     switch(tx)
     547                 :            :     {
     548                 :            :       case t_INT:
     549                 :            :       case t_COMPLEX:
     550                 :            :       case t_POL:
     551                 :         15 :       case t_RFRAC: return prime_fact(x);
     552                 :          0 :       default: pari_err_TYPE("factor",x);
     553                 :            :     }
     554                 :      26835 :   av = avma;
     555   [ +  +  +  +  :      26835 :   switch(tx)
                   +  - ]
     556                 :            :   {
     557                 :      26420 :     case t_INT: return Z_factor(x);
     558                 :            : 
     559                 :            :     case t_FRAC:
     560                 :        205 :       p1 = Z_factor(gel(x,1));
     561                 :        205 :       p2 = Z_factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
     562                 :        205 :       return gerepilecopy(av, merge_factor_i(p1,p2));
     563                 :            : 
     564                 :            :     case t_POL:
     565                 :        175 :       tx=RgX_type(x,&p,&pol,&pa);
     566   [ +  +  +  +  :        175 :       switch(tx)
             +  +  +  + ]
     567                 :            :       {
     568                 :          5 :         case 0: pari_err_IMPL("factor for general polynomials");
     569                 :         55 :         case t_INT: return QX_factor(x);
     570                 :          5 :         case t_INTMOD: return factmod(x,p);
     571                 :            : 
     572                 :          5 :         case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2;
     573                 :          5 :           av = avma; p1 = roots(x,pa); tetpil = avma;
     574                 :          5 :           p1 = deg1_from_roots(p1, varn(x));
     575                 :          5 :           gel(y,1) = gerepile(av,tetpil,p1);
     576                 :          5 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     577                 :            : 
     578                 :         10 :         case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
     579                 :         10 :           av=avma; p1=cleanroots(x,pa); tetpil=avma;
     580         [ +  + ]:         25 :           for(r1=1; r1<lx; r1++)
     581         [ +  + ]:         20 :             if (typ(gel(p1,r1)) == t_COMPLEX) break;
     582                 :         10 :           lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     583         [ +  + ]:         25 :           for(i=1; i<r1; i++)
     584                 :         15 :             gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     585         [ +  + ]:         15 :           for(   ; i<lx; i++)
     586                 :            :           {
     587                 :          5 :             GEN a = gel(p1,2*i-r1);
     588                 :          5 :             p = cgetg(5, t_POL); gel(p2,i) = p;
     589                 :          5 :             p[1] = x[1];
     590                 :          5 :             gel(p,2) = gnorm(a);
     591                 :          5 :             gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     592                 :          5 :             gel(p,4) = gen_1;
     593                 :            :           }
     594                 :         10 :           gel(y,1) = gerepile(av,tetpil,p2);
     595                 :         10 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     596                 :            : 
     597                 :         10 :         case t_PADIC: return factorpadic(x,p,pa);
     598                 :            : 
     599                 :         25 :         case t_FFELT: return FFX_factor(x,pol);
     600                 :            : 
     601                 :            :         default:
     602                 :            :         {
     603                 :            :           GEN w;
     604                 :            :           long killv, t1, t2;
     605                 :         60 :           x = leafcopy(x); lx=lg(x);
     606                 :         60 :           pol = leafcopy(pol);
     607                 :         60 :           v = pari_var_next_temp();
     608         [ +  + ]:        240 :           for(i=2; i<lx; i++)
     609                 :            :           {
     610                 :        180 :             p1 = gel(x,i);
     611      [ +  +  + ]:        180 :             switch(typ(p1))
     612                 :            :             {
     613                 :         20 :               case t_QUAD: p1++;
     614                 :            :               case t_COMPLEX:
     615                 :         35 :                 gel(x,i) = mkpolmod(deg1pol_shallow(gel(p1,2), gel(p1,1), v), pol);
     616                 :            :             }
     617                 :            :           }
     618                 :         60 :           killv = (avma != (pari_sp)pol);
     619         [ +  + ]:         60 :           if (killv) setvarn(pol, fetch_var());
     620                 :         60 :           RgX_type_decode(tx, &t1, &t2);
     621      [ +  +  + ]:         60 :           switch (t2)
     622                 :            :           {
     623                 :         20 :             case t_INT: p1 = polfnf(x,pol); break;
     624                 :            :             case t_INTMOD:
     625                 :         25 :               pol = RgX_to_FpX(pol, p);
     626 [ +  + ][ +  + ]:         25 :               if (FpX_is_squarefree(pol,p) && FpX_nbfact(pol, p) == 1)
     627                 :            :               {
     628                 :         15 :                 p1 = factorff(x,p,pol); break;
     629                 :            :               }
     630                 :            :             /*fall through*/
     631                 :         25 :             default: pari_err_IMPL("factor for general polynomial");
     632                 :          0 :               return NULL; /* not reached */
     633                 :            :           }
     634   [ +  +  +  - ]:         35 :           switch (t1)
     635                 :            :           {
     636                 :            :             case t_POLMOD:
     637         [ -  + ]:         15 :               if (killv) (void)delete_var();
     638                 :         15 :               return gerepileupto(av,p1);
     639                 :         10 :             case t_COMPLEX: w = gen_I(); break;
     640                 :         10 :             case t_QUAD: w = mkquad(pol,gen_0,gen_1);
     641                 :         10 :               break;
     642                 :          0 :             default: pari_err_IMPL("factor for general polynomial");
     643                 :          0 :               return NULL; /* not reached */
     644                 :            :           }
     645                 :         20 :           p2=gel(p1,1);
     646         [ +  + ]:         50 :           for(i=1; i<lg(p2); i++)
     647                 :            :           {
     648                 :         30 :             GEN P = gel(p2,i);
     649                 :            :             long j;
     650         [ +  + ]:        100 :             for(j=2; j<lg(P); j++)
     651                 :            :             {
     652                 :         70 :               GEN p = gel(P,j);
     653         [ +  + ]:         70 :               if(typ(p)==t_POLMOD) gel(P,j) = gsubst(gel(p,2),v,w);
     654                 :            :             }
     655                 :            :           }
     656         [ +  - ]:         20 :           if (killv) (void)delete_var();
     657                 :         35 :           return gerepilecopy(av, p1);
     658                 :            :         }
     659                 :            :       }
     660                 :            :     case t_RFRAC: {
     661                 :          5 :       GEN a = gel(x,1), b = gel(x,2);
     662                 :          5 :       y = factor(b); gel(y,2) = gneg_i(gel(y,2));
     663 [ +  - ][ +  - ]:          5 :       if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a), y);
     664                 :          5 :       return gerepilecopy(av, y);
     665                 :            :     }
     666                 :            : 
     667                 :            :     case t_COMPLEX:
     668                 :         30 :       y = gauss_factor(x);
     669         [ +  - ]:         30 :       if (y) return y;
     670                 :            :       /* fall through */
     671                 :            :   }
     672                 :          0 :   pari_err_TYPE("factor",x);
     673                 :      26820 :   return NULL; /* not reached */
     674                 :            : }
     675                 :            : 
     676                 :            : /*******************************************************************/
     677                 :            : /*                                                                 */
     678                 :            : /*                     ROOTS --> MONIC POLYNOMIAL                  */
     679                 :            : /*                                                                 */
     680                 :            : /*******************************************************************/
     681                 :            : static GEN
     682                 :     121292 : normalized_mul(GEN x, GEN y)
     683                 :            : {
     684                 :     121292 :   long a = gel(x,1)[1], b = gel(y,1)[1];
     685                 :     121292 :   return mkvec2(mkvecsmall(a + b),
     686                 :     242584 :                 RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
     687                 :            : }
     688                 :            : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
     689                 :            : static GEN
     690                 :      24419 : normalized_to_RgX(GEN L)
     691                 :            : {
     692                 :      24419 :   long i, a = gel(L,1)[1];
     693                 :      24419 :   GEN A = gel(L,2);
     694                 :      24419 :   GEN z = cgetg(a + 3, t_POL);
     695                 :      24419 :   z[1] = evalsigne(1) | evalvarn(varn(A));
     696         [ +  + ]:     307820 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
     697         [ +  + ]:      24429 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
     698                 :      24419 :   gel(z,i) = gen_1; return z;
     699                 :            : }
     700                 :            : 
     701                 :            : /* compute prod (x - a[i]) */
     702                 :            : GEN
     703                 :       7279 : roots_to_pol(GEN a, long v)
     704                 :            : {
     705                 :       7279 :   pari_sp av = avma;
     706                 :       7279 :   long i, k, lx = lg(a);
     707                 :            :   GEN L;
     708         [ -  + ]:       7279 :   if (lx == 1) return pol_1(v);
     709                 :       7279 :   L = cgetg(lx, t_VEC);
     710         [ +  + ]:      34326 :   for (k=1,i=1; i<lx-1; i+=2)
     711                 :            :   {
     712                 :      27047 :     GEN s = gel(a,i), t = gel(a,i+1);
     713                 :      27047 :     GEN x0 = gmul(s,t);
     714                 :      27047 :     GEN x1 = gneg(gadd(s,t));
     715                 :      27047 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     716                 :            :   }
     717         [ +  + ]:       7279 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
     718                 :       3346 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
     719                 :       7279 :   setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
     720                 :       7279 :   return gerepileupto(av, normalized_to_RgX(L));
     721                 :            : }
     722                 :            : 
     723                 :            : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
     724                 :            : GEN
     725                 :      17140 : roots_to_pol_r1(GEN a, long v, long r1)
     726                 :            : {
     727                 :      17140 :   pari_sp av = avma;
     728                 :      17140 :   long i, k, lx = lg(a);
     729                 :            :   GEN L;
     730         [ -  + ]:      17140 :   if (lx == 1) return pol_1(v);
     731                 :      17140 :   L = cgetg(lx, t_VEC);
     732         [ +  + ]:      87197 :   for (k=1,i=1; i<r1; i+=2)
     733                 :            :   {
     734                 :      70057 :     GEN s = gel(a,i), t = gel(a,i+1);
     735                 :      70057 :     GEN x0 = gmul(s,t);
     736                 :      70057 :     GEN x1 = gneg(gadd(s,t));
     737                 :      70057 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     738                 :            :   }
     739         [ +  + ]:      17140 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
     740                 :       4665 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
     741         [ +  + ]:      57736 :   for (i=r1+1; i<lx; i++)
     742                 :            :   {
     743                 :      40596 :     GEN s = gel(a,i);
     744                 :      40596 :     GEN x0 = gnorm(s);
     745                 :      40596 :     GEN x1 = gneg(gtrace(s));
     746                 :      40596 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     747                 :            :   }
     748                 :      17140 :   setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
     749                 :      17140 :   return gerepileupto(av, normalized_to_RgX(L));
     750                 :            : }
     751                 :            : 
     752                 :            : GEN
     753                 :     573895 : divide_conquer_assoc(GEN x, void *data, GEN (*mul)(void *,GEN,GEN))
     754                 :            : {
     755                 :            :   pari_sp ltop, lim;
     756                 :     573895 :   long i,k,lx = lg(x);
     757                 :            : 
     758         [ +  + ]:     573895 :   if (lx == 1) return gen_1;
     759         [ +  + ]:     573490 :   if (lx == 2) return gcopy(gel(x,1));
     760                 :     527343 :   x = leafcopy(x); k = lx;
     761                 :     527343 :   ltop=avma; lim = stack_lim(ltop,1);
     762         [ +  + ]:    1484748 :   while (k > 2)
     763                 :            :   {
     764         [ -  + ]:     957405 :     if (DEBUGLEVEL>7)
     765                 :          0 :       err_printf("prod: remaining objects %ld\n",k-1);
     766                 :     957405 :     lx = k; k = 1;
     767         [ +  + ]:    2947358 :     for (i=1; i<lx-1; i+=2)
     768                 :    1989953 :       gel(x,k++) = mul(data,gel(x,i),gel(x,i+1));
     769         [ +  + ]:     957405 :     if (i < lx) gel(x,k++) = gel(x,i);
     770         [ +  + ]:     957405 :     if (low_stack(lim,stack_lim(ltop,1)))
     771                 :          6 :       gerepilecoeffs(ltop,x+1,k-1);
     772                 :            :   }
     773                 :     573895 :   return gel(x,1);
     774                 :            : }
     775                 :            : 
     776                 :            : static GEN
     777                 :    1217687 : _domul(void *data, GEN x, GEN y)
     778                 :            : {
     779                 :    1217687 :   GEN (*mul)(GEN,GEN)=(GEN (*)(GEN,GEN)) data;
     780                 :    1217687 :   return mul(x,y);
     781                 :            : }
     782                 :            : GEN
     783                 :     225675 : divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
     784                 :     225675 : { return divide_conquer_assoc(x, (void *)mul, _domul); }
     785                 :            : 
     786                 :            : /*******************************************************************/
     787                 :            : /*                                                                 */
     788                 :            : /*                          FACTORBACK                             */
     789                 :            : /*                                                                 */
     790                 :            : /*******************************************************************/
     791                 :            : static GEN
     792                 :          0 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
     793                 :            : static GEN
     794                 :         40 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
     795                 :            : static GEN
     796                 :          5 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
     797                 :            : static GEN
     798                 :         70 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
     799                 :            : static GEN
     800                 :       1511 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
     801                 :            : static GEN
     802                 :       2166 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
     803                 :            : static GEN
     804                 :     737115 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
     805                 :            : static GEN
     806                 :    1071945 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
     807                 :            : 
     808                 :            : #if 0
     809                 :            : static GEN
     810                 :            : _ellmul(void *ell, GEN x, GEN y) { return elladd((GEN) ell, x, y); }
     811                 :            : static GEN
     812                 :            : _ellpow(void *ell, GEN x, GEN n) { return ellmul((GEN) ell, x, n); }
     813                 :            : #endif
     814                 :            : 
     815                 :            : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
     816                 :            : GEN
     817                 :     336090 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
     818                 :            :                GEN (*_pow)(void*,GEN,GEN), void *data)
     819                 :            : {
     820                 :     336090 :   pari_sp av = avma;
     821                 :            :   long k, l, lx;
     822                 :            :   GEN p,x;
     823                 :            : 
     824         [ +  + ]:     336090 :   if (e) /* supplied vector of exponents */
     825                 :     324605 :     p = L;
     826                 :            :   else
     827                 :            :   {
     828      [ +  +  - ]:      11485 :     switch(typ(L)) {
     829                 :            :       case t_VEC:
     830                 :            :       case t_COL: /* product of the L[i] */
     831                 :         35 :         return gerepileupto(av, divide_conquer_assoc(L, data, _mul));
     832                 :            :       case t_MAT: /* genuine factorization */
     833                 :      11450 :         l = lg(L);
     834         [ -  + ]:      11450 :         if (l == 1) return gen_1;
     835         [ +  + ]:      11450 :         if (l == 3) break;
     836                 :            :         /*fall through*/
     837                 :            :       default:
     838                 :          5 :         pari_err_TYPE("factorback [not a factorization]", L);
     839                 :            :     }
     840                 :      11445 :     p = gel(L,1);
     841                 :      11445 :     e = gel(L,2);
     842                 :            :   }
     843                 :            :   /* p = elts, e = expo */
     844                 :     336050 :   lx = lg(p);
     845                 :            :   /* check whether e is an integral vector of correct length */
     846 [ +  - ][ +  - ]:     336050 :   if (!is_vec_t(typ(e)) || lx != lg(e) || !RgV_is_ZV(e))
                 [ +  + ]
     847                 :          5 :     pari_err_TYPE("factorback [not an exponent vector]", e);
     848         [ +  + ]:     336045 :   if (lx == 1) return gen_1;
     849                 :     336025 :   x = cgetg(lx,t_VEC);
     850         [ +  + ]:    1412425 :   for (l=1,k=1; k<lx; k++)
     851         [ +  + ]:    1076400 :     if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
     852                 :     336025 :   x[0] = evaltyp(t_VEC) | _evallg(l);
     853                 :     336080 :   return gerepileupto(av, divide_conquer_assoc(x, data, _mul));
     854                 :            : }
     855                 :            : 
     856                 :            : GEN
     857                 :        105 : idealfactorback(GEN nf, GEN L, GEN e, int red)
     858                 :            : {
     859                 :        105 :   nf = checknf(nf);
     860         [ +  + ]:        105 :   if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
     861                 :        105 :   else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
     862                 :            : }
     863                 :            : 
     864                 :            : GEN
     865                 :        915 : nffactorback(GEN nf, GEN L, GEN e)
     866                 :        915 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
     867                 :            : 
     868                 :            : GEN
     869                 :     335070 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
     870                 :            : GEN
     871                 :      11260 : factorback(GEN fa) { return factorback2(fa, NULL); }
     872                 :            : 
     873                 :            : static int
     874                 :         20 : RgX_is_irred_i(GEN x)
     875                 :            : {
     876                 :            :   GEN y, p, pol;
     877                 :         20 :   long l = lg(x), pa;
     878                 :            : 
     879 [ +  - ][ -  + ]:         20 :   if (!signe(x) || l <= 3) return 0;
     880   [ +  -  -  + ]:         20 :   switch(RgX_type(x,&p,&pol,&pa))
     881                 :            :   {
     882                 :         10 :     case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
     883                 :          0 :     case t_COMPLEX: return l == 4;
     884                 :            :     case t_REAL:
     885         [ #  # ]:          0 :       if (l == 4) return 1;
     886         [ #  # ]:          0 :       if (l > 5) return 0;
     887                 :          0 :       return gsigne(RgX_disc(x)) > 0;
     888                 :            :   }
     889                 :         10 :   y = factor(x);
     890                 :         20 :   return (lg(gcoeff(y,1,1))==l);
     891                 :            : }
     892                 :            : static int
     893                 :         20 : RgX_is_irred(GEN x)
     894                 :            : {
     895                 :         20 :   pari_sp av = avma;
     896                 :         20 :   int r = RgX_is_irred_i(x);
     897                 :         20 :   avma = av; return r;
     898                 :            : }
     899                 :            : long
     900                 :         20 : isirreducible(GEN x)
     901                 :            : {
     902      [ -  +  - ]:         20 :   switch(typ(x))
     903                 :            :   {
     904                 :          0 :     case t_INT: case t_REAL: case t_FRAC: return 0;
     905                 :         20 :     case t_POL: return RgX_is_irred(x);
     906                 :            :   }
     907                 :          0 :   pari_err_TYPE("isirreducible",x);
     908                 :         20 :   return 0;
     909                 :            : }
     910                 :            : 
     911                 :            : /*******************************************************************/
     912                 :            : /*                                                                 */
     913                 :            : /*                         GENERIC GCD                             */
     914                 :            : /*                                                                 */
     915                 :            : /*******************************************************************/
     916                 :            : /* x is a COMPLEX or a QUAD */
     917                 :            : static GEN
     918                 :        485 : triv_cont_gcd(GEN x, GEN y)
     919                 :            : {
     920                 :        485 :   pari_sp av = avma;
     921                 :            :   GEN c;
     922         [ +  + ]:        485 :   if (typ(x)==t_COMPLEX)
     923                 :            :   {
     924                 :         25 :     GEN a = gel(x,1), b = gel(x,2);
     925 [ +  + ][ -  + ]:         25 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
     926                 :         15 :     c = ggcd(a,b);
     927                 :            :   }
     928                 :            :   else
     929                 :        460 :     c = ggcd(gel(x,2),gel(x,3));
     930                 :        485 :   return gerepileupto(av, ggcd(c,y));
     931                 :            : }
     932                 :            : 
     933                 :            : /* y is a PADIC, x a rational number or an INTMOD */
     934                 :            : static GEN
     935                 :         65 : padic_gcd(GEN x, GEN y)
     936                 :            : {
     937                 :         65 :   GEN p = gel(y,2);
     938                 :         65 :   long v = gvaluation(x,p), w = valp(y);
     939         [ -  + ]:         65 :   if (w < v) v = w;
     940                 :         65 :   return powis(p, v);
     941                 :            : }
     942                 :            : 
     943                 :            : /* x,y in Z[i], at least one of which is t_COMPLEX */
     944                 :            : static GEN
     945                 :         35 : gauss_gcd(GEN x, GEN y)
     946                 :            : {
     947                 :         35 :   pari_sp av = avma;
     948                 :            :   GEN dx, dy;
     949                 :         35 :   dx = denom(x); x = gmul(x, dx);
     950                 :         35 :   dy = denom(y); y = gmul(y, dy);
     951         [ +  + ]:         65 :   while (!gequal0(y))
     952                 :            :   {
     953                 :         30 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
     954                 :         30 :     x = y; y = z;
     955                 :            :   }
     956                 :         35 :   x = gauss_normal(x);
     957         [ +  + ]:         35 :   if (typ(x) == t_COMPLEX)
     958                 :            :   {
     959         [ -  + ]:         25 :     if      (gequal0(gel(x,2))) x = gel(x,1);
     960         [ +  + ]:         25 :     else if (gequal0(gel(x,1))) x = gel(x,2);
     961                 :            :   }
     962                 :         35 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
     963                 :            : }
     964                 :            : 
     965                 :            : static int
     966                 :         75 : is_rational(GEN x) { long t = typ(x); return is_rational_t(t); }
     967                 :            : static int
     968                 :         40 : c_is_rational(GEN x)
     969                 :            : {
     970 [ +  + ][ +  + ]:         40 :   return (is_rational(gel(x,1)) && is_rational(gel(x,2)));
     971                 :            : }
     972                 :            : static GEN
     973                 :         20 : c_zero_gcd(GEN c)
     974                 :            : {
     975                 :         20 :   GEN x = gel(c,1), y = gel(c,2);
     976                 :         20 :   long tx = typ(x), ty = typ(y);
     977 [ +  + ][ -  + ]:         20 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
     978 [ +  - ][ +  + ]:         15 :   if (tx == t_PADIC || tx == t_INTMOD
     979 [ +  - ][ -  + ]:         15 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
     980                 :         20 :   return gauss_gcd(c, gen_0);
     981                 :            : }
     982                 :            : 
     983                 :            : /* y == 0 */
     984                 :            : static GEN
     985                 :    5915834 : zero_gcd(GEN x, long tx)
     986                 :            : {
     987                 :            :   pari_sp av;
     988   [ +  +  +  +  :    5915834 :   switch(tx)
          +  +  +  +  +  
                      + ]
     989                 :            :   {
     990                 :      22959 :     case t_INT: return absi(x);
     991                 :       6465 :     case t_FRAC: return absfrac(x);
     992                 :         20 :     case t_COMPLEX: return c_zero_gcd(x);
     993                 :         30 :     case t_REAL: return gen_1;
     994                 :         25 :     case t_PADIC: return powis(gel(x,2), valp(x));
     995                 :        380 :     case t_SER: return monomial(gen_1, valp(x), varn(x));
     996                 :            :     case t_POLMOD: {
     997                 :       5675 :       GEN d = gel(x,2);
     998 [ +  + ][ +  - ]:       5675 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
     999         [ -  + ]:        200 :       return isinexact(d)? zero_gcd(d, typ(d)): gcopy(d);
    1000                 :            :     }
    1001                 :            :     case t_POL:
    1002         [ +  - ]:    5723694 :       if (!isinexact(x)) break;
    1003                 :          0 :       av = avma;
    1004                 :          0 :       return gerepileupto(av,
    1005                 :          0 :         monomialcopy(content(x), RgX_val(x), varn(x))
    1006                 :            :       );
    1007                 :            : 
    1008                 :            :     case t_RFRAC:
    1009         [ +  - ]:     156351 :       if (!isinexact(x)) break;
    1010                 :          0 :       av = avma;
    1011                 :          0 :       return gerepileupto(av,
    1012                 :          0 :         gdiv(zero_gcd(gel(x,1), typ(gel(x,1))), gel(x,2))
    1013                 :            :       );
    1014                 :            :   }
    1015                 :    5915834 :   return gcopy(x);
    1016                 :            : }
    1017                 :            : 
    1018                 :            : /* tx = t_RFRAC, y considered as constant */
    1019                 :            : static GEN
    1020                 :     627638 : cont_gcd_rfrac(GEN x, GEN y)
    1021                 :            : {
    1022                 :     627638 :   pari_sp av = avma;
    1023                 :     627638 :   GEN cx; x = primitive_part(x, &cx);
    1024         [ +  + ]:     627638 :   return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
    1025                 :            : }
    1026                 :            : /* tx = t_POL, y considered as constant */
    1027                 :            : static GEN
    1028                 :    1627257 : cont_gcd_pol(GEN x, GEN y)
    1029                 :            : {
    1030                 :    1627257 :   pari_sp av = avma;
    1031                 :    1627257 :   return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
    1032                 :            : }
    1033                 :            : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
    1034                 :            : static GEN
    1035                 :       5445 : cont_gcd_gen(GEN x, GEN y)
    1036                 :            : {
    1037                 :       5445 :   pari_sp av = avma;
    1038                 :       5445 :   return gerepileupto(av, ggcd(content(x),y));
    1039                 :            : }
    1040                 :            : /* !is_const(tx), y considered as constant */
    1041                 :            : static GEN
    1042                 :    2217502 : cont_gcd(GEN x, long tx, GEN y)
    1043                 :            : {
    1044      [ +  +  + ]:    2217502 :   switch(tx)
    1045                 :            :   {
    1046                 :     584805 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1047                 :    1627252 :     case t_POL: return cont_gcd_pol(x,y);
    1048                 :    2217502 :     default: return cont_gcd_gen(x,y);
    1049                 :            :   }
    1050                 :            : }
    1051                 :            : static GEN
    1052                 :     269833 : gcdiq(GEN x, GEN y)
    1053                 :            : {
    1054                 :            :   GEN z;
    1055         [ +  + ]:     269833 :   if (!signe(x)) return Q_abs(y);
    1056                 :      97946 :   z = cgetg(3,t_FRAC);
    1057                 :      97946 :   gel(z,1) = gcdii(x,gel(y,1));
    1058                 :      97946 :   gel(z,2) = icopy(gel(y,2));
    1059                 :     269833 :   return z;
    1060                 :            : }
    1061                 :            : static GEN
    1062                 :     582648 : gcdqq(GEN x, GEN y)
    1063                 :            : {
    1064                 :     582648 :   GEN z = cgetg(3,t_FRAC);
    1065                 :     582648 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1066                 :     582648 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1067                 :     582648 :   return z;
    1068                 :            : }
    1069                 :            : /* assume x,y t_INT or t_FRAC */
    1070                 :            : GEN
    1071                 :   40259878 : Q_gcd(GEN x, GEN y)
    1072                 :            : {
    1073                 :   40259878 :   long tx = typ(x), ty = typ(y);
    1074         [ +  + ]:   40259878 :   if (tx == t_INT)
    1075         [ +  + ]:   39578823 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1076                 :            :   else
    1077         [ +  + ]:   40259878 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1078                 :            : }
    1079                 :            : 
    1080                 :            : GEN
    1081                 :   15723912 : ggcd(GEN x, GEN y)
    1082                 :            : {
    1083                 :   15723912 :   long l, i, vx, vy, tx = typ(x), ty = typ(y);
    1084                 :            :   pari_sp av, tetpil;
    1085                 :            :   GEN p1,z;
    1086                 :            : 
    1087 [ +  - ][ -  + ]:   15723912 :   if (is_noncalc_t(tx) || is_noncalc_t(ty)) pari_err_TYPE2("gcd",x,y);
    1088         [ -  + ]:   15723912 :   if (is_matvec_t(ty))
    1089                 :            :   {
    1090                 :          0 :     z = cgetg_copy(y, &l);
    1091         [ #  # ]:          0 :     for (i=1; i<l; i++) gel(z,i) = ggcd(x,gel(y,i));
    1092                 :          0 :     return z;
    1093                 :            :   }
    1094         [ -  + ]:   15723912 :   if (is_matvec_t(tx))
    1095                 :            :   {
    1096                 :          0 :     z = cgetg_copy(x, &l);
    1097         [ #  # ]:          0 :     for (i=1; i<l; i++) gel(z,i) = ggcd(gel(x,i),y);
    1098                 :          0 :     return z;
    1099                 :            :   }
    1100         [ +  + ]:   15723912 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1101                 :            :   /* tx <= ty */
    1102         [ +  + ]:   15723912 :   if (isrationalzero(x)) return zero_gcd(y, ty);
    1103         [ +  + ]:   13512578 :   if (isrationalzero(y)) return zero_gcd(x, tx);
    1104         [ +  + ]:    9814558 :   if (is_const_t(tx))
    1105                 :            :   {
    1106 [ +  + ][ +  +  :    4782073 :     if (ty == tx) switch(tx)
          +  +  +  +  +  
                      - ]
    1107                 :            :     {
    1108                 :            :       case t_INT:
    1109                 :    2334173 :         return gcdii(x,y);
    1110                 :            : 
    1111                 :     131275 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1112         [ +  + ]:     131275 :         if (equalii(gel(x,1),gel(y,1)))
    1113                 :     131270 :           gel(z,1) = icopy(gel(x,1));
    1114                 :            :         else
    1115                 :          5 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1116         [ -  + ]:     131275 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1117                 :            :         else
    1118                 :            :         {
    1119                 :     131275 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1120         [ +  + ]:     131275 :           if (!is_pm1(p1))
    1121                 :            :           {
    1122                 :         85 :             p1 = gcdii(p1,gel(y,2));
    1123         [ +  + ]:         85 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1124                 :         80 :             else p1 = gerepileuptoint(av, p1);
    1125                 :            :           }
    1126                 :     131275 :           gel(z,2) = p1;
    1127                 :            :         }
    1128                 :     131275 :         return z;
    1129                 :            : 
    1130                 :            :       case t_FRAC:
    1131                 :      72075 :         return gcdqq(x,y);
    1132                 :            : 
    1133                 :            :       case t_FFELT:
    1134         [ -  + ]:      18180 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1135 [ -  + ][ #  # ]:      18180 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1136                 :            : 
    1137                 :            :       case t_COMPLEX:
    1138 [ +  + ][ +  - ]:         10 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1139                 :          5 :         return triv_cont_gcd(y,x);
    1140                 :            : 
    1141                 :            :       case t_PADIC:
    1142         [ +  + ]:         10 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1143                 :          5 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1144                 :            : 
    1145                 :            :       case t_QUAD:
    1146                 :        185 :         av=avma; p1=gdiv(x,y);
    1147         [ +  + ]:        185 :         if (gequal0(gel(p1,3)))
    1148                 :            :         {
    1149                 :         10 :           p1=gel(p1,2);
    1150         [ +  + ]:         10 :           if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
    1151                 :          5 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1152                 :            :         }
    1153 [ +  + ][ +  + ]:        175 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
    1154                 :        165 :         p1 = ginv(p1); avma=av;
    1155 [ +  + ][ +  - ]:        165 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1156                 :        155 :         return triv_cont_gcd(y,x);
    1157                 :            : 
    1158                 :          0 :       default: return gen_1; /* t_REAL */
    1159                 :            :     }
    1160 [ +  + ][ +  +  :    2226165 :     if (is_const_t(ty)) switch(tx)
          +  +  +  +  +  
                      - ]
    1161                 :            :     {
    1162                 :            :       case t_INT:
    1163   [ +  +  +  +  :      26813 :         switch(ty)
             +  +  +  - ]
    1164                 :            :         {
    1165                 :         70 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1166                 :         70 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1167                 :         70 :             p1 = gcdii(gel(y,1),gel(y,2));
    1168         [ +  + ]:         70 :             if (!is_pm1(p1)) {
    1169                 :         25 :               p1 = gcdii(x,p1);
    1170         [ +  + ]:         25 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1171                 :            :               else
    1172                 :         20 :                 p1 = gerepileuptoint(av, p1);
    1173                 :            :             }
    1174                 :         70 :             gel(z,2) = p1; return z;
    1175                 :            : 
    1176                 :         25 :           case t_REAL: return gen_1;
    1177                 :            : 
    1178                 :            :           case t_FRAC:
    1179                 :      24313 :             return gcdiq(x,y);
    1180                 :            : 
    1181                 :            :           case t_COMPLEX:
    1182         [ +  - ]:         15 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1183                 :          0 :             return triv_cont_gcd(y,x);
    1184                 :            : 
    1185                 :            :           case t_FFELT:
    1186         [ +  + ]:       2150 :             if (!FF_equal0(y)) return FF_1(y);
    1187         [ +  + ]:        295 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1188                 :            : 
    1189                 :            :           case t_PADIC:
    1190                 :         55 :             return padic_gcd(x,y);
    1191                 :            : 
    1192                 :            :           case t_QUAD:
    1193                 :        185 :             return triv_cont_gcd(y,x);
    1194                 :            :           default:
    1195                 :          0 :             pari_err_TYPE2("gcd",x,y);
    1196                 :            :         }
    1197                 :            : 
    1198                 :            :       case t_REAL:
    1199         [ +  + ]:         35 :         switch(ty)
    1200                 :            :         {
    1201                 :            :           case t_INTMOD:
    1202                 :            :           case t_FFELT:
    1203                 :         15 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1204                 :         20 :           default: return gen_1;
    1205                 :            :         }
    1206                 :            : 
    1207                 :            :       case t_INTMOD:
    1208   [ +  +  +  +  :         50 :         switch(ty)
                      - ]
    1209                 :            :         {
    1210                 :            :           case t_FRAC:
    1211                 :         10 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
    1212         [ +  + ]:         10 :             if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
    1213                 :          5 :             return ggcd(gel(y,1), x);
    1214                 :            : 
    1215                 :            :           case t_FFELT:
    1216                 :            :           {
    1217                 :         25 :             GEN p = gel(y,4);
    1218         [ +  + ]:         25 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1219         [ +  + ]:         20 :             if (!FF_equal0(y)) return FF_1(y);
    1220         [ -  + ]:          5 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1221                 :            :           }
    1222                 :            : 
    1223                 :            :           case t_COMPLEX: case t_QUAD:
    1224                 :         10 :             return triv_cont_gcd(y,x);
    1225                 :            : 
    1226                 :            :           case t_PADIC:
    1227                 :          5 :             return padic_gcd(x,y);
    1228                 :            : 
    1229                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1230                 :            :         }
    1231                 :            : 
    1232                 :            :       case t_FRAC:
    1233   [ +  +  +  +  :        170 :         switch(ty)
                      - ]
    1234                 :            :         {
    1235                 :            :           case t_COMPLEX:
    1236         [ +  + ]:         10 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1237                 :            :           case t_QUAD:
    1238                 :        115 :             return triv_cont_gcd(y,x);
    1239                 :            :           case t_FFELT:
    1240                 :            :           {
    1241                 :         45 :             GEN p = gel(y,4);
    1242         [ +  + ]:         45 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1243         [ +  + ]:         25 :             if (!FF_equal0(y)) return FF_1(y);
    1244         [ +  + ]:         10 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1245                 :            :           }
    1246                 :            : 
    1247                 :            :           case t_PADIC:
    1248                 :          5 :             return padic_gcd(x,y);
    1249                 :            : 
    1250                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1251                 :            :         }
    1252                 :            :       case t_FFELT:
    1253         [ +  + ]:         75 :         switch(ty)
    1254                 :            :         {
    1255                 :            :           case t_PADIC:
    1256                 :            :           {
    1257                 :         45 :             GEN p = gel(y,2);
    1258                 :         45 :             long v = valp(y);
    1259 [ +  + ][ +  + ]:         45 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1260 [ +  + ][ +  + ]:         20 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1261                 :            :           }
    1262                 :         30 :           default: pari_err_TYPE2("gcd",x,y);
    1263                 :            :         }
    1264                 :            : 
    1265                 :            :       case t_COMPLEX:
    1266         [ +  - ]:         10 :         switch(ty)
    1267                 :            :         {
    1268                 :            :           case t_PADIC:
    1269                 :         10 :           case t_QUAD: return triv_cont_gcd(x,y);
    1270                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1271                 :            :         }
    1272                 :            : 
    1273                 :            :       case t_PADIC:
    1274         [ +  - ]:          5 :         switch(ty)
    1275                 :            :         {
    1276                 :          5 :           case t_QUAD: return triv_cont_gcd(y,x);
    1277                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1278                 :            :         }
    1279                 :            : 
    1280                 :          0 :       default: return gen_1; /* tx = t_REAL */
    1281                 :            :     }
    1282                 :    2199007 :     return cont_gcd(y,ty, x);
    1283                 :            :   }
    1284                 :            : 
    1285         [ +  + ]:    5032485 :   if (tx == t_POLMOD)
    1286                 :            :   {
    1287         [ +  + ]:       5075 :     if (ty == t_POLMOD)
    1288                 :            :     {
    1289                 :       4865 :       GEN T = gel(x,1);
    1290                 :       4865 :       z = cgetg(3,t_POLMOD);
    1291         [ +  - ]:       4865 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1292                 :       4865 :       gel(z,1) = T;
    1293         [ -  + ]:       4865 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1294                 :            :       else
    1295                 :            :       {
    1296                 :            :         GEN X, Y, d;
    1297                 :       4865 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1298                 :       4865 :         d = ggcd(content(X), content(Y));
    1299         [ +  + ]:       4865 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1300                 :       4865 :         p1 = ggcd(T, X);
    1301                 :       4865 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1302                 :            :       }
    1303                 :       4865 :       return z;
    1304                 :            :     }
    1305                 :        210 :     vx = varn(gel(x,1));
    1306      [ +  +  - ]:        210 :     switch(ty)
    1307                 :            :     {
    1308                 :            :       case t_POL:
    1309                 :        190 :         vy = varn(y);
    1310         [ +  + ]:        190 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1311                 :        185 :         z = cgetg(3,t_POLMOD);
    1312                 :        185 :         gel(z,1) = RgX_copy(gel(x,1));
    1313                 :        185 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1314                 :        185 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1315                 :        185 :         return z;
    1316                 :            : 
    1317                 :            :       case t_RFRAC:
    1318                 :         20 :         vy = varn(gel(y,2));
    1319         [ -  + ]:         20 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1320                 :         20 :         av = avma;
    1321                 :         20 :         p1 = ggcd(gel(x,1),gel(y,2));
    1322         [ +  + ]:         20 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1323                 :         15 :         avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1324                 :            :     }
    1325                 :            :   }
    1326                 :            : 
    1327                 :    5027410 :   vx = gvar(x);
    1328                 :    5027410 :   vy = gvar(y);
    1329         [ +  + ]:    5027410 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1330         [ +  + ]:    5019085 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1331                 :            : 
    1332                 :            :   /* vx = vy: same main variable */
    1333   [ +  +  +  - ]:    5008915 :   switch(tx)
    1334                 :            :   {
    1335                 :            :     case t_POL:
    1336   [ +  +  +  - ]:    4885945 :       switch(ty)
    1337                 :            :       {
    1338                 :    4838962 :         case t_POL: return RgX_gcd(x,y);
    1339                 :            :         case t_SER:
    1340                 :       4150 :           z = ggcd(content(x), content(y));
    1341                 :       4150 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1342                 :      42833 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1343                 :            :       }
    1344                 :          0 :       break;
    1345                 :            : 
    1346                 :            :     case t_SER:
    1347                 :         10 :       z = ggcd(content(x), content(y));
    1348      [ +  +  - ]:         10 :       switch(ty)
    1349                 :            :       {
    1350                 :          5 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1351                 :          5 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1352                 :            :       }
    1353                 :          0 :       break;
    1354                 :            : 
    1355                 :            :     case t_RFRAC:
    1356                 :            :     {
    1357                 :     122960 :       GEN xd = gel(x,2), yd = gel(y,2);
    1358         [ -  + ]:     122960 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1359                 :     122960 :       z = cgetg(3,t_RFRAC); av = avma;
    1360                 :     122960 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1361                 :     122960 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1362                 :            :     }
    1363                 :            :   }
    1364                 :          0 :   pari_err_TYPE2("gcd",x,y);
    1365                 :   15723807 :   return NULL; /* not reached */
    1366                 :            : }
    1367                 :            : GEN
    1368         [ +  - ]:       1005 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
    1369                 :            : 
    1370                 :            : /* x a t_VEC,t_COL or t_MAT */
    1371                 :            : static GEN
    1372                 :          0 : vec_lcm(GEN x)
    1373                 :            : {
    1374         [ #  # ]:          0 :   if (typ(x) == t_MAT)
    1375                 :            :   {
    1376                 :          0 :     long i, l = lg(x);
    1377                 :          0 :     GEN z = cgetg(l, t_VEC);
    1378         [ #  # ]:          0 :     for (i = 1; i < l; i++) gel(z,i) = glcm0(gel(x,i), NULL);
    1379                 :          0 :     x = z;
    1380                 :            :   }
    1381                 :          0 :   return glcm0(x, NULL);
    1382                 :            : }
    1383                 :            : static GEN
    1384                 :       2040 : scal_lcm(GEN x, GEN y)
    1385                 :            : {
    1386                 :       2040 :   pari_sp av = avma;
    1387                 :       2040 :   long tx = typ(x), ty = typ(y);
    1388         [ -  + ]:       2040 :   if (is_matvec_t(tx)) x = vec_lcm(x);
    1389         [ -  + ]:       2040 :   if (is_matvec_t(ty)) y = vec_lcm(y);
    1390                 :       2040 :   return gerepileupto(av, glcm(x, y));
    1391                 :            : }
    1392                 :            : 
    1393                 :            : static GEN
    1394                 :         35 : fix_lcm(GEN x)
    1395                 :            : {
    1396                 :            :   GEN t;
    1397      [ -  +  - ]:         35 :   switch(typ(x))
    1398                 :            :   {
    1399         [ #  # ]:          0 :     case t_INT: if (signe(x)<0) x = negi(x);
    1400                 :          0 :       break;
    1401                 :            :     case t_POL:
    1402         [ -  + ]:         35 :       if (lg(x) <= 2) break;
    1403                 :         35 :       t = leading_term(x);
    1404 [ +  - ][ -  + ]:         35 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1405                 :            :   }
    1406                 :         35 :   return x;
    1407                 :            : }
    1408                 :            : 
    1409                 :            : GEN
    1410                 :       2035 : glcm0(GEN x, GEN y) {
    1411 [ +  + ][ -  + ]:       2035 :   if (!y && lg(x) == 2)
    1412                 :            :   {
    1413                 :          0 :     long tx = typ(x);
    1414         [ #  # ]:          0 :     if (is_vec_t(tx))
    1415                 :            :     {
    1416                 :          0 :       x = gel(x,1);
    1417                 :          0 :       tx = typ(x);
    1418         [ #  # ]:          0 :       return is_matvec_t(tx)? vec_lcm(x): fix_lcm(x);
    1419                 :            :     }
    1420                 :            :   }
    1421                 :       2035 :   return gassoc_proto(scal_lcm,x,y);
    1422                 :            : }
    1423                 :            : 
    1424                 :            : GEN
    1425                 :       2085 : glcm(GEN x, GEN y)
    1426                 :            : {
    1427                 :            :   long tx, ty, i, l;
    1428                 :            :   pari_sp av;
    1429                 :            :   GEN p1, z;
    1430                 :            : 
    1431                 :       2085 :   ty = typ(y);
    1432         [ -  + ]:       2085 :   if (is_matvec_t(ty))
    1433                 :            :   {
    1434                 :          0 :     z = cgetg_copy(y, &l);
    1435         [ #  # ]:          0 :     for (i=1; i<l; i++) gel(z,i) = glcm(x,gel(y,i));
    1436                 :          0 :     return z;
    1437                 :            :   }
    1438                 :       2085 :   tx = typ(x);
    1439         [ -  + ]:       2085 :   if (is_matvec_t(tx))
    1440                 :            :   {
    1441                 :          0 :     z = cgetg_copy(x, &l);
    1442         [ #  # ]:          0 :     for (i=1; i<l; i++) gel(z,i) = glcm(gel(x,i),y);
    1443                 :          0 :     return z;
    1444                 :            :   }
    1445 [ +  + ][ +  - ]:       2085 :   if (tx==t_INT && ty==t_INT) return lcmii(x,y);
    1446         [ -  + ]:         35 :   if (gequal0(x)) return gen_0;
    1447                 :            : 
    1448                 :         35 :   av = avma;
    1449         [ +  - ]:         35 :   p1 = ggcd(x,y); if (!gequal1(p1)) y = gdiv(y,p1);
    1450                 :       2085 :   return gerepileupto(av, fix_lcm(gmul(x,y)));
    1451                 :            : }
    1452                 :            : 
    1453                 :            : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
    1454                 :            : static int
    1455                 :      27655 : pol_approx0(GEN r, GEN x, int exact)
    1456                 :            : {
    1457                 :            :   long i, lx,lr;
    1458         [ +  - ]:      27655 :   if (exact) return gequal0(r);
    1459                 :          0 :   lx = lg(x);
    1460         [ #  # ]:          0 :   lr = lg(r); if (lr < lx) lx = lr;
    1461         [ #  # ]:          0 :   for (i=2; i<lx; i++)
    1462         [ #  # ]:          0 :     if (!approx_0(gel(r,i), gel(x,i))) return 0;
    1463                 :      27655 :   return 1;
    1464                 :            : }
    1465                 :            : 
    1466                 :            : GEN
    1467                 :       6200 : RgX_gcd_simple(GEN x, GEN y)
    1468                 :            : {
    1469                 :       6200 :   pari_sp av1, av = avma, lim = stack_lim(av, 1);
    1470                 :       6200 :   GEN r, yorig = y;
    1471 [ +  - ][ +  - ]:       6200 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1472                 :            : 
    1473                 :            :   for(;;)
    1474                 :            :   {
    1475                 :      27655 :     av1 = avma; r = RgX_rem(x,y);
    1476         [ +  + ]:      27655 :     if (pol_approx0(r, x, exact))
    1477                 :            :     {
    1478                 :       6200 :       avma = av1;
    1479         [ +  + ]:       6200 :       if (y == yorig) return RgX_copy(y);
    1480                 :       5885 :       y = normalizepol_approx(y, lg(y));
    1481         [ +  + ]:       5885 :       if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
    1482                 :         60 :       return gerepileupto(av,y);
    1483                 :            :     }
    1484                 :      21455 :     x = y; y = r;
    1485         [ -  + ]:      21455 :     if (low_stack(lim,stack_lim(av,1))) {
    1486         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1487                 :          0 :       gerepileall(av,2, &x,&y);
    1488                 :            :     }
    1489                 :      27655 :   }
    1490                 :            : }
    1491                 :            : GEN
    1492                 :          0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
    1493                 :            : {
    1494                 :          0 :   pari_sp av = avma;
    1495                 :            :   GEN q, r, d, d1, u, v, v1;
    1496 [ #  # ][ #  # ]:          0 :   int exact = !(isinexactreal(a) || isinexactreal(b));
    1497                 :            : 
    1498                 :          0 :   d = a; d1 = b; v = gen_0; v1 = gen_1;
    1499                 :            :   for(;;)
    1500                 :            :   {
    1501         [ #  # ]:          0 :     if (pol_approx0(d1, a, exact)) break;
    1502                 :          0 :     q = poldivrem(d,d1, &r);
    1503                 :          0 :     v = gsub(v, gmul(q,v1));
    1504                 :          0 :     u=v; v=v1; v1=u;
    1505                 :          0 :     u=r; d=d1; d1=u;
    1506                 :          0 :   }
    1507                 :          0 :   u = gsub(d, gmul(b,v));
    1508                 :          0 :   u = RgX_div(u,a);
    1509                 :            : 
    1510                 :          0 :   gerepileall(av, 3, &u,&v,&d);
    1511                 :          0 :   *pu = u;
    1512                 :          0 :   *pv = v; return d;
    1513                 :            : }
    1514                 :            : /*******************************************************************/
    1515                 :            : /*                                                                 */
    1516                 :            : /*                    CONTENT / PRIMITIVE PART                     */
    1517                 :            : /*                                                                 */
    1518                 :            : /*******************************************************************/
    1519                 :            : 
    1520                 :            : GEN
    1521                 :   53898874 : content(GEN x)
    1522                 :            : {
    1523                 :   53898874 :   long lx, i, t, tx = typ(x);
    1524                 :   53898874 :   pari_sp av = avma;
    1525                 :            :   GEN c;
    1526                 :            : 
    1527         [ +  + ]:   53898874 :   if (is_scalar_t(tx)) return zero_gcd(x, tx);
    1528   [ +  +  +  +  :   53892819 :   switch(tx)
                   +  - ]
    1529                 :            :   {
    1530                 :            :     case t_RFRAC:
    1531                 :            :     {
    1532                 :     627648 :       GEN n = gel(x,1), d = gel(x,2);
    1533                 :            :       /* -- varncmp(vn, vd) < 0 can't happen
    1534                 :            :        * -- if n is POLMOD, its main variable (in the sense of gvar2)
    1535                 :            :        *    has lower priority than denominator */
    1536 [ +  - ][ +  + ]:     627648 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1537         [ -  + ]:     513843 :         n = isinexact(n)? zero_gcd(n, typ(n)): gcopy(n);
    1538                 :            :       else
    1539                 :     113805 :         n = content(n);
    1540                 :     627648 :       return gerepileupto(av, gdiv(n, content(d)));
    1541                 :            :     }
    1542                 :            : 
    1543                 :            :     case t_VEC: case t_COL:
    1544         [ -  + ]:    5454124 :       lx = lg(x); if (lx==1) return gen_1;
    1545                 :    5454124 :       break;
    1546                 :            : 
    1547                 :            :     case t_MAT:
    1548                 :            :     {
    1549                 :            :       long hx, j;
    1550                 :        145 :       lx = lg(x);
    1551         [ -  + ]:        145 :       if (lx == 1) return gen_1;
    1552                 :        145 :       hx = lgcols(x);
    1553         [ -  + ]:        145 :       if (hx == 1) return gen_1;
    1554         [ -  + ]:        145 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1555         [ -  + ]:        145 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1556                 :        145 :       c = content(gel(x,1));
    1557         [ +  + ]:        290 :       for (j=2; j<lx; j++)
    1558         [ +  + ]:        435 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1559 [ +  - ][ -  + ]:        145 :       if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
    1560                 :        145 :       return gerepileupto(av,c);
    1561                 :            :     }
    1562                 :            : 
    1563                 :            :     case t_POL: case t_SER:
    1564         [ +  + ]:   47810892 :       lx = lg(x); if (lx == 2) return gen_0;
    1565                 :   47809022 :       break;
    1566                 :            :     case t_QFR: case t_QFI:
    1567                 :         10 :       lx = 4; break;
    1568                 :            : 
    1569                 :          0 :     default: pari_err_TYPE("content",x);
    1570                 :          0 :       return NULL; /* not reached */
    1571                 :            :   }
    1572         [ +  + ]:  151663711 :   for (i=lontyp[tx]; i<lx; i++)
    1573         [ +  + ]:  107765137 :     if (typ(gel(x,i)) != t_INT) break;
    1574                 :   53263156 :   lx--; c = gel(x,lx);
    1575         [ -  + ]:   53263156 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1576         [ +  + ]:   53263156 :   if (i > lx)
    1577                 :            :   { /* integer coeffs */
    1578         [ +  + ]:   44345156 :     while (lx-- > lontyp[tx])
    1579                 :            :     {
    1580                 :   43147791 :       c = gcdii(c, gel(x,lx));
    1581         [ +  + ]:   43147791 :       if (is_pm1(c)) { avma=av; return gen_1; }
    1582                 :            :     }
    1583                 :            :   }
    1584                 :            :   else
    1585                 :            :   {
    1586         [ +  + ]:    9364582 :     if (isinexact(c)) c = zero_gcd(c, typ(c));
    1587         [ +  + ]:   22681014 :     while (lx-- > lontyp[tx])
    1588                 :            :     {
    1589                 :   13316432 :       GEN d = gel(x,lx);
    1590         [ -  + ]:   13316432 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1591                 :   13316432 :       c = ggcd(c, d);
    1592                 :            :     }
    1593         [ -  + ]:    9364582 :     if (isinexact(c)) { avma=av; return gen_1; }
    1594                 :            :   }
    1595      [ +  -  + ]:   10561947 :   switch(typ(c))
    1596                 :            :   {
    1597                 :            :     case t_INT:
    1598         [ +  + ]:    1200530 :       if (signe(c) < 0) c = negi(c);
    1599                 :    1200530 :       break;
    1600                 :            :     case t_VEC: case t_COL: case t_MAT:
    1601                 :          0 :       pari_err_TYPE("content",x);
    1602                 :            :   }
    1603                 :            : 
    1604         [ +  + ]:   53898874 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1605                 :            : }
    1606                 :            : 
    1607                 :            : GEN
    1608                 :     973850 : primitive_part(GEN x, GEN *ptc)
    1609                 :            : {
    1610                 :     973850 :   pari_sp av = avma;
    1611                 :     973850 :   GEN c = content(x);
    1612         [ +  + ]:     973850 :   if (gequal1(c)) { avma = av; c = NULL; }
    1613         [ +  + ]:      78432 :   else if (!gequal0(c)) x = gdiv(x,c);
    1614         [ +  + ]:     973850 :   if (ptc) *ptc = c;
    1615                 :     973850 :   return x;
    1616                 :            : }
    1617                 :            : GEN
    1618                 :       1810 : primpart(GEN x) { return primitive_part(x, NULL); }
    1619                 :            : 
    1620                 :            : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1621                 :            :  * of Q(x2,...,xn)[x1] */
    1622                 :            : GEN
    1623                 :   57036423 : Q_content(GEN x)
    1624                 :            : {
    1625                 :            :   long i, l;
    1626                 :            :   GEN d;
    1627                 :            :   pari_sp av;
    1628                 :            : 
    1629   [ +  +  +  +  :   57036423 :   switch(typ(x))
                -  +  - ]
    1630                 :            :   {
    1631                 :   44812586 :     case t_INT:  return absi(x);
    1632                 :     627368 :     case t_FRAC: return absfrac(x);
    1633                 :            : 
    1634                 :            :     case t_VEC: case t_COL: case t_MAT:
    1635         [ +  + ]:    5684988 :       l = lg(x); if (l == 1) return gen_1;
    1636                 :    5684968 :       av = avma; d = Q_content(gel(x,1));
    1637         [ +  + ]:   26892461 :       for (i=2; i<l; i++)
    1638                 :            :       {
    1639                 :   21207493 :         d = Q_gcd(d, Q_content(gel(x,i)));
    1640         [ +  + ]:   21207493 :         if ((i & 255) == 0) d = gerepileupto(av, d);
    1641                 :            :       }
    1642                 :    5684968 :       return gerepileupto(av, d);
    1643                 :            : 
    1644                 :            :     case t_POL:
    1645         [ +  + ]:    5911466 :       l = lg(x); if (l == 2) return gen_0;
    1646                 :    5910441 :       av = avma; d = Q_content(gel(x,2));
    1647         [ +  + ]:   22563345 :       for (i=3; i<l; i++) d = Q_gcd(d, Q_content(gel(x,i)));
    1648                 :    5910441 :       return gerepileupto(av, d);
    1649                 :          0 :     case t_POLMOD: return Q_content(gel(x,2));
    1650                 :         15 :     case t_COMPLEX: return Q_gcd(Q_content(gel(x,1)), Q_content(gel(x,2)));
    1651                 :            :   }
    1652                 :          0 :   pari_err_TYPE("Q_content",x);
    1653                 :   57036423 :   return NULL; /* not reached */
    1654                 :            : }
    1655                 :            : 
    1656                 :            : GEN
    1657                 :       1095 : ZX_content(GEN x)
    1658                 :            : {
    1659                 :       1095 :   long i, l = lg(x);
    1660                 :            :   GEN d;
    1661                 :            :   pari_sp av;
    1662                 :            : 
    1663         [ -  + ]:       1095 :   if (l == 2) return gen_0;
    1664                 :       1095 :   d = gel(x,2);
    1665         [ -  + ]:       1095 :   if (l == 3) return absi(d);
    1666                 :       1095 :   av = avma;
    1667 [ +  - ][ +  + ]:      17515 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    1668         [ -  + ]:       1095 :   if (signe(d) < 0) d = absi(d);
    1669                 :       1095 :   return gerepileuptoint(av, d);
    1670                 :            : }
    1671                 :            : 
    1672                 :            : /* NOT MEMORY CLEAN (because of t_FRAC).
    1673                 :            :  * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
    1674                 :            :  * of Q(x2,...,xn)[x1] */
    1675                 :            : GEN
    1676                 :    7609614 : Q_denom(GEN x)
    1677                 :            : {
    1678                 :            :   long i, l;
    1679                 :            :   GEN d, D;
    1680                 :            :   pari_sp av;
    1681                 :            : 
    1682   [ +  +  +  +  :    7609614 :   switch(typ(x))
                      - ]
    1683                 :            :   {
    1684                 :    4322637 :     case t_INT: return gen_1;
    1685                 :    1968533 :     case t_FRAC: return gel(x,2);
    1686                 :            : 
    1687                 :            :     case t_VEC: case t_COL: case t_MAT:
    1688         [ +  + ]:    1140121 :       l = lg(x); if (l == 1) return gen_1;
    1689                 :    1139446 :       av = avma; d = Q_denom(gel(x,1));
    1690         [ +  + ]:    5186120 :       for (i=2; i<l; i++)
    1691                 :            :       {
    1692                 :    4046674 :         D = Q_denom(gel(x,i));
    1693         [ +  + ]:    4046674 :         if (D != gen_1) d = lcmii(d, D);
    1694         [ -  + ]:    4046674 :         if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1695                 :            :       }
    1696                 :    1139446 :       return gerepileuptoint(av, d);
    1697                 :            : 
    1698                 :            :     case t_POL:
    1699         [ +  + ]:     178323 :       l = lg(x); if (l == 2) return gen_1;
    1700                 :     178283 :       av = avma; d = Q_denom(gel(x,2));
    1701         [ +  + ]:     801544 :       for (i=3; i<l; i++)
    1702                 :            :       {
    1703                 :     623261 :         D = Q_denom(gel(x,i));
    1704         [ +  + ]:     623261 :         if (D != gen_1) d = lcmii(d, D);
    1705                 :            :       }
    1706                 :     178283 :       return gerepileuptoint(av, d);
    1707                 :            :   }
    1708                 :          0 :   pari_err_TYPE("Q_denom",x);
    1709                 :    7609614 :   return NULL; /* not reached */
    1710                 :            : }
    1711                 :            : 
    1712                 :            : GEN
    1713                 :     957623 : Q_remove_denom(GEN x, GEN *ptd)
    1714                 :            : {
    1715                 :     957623 :   GEN d = Q_denom(x);
    1716         [ +  + ]:     957623 :   if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
    1717         [ +  + ]:     957623 :   if (ptd) *ptd = d;
    1718                 :     957623 :   return x;
    1719                 :            : }
    1720                 :            : 
    1721                 :            : /* return y = x * d, assuming x rational, and d,y integral */
    1722                 :            : GEN
    1723                 :    6018966 : Q_muli_to_int(GEN x, GEN d)
    1724                 :            : {
    1725                 :            :   long i, l;
    1726                 :            :   GEN y, xn, xd;
    1727                 :            :   pari_sp av;
    1728                 :            : 
    1729         [ -  + ]:    6018966 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    1730   [ +  +  +  +  :    6018966 :   switch (typ(x))
                   -  - ]
    1731                 :            :   {
    1732                 :            :     case t_INT:
    1733                 :    2279945 :       return mulii(x,d);
    1734                 :            : 
    1735                 :            :     case t_FRAC:
    1736                 :    2689757 :       xn = gel(x,1);
    1737                 :    2689757 :       xd = gel(x,2); av = avma;
    1738                 :    2689757 :       y = mulii(xn, diviiexact(d, xd));
    1739                 :    2689757 :       return gerepileuptoint(av, y);
    1740                 :            : 
    1741                 :            :     case t_VEC: case t_COL: case t_MAT:
    1742                 :     695963 :       y = cgetg_copy(x, &l);
    1743         [ +  + ]:    4166770 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1744                 :     695963 :       return y;
    1745                 :            : 
    1746                 :            :     case t_POL:
    1747                 :     353301 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1748         [ +  + ]:    2272271 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1749                 :     353301 :       return y;
    1750                 :            : 
    1751                 :            :     case t_POLMOD:
    1752                 :          0 :       y = cgetg(3, t_POLMOD);
    1753                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1754                 :          0 :       gel(y,2) = Q_muli_to_int(gel(x,2), d);
    1755                 :          0 :       return y;
    1756                 :            :   }
    1757                 :          0 :   pari_err_TYPE("Q_muli_to_int",x);
    1758                 :    6018966 :   return NULL; /* not reached */
    1759                 :            : }
    1760                 :            : 
    1761                 :            : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    1762                 :            : static GEN
    1763                 :      17522 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    1764                 :            : {
    1765                 :            :   long i, l;
    1766                 :            :   GEN y, xn, xd;
    1767                 :            :   pari_sp av;
    1768                 :            : 
    1769   [ +  +  +  +  :      17522 :   switch(typ(x))
                   -  - ]
    1770                 :            :   {
    1771                 :            :     case t_INT:
    1772                 :       1820 :       av = avma; y = diviiexact(x,d);
    1773                 :       1820 :       return gerepileuptoint(av, mulii(y,n));
    1774                 :            : 
    1775                 :            :     case t_FRAC:
    1776                 :      11205 :       xn = gel(x,1);
    1777                 :      11205 :       xd = gel(x,2); av = avma;
    1778                 :      11205 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    1779                 :      11205 :       return gerepileuptoint(av, y);
    1780                 :            : 
    1781                 :            :     case t_VEC: case t_COL: case t_MAT:
    1782                 :       2462 :       y = cgetg_copy(x, &l);
    1783         [ +  + ]:       9293 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1784                 :       2462 :       return y;
    1785                 :            : 
    1786                 :            :     case t_POL:
    1787                 :       2035 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1788         [ +  + ]:       9420 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1789                 :       2035 :       return y;
    1790                 :            : 
    1791                 :            :     case t_POLMOD:
    1792                 :          0 :       y = cgetg(3, t_POLMOD);
    1793                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1794                 :          0 :       gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
    1795                 :          0 :       return y;
    1796                 :            :   }
    1797                 :          0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    1798                 :      17522 :   return NULL; /* not reached */
    1799                 :            : }
    1800                 :            : 
    1801                 :            : /* return x / d. x: rational; d,result: integral. */
    1802                 :            : static GEN
    1803                 :    8286713 : Q_divi_to_int(GEN x, GEN d)
    1804                 :            : {
    1805                 :            :   long i, l;
    1806                 :            :   GEN y;
    1807                 :            : 
    1808   [ +  +  +  -  :    8286713 :   switch(typ(x))
                      - ]
    1809                 :            :   {
    1810                 :            :     case t_INT:
    1811                 :    6824415 :       return diviiexact(x,d);
    1812                 :            : 
    1813                 :            :     case t_VEC: case t_COL: case t_MAT:
    1814                 :     644533 :       y = cgetg_copy(x, &l);
    1815         [ +  + ]:    4384651 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1816                 :     644533 :       return y;
    1817                 :            : 
    1818                 :            :     case t_POL:
    1819                 :     817765 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1820         [ +  + ]:    4270348 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1821                 :     817765 :       return y;
    1822                 :            : 
    1823                 :            :     case t_POLMOD:
    1824                 :          0 :       y = cgetg(3, t_POLMOD);
    1825                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1826                 :          0 :       gel(y,2) = Q_divi_to_int(gel(x,2), d);
    1827                 :          0 :       return y;
    1828                 :            :   }
    1829                 :          0 :   pari_err_TYPE("Q_divi_to_int",x);
    1830                 :    8286713 :   return NULL; /* not reached */
    1831                 :            : }
    1832                 :            : /* c t_FRAC */
    1833                 :            : static GEN
    1834                 :     116970 : Q_divq_to_int(GEN x, GEN c)
    1835                 :            : {
    1836                 :     116970 :   GEN n = gel(c,1), d = gel(c,2);
    1837         [ +  + ]:     116970 :   if (is_pm1(n)) {
    1838                 :     113664 :     GEN y = Q_muli_to_int(x,d);
    1839         [ -  + ]:     113664 :     if (signe(n) < 0) y = gneg(y);
    1840                 :     113664 :     return y;
    1841                 :            :   }
    1842                 :     116970 :   return Q_divmuli_to_int(x, n,d);
    1843                 :            : }
    1844                 :            : 
    1845                 :            : /* return y = x / c, assuming x,c rational, and y integral */
    1846                 :            : GEN
    1847                 :       1775 : Q_div_to_int(GEN x, GEN c)
    1848                 :            : {
    1849      [ +  +  - ]:       1775 :   switch(typ(c))
    1850                 :            :   {
    1851                 :       1715 :     case t_INT:  return Q_divi_to_int(x, c);
    1852                 :         60 :     case t_FRAC: return Q_divq_to_int(x, c);
    1853                 :            :   }
    1854                 :          0 :   pari_err_TYPE("Q_div_to_int",c);
    1855                 :       1775 :   return NULL; /* not reached */
    1856                 :            : }
    1857                 :            : /* return y = x * c, assuming x,c rational, and y integral */
    1858                 :            : GEN
    1859                 :          0 : Q_mul_to_int(GEN x, GEN c)
    1860                 :            : {
    1861                 :            :   GEN d, n;
    1862      [ #  #  # ]:          0 :   switch(typ(c))
    1863                 :            :   {
    1864                 :          0 :     case t_INT: return Q_muli_to_int(x, c);
    1865                 :            :     case t_FRAC:
    1866                 :          0 :       n = gel(c,1);
    1867                 :          0 :       d = gel(c,2);
    1868                 :          0 :       return Q_divmuli_to_int(x, d,n);
    1869                 :            :   }
    1870                 :          0 :   pari_err_TYPE("Q_mul_to_int",c);
    1871                 :          0 :   return NULL; /* not reached */
    1872                 :            : }
    1873                 :            : 
    1874                 :            : GEN
    1875                 :    7411855 : Q_primitive_part(GEN x, GEN *ptc)
    1876                 :            : {
    1877                 :    7411855 :   pari_sp av = avma;
    1878                 :    7411855 :   GEN c = Q_content(x);
    1879         [ +  + ]:    7411855 :   if (typ(c) == t_INT)
    1880                 :            :   {
    1881         [ +  + ]:    7294945 :     if (is_pm1(c)) { avma = av; c = NULL; }
    1882         [ +  + ]:    1092332 :     else if (signe(c)) x = Q_divi_to_int(x, c);
    1883                 :            :   }
    1884                 :     116910 :   else x = Q_divq_to_int(x, c);
    1885         [ +  + ]:    7411855 :   if (ptc) *ptc = c;
    1886                 :    7411855 :   return x;
    1887                 :            : }
    1888                 :            : GEN
    1889                 :     904095 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    1890                 :            : 
    1891                 :            : /*******************************************************************/
    1892                 :            : /*                                                                 */
    1893                 :            : /*                           SUBRESULTANT                          */
    1894                 :            : /*                                                                 */
    1895                 :            : /*******************************************************************/
    1896                 :            : /* for internal use */
    1897                 :            : GEN
    1898                 :    1285979 : gdivexact(GEN x, GEN y)
    1899                 :            : {
    1900                 :            :   long i,lx;
    1901                 :            :   GEN z;
    1902         [ +  + ]:    1285979 :   if (gequal1(y)) return x;
    1903   [ +  +  +  +  :    1272414 :   switch(typ(x))
                      + ]
    1904                 :            :   {
    1905                 :            :     case t_INT:
    1906         [ +  + ]:    1086021 :       if (typ(y)==t_INT) return diviiexact(x,y);
    1907         [ +  - ]:        240 :       if (!signe(x)) return gen_0;
    1908                 :          0 :       break;
    1909                 :            :     case t_INTMOD:
    1910                 :       2208 :     case t_POLMOD: return gmul(x,ginv(y));
    1911                 :            :     case t_POL:
    1912      [ -  +  + ]:     183800 :       switch(typ(y))
    1913                 :            :       {
    1914                 :            :         case t_INTMOD:
    1915                 :          0 :         case t_POLMOD: return gmul(x,ginv(y));
    1916                 :            :         case t_POL: { /* not stack-clean */
    1917                 :            :           long v;
    1918         [ +  + ]:      15675 :           if (varn(x)!=varn(y)) break;
    1919                 :      15190 :           v = RgX_valrem(y,&y);
    1920         [ +  + ]:      15190 :           if (v) x = RgX_shift_shallow(x,-v);
    1921         [ +  + ]:      15190 :           if (!degpol(y)) { y = gel(y,2); break; }
    1922                 :      13910 :           return RgX_div(x,y);
    1923                 :            :         }
    1924                 :            :       }
    1925                 :     169890 :       return RgX_Rg_divexact(x, y);
    1926                 :            :     case t_VEC: case t_COL: case t_MAT:
    1927                 :          5 :       lx = lg(x); z = new_chunk(lx);
    1928         [ +  + ]:         35 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1929                 :          5 :       z[0] = x[0]; return z;
    1930                 :            :   }
    1931         [ -  + ]:        380 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    1932                 :    1285979 :   return gdiv(x,y);
    1933                 :            : }
    1934                 :            : 
    1935                 :            : static GEN
    1936                 :      26937 : init_resultant(GEN x, GEN y)
    1937                 :            : {
    1938                 :      26937 :   long tx = typ(x), ty = typ(y), vx, vy;
    1939 [ +  - ][ -  + ]:      26937 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    1940                 :            :   {
    1941 [ #  # ][ #  # ]:          0 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    1942         [ #  # ]:          0 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    1943         [ #  # ]:          0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    1944                 :          0 :     return gen_1;
    1945                 :            :   }
    1946         [ -  + ]:      26937 :   if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
    1947         [ -  + ]:      26937 :   if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
    1948 [ +  - ][ +  + ]:      26937 :   if (!signe(x) || !signe(y)) return gmul(RgX_get_0(x),RgX_get_0(y)); /*type*/
    1949                 :      26887 :   vx = varn(x);
    1950         [ +  - ]:      26887 :   vy = varn(y); if (vx == vy) return NULL;
    1951         [ #  # ]:      26937 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    1952                 :            : }
    1953                 :            : 
    1954                 :            : static long
    1955                 :      55989 : RgX_simpletype(GEN x)
    1956                 :            : {
    1957                 :      55989 :   long T = t_INT, i, lx = lg(x);
    1958         [ +  + ]:     361746 :   for (i = 2; i < lx; i++)
    1959                 :            :   {
    1960                 :     305767 :     GEN c = gel(x,i);
    1961                 :     305767 :     long tc = typ(c);
    1962      [ +  +  + ]:     305767 :     switch(tc) {
    1963                 :            :       case t_INT:
    1964                 :     273319 :         break;
    1965                 :            :       case t_FRAC:
    1966         [ +  + ]:      15986 :         if (T == t_INT) T = t_FRAC;
    1967                 :      15986 :         break;
    1968                 :            :       default:
    1969         [ +  + ]:      16462 :         if (isinexact(c)) return t_REAL;
    1970                 :      16452 :         T = 0; break;
    1971                 :            :     }
    1972                 :            :   }
    1973                 :      55989 :   return T;
    1974                 :            : }
    1975                 :            : 
    1976                 :            : /* x an RgX, y a scalar */
    1977                 :            : static GEN
    1978                 :          0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    1979                 :            : {
    1980                 :          0 :   *V = gpowgs(y,degpol(x)-1);
    1981                 :          0 :   *U = gen_0; return gmul(y, *V);
    1982                 :            : }
    1983                 :            : 
    1984                 :            : /* return 0 if the subresultant chain can be interrupted.
    1985                 :            :  * Set u = NULL if the resultant is 0. */
    1986                 :            : static int
    1987                 :     127884 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    1988                 :            : {
    1989                 :     127884 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    1990                 :            :   long du, dv, dr, degq;
    1991                 :            : 
    1992         [ +  + ]:     127884 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    1993                 :     127754 :   du = degpol(*u);
    1994                 :     127754 :   dv = degpol(*v);
    1995                 :     127754 :   degq = du - dv;
    1996         [ +  + ]:     127754 :   if (*um1 == gen_1)
    1997                 :      79459 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    1998         [ +  + ]:      48295 :   else if (*um1 == gen_0)
    1999                 :      22438 :     u0 = gen_0;
    2000                 :            :   else /* except in those 2 cases, um1 is an RgX */
    2001                 :      25857 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2002                 :            : 
    2003         [ +  + ]:     127754 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2004                 :      79459 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2005                 :            :   else
    2006                 :      48295 :     u0 = gsub(u0, gmul(q,*uze));
    2007                 :            : 
    2008                 :     127754 :   *um1 = *uze;
    2009                 :     127754 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2010                 :            : 
    2011                 :     127754 :   *u = *v; c = *g; *g  = leading_term(*u);
    2012      [ +  +  + ]:     127754 :   switch(degq)
    2013                 :            :   {
    2014                 :         20 :     case 0: break;
    2015                 :            :     case 1:
    2016                 :     112476 :       c = gmul(*h,c); *h = *g; break;
    2017                 :            :     default:
    2018                 :      15258 :       c = gmul(gpowgs(*h,degq), c);
    2019                 :      15258 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2020                 :            :   }
    2021                 :     127754 :   *v  = RgX_Rg_divexact(r,c);
    2022                 :     127754 :   *uze= RgX_Rg_divexact(*uze,c);
    2023         [ +  + ]:     127754 :   if (both_odd(du, dv)) *signh = -*signh;
    2024                 :     127884 :   return (dr > 3);
    2025                 :            : }
    2026                 :            : 
    2027                 :            : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2028                 :            : static GEN
    2029                 :      79614 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2030                 :            : {
    2031                 :            :   pari_sp av, av2, lim;
    2032                 :      79614 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2033                 :            :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2034                 :            : 
    2035         [ -  + ]:      79614 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2036         [ -  + ]:      79614 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2037 [ +  + ][ -  + ]:      79614 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2038         [ -  + ]:      79609 :   if (tx != t_POL) {
    2039         [ #  # ]:          0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2040                 :          0 :     return scalar_res(y,x,V,U);
    2041                 :            :   }
    2042         [ -  + ]:      79609 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2043         [ -  + ]:      79609 :   if (varn(x) != varn(y))
    2044                 :          0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2045         [ #  # ]:          0 :                                         : scalar_res(y,x,V,U);
    2046                 :      79609 :   dx = degpol(x); dy = degpol(y); signh = 1;
    2047         [ +  + ]:      79609 :   if (dx < dy)
    2048                 :            :   {
    2049                 :      79239 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2050         [ +  + ]:      79239 :     if (both_odd(dx, dy)) signh = -signh;
    2051                 :            :   }
    2052         [ +  + ]:      79609 :   if (dy == 0)
    2053                 :            :   {
    2054                 :        215 :     *V = gpowgs(gel(y,2),dx-1);
    2055                 :        215 :     *U = gen_0; return gmul(*V,gel(y,2));
    2056                 :            :   }
    2057                 :      79394 :   av = avma;
    2058                 :      79394 :   u = x = primitive_part(x, &cu);
    2059                 :      79394 :   v = y = primitive_part(y, &cv);
    2060                 :      79394 :   g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
    2061                 :      79394 :   um1 = gen_1; uze = gen_0;
    2062                 :            :   for(;;)
    2063                 :            :   {
    2064         [ +  + ]:     127639 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2065         [ -  + ]:      48245 :     if (low_stack(lim,stack_lim(av2,1)))
    2066                 :            :     {
    2067         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2068                 :          0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2069                 :            :     }
    2070                 :      48245 :   }
    2071                 :            :   /* uze an RgX */
    2072         [ -  + ]:      79394 :   if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
    2073                 :      79394 :   z = gel(v,2); du = degpol(u);
    2074         [ +  + ]:      79394 :   if (du > 1)
    2075                 :            :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2076                 :       2499 :     p1 = gpowgs(gdiv(z,h),du-1);
    2077                 :       2499 :     z = gmul(z,p1);
    2078                 :       2499 :     uze = RgX_Rg_mul(uze, p1);
    2079                 :            :   }
    2080         [ +  + ]:      79394 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2081                 :            : 
    2082                 :      79394 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2083         [ -  + ]:      79394 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2084                 :            :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2085                 :      79394 :   p1 = gen_1;
    2086         [ +  + ]:      79394 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2087         [ +  + ]:      79394 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2088         [ +  + ]:      79394 :   cu = cu? gdiv(p1,cu): p1;
    2089         [ +  + ]:      79394 :   cv = cv? gdiv(p1,cv): p1;
    2090                 :      79394 :   z = gmul(z,p1);
    2091                 :      79394 :   *U = RgX_Rg_mul(uze,cu);
    2092                 :      79394 :   *V = RgX_Rg_mul(vze,cv);
    2093                 :      79614 :   return z;
    2094                 :            : }
    2095                 :            : GEN
    2096                 :          0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2097                 :            : {
    2098                 :          0 :   pari_sp av = avma;
    2099                 :          0 :   GEN z = subresext_i(x, y, U, V);
    2100                 :          0 :   gerepileall(av, 3, &z, U, V);
    2101                 :          0 :   return z;
    2102                 :            : }
    2103                 :            : 
    2104                 :            : static GEN
    2105                 :         20 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2106                 :            : {
    2107                 :         20 :   GEN x=content(y);
    2108                 :         20 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2109                 :            : }
    2110                 :            : 
    2111                 :            : static int
    2112                 :       8020 : must_negate(GEN x)
    2113                 :            : {
    2114                 :       8020 :   GEN t = leading_term(x);
    2115      [ +  -  + ]:       8020 :   switch(typ(t))
    2116                 :            :   {
    2117                 :            :     case t_INT: case t_REAL:
    2118                 :       6465 :       return (signe(t) < 0);
    2119                 :            :     case t_FRAC:
    2120                 :          0 :       return (signe(gel(t,1)) < 0);
    2121                 :            :   }
    2122                 :       8020 :   return 0;
    2123                 :            : }
    2124                 :            : 
    2125                 :            : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2126                 :            : GEN
    2127                 :        310 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2128                 :            : {
    2129                 :            :   pari_sp av, av2, tetpil, lim;
    2130                 :            :   long signh; /* junk */
    2131                 :        310 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2132                 :            :   GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2133                 :            : 
    2134         [ -  + ]:        310 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2135         [ -  + ]:        310 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2136         [ -  + ]:        310 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2137                 :        310 :   vx=varn(x);
    2138         [ +  + ]:        310 :   if (!signe(x))
    2139                 :            :   {
    2140         [ +  + ]:         10 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2141                 :          5 :     *U = pol_0(vx); *V = pol_0(vx);
    2142                 :          5 :     return pol_0(vx);
    2143                 :            :   }
    2144         [ +  + ]:        300 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2145                 :        285 :   dx = degpol(x); dy = degpol(y);
    2146         [ +  + ]:        285 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2147         [ +  + ]:        285 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2148                 :            : 
    2149                 :        155 :   av = avma;
    2150                 :        155 :   u = x = primitive_part(x, &cu);
    2151                 :        155 :   v = y = primitive_part(y, &cv);
    2152                 :        155 :   g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
    2153                 :        155 :   um1 = gen_1; uze = gen_0;
    2154                 :            :   for(;;)
    2155                 :            :   {
    2156         [ +  + ]:        165 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2157         [ -  + ]:         10 :     if (low_stack(lim,stack_lim(av2,1)))
    2158                 :            :     {
    2159         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2160                 :          0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2161                 :            :     }
    2162                 :         10 :   }
    2163         [ +  + ]:        155 :   if (uze != gen_0) {
    2164                 :            :     GEN r;
    2165                 :         25 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2166         [ -  + ]:         25 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2167         [ -  + ]:         25 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2168         [ -  + ]:         25 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2169                 :         25 :     p1 = ginv(content(v));
    2170                 :            :   }
    2171                 :            :   else /* y | x */
    2172                 :            :   {
    2173         [ +  + ]:        130 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2174                 :        130 :     uze = pol_0(vx);
    2175                 :        130 :     p1 = gen_1;
    2176                 :            :   }
    2177         [ +  + ]:        155 :   if (must_negate(v)) p1 = gneg(p1);
    2178                 :        155 :   tetpil = avma;
    2179                 :        155 :   z = RgX_Rg_mul(v,p1);
    2180                 :        155 :   *U = RgX_Rg_mul(uze,p1);
    2181                 :        155 :   *V = RgX_Rg_mul(vze,p1);
    2182                 :        155 :   gptr[0] = &z;
    2183                 :        155 :   gptr[1] = U;
    2184                 :        155 :   gptr[2] = V;
    2185                 :        310 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2186                 :            : }
    2187                 :            : 
    2188                 :            : int
    2189                 :         40 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2190                 :            : {
    2191                 :         40 :   pari_sp av = avma, av2, tetpil, lim;
    2192                 :            :   long signh; /* junk */
    2193                 :            :   long vx;
    2194                 :            :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2195                 :            : 
    2196         [ -  + ]:         40 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2197         [ -  + ]:         40 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2198         [ -  + ]:         40 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2199         [ -  + ]:         40 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2200         [ -  + ]:         40 :   if (!signe(T)) {
    2201         [ #  # ]:          0 :     if (degpol(x) <= amax) {
    2202                 :          0 :       *P = RgX_copy(x);
    2203                 :          0 :       *Q = pol_1(varn(x));
    2204                 :          0 :       return 1;
    2205                 :            :     }
    2206                 :          0 :     return 0;
    2207                 :            :   }
    2208         [ -  + ]:         40 :   if (amax+bmax >= degpol(T))
    2209                 :          0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2210                 :            :                     mkvec3(stoi(amax), stoi(bmax), T));
    2211                 :         40 :   vx = varn(T);
    2212                 :         40 :   u = x = primitive_part(x, &cu);
    2213                 :         40 :   v = T = primitive_part(T, &cv);
    2214                 :         40 :   g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
    2215                 :         40 :   um1 = gen_1; uze = gen_0;
    2216                 :            :   for(;;)
    2217                 :            :   {
    2218                 :         80 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2219 [ +  - ][ +  - ]:         80 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
                 [ -  + ]
    2220 [ +  - ][ +  + ]:         80 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2221         [ -  + ]:         40 :     if (low_stack(lim,stack_lim(av2,1)))
    2222                 :            :     {
    2223         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2224                 :          0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2225                 :            :     }
    2226                 :         40 :   }
    2227         [ -  + ]:         40 :   if (uze == gen_0)
    2228                 :            :   {
    2229                 :          0 :     avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    2230                 :          0 :     return 1;
    2231                 :            :   }
    2232         [ +  + ]:         40 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2233                 :         40 :   p1 = ginv(content(v));
    2234         [ -  + ]:         40 :   if (must_negate(v)) p1 = gneg(p1);
    2235                 :         40 :   tetpil = avma;
    2236                 :         40 :   *P = RgX_Rg_mul(v,p1);
    2237                 :         40 :   *Q = RgX_Rg_mul(uze,p1);
    2238                 :         40 :   gptr[0] = P;
    2239                 :         40 :   gptr[1] = Q;
    2240                 :         40 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2241                 :            : }
    2242                 :            : 
    2243                 :            : /*******************************************************************/
    2244                 :            : /*                                                                 */
    2245                 :            : /*                RESULTANT USING DUCOS VARIANT                    */
    2246                 :            : /*                                                                 */
    2247                 :            : /*******************************************************************/
    2248                 :            : /* x^n / y^(n-1), assume n > 0 */
    2249                 :            : static GEN
    2250                 :      11442 : Lazard(GEN x, GEN y, long n)
    2251                 :            : {
    2252                 :            :   long a;
    2253                 :            :   GEN c;
    2254                 :            : 
    2255         [ +  + ]:      11442 :   if (n == 1) return x;
    2256                 :        558 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2257                 :        558 :   c=x; n-=a;
    2258         [ +  + ]:       1226 :   while (a>1)
    2259                 :            :   {
    2260                 :        668 :     a>>=1; c=gdivexact(gsqr(c),y);
    2261         [ +  + ]:        668 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2262                 :            :   }
    2263                 :      11442 :   return c;
    2264                 :            : }
    2265                 :            : 
    2266                 :            : /* F (x/y)^(n-1), assume n >= 1 */
    2267                 :            : static GEN
    2268                 :      14142 : Lazard2(GEN F, GEN x, GEN y, long n)
    2269                 :            : {
    2270         [ +  + ]:      14142 :   if (n == 1) return F;
    2271                 :      14142 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2272                 :            : }
    2273                 :            : 
    2274                 :            : static GEN
    2275                 :      14142 : RgX_neg_i(GEN x, long lx)
    2276                 :            : {
    2277                 :            :   long i;
    2278                 :      14142 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2279         [ +  + ]:      40419 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2280                 :      14142 :   return y;
    2281                 :            : }
    2282                 :            : static GEN
    2283                 :      41211 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2284                 :            : {
    2285                 :            :   long i;
    2286                 :            :   GEN z;
    2287         [ +  + ]:      41211 :   if (isrationalzero(x)) return pol_0(varn(y));
    2288                 :      41206 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2289         [ +  + ]:     115827 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2290                 :      41211 :   return z;
    2291                 :            : }
    2292                 :            : static long
    2293                 :      40591 : reductum_lg(GEN x, long lx)
    2294                 :            : {
    2295                 :      40591 :   long i = lx-2;
    2296 [ +  + ][ +  + ]:      44566 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2297                 :      40591 :   return i+1;
    2298                 :            : }
    2299                 :            : 
    2300                 :            : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2301                 :            :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2302                 :            : static GEN
    2303                 :      14142 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2304                 :            : {
    2305                 :      14142 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_term(Z);
    2306                 :            :   long p, q, j, lP, lQ;
    2307                 :            :   pari_sp av, lim;
    2308                 :            : 
    2309                 :      14142 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2310                 :      14142 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2311                 :            :   /* p > q. Very often p - 1 = q */
    2312                 :      14142 :   av = avma; lim = stack_lim(av,1);
    2313                 :            :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2314                 :      14142 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2315                 :            : 
    2316         [ +  + ]:      14142 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2317         [ +  + ]:      15157 :   for (j = q+1; j < p; j++)
    2318                 :            :   {
    2319         [ +  + ]:       1015 :     if (degpol(H) == q-1)
    2320                 :            :     { /* h0 = coeff of degree q-1 = leading coeff */
    2321                 :        715 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2322                 :        715 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2323                 :            :     }
    2324                 :            :     else
    2325                 :        300 :       H = RgX_shift_shallow(H, 1);
    2326         [ +  + ]:       1015 :     if (j+2 < lP)
    2327                 :            :     {
    2328                 :        630 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2329         [ +  - ]:        630 :       A = A? RgX_add(A, TMP): TMP;
    2330                 :            :     }
    2331         [ +  + ]:       1015 :     if (low_stack(lim,stack_lim(av,1)))
    2332                 :            :     {
    2333         [ -  + ]:        101 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2334         [ +  - ]:        101 :       gerepileall(av,A?2:1,&H,&A);
    2335                 :            :     }
    2336                 :            :   }
    2337         [ +  + ]:      14142 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2338                 :      14142 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2339         [ +  + ]:      14142 :   A = A? RgX_add(A, TMP): TMP;
    2340                 :      14142 :   A = RgX_Rg_divexact(A, p0);
    2341         [ +  + ]:      14142 :   if (degpol(H) == q-1)
    2342                 :            :   {
    2343                 :      14047 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2344                 :      14047 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2345                 :            :   }
    2346                 :            :   else
    2347                 :         95 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2348                 :      14142 :   return RgX_Rg_divexact(A, s);
    2349                 :            : }
    2350                 :            : 
    2351                 :            : /* Ducos's subresultant */
    2352                 :            : GEN
    2353                 :      11467 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2354                 :            : {
    2355                 :            :   pari_sp av, av2, lim;
    2356                 :      11467 :   long dP, dQ, delta, sig = 1;
    2357                 :            :   GEN cP, cQ, Z, s;
    2358                 :            : 
    2359                 :      11467 :   dP = degpol(P);
    2360                 :      11467 :   dQ = degpol(Q); delta = dP - dQ;
    2361         [ +  + ]:      11467 :   if (delta < 0)
    2362                 :            :   {
    2363         [ +  + ]:       1005 :     if (both_odd(dP, dQ)) sig = -1;
    2364                 :       1005 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2365                 :            :   }
    2366         [ +  + ]:      11467 :   if (sol) *sol = gen_0;
    2367                 :      11467 :   av = avma;
    2368         [ +  + ]:      11467 :   if (dQ <= 0)
    2369                 :            :   {
    2370         [ -  + ]:        205 :     if (dQ < 0) return RgX_get_0(P);
    2371                 :        205 :     s = gpowgs(gel(Q,2), dP);
    2372         [ -  + ]:        205 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2373                 :        205 :     return s;
    2374                 :            :   }
    2375                 :      11262 :   P = primitive_part(P, &cP);
    2376                 :      11262 :   Q = primitive_part(Q, &cQ);
    2377                 :      11262 :   av2 = avma; lim = stack_lim(av2,1);
    2378                 :      11262 :   s = gpowgs(leading_term(Q),delta);
    2379         [ +  + ]:      11262 :   if (both_odd(dP, dQ)) sig = -sig;
    2380                 :      11262 :   Z = Q;
    2381                 :      11262 :   Q = RgX_pseudorem(P, Q);
    2382                 :      11262 :   P = Z;
    2383         [ +  + ]:      25404 :   while(degpol(Q) > 0)
    2384                 :            :   {
    2385                 :      14142 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2386                 :      14142 :     Z = Lazard2(Q, leading_term(Q), s, delta);
    2387         [ +  + ]:      14142 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2388                 :      14142 :     Q = nextSousResultant(P, Q, Z, s);
    2389                 :      14142 :     P = Z;
    2390         [ +  + ]:      14142 :     if (low_stack(lim,stack_lim(av,1)))
    2391                 :            :     {
    2392         [ -  + ]:          9 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2393                 :          9 :       gerepileall(av2,2,&P,&Q);
    2394                 :            :     }
    2395                 :      14142 :     s = leading_term(P);
    2396                 :            :   }
    2397         [ -  + ]:      11262 :   if (!signe(Q)) { avma = av; return RgX_get_0(Q); }
    2398                 :      11262 :   av2 = avma;
    2399                 :      11262 :   s = Lazard(leading_term(Q), s, degpol(P));
    2400         [ +  + ]:      11262 :   if (sig == -1) s = gneg(s);
    2401         [ -  + ]:      11262 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2402         [ +  + ]:      11262 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2403         [ +  + ]:      11262 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2404         [ +  + ]:      11467 :   return (avma == av2)? gerepilecopy(av, s): gerepileupto(av, s);
    2405                 :            : }
    2406                 :            : /* Return resultant(P,Q). If sol != NULL: set *sol to the last non-constant
    2407                 :            :  * polynomial in the prs IF the sequence was computed, and gen_0 otherwise.
    2408                 :            :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    2409                 :            :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    2410                 :            :  * in the "generic" case. */
    2411                 :            : GEN
    2412                 :      26917 : resultant_all(GEN P, GEN Q, GEN *sol)
    2413                 :            : {
    2414                 :            :   long TP, TQ;
    2415                 :            :   GEN s;
    2416                 :            : 
    2417         [ +  + ]:      26917 :   if (sol) *sol = gen_0;
    2418         [ +  + ]:      26917 :   if ((s = init_resultant(P,Q))) return s;
    2419 [ +  - ][ -  + ]:      26867 :   if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    2420                 :          0 :     return resultant2(P,Q); /* inexact */
    2421 [ +  + ][ +  + ]:      26867 :   if (TP && TQ) /* rational */
    2422                 :            :   {
    2423 [ +  - ][ +  + ]:      15475 :     if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    2424                 :       6832 :     return QX_resultant(P,Q);
    2425                 :            :   }
    2426                 :      26917 :   return RgX_resultant_all(P, Q, sol);
    2427                 :            : }
    2428                 :            : 
    2429                 :            : /*******************************************************************/
    2430                 :            : /*                                                                 */
    2431                 :            : /*               RESULTANT USING SYLVESTER MATRIX                  */
    2432                 :            : /*                                                                 */
    2433                 :            : /*******************************************************************/
    2434                 :            : static GEN
    2435                 :          0 : _pol_0(void)
    2436                 :            : {
    2437                 :          0 :   GEN x = cgetg(3,t_POL);
    2438                 :          0 :   x[1] = 0;
    2439                 :          0 :   gel(x,2) = gen_0; return x;
    2440                 :            : }
    2441                 :            : 
    2442                 :            : static GEN
    2443                 :        140 : sylvester_col(GEN x, long j, long d, long D)
    2444                 :            : {
    2445                 :        140 :   GEN c = cgetg(d+1,t_COL);
    2446                 :            :   long i;
    2447         [ +  + ]:        250 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    2448         [ +  + ]:        620 :   for (   ; i<=D; i++) gel(c,i) = gel(x,D-i+2);
    2449         [ +  + ]:        250 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    2450                 :        140 :   return c;
    2451                 :            : }
    2452                 :            : 
    2453                 :            : GEN
    2454                 :         30 : sylvestermatrix_i(GEN x, GEN y)
    2455                 :            : {
    2456                 :            :   long j,d,dx,dy;
    2457                 :            :   GEN M;
    2458                 :            : 
    2459         [ -  + ]:         30 :   dx = degpol(x); if (dx < 0) { dx = 0; x = _pol_0(); }
    2460         [ -  + ]:         30 :   dy = degpol(y); if (dy < 0) { dy = 0; y = _pol_0(); }
    2461                 :         30 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    2462         [ +  + ]:         90 :   for (j=1; j<=dy; j++) gel(M,j)    = sylvester_col(x,j,d,j+dx);
    2463         [ +  + ]:        110 :   for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
    2464                 :         30 :   return M;
    2465                 :            : }
    2466                 :            : 
    2467                 :            : GEN
    2468                 :         10 : sylvestermatrix(GEN x, GEN y)
    2469                 :            : {
    2470                 :            :   long i,j,lx;
    2471         [ -  + ]:         10 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    2472         [ -  + ]:         10 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    2473         [ -  + ]:         10 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    2474                 :         10 :   x = sylvestermatrix_i(x,y); lx = lg(x);
    2475         [ +  + ]:         40 :   for (i=1; i<lx; i++)
    2476         [ +  + ]:        120 :     for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
    2477                 :         10 :   return x;
    2478                 :            : }
    2479                 :            : 
    2480                 :            : GEN
    2481                 :         20 : resultant2(GEN x, GEN y)
    2482                 :            : {
    2483                 :            :   pari_sp av;
    2484                 :            :   GEN r;
    2485         [ -  + ]:         20 :   if ((r = init_resultant(x,y))) return r;
    2486                 :         20 :   av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
    2487                 :            : }
    2488                 :            : 
    2489                 :            : /* Let vx = main variable of x. Return a polynomial in variable 0:
    2490                 :            :  * if vx is 0 and v != 0, set *mx = 1 and replace vx by pol_x(MAXVARN)
    2491                 :            :  * if vx = v, copy x, set its main variable to 0 and return
    2492                 :            :  * if vx < v, return subst(x, v, pol_x(0))
    2493                 :            :  * if vx > v, return scalarpol(x, 0) */
    2494                 :            : static GEN
    2495                 :         50 : fix_pol(GEN x, long v, long *mx)
    2496                 :            : {
    2497                 :            :   long vx;
    2498         [ -  + ]:         50 :   if (typ(x) != t_POL) return x;
    2499                 :         50 :   vx = varn(x);
    2500         [ +  + ]:         50 :   if (v == vx)
    2501                 :            :   {
    2502         [ +  + ]:         20 :     if (v) { x = leafcopy(x); setvarn(x, 0); }
    2503                 :         20 :     return x;
    2504                 :            :   }
    2505         [ +  + ]:         30 :   if (!vx)
    2506                 :            :   {
    2507                 :         20 :     *mx = 1;
    2508                 :         20 :     x = poleval(x, pol_x(MAXVARN));
    2509                 :         20 :     vx = varn(x);
    2510         [ +  + ]:         20 :     if (v == vx) { setvarn(x, 0); return x; }
    2511                 :            :   }
    2512         [ +  + ]:         20 :   if (varncmp(v, vx) > 0)
    2513                 :            :   {
    2514                 :         10 :     x = gsubst(x,v,pol_x(0));
    2515 [ +  - ][ -  + ]:         10 :     if (typ(x) == t_POL && varn(x) == 0) return x;
    2516                 :            :   }
    2517                 :         50 :   return scalarpol_shallow(x, 0);
    2518                 :            : }
    2519                 :            : 
    2520                 :            : /* resultant of x and y with respect to variable v, or with respect to their
    2521                 :            :  * main variable if v < 0. */
    2522                 :            : GEN
    2523                 :        165 : polresultant0(GEN x, GEN y, long v, long flag)
    2524                 :            : {
    2525                 :        165 :   long m = 0;
    2526                 :        165 :   pari_sp av = avma;
    2527                 :            : 
    2528         [ +  + ]:        165 :   if (v >= 0)
    2529                 :            :   {
    2530                 :         15 :     x = fix_pol(x,v, &m);
    2531                 :         15 :     y = fix_pol(y,v, &m);
    2532                 :            :   }
    2533      [ +  +  - ]:        165 :   switch(flag)
    2534                 :            :   {
    2535                 :            :     case 2:
    2536                 :        160 :     case 0: x=resultant_all(x,y,NULL); break;
    2537                 :          5 :     case 1: x=resultant2(x,y); break;
    2538                 :          0 :     default: pari_err_FLAG("polresultant");
    2539                 :            :   }
    2540         [ +  + ]:        165 :   if (m) x = gsubst(x,MAXVARN,pol_x(0));
    2541                 :        165 :   return gerepileupto(av,x);
    2542                 :            : }
    2543                 :            : GEN
    2544                 :        360 : polresultantext0(GEN x, GEN y, long v)
    2545                 :            : {
    2546                 :            :   GEN R, U, V;
    2547                 :        360 :   long m = 0;
    2548                 :        360 :   pari_sp av = avma;
    2549                 :            : 
    2550         [ +  + ]:        360 :   if (v >= 0)
    2551                 :            :   {
    2552                 :         10 :     x = fix_pol(x,v, &m);
    2553                 :         10 :     y = fix_pol(y,v, &m);
    2554                 :            :   }
    2555                 :        360 :   R = subresext_i(x,y, &U,&V);
    2556         [ +  + ]:        360 :   if (m)
    2557                 :            :   {
    2558                 :          5 :     U = gsubst(gsubst(U, 0, pol_x(v)), MAXVARN, pol_x(0));
    2559                 :          5 :     V = gsubst(gsubst(V, 0, pol_x(v)), MAXVARN, pol_x(0));
    2560                 :          5 :     R = gsubst(R,MAXVARN,pol_x(0));
    2561                 :            :   }
    2562         [ +  + ]:        355 :   else if (v >= 0)
    2563                 :            :   {
    2564 [ +  - ][ +  - ]:          5 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    2565 [ +  - ][ +  - ]:          5 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    2566                 :            :   }
    2567                 :        360 :   return gerepilecopy(av, mkvec3(U,V,R));
    2568                 :            : }
    2569                 :            : GEN
    2570                 :        340 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    2571                 :            : 
    2572                 :            : /*******************************************************************/
    2573                 :            : /*                                                                 */
    2574                 :            : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    2575                 :            : /*                                                                 */
    2576                 :            : /*******************************************************************/
    2577                 :            : 
    2578                 :            : /* (v - x)^d */
    2579                 :            : static GEN
    2580                 :         60 : caract_const(pari_sp av, GEN x, long v, long d)
    2581                 :         60 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    2582                 :            : 
    2583                 :            : /* return caract(Mod(x,T)) in variable v */
    2584                 :            : GEN
    2585                 :       8902 : RgXQ_charpoly(GEN x, GEN T, long v)
    2586                 :            : {
    2587                 :       8902 :   pari_sp av = avma;
    2588                 :       8902 :   long d = degpol(T), dx, vx, vp;
    2589                 :            :   GEN ch, L;
    2590                 :            : 
    2591         [ -  + ]:       8902 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    2592                 :       8902 :   vx = varn(x);
    2593                 :       8902 :   vp = varn(T);
    2594         [ -  + ]:       8902 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    2595         [ -  + ]:       8902 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    2596                 :       8902 :   dx = degpol(x);
    2597         [ +  + ]:       8902 :   if (dx <= 0)
    2598         [ -  + ]:         10 :     return dx? monomial(gen_1, d, v): caract_const(av, gel(x,2), v, d);
    2599                 :            : 
    2600                 :       8892 :   x = RgX_neg(x);
    2601         [ -  + ]:       8892 :   if (varn(x) == MAXVARN) { setvarn(x, 0); T = leafcopy(T); setvarn(T, 0); }
    2602                 :       8892 :   gel(x,2) = gadd(gel(x,2), pol_x(MAXVARN));
    2603                 :       8892 :   ch = resultant_all(T, x, NULL);
    2604         [ +  - ]:       8892 :   if (v != MAXVARN)
    2605                 :            :   {
    2606 [ +  + ][ +  + ]:       8892 :     if (typ(ch) == t_POL && varn(ch) == MAXVARN)
    2607                 :       8827 :       setvarn(ch, v);
    2608                 :            :     else
    2609                 :         65 :       ch = gsubst(ch, MAXVARN, pol_x(v));
    2610                 :            :   }
    2611                 :            :   /* test for silly input: x mod (deg 0 polynomial) */
    2612         [ -  + ]:       8892 :   if (typ(ch) != t_POL) { avma = av; return pol_1(v); }
    2613                 :            : 
    2614                 :       8892 :   L = leading_term(ch);
    2615         [ -  + ]:       8892 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    2616                 :       8902 :   return gerepileupto(av, ch);
    2617                 :            : }
    2618                 :            : 
    2619                 :            : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    2620                 :            :  * algebra nf[t]/(Q(t)) */
    2621                 :            : GEN
    2622                 :        160 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    2623                 :            : {
    2624                 :        160 :   const char *f = "rnfcharpoly";
    2625                 :        160 :   long dQ = degpol(Q);
    2626                 :        160 :   pari_sp av = avma;
    2627                 :            :   GEN T;
    2628                 :            : 
    2629         [ +  - ]:        160 :   if (v < 0) v = 0;
    2630                 :        160 :   nf = checknf(nf); T = nf_get_pol(nf);
    2631                 :        160 :   Q = RgX_nffix(f, T,Q,0);
    2632   [ +  +  +  + ]:        160 :   switch(typ(x))
    2633                 :            :   {
    2634                 :            :     case t_INT:
    2635                 :         20 :     case t_FRAC: return caract_const(av, x, v, dQ);
    2636                 :            :     case t_POLMOD:
    2637                 :         65 :       x = polmod_nffix2(f,T,Q, x,0);
    2638                 :         35 :       break;
    2639                 :            :     case t_POL:
    2640         [ +  + ]:         40 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    2641                 :         30 :       break;
    2642                 :         35 :     default: pari_err_TYPE(f,x);
    2643                 :            :   }
    2644         [ +  + ]:         65 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    2645                 :            :   /* x a t_POL in variable vQ */
    2646         [ +  + ]:         35 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    2647         [ -  + ]:         35 :   if (dQ <= 1) return caract_const(av, constant_term(x), v, 1);
    2648                 :         85 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    2649                 :            : }
    2650                 :            : 
    2651                 :            : /*******************************************************************/
    2652                 :            : /*                                                                 */
    2653                 :            : /*                  GCD USING SUBRESULTANT                         */
    2654                 :            : /*                                                                 */
    2655                 :            : /*******************************************************************/
    2656                 :            : static int inexact(GEN x, int *simple, int *rational);
    2657                 :            : static int
    2658                 :    4835597 : isinexactall(GEN x, int *simple, int *rational)
    2659                 :            : {
    2660                 :    4835597 :   long i, lx = lg(x);
    2661         [ +  + ]:   22770508 :   for (i=2; i<lx; i++)
    2662         [ -  + ]:   17934911 :     if (inexact(gel(x,i), simple, rational)) return 1;
    2663                 :    4835597 :   return 0;
    2664                 :            : }
    2665                 :            : /* return 1 if coeff explosion is not possible */
    2666                 :            : static int
    2667                 :   17934951 : inexact(GEN x, int *simple, int *rational)
    2668                 :            : {
    2669                 :   17934951 :   int junk = 0;
    2670   [ +  -  +  -  :   17934951 :   switch(typ(x))
             -  +  +  +  
                      - ]
    2671                 :            :   {
    2672                 :   17874706 :     case t_INT: case t_FRAC: return 0;
    2673                 :            : 
    2674                 :          0 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    2675                 :            : 
    2676                 :            :     case t_INTMOD:
    2677                 :            :     case t_FFELT:
    2678                 :      40760 :       *rational = 0;
    2679         [ +  + ]:      40760 :       if (!*simple) *simple = 1;
    2680                 :      40760 :       return 0;
    2681                 :            : 
    2682                 :            :     case t_COMPLEX:
    2683                 :          0 :       *rational = 0;
    2684                 :          0 :       return inexact(gel(x,1), simple, rational)
    2685 [ #  # ][ #  # ]:          0 :           || inexact(gel(x,2), simple, rational);
    2686                 :            :     case t_QUAD:
    2687                 :          0 :       *rational = *simple = 0;
    2688                 :          0 :       return inexact(gel(x,2), &junk, rational)
    2689 [ #  # ][ #  # ]:          0 :           || inexact(gel(x,3), &junk, rational);
    2690                 :            : 
    2691                 :            :     case t_POLMOD:
    2692                 :       2795 :       *rational = 0;
    2693                 :       2795 :       return isinexactall(gel(x,1), simple, rational);
    2694                 :            :     case t_POL:
    2695                 :      16670 :       *rational = 0;
    2696                 :      16670 :       *simple = -1;
    2697                 :      16670 :       return isinexactall(x, &junk, rational);
    2698                 :            :     case t_RFRAC:
    2699                 :         20 :       *rational = 0;
    2700                 :         20 :       *simple = -1;
    2701                 :         40 :       return inexact(gel(x,1), &junk, rational)
    2702 [ +  - ][ -  + ]:         20 :           || inexact(gel(x,2), &junk, rational);
    2703                 :            :   }
    2704                 :          0 :   *rational = 0;
    2705                 :   17934951 :   *simple = -1; return 0;
    2706                 :            : }
    2707                 :            : 
    2708                 :            : /* x monomial, y t_POL in the same variable */
    2709                 :            : static GEN
    2710                 :    5296244 : gcdmonome(GEN x, GEN y)
    2711                 :            : {
    2712                 :    5296244 :   pari_sp av = avma;
    2713                 :    5296244 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    2714                 :    5296244 :   long i, l = lg(y);
    2715                 :    5296244 :   GEN t, v = cgetg(l, t_VEC);
    2716                 :    5296244 :   gel(v,1) = gel(x,dx+2);
    2717         [ +  + ]:   11757730 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    2718                 :    5296244 :   t = content(v); /* gcd(lc(x), cont(y)) */
    2719                 :    5296244 :   t = simplify_shallow(t);
    2720         [ +  + ]:    5296244 :   if (dx < e) e = dx;
    2721                 :    5296244 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    2722                 :            : }
    2723                 :            : 
    2724                 :            : /* x, y are t_POL in the same variable */
    2725                 :            : GEN
    2726                 :    7704525 : RgX_gcd(GEN x, GEN y)
    2727                 :            : {
    2728                 :            :   long dx, dy;
    2729                 :            :   pari_sp av, av1, lim;
    2730                 :            :   GEN d, g, h, p1, p2, u, v;
    2731                 :    7704525 :   int simple = 0, rational = 1;
    2732                 :            : 
    2733         [ +  + ]:    7704525 :   if (isexactzero(y)) return RgX_copy(x);
    2734         [ +  + ]:    7704325 :   if (isexactzero(x)) return RgX_copy(y);
    2735         [ +  + ]:    7704310 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    2736         [ +  + ]:    2878631 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    2737 [ +  - ][ -  + ]:    2408066 :   if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
    2738                 :            :   {
    2739                 :          0 :     av = avma; u = ggcd(content(x), content(y));
    2740                 :          0 :     return gerepileupto(av, scalarpol(u, varn(x)));
    2741                 :            :   }
    2742         [ +  + ]:    2408066 :   if (rational) return QX_gcd(x,y); /* Q[X] */
    2743                 :            : 
    2744                 :       8685 :   av = avma;
    2745         [ +  + ]:       8685 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    2746                 :            :   else
    2747                 :            :   {
    2748                 :       2485 :     dx = lg(x); dy = lg(y);
    2749         [ +  + ]:       2485 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    2750         [ -  + ]:       2485 :     if (dy==3)
    2751                 :            :     {
    2752                 :          0 :       d = ggcd(gel(y,2), content(x));
    2753                 :          0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    2754                 :            :     }
    2755         [ +  + ]:       2485 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    2756         [ +  + ]:       2485 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    2757                 :       2485 :     d = ggcd(p1,p2);
    2758                 :       2485 :     av1 = avma; lim = stack_lim(av1,1);
    2759                 :       2485 :     g = h = gen_1;
    2760                 :            :     for(;;)
    2761                 :            :     {
    2762                 :       3540 :       GEN r = RgX_pseudorem(u,v);
    2763                 :       3540 :       long degq, du, dv, dr = lg(r);
    2764                 :            : 
    2765         [ +  + ]:       3540 :       if (!signe(r)) break;
    2766         [ +  + ]:       1915 :       if (dr <= 3)
    2767                 :            :       {
    2768                 :        860 :         avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
    2769                 :            :       }
    2770         [ -  + ]:       1055 :       if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
    2771                 :       1055 :       du = lg(u); dv = lg(v); degq = du-dv;
    2772                 :       1055 :       u = v; p1 = g; g = leading_term(u);
    2773      [ +  +  + ]:       1055 :       switch(degq)
    2774                 :            :       {
    2775                 :         40 :         case 0: break;
    2776                 :            :         case 1:
    2777                 :        895 :           p1 = gmul(h,p1); h = g; break;
    2778                 :            :         default:
    2779                 :        120 :           p1 = gmul(gpowgs(h,degq), p1);
    2780                 :        120 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    2781                 :            :       }
    2782                 :       1055 :       v = RgX_Rg_div(r,p1);
    2783         [ -  + ]:       1055 :       if (low_stack(lim, stack_lim(av1,1)))
    2784                 :            :       {
    2785         [ #  # ]:          0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
    2786                 :          0 :         gerepileall(av1,4, &u,&v,&g,&h);
    2787                 :            :       }
    2788                 :       1055 :     }
    2789                 :       1625 :     x = RgX_Rg_mul(primpart(v), d);
    2790                 :            :   }
    2791         [ +  + ]:       7825 :   if (must_negate(x)) x = RgX_neg(x);
    2792                 :    7704525 :   return gerepileupto(av,x);
    2793                 :            : }
    2794                 :            : 
    2795                 :            : static GEN
    2796                 :       2355 : RgX_disc_aux(GEN x)
    2797                 :            : {
    2798                 :       2355 :   long dx = degpol(x), Tx;
    2799                 :            :   GEN D, L, y, p;
    2800 [ +  - ][ -  + ]:       2355 :   if (!signe(x) || !dx) return RgX_get_0(x);
    2801         [ +  + ]:       2355 :   if (dx == 1) return RgX_get_1(x);
    2802         [ +  + ]:       2345 :   if (dx == 2) {
    2803                 :         90 :     GEN a = gel(x,4), b = gel(x,3), c = gel(x,2);
    2804                 :         90 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    2805                 :            :   }
    2806                 :       2255 :   Tx = RgX_simpletype(x);
    2807         [ +  + ]:       2255 :   if (Tx == t_INT) return ZX_disc(x);
    2808         [ -  + ]:         95 :   if (Tx == t_FRAC) return QX_disc(x);
    2809                 :         95 :   p = NULL;
    2810 [ +  + ][ +  - ]:         95 :   if (RgX_is_FpX(x, &p) && p)
    2811                 :         10 :     return Fp_to_mod(FpX_disc(RgX_to_FpX(x,p), p), p);
    2812                 :            : 
    2813                 :         85 :   y = RgX_deriv(x);
    2814         [ -  + ]:         85 :   if (!signe(y)) return RgX_get_0(y);
    2815         [ +  + ]:         85 :   if (Tx == t_REAL)
    2816                 :         10 :     D = resultant2(x,y);
    2817                 :            :   else
    2818                 :            :   {
    2819                 :         75 :     D = RgX_resultant_all(x, y, NULL);
    2820         [ -  + ]:         75 :     if (D == gen_0) return RgX_get_0(y);
    2821                 :            :   }
    2822         [ +  + ]:         85 :   L = leading_term(x); if (!gequal1(L)) D = gdiv(D,L);
    2823         [ +  + ]:         85 :   if (dx & 2) D = gneg(D);
    2824                 :       2355 :   return D;
    2825                 :            : }
    2826                 :            : GEN
    2827                 :        195 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
    2828                 :            : 
    2829                 :            : GEN
    2830                 :       2160 : poldisc0(GEN x, long v)
    2831                 :            : {
    2832                 :            :   long i;
    2833                 :            :   pari_sp av;
    2834                 :            :   GEN z, D;
    2835                 :            : 
    2836   [ +  -  -  -  :       2160 :   switch(typ(x))
                -  -  - ]
    2837                 :            :   {
    2838                 :            :     case t_POL:
    2839                 :       2160 :       av = avma; i = 0;
    2840 [ -  + ][ #  # ]:       2160 :       if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
    2841                 :       2160 :       D = RgX_disc_aux(x);
    2842         [ -  + ]:       2160 :       if (i) D = gsubst(D, MAXVARN, pol_x(0));
    2843                 :       2160 :       return gerepileupto(av, D);
    2844                 :            : 
    2845                 :            :     case t_COMPLEX:
    2846                 :          0 :       return utoineg(4);
    2847                 :            : 
    2848                 :            :     case t_QUAD:
    2849                 :          0 :       return quad_disc(x);
    2850                 :            :     case t_POLMOD:
    2851                 :          0 :       return poldisc0(gel(x,1), v);
    2852                 :            : 
    2853                 :            :     case t_QFR: case t_QFI:
    2854                 :          0 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    2855                 :            : 
    2856                 :            :     case t_VEC: case t_COL: case t_MAT:
    2857                 :          0 :       z = cgetg_copy(x, &i);
    2858         [ #  # ]:          0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    2859                 :          0 :       return z;
    2860                 :            :   }
    2861                 :          0 :   pari_err_TYPE("poldisc",x);
    2862                 :       2160 :   return NULL; /* not reached */
    2863                 :            : }
    2864                 :            : 
    2865                 :            : GEN
    2866                 :         10 : reduceddiscsmith(GEN x)
    2867                 :            : {
    2868                 :         10 :   long j, n = degpol(x);
    2869                 :         10 :   pari_sp av = avma;
    2870                 :            :   GEN xp, M;
    2871                 :            : 
    2872         [ -  + ]:         10 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    2873         [ -  + ]:         10 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    2874                 :         10 :   RgX_check_ZX(x,"poldiscreduced");
    2875         [ -  + ]:         10 :   if (!gequal1(gel(x,n+2)))
    2876                 :          0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    2877                 :         10 :   M = cgetg(n+1,t_MAT);
    2878                 :         10 :   xp = ZX_deriv(x);
    2879         [ +  + ]:         40 :   for (j=1; j<=n; j++)
    2880                 :            :   {
    2881                 :         30 :     gel(M,j) = RgX_to_RgC(xp, n);
    2882         [ +  + ]:         30 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    2883                 :            :   }
    2884                 :         10 :   return gerepileupto(av, ZM_snf(M));
    2885                 :            : }
    2886                 :            : 
    2887                 :            : /***********************************************************************/
    2888                 :            : /**                                                                   **/
    2889                 :            : /**                       STURM ALGORITHM                             **/
    2890                 :            : /**              (number of real roots of x in [a,b])                 **/
    2891                 :            : /**                                                                   **/
    2892                 :            : /***********************************************************************/
    2893                 :            : static int
    2894                 :        135 : exact_sturm(GEN a)
    2895                 :            : {
    2896         [ +  + ]:        135 :   switch(typ(a))
    2897                 :            :   {
    2898                 :        130 :     case t_INT: case t_FRAC: case t_INFINITY: return 1;
    2899                 :        135 :     default: return 0;
    2900                 :            :   }
    2901                 :            : }
    2902                 :            : 
    2903                 :            : /* Deprecated: support the old format: if a (resp. b) is NULL, set it
    2904                 :            :  * to -oo resp. +oo). ZX_sturmpart() should be preferred  */
    2905                 :            : static long
    2906                 :         95 : sturmpart_i(GEN x, GEN a, GEN b)
    2907                 :            : {
    2908                 :            :   long sl, sr, s, t, r1;
    2909                 :         95 :   pari_sp av = avma, lim = stack_lim(av, 1);
    2910                 :            :   int integral;
    2911                 :            :   GEN g,h,u,v;
    2912                 :            : 
    2913         [ -  + ]:         95 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    2914                 :         95 :   t = typ(x);
    2915         [ -  + ]:         95 :   if (t != t_POL)
    2916                 :            :   {
    2917 [ #  # ][ #  # ]:          0 :     if (t == t_INT || t == t_REAL || t == t_FRAC) return 0;
                 [ #  # ]
    2918                 :          0 :     pari_err_TYPE("sturm",x);
    2919                 :            :   }
    2920         [ -  + ]:         95 :   s=lg(x); if (s==3) return 0;
    2921                 :         95 :   u = primpart(x);
    2922                 :         95 :   integral = RgX_is_ZX(u);
    2923 [ +  + ][ +  + ]:         95 :   if (!b && a && typ(a) == t_VEC && lg(a) == 3)
         [ +  - ][ +  - ]
    2924                 :            :   { /* new format */
    2925 [ +  - ][ +  - ]:         50 :     if (integral && exact_sturm(gel(a,1)) && exact_sturm(gel(a,2)))
                 [ +  - ]
    2926                 :         50 :       return ZX_sturmpart(u, a);
    2927                 :            :     /* but can't use new function; convert to old form */
    2928                 :          0 :     integral = 0;
    2929                 :          0 :     b = gel(a,2);
    2930         [ #  # ]:          0 :     if (typ(b) == t_INFINITY)
    2931                 :            :     {
    2932         [ #  # ]:          0 :       if (inf_get_sign(b) < 0) return 0;
    2933                 :          0 :       b = NULL;
    2934                 :            :     }
    2935                 :          0 :     a = gel(a,1);
    2936         [ #  # ]:          0 :     if (typ(a) == t_INFINITY)
    2937                 :            :     {
    2938         [ #  # ]:          0 :       if (inf_get_sign(a) > 0) return 0;
    2939                 :          0 :       a = NULL;
    2940                 :            :     }
    2941                 :            :   }
    2942         [ +  + ]:         45 :   if (integral)
    2943                 :            :   {
    2944         [ +  + ]:         20 :     if (!a) a = mkmoo();
    2945         [ +  + ]:         20 :     if (!b) b = mkoo();
    2946 [ +  + ][ +  - ]:         20 :     if (exact_sturm(a) && exact_sturm(b)) return ZX_sturmpart(u, mkvec2(a,b));
    2947                 :            :   }
    2948                 :            :   /* legacy code: should only be used if we have a t_REAL somewhere; and even
    2949                 :            :    * then, the calling program should be changed */
    2950                 :         30 :   sl = gsigne(leading_term(u));
    2951 [ +  + ][ -  + ]:         30 :   t = a? gsigne(poleval(u,a)): (odd(s)? sl: -sl);
    2952         [ -  + ]:         30 :   if (s==4)
    2953                 :            :   {
    2954         [ #  # ]:          0 :     if (t == 0) return 1;
    2955         [ #  # ]:          0 :     s = b? gsigne(poleval(u,b)):  sl;
    2956                 :          0 :     return (s == t)? 0: 1;
    2957                 :            :   }
    2958         [ +  + ]:         30 :   s = b? gsigne(poleval(u,b)): sl;
    2959                 :         30 :   r1= (t == 0)? 1: 0;
    2960                 :         30 :   v = primpart(RgX_deriv(x));
    2961         [ +  + ]:         30 :   sr = b? gsigne(poleval(v,b)): s;
    2962         [ +  - ]:         30 :   if (sr)
    2963                 :            :   {
    2964         [ +  + ]:         30 :     if (!s) s=sr;
    2965         [ -  + ]:         25 :     else if (sr!=s) { s= -s; r1--; }
    2966                 :            :   }
    2967         [ +  + ]:         30 :   sr = a? gsigne(poleval(v,a)): -t;
    2968         [ +  - ]:         30 :   if (sr)
    2969                 :            :   {
    2970         [ +  + ]:         30 :     if (!t) t=sr;
    2971         [ +  + ]:         25 :     else if (sr!=t) { t= -t; r1++; }
    2972                 :            :   }
    2973                 :         30 :   g=gen_1; h=gen_1;
    2974                 :            :   for(;;)
    2975                 :            :   {
    2976                 :        195 :     GEN p1, r = RgX_pseudorem(u,v);
    2977                 :        195 :     long du=lg(u), dv=lg(v), dr=lg(r), degq=du-dv;
    2978                 :            : 
    2979         [ +  + ]:        195 :     if (dr<=2) pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
    2980 [ +  + ][ +  - ]:        190 :     if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
    2981                 :        190 :     sl = gsigne(gel(r,dr-1));
    2982         [ +  + ]:        190 :     sr = b? gsigne(poleval(r,b)): sl;
    2983         [ +  - ]:        190 :     if (sr)
    2984                 :            :     {
    2985         [ -  + ]:        190 :       if (!s) s=sr;
    2986         [ +  + ]:        190 :       else if (sr!=s) { s= -s; r1--; }
    2987                 :            :     }
    2988 [ +  + ][ +  + ]:        190 :     sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
    2989         [ +  + ]:        190 :     if (sr)
    2990                 :            :     {
    2991         [ -  + ]:        185 :       if (!t) t=sr;
    2992         [ +  + ]:        185 :       else if (sr!=t) { t= -t; r1++; }
    2993                 :            :     }
    2994         [ +  + ]:        190 :     if (dr==3) return r1;
    2995                 :            : 
    2996                 :        165 :     u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
    2997      [ -  +  - ]:        165 :     switch(degq)
    2998                 :            :     {
    2999                 :          0 :       case 0: break;
    3000                 :            :       case 1:
    3001                 :        165 :         p1 = gmul(h,p1); h = g; break;
    3002                 :            :       default:
    3003                 :          0 :         p1 = gmul(gpowgs(h,degq),p1);
    3004                 :          0 :         h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
    3005                 :            :     }
    3006                 :        165 :     v = RgX_Rg_divexact(r,p1);
    3007         [ -  + ]:        165 :     if (low_stack(lim,stack_lim(av,1)))
    3008                 :            :     {
    3009         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polsturm, dr = %ld",dr);
    3010                 :          0 :       gerepileall(av,4,&u,&v,&g,&h);
    3011                 :            :     }
    3012                 :        255 :   }
    3013                 :            : }
    3014                 :            : long
    3015                 :         95 : sturmpart(GEN x, GEN a, GEN b)
    3016                 :            : {
    3017                 :         95 :   pari_sp av = avma;
    3018                 :         95 :   long r = sturmpart_i(x,a,b);
    3019                 :         90 :   avma = av; return r;
    3020                 :            : }
    3021                 :            : long
    3022                 :          0 : RgX_sturmpart(GEN x, GEN ab) { return sturmpart(x, ab, NULL); }
    3023                 :            : 
    3024                 :            : /***********************************************************************/
    3025                 :            : /**                                                                   **/
    3026                 :            : /**                        GENERIC EXTENDED GCD                       **/
    3027                 :            : /**                                                                   **/
    3028                 :            : /***********************************************************************/
    3029                 :            : /* assume typ(x) = typ(y) = t_POL */
    3030                 :            : GEN
    3031                 :      79254 : RgXQ_inv(GEN x, GEN y)
    3032                 :            : {
    3033                 :      79254 :   long vx=varn(x), vy=varn(y);
    3034                 :            :   pari_sp av;
    3035                 :            :   GEN u, v, d;
    3036                 :            : 
    3037         [ -  + ]:      79254 :   while (vx != vy)
    3038                 :            :   {
    3039         [ #  # ]:          0 :     if (varncmp(vx,vy) > 0)
    3040                 :            :     {
    3041         [ #  # ]:          0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3042                 :          0 :       return scalarpol(d, vy);
    3043                 :            :     }
    3044         [ #  # ]:          0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3045                 :          0 :     x = gel(x,2); vx = gvar(x);
    3046                 :            :   }
    3047                 :      79254 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3048         [ +  + ]:      79254 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3049                 :      79249 :   d = gdiv(u,d);
    3050 [ +  + ][ -  + ]:      79249 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3051                 :      79249 :   return gerepileupto(av, d);
    3052                 :            : }
    3053                 :            : 
    3054                 :            : /*Assume x is a polynomial and y is not */
    3055                 :            : static GEN
    3056                 :         80 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3057                 :            : {
    3058                 :         80 :   long vx = varn(x);
    3059                 :         80 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3060 [ +  + ][ +  + ]:         80 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3061         [ +  + ]:         60 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3062                 :         80 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3063                 :            : }
    3064                 :            : /* Assume x==0, y!=0 */
    3065                 :            : static GEN
    3066                 :         45 : zero_bezout(GEN y, GEN *U, GEN *V)
    3067                 :            : {
    3068                 :         45 :   *U=gen_0; *V = ginv(y); return gen_1;
    3069                 :            : }
    3070                 :            : 
    3071                 :            : GEN
    3072                 :        200 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3073                 :            : {
    3074                 :        200 :   long tx=typ(x), ty=typ(y), vx;
    3075 [ +  + ][ +  + ]:        200 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3076         [ +  + ]:        170 :   if (tx != t_POL)
    3077                 :            :   {
    3078         [ +  + ]:        100 :     if (ty == t_POL)
    3079                 :         40 :       return scalar_bezout(y,x,v,u);
    3080                 :            :     else
    3081                 :            :     {
    3082                 :         60 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3083 [ +  + ][ +  + ]:         60 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3084         [ +  + ]:         45 :       if (xis0) return zero_bezout(y,u,v);
    3085                 :         30 :       else return zero_bezout(x,v,u);
    3086                 :            :     }
    3087                 :            :   }
    3088         [ +  + ]:         70 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3089                 :         30 :   vx = varn(x);
    3090         [ -  + ]:         30 :   if (vx != varn(y))
    3091                 :          0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3092         [ #  # ]:          0 :                                    : scalar_bezout(y,x,v,u);
    3093                 :        200 :   return RgX_extgcd(x,y,u,v);
    3094                 :            : }
    3095                 :            : 
    3096                 :            : GEN
    3097                 :        200 : gcdext0(GEN x, GEN y)
    3098                 :            : {
    3099                 :        200 :   GEN z=cgetg(4,t_VEC);
    3100                 :        200 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3101                 :        200 :   return z;
    3102                 :            : }
    3103                 :            : 
    3104                 :            : /*******************************************************************/
    3105                 :            : /*                                                                 */
    3106                 :            : /*                    GENERIC (modular) INVERSE                    */
    3107                 :            : /*                                                                 */
    3108                 :            : /*******************************************************************/
    3109                 :            : 
    3110                 :            : GEN
    3111                 :       1755 : ginvmod(GEN x, GEN y)
    3112                 :            : {
    3113                 :       1755 :   long tx=typ(x);
    3114                 :            : 
    3115      [ +  -  - ]:       1755 :   switch(typ(y))
    3116                 :            :   {
    3117                 :            :     case t_POL:
    3118         [ +  + ]:       1755 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3119         [ +  - ]:       1420 :       if (is_scalar_t(tx)) return ginv(x);
    3120                 :          0 :       break;
    3121                 :            :     case t_INT:
    3122         [ #  # ]:          0 :       if (tx==t_INT) return Fp_inv(x,y);
    3123         [ #  # ]:          0 :       if (tx==t_POL) return gen_0;
    3124                 :            :   }
    3125                 :          0 :   pari_err_TYPE2("ginvmod",x,y);
    3126                 :       1750 :   return NULL; /* not reached */
    3127                 :            : }
    3128                 :            : 
    3129                 :            : /***********************************************************************/
    3130                 :            : /**                                                                   **/
    3131                 :            : /**                         NEWTON POLYGON                            **/
    3132                 :            : /**                                                                   **/
    3133                 :            : /***********************************************************************/
    3134                 :            : 
    3135                 :            : /* assume leading coeff of x is non-zero */
    3136                 :            : GEN
    3137                 :         10 : newtonpoly(GEN x, GEN p)
    3138                 :            : {
    3139                 :            :   GEN y;
    3140                 :            :   long n,ind,a,b,c,u1,u2,r1,r2;
    3141                 :         10 :   long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
    3142                 :            : 
    3143         [ -  + ]:         10 :   if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
    3144         [ -  + ]:         10 :   n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3145                 :         10 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3146                 :         10 :   vval = (long *) pari_malloc(sizeof(long)*(n+1));
    3147         [ +  + ]:         60 :   for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
    3148         [ +  - ]:         10 :   for (a=0, ind=1; a<n; a++)
    3149                 :            :   {
    3150         [ +  - ]:         10 :     if (vval[a] != LONG_MAX) break;
    3151                 :          0 :     gel(y,ind++) = utoipos(LONG_MAX);
    3152                 :            :   }
    3153         [ +  + ]:         30 :   for (b=a+1; b<=n; a=b, b=a+1)
    3154                 :            :   {
    3155         [ -  + ]:         20 :     while (vval[b] == LONG_MAX) b++;
    3156                 :         20 :     u1=vval[a]-vval[b]; u2=b-a;
    3157         [ +  + ]:         70 :     for (c=b+1; c<=n; c++)
    3158                 :            :     {
    3159         [ -  + ]:         50 :       if (vval[c] == LONG_MAX) continue;
    3160                 :         50 :       r1=vval[a]-vval[c]; r2=c-a;
    3161         [ +  + ]:         50 :       if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
    3162                 :            :     }
    3163         [ +  + ]:         60 :     while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
    3164                 :            :   }
    3165                 :         10 :   pari_free(vval); return y;
    3166                 :            : }

Generated by: LCOV version 1.9