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

Generated by: LCOV version 1.9