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 17242-6994acc) Lines: 1505 1730 87.0 %
Date: 2014-12-22 Functions: 123 133 92.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1150 1590 72.3 %

           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                 :      19684 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
      32                 :            : {
      33                 :      19684 :   long dP=degpol(P), i, k, m;
      34                 :            :   pari_sp av1, av2;
      35                 :            :   GEN s,y,P_lead;
      36                 :            : 
      37         [ -  + ]:      19684 :   if (n<0) pari_err_IMPL("polsym of a negative n");
      38         [ -  + ]:      19684 :   if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
      39         [ -  + ]:      19684 :   if (!signe(P)) pari_err_ROOTS0("polsym");
      40                 :      19684 :   y = cgetg(n+2,t_COL);
      41         [ +  + ]:      19684 :   if (y0)
      42                 :            :   {
      43         [ -  + ]:      10675 :     if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
      44                 :      10675 :     m = lg(y0)-1;
      45         [ +  + ]:      48363 :     for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
      46                 :            :   }
      47                 :            :   else
      48                 :            :   {
      49                 :       9009 :     m = 1;
      50                 :       9009 :     gel(y,1) = stoi(dP);
      51                 :            :   }
      52                 :      19684 :   P += 2; /* strip codewords */
      53                 :            : 
      54         [ +  + ]:      19684 :   P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
      55         [ +  + ]:      19684 :   if (P_lead)
      56                 :            :   {
      57         [ -  + ]:          7 :     if (N) P_lead = Fq_inv(P_lead,T,N);
      58         [ -  + ]:          7 :     else if (T) P_lead = QXQ_inv(P_lead,T);
      59                 :            :   }
      60         [ +  + ]:      53004 :   for (k=m; k<=n; k++)
      61                 :            :   {
      62         [ +  + ]:      33320 :     av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
      63 [ +  + ][ +  + ]:     168574 :     for (i=1; i<k && i<=dP; i++)
      64                 :     135254 :       s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
      65         [ +  + ]:      33320 :     if (N)
      66                 :            :     {
      67                 :      15764 :       s = Fq_red(s, T, N);
      68         [ -  + ]:      15764 :       if (P_lead) s = Fq_mul(s, P_lead, T, N);
      69                 :            :     }
      70         [ +  + ]:      17556 :     else if (T)
      71                 :            :     {
      72                 :       1148 :       s = grem(s, T);
      73         [ -  + ]:       1148 :       if (P_lead) s = grem(gmul(s, P_lead), T);
      74                 :            :     }
      75                 :            :     else
      76         [ +  + ]:      16408 :       if (P_lead) s = gdiv(s, P_lead);
      77                 :      33320 :     av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
      78                 :            :   }
      79                 :      19684 :   return y;
      80                 :            : }
      81                 :            : 
      82                 :            : GEN
      83                 :       4802 : polsym(GEN x, long n)
      84                 :            : {
      85                 :       4802 :   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                 :   92882100 : centermodii(GEN x, GEN p, GEN po2)
      91                 :            : {
      92                 :   92882100 :   GEN y = remii(x, p);
      93   [ +  +  +  - ]:   92882100 :   switch(signe(y))
      94                 :            :   {
      95                 :   25514783 :     case 0: break;
      96 [ +  + ][ +  + ]:   45302880 :     case 1: if (po2 && absi_cmp(y,po2) > 0) y = subii(y, p);
      97                 :   45302880 :       break;
      98 [ +  + ][ +  + ]:   22064437 :     case -1: if (!po2 || absi_cmp(y,po2) > 0) y = addii(y, p);
      99                 :   22064437 :       break;
     100                 :            :   }
     101                 :   92882100 :   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                 :    6183655 : centermod_i(GEN x, GEN p, GEN ps2)
     115                 :            : {
     116                 :            :   long i, lx;
     117                 :            :   pari_sp av;
     118                 :            :   GEN y;
     119                 :            : 
     120         [ +  + ]:    6183655 :   if (!ps2) ps2 = shifti(p,-1);
     121   [ +  +  +  +  :    6183655 :   switch(typ(x))
                   -  - ]
     122                 :            :   {
     123                 :    3634974 :     case t_INT: return centermodii(x,p,ps2);
     124                 :            : 
     125                 :    1931832 :     case t_POL: lx = lg(x);
     126                 :    1931832 :       y = cgetg(lx,t_POL); y[1] = x[1];
     127         [ +  + ]:   13766354 :       for (i=2; i<lx; i++)
     128                 :            :       {
     129                 :   11834522 :         av = avma;
     130                 :   11834522 :         gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
     131                 :            :       }
     132                 :    1931832 :       return normalizepol_lg(y, lx);
     133                 :            : 
     134                 :     616009 :     case t_COL: lx = lg(x);
     135                 :     616009 :       y = cgetg(lx,t_COL);
     136         [ +  + ]:    3081564 :       for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
     137                 :     616009 :       return y;
     138                 :            : 
     139                 :        840 :     case t_MAT: lx = lg(x);
     140                 :        840 :       y = cgetg(lx,t_MAT);
     141         [ +  + ]:      11109 :       for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
     142                 :        840 :       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                 :    6183655 :   return x;
     153                 :            : }
     154                 :            : 
     155                 :            : GEN
     156                 :    4558577 : 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                 :        168 : code(long t1, long t2) { return (t1 << tsh) | t2; }
     171                 :            : void
     172                 :        133 : RgX_type_decode(long x, long *t1, long *t2)
     173                 :            : {
     174                 :        133 :   *t1 = x >> tsh;
     175                 :        133 :   *t2 = (x & ((1L<<tsh)-1));
     176                 :        133 : }
     177                 :            : int
     178                 :     295248 : RgX_type_is_composite(long t) { return t >= tsh; }
     179                 :            : 
     180                 :            : long
     181                 :     298321 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
     182                 :            : {
     183                 :     298321 :   long t[] = {0,0,0,0,0,0,0,0,0,0};
     184                 :     298321 :   long tx = typ(x), lx, i, j, t2 = 0;
     185                 :     298321 :   GEN ff = NULL;
     186                 :     298321 :   *p = *pol = NULL; *pa = LONG_MAX;
     187         [ +  + ]:     298321 :   if (is_scalar_t(tx))
     188                 :            :   {
     189         [ -  + ]:        224 :     if (tx == t_POLMOD) return 0;
     190                 :        224 :     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                 :     298321 :   lx = lg(x);
     204         [ +  + ]:    1417975 :   for (i=2; i<lx; i++)
     205                 :            :   {
     206                 :    1126381 :     GEN c = gel(x,i);
     207   [ +  +  +  +  :    1126381 :     switch(typ(c))
             +  +  +  +  
                      + ]
     208                 :            :     {
     209                 :            :       case t_INT: case t_FRAC:
     210                 :    1115867 :         break;
     211                 :            :       case t_REAL:
     212         [ +  + ]:       1806 :         update_prec(precision(c), pa);
     213                 :       1806 :         t[2]=1; break;
     214                 :            :       case t_INTMOD:
     215 [ +  + ][ -  + ]:        175 :         assign_or_fail(gel(c,1),p);
     216                 :        175 :         t[3]=1; break;
     217                 :            :       case t_FFELT:
     218 [ +  + ][ -  + ]:        427 :         if (!ff) ff=c; else if (!FF_samefield(c,ff)) return 0;
     219 [ +  + ][ -  + ]:        427 :         assign_or_fail(FF_p_i(c),p);
     220                 :        427 :         t[5]=1; break;
     221                 :            :       case t_COMPLEX:
     222         [ +  + ]:         91 :         for (j=1; j<=2; j++)
     223                 :            :         {
     224                 :         63 :           GEN d = gel(c,j);
     225   [ +  +  +  +  :         63 :           switch(typ(d))
                      - ]
     226                 :            :           {
     227                 :            :             case t_INT: case t_FRAC:
     228                 :         35 :               t[4]=1; break;
     229                 :            :             case t_REAL:
     230         [ +  - ]:          7 :               update_prec(precision(d), pa);
     231                 :          7 :               t[6]=1; break;
     232                 :            :             case t_INTMOD:
     233 [ +  - ][ #  # ]:         14 :               assign_or_fail(gel(d,1),p);
     234 [ +  - ][ +  + ]:         14 :               if (!signe(*p) || mod4(*p) != 3) return 0;
     235         [ +  - ]:          7 :               if (!t2) t2 = t_COMPLEX;
     236                 :          7 :               t[3]=1; break;
     237                 :            :             case t_PADIC:
     238         [ +  - ]:          7 :               update_prec(precp(d)+valp(d), pa);
     239 [ +  - ][ #  # ]:          7 :               assign_or_fail(gel(d,2),p);
     240         [ +  - ]:          7 :               if (!t2) t2 = t_COMPLEX;
     241                 :          7 :               t[7]=1; break;
     242                 :          0 :             default: return 0;
     243                 :            :           }
     244                 :            :         }
     245 [ +  + ][ +  - ]:         28 :         if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
                 [ #  # ]
     246                 :         28 :         break;
     247                 :            :       case t_PADIC:
     248         [ +  - ]:         28 :         update_prec(precp(c)+valp(c), pa);
     249 [ +  - ][ #  # ]:         28 :         assign_or_fail(gel(c,2),p);
     250                 :         28 :         t[7]=1; break;
     251                 :            :       case t_QUAD:
     252 [ +  - ][ #  # ]:         28 :         assign_or_fail(gel(c,1),pol);
     253         [ +  + ]:         84 :         for (j=2; j<=3; j++)
     254                 :            :         {
     255                 :         56 :           GEN d = gel(c,j);
     256   [ +  +  +  - ]:         56 :           switch(typ(d))
     257                 :            :           {
     258                 :            :             case t_INT: case t_FRAC:
     259                 :         21 :               t[8]=1; break;
     260                 :            :             case t_INTMOD:
     261 [ +  + ][ -  + ]:         28 :               assign_or_fail(gel(d,1),p);
     262         [ +  - ]:         28 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     263                 :         28 :               t[3]=1; break;
     264                 :            :             case t_PADIC:
     265         [ +  - ]:          7 :               update_prec(precp(d)+valp(d), pa);
     266 [ +  - ][ #  # ]:          7 :               assign_or_fail(gel(d,2),p);
     267         [ +  - ]:          7 :               if (t2 != t_POLMOD) t2 = t_QUAD;
     268                 :          7 :               t[7]=1; break;
     269                 :          0 :             default: return 0;
     270                 :            :           }
     271                 :            :         }
     272                 :         28 :         break;
     273                 :            :       case t_POLMOD:
     274 [ +  + ][ -  + ]:       1295 :         assign_or_fail(gel(c,1),pol);
     275         [ +  + ]:       3885 :         for (j=1; j<=2; j++)
     276                 :            :         {
     277                 :       2590 :           GEN pbis = NULL, polbis = NULL;
     278                 :            :           long pabis;
     279   [ +  +  +  - ]:       2590 :           switch(RgX_type(gel(c,j),&pbis,&polbis,&pabis))
     280                 :            :           {
     281                 :       2569 :             case t_INT: t[9]=1; break;
     282                 :         14 :             case t_INTMOD: t[3]=1; t2 = t_POLMOD; break;
     283         [ +  - ]:          7 :             case t_PADIC: t[7]=1; t2 = t_POLMOD; update_prec(pabis,pa); break;
     284                 :          0 :             default: return 0;
     285                 :            :           }
     286 [ +  + ][ +  - ]:       2590 :           if (pbis) assign_or_fail(pbis,p);
                 [ #  # ]
     287 [ -  + ][ #  # ]:       2590 :           if (polbis) assign_or_fail(polbis,pol);
                 [ #  # ]
     288                 :            :         }
     289                 :       1295 :         break;
     290                 :       6720 :       default: return 0;
     291                 :            :     }
     292                 :            :   }
     293         [ +  + ]:     291594 :   if (t[5]) /* ffelt */
     294                 :            :   {
     295 [ +  - ][ +  - ]:        105 :     if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
         [ +  - ][ +  - ]
         [ +  - ][ -  + ]
     296                 :        105 :     *pol=ff; return t_FFELT;
     297                 :            :   }
     298         [ +  + ]:     291489 :   if (t[6]) /* inexact, complex */
     299                 :            :   {
     300 [ +  - ][ +  - ]:          7 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
         [ +  - ][ -  + ]
     301                 :          7 :     return t_COMPLEX;
     302                 :            :   }
     303         [ +  + ]:     291482 :   if (t[2]) /* inexact, real */
     304                 :            :   {
     305 [ +  - ][ +  - ]:        266 :     if (t2 ||t[3]||t[7]||t[9]) return 0;
         [ +  - ][ -  + ]
     306         [ -  + ]:        266 :     return t[4]?t_COMPLEX:t_REAL;
     307                 :            :   }
     308         [ +  + ]:     291216 :   if (t2) /* polmod/quad/complex of intmod/padic */
     309                 :            :   {
     310         [ +  + ]:         56 :     if (t[3]) return code(t2,t_INTMOD);
     311         [ +  - ]:         21 :     if (t[7]) return code(t2,t_PADIC);
     312                 :            :   }
     313         [ +  + ]:     291160 :   if (t[9]) return code(t_POLMOD,t_INT);
     314         [ +  + ]:     291062 :   if (t[8]) return code(t_QUAD,t_INT);
     315         [ +  + ]:     291055 :   if (t[4]) return code(t_COMPLEX,t_INT);
     316         [ +  + ]:     291048 :   if (t[3]) return t_INTMOD;
     317         [ +  + ]:     290971 :   if (t[7]) return t_PADIC;
     318                 :     298321 :   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                 :      36183 : gp_factor0(GEN x, GEN flag)
     330                 :            : {
     331                 :            :   ulong B;
     332         [ +  + ]:      36183 :   if (!flag) return factor(x);
     333 [ +  - ][ -  + ]:         35 :   if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
     334      [ +  +  - ]:         35 :   switch(lgefint(flag))
     335                 :            :   {
     336                 :          7 :     case 2: B = 0; break;
     337                 :         28 :     case 3: B = flag[2]; break;
     338                 :          0 :     default: pari_err_OVERFLOW("factor [large prime bound]");
     339                 :          0 :              return NULL; /*not reached*/
     340                 :            :   }
     341                 :      36141 :   return boundfact(x, B);
     342                 :            : }
     343                 :            : 
     344                 :            : GEN
     345                 :      38563 : deg1_from_roots(GEN L, long v)
     346                 :            : {
     347                 :      38563 :   long i, l = lg(L);
     348                 :      38563 :   GEN z = cgetg(l,t_COL);
     349         [ +  + ]:      77658 :   for (i=1; i<l; i++)
     350                 :      39095 :     gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
     351                 :      38563 :   return z;
     352                 :            : }
     353                 :            : GEN
     354                 :        609 : roots_from_deg1(GEN x)
     355                 :            : {
     356                 :        609 :   long i,l = lg(x);
     357                 :        609 :   GEN r = cgetg(l,t_VEC);
     358         [ +  + ]:       8813 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
     359                 :        609 :   return r;
     360                 :            : }
     361                 :            : 
     362                 :            : static GEN
     363                 :         35 : gauss_factor_p(GEN p)
     364                 :            : {
     365                 :         35 :   GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
     366                 :         35 :   return mkcomplex(a, b);
     367                 :            : }
     368                 :            : 
     369                 :            : static GEN
     370                 :         42 : gauss_primpart(GEN x, GEN *c)
     371                 :            : {
     372                 :         42 :   GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
     373         [ -  + ]:         42 :   *c = n; if (n == gen_1) return x;
     374                 :         42 :   retmkcomplex(diviiexact(a,n), diviiexact(b,n));
     375                 :            : }
     376                 :            : 
     377                 :            : static GEN
     378                 :         70 : gauss_primpart_try(GEN x, GEN c)
     379                 :            : {
     380                 :            :   GEN r, y;
     381         [ +  + ]:         70 :   if (typ(x) == t_INT)
     382                 :            :   {
     383         [ -  + ]:         42 :     y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
     384                 :            :   }
     385                 :            :   else
     386                 :            :   {
     387                 :         28 :     GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
     388         [ +  + ]:         28 :     gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
     389         [ -  + ]:         14 :     gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
     390                 :            :   }
     391                 :         70 :   return y;
     392                 :            : }
     393                 :            : 
     394                 :            : static int
     395                 :         77 : gauss_cmp(GEN x, GEN y)
     396                 :            : {
     397                 :            :   int v;
     398         [ -  + ]:         77 :   if (typ(x) != t_COMPLEX)
     399         [ #  # ]:          0 :     return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
     400         [ +  + ]:         77 :   if (typ(y) != t_COMPLEX) return 1;
     401                 :         49 :   v = cmpii(gel(x,2), gel(y,2));
     402         [ +  + ]:         49 :   if (v) return v;
     403                 :         77 :   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                 :        112 : gauss_normal(GEN x)
     409                 :            : {
     410 [ +  + ][ +  - ]:        112 :   if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
     411         [ -  + ]:        105 :   if (signe(gel(x,1)) < 0) x = gneg(x);
     412         [ +  + ]:        105 :   if (signe(gel(x,2)) < 0) x = mulcxI(x);
     413                 :        112 :   return x;
     414                 :            : }
     415                 :            : 
     416                 :            : static GEN
     417                 :         42 : gauss_factor(GEN x)
     418                 :            : {
     419                 :         42 :   pari_sp av = avma;
     420                 :         42 :   GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
     421                 :         42 :   long t1 = typ(a);
     422                 :         42 :   long t2 = typ(b), i, j, l, exp = 0;
     423         [ +  + ]:         42 :   if (t1 == t_FRAC) d = gel(a,2);
     424         [ +  + ]:         42 :   if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
     425         [ +  + ]:         42 :   if (d == gen_1) y = x;
     426                 :            :   else
     427                 :            :   {
     428                 :         14 :     y = gmul(x, d);
     429                 :         14 :     a = gel(y,1); t1 = typ(a);
     430                 :         14 :     b = gel(y,2); t2 = typ(b);
     431                 :            :   }
     432 [ +  - ][ -  + ]:         42 :   if (t1 != t_INT || t2 != t_INT) return NULL;
     433                 :         42 :   y = gauss_primpart(y, &n);
     434                 :         42 :   fa = factor(cxnorm(y));
     435                 :         42 :   P = gel(fa,1);
     436                 :         42 :   E = gel(fa,2); l = lg(P);
     437                 :         42 :   P2 = cgetg(l, t_COL);
     438                 :         42 :   E2 = cgetg(l, t_COL);
     439         [ +  + ]:         98 :   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                 :         56 :     GEN p = gel(P,i), w, w2, t, we, pe;
     442                 :         56 :     long v, e = itos(gel(E,i));
     443                 :         56 :     int is2 = equaliu(p, 2);
     444         [ +  + ]:         56 :     w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     445                 :         56 :     w2 = gauss_normal( gconj(w) );
     446                 :            :     /* w * w2 * I^3 = p, w2 = gconj(w) * I */
     447                 :         56 :     pe = powiu(p, e);
     448                 :         56 :     we = gpowgs(w, e);
     449                 :         56 :     t = gauss_primpart_try( gmul(y, gconj(we)), pe );
     450         [ +  + ]:         56 :     if (t) y = t; /* y /= w^e */
     451                 :            :     else {
     452                 :            :       /* y /= conj(w)^e, should be y /= w2^e */
     453                 :         14 :       y = gauss_primpart_try( gmul(y, we), pe );
     454                 :         14 :       swap(w, w2); exp -= e; /* += 3*e mod 4 */
     455                 :            :     }
     456                 :         56 :     gel(P,i) = w;
     457                 :         56 :     v = Z_pvalrem(n, p, &n);
     458         [ +  + ]:         56 :     if (v) {
     459                 :          7 :       exp -= v; /* += 3*v mod 4 */
     460         [ +  - ]:          7 :       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                 :          7 :       gel(E,i) = stoi(e + v);
     466                 :            :     }
     467                 :         56 :     v = Z_pvalrem(d, p, &d);
     468         [ +  + ]:         56 :     if (v) {
     469                 :          7 :       exp += v; /* -= 3*v mod 4 */
     470         [ -  + ]:          7 :       if (is2) v <<= 1; /* 2 is ramified */
     471                 :            :       else {
     472                 :          7 :         gel(P2,j) = w2;
     473                 :          7 :         gel(E2,j) = utoineg(v); j++;
     474                 :            :       }
     475                 :          7 :       gel(E,i) = stoi(e - v);
     476                 :            :     }
     477                 :         56 :     exp &= 3;
     478                 :            :   }
     479         [ +  + ]:         42 :   if (j > 1) {
     480                 :          7 :     long k = 1;
     481                 :          7 :     GEN P1 = cgetg(l, t_COL);
     482                 :          7 :     GEN E1 = cgetg(l, t_COL);
     483                 :            :     /* remove factors with exponent 0 */
     484         [ +  + ]:         14 :     for (i = 1; i < l; i++)
     485         [ -  + ]:          7 :       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                 :          7 :     setlg(P1, k); setlg(E1, k);
     492                 :          7 :     setlg(P2, j); setlg(E2, j);
     493                 :          7 :     fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
     494                 :            :   }
     495 [ +  + ][ +  + ]:         42 :   if (!is_pm1(n) || !is_pm1(d))
     496                 :            :   {
     497                 :         21 :     GEN Fa = factor(gdiv(n, d));
     498                 :         21 :     P = gel(Fa,1); l = lg(P);
     499                 :         21 :     E = gel(Fa,2);
     500         [ +  + ]:         49 :     for (i = 1; i < l; i++)
     501                 :            :     {
     502                 :         28 :       GEN w, p = gel(P,i);
     503                 :            :       long e;
     504                 :            :       int is2;
     505      [ +  +  + ]:         28 :       switch(mod4(p))
     506                 :            :       {
     507                 :         14 :         case 3: continue;
     508                 :          7 :         case 2: is2 = 1; break;
     509                 :          7 :         default:is2 = 0; break;
     510                 :            :       }
     511                 :         14 :       e = itos(gel(E,i));
     512         [ +  + ]:         14 :       w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
     513                 :         14 :       gel(P,i) = w;
     514         [ +  + ]:         14 :       if (is2)
     515                 :          7 :         gel(E,i) = stoi(e << 1);
     516                 :            :       else
     517                 :            :       {
     518                 :          7 :         P = shallowconcat(P, gauss_normal( gconj(w) ));
     519                 :          7 :         E = shallowconcat(E, gel(E,i));
     520                 :            :       }
     521                 :         14 :       exp -= e; /* += 3*e mod 4 */
     522                 :         14 :       exp &= 3;
     523                 :            :     }
     524                 :         21 :     gel(Fa,1) = P;
     525                 :         21 :     gel(Fa,2) = E;
     526                 :         21 :     fa = famat_mul_shallow(fa, Fa);
     527                 :            :   }
     528                 :         42 :   fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
     529                 :            : 
     530                 :         42 :   y = gmul(y, powIs(exp));
     531         [ +  + ]:         42 :   if (!gequal1(y)) {
     532                 :         35 :     gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
     533                 :         35 :     gel(fa,2) = shallowconcat(gen_1,    gel(fa,2));
     534                 :            :   }
     535                 :         42 :   return gerepilecopy(av, fa);
     536                 :            : }
     537                 :            : 
     538                 :            : GEN
     539                 :      38164 : factor(GEN x)
     540                 :            : {
     541                 :      38164 :   long tx=typ(x), lx, i, pa, v, r1;
     542                 :            :   pari_sp av, tetpil;
     543                 :            :   GEN  y, p, p1, p2, pol;
     544                 :            : 
     545         [ +  + ]:      38164 :   if (gequal0(x))
     546         [ +  - ]:         21 :     switch(tx)
     547                 :            :     {
     548                 :            :       case t_INT:
     549                 :            :       case t_COMPLEX:
     550                 :            :       case t_POL:
     551                 :         21 :       case t_RFRAC: return prime_fact(x);
     552                 :          0 :       default: pari_err_TYPE("factor",x);
     553                 :            :     }
     554                 :      38143 :   av = avma;
     555   [ +  +  +  +  :      38143 :   switch(tx)
                   +  - ]
     556                 :            :   {
     557                 :      37450 :     case t_INT: return Z_factor(x);
     558                 :            : 
     559                 :            :     case t_FRAC:
     560                 :        294 :       p1 = Z_factor(gel(x,1));
     561                 :        294 :       p2 = Z_factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
     562                 :        294 :       return gerepilecopy(av, merge_factor_i(p1,p2));
     563                 :            : 
     564                 :            :     case t_POL:
     565                 :        350 :       tx=RgX_type(x,&p,&pol,&pa);
     566   [ +  +  +  +  :        350 :       switch(tx)
             +  +  +  + ]
     567                 :            :       {
     568                 :          7 :         case 0: pari_err_IMPL("factor for general polynomials");
     569                 :        140 :         case t_INT: return QX_factor(x);
     570                 :          7 :         case t_INTMOD: return factmod(x,p);
     571                 :            : 
     572                 :          7 :         case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2;
     573                 :          7 :           av = avma; p1 = roots(x,pa); tetpil = avma;
     574                 :          7 :           p1 = deg1_from_roots(p1, varn(x));
     575                 :          7 :           gel(y,1) = gerepile(av,tetpil,p1);
     576                 :          7 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     577                 :            : 
     578                 :         14 :         case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
     579                 :         14 :           av=avma; p1=cleanroots(x,pa); tetpil=avma;
     580         [ +  + ]:         35 :           for(r1=1; r1<lx; r1++)
     581         [ +  + ]:         28 :             if (typ(gel(p1,r1)) == t_COMPLEX) break;
     582                 :         14 :           lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
     583         [ +  + ]:         35 :           for(i=1; i<r1; i++)
     584                 :         21 :             gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
     585         [ +  + ]:         21 :           for(   ; i<lx; i++)
     586                 :            :           {
     587                 :          7 :             GEN a = gel(p1,2*i-r1);
     588                 :          7 :             p = cgetg(5, t_POL); gel(p2,i) = p;
     589                 :          7 :             p[1] = x[1];
     590                 :          7 :             gel(p,2) = gnorm(a);
     591                 :          7 :             gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
     592                 :          7 :             gel(p,4) = gen_1;
     593                 :            :           }
     594                 :         14 :           gel(y,1) = gerepile(av,tetpil,p2);
     595                 :         14 :           gel(y,2) = const_col(lx-1, gen_1); return y;
     596                 :            : 
     597                 :         14 :         case t_PADIC: return factorpadic(x,p,pa);
     598                 :            : 
     599                 :         35 :         case t_FFELT: return FFX_factor(x,pol);
     600                 :            : 
     601                 :            :         default:
     602                 :            :         {
     603                 :            :           GEN w;
     604                 :            :           long killv, t1, t2;
     605                 :        126 :           x = leafcopy(x); lx=lg(x);
     606                 :        126 :           pol = leafcopy(pol);
     607                 :        126 :           v = pari_var_next_temp();
     608         [ +  + ]:       1001 :           for(i=2; i<lx; i++)
     609                 :            :           {
     610                 :        875 :             p1 = gel(x,i);
     611      [ +  +  + ]:        875 :             switch(typ(p1))
     612                 :            :             {
     613                 :         28 :               case t_QUAD: p1++;
     614                 :            :               case t_COMPLEX:
     615                 :         49 :                 gel(x,i) = mkpolmod(deg1pol_shallow(gel(p1,2), gel(p1,1), v), pol);
     616                 :            :             }
     617                 :            :           }
     618                 :        126 :           killv = (avma != (pari_sp)pol);
     619         [ +  + ]:        126 :           if (killv) setvarn(pol, fetch_var());
     620                 :        126 :           RgX_type_decode(tx, &t1, &t2);
     621      [ +  +  + ]:        126 :           switch (t2)
     622                 :            :           {
     623                 :         70 :             case t_INT: p1 = polfnf(x,pol); break;
     624                 :            :             case t_INTMOD:
     625                 :         35 :               pol = RgX_to_FpX(pol, p);
     626 [ +  + ][ +  + ]:         35 :               if (FpX_is_squarefree(pol,p) && FpX_nbfact(pol, p) == 1)
     627                 :            :               {
     628                 :         21 :                 p1 = factorff(x,p,pol); break;
     629                 :            :               }
     630                 :            :             /*fall through*/
     631                 :         35 :             default: pari_err_IMPL("factor for general polynomial");
     632                 :          0 :               return NULL; /* not reached */
     633                 :            :           }
     634   [ +  +  +  - ]:         91 :           switch (t1)
     635                 :            :           {
     636                 :            :             case t_POLMOD:
     637         [ -  + ]:         63 :               if (killv) (void)delete_var();
     638                 :         63 :               return gerepileupto(av,p1);
     639                 :         14 :             case t_COMPLEX: w = gen_I(); break;
     640                 :         14 :             case t_QUAD: w = mkquad(pol,gen_0,gen_1);
     641                 :         14 :               break;
     642                 :          0 :             default: pari_err_IMPL("factor for general polynomial");
     643                 :          0 :               return NULL; /* not reached */
     644                 :            :           }
     645                 :         28 :           p2=gel(p1,1);
     646         [ +  + ]:         70 :           for(i=1; i<lg(p2); i++)
     647                 :            :           {
     648                 :         42 :             GEN P = gel(p2,i);
     649                 :            :             long j;
     650         [ +  + ]:        140 :             for(j=2; j<lg(P); j++)
     651                 :            :             {
     652                 :         98 :               GEN p = gel(P,j);
     653         [ +  + ]:         98 :               if(typ(p)==t_POLMOD) gel(P,j) = gsubst(gel(p,2),v,w);
     654                 :            :             }
     655                 :            :           }
     656         [ +  - ]:         28 :           if (killv) (void)delete_var();
     657                 :         91 :           return gerepilecopy(av, p1);
     658                 :            :         }
     659                 :            :       }
     660                 :            :     case t_RFRAC: {
     661                 :          7 :       GEN a = gel(x,1), b = gel(x,2);
     662                 :          7 :       y = factor(b); gel(y,2) = gneg_i(gel(y,2));
     663 [ +  - ][ +  - ]:          7 :       if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a), y);
     664                 :          7 :       return gerepilecopy(av, y);
     665                 :            :     }
     666                 :            : 
     667                 :            :     case t_COMPLEX:
     668                 :         42 :       y = gauss_factor(x);
     669         [ +  - ]:         42 :       if (y) return y;
     670                 :            :       /* fall through */
     671                 :            :   }
     672                 :          0 :   pari_err_TYPE("factor",x);
     673                 :      38122 :   return NULL; /* not reached */
     674                 :            : }
     675                 :            : 
     676                 :            : /*******************************************************************/
     677                 :            : /*                                                                 */
     678                 :            : /*                     ROOTS --> MONIC POLYNOMIAL                  */
     679                 :            : /*                                                                 */
     680                 :            : /*******************************************************************/
     681                 :            : static GEN
     682                 :     170686 : normalized_mul(GEN x, GEN y)
     683                 :            : {
     684                 :     170686 :   long a = gel(x,1)[1], b = gel(y,1)[1];
     685                 :     170686 :   return mkvec2(mkvecsmall(a + b),
     686                 :     341372 :                 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                 :      35355 : normalized_to_RgX(GEN L)
     691                 :            : {
     692                 :      35355 :   long i, a = gel(L,1)[1];
     693                 :      35355 :   GEN A = gel(L,2);
     694                 :      35355 :   GEN z = cgetg(a + 3, t_POL);
     695                 :      35355 :   z[1] = evalsigne(1) | evalvarn(varn(A));
     696         [ +  + ]:     436431 :   for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
     697         [ +  + ]:      35362 :   for (     ; i < a+2;   i++) gel(z,i) = gen_0;
     698                 :      35355 :   gel(z,i) = gen_1; return z;
     699                 :            : }
     700                 :            : 
     701                 :            : /* compute prod (x - a[i]) */
     702                 :            : GEN
     703                 :      11760 : roots_to_pol(GEN a, long v)
     704                 :            : {
     705                 :      11760 :   pari_sp av = avma;
     706                 :      11760 :   long i, k, lx = lg(a);
     707                 :            :   GEN L;
     708         [ -  + ]:      11760 :   if (lx == 1) return pol_1(v);
     709                 :      11760 :   L = cgetg(lx, t_VEC);
     710         [ +  + ]:      52059 :   for (k=1,i=1; i<lx-1; i+=2)
     711                 :            :   {
     712                 :      40299 :     GEN s = gel(a,i), t = gel(a,i+1);
     713                 :      40299 :     GEN x0 = gmul(s,t);
     714                 :      40299 :     GEN x1 = gneg(gadd(s,t));
     715                 :      40299 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     716                 :            :   }
     717         [ +  + ]:      11760 :   if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
     718                 :       5091 :                                   scalarpol_shallow(gneg(gel(a,i)), v));
     719                 :      11760 :   setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
     720                 :      11760 :   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                 :      23595 : roots_to_pol_r1(GEN a, long v, long r1)
     726                 :            : {
     727                 :      23595 :   pari_sp av = avma;
     728                 :      23595 :   long i, k, lx = lg(a);
     729                 :            :   GEN L;
     730         [ -  + ]:      23595 :   if (lx == 1) return pol_1(v);
     731                 :      23595 :   L = cgetg(lx, t_VEC);
     732         [ +  + ]:     120869 :   for (k=1,i=1; i<r1; i+=2)
     733                 :            :   {
     734                 :      97274 :     GEN s = gel(a,i), t = gel(a,i+1);
     735                 :      97274 :     GEN x0 = gmul(s,t);
     736                 :      97274 :     GEN x1 = gneg(gadd(s,t));
     737                 :      97274 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     738                 :            :   }
     739         [ +  + ]:      23595 :   if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
     740                 :       5908 :                                     scalarpol_shallow(gneg(gel(a,i)), v));
     741         [ +  + ]:      81064 :   for (i=r1+1; i<lx; i++)
     742                 :            :   {
     743                 :      57469 :     GEN s = gel(a,i);
     744                 :      57469 :     GEN x0 = gnorm(s);
     745                 :      57469 :     GEN x1 = gneg(gtrace(s));
     746                 :      57469 :     gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
     747                 :            :   }
     748                 :      23595 :   setlg(L, k); L = divide_conquer_prod(L, normalized_mul);
     749                 :      23595 :   return gerepileupto(av, normalized_to_RgX(L));
     750                 :            : }
     751                 :            : 
     752                 :            : GEN
     753                 :     930554 : divide_conquer_assoc(GEN x, void *data, GEN (*mul)(void *,GEN,GEN))
     754                 :            : {
     755                 :            :   pari_sp ltop;
     756                 :     930554 :   long i,k,lx = lg(x);
     757                 :            : 
     758         [ +  + ]:     930554 :   if (lx == 1) return gen_1;
     759         [ +  + ]:     929924 :   if (lx == 2) return gcopy(gel(x,1));
     760                 :     869318 :   x = leafcopy(x); k = lx;
     761                 :     869318 :   ltop=avma;
     762         [ +  + ]:    2686016 :   while (k > 2)
     763                 :            :   {
     764         [ -  + ]:    1816698 :     if (DEBUGLEVEL>7)
     765                 :          0 :       err_printf("prod: remaining objects %ld\n",k-1);
     766                 :    1816698 :     lx = k; k = 1;
     767         [ +  + ]:    5747391 :     for (i=1; i<lx-1; i+=2)
     768                 :    3930693 :       gel(x,k++) = mul(data,gel(x,i),gel(x,i+1));
     769         [ +  + ]:    1816698 :     if (i < lx) gel(x,k++) = gel(x,i);
     770         [ +  + ]:    1816698 :     if (gc_needed(ltop,1))
     771                 :          8 :       gerepilecoeffs(ltop,x+1,k-1);
     772                 :            :   }
     773                 :     930554 :   return gel(x,1);
     774                 :            : }
     775                 :            : 
     776                 :            : static GEN
     777                 :    1740706 : _domul(void *data, GEN x, GEN y)
     778                 :            : {
     779                 :    1740706 :   GEN (*mul)(GEN,GEN)=(GEN (*)(GEN,GEN)) data;
     780                 :    1740706 :   return mul(x,y);
     781                 :            : }
     782                 :            : GEN
     783                 :     312614 : divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
     784                 :     312614 : { 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                 :         56 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
     795                 :            : static GEN
     796                 :         91 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
     797                 :            : static GEN
     798                 :        609 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
     799                 :            : static GEN
     800                 :       2346 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
     801                 :            : static GEN
     802                 :       3445 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
     803                 :            : static GEN
     804                 :    1032269 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
     805                 :            : static GEN
     806                 :    1503033 : 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                 :     474488 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
     818                 :            :                GEN (*_pow)(void*,GEN,GEN), void *data)
     819                 :            : {
     820                 :     474488 :   pari_sp av = avma;
     821                 :            :   long k, l, lx;
     822                 :            :   GEN p,x;
     823                 :            : 
     824         [ +  + ]:     474488 :   if (e) /* supplied vector of exponents */
     825                 :     455938 :     p = L;
     826                 :            :   else
     827                 :            :   {
     828      [ +  +  - ]:      18550 :     switch(typ(L)) {
     829                 :            :       case t_VEC:
     830                 :            :       case t_COL: /* product of the L[i] */
     831                 :         49 :         return gerepileupto(av, divide_conquer_assoc(L, data, _mul));
     832                 :            :       case t_MAT: /* genuine factorization */
     833                 :      18501 :         l = lg(L);
     834         [ -  + ]:      18501 :         if (l == 1) return gen_1;
     835         [ +  + ]:      18501 :         if (l == 3) break;
     836                 :            :         /*fall through*/
     837                 :            :       default:
     838                 :          7 :         pari_err_TYPE("factorback [not a factorization]", L);
     839                 :            :     }
     840                 :      18494 :     p = gel(L,1);
     841                 :      18494 :     e = gel(L,2);
     842                 :            :   }
     843                 :            :   /* p = elts, e = expo */
     844                 :     474432 :   lx = lg(p);
     845                 :            :   /* check whether e is an integral vector of correct length */
     846 [ +  - ][ +  - ]:     474432 :   if (!is_vec_t(typ(e)) || lx != lg(e) || !RgV_is_ZV(e))
                 [ +  + ]
     847                 :          7 :     pari_err_TYPE("factorback [not an exponent vector]", e);
     848         [ +  + ]:     474425 :   if (lx == 1) return gen_1;
     849                 :     474383 :   x = cgetg(lx,t_VEC);
     850         [ +  + ]:    1986432 :   for (l=1,k=1; k<lx; k++)
     851         [ +  + ]:    1512049 :     if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
     852                 :     474383 :   x[0] = evaltyp(t_VEC) | _evallg(l);
     853                 :     474474 :   return gerepileupto(av, divide_conquer_assoc(x, data, _mul));
     854                 :            : }
     855                 :            : 
     856                 :            : GEN
     857                 :        574 : idealfactorback(GEN nf, GEN L, GEN e, int red)
     858                 :            : {
     859                 :        574 :   nf = checknf(nf);
     860         [ +  + ]:        574 :   if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
     861                 :        574 :   else     return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
     862                 :            : }
     863                 :            : 
     864                 :            : GEN
     865                 :       1547 : nffactorback(GEN nf, GEN L, GEN e)
     866                 :       1547 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
     867                 :            : 
     868                 :            : GEN
     869                 :     471093 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
     870                 :            : GEN
     871                 :      17759 : factorback(GEN fa) { return factorback2(fa, NULL); }
     872                 :            : 
     873                 :            : static int
     874                 :         63 : RgX_is_irred_i(GEN x)
     875                 :            : {
     876                 :            :   GEN y, p, pol;
     877                 :         63 :   long l = lg(x), pa;
     878                 :            : 
     879 [ +  - ][ -  + ]:         63 :   if (!signe(x) || l <= 3) return 0;
     880   [ +  -  -  + ]:         63 :   switch(RgX_type(x,&p,&pol,&pa))
     881                 :            :   {
     882                 :         14 :     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                 :         49 :   y = factor(x);
     890                 :         63 :   return (lg(gcoeff(y,1,1))==l);
     891                 :            : }
     892                 :            : static int
     893                 :         63 : RgX_is_irred(GEN x)
     894                 :            : {
     895                 :         63 :   pari_sp av = avma;
     896                 :         63 :   int r = RgX_is_irred_i(x);
     897                 :         63 :   avma = av; return r;
     898                 :            : }
     899                 :            : long
     900                 :         63 : isirreducible(GEN x)
     901                 :            : {
     902      [ -  +  - ]:         63 :   switch(typ(x))
     903                 :            :   {
     904                 :          0 :     case t_INT: case t_REAL: case t_FRAC: return 0;
     905                 :         63 :     case t_POL: return RgX_is_irred(x);
     906                 :            :   }
     907                 :          0 :   pari_err_TYPE("isirreducible",x);
     908                 :         63 :   return 0;
     909                 :            : }
     910                 :            : 
     911                 :            : /*******************************************************************/
     912                 :            : /*                                                                 */
     913                 :            : /*                         GENERIC GCD                             */
     914                 :            : /*                                                                 */
     915                 :            : /*******************************************************************/
     916                 :            : /* x is a COMPLEX or a QUAD */
     917                 :            : static GEN
     918                 :        357 : triv_cont_gcd(GEN x, GEN y)
     919                 :            : {
     920                 :        357 :   pari_sp av = avma;
     921                 :            :   GEN c;
     922         [ +  + ]:        357 :   if (typ(x)==t_COMPLEX)
     923                 :            :   {
     924                 :         35 :     GEN a = gel(x,1), b = gel(x,2);
     925 [ +  + ][ -  + ]:         35 :     if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
     926                 :         21 :     c = ggcd(a,b);
     927                 :            :   }
     928                 :            :   else
     929                 :        322 :     c = ggcd(gel(x,2),gel(x,3));
     930                 :        357 :   return gerepileupto(av, ggcd(c,y));
     931                 :            : }
     932                 :            : 
     933                 :            : /* y is a PADIC, x a rational number or an INTMOD */
     934                 :            : static GEN
     935                 :         91 : padic_gcd(GEN x, GEN y)
     936                 :            : {
     937                 :         91 :   GEN p = gel(y,2);
     938                 :         91 :   long v = gvaluation(x,p), w = valp(y);
     939         [ -  + ]:         91 :   if (w < v) v = w;
     940                 :         91 :   return powis(p, v);
     941                 :            : }
     942                 :            : 
     943                 :            : /* x,y in Z[i], at least one of which is t_COMPLEX */
     944                 :            : static GEN
     945                 :         49 : gauss_gcd(GEN x, GEN y)
     946                 :            : {
     947                 :         49 :   pari_sp av = avma;
     948                 :            :   GEN dx, dy;
     949                 :         49 :   dx = denom(x); x = gmul(x, dx);
     950                 :         49 :   dy = denom(y); y = gmul(y, dy);
     951         [ +  + ]:         91 :   while (!gequal0(y))
     952                 :            :   {
     953                 :         42 :     GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
     954                 :         42 :     x = y; y = z;
     955                 :            :   }
     956                 :         49 :   x = gauss_normal(x);
     957         [ +  + ]:         49 :   if (typ(x) == t_COMPLEX)
     958                 :            :   {
     959         [ -  + ]:         35 :     if      (gequal0(gel(x,2))) x = gel(x,1);
     960         [ +  + ]:         35 :     else if (gequal0(gel(x,1))) x = gel(x,2);
     961                 :            :   }
     962                 :         49 :   return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
     963                 :            : }
     964                 :            : 
     965                 :            : static int
     966                 :        105 : is_rational(GEN x) { long t = typ(x); return is_rational_t(t); }
     967                 :            : static int
     968                 :         56 : c_is_rational(GEN x)
     969                 :            : {
     970 [ +  + ][ +  + ]:         56 :   return (is_rational(gel(x,1)) && is_rational(gel(x,2)));
     971                 :            : }
     972                 :            : static GEN
     973                 :         28 : c_zero_gcd(GEN c)
     974                 :            : {
     975                 :         28 :   GEN x = gel(c,1), y = gel(c,2);
     976                 :         28 :   long tx = typ(x), ty = typ(y);
     977 [ +  + ][ -  + ]:         28 :   if (tx == t_REAL || ty == t_REAL) return gen_1;
     978 [ +  - ][ +  + ]:         21 :   if (tx == t_PADIC || tx == t_INTMOD
     979 [ +  - ][ -  + ]:         21 :    || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
     980                 :         28 :   return gauss_gcd(c, gen_0);
     981                 :            : }
     982                 :            : 
     983                 :            : /* y == 0 */
     984                 :            : static GEN
     985                 :    8160872 : zero_gcd(GEN x, long tx)
     986                 :            : {
     987                 :            :   pari_sp av;
     988   [ +  +  +  +  :    8160872 :   switch(tx)
          +  +  +  +  +  
                      + ]
     989                 :            :   {
     990                 :      32473 :     case t_INT: return absi(x);
     991                 :       9473 :     case t_FRAC: return absfrac(x);
     992                 :         28 :     case t_COMPLEX: return c_zero_gcd(x);
     993                 :         42 :     case t_REAL: return gen_1;
     994                 :         28 :     case t_PADIC: return powis(gel(x,2), valp(x));
     995                 :        210 :     case t_SER: return monomial(gen_1, valp(x), varn(x));
     996                 :            :     case t_POLMOD: {
     997                 :       9856 :       GEN d = gel(x,2);
     998 [ +  + ][ +  - ]:       9856 :       if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
     999         [ -  + ]:        609 :       return isinexact(d)? zero_gcd(d, typ(d)): gcopy(d);
    1000                 :            :     }
    1001                 :            :     case t_POL:
    1002         [ +  - ]:    7891513 :       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         [ +  - ]:     217025 :       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                 :    8160872 :   return gcopy(x);
    1016                 :            : }
    1017                 :            : 
    1018                 :            : /* tx = t_RFRAC, y considered as constant */
    1019                 :            : static GEN
    1020                 :     864920 : cont_gcd_rfrac(GEN x, GEN y)
    1021                 :            : {
    1022                 :     864920 :   pari_sp av = avma;
    1023                 :     864920 :   GEN cx; x = primitive_part(x, &cx);
    1024         [ +  + ]:     864920 :   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                 :    2231705 : cont_gcd_pol(GEN x, GEN y)
    1029                 :            : {
    1030                 :    2231705 :   pari_sp av = avma;
    1031                 :    2231705 :   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                 :       9177 : cont_gcd_gen(GEN x, GEN y)
    1036                 :            : {
    1037                 :       9177 :   pari_sp av = avma;
    1038                 :       9177 :   return gerepileupto(av, ggcd(content(x),y));
    1039                 :            : }
    1040                 :            : /* !is_const(tx), y considered as constant */
    1041                 :            : static GEN
    1042                 :    3048994 : cont_gcd(GEN x, long tx, GEN y)
    1043                 :            : {
    1044      [ +  +  + ]:    3048994 :   switch(tx)
    1045                 :            :   {
    1046                 :     808161 :     case t_RFRAC: return cont_gcd_rfrac(x,y);
    1047                 :    2231656 :     case t_POL: return cont_gcd_pol(x,y);
    1048                 :    3048994 :     default: return cont_gcd_gen(x,y);
    1049                 :            :   }
    1050                 :            : }
    1051                 :            : static GEN
    1052                 :     438816 : gcdiq(GEN x, GEN y)
    1053                 :            : {
    1054                 :            :   GEN z;
    1055         [ +  + ]:     438816 :   if (!signe(x)) return Q_abs(y);
    1056                 :     138246 :   z = cgetg(3,t_FRAC);
    1057                 :     138246 :   gel(z,1) = gcdii(x,gel(y,1));
    1058                 :     138246 :   gel(z,2) = icopy(gel(y,2));
    1059                 :     438816 :   return z;
    1060                 :            : }
    1061                 :            : static GEN
    1062                 :     935195 : gcdqq(GEN x, GEN y)
    1063                 :            : {
    1064                 :     935195 :   GEN z = cgetg(3,t_FRAC);
    1065                 :     935195 :   gel(z,1) = gcdii(gel(x,1), gel(y,1));
    1066                 :     935195 :   gel(z,2) = lcmii(gel(x,2), gel(y,2));
    1067                 :     935195 :   return z;
    1068                 :            : }
    1069                 :            : /* assume x,y t_INT or t_FRAC */
    1070                 :            : GEN
    1071                 :   60102157 : Q_gcd(GEN x, GEN y)
    1072                 :            : {
    1073                 :   60102157 :   long tx = typ(x), ty = typ(y);
    1074         [ +  + ]:   60102157 :   if (tx == t_INT)
    1075         [ +  + ]:   58989741 :   { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
    1076                 :            :   else
    1077         [ +  + ]:   60102157 :   { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
    1078                 :            : }
    1079                 :            : 
    1080                 :            : GEN
    1081                 :   21495648 : ggcd(GEN x, GEN y)
    1082                 :            : {
    1083                 :   21495648 :   long l, i, vx, vy, tx = typ(x), ty = typ(y);
    1084                 :            :   pari_sp av, tetpil;
    1085                 :            :   GEN p1,z;
    1086                 :            : 
    1087 [ +  - ][ -  + ]:   21495648 :   if (is_noncalc_t(tx) || is_noncalc_t(ty)) pari_err_TYPE2("gcd",x,y);
    1088         [ -  + ]:   21495648 :   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         [ -  + ]:   21495648 :   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         [ +  + ]:   21495648 :   if (tx>ty) { swap(x,y); lswap(tx,ty); }
    1101                 :            :   /* tx <= ty */
    1102         [ +  + ]:   21495648 :   if (isrationalzero(x)) return zero_gcd(y, ty);
    1103         [ +  + ]:   18467031 :   if (isrationalzero(y)) return zero_gcd(x, tx);
    1104         [ +  + ]:   13345794 :   if (is_const_t(tx))
    1105                 :            :   {
    1106 [ +  + ][ +  +  :    6404398 :     if (ty == tx) switch(tx)
          +  +  +  +  +  
                      - ]
    1107                 :            :     {
    1108                 :            :       case t_INT:
    1109                 :    3197808 :         return gcdii(x,y);
    1110                 :            : 
    1111                 :       2583 :       case t_INTMOD: z=cgetg(3,t_INTMOD);
    1112         [ +  + ]:       2583 :         if (equalii(gel(x,1),gel(y,1)))
    1113                 :       2576 :           gel(z,1) = icopy(gel(x,1));
    1114                 :            :         else
    1115                 :          7 :           gel(z,1) = gcdii(gel(x,1),gel(y,1));
    1116         [ -  + ]:       2583 :         if (gequal1(gel(z,1))) gel(z,2) = gen_0;
    1117                 :            :         else
    1118                 :            :         {
    1119                 :       2583 :           av = avma; p1 = gcdii(gel(z,1),gel(x,2));
    1120         [ +  + ]:       2583 :           if (!is_pm1(p1))
    1121                 :            :           {
    1122                 :        112 :             p1 = gcdii(p1,gel(y,2));
    1123         [ +  + ]:        112 :             if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1124                 :        105 :             else p1 = gerepileuptoint(av, p1);
    1125                 :            :           }
    1126                 :       2583 :           gel(z,2) = p1;
    1127                 :            :         }
    1128                 :       2583 :         return z;
    1129                 :            : 
    1130                 :            :       case t_FRAC:
    1131                 :     116663 :         return gcdqq(x,y);
    1132                 :            : 
    1133                 :            :       case t_FFELT:
    1134         [ -  + ]:      25452 :         if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
    1135 [ -  + ][ #  # ]:      25452 :         return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
    1136                 :            : 
    1137                 :            :       case t_COMPLEX:
    1138 [ +  + ][ +  - ]:         14 :         if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
    1139                 :          7 :         return triv_cont_gcd(y,x);
    1140                 :            : 
    1141                 :            :       case t_PADIC:
    1142         [ +  + ]:         14 :         if (!equalii(gel(x,2),gel(y,2))) return gen_1;
    1143                 :          7 :         return powis(gel(y,2), minss(valp(x), valp(y)));
    1144                 :            : 
    1145                 :            :       case t_QUAD:
    1146                 :        133 :         av=avma; p1=gdiv(x,y);
    1147         [ +  + ]:        133 :         if (gequal0(gel(p1,3)))
    1148                 :            :         {
    1149                 :         14 :           p1=gel(p1,2);
    1150         [ +  + ]:         14 :           if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
    1151                 :          7 :           tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
    1152                 :            :         }
    1153 [ +  + ][ +  + ]:        119 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
    1154                 :        112 :         p1 = ginv(p1); avma=av;
    1155 [ +  + ][ +  - ]:        112 :         if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
    1156                 :        105 :         return triv_cont_gcd(y,x);
    1157                 :            : 
    1158                 :          0 :       default: return gen_1; /* t_REAL */
    1159                 :            :     }
    1160 [ +  + ][ +  +  :    3061731 :     if (is_const_t(ty)) switch(tx)
          +  +  +  +  +  
                      - ]
    1161                 :            :     {
    1162                 :            :       case t_INT:
    1163   [ +  +  +  +  :      38483 :         switch(ty)
             +  +  +  - ]
    1164                 :            :         {
    1165                 :         56 :           case t_INTMOD: z = cgetg(3,t_INTMOD);
    1166                 :         56 :             gel(z,1) = icopy(gel(y,1)); av = avma;
    1167                 :         56 :             p1 = gcdii(gel(y,1),gel(y,2));
    1168         [ +  + ]:         56 :             if (!is_pm1(p1)) {
    1169                 :         35 :               p1 = gcdii(x,p1);
    1170         [ +  + ]:         35 :               if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
    1171                 :            :               else
    1172                 :         28 :                 p1 = gerepileuptoint(av, p1);
    1173                 :            :             }
    1174                 :         56 :             gel(z,2) = p1; return z;
    1175                 :            : 
    1176                 :         35 :           case t_REAL: return gen_1;
    1177                 :            : 
    1178                 :            :           case t_FRAC:
    1179                 :      35158 :             return gcdiq(x,y);
    1180                 :            : 
    1181                 :            :           case t_COMPLEX:
    1182         [ +  - ]:         21 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1183                 :          0 :             return triv_cont_gcd(y,x);
    1184                 :            : 
    1185                 :            :           case t_FFELT:
    1186         [ +  + ]:       3010 :             if (!FF_equal0(y)) return FF_1(y);
    1187         [ +  + ]:        413 :             return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
    1188                 :            : 
    1189                 :            :           case t_PADIC:
    1190                 :         77 :             return padic_gcd(x,y);
    1191                 :            : 
    1192                 :            :           case t_QUAD:
    1193                 :        126 :             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                 :         21 :           case t_PADIC: pari_err_TYPE2("gcd",x,y);
    1204                 :         14 :           default: return gen_1;
    1205                 :            :         }
    1206                 :            : 
    1207                 :            :       case t_INTMOD:
    1208   [ +  +  +  +  :         70 :         switch(ty)
                      - ]
    1209                 :            :         {
    1210                 :            :           case t_FRAC:
    1211                 :         14 :             av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
    1212         [ +  + ]:         14 :             if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
    1213                 :          7 :             return ggcd(gel(y,1), x);
    1214                 :            : 
    1215                 :            :           case t_FFELT:
    1216                 :            :           {
    1217                 :         35 :             GEN p = gel(y,4);
    1218         [ +  + ]:         35 :             if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
    1219         [ +  + ]:         28 :             if (!FF_equal0(y)) return FF_1(y);
    1220         [ -  + ]:          7 :             return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
    1221                 :            :           }
    1222                 :            : 
    1223                 :            :           case t_COMPLEX: case t_QUAD:
    1224                 :         14 :             return triv_cont_gcd(y,x);
    1225                 :            : 
    1226                 :            :           case t_PADIC:
    1227                 :          7 :             return padic_gcd(x,y);
    1228                 :            : 
    1229                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1230                 :            :         }
    1231                 :            : 
    1232                 :            :       case t_FRAC:
    1233   [ +  +  +  +  :        161 :         switch(ty)
                      - ]
    1234                 :            :         {
    1235                 :            :           case t_COMPLEX:
    1236         [ +  + ]:         14 :             if (c_is_rational(y)) return gauss_gcd(x,y);
    1237                 :            :           case t_QUAD:
    1238                 :         84 :             return triv_cont_gcd(y,x);
    1239                 :            :           case t_FFELT:
    1240                 :            :           {
    1241                 :         63 :             GEN p = gel(y,4);
    1242         [ +  + ]:         63 :             if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
    1243         [ +  + ]:         35 :             if (!FF_equal0(y)) return FF_1(y);
    1244         [ +  + ]:         14 :             return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
    1245                 :            :           }
    1246                 :            : 
    1247                 :            :           case t_PADIC:
    1248                 :          7 :             return padic_gcd(x,y);
    1249                 :            : 
    1250                 :          0 :           default: pari_err_TYPE2("gcd",x,y);
    1251                 :            :         }
    1252                 :            :       case t_FFELT:
    1253         [ +  + ]:        105 :         switch(ty)
    1254                 :            :         {
    1255                 :            :           case t_PADIC:
    1256                 :            :           {
    1257                 :         63 :             GEN p = gel(y,2);
    1258                 :         63 :             long v = valp(y);
    1259 [ +  + ][ +  + ]:         63 :             if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
    1260 [ +  + ][ +  + ]:         28 :             return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
    1261                 :            :           }
    1262                 :         42 :           default: pari_err_TYPE2("gcd",x,y);
    1263                 :            :         }
    1264                 :            : 
    1265                 :            :       case t_COMPLEX:
    1266         [ +  - ]:         14 :         switch(ty)
    1267                 :            :         {
    1268                 :            :           case t_PADIC:
    1269                 :         14 :           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         [ +  - ]:          7 :         switch(ty)
    1275                 :            :         {
    1276                 :          7 :           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                 :    3022856 :     return cont_gcd(y,ty, x);
    1283                 :            :   }
    1284                 :            : 
    1285         [ +  + ]:    6941396 :   if (tx == t_POLMOD)
    1286                 :            :   {
    1287         [ +  + ]:       7721 :     if (ty == t_POLMOD)
    1288                 :            :     {
    1289                 :       7385 :       GEN T = gel(x,1);
    1290                 :       7385 :       z = cgetg(3,t_POLMOD);
    1291         [ +  - ]:       7385 :       T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
    1292                 :       7385 :       gel(z,1) = T;
    1293         [ -  + ]:       7385 :       if (degpol(T) <= 0) gel(z,2) = gen_0;
    1294                 :            :       else
    1295                 :            :       {
    1296                 :            :         GEN X, Y, d;
    1297                 :       7385 :         av = avma; X = gel(x,2); Y = gel(y,2);
    1298                 :       7385 :         d = ggcd(content(X), content(Y));
    1299         [ +  + ]:       7385 :         if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
    1300                 :       7385 :         p1 = ggcd(T, X);
    1301                 :       7385 :         gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
    1302                 :            :       }
    1303                 :       7385 :       return z;
    1304                 :            :     }
    1305                 :        336 :     vx = varn(gel(x,1));
    1306      [ +  +  - ]:        336 :     switch(ty)
    1307                 :            :     {
    1308                 :            :       case t_POL:
    1309                 :        308 :         vy = varn(y);
    1310         [ +  + ]:        308 :         if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
    1311                 :        259 :         z = cgetg(3,t_POLMOD);
    1312                 :        259 :         gel(z,1) = RgX_copy(gel(x,1));
    1313                 :        259 :         av = avma; p1 = ggcd(gel(x,1),gel(x,2));
    1314                 :        259 :         gel(z,2) = gerepileupto(av, ggcd(p1,y));
    1315                 :        259 :         return z;
    1316                 :            : 
    1317                 :            :       case t_RFRAC:
    1318                 :         28 :         vy = varn(gel(y,2));
    1319         [ -  + ]:         28 :         if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
    1320                 :         28 :         av = avma;
    1321                 :         28 :         p1 = ggcd(gel(x,1),gel(y,2));
    1322         [ +  + ]:         28 :         if (degpol(p1)) pari_err_OP("gcd",x,y);
    1323                 :         21 :         avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
    1324                 :            :     }
    1325                 :            :   }
    1326                 :            : 
    1327                 :    6933675 :   vx = gvar(x);
    1328                 :    6933675 :   vy = gvar(y);
    1329         [ +  + ]:    6933675 :   if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
    1330         [ +  + ]:    6922062 :   if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
    1331                 :            : 
    1332                 :            :   /* vx = vy: same main variable */
    1333   [ +  +  +  - ]:    6907537 :   switch(tx)
    1334                 :            :   {
    1335                 :            :     case t_POL:
    1336   [ +  +  +  - ]:    6746575 :       switch(ty)
    1337                 :            :       {
    1338                 :    6689802 :         case t_POL: return RgX_gcd(x,y);
    1339                 :            :         case t_SER:
    1340                 :         14 :           z = ggcd(content(x), content(y));
    1341                 :         14 :           return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
    1342                 :      56759 :         case t_RFRAC: return cont_gcd_rfrac(y, x);
    1343                 :            :       }
    1344                 :          0 :       break;
    1345                 :            : 
    1346                 :            :     case t_SER:
    1347                 :         14 :       z = ggcd(content(x), content(y));
    1348      [ +  +  - ]:         14 :       switch(ty)
    1349                 :            :       {
    1350                 :          7 :         case t_SER:   return monomialcopy(z, minss(valp(x),valp(y)), vx);
    1351                 :          7 :         case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
    1352                 :            :       }
    1353                 :          0 :       break;
    1354                 :            : 
    1355                 :            :     case t_RFRAC:
    1356                 :            :     {
    1357                 :     160948 :       GEN xd = gel(x,2), yd = gel(y,2);
    1358         [ -  + ]:     160948 :       if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
    1359                 :     160948 :       z = cgetg(3,t_RFRAC); av = avma;
    1360                 :     160948 :       gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
    1361                 :     160948 :       gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
    1362                 :            :     }
    1363                 :            :   }
    1364                 :          0 :   pari_err_TYPE2("gcd",x,y);
    1365                 :   21495501 :   return NULL; /* not reached */
    1366                 :            : }
    1367                 :            : GEN
    1368         [ +  - ]:       1407 : 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                 :       2856 : scal_lcm(GEN x, GEN y)
    1385                 :            : {
    1386                 :       2856 :   pari_sp av = avma;
    1387                 :       2856 :   long tx = typ(x), ty = typ(y);
    1388         [ -  + ]:       2856 :   if (is_matvec_t(tx)) x = vec_lcm(x);
    1389         [ -  + ]:       2856 :   if (is_matvec_t(ty)) y = vec_lcm(y);
    1390                 :       2856 :   return gerepileupto(av, glcm(x, y));
    1391                 :            : }
    1392                 :            : 
    1393                 :            : static GEN
    1394                 :         49 : fix_lcm(GEN x)
    1395                 :            : {
    1396                 :            :   GEN t;
    1397      [ -  +  - ]:         49 :   switch(typ(x))
    1398                 :            :   {
    1399         [ #  # ]:          0 :     case t_INT: if (signe(x)<0) x = negi(x);
    1400                 :          0 :       break;
    1401                 :            :     case t_POL:
    1402         [ -  + ]:         49 :       if (lg(x) <= 2) break;
    1403                 :         49 :       t = leading_term(x);
    1404 [ +  - ][ -  + ]:         49 :       if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
    1405                 :            :   }
    1406                 :         49 :   return x;
    1407                 :            : }
    1408                 :            : 
    1409                 :            : GEN
    1410                 :       2849 : glcm0(GEN x, GEN y) {
    1411 [ +  + ][ -  + ]:       2849 :   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                 :       2849 :   return gassoc_proto(scal_lcm,x,y);
    1422                 :            : }
    1423                 :            : 
    1424                 :            : GEN
    1425                 :       2912 : glcm(GEN x, GEN y)
    1426                 :            : {
    1427                 :            :   long tx, ty, i, l;
    1428                 :            :   pari_sp av;
    1429                 :            :   GEN p1, z;
    1430                 :            : 
    1431                 :       2912 :   ty = typ(y);
    1432         [ -  + ]:       2912 :   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                 :       2912 :   tx = typ(x);
    1439         [ -  + ]:       2912 :   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 [ +  + ][ +  - ]:       2912 :   if (tx==t_INT && ty==t_INT) return lcmii(x,y);
    1446         [ -  + ]:         49 :   if (gequal0(x)) return gen_0;
    1447                 :            : 
    1448                 :         49 :   av = avma;
    1449         [ +  - ]:         49 :   p1 = ggcd(x,y); if (!gequal1(p1)) y = gdiv(y,p1);
    1450                 :       2912 :   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                 :       2849 : pol_approx0(GEN r, GEN x, int exact)
    1456                 :            : {
    1457                 :            :   long i, lx,lr;
    1458         [ +  - ]:       2849 :   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                 :       2849 :   return 1;
    1464                 :            : }
    1465                 :            : 
    1466                 :            : GEN
    1467                 :       1127 : RgX_gcd_simple(GEN x, GEN y)
    1468                 :            : {
    1469                 :       1127 :   pari_sp av1, av = avma;
    1470                 :       1127 :   GEN r, yorig = y;
    1471 [ +  - ][ +  - ]:       1127 :   int exact = !(isinexactreal(x) || isinexactreal(y));
    1472                 :            : 
    1473                 :            :   for(;;)
    1474                 :            :   {
    1475                 :       2849 :     av1 = avma; r = RgX_rem(x,y);
    1476         [ +  + ]:       2849 :     if (pol_approx0(r, x, exact))
    1477                 :            :     {
    1478                 :       1127 :       avma = av1;
    1479         [ +  + ]:       1127 :       if (y == yorig) return RgX_copy(y);
    1480                 :        910 :       y = normalizepol_approx(y, lg(y));
    1481         [ +  + ]:        910 :       if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
    1482                 :         84 :       return gerepileupto(av,y);
    1483                 :            :     }
    1484                 :       1722 :     x = y; y = r;
    1485         [ -  + ]:       1722 :     if (gc_needed(av,1)) {
    1486         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
    1487                 :          0 :       gerepileall(av,2, &x,&y);
    1488                 :            :     }
    1489                 :       2849 :   }
    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                 :   74423306 : content(GEN x)
    1522                 :            : {
    1523                 :   74423306 :   long lx, i, t, tx = typ(x);
    1524                 :   74423306 :   pari_sp av = avma;
    1525                 :            :   GEN c;
    1526                 :            : 
    1527         [ +  + ]:   74423306 :   if (is_scalar_t(tx)) return zero_gcd(x, tx);
    1528   [ +  +  +  +  :   74412554 :   switch(tx)
                   +  - ]
    1529                 :            :   {
    1530                 :            :     case t_RFRAC:
    1531                 :            :     {
    1532                 :     864934 :       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 [ +  - ][ +  + ]:     864934 :       if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
    1537         [ -  + ]:     709211 :         n = isinexact(n)? zero_gcd(n, typ(n)): gcopy(n);
    1538                 :            :       else
    1539                 :     155723 :         n = content(n);
    1540                 :     864934 :       return gerepileupto(av, gdiv(n, content(d)));
    1541                 :            :     }
    1542                 :            : 
    1543                 :            :     case t_VEC: case t_COL:
    1544         [ -  + ]:    7395367 :       lx = lg(x); if (lx==1) return gen_1;
    1545                 :    7395367 :       break;
    1546                 :            : 
    1547                 :            :     case t_MAT:
    1548                 :            :     {
    1549                 :            :       long hx, j;
    1550                 :        231 :       lx = lg(x);
    1551         [ -  + ]:        231 :       if (lx == 1) return gen_1;
    1552                 :        231 :       hx = lgcols(x);
    1553         [ -  + ]:        231 :       if (hx == 1) return gen_1;
    1554         [ -  + ]:        231 :       if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
    1555         [ -  + ]:        231 :       if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
    1556                 :        231 :       c = content(gel(x,1));
    1557         [ +  + ]:        462 :       for (j=2; j<lx; j++)
    1558         [ +  + ]:        693 :         for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
    1559 [ +  - ][ -  + ]:        231 :       if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
    1560                 :        231 :       return gerepileupto(av,c);
    1561                 :            :     }
    1562                 :            : 
    1563                 :            :     case t_POL: case t_SER:
    1564         [ +  + ]:   66152008 :       lx = lg(x); if (lx == 2) return gen_0;
    1565                 :   66151994 :       break;
    1566                 :            :     case t_QFR: case t_QFI:
    1567                 :         14 :       lx = 4; break;
    1568                 :            : 
    1569                 :          0 :     default: pari_err_TYPE("content",x);
    1570                 :          0 :       return NULL; /* not reached */
    1571                 :            :   }
    1572         [ +  + ]:  209693878 :   for (i=lontyp[tx]; i<lx; i++)
    1573         [ +  + ]:  148972229 :     if (typ(gel(x,i)) != t_INT) break;
    1574                 :   73547375 :   lx--; c = gel(x,lx);
    1575         [ -  + ]:   73547375 :   t = typ(c); if (is_matvec_t(t)) c = content(c);
    1576         [ +  + ]:   73547375 :   if (i > lx)
    1577                 :            :   { /* integer coeffs */
    1578         [ +  + ]:   61430538 :     while (lx-- > lontyp[tx])
    1579                 :            :     {
    1580                 :   59802345 :       c = gcdii(c, gel(x,lx));
    1581         [ +  + ]:   59802345 :       if (is_pm1(c)) { avma=av; return gen_1; }
    1582                 :            :     }
    1583                 :            :   }
    1584                 :            :   else
    1585                 :            :   {
    1586         [ +  + ]:   12825726 :     if (isinexact(c)) c = zero_gcd(c, typ(c));
    1587         [ +  + ]:   31025693 :     while (lx-- > lontyp[tx])
    1588                 :            :     {
    1589                 :   18199967 :       GEN d = gel(x,lx);
    1590         [ -  + ]:   18199967 :       t = typ(d); if (is_matvec_t(t)) d = content(d);
    1591                 :   18199967 :       c = ggcd(c, d);
    1592                 :            :     }
    1593         [ -  + ]:   12825726 :     if (isinexact(c)) { avma=av; return gen_1; }
    1594                 :            :   }
    1595      [ +  -  + ]:   14453919 :   switch(typ(c))
    1596                 :            :   {
    1597                 :            :     case t_INT:
    1598         [ +  + ]:    1632946 :       if (signe(c) < 0) c = negi(c);
    1599                 :    1632946 :       break;
    1600                 :            :     case t_VEC: case t_COL: case t_MAT:
    1601                 :          0 :       pari_err_TYPE("content",x);
    1602                 :            :   }
    1603                 :            : 
    1604         [ +  + ]:   74423306 :   return av==avma? gcopy(c): gerepileupto(av,c);
    1605                 :            : }
    1606                 :            : 
    1607                 :            : GEN
    1608                 :    1357848 : primitive_part(GEN x, GEN *ptc)
    1609                 :            : {
    1610                 :    1357848 :   pari_sp av = avma;
    1611                 :    1357848 :   GEN c = content(x);
    1612         [ +  + ]:    1357848 :   if (gequal1(c)) { avma = av; c = NULL; }
    1613         [ +  + ]:     116001 :   else if (!gequal0(c)) x = gdiv(x,c);
    1614         [ +  + ]:    1357848 :   if (ptc) *ptc = c;
    1615                 :    1357848 :   return x;
    1616                 :            : }
    1617                 :            : GEN
    1618                 :       2716 : 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                 :   83028774 : Q_content(GEN x)
    1624                 :            : {
    1625                 :            :   long i, l;
    1626                 :            :   GEN d;
    1627                 :            :   pari_sp av;
    1628                 :            : 
    1629   [ +  +  +  +  :   83028774 :   switch(typ(x))
                -  +  - ]
    1630                 :            :   {
    1631                 :   66223139 :     case t_INT:  return absi(x);
    1632                 :     976415 :     case t_FRAC: return absfrac(x);
    1633                 :            : 
    1634                 :            :     case t_VEC: case t_COL: case t_MAT:
    1635         [ +  + ]:    7542360 :       l = lg(x); if (l == 1) return gen_1;
    1636                 :    7542241 :       av = avma; d = Q_content(gel(x,1));
    1637         [ +  + ]:   40733617 :       for (i=2; i<l; i++)
    1638                 :            :       {
    1639                 :   33191376 :         d = Q_gcd(d, Q_content(gel(x,i)));
    1640         [ +  + ]:   33191376 :         if ((i & 255) == 0) d = gerepileupto(av, d);
    1641                 :            :       }
    1642                 :    7542241 :       return gerepileupto(av, d);
    1643                 :            : 
    1644                 :            :     case t_POL:
    1645         [ +  + ]:    8286839 :       l = lg(x); if (l == 2) return gen_0;
    1646                 :    8284900 :       av = avma; d = Q_content(gel(x,2));
    1647         [ +  + ]:   31838634 :       for (i=3; i<l; i++) d = Q_gcd(d, Q_content(gel(x,i)));
    1648                 :    8284900 :       return gerepileupto(av, d);
    1649                 :          0 :     case t_POLMOD: return Q_content(gel(x,2));
    1650                 :         21 :     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                 :   83028774 :   return NULL; /* not reached */
    1654                 :            : }
    1655                 :            : 
    1656                 :            : GEN
    1657                 :       5852 : ZX_content(GEN x)
    1658                 :            : {
    1659                 :       5852 :   long i, l = lg(x);
    1660                 :            :   GEN d;
    1661                 :            :   pari_sp av;
    1662                 :            : 
    1663         [ -  + ]:       5852 :   if (l == 2) return gen_0;
    1664                 :       5852 :   d = gel(x,2);
    1665         [ +  + ]:       5852 :   if (l == 3) return absi(d);
    1666                 :       3682 :   av = avma;
    1667 [ +  + ][ +  + ]:      30905 :   for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
    1668         [ -  + ]:       3682 :   if (signe(d) < 0) d = absi(d);
    1669                 :       5852 :   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                 :   11349277 : Q_denom(GEN x)
    1677                 :            : {
    1678                 :            :   long i, l;
    1679                 :            :   GEN d, D;
    1680                 :            :   pari_sp av;
    1681                 :            : 
    1682   [ +  +  +  +  :   11349277 :   switch(typ(x))
                   +  - ]
    1683                 :            :   {
    1684                 :    6659918 :     case t_INT: return gen_1;
    1685                 :    2759793 :     case t_FRAC: return gel(x,2);
    1686                 :            : 
    1687                 :            :     case t_VEC: case t_COL: case t_MAT:
    1688         [ +  + ]:    1625034 :       l = lg(x); if (l == 1) return gen_1;
    1689                 :    1623746 :       av = avma; d = Q_denom(gel(x,1));
    1690         [ +  + ]:    7740695 :       for (i=2; i<l; i++)
    1691                 :            :       {
    1692                 :    6116949 :         D = Q_denom(gel(x,i));
    1693         [ +  + ]:    6116949 :         if (D != gen_1) d = lcmii(d, D);
    1694         [ -  + ]:    6116949 :         if ((i & 255) == 0) d = gerepileuptoint(av, d);
    1695                 :            :       }
    1696                 :    1623746 :       return gerepileuptoint(av, d);
    1697                 :            : 
    1698                 :            :     case t_POL:
    1699         [ +  + ]:     303734 :       l = lg(x); if (l == 2) return gen_1;
    1700                 :     303244 :       av = avma; d = Q_denom(gel(x,2));
    1701         [ +  + ]:    1343165 :       for (i=3; i<l; i++)
    1702                 :            :       {
    1703                 :    1039921 :         D = Q_denom(gel(x,i));
    1704         [ +  + ]:    1039921 :         if (D != gen_1) d = lcmii(d, D);
    1705                 :            :       }
    1706                 :     303244 :       return gerepileuptoint(av, d);
    1707                 :            : 
    1708                 :        798 :     case t_POLMOD: return Q_denom(gel(x,2));
    1709                 :            :   }
    1710                 :          0 :   pari_err_TYPE("Q_denom",x);
    1711                 :   11349277 :   return NULL; /* not reached */
    1712                 :            : }
    1713                 :            : 
    1714                 :            : GEN
    1715                 :    1322622 : Q_remove_denom(GEN x, GEN *ptd)
    1716                 :            : {
    1717                 :    1322622 :   GEN d = Q_denom(x);
    1718         [ +  + ]:    1322622 :   if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
    1719         [ +  + ]:    1322622 :   if (ptd) *ptd = d;
    1720                 :    1322622 :   return x;
    1721                 :            : }
    1722                 :            : 
    1723                 :            : /* return y = x * d, assuming x rational, and d,y integral */
    1724                 :            : GEN
    1725                 :    8950571 : Q_muli_to_int(GEN x, GEN d)
    1726                 :            : {
    1727                 :            :   long i, l;
    1728                 :            :   GEN y, xn, xd;
    1729                 :            :   pari_sp av;
    1730                 :            : 
    1731         [ -  + ]:    8950571 :   if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
    1732   [ +  +  +  +  :    8950571 :   switch (typ(x))
                   -  - ]
    1733                 :            :   {
    1734                 :            :     case t_INT:
    1735                 :    3596658 :       return mulii(x,d);
    1736                 :            : 
    1737                 :            :     case t_FRAC:
    1738                 :    3869994 :       xn = gel(x,1);
    1739                 :    3869994 :       xd = gel(x,2); av = avma;
    1740                 :    3869994 :       y = mulii(xn, diviiexact(d, xd));
    1741                 :    3869994 :       return gerepileuptoint(av, y);
    1742                 :            : 
    1743                 :            :     case t_VEC: case t_COL: case t_MAT:
    1744                 :     950043 :       y = cgetg_copy(x, &l);
    1745         [ +  + ]:    6058927 :       for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1746                 :     950043 :       return y;
    1747                 :            : 
    1748                 :            :     case t_POL:
    1749                 :     533876 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1750         [ +  + ]:    3544771 :       for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
    1751                 :     533876 :       return y;
    1752                 :            : 
    1753                 :            :     case t_POLMOD:
    1754                 :          0 :       y = cgetg(3, t_POLMOD);
    1755                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1756                 :          0 :       gel(y,2) = Q_muli_to_int(gel(x,2), d);
    1757                 :          0 :       return y;
    1758                 :            :   }
    1759                 :          0 :   pari_err_TYPE("Q_muli_to_int",x);
    1760                 :    8950571 :   return NULL; /* not reached */
    1761                 :            : }
    1762                 :            : 
    1763                 :            : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
    1764                 :            : static GEN
    1765                 :      19903 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
    1766                 :            : {
    1767                 :            :   long i, l;
    1768                 :            :   GEN y, xn, xd;
    1769                 :            :   pari_sp av;
    1770                 :            : 
    1771   [ +  +  +  +  :      19903 :   switch(typ(x))
                   -  - ]
    1772                 :            :   {
    1773                 :            :     case t_INT:
    1774                 :       2730 :       av = avma; y = diviiexact(x,d);
    1775                 :       2730 :       return gerepileuptoint(av, mulii(y,n));
    1776                 :            : 
    1777                 :            :     case t_FRAC:
    1778                 :      11283 :       xn = gel(x,1);
    1779                 :      11283 :       xd = gel(x,2); av = avma;
    1780                 :      11283 :       y = mulii(diviiexact(xn, d), diviiexact(n, xd));
    1781                 :      11283 :       return gerepileuptoint(av, y);
    1782                 :            : 
    1783                 :            :     case t_VEC: case t_COL: case t_MAT:
    1784                 :       4266 :       y = cgetg_copy(x, &l);
    1785         [ +  + ]:      15242 :       for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1786                 :       4266 :       return y;
    1787                 :            : 
    1788                 :            :     case t_POL:
    1789                 :       1624 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1790         [ +  + ]:       6706 :       for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
    1791                 :       1624 :       return y;
    1792                 :            : 
    1793                 :            :     case t_POLMOD:
    1794                 :          0 :       y = cgetg(3, t_POLMOD);
    1795                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1796                 :          0 :       gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
    1797                 :          0 :       return y;
    1798                 :            :   }
    1799                 :          0 :   pari_err_TYPE("Q_divmuli_to_int",x);
    1800                 :      19903 :   return NULL; /* not reached */
    1801                 :            : }
    1802                 :            : 
    1803                 :            : /* return x / d. x: rational; d,result: integral. */
    1804                 :            : static GEN
    1805                 :   12249380 : Q_divi_to_int(GEN x, GEN d)
    1806                 :            : {
    1807                 :            :   long i, l;
    1808                 :            :   GEN y;
    1809                 :            : 
    1810   [ +  +  +  -  :   12249380 :   switch(typ(x))
                      - ]
    1811                 :            :   {
    1812                 :            :     case t_INT:
    1813                 :   10111565 :       return diviiexact(x,d);
    1814                 :            : 
    1815                 :            :     case t_VEC: case t_COL: case t_MAT:
    1816                 :     991044 :       y = cgetg_copy(x, &l);
    1817         [ +  + ]:    6799743 :       for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1818                 :     991044 :       return y;
    1819                 :            : 
    1820                 :            :     case t_POL:
    1821                 :    1146771 :       y = cgetg_copy(x, &l); y[1] = x[1];
    1822         [ +  + ]:    6015558 :       for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
    1823                 :    1146771 :       return y;
    1824                 :            : 
    1825                 :            :     case t_POLMOD:
    1826                 :          0 :       y = cgetg(3, t_POLMOD);
    1827                 :          0 :       gel(y,1) = RgX_copy(gel(x,1));
    1828                 :          0 :       gel(y,2) = Q_divi_to_int(gel(x,2), d);
    1829                 :          0 :       return y;
    1830                 :            :   }
    1831                 :          0 :   pari_err_TYPE("Q_divi_to_int",x);
    1832                 :   12249380 :   return NULL; /* not reached */
    1833                 :            : }
    1834                 :            : /* c t_FRAC */
    1835                 :            : static GEN
    1836                 :     157960 : Q_divq_to_int(GEN x, GEN c)
    1837                 :            : {
    1838                 :     157960 :   GEN n = gel(c,1), d = gel(c,2);
    1839         [ +  + ]:     157960 :   if (is_pm1(n)) {
    1840                 :     154115 :     GEN y = Q_muli_to_int(x,d);
    1841         [ -  + ]:     154115 :     if (signe(n) < 0) y = gneg(y);
    1842                 :     154115 :     return y;
    1843                 :            :   }
    1844                 :     157960 :   return Q_divmuli_to_int(x, n,d);
    1845                 :            : }
    1846                 :            : 
    1847                 :            : /* return y = x / c, assuming x,c rational, and y integral */
    1848                 :            : GEN
    1849                 :       3668 : Q_div_to_int(GEN x, GEN c)
    1850                 :            : {
    1851      [ +  +  - ]:       3668 :   switch(typ(c))
    1852                 :            :   {
    1853                 :       3626 :     case t_INT:  return Q_divi_to_int(x, c);
    1854                 :         42 :     case t_FRAC: return Q_divq_to_int(x, c);
    1855                 :            :   }
    1856                 :          0 :   pari_err_TYPE("Q_div_to_int",c);
    1857                 :       3668 :   return NULL; /* not reached */
    1858                 :            : }
    1859                 :            : /* return y = x * c, assuming x,c rational, and y integral */
    1860                 :            : GEN
    1861                 :          0 : Q_mul_to_int(GEN x, GEN c)
    1862                 :            : {
    1863                 :            :   GEN d, n;
    1864      [ #  #  # ]:          0 :   switch(typ(c))
    1865                 :            :   {
    1866                 :          0 :     case t_INT: return Q_muli_to_int(x, c);
    1867                 :            :     case t_FRAC:
    1868                 :          0 :       n = gel(c,1);
    1869                 :          0 :       d = gel(c,2);
    1870                 :          0 :       return Q_divmuli_to_int(x, d,n);
    1871                 :            :   }
    1872                 :          0 :   pari_err_TYPE("Q_mul_to_int",c);
    1873                 :          0 :   return NULL; /* not reached */
    1874                 :            : }
    1875                 :            : 
    1876                 :            : GEN
    1877                 :   10279529 : Q_primitive_part(GEN x, GEN *ptc)
    1878                 :            : {
    1879                 :   10279529 :   pari_sp av = avma;
    1880                 :   10279529 :   GEN c = Q_content(x);
    1881         [ +  + ]:   10279529 :   if (typ(c) == t_INT)
    1882                 :            :   {
    1883         [ +  + ]:   10121611 :     if (is_pm1(c)) { avma = av; c = NULL; }
    1884         [ +  + ]:    1568352 :     else if (signe(c)) x = Q_divi_to_int(x, c);
    1885                 :            :   }
    1886                 :     157918 :   else x = Q_divq_to_int(x, c);
    1887         [ +  + ]:   10279529 :   if (ptc) *ptc = c;
    1888                 :   10279529 :   return x;
    1889                 :            : }
    1890                 :            : GEN
    1891                 :    1174640 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
    1892                 :            : 
    1893                 :            : /*******************************************************************/
    1894                 :            : /*                                                                 */
    1895                 :            : /*                           SUBRESULTANT                          */
    1896                 :            : /*                                                                 */
    1897                 :            : /*******************************************************************/
    1898                 :            : /* for internal use */
    1899                 :            : GEN
    1900                 :    2747891 : gdivexact(GEN x, GEN y)
    1901                 :            : {
    1902                 :            :   long i,lx;
    1903                 :            :   GEN z;
    1904         [ +  + ]:    2747891 :   if (gequal1(y)) return x;
    1905   [ +  +  +  +  :    2727787 :   switch(typ(x))
                      + ]
    1906                 :            :   {
    1907                 :            :     case t_INT:
    1908         [ +  + ]:    2180188 :       if (typ(y)==t_INT) return diviiexact(x,y);
    1909         [ +  - ]:        224 :       if (!signe(x)) return gen_0;
    1910                 :          0 :       break;
    1911                 :            :     case t_INTMOD:
    1912                 :      12687 :     case t_POLMOD: return gmul(x,ginv(y));
    1913                 :            :     case t_POL:
    1914      [ -  +  + ]:     534639 :       switch(typ(y))
    1915                 :            :       {
    1916                 :            :         case t_INTMOD:
    1917                 :          0 :         case t_POLMOD: return gmul(x,ginv(y));
    1918                 :            :         case t_POL: { /* not stack-clean */
    1919                 :            :           long v;
    1920         [ +  + ]:      24864 :           if (varn(x)!=varn(y)) break;
    1921                 :      24017 :           v = RgX_valrem(y,&y);
    1922         [ +  + ]:      24017 :           if (v) x = RgX_shift_shallow(x,-v);
    1923         [ +  + ]:      24017 :           if (!degpol(y)) { y = gel(y,2); break; }
    1924                 :      20545 :           return RgX_div(x,y);
    1925                 :            :         }
    1926                 :            :       }
    1927                 :     514094 :       return RgX_Rg_divexact(x, y);
    1928                 :            :     case t_VEC: case t_COL: case t_MAT:
    1929                 :          7 :       lx = lg(x); z = new_chunk(lx);
    1930         [ +  + ]:         49 :       for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
    1931                 :          7 :       z[0] = x[0]; return z;
    1932                 :            :   }
    1933         [ -  + ]:        266 :   if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
    1934                 :    2747891 :   return gdiv(x,y);
    1935                 :            : }
    1936                 :            : 
    1937                 :            : static GEN
    1938                 :      39367 : init_resultant(GEN x, GEN y)
    1939                 :            : {
    1940                 :      39367 :   long tx = typ(x), ty = typ(y), vx, vy;
    1941 [ +  - ][ -  + ]:      39367 :   if (is_scalar_t(tx) || is_scalar_t(ty))
    1942                 :            :   {
    1943 [ #  # ][ #  # ]:          0 :     if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
    1944         [ #  # ]:          0 :     if (tx==t_POL) return gpowgs(y, degpol(x));
    1945         [ #  # ]:          0 :     if (ty==t_POL) return gpowgs(x, degpol(y));
    1946                 :          0 :     return gen_1;
    1947                 :            :   }
    1948         [ -  + ]:      39367 :   if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
    1949         [ -  + ]:      39367 :   if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
    1950 [ +  - ][ +  + ]:      39367 :   if (!signe(x) || !signe(y)) return gmul(RgX_get_0(x),RgX_get_0(y)); /*type*/
    1951                 :      39290 :   vx = varn(x);
    1952         [ +  - ]:      39290 :   vy = varn(y); if (vx == vy) return NULL;
    1953         [ #  # ]:      39367 :   return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
    1954                 :            : }
    1955                 :            : 
    1956                 :            : static long
    1957                 :      81814 : RgX_simpletype(GEN x)
    1958                 :            : {
    1959                 :      81814 :   long T = t_INT, i, lx = lg(x);
    1960         [ +  + ]:     525012 :   for (i = 2; i < lx; i++)
    1961                 :            :   {
    1962                 :     443212 :     GEN c = gel(x,i);
    1963                 :     443212 :     long tc = typ(c);
    1964      [ +  +  + ]:     443212 :     switch(tc) {
    1965                 :            :       case t_INT:
    1966                 :     393266 :         break;
    1967                 :            :       case t_FRAC:
    1968         [ +  + ]:      25320 :         if (T == t_INT) T = t_FRAC;
    1969                 :      25320 :         break;
    1970                 :            :       default:
    1971         [ +  + ]:      24626 :         if (isinexact(c)) return t_REAL;
    1972                 :      24612 :         T = 0; break;
    1973                 :            :     }
    1974                 :            :   }
    1975                 :      81814 :   return T;
    1976                 :            : }
    1977                 :            : 
    1978                 :            : /* x an RgX, y a scalar */
    1979                 :            : static GEN
    1980                 :          0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
    1981                 :            : {
    1982                 :          0 :   *V = gpowgs(y,degpol(x)-1);
    1983                 :          0 :   *U = gen_0; return gmul(y, *V);
    1984                 :            : }
    1985                 :            : 
    1986                 :            : /* return 0 if the subresultant chain can be interrupted.
    1987                 :            :  * Set u = NULL if the resultant is 0. */
    1988                 :            : static int
    1989                 :     202199 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
    1990                 :            : {
    1991                 :     202199 :   GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
    1992                 :            :   long du, dv, dr, degq;
    1993                 :            : 
    1994         [ +  + ]:     202199 :   dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
    1995                 :     202066 :   du = degpol(*u);
    1996                 :     202066 :   dv = degpol(*v);
    1997                 :     202066 :   degq = du - dv;
    1998         [ +  + ]:     202066 :   if (*um1 == gen_1)
    1999                 :     119204 :     u0 = gpowgs(gel(*v,dv+2),degq+1);
    2000         [ +  + ]:      82862 :   else if (*um1 == gen_0)
    2001                 :      37453 :     u0 = gen_0;
    2002                 :            :   else /* except in those 2 cases, um1 is an RgX */
    2003                 :      45409 :     u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
    2004                 :            : 
    2005         [ +  + ]:     202066 :   if (*uze == gen_0) /* except in that case, uze is an RgX */
    2006                 :     119204 :     u0 = scalarpol(u0, varn(*u)); /* now an RgX */
    2007                 :            :   else
    2008                 :      82862 :     u0 = gsub(u0, gmul(q,*uze));
    2009                 :            : 
    2010                 :     202066 :   *um1 = *uze;
    2011                 :     202066 :   *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
    2012                 :            : 
    2013                 :     202066 :   *u = *v; c = *g; *g  = leading_term(*u);
    2014      [ +  +  + ]:     202066 :   switch(degq)
    2015                 :            :   {
    2016                 :         49 :     case 0: break;
    2017                 :            :     case 1:
    2018                 :     179766 :       c = gmul(*h,c); *h = *g; break;
    2019                 :            :     default:
    2020                 :      22251 :       c = gmul(gpowgs(*h,degq), c);
    2021                 :      22251 :       *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
    2022                 :            :   }
    2023                 :     202066 :   *v  = RgX_Rg_divexact(r,c);
    2024                 :     202066 :   *uze= RgX_Rg_divexact(*uze,c);
    2025         [ +  + ]:     202066 :   if (both_odd(du, dv)) *signh = -*signh;
    2026                 :     202199 :   return (dr > 3);
    2027                 :            : }
    2028                 :            : 
    2029                 :            : /* compute U, V s.t Ux + Vy = resultant(x,y) */
    2030                 :            : static GEN
    2031                 :     119624 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
    2032                 :            : {
    2033                 :            :   pari_sp av, av2;
    2034                 :     119624 :   long dx, dy, du, signh, tx = typ(x), ty = typ(y);
    2035                 :            :   GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
    2036                 :            : 
    2037         [ -  + ]:     119624 :   if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
    2038         [ -  + ]:     119624 :   if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
    2039 [ +  + ][ -  + ]:     119624 :   if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
    2040         [ -  + ]:     119617 :   if (tx != t_POL) {
    2041         [ #  # ]:          0 :     if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
    2042                 :          0 :     return scalar_res(y,x,V,U);
    2043                 :            :   }
    2044         [ -  + ]:     119617 :   if (ty != t_POL) return scalar_res(x,y,U,V);
    2045         [ -  + ]:     119617 :   if (varn(x) != varn(y))
    2046                 :          0 :     return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
    2047         [ #  # ]:          0 :                                         : scalar_res(y,x,V,U);
    2048                 :     119617 :   dx = degpol(x); dy = degpol(y); signh = 1;
    2049         [ +  + ]:     119617 :   if (dx < dy)
    2050                 :            :   {
    2051                 :     119057 :     pswap(U,V); lswap(dx,dy); swap(x,y);
    2052         [ +  + ]:     119057 :     if (both_odd(dx, dy)) signh = -signh;
    2053                 :            :   }
    2054         [ +  + ]:     119617 :   if (dy == 0)
    2055                 :            :   {
    2056                 :        518 :     *V = gpowgs(gel(y,2),dx-1);
    2057                 :        518 :     *U = gen_0; return gmul(*V,gel(y,2));
    2058                 :            :   }
    2059                 :     119099 :   av = avma;
    2060                 :     119099 :   u = x = primitive_part(x, &cu);
    2061                 :     119099 :   v = y = primitive_part(y, &cv);
    2062                 :     119099 :   g = h = gen_1; av2 = avma;
    2063                 :     119099 :   um1 = gen_1; uze = gen_0;
    2064                 :            :   for(;;)
    2065                 :            :   {
    2066         [ +  + ]:     201891 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2067         [ -  + ]:      82792 :     if (gc_needed(av2,1))
    2068                 :            :     {
    2069         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
    2070                 :          0 :       gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
    2071                 :            :     }
    2072                 :      82792 :   }
    2073                 :            :   /* uze an RgX */
    2074         [ +  + ]:     119099 :   if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
    2075                 :     119085 :   z = gel(v,2); du = degpol(u);
    2076         [ +  + ]:     119085 :   if (du > 1)
    2077                 :            :   { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
    2078                 :       3716 :     p1 = gpowgs(gdiv(z,h),du-1);
    2079                 :       3716 :     z = gmul(z,p1);
    2080                 :       3716 :     uze = RgX_Rg_mul(uze, p1);
    2081                 :            :   }
    2082         [ +  + ]:     119085 :   if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
    2083                 :            : 
    2084                 :     119085 :   vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
    2085         [ -  + ]:     119085 :   if (signe(r)) pari_warn(warner,"inexact computation in subresext");
    2086                 :            :   /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
    2087                 :     119085 :   p1 = gen_1;
    2088         [ +  + ]:     119085 :   if (cu) p1 = gmul(p1, gpowgs(cu,dy));
    2089         [ +  + ]:     119085 :   if (cv) p1 = gmul(p1, gpowgs(cv,dx));
    2090         [ +  + ]:     119085 :   cu = cu? gdiv(p1,cu): p1;
    2091         [ +  + ]:     119085 :   cv = cv? gdiv(p1,cv): p1;
    2092                 :     119085 :   z = gmul(z,p1);
    2093                 :     119085 :   *U = RgX_Rg_mul(uze,cu);
    2094                 :     119085 :   *V = RgX_Rg_mul(vze,cv);
    2095                 :     119624 :   return z;
    2096                 :            : }
    2097                 :            : GEN
    2098                 :          0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
    2099                 :            : {
    2100                 :          0 :   pari_sp av = avma;
    2101                 :          0 :   GEN z = subresext_i(x, y, U, V);
    2102                 :          0 :   gerepileall(av, 3, &z, U, V);
    2103                 :          0 :   return z;
    2104                 :            : }
    2105                 :            : 
    2106                 :            : static GEN
    2107                 :         28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
    2108                 :            : {
    2109                 :         28 :   GEN x=content(y);
    2110                 :         28 :   *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
    2111                 :            : }
    2112                 :            : 
    2113                 :            : static int
    2114                 :       3647 : must_negate(GEN x)
    2115                 :            : {
    2116                 :       3647 :   GEN t = leading_term(x);
    2117      [ +  -  + ]:       3647 :   switch(typ(t))
    2118                 :            :   {
    2119                 :            :     case t_INT: case t_REAL:
    2120                 :       1484 :       return (signe(t) < 0);
    2121                 :            :     case t_FRAC:
    2122                 :          0 :       return (signe(gel(t,1)) < 0);
    2123                 :            :   }
    2124                 :       3647 :   return 0;
    2125                 :            : }
    2126                 :            : 
    2127                 :            : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
    2128                 :            : GEN
    2129                 :        343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
    2130                 :            : {
    2131                 :            :   pari_sp av, av2, tetpil;
    2132                 :            :   long signh; /* junk */
    2133                 :        343 :   long dx, dy, vx, tx = typ(x), ty = typ(y);
    2134                 :            :   GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
    2135                 :            : 
    2136         [ -  + ]:        343 :   if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
    2137         [ -  + ]:        343 :   if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
    2138         [ -  + ]:        343 :   if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
    2139                 :        343 :   vx=varn(x);
    2140         [ +  + ]:        343 :   if (!signe(x))
    2141                 :            :   {
    2142         [ +  + ]:         14 :     if (signe(y)) return zero_extgcd(y,U,V,vx);
    2143                 :          7 :     *U = pol_0(vx); *V = pol_0(vx);
    2144                 :          7 :     return pol_0(vx);
    2145                 :            :   }
    2146         [ +  + ]:        329 :   if (!signe(y)) return zero_extgcd(x,V,U,vx);
    2147                 :        308 :   dx = degpol(x); dy = degpol(y);
    2148         [ +  + ]:        308 :   if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
    2149         [ +  + ]:        308 :   if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
    2150                 :            : 
    2151                 :        182 :   av = avma;
    2152                 :        182 :   u = x = primitive_part(x, &cu);
    2153                 :        182 :   v = y = primitive_part(y, &cv);
    2154                 :        182 :   g = h = gen_1; av2 = avma;
    2155                 :        182 :   um1 = gen_1; uze = gen_0;
    2156                 :            :   for(;;)
    2157                 :            :   {
    2158         [ +  + ]:        196 :     if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
    2159         [ -  + ]:         14 :     if (gc_needed(av2,1))
    2160                 :            :     {
    2161         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
    2162                 :          0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2163                 :            :     }
    2164                 :         14 :   }
    2165         [ +  + ]:        182 :   if (uze != gen_0) {
    2166                 :            :     GEN r;
    2167                 :         63 :     vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
    2168         [ -  + ]:         63 :     if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
    2169         [ -  + ]:         63 :     if (cu) uze = RgX_Rg_div(uze,cu);
    2170         [ -  + ]:         63 :     if (cv) vze = RgX_Rg_div(vze,cv);
    2171                 :         63 :     p1 = ginv(content(v));
    2172                 :            :   }
    2173                 :            :   else /* y | x */
    2174                 :            :   {
    2175         [ +  + ]:        119 :     vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
    2176                 :        119 :     uze = pol_0(vx);
    2177                 :        119 :     p1 = gen_1;
    2178                 :            :   }
    2179         [ +  + ]:        182 :   if (must_negate(v)) p1 = gneg(p1);
    2180                 :        182 :   tetpil = avma;
    2181                 :        182 :   z = RgX_Rg_mul(v,p1);
    2182                 :        182 :   *U = RgX_Rg_mul(uze,p1);
    2183                 :        182 :   *V = RgX_Rg_mul(vze,p1);
    2184                 :        182 :   gptr[0] = &z;
    2185                 :        182 :   gptr[1] = U;
    2186                 :        182 :   gptr[2] = V;
    2187                 :        343 :   gerepilemanysp(av,tetpil,gptr,3); return z;
    2188                 :            : }
    2189                 :            : 
    2190                 :            : int
    2191                 :         56 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
    2192                 :            : {
    2193                 :         56 :   pari_sp av = avma, av2, tetpil;
    2194                 :            :   long signh; /* junk */
    2195                 :            :   long vx;
    2196                 :            :   GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
    2197                 :            : 
    2198         [ -  + ]:         56 :   if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
    2199         [ -  + ]:         56 :   if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
    2200         [ -  + ]:         56 :   if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
    2201         [ -  + ]:         56 :   if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
    2202         [ -  + ]:         56 :   if (!signe(T)) {
    2203         [ #  # ]:          0 :     if (degpol(x) <= amax) {
    2204                 :          0 :       *P = RgX_copy(x);
    2205                 :          0 :       *Q = pol_1(varn(x));
    2206                 :          0 :       return 1;
    2207                 :            :     }
    2208                 :          0 :     return 0;
    2209                 :            :   }
    2210         [ -  + ]:         56 :   if (amax+bmax >= degpol(T))
    2211                 :          0 :     pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
    2212                 :            :                     mkvec3(stoi(amax), stoi(bmax), T));
    2213                 :         56 :   vx = varn(T);
    2214                 :         56 :   u = x = primitive_part(x, &cu);
    2215                 :         56 :   v = T = primitive_part(T, &cv);
    2216                 :         56 :   g = h = gen_1; av2 = avma;
    2217                 :         56 :   um1 = gen_1; uze = gen_0;
    2218                 :            :   for(;;)
    2219                 :            :   {
    2220                 :        112 :     (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
    2221 [ +  - ][ +  - ]:        112 :     if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
                 [ -  + ]
    2222 [ +  - ][ +  + ]:        112 :     if (typ(v)!=t_POL || degpol(v)<=amax) break;
    2223         [ -  + ]:         56 :     if (gc_needed(av2,1))
    2224                 :            :     {
    2225         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
    2226                 :          0 :       gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
    2227                 :            :     }
    2228                 :         56 :   }
    2229         [ -  + ]:         56 :   if (uze == gen_0)
    2230                 :            :   {
    2231                 :          0 :     avma = av; *P = pol_0(vx); *Q = pol_1(vx);
    2232                 :          0 :     return 1;
    2233                 :            :   }
    2234         [ +  + ]:         56 :   if (cu) uze = RgX_Rg_div(uze,cu);
    2235                 :         56 :   p1 = ginv(content(v));
    2236         [ -  + ]:         56 :   if (must_negate(v)) p1 = gneg(p1);
    2237                 :         56 :   tetpil = avma;
    2238                 :         56 :   *P = RgX_Rg_mul(v,p1);
    2239                 :         56 :   *Q = RgX_Rg_mul(uze,p1);
    2240                 :         56 :   gptr[0] = P;
    2241                 :         56 :   gptr[1] = Q;
    2242                 :         56 :   gerepilemanysp(av,tetpil,gptr,2); return 1;
    2243                 :            : }
    2244                 :            : 
    2245                 :            : /*******************************************************************/
    2246                 :            : /*                                                                 */
    2247                 :            : /*                RESULTANT USING DUCOS VARIANT                    */
    2248                 :            : /*                                                                 */
    2249                 :            : /*******************************************************************/
    2250                 :            : /* x^n / y^(n-1), assume n > 0 */
    2251                 :            : static GEN
    2252                 :      16653 : Lazard(GEN x, GEN y, long n)
    2253                 :            : {
    2254                 :            :   long a;
    2255                 :            :   GEN c;
    2256                 :            : 
    2257         [ +  + ]:      16653 :   if (n == 1) return x;
    2258                 :        794 :   a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
    2259                 :        794 :   c=x; n-=a;
    2260         [ +  + ]:       1756 :   while (a>1)
    2261                 :            :   {
    2262                 :        962 :     a>>=1; c=gdivexact(gsqr(c),y);
    2263         [ +  + ]:        962 :     if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
    2264                 :            :   }
    2265                 :      16653 :   return c;
    2266                 :            : }
    2267                 :            : 
    2268                 :            : /* F (x/y)^(n-1), assume n >= 1 */
    2269                 :            : static GEN
    2270                 :      20412 : Lazard2(GEN F, GEN x, GEN y, long n)
    2271                 :            : {
    2272         [ +  + ]:      20412 :   if (n == 1) return F;
    2273                 :      20412 :   return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
    2274                 :            : }
    2275                 :            : 
    2276                 :            : static GEN
    2277                 :      20412 : RgX_neg_i(GEN x, long lx)
    2278                 :            : {
    2279                 :            :   long i;
    2280                 :      20412 :   GEN y = cgetg(lx, t_POL); y[1] = x[1];
    2281         [ +  + ]:      58184 :   for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
    2282                 :      20412 :   return y;
    2283                 :            : }
    2284                 :            : static GEN
    2285                 :      59409 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
    2286                 :            : {
    2287                 :            :   long i;
    2288                 :            :   GEN z;
    2289         [ +  + ]:      59409 :   if (isrationalzero(x)) return pol_0(varn(y));
    2290                 :      59402 :   z = cgetg(ly,t_POL); z[1] = y[1];
    2291         [ +  + ]:     166670 :   for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
    2292                 :      59409 :   return z;
    2293                 :            : }
    2294                 :            : static long
    2295                 :      58541 : reductum_lg(GEN x, long lx)
    2296                 :            : {
    2297                 :      58541 :   long i = lx-2;
    2298 [ +  + ][ +  + ]:      64218 :   while (i > 1 && gequal0(gel(x,i))) i--;
    2299                 :      58541 :   return i+1;
    2300                 :            : }
    2301                 :            : 
    2302                 :            : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
    2303                 :            :  * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
    2304                 :            : static GEN
    2305                 :      20412 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
    2306                 :            : {
    2307                 :      20412 :   GEN p0, q0, h0, TMP, H, A, z0 = leading_term(Z);
    2308                 :            :   long p, q, j, lP, lQ;
    2309                 :            :   pari_sp av;
    2310                 :            : 
    2311                 :      20412 :   p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
    2312                 :      20412 :   q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
    2313                 :            :   /* p > q. Very often p - 1 = q */
    2314                 :      20412 :   av = avma;
    2315                 :            :   /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
    2316                 :      20412 :   H = RgX_neg_i(Z, lQ); /* deg H < q */
    2317                 :            : 
    2318         [ +  + ]:      20412 :   A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
    2319         [ +  + ]:      21826 :   for (j = q+1; j < p; j++)
    2320                 :            :   {
    2321         [ +  + ]:       1414 :     if (degpol(H) == q-1)
    2322                 :            :     { /* h0 = coeff of degree q-1 = leading coeff */
    2323                 :        994 :       h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
    2324                 :        994 :       H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
    2325                 :            :     }
    2326                 :            :     else
    2327                 :        420 :       H = RgX_shift_shallow(H, 1);
    2328         [ +  + ]:       1414 :     if (j+2 < lP)
    2329                 :            :     {
    2330                 :        903 :       TMP = RgX_Rg_mul(H, gel(P,j+2));
    2331         [ +  - ]:        903 :       A = A? RgX_add(A, TMP): TMP;
    2332                 :            :     }
    2333         [ +  + ]:       1414 :     if (gc_needed(av,1))
    2334                 :            :     {
    2335         [ -  + ]:        147 :       if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
    2336         [ +  - ]:        147 :       gerepileall(av,A?2:1,&H,&A);
    2337                 :            :     }
    2338                 :            :   }
    2339         [ +  + ]:      20412 :   if (q+2 < lP) lP = reductum_lg(P, q+3);
    2340                 :      20412 :   TMP = RgX_Rg_mul_i(P, z0, lP);
    2341         [ +  + ]:      20412 :   A = A? RgX_add(A, TMP): TMP;
    2342                 :      20412 :   A = RgX_Rg_divexact(A, p0);
    2343         [ +  + ]:      20412 :   if (degpol(H) == q-1)
    2344                 :            :   {
    2345                 :      20286 :     h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
    2346                 :      20286 :     A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
    2347                 :            :   }
    2348                 :            :   else
    2349                 :        126 :     A = RgX_Rg_mul(addshift(H,A), q0);
    2350                 :      20412 :   return RgX_Rg_divexact(A, s);
    2351                 :            : }
    2352                 :            : 
    2353                 :            : /* Ducos's subresultant */
    2354                 :            : GEN
    2355                 :      16912 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
    2356                 :            : {
    2357                 :            :   pari_sp av, av2;
    2358                 :      16912 :   long dP, dQ, delta, sig = 1;
    2359                 :            :   GEN cP, cQ, Z, s;
    2360                 :            : 
    2361                 :      16912 :   dP = degpol(P);
    2362                 :      16912 :   dQ = degpol(Q); delta = dP - dQ;
    2363         [ +  + ]:      16912 :   if (delta < 0)
    2364                 :            :   {
    2365         [ +  + ]:       1701 :     if (both_odd(dP, dQ)) sig = -1;
    2366                 :       1701 :     swap(P,Q); lswap(dP, dQ); delta = -delta;
    2367                 :            :   }
    2368         [ +  + ]:      16912 :   if (sol) *sol = gen_0;
    2369                 :      16912 :   av = avma;
    2370         [ +  + ]:      16912 :   if (dQ <= 0)
    2371                 :            :   {
    2372         [ -  + ]:        511 :     if (dQ < 0) return RgX_get_0(P);
    2373                 :        511 :     s = gpowgs(gel(Q,2), dP);
    2374         [ -  + ]:        511 :     if (sig == -1) s = gerepileupto(av, gneg(s));
    2375                 :        511 :     return s;
    2376                 :            :   }
    2377                 :      16401 :   P = primitive_part(P, &cP);
    2378                 :      16401 :   Q = primitive_part(Q, &cQ);
    2379                 :      16401 :   av2 = avma;
    2380                 :      16401 :   s = gpowgs(leading_term(Q),delta);
    2381         [ +  + ]:      16401 :   if (both_odd(dP, dQ)) sig = -sig;
    2382                 :      16401 :   Z = Q;
    2383                 :      16401 :   Q = RgX_pseudorem(P, Q);
    2384                 :      16401 :   P = Z;
    2385         [ +  + ]:      36813 :   while(degpol(Q) > 0)
    2386                 :            :   {
    2387                 :      20412 :     delta = degpol(P) - degpol(Q); /* > 0 */
    2388                 :      20412 :     Z = Lazard2(Q, leading_term(Q), s, delta);
    2389         [ +  + ]:      20412 :     if (both_odd(degpol(P), degpol(Q))) sig = -sig;
    2390                 :      20412 :     Q = nextSousResultant(P, Q, Z, s);
    2391                 :      20412 :     P = Z;
    2392         [ +  + ]:      20412 :     if (gc_needed(av,1))
    2393                 :            :     {
    2394         [ -  + ]:         13 :       if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
    2395                 :         13 :       gerepileall(av2,2,&P,&Q);
    2396                 :            :     }
    2397                 :      20412 :     s = leading_term(P);
    2398                 :            :   }
    2399         [ -  + ]:      16401 :   if (!signe(Q)) { avma = av; return RgX_get_0(Q); }
    2400                 :      16401 :   s = Lazard(leading_term(Q), s, degpol(P));
    2401         [ +  + ]:      16401 :   if (sig == -1) s = gneg(s);
    2402         [ +  + ]:      16401 :   if (cP) s = gmul(s, gpowgs(cP,dQ));
    2403         [ +  + ]:      16401 :   if (cQ) s = gmul(s, gpowgs(cQ,dP));
    2404         [ +  + ]:      16401 :   if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
    2405                 :      16912 :   return gerepilecopy(av, s);
    2406                 :            : }
    2407                 :            : /* Return resultant(P,Q). If sol != NULL: set *sol to the last non-constant
    2408                 :            :  * polynomial in the prs IF the sequence was computed, and gen_0 otherwise.
    2409                 :            :  * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
    2410                 :            :  * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
    2411                 :            :  * in the "generic" case. */
    2412                 :            : GEN
    2413                 :      39346 : resultant_all(GEN P, GEN Q, GEN *sol)
    2414                 :            : {
    2415                 :            :   long TP, TQ;
    2416                 :            :   GEN s;
    2417                 :            : 
    2418         [ +  + ]:      39346 :   if (sol) *sol = gen_0;
    2419         [ +  + ]:      39346 :   if ((s = init_resultant(P,Q))) return s;
    2420 [ +  - ][ -  + ]:      39269 :   if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
    2421                 :          0 :     return resultant2(P,Q); /* inexact */
    2422 [ +  + ][ +  + ]:      39269 :   if (TP && TQ) /* rational */
    2423                 :            :   {
    2424 [ +  - ][ +  + ]:      22588 :     if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
    2425                 :      10125 :     return QX_resultant(P,Q);
    2426                 :            :   }
    2427                 :      39346 :   return RgX_resultant_all(P, Q, sol);
    2428                 :            : }
    2429                 :            : 
    2430                 :            : /*******************************************************************/
    2431                 :            : /*                                                                 */
    2432                 :            : /*               RESULTANT USING SYLVESTER MATRIX                  */
    2433                 :            : /*                                                                 */
    2434                 :            : /*******************************************************************/
    2435                 :            : static GEN
    2436                 :          0 : _pol_0(void)
    2437                 :            : {
    2438                 :          0 :   GEN x = cgetg(3,t_POL);
    2439                 :          0 :   x[1] = 0;
    2440                 :          0 :   gel(x,2) = gen_0; return x;
    2441                 :            : }
    2442                 :            : 
    2443                 :            : static GEN
    2444                 :        133 : sylvester_col(GEN x, long j, long d, long D)
    2445                 :            : {
    2446                 :        133 :   GEN c = cgetg(d+1,t_COL);
    2447                 :            :   long i;
    2448         [ +  + ]:        238 :   for (i=1; i< j; i++) gel(c,i) = gen_0;
    2449         [ +  + ]:        588 :   for (   ; i<=D; i++) gel(c,i) = gel(x,D-i+2);
    2450         [ +  + ]:        238 :   for (   ; i<=d; i++) gel(c,i) = gen_0;
    2451                 :        133 :   return c;
    2452                 :            : }
    2453                 :            : 
    2454                 :            : GEN
    2455                 :         28 : sylvestermatrix_i(GEN x, GEN y)
    2456                 :            : {
    2457                 :            :   long j,d,dx,dy;
    2458                 :            :   GEN M;
    2459                 :            : 
    2460         [ -  + ]:         28 :   dx = degpol(x); if (dx < 0) { dx = 0; x = _pol_0(); }
    2461         [ -  + ]:         28 :   dy = degpol(y); if (dy < 0) { dy = 0; y = _pol_0(); }
    2462                 :         28 :   d = dx+dy; M = cgetg(d+1,t_MAT);
    2463         [ +  + ]:         84 :   for (j=1; j<=dy; j++) gel(M,j)    = sylvester_col(x,j,d,j+dx);
    2464         [ +  + ]:        105 :   for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
    2465                 :         28 :   return M;
    2466                 :            : }
    2467                 :            : 
    2468                 :            : GEN
    2469                 :          7 : sylvestermatrix(GEN x, GEN y)
    2470                 :            : {
    2471                 :            :   long i,j,lx;
    2472         [ -  + ]:          7 :   if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
    2473         [ -  + ]:          7 :   if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
    2474         [ -  + ]:          7 :   if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
    2475                 :          7 :   x = sylvestermatrix_i(x,y); lx = lg(x);
    2476         [ +  + ]:         28 :   for (i=1; i<lx; i++)
    2477         [ +  + ]:         84 :     for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
    2478                 :          7 :   return x;
    2479                 :            : }
    2480                 :            : 
    2481                 :            : GEN
    2482                 :         21 : resultant2(GEN x, GEN y)
    2483                 :            : {
    2484                 :            :   pari_sp av;
    2485                 :            :   GEN r;
    2486         [ -  + ]:         21 :   if ((r = init_resultant(x,y))) return r;
    2487                 :         21 :   av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
    2488                 :            : }
    2489                 :            : 
    2490                 :            : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
    2491                 :            :  * if vx <= v, return subst(x, v, pol_x(v0))
    2492                 :            :  * if vx >  v, return scalarpol(x, v0) */
    2493                 :            : static GEN
    2494                 :         84 : fix_pol(GEN x, long v, long v0)
    2495                 :            : {
    2496                 :            :   long vx;
    2497         [ -  + ]:         84 :   if (typ(x) != t_POL) return x;
    2498                 :         84 :   vx = varn(x);
    2499         [ +  + ]:         84 :   if (v == vx)
    2500                 :            :   {
    2501         [ +  + ]:         28 :     if (v) { x = leafcopy(x); setvarn(x, v0); }
    2502                 :         28 :     return x;
    2503                 :            :   }
    2504         [ +  - ]:         56 :   if (varncmp(v, vx) > 0)
    2505                 :            :   {
    2506                 :         56 :     x = gsubst(x,v,pol_x(v0));
    2507 [ +  - ][ +  + ]:         56 :     if (typ(x) == t_POL && varn(x) == v0) return x;
    2508                 :            :   }
    2509                 :         84 :   return scalarpol_shallow(x, v0);
    2510                 :            : }
    2511                 :            : 
    2512                 :            : /* resultant of x and y with respect to variable v, or with respect to their
    2513                 :            :  * main variable if v < 0. */
    2514                 :            : GEN
    2515                 :        238 : polresultant0(GEN x, GEN y, long v, long flag)
    2516                 :            : {
    2517                 :        238 :   long v0 = 0;
    2518                 :        238 :   pari_sp av = avma;
    2519                 :            : 
    2520         [ +  + ]:        238 :   if (v >= 0)
    2521                 :            :   {
    2522                 :         28 :     v0 = fetch_var_higher();
    2523                 :         28 :     x = fix_pol(x,v, v0);
    2524                 :         28 :     y = fix_pol(y,v, v0);
    2525                 :            :   }
    2526      [ +  +  - ]:        238 :   switch(flag)
    2527                 :            :   {
    2528                 :            :     case 2:
    2529                 :        231 :     case 0: x=resultant_all(x,y,NULL); break;
    2530                 :          7 :     case 1: x=resultant2(x,y); break;
    2531                 :          0 :     default: pari_err_FLAG("polresultant");
    2532                 :            :   }
    2533         [ +  + ]:        238 :   if (v >= 0) (void)delete_var();
    2534                 :        238 :   return gerepileupto(av,x);
    2535                 :            : }
    2536                 :            : GEN
    2537                 :        532 : polresultantext0(GEN x, GEN y, long v)
    2538                 :            : {
    2539                 :            :   GEN R, U, V;
    2540                 :        532 :   long v0 = 0;
    2541                 :        532 :   pari_sp av = avma;
    2542                 :            : 
    2543         [ +  + ]:        532 :   if (v >= 0)
    2544                 :            :   {
    2545                 :         14 :     v0 = fetch_var_higher();
    2546                 :         14 :     x = fix_pol(x,v, v0);
    2547                 :         14 :     y = fix_pol(y,v, v0);
    2548                 :            :   }
    2549                 :        532 :   R = subresext_i(x,y, &U,&V);
    2550         [ +  + ]:        532 :   if (v >= 0)
    2551                 :            :   {
    2552                 :         14 :     (void)delete_var();
    2553 [ +  - ][ +  - ]:         14 :     if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
    2554 [ +  - ][ +  - ]:         14 :     if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
    2555                 :            :   }
    2556                 :        532 :   return gerepilecopy(av, mkvec3(U,V,R));
    2557                 :            : }
    2558                 :            : GEN
    2559                 :        504 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
    2560                 :            : 
    2561                 :            : /*******************************************************************/
    2562                 :            : /*                                                                 */
    2563                 :            : /*             CHARACTERISTIC POLYNOMIAL USING RESULTANT           */
    2564                 :            : /*                                                                 */
    2565                 :            : /*******************************************************************/
    2566                 :            : 
    2567                 :            : /* (v - x)^d */
    2568                 :            : static GEN
    2569                 :         84 : caract_const(pari_sp av, GEN x, long v, long d)
    2570                 :         84 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
    2571                 :            : 
    2572                 :            : /* return caract(Mod(x,T)) in variable v */
    2573                 :            : GEN
    2574                 :      12565 : RgXQ_charpoly(GEN x, GEN T, long v)
    2575                 :            : {
    2576                 :      12565 :   pari_sp av = avma;
    2577                 :      12565 :   long d = degpol(T), dx, vx, vp, v0;
    2578                 :            :   GEN ch, L;
    2579                 :            : 
    2580         [ -  + ]:      12565 :   if (typ(x) != t_POL) return caract_const(av, x, v, d);
    2581                 :      12565 :   vx = varn(x);
    2582                 :      12565 :   vp = varn(T);
    2583         [ -  + ]:      12565 :   if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
    2584         [ -  + ]:      12565 :   if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
    2585                 :      12565 :   dx = degpol(x);
    2586         [ +  + ]:      12565 :   if (dx <= 0)
    2587         [ -  + ]:         14 :     return dx? monomial(gen_1, d, v): caract_const(av, gel(x,2), v, d);
    2588                 :            : 
    2589                 :      12551 :   v0 = fetch_var_higher();
    2590                 :      12551 :   x = RgX_neg(x);
    2591                 :      12551 :   gel(x,2) = gadd(gel(x,2), pol_x(v));
    2592                 :      12551 :   setvarn(x, v0);
    2593                 :      12551 :   T = leafcopy(T); setvarn(T, v0);
    2594                 :      12551 :   ch = resultant_all(T, x, NULL);
    2595                 :      12551 :   (void)delete_var();
    2596                 :            :   /* test for silly input: x mod (deg 0 polynomial) */
    2597         [ -  + ]:      12551 :   if (typ(ch) != t_POL) { avma = av; return pol_1(v); }
    2598                 :            : 
    2599                 :      12551 :   L = leading_term(ch);
    2600         [ -  + ]:      12551 :   if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
    2601                 :      12565 :   return gerepileupto(av, ch);
    2602                 :            : }
    2603                 :            : 
    2604                 :            : /* characteristic polynomial (in v) of x over nf, where x is an element of the
    2605                 :            :  * algebra nf[t]/(Q(t)) */
    2606                 :            : GEN
    2607                 :        224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
    2608                 :            : {
    2609                 :        224 :   const char *f = "rnfcharpoly";
    2610                 :        224 :   long dQ = degpol(Q);
    2611                 :        224 :   pari_sp av = avma;
    2612                 :            :   GEN T;
    2613                 :            : 
    2614         [ +  - ]:        224 :   if (v < 0) v = 0;
    2615                 :        224 :   nf = checknf(nf); T = nf_get_pol(nf);
    2616                 :        224 :   Q = RgX_nffix(f, T,Q,0);
    2617   [ +  +  +  + ]:        224 :   switch(typ(x))
    2618                 :            :   {
    2619                 :            :     case t_INT:
    2620                 :         28 :     case t_FRAC: return caract_const(av, x, v, dQ);
    2621                 :            :     case t_POLMOD:
    2622                 :         91 :       x = polmod_nffix2(f,T,Q, x,0);
    2623                 :         49 :       break;
    2624                 :            :     case t_POL:
    2625         [ +  + ]:         56 :       x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
    2626                 :         42 :       break;
    2627                 :         49 :     default: pari_err_TYPE(f,x);
    2628                 :            :   }
    2629         [ +  + ]:         91 :   if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
    2630                 :            :   /* x a t_POL in variable vQ */
    2631         [ +  + ]:         49 :   if (degpol(x) >= dQ) x = RgX_rem(x, Q);
    2632         [ -  + ]:         49 :   if (dQ <= 1) return caract_const(av, constant_term(x), v, 1);
    2633                 :        119 :   return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
    2634                 :            : }
    2635                 :            : 
    2636                 :            : /*******************************************************************/
    2637                 :            : /*                                                                 */
    2638                 :            : /*                  GCD USING SUBRESULTANT                         */
    2639                 :            : /*                                                                 */
    2640                 :            : /*******************************************************************/
    2641                 :            : static int inexact(GEN x, int *simple, int *rational);
    2642                 :            : static int
    2643                 :    6749577 : isinexactall(GEN x, int *simple, int *rational)
    2644                 :            : {
    2645                 :    6749577 :   long i, lx = lg(x);
    2646         [ +  + ]:   31761052 :   for (i=2; i<lx; i++)
    2647         [ -  + ]:   25011475 :     if (inexact(gel(x,i), simple, rational)) return 1;
    2648                 :    6749577 :   return 0;
    2649                 :            : }
    2650                 :            : /* return 1 if coeff explosion is not possible */
    2651                 :            : static int
    2652                 :   25011531 : inexact(GEN x, int *simple, int *rational)
    2653                 :            : {
    2654                 :   25011531 :   int junk = 0;
    2655   [ +  -  +  -  :   25011531 :   switch(typ(x))
             -  +  +  +  
                      - ]
    2656                 :            :   {
    2657                 :   24972149 :     case t_INT: case t_FRAC: return 0;
    2658                 :            : 
    2659                 :          0 :     case t_REAL: case t_PADIC: case t_SER: return 1;
    2660                 :            : 
    2661                 :            :     case t_INTMOD:
    2662                 :            :     case t_FFELT:
    2663                 :      12271 :       *rational = 0;
    2664         [ +  + ]:      12271 :       if (!*simple) *simple = 1;
    2665                 :      12271 :       return 0;
    2666                 :            : 
    2667                 :            :     case t_COMPLEX:
    2668                 :          0 :       *rational = 0;
    2669                 :          0 :       return inexact(gel(x,1), simple, rational)
    2670 [ #  # ][ #  # ]:          0 :           || inexact(gel(x,2), simple, rational);
    2671                 :            :     case t_QUAD:
    2672                 :          0 :       *rational = *simple = 0;
    2673                 :          0 :       return inexact(gel(x,2), &junk, rational)
    2674 [ #  # ][ #  # ]:          0 :           || inexact(gel(x,3), &junk, rational);
    2675                 :            : 
    2676                 :            :     case t_POLMOD:
    2677                 :       3787 :       *rational = 0;
    2678                 :       3787 :       return isinexactall(gel(x,1), simple, rational);
    2679                 :            :     case t_POL:
    2680                 :      23296 :       *rational = 0;
    2681                 :      23296 :       *simple = -1;
    2682                 :      23296 :       return isinexactall(x, &junk, rational);
    2683                 :            :     case t_RFRAC:
    2684                 :         28 :       *rational = 0;
    2685                 :         28 :       *simple = -1;
    2686                 :         56 :       return inexact(gel(x,1), &junk, rational)
    2687 [ +  - ][ -  + ]:         28 :           || inexact(gel(x,2), &junk, rational);
    2688                 :            :   }
    2689                 :          0 :   *rational = 0;
    2690                 :   25011531 :   *simple = -1; return 0;
    2691                 :            : }
    2692                 :            : 
    2693                 :            : /* x monomial, y t_POL in the same variable */
    2694                 :            : static GEN
    2695                 :    7183316 : gcdmonome(GEN x, GEN y)
    2696                 :            : {
    2697                 :    7183316 :   pari_sp av = avma;
    2698                 :    7183316 :   long dx = degpol(x), e = RgX_valrem(y, &y);
    2699                 :    7183316 :   long i, l = lg(y);
    2700                 :    7183316 :   GEN t, v = cgetg(l, t_VEC);
    2701                 :    7183316 :   gel(v,1) = gel(x,dx+2);
    2702         [ +  + ]:   15977253 :   for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
    2703                 :    7183316 :   t = content(v); /* gcd(lc(x), cont(y)) */
    2704                 :    7183316 :   t = simplify_shallow(t);
    2705         [ +  + ]:    7183316 :   if (dx < e) e = dx;
    2706                 :    7183316 :   return gerepileupto(av, monomialcopy(t, e, varn(x)));
    2707                 :            : }
    2708                 :            : 
    2709                 :            : /* x, y are t_POL in the same variable */
    2710                 :            : GEN
    2711                 :   10544864 : RgX_gcd(GEN x, GEN y)
    2712                 :            : {
    2713                 :            :   long dx, dy;
    2714                 :            :   pari_sp av, av1;
    2715                 :            :   GEN d, g, h, p1, p2, u, v;
    2716                 :   10544864 :   int simple = 0, rational = 1;
    2717                 :            : 
    2718         [ +  + ]:   10544864 :   if (isexactzero(y)) return RgX_copy(x);
    2719         [ +  + ]:   10544584 :   if (isexactzero(x)) return RgX_copy(y);
    2720         [ +  + ]:   10544563 :   if (RgX_is_monomial(x)) return gcdmonome(x,y);
    2721         [ +  + ]:    4016142 :   if (RgX_is_monomial(y)) return gcdmonome(y,x);
    2722 [ +  - ][ -  + ]:    3361247 :   if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
    2723                 :            :   {
    2724                 :          0 :     av = avma; u = ggcd(content(x), content(y));
    2725                 :          0 :     return gerepileupto(av, scalarpol(u, varn(x)));
    2726                 :            :   }
    2727         [ +  + ]:    3361247 :   if (rational) return QX_gcd(x,y); /* Q[X] */
    2728                 :            : 
    2729                 :       4585 :   av = avma;
    2730         [ +  + ]:       4585 :   if (simple > 0) x = RgX_gcd_simple(x,y);
    2731                 :            :   else
    2732                 :            :   {
    2733                 :       3458 :     dx = lg(x); dy = lg(y);
    2734         [ +  + ]:       3458 :     if (dx < dy) { swap(x,y); lswap(dx,dy); }
    2735         [ -  + ]:       3458 :     if (dy==3)
    2736                 :            :     {
    2737                 :          0 :       d = ggcd(gel(y,2), content(x));
    2738                 :          0 :       return gerepileupto(av, scalarpol(d, varn(x)));
    2739                 :            :     }
    2740         [ +  + ]:       3458 :     u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
    2741         [ +  + ]:       3458 :     v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
    2742                 :       3458 :     d = ggcd(p1,p2);
    2743                 :       3458 :     av1 = avma;
    2744                 :       3458 :     g = h = gen_1;
    2745                 :            :     for(;;)
    2746                 :            :     {
    2747                 :       4879 :       GEN r = RgX_pseudorem(u,v);
    2748                 :       4879 :       long degq, du, dv, dr = lg(r);
    2749                 :            : 
    2750         [ +  + ]:       4879 :       if (!signe(r)) break;
    2751         [ +  + ]:       2597 :       if (dr <= 3)
    2752                 :            :       {
    2753                 :       1176 :         avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
    2754                 :            :       }
    2755         [ -  + ]:       1421 :       if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
    2756                 :       1421 :       du = lg(u); dv = lg(v); degq = du-dv;
    2757                 :       1421 :       u = v; p1 = g; g = leading_term(u);
    2758      [ +  +  + ]:       1421 :       switch(degq)
    2759                 :            :       {
    2760                 :         56 :         case 0: break;
    2761                 :            :         case 1:
    2762                 :       1204 :           p1 = gmul(h,p1); h = g; break;
    2763                 :            :         default:
    2764                 :        161 :           p1 = gmul(gpowgs(h,degq), p1);
    2765                 :        161 :           h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
    2766                 :            :       }
    2767                 :       1421 :       v = RgX_Rg_div(r,p1);
    2768         [ -  + ]:       1421 :       if (gc_needed(av1,1))
    2769                 :            :       {
    2770         [ #  # ]:          0 :         if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
    2771                 :          0 :         gerepileall(av1,4, &u,&v,&g,&h);
    2772                 :            :       }
    2773                 :       1421 :     }
    2774                 :       2282 :     x = RgX_Rg_mul(primpart(v), d);
    2775                 :            :   }
    2776         [ +  + ]:       3409 :   if (must_negate(x)) x = RgX_neg(x);
    2777                 :   10544864 :   return gerepileupto(av,x);
    2778                 :            : }
    2779                 :            : 
    2780                 :            : static GEN
    2781                 :       3584 : RgX_disc_aux(GEN x)
    2782                 :            : {
    2783                 :       3584 :   long dx = degpol(x), Tx;
    2784                 :            :   GEN D, L, y, p;
    2785 [ +  - ][ -  + ]:       3584 :   if (!signe(x) || !dx) return RgX_get_0(x);
    2786         [ +  + ]:       3584 :   if (dx == 1) return RgX_get_1(x);
    2787         [ +  + ]:       3556 :   if (dx == 2) {
    2788                 :        280 :     GEN a = gel(x,4), b = gel(x,3), c = gel(x,2);
    2789                 :        280 :     return gsub(gsqr(b), gmul2n(gmul(a,c),2));
    2790                 :            :   }
    2791                 :       3276 :   Tx = RgX_simpletype(x);
    2792         [ +  + ]:       3276 :   if (Tx == t_INT) return ZX_disc(x);
    2793         [ -  + ]:        161 :   if (Tx == t_FRAC) return QX_disc(x);
    2794                 :        161 :   p = NULL;
    2795 [ +  + ][ +  - ]:        161 :   if (RgX_is_FpX(x, &p) && p)
    2796                 :         14 :     return Fp_to_mod(FpX_disc(RgX_to_FpX(x,p), p), p);
    2797                 :            : 
    2798                 :        147 :   y = RgX_deriv(x);
    2799         [ -  + ]:        147 :   if (!signe(y)) return RgX_get_0(y);
    2800         [ +  + ]:        147 :   if (Tx == t_REAL)
    2801                 :         14 :     D = resultant2(x,y);
    2802                 :            :   else
    2803                 :            :   {
    2804                 :        133 :     D = RgX_resultant_all(x, y, NULL);
    2805         [ -  + ]:        133 :     if (D == gen_0) return RgX_get_0(y);
    2806                 :            :   }
    2807         [ +  + ]:        147 :   L = leading_term(x); if (!gequal1(L)) D = gdiv(D,L);
    2808         [ +  + ]:        147 :   if (dx & 2) D = gneg(D);
    2809                 :       3584 :   return D;
    2810                 :            : }
    2811                 :            : GEN
    2812                 :        567 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
    2813                 :            : 
    2814                 :            : GEN
    2815                 :       3024 : poldisc0(GEN x, long v)
    2816                 :            : {
    2817                 :            :   pari_sp av;
    2818   [ +  -  -  -  :       3024 :   switch(typ(x))
                +  -  - ]
    2819                 :            :   {
    2820                 :            :     case t_POL:
    2821                 :            :     {
    2822                 :            :       GEN D;
    2823                 :       3017 :       long v0 = -1;
    2824                 :       3017 :       av = avma;
    2825 [ -  + ][ #  # ]:       3017 :       if (v >= 0 && v != varn(x))
    2826                 :            :       {
    2827                 :          0 :         v0 = fetch_var_higher();
    2828                 :          0 :         x = fix_pol(x,v, v0);
    2829                 :            :       }
    2830                 :       3017 :       D = RgX_disc_aux(x);
    2831         [ -  + ]:       3017 :       if (v0 >= 0) (void)delete_var();
    2832                 :       3017 :       return gerepileupto(av, D);
    2833                 :            :     }
    2834                 :            : 
    2835                 :            :     case t_COMPLEX:
    2836                 :          0 :       return utoineg(4);
    2837                 :            : 
    2838                 :            :     case t_QUAD:
    2839                 :          0 :       return quad_disc(x);
    2840                 :            :     case t_POLMOD:
    2841                 :          0 :       return poldisc0(gel(x,1), v);
    2842                 :            : 
    2843                 :            :     case t_QFR: case t_QFI:
    2844                 :          7 :       av = avma; return gerepileuptoint(av, qfb_disc(x));
    2845                 :            : 
    2846                 :            :     case t_VEC: case t_COL: case t_MAT:
    2847                 :            :     {
    2848                 :            :       long i;
    2849                 :          0 :       GEN z = cgetg_copy(x, &i);
    2850         [ #  # ]:          0 :       for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
    2851                 :          0 :       return z;
    2852                 :            :     }
    2853                 :            :   }
    2854                 :          0 :   pari_err_TYPE("poldisc",x);
    2855                 :       3024 :   return NULL; /* not reached */
    2856                 :            : }
    2857                 :            : 
    2858                 :            : GEN
    2859                 :          7 : reduceddiscsmith(GEN x)
    2860                 :            : {
    2861                 :          7 :   long j, n = degpol(x);
    2862                 :          7 :   pari_sp av = avma;
    2863                 :            :   GEN xp, M;
    2864                 :            : 
    2865         [ -  + ]:          7 :   if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
    2866         [ -  + ]:          7 :   if (n<=0) pari_err_CONSTPOL("poldiscreduced");
    2867                 :          7 :   RgX_check_ZX(x,"poldiscreduced");
    2868         [ -  + ]:          7 :   if (!gequal1(gel(x,n+2)))
    2869                 :          0 :     pari_err_IMPL("non-monic polynomial in poldiscreduced");
    2870                 :          7 :   M = cgetg(n+1,t_MAT);
    2871                 :          7 :   xp = ZX_deriv(x);
    2872         [ +  + ]:         28 :   for (j=1; j<=n; j++)
    2873                 :            :   {
    2874                 :         21 :     gel(M,j) = RgX_to_RgC(xp, n);
    2875         [ +  + ]:         21 :     if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
    2876                 :            :   }
    2877                 :          7 :   return gerepileupto(av, ZM_snf(M));
    2878                 :            : }
    2879                 :            : 
    2880                 :            : /***********************************************************************/
    2881                 :            : /**                                                                   **/
    2882                 :            : /**                       STURM ALGORITHM                             **/
    2883                 :            : /**              (number of real roots of x in [a,b])                 **/
    2884                 :            : /**                                                                   **/
    2885                 :            : /***********************************************************************/
    2886                 :            : static int
    2887                 :        651 : exact_sturm(GEN a)
    2888                 :            : {
    2889         [ +  + ]:        651 :   switch(typ(a))
    2890                 :            :   {
    2891                 :        644 :     case t_INT: case t_FRAC: case t_INFINITY: return 1;
    2892                 :        651 :     default: return 0;
    2893                 :            :   }
    2894                 :            : }
    2895                 :            : 
    2896                 :            : /* Deprecated: support the old format: if a (resp. b) is NULL, set it
    2897                 :            :  * to -oo resp. +oo). ZX_sturmpart() should be preferred  */
    2898                 :            : static long
    2899                 :        350 : sturmpart_i(GEN x, GEN a, GEN b)
    2900                 :            : {
    2901                 :            :   long sl, sr, s, t, r1;
    2902                 :        350 :   pari_sp av = avma;
    2903                 :            :   int integral;
    2904                 :            :   GEN g,h,u,v;
    2905                 :            : 
    2906         [ -  + ]:        350 :   if (gequal0(x)) pari_err_ROOTS0("sturm");
    2907                 :        350 :   t = typ(x);
    2908         [ -  + ]:        350 :   if (t != t_POL)
    2909                 :            :   {
    2910 [ #  # ][ #  # ]:          0 :     if (t == t_INT || t == t_REAL || t == t_FRAC) return 0;
                 [ #  # ]
    2911                 :          0 :     pari_err_TYPE("sturm",x);
    2912                 :            :   }
    2913         [ -  + ]:        350 :   s=lg(x); if (s==3) return 0;
    2914                 :        350 :   u = primpart(x);
    2915                 :        350 :   integral = RgX_is_ZX(u);
    2916 [ +  + ][ +  + ]:        350 :   if (!b && a && typ(a) == t_VEC && lg(a) == 3)
         [ +  - ][ +  - ]
    2917                 :            :   { /* new format */
    2918 [ +  - ][ +  - ]:         70 :     if (integral && exact_sturm(gel(a,1)) && exact_sturm(gel(a,2)))
                 [ +  - ]
    2919                 :            :     {
    2920         [ -  + ]:         70 :       if (!ZX_is_squarefree(u))
    2921                 :          0 :         pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,u);
    2922                 :         70 :       return ZX_sturmpart(u, a);
    2923                 :            :     }
    2924                 :            :     /* but can't use new function; convert to old form */
    2925                 :          0 :     integral = 0;
    2926                 :          0 :     b = gel(a,2);
    2927         [ #  # ]:          0 :     if (typ(b) == t_INFINITY)
    2928                 :            :     {
    2929         [ #  # ]:          0 :       if (inf_get_sign(b) < 0) return 0;
    2930                 :          0 :       b = NULL;
    2931                 :            :     }
    2932                 :          0 :     a = gel(a,1);
    2933         [ #  # ]:          0 :     if (typ(a) == t_INFINITY)
    2934                 :            :     {
    2935         [ #  # ]:          0 :       if (inf_get_sign(a) > 0) return 0;
    2936                 :          0 :       a = NULL;
    2937                 :            :     }
    2938                 :            :   }
    2939         [ +  + ]:        280 :   if (integral)
    2940                 :            :   {
    2941         [ +  + ]:        259 :     if (!a) a = mkmoo();
    2942         [ +  + ]:        259 :     if (!b) b = mkoo();
    2943 [ +  + ][ +  - ]:        259 :     if (exact_sturm(a) && exact_sturm(b))
    2944                 :            :     {
    2945         [ +  + ]:        252 :       if (!ZX_is_squarefree(u))
    2946                 :          7 :         pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,u);
    2947                 :        245 :       return ZX_sturmpart(u, mkvec2(a,b));
    2948                 :            :     }
    2949                 :            :   }
    2950                 :            :   /* legacy code: should only be used if we have a t_REAL somewhere; and even
    2951                 :            :    * then, the calling program should be changed */
    2952                 :         28 :   sl = gsigne(leading_term(u));
    2953 [ +  + ][ -  + ]:         28 :   t = a? gsigne(poleval(u,a)): (odd(s)? sl: -sl);
    2954         [ -  + ]:         28 :   if (s==4)
    2955                 :            :   {
    2956         [ #  # ]:          0 :     if (t == 0) return 1;
    2957         [ #  # ]:          0 :     s = b? gsigne(poleval(u,b)):  sl;
    2958                 :          0 :     return (s == t)? 0: 1;
    2959                 :            :   }
    2960         [ +  + ]:         28 :   s = b? gsigne(poleval(u,b)): sl;
    2961                 :         28 :   r1= (t == 0)? 1: 0;
    2962                 :         28 :   v = primpart(RgX_deriv(x));
    2963         [ +  + ]:         28 :   sr = b? gsigne(poleval(v,b)): s;
    2964         [ +  - ]:         28 :   if (sr)
    2965                 :            :   {
    2966         [ +  + ]:         28 :     if (!s) s=sr;
    2967         [ -  + ]:         21 :     else if (sr!=s) { s= -s; r1--; }
    2968                 :            :   }
    2969         [ +  + ]:         28 :   sr = a? gsigne(poleval(v,a)): -t;
    2970         [ +  - ]:         28 :   if (sr)
    2971                 :            :   {
    2972         [ +  + ]:         28 :     if (!t) t=sr;
    2973         [ +  + ]:         21 :     else if (sr!=t) { t= -t; r1++; }
    2974                 :            :   }
    2975                 :         28 :   g=gen_1; h=gen_1;
    2976                 :            :   for(;;)
    2977                 :            :   {
    2978                 :        147 :     GEN p1, r = RgX_pseudorem(u,v);
    2979                 :        147 :     long du=lg(u), dv=lg(v), dr=lg(r), degq=du-dv;
    2980                 :            : 
    2981         [ +  + ]:        147 :     if (dr<=2) pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
    2982 [ +  + ][ +  - ]:        140 :     if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
    2983                 :        140 :     sl = gsigne(gel(r,dr-1));
    2984         [ +  + ]:        140 :     sr = b? gsigne(poleval(r,b)): sl;
    2985         [ +  - ]:        140 :     if (sr)
    2986                 :            :     {
    2987         [ -  + ]:        140 :       if (!s) s=sr;
    2988         [ +  + ]:        140 :       else if (sr!=s) { s= -s; r1--; }
    2989                 :            :     }
    2990 [ +  + ][ +  + ]:        140 :     sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
    2991         [ +  + ]:        140 :     if (sr)
    2992                 :            :     {
    2993         [ -  + ]:        133 :       if (!t) t=sr;
    2994         [ +  + ]:        133 :       else if (sr!=t) { t= -t; r1++; }
    2995                 :            :     }
    2996         [ +  + ]:        140 :     if (dr==3) return r1;
    2997                 :            : 
    2998                 :        119 :     u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
    2999      [ -  +  - ]:        119 :     switch(degq)
    3000                 :            :     {
    3001                 :          0 :       case 0: break;
    3002                 :            :       case 1:
    3003                 :        119 :         p1 = gmul(h,p1); h = g; break;
    3004                 :            :       default:
    3005                 :          0 :         p1 = gmul(gpowgs(h,degq),p1);
    3006                 :          0 :         h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
    3007                 :            :     }
    3008                 :        119 :     v = RgX_Rg_divexact(r,p1);
    3009         [ -  + ]:        119 :     if (gc_needed(av,1))
    3010                 :            :     {
    3011         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polsturm, dr = %ld",dr);
    3012                 :          0 :       gerepileall(av,4,&u,&v,&g,&h);
    3013                 :            :     }
    3014                 :        455 :   }
    3015                 :            : }
    3016                 :            : long
    3017                 :        350 : sturmpart(GEN x, GEN a, GEN b)
    3018                 :            : {
    3019                 :        350 :   pari_sp av = avma;
    3020                 :        350 :   long r = sturmpart_i(x,a,b);
    3021                 :        336 :   avma = av; return r;
    3022                 :            : }
    3023                 :            : long
    3024                 :          0 : RgX_sturmpart(GEN x, GEN ab) { return sturmpart(x, ab, NULL); }
    3025                 :            : 
    3026                 :            : /***********************************************************************/
    3027                 :            : /**                                                                   **/
    3028                 :            : /**                        GENERIC EXTENDED GCD                       **/
    3029                 :            : /**                                                                   **/
    3030                 :            : /***********************************************************************/
    3031                 :            : /* assume typ(x) = typ(y) = t_POL */
    3032                 :            : GEN
    3033                 :     119092 : RgXQ_inv(GEN x, GEN y)
    3034                 :            : {
    3035                 :     119092 :   long vx=varn(x), vy=varn(y);
    3036                 :            :   pari_sp av;
    3037                 :            :   GEN u, v, d;
    3038                 :            : 
    3039         [ -  + ]:     119092 :   while (vx != vy)
    3040                 :            :   {
    3041         [ #  # ]:          0 :     if (varncmp(vx,vy) > 0)
    3042                 :            :     {
    3043         [ #  # ]:          0 :       d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
    3044                 :          0 :       return scalarpol(d, vy);
    3045                 :            :     }
    3046         [ #  # ]:          0 :     if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3047                 :          0 :     x = gel(x,2); vx = gvar(x);
    3048                 :            :   }
    3049                 :     119092 :   av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
    3050         [ +  + ]:     119092 :   if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
    3051                 :     119085 :   d = gdiv(u,d);
    3052 [ +  + ][ -  + ]:     119085 :   if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
    3053                 :     119085 :   return gerepileupto(av, d);
    3054                 :            : }
    3055                 :            : 
    3056                 :            : /*Assume x is a polynomial and y is not */
    3057                 :            : static GEN
    3058                 :        112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
    3059                 :            : {
    3060                 :        112 :   long vx = varn(x);
    3061                 :        112 :   int xis0 = signe(x)==0, yis0 = gequal0(y);
    3062 [ +  + ][ +  + ]:        112 :   if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
    3063         [ +  + ]:         84 :   if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
    3064                 :        112 :   *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
    3065                 :            : }
    3066                 :            : /* Assume x==0, y!=0 */
    3067                 :            : static GEN
    3068                 :         63 : zero_bezout(GEN y, GEN *U, GEN *V)
    3069                 :            : {
    3070                 :         63 :   *U=gen_0; *V = ginv(y); return gen_1;
    3071                 :            : }
    3072                 :            : 
    3073                 :            : GEN
    3074                 :        273 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
    3075                 :            : {
    3076                 :        273 :   long tx=typ(x), ty=typ(y), vx;
    3077 [ +  + ][ +  + ]:        273 :   if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
    3078         [ +  + ]:        238 :   if (tx != t_POL)
    3079                 :            :   {
    3080         [ +  + ]:        140 :     if (ty == t_POL)
    3081                 :         56 :       return scalar_bezout(y,x,v,u);
    3082                 :            :     else
    3083                 :            :     {
    3084                 :         84 :       int xis0 = gequal0(x), yis0 = gequal0(y);
    3085 [ +  + ][ +  + ]:         84 :       if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
    3086         [ +  + ]:         63 :       if (xis0) return zero_bezout(y,u,v);
    3087                 :         42 :       else return zero_bezout(x,v,u);
    3088                 :            :     }
    3089                 :            :   }
    3090         [ +  + ]:         98 :   else if (ty != t_POL) return scalar_bezout(x,y,u,v);
    3091                 :         42 :   vx = varn(x);
    3092         [ -  + ]:         42 :   if (vx != varn(y))
    3093                 :          0 :     return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
    3094         [ #  # ]:          0 :                                    : scalar_bezout(y,x,v,u);
    3095                 :        273 :   return RgX_extgcd(x,y,u,v);
    3096                 :            : }
    3097                 :            : 
    3098                 :            : GEN
    3099                 :        273 : gcdext0(GEN x, GEN y)
    3100                 :            : {
    3101                 :        273 :   GEN z=cgetg(4,t_VEC);
    3102                 :        273 :   gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
    3103                 :        273 :   return z;
    3104                 :            : }
    3105                 :            : 
    3106                 :            : /*******************************************************************/
    3107                 :            : /*                                                                 */
    3108                 :            : /*                    GENERIC (modular) INVERSE                    */
    3109                 :            : /*                                                                 */
    3110                 :            : /*******************************************************************/
    3111                 :            : 
    3112                 :            : GEN
    3113                 :       4459 : ginvmod(GEN x, GEN y)
    3114                 :            : {
    3115                 :       4459 :   long tx=typ(x);
    3116                 :            : 
    3117      [ +  -  - ]:       4459 :   switch(typ(y))
    3118                 :            :   {
    3119                 :            :     case t_POL:
    3120         [ +  + ]:       4459 :       if (tx==t_POL) return RgXQ_inv(x,y);
    3121         [ +  - ]:       3759 :       if (is_scalar_t(tx)) return ginv(x);
    3122                 :          0 :       break;
    3123                 :            :     case t_INT:
    3124         [ #  # ]:          0 :       if (tx==t_INT) return Fp_inv(x,y);
    3125         [ #  # ]:          0 :       if (tx==t_POL) return gen_0;
    3126                 :            :   }
    3127                 :          0 :   pari_err_TYPE2("ginvmod",x,y);
    3128                 :       4452 :   return NULL; /* not reached */
    3129                 :            : }
    3130                 :            : 
    3131                 :            : /***********************************************************************/
    3132                 :            : /**                                                                   **/
    3133                 :            : /**                         NEWTON POLYGON                            **/
    3134                 :            : /**                                                                   **/
    3135                 :            : /***********************************************************************/
    3136                 :            : 
    3137                 :            : /* assume leading coeff of x is non-zero */
    3138                 :            : GEN
    3139                 :          7 : newtonpoly(GEN x, GEN p)
    3140                 :            : {
    3141                 :            :   GEN y;
    3142                 :            :   long n,ind,a,b,c,u1,u2,r1,r2;
    3143                 :          7 :   long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
    3144                 :            : 
    3145         [ -  + ]:          7 :   if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
    3146         [ -  + ]:          7 :   n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
    3147                 :          7 :   y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
    3148                 :          7 :   vval = (long *) pari_malloc(sizeof(long)*(n+1));
    3149         [ +  + ]:         42 :   for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
    3150         [ +  - ]:          7 :   for (a=0, ind=1; a<n; a++)
    3151                 :            :   {
    3152         [ +  - ]:          7 :     if (vval[a] != LONG_MAX) break;
    3153                 :          0 :     gel(y,ind++) = utoipos(LONG_MAX);
    3154                 :            :   }
    3155         [ +  + ]:         21 :   for (b=a+1; b<=n; a=b, b=a+1)
    3156                 :            :   {
    3157         [ -  + ]:         14 :     while (vval[b] == LONG_MAX) b++;
    3158                 :         14 :     u1=vval[a]-vval[b]; u2=b-a;
    3159         [ +  + ]:         49 :     for (c=b+1; c<=n; c++)
    3160                 :            :     {
    3161         [ -  + ]:         35 :       if (vval[c] == LONG_MAX) continue;
    3162                 :         35 :       r1=vval[a]-vval[c]; r2=c-a;
    3163         [ +  + ]:         35 :       if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
    3164                 :            :     }
    3165         [ +  + ]:         42 :     while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
    3166                 :            :   }
    3167                 :          7 :   pari_free(vval); return y;
    3168                 :            : }

Generated by: LCOV version 1.9