Code coverage tests

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

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

The target is to exceed 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 - bibli2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 1159 1219 95.1 %
Date: 2020-09-18 06:10:04 Functions: 108 114 94.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*******************************************************************/
      18             : /**                                                               **/
      19             : /**                      SPECIAL POLYNOMIALS                      **/
      20             : /**                                                               **/
      21             : /*******************************************************************/
      22             : /* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
      23             :  * T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
      24             :  *   where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
      25             :  *   and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
      26             : GEN
      27        2156 : polchebyshev1(long n, long v) /* Assume 4*n < LONG_MAX */
      28             : {
      29             :   long k, l;
      30             :   pari_sp av;
      31             :   GEN q,a,r;
      32             : 
      33        2156 :   if (v<0) v = 0;
      34             :   /* polchebyshev(-n,1) = polchebyshev(n,1) */
      35        2156 :   if (n < 0) n = -n;
      36        2156 :   if (n==0) return pol_1(v);
      37        2135 :   if (n==1) return pol_x(v);
      38             : 
      39        2093 :   q = cgetg(n+3, t_POL); r = q + n+2;
      40        2093 :   a = int2n(n-1);
      41        2093 :   gel(r--,0) = a;
      42        2093 :   gel(r--,0) = gen_0;
      43       31955 :   for (k=1,l=n; l>1; k++,l-=2)
      44             :   {
      45       29862 :     av = avma;
      46       29862 :     a = diviuuexact(muluui(l, l-1, a), 4*k, n-k);
      47       29862 :     togglesign(a); a = gerepileuptoint(av, a);
      48       29862 :     gel(r--,0) = a;
      49       29862 :     gel(r--,0) = gen_0;
      50             :   }
      51        2093 :   q[1] = evalsigne(1) | evalvarn(v);
      52        2093 :   return q;
      53             : }
      54             : static void
      55          70 : polchebyshev1_eval_aux(long n, GEN x, GEN *pt1, GEN *pt2)
      56             : {
      57             :   GEN t1, t2, b;
      58          70 :   if (n == 1) { *pt1 = gen_1; *pt2 = x; return; }
      59          56 :   if (n == 0) { *pt1 = x; *pt2 = gen_1; return; }
      60          56 :   polchebyshev1_eval_aux((n+1) >> 1, x, &t1, &t2);
      61          56 :   b = gsub(gmul(gmul2n(t1,1), t2), x);
      62          56 :   if (odd(n)) { *pt1 = gadd(gmul2n(gsqr(t1), 1), gen_m1); *pt2 = b; }
      63          42 :   else        { *pt1 = b; *pt2 = gadd(gmul2n(gsqr(t2), 1), gen_m1); }
      64             : }
      65             : static GEN
      66          14 : polchebyshev1_eval(long n, GEN x)
      67             : {
      68             :   GEN t1, t2;
      69             :   long i, v;
      70             :   pari_sp av;
      71             : 
      72          14 :   if (n < 0) n = -n;
      73          14 :   if (n==0) return gen_1;
      74          14 :   if (n==1) return gcopy(x);
      75          14 :   av = avma;
      76          14 :   v = u_lvalrem(n, 2, (ulong*)&n);
      77          14 :   polchebyshev1_eval_aux((n+1)>>1, x, &t1, &t2);
      78          14 :   if (n != 1) t2 = gsub(gmul(gmul2n(t1,1), t2), x);
      79          35 :   for (i = 1; i <= v; i++) t2 = gadd(gmul2n(gsqr(t2), 1), gen_m1);
      80          14 :   return gerepileupto(av, t2);
      81             : }
      82             : 
      83             : /* Chebychev  polynomial of the second kind U(n,x): the coefficient in front of
      84             :  * x^(n-2*m) is (-1)^m * 2^(n-2m)*(n-m)!/m!/(n-2m)!  for m=0,1,...,n/2 */
      85             : GEN
      86        2135 : polchebyshev2(long n, long v)
      87             : {
      88             :   pari_sp av;
      89             :   GEN q, a, r;
      90             :   long m;
      91        2135 :   int neg = 0;
      92             : 
      93        2135 :   if (v<0) v = 0;
      94             :   /* polchebyshev(-n,2) = -polchebyshev(n-2,2) */
      95        2135 :   if (n < 0) {
      96        1050 :     if (n == -1) return zeropol(v);
      97        1029 :     neg = 1; n = -n-2;
      98             :   }
      99        2114 :   if (n==0) return neg ? scalar_ZX_shallow(gen_m1, v): pol_1(v);
     100             : 
     101        2072 :   q = cgetg(n+3, t_POL); r = q + n+2;
     102        2072 :   a = int2n(n);
     103        2072 :   if (neg) togglesign(a);
     104        2072 :   gel(r--,0) = a;
     105        2072 :   gel(r--,0) = gen_0;
     106       30807 :   for (m=1; 2*m<= n; m++)
     107             :   {
     108       28735 :     av = avma;
     109       28735 :     a = diviuuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m, n-m+1);
     110       28735 :     togglesign(a); a = gerepileuptoint(av, a);
     111       28735 :     gel(r--,0) = a;
     112       28735 :     gel(r--,0) = gen_0;
     113             :   }
     114        2072 :   q[1] = evalsigne(1) | evalvarn(v);
     115        2072 :   return q;
     116             : }
     117             : static void
     118          91 : polchebyshev2_eval_aux(long n, GEN x, GEN *pu1, GEN *pu2)
     119             : {
     120             :   GEN u1, u2, u, mu1;
     121          91 :   if (n == 1) { *pu1 = gen_1; *pu2 = gmul2n(x,1); return; }
     122          70 :   if (n == 0) { *pu1 = gen_0; *pu2 = gen_1; return; }
     123          70 :   polchebyshev2_eval_aux(n >> 1, x, &u1, &u2);
     124          70 :   mu1 = gneg(u1);
     125          70 :   u = gmul(gadd(u2,u1), gadd(u2,mu1));
     126          70 :   if (odd(n)) { *pu1 = u; *pu2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1)); }
     127          35 :   else        { *pu2 = u; *pu1 = gmul(gmul2n(u1,1), gadd(u2, gmul(x,mu1))); }
     128             : }
     129             : static GEN
     130          35 : polchebyshev2_eval(long n, GEN x)
     131             : {
     132             :   GEN u1, u2, mu1;
     133          35 :   long neg = 0;
     134             :   pari_sp av;
     135             : 
     136          35 :   if (n < 0) {
     137          14 :     if (n == -1) return gen_0;
     138           7 :     neg = 1; n = -n-2;
     139             :   }
     140          28 :   if (n==0) return neg ? gen_m1: gen_1;
     141          21 :   av = avma;
     142          21 :   polchebyshev2_eval_aux(n>>1, x, &u1, &u2);
     143          21 :   mu1 = gneg(u1);
     144          21 :   if (odd(n)) u2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1));
     145          14 :   else        u2 = gmul(gadd(u2,u1), gadd(u2,mu1));
     146          21 :   if (neg) u2 = gneg(u2);
     147          21 :   return gerepileupto(av, u2);
     148             : }
     149             : 
     150             : GEN
     151        4284 : polchebyshev(long n, long kind, long v)
     152             : {
     153        4284 :   switch (kind)
     154             :   {
     155        2149 :     case 1: return polchebyshev1(n, v);
     156        2135 :     case 2: return polchebyshev2(n, v);
     157           0 :     default: pari_err_FLAG("polchebyshev");
     158             :   }
     159             :   return NULL; /* LCOV_EXCL_LINE */
     160             : }
     161             : GEN
     162        4333 : polchebyshev_eval(long n, long kind, GEN x)
     163             : {
     164        4333 :   if (!x) return polchebyshev(n, kind, 0);
     165          63 :   if (gequalX(x)) return polchebyshev(n, kind, varn(x));
     166          49 :   switch (kind)
     167             :   {
     168          14 :     case 1: return polchebyshev1_eval(n, x);
     169          35 :     case 2: return polchebyshev2_eval(n, x);
     170           0 :     default: pari_err_FLAG("polchebyshev");
     171             :   }
     172             :   return NULL; /* LCOV_EXCL_LINE */
     173             : }
     174             : 
     175             : /* Hermite polynomial H(n,x):  H(n+1) = 2x H(n) - 2n H(n-1)
     176             :  * The coefficient in front of x^(n-2*m) is
     177             :  * (-1)^m * n! * 2^(n-2m)/m!/(n-2m)!  for m=0,1,...,n/2.. */
     178             : GEN
     179        1442 : polhermite(long n, long v)
     180             : {
     181             :   long m;
     182             :   pari_sp av;
     183             :   GEN q,a,r;
     184             : 
     185        1442 :   if (v<0) v = 0;
     186        1442 :   if (n==0) return pol_1(v);
     187             : 
     188        1435 :   q = cgetg(n+3, t_POL); r = q + n+2;
     189        1435 :   a = int2n(n);
     190        1435 :   gel(r--,0) = a;
     191        1435 :   gel(r--,0) = gen_0;
     192       40327 :   for (m=1; 2*m<= n; m++)
     193             :   {
     194       38892 :     av = avma;
     195       38892 :     a = diviuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m);
     196       38892 :     togglesign(a);
     197       38892 :     gel(r--,0) = a = gerepileuptoint(av, a);
     198       38892 :     gel(r--,0) = gen_0;
     199             :   }
     200        1435 :   q[1] = evalsigne(1) | evalvarn(v);
     201        1435 :   return q;
     202             : }
     203             : static void
     204          21 : err_hermite(long n)
     205          21 : { pari_err_DOMAIN("polhermite", "degree", "<", gen_0, stoi(n)); }
     206             : GEN
     207        1477 : polhermite_eval0(long n, GEN x, long flag)
     208             : {
     209             :   long i;
     210             :   pari_sp av, av2;
     211             :   GEN x2, u, v;
     212             : 
     213        1477 :   if (n < 0) err_hermite(n);
     214        1470 :   if (!x || gequalX(x))
     215             :   {
     216        1442 :     long v = x? varn(x): 0;
     217        1442 :     if (flag)
     218             :     {
     219          14 :       if (!n) err_hermite(-1);
     220           7 :       retmkvec2(polhermite(n-1,v),polhermite(n,v));
     221             :     }
     222        1428 :     return polhermite(n, v);
     223             :   }
     224          28 :   if (n==0)
     225             :   {
     226           7 :     if (flag) err_hermite(-1);
     227           0 :     return gen_1;
     228             :   }
     229          21 :   if (n==1)
     230             :   {
     231           0 :     if (flag) retmkvec2(gen_1, gmul2n(x,1));
     232           0 :     return gmul2n(x,1);
     233             :   }
     234          21 :   av = avma; x2 = gmul2n(x,1); v = gen_1; u = x2;
     235          21 :   av2= avma;
     236        7070 :   for (i=1; i<n; i++)
     237             :   { /* u = H_i(x), v = H_{i-1}(x), compute t = H_{i+1}(x) */
     238             :     GEN t;
     239        7049 :     if ((i & 0xff) == 0) gerepileall(av2,2,&u, &v);
     240        7049 :     t = gsub(gmul(x2, u), gmulsg(2*i,v));
     241        7049 :     v = u; u = t;
     242             :   }
     243          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     244          14 :   return gerepileupto(av, u);
     245             : }
     246             : GEN
     247           0 : polhermite_eval(long n, GEN x) { return polhermite_eval0(n, x, 0); }
     248             : 
     249             : /* Legendre polynomial
     250             :  * L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1)
     251             :  * L(n) = 2^-n sum_{k=0}^{n/2} a_k x^(n-2k)
     252             :  *   where a_k = (-1)^k (2n-2k)! / k! (n-k)! (n-2k)! is an integer
     253             :  *   and a_0 = binom(2n,n), a_k / a_{k-1} = - (n-2k+1)(n-2k+2) / 2k (2n-2k+1) */
     254             : GEN
     255        2163 : pollegendre(long n, long v)
     256             : {
     257             :   long k, l;
     258             :   pari_sp av;
     259             :   GEN a, r, q;
     260             : 
     261        2163 :   if (v<0) v = 0;
     262             :   /* pollegendre(-n) = pollegendre(n-1) */
     263        2163 :   if (n < 0) n = -n-1;
     264        2163 :   if (n==0) return pol_1(v);
     265        2121 :   if (n==1) return pol_x(v);
     266             : 
     267        2079 :   av = avma;
     268        2079 :   q = cgetg(n+3, t_POL); r = q + n+2;
     269        2079 :   gel(r--,0) = a = binomialuu(n<<1,n);
     270        2079 :   gel(r--,0) = gen_0;
     271       31423 :   for (k=1,l=n; l>1; k++,l-=2)
     272             :   { /* l = n-2*k+2 */
     273       29344 :     av = avma;
     274       29344 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     275       29344 :     togglesign(a); a = gerepileuptoint(av, a);
     276       29344 :     gel(r--,0) = a;
     277       29344 :     gel(r--,0) = gen_0;
     278             :   }
     279        2079 :   q[1] = evalsigne(1) | evalvarn(v);
     280        2079 :   return gerepileupto(av, gmul2n(q,-n));
     281             : }
     282             : /* q such that Ln * 2^n = q(x^2) [n even] or x q(x^2) [n odd] */
     283             : GEN
     284           0 : pollegendre_reduced(long n, long v)
     285             : {
     286             :   long k, l, N;
     287             :   pari_sp av;
     288             :   GEN a, r, q;
     289             : 
     290           0 :   if (v<0) v = 0;
     291             :   /* pollegendre(-n) = pollegendre(n-1) */
     292           0 :   if (n < 0) n = -n-1;
     293           0 :   if (n<=1) return n? scalarpol_shallow(gen_2,v): pol_1(v);
     294             : 
     295           0 :   N = n >> 1;
     296           0 :   q = cgetg(N+3, t_POL); r = q + N+2;
     297           0 :   gel(r--,0) = a = binomialuu(n<<1,n);
     298           0 :   for (k=1,l=n; l>1; k++,l-=2)
     299             :   { /* l = n-2*k+2 */
     300           0 :     av = avma;
     301           0 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     302           0 :     togglesign(a);
     303           0 :     gel(r--,0) = a = gerepileuptoint(av, a);
     304             :   }
     305           0 :   q[1] = evalsigne(1) | evalvarn(v);
     306           0 :   return q;
     307             : }
     308             : 
     309             : GEN
     310        2177 : pollegendre_eval0(long n, GEN x, long flag)
     311             : {
     312             :   pari_sp av;
     313             :   GEN u, v;
     314             :   long i;
     315             : 
     316        2177 :   if (n < 0) n = -n-1; /* L(-n) = L(n-1) */
     317             :   /* n >= 0 */
     318        2177 :   if (flag && flag != 1) pari_err_FLAG("pollegendre");
     319        2177 :   if (!x || gequalX(x))
     320             :   {
     321        2156 :     long v = x? varn(x): 0;
     322        2156 :     if (flag) retmkvec2(pollegendre(n-1,v), pollegendre(n,v));
     323        2149 :     return pollegendre(n, v);
     324             :   }
     325          21 :   if (n==0)
     326             :   {
     327           0 :     if (flag) retmkvec2(gen_1, gcopy(x));
     328           0 :     return gen_1;
     329             :   }
     330          21 :   if (n==1)
     331             :   {
     332           0 :     if (flag) retmkvec2(gcopy(x), gen_1);
     333           0 :     return gcopy(x);
     334             :   }
     335          21 :   av = avma; v = gen_1; u = x;
     336        7070 :   for (i=1; i<n; i++)
     337             :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     338             :     GEN t;
     339        7049 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     340        7049 :     t = gdivgs(gsub(gmul(gmulsg(2*i+1,x), u), gmulsg(i,v)), i+1);
     341        7049 :     v = u; u = t;
     342             :   }
     343          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     344          14 :   return gerepileupto(av, u);
     345             : }
     346             : GEN
     347           0 : pollegendre_eval(long n, GEN x) { return pollegendre_eval0(n, x, 0); }
     348             : 
     349             : /* Laguerre polynomial
     350             :  * L0^a = 1; L1^a = -X+a+1;
     351             :  * (n+1)*L^a(n+1) = (-X+(2*n+a+1))*L^a(n) - (n+a)*L^a(n-1)
     352             :  * L^a(n) = sum_{k=0}^n (-1)^k * binom(n+a,n-k) * x^k/k! */
     353             : GEN
     354        2128 : pollaguerre(long n, GEN a, long v)
     355             : {
     356        2128 :   pari_sp av = avma;
     357        2128 :   GEN L = cgetg(n+3, t_POL), c1 = gen_1, c2 = mpfact(n);
     358             :   long i;
     359             : 
     360        2128 :   L[1] = evalsigne(1) | evalvarn(v);
     361        2128 :   if (odd(n)) togglesign_safe(&c2);
     362      117404 :   for (i = n; i >= 0; i--)
     363             :   {
     364      115276 :     gel(L, i+2) = gdiv(c1, c2);
     365      115276 :     if (i)
     366             :     {
     367      113148 :       c2 = divis(c2,-i);
     368      113148 :       c1 = gdivgs(gmul(c1, gaddsg(i,a)), n+1-i);
     369             :     }
     370             :   }
     371        2128 :   return gerepilecopy(av, L);
     372             : }
     373             : static void
     374          21 : err_lag(long n)
     375          21 : { pari_err_DOMAIN("pollaguerre", "degree", "<", gen_0, stoi(n)); }
     376             : GEN
     377        2163 : pollaguerre_eval0(long n, GEN a, GEN x, long flag)
     378             : {
     379        2163 :   pari_sp av = avma;
     380             :   long i;
     381             :   GEN v, u;
     382             : 
     383        2163 :   if (n < 0) err_lag(n);
     384        2156 :   if (flag && flag != 1) pari_err_FLAG("pollaguerre");
     385        2156 :   if (!a) a = gen_0;
     386        2156 :   if (!x || gequalX(x))
     387             :   {
     388        2128 :     long v = x? varn(x): 0;
     389        2128 :     if (flag)
     390             :     {
     391          14 :       if (!n) err_lag(-1);
     392           7 :       retmkvec2(pollaguerre(n-1,a,v), pollaguerre(n,a,v));
     393             :     }
     394        2114 :     return pollaguerre(n,a,v);
     395             :   }
     396          28 :   if (n==0)
     397             :   {
     398           7 :     if (flag) err_lag(-1);
     399           0 :     return gen_1;
     400             :   }
     401          21 :   if (n==1)
     402             :   {
     403           0 :     if (flag) retmkvec2(gsub(gaddgs(a,1),x), gen_1);
     404           0 :     return gsub(gaddgs(a,1),x);
     405             :   }
     406          21 :   av = avma; v = gen_1; u = gsub(gaddgs(a,1),x);
     407        7070 :   for (i=1; i<n; i++)
     408             :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     409             :     GEN t;
     410        7049 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     411        7049 :     t = gdivgs(gsub(gmul(gsub(gaddsg(2*i+1,a),x), u), gmul(gaddsg(i,a),v)), i+1);
     412        7049 :     v = u; u = t;
     413             :   }
     414          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     415          14 :   return gerepileupto(av, u);
     416             : }
     417             : GEN
     418           0 : pollaguerre_eval(long n, GEN x, GEN a) { return pollaguerre_eval0(n, x, a, 0); }
     419             : 
     420             : /* polcyclo(p) = X^(p-1) + ... + 1 */
     421             : static GEN
     422      337995 : polcyclo_prime(long p, long v)
     423             : {
     424      337995 :   GEN T = cgetg(p+2, t_POL);
     425             :   long i;
     426      337995 :   T[1] = evalsigne(1) | evalvarn(v);
     427     1239652 :   for (i = 2; i < p+2; i++) gel(T,i) = gen_1;
     428      337995 :   return T;
     429             : }
     430             : 
     431             : /* cyclotomic polynomial */
     432             : GEN
     433      453242 : polcyclo(long n, long v)
     434             : {
     435             :   long s, q, i, l;
     436      453242 :   pari_sp av=avma;
     437             :   GEN T, P;
     438             : 
     439      453242 :   if (v<0) v = 0;
     440      453242 :   if (n < 3)
     441      115247 :     switch(n)
     442             :     {
     443       27482 :       case 1: return deg1pol_shallow(gen_1, gen_m1, v);
     444       87765 :       case 2: return deg1pol_shallow(gen_1, gen_1, v);
     445           0 :       default: pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     446             :     }
     447      337995 :   P = gel(factoru(n), 1); l = lg(P);
     448      337995 :   s = P[1]; T = polcyclo_prime(s, v);
     449      543382 :   for (i = 2; i < l; i++)
     450             :   { /* Phi_{np}(X) = Phi_n(X^p) / Phi_n(X) */
     451      205387 :     s *= P[i];
     452      205387 :     T = RgX_div(RgX_inflate(T, P[i]), T);
     453             :   }
     454             :   /* s = squarefree part of n */
     455      337995 :   q = n / s;
     456      337995 :   if (q == 1) return gerepileupto(av, T);
     457      178866 :   return gerepilecopy(av, RgX_inflate(T,q));
     458             : }
     459             : 
     460             : /* cyclotomic polynomial */
     461             : GEN
     462        8799 : polcyclo_eval(long n, GEN x)
     463             : {
     464        8799 :   pari_sp av= avma;
     465             :   GEN P, md, xd, yneg, ypos;
     466             :   long l, s, i, j, q, tx;
     467        8799 :   long root_of_1 = 0;
     468             : 
     469        8799 :   if (!x) return polcyclo(n, 0);
     470        7287 :   tx = typ(x);
     471        7287 :   if (gequalX(x)) return polcyclo(n, varn(x));
     472        6706 :   if (n <= 0) pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     473        6706 :   if (n == 1) return gsubgs(x, 1);
     474        6706 :   if (tx == t_INT && !signe(x)) return gen_1;
     475        7259 :   while ((n & 3) == 0) { n >>= 1; x = gsqr(x); } /* Phi_4n(x) = Phi_2n(x^2) */
     476             :   /* n not divisible by 4 */
     477        6706 :   if (n == 2) return gerepileupto(av, gaddgs(x,1));
     478         875 :   if (!odd(n)) { n >>= 1; x = gneg(x); } /* Phi_2n(x) = Phi_n(-x) for n>1 odd */
     479             :   /* n odd > 2.  s largest squarefree divisor of n */
     480         875 :   P = gel(factoru(n), 1); s = zv_prod(P);
     481             :   /* replace n by largest squarefree divisor */
     482         875 :   q = n/s; if (q != 1) { x = gpowgs(x, q); n = s; }
     483         875 :   l = lg(P)-1;
     484             :   /* n squarefree odd > 2, l distinct prime divisors. Now handle x = 1 or -1 */
     485         875 :   if (tx == t_INT) { /* shortcut */
     486         805 :     if (is_pm1(x))
     487             :     {
     488          56 :       set_avma(av);
     489          56 :       if (signe(x) > 0 && l == 1) return utoipos(P[1]);
     490          35 :       return gen_1;
     491             :     }
     492             :   } else {
     493          70 :     if (gequal1(x))
     494             :     { /* n is prime, return n; multiply by x to keep the type */
     495          14 :       if (l == 1) return gerepileupto(av, gmulgs(x,n));
     496           7 :       return gerepilecopy(av, x); /* else 1 */
     497             :     }
     498          56 :     if (gequalm1(x)) return gerepileupto(av, gneg(x)); /* -1 */
     499             :   }
     500             :   /* Heuristic: evaluation will probably not improve things */
     501         798 :   if (tx == t_POL || tx == t_MAT || lg(x) > n)
     502          17 :     return gerepileupto(av, poleval(polcyclo(n,0), x));
     503             : 
     504         781 :   xd = cgetg((1L<<l) + 1, t_VEC); /* the x^d, where d | n */
     505         781 :   md = cgetg((1L<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
     506         781 :   gel(xd, 1) = x;
     507         781 :   md[1] = 1;
     508             :   /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d).
     509             :    * If x has exact order D, n = Dq, then the result is 0 if q = 1. Otherwise
     510             :    * the factors with x^d-1, D|d are omitted and we multiply at the end by
     511             :    *   prod_{d | q} d^mu(q/d) = q if prime, 1 otherwise */
     512             :   /* We store the factors with mu(d)= 1 (resp.-1) in ypos (resp yneg).
     513             :    * At the end we return ypos/yneg if mu(n)=1 and yneg/ypos if mu(n)=-1 */
     514         781 :   ypos = gsubgs(x,1);
     515         781 :   yneg = gen_1;
     516        1660 :   for (i = 1; i <= l; i++)
     517             :   {
     518         879 :     long ti = 1L<<(i-1), p = P[i];
     519        1870 :     for (j = 1; j <= ti; j++) {
     520         991 :       GEN X = gpowgs(gel(xd,j), p), t = gsubgs(X,1);
     521         991 :       gel(xd,ti+j) = X;
     522         991 :       md[ti+j] = -md[j];
     523         991 :       if (gequal0(t))
     524             :       { /* x^d = 1; root_of_1 := the smallest index ti+j such that X == 1
     525             :         * (whose bits code d: bit i-1 is set iff P[i] | d). If no such index
     526             :         * exists, then root_of_1 remains 0. Do not multiply with X-1 if X = 1,
     527             :         * we handle these factors at the end */
     528          28 :         if (!root_of_1) root_of_1 = ti+j;
     529             :       }
     530             :       else
     531             :       {
     532         963 :         if (md[ti+j] == 1) ypos = gmul(ypos, t);
     533         872 :         else               yneg = gmul(yneg, t);
     534             :       }
     535             :     }
     536             :   }
     537         781 :   ypos = odd(l)? gdiv(yneg,ypos): gdiv(ypos,yneg);
     538         781 :   if (root_of_1)
     539             :   {
     540          21 :     GEN X = gel(xd,(1<<l)); /* = x^n = 1 */
     541          21 :     long bitmask_q = (1<<l) - root_of_1;
     542             :     /* bitmask_q encodes q = n/d: bit (i-1) is 1 iff P[i] | q */
     543             : 
     544             :     /* x is a root of unity.  If bitmask_q = 0, then x was a primitive n-th
     545             :      * root of 1 and the result is zero. Return X - 1 to preserve type. */
     546          21 :     if (!bitmask_q) return gerepileupto(av, gsubgs(X, 1));
     547             :     /* x is a primitive d-th root of unity, where d|n and d<n: we
     548             :      * must multiply ypos by if(isprime(n/d), n/d, 1) */
     549           7 :     ypos = gmul(ypos, X); /* multiply by X = 1 to preserve type */
     550             :     /* If bitmask_q = 1<<(i-1) for some i <= l, then q == P[i] and we multiply
     551             :      * by P[i]; otherwise q is composite and nothing more needs to be done */
     552           7 :     if (!(bitmask_q & (bitmask_q-1))) /* detects power of 2, since bitmask!=0 */
     553             :     {
     554           7 :       i = vals(bitmask_q)+1; /* q = P[i] */
     555           7 :       ypos = gmulgs(ypos, P[i]);
     556             :     }
     557             :   }
     558         767 :   return gerepileupto(av, ypos);
     559             : }
     560             : /********************************************************************/
     561             : /**                                                                **/
     562             : /**                  HILBERT & PASCAL MATRICES                     **/
     563             : /**                                                                **/
     564             : /********************************************************************/
     565             : GEN
     566         133 : mathilbert(long n) /* Hilbert matrix of order n */
     567             : {
     568             :   long i,j;
     569             :   GEN p;
     570             : 
     571         133 :   if (n < 0) pari_err_DOMAIN("mathilbert", "dimension", "<", gen_0, stoi(n));
     572         133 :   p = cgetg(n+1,t_MAT);
     573        1120 :   for (j=1; j<=n; j++)
     574             :   {
     575         987 :     gel(p,j) = cgetg(n+1,t_COL);
     576       16583 :     for (i=1+(j==1); i<=n; i++)
     577       15596 :       gcoeff(p,i,j) = mkfrac(gen_1, utoipos(i+j-1));
     578             :   }
     579         133 :   if (n) gcoeff(p,1,1) = gen_1;
     580         133 :   return p;
     581             : }
     582             : 
     583             : /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
     584             : GEN
     585        2968 : matqpascal(long n, GEN q)
     586             : {
     587             :   long i, j, I;
     588        2968 :   pari_sp av = avma;
     589        2968 :   GEN m, qpow = NULL; /* gcc -Wall */
     590             : 
     591        2968 :   if (n < -1)  pari_err_DOMAIN("matpascal", "n", "<", gen_m1, stoi(n));
     592        2968 :   n++; m = cgetg(n+1,t_MAT);
     593       34797 :   for (j=1; j<=n; j++) gel(m,j) = cgetg(n+1,t_COL);
     594        2968 :   if (q)
     595             :   {
     596          42 :     I = (n+1)/2;
     597          42 :     if (I > 1) { qpow = new_chunk(I+1); gel(qpow,2)=q; }
     598          84 :     for (j=3; j<=I; j++) gel(qpow,j) = gmul(q, gel(qpow,j-1));
     599             :   }
     600       34797 :   for (i=1; i<=n; i++)
     601             :   {
     602       31829 :     I = (i+1)/2; gcoeff(m,i,1)= gen_1;
     603       31829 :     if (q)
     604             :     {
     605         483 :       for (j=2; j<=I; j++)
     606         238 :         gcoeff(m,i,j) = gadd(gmul(gel(qpow,j),gcoeff(m,i-1,j)),
     607         238 :                              gcoeff(m,i-1,j-1));
     608             :     }
     609             :     else
     610             :     {
     611      991893 :       for (j=2; j<=I; j++)
     612      960309 :         gcoeff(m,i,j) = addii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
     613             :     }
     614     1007678 :     for (   ; j<=i; j++) gcoeff(m,i,j) = gcoeff(m,i,i+1-j);
     615     1968225 :     for (   ; j<=n; j++) gcoeff(m,i,j) = gen_0;
     616             :   }
     617        2968 :   return gerepilecopy(av, m);
     618             : }
     619             : 
     620             : GEN
     621          56 : eulerianpol(long N, long v)
     622             : {
     623          56 :   pari_sp av = avma;
     624          56 :   long n, n2, k = 0;
     625             :   GEN A;
     626          56 :   if (v < 0) v = 0;
     627          56 :   if (N <= 0) pari_err_DOMAIN("eulerianpol", "index", "<=", gen_0, stoi(N));
     628          49 :   if (N == 1) return pol_1(v);
     629          42 :   if (N == 2) return deg1pol_shallow(gen_1, gen_1, v);
     630          35 :   A = cgetg(N+1, t_VEC);
     631          35 :   gel(A,1) = gen_1; gel(A,2) = gen_1; /* A_2 = x+1 */
     632         567 :   for (n = 3; n <= N; n++)
     633             :   { /* A(n,k) = (n-k)A(n-1,k-1) + (k+1)A(n-1,k) */
     634         532 :     n2 = n >> 1;
     635         532 :     if (odd(n)) gel(A,n2+1) = mului(n+1, gel(A,n2));
     636        8652 :     for (k = n2-1; k; k--)
     637        8120 :       gel(A,k+1) = addii(mului(n-k, gel(A,k)), mului(k+1, gel(A,k+1)));
     638         532 :     if (gc_needed(av,1))
     639             :     {
     640           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"eulerianpol, %ld/%ld",n,N);
     641           0 :       for (k = odd(n)? n2+1: n2; k < N; k++) gel(A,k+1) = gen_0;
     642           0 :       A = gerepilecopy(av, A);
     643             :     }
     644             :   }
     645          35 :   k = N >> 1; if (odd(N)) k++;
     646         329 :   for (; k < N; k++) gel(A,k+1) = gel(A, N-k);
     647          35 :   return gerepilecopy(av, RgV_to_RgX(A, v));
     648             : }
     649             : 
     650             : /******************************************************************/
     651             : /**                                                              **/
     652             : /**                       PRECISION CHANGES                      **/
     653             : /**                                                              **/
     654             : /******************************************************************/
     655             : 
     656             : GEN
     657          84 : gprec(GEN x, long d)
     658             : {
     659          84 :   pari_sp av = avma;
     660          84 :   if (d <= 0) pari_err_DOMAIN("gprec", "precision", "<=", gen_0, stoi(d));
     661          84 :   return gerepilecopy(av, gprec_w(x, ndec2prec(d)));
     662             : }
     663             : 
     664             : /* not GC-safe; precision given in word length (including codewords) */
     665             : GEN
     666     3375784 : gprec_w(GEN x, long pr)
     667             : {
     668             :   long lx, i;
     669             :   GEN y;
     670             : 
     671     3375784 :   switch(typ(x))
     672             :   {
     673     2611143 :     case t_REAL:
     674     2611143 :       if (signe(x)) return realprec(x) != pr? rtor(x,pr): x;
     675       26236 :       i = -prec2nbits(pr);
     676       26236 :       if (i < expo(x)) return real_0_bit(i);
     677       25294 :       y = cgetr(2); y[1] = x[1]; return y;
     678      435246 :     case t_COMPLEX:
     679      435246 :       y = cgetg(3, t_COMPLEX);
     680      435246 :       gel(y,1) = gprec_w(gel(x,1),pr);
     681      435246 :       gel(y,2) = gprec_w(gel(x,2),pr);
     682      435246 :       break;
     683      102400 :    case t_POL: case t_SER:
     684      102400 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     685      790421 :       for (i=2; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     686      102400 :       break;
     687       75519 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     688       75519 :       y = cgetg_copy(x, &lx);
     689      347166 :       for (i=1; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     690       75519 :       break;
     691      151476 :     default: return x;
     692             :   }
     693      613165 :   return y;
     694             : }
     695             : /* not GC-safe; precision given in word length (including codewords) */
     696             : GEN
     697     4169268 : gprec_wensure(GEN x, long pr)
     698             : {
     699             :   long lx, i;
     700             :   GEN y;
     701             : 
     702     4169268 :   switch(typ(x))
     703             :   {
     704     3613895 :     case t_REAL:
     705     3613895 :       if (signe(x)) return realprec(x) < pr? rtor(x,pr): x;
     706        9300 :       i = -prec2nbits(pr);
     707        9300 :       if (i < expo(x)) return real_0_bit(i);
     708        7004 :       y = cgetr(2); y[1] = x[1]; return y;
     709      130566 :     case t_COMPLEX:
     710      130566 :       y = cgetg(3, t_COMPLEX);
     711      130566 :       gel(y,1) = gprec_wensure(gel(x,1),pr);
     712      130566 :       gel(y,2) = gprec_wensure(gel(x,2),pr);
     713      130566 :       break;
     714       49784 :    case t_POL: case t_SER:
     715       49784 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     716      868336 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wensure(gel(x,i),pr);
     717       49784 :       break;
     718       82605 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     719       82605 :       y = cgetg_copy(x, &lx);
     720     1298434 :       for (i=1; i<lx; i++) gel(y,i) = gprec_wensure(gel(x,i),pr);
     721       82605 :       break;
     722      292418 :     default: return x;
     723             :   }
     724      262955 :   return y;
     725             : }
     726             : 
     727             : /* not GC-safe; precision given in word length (including codewords),
     728             :  * truncate mantissa to precision 'pr' but never increase it */
     729             : GEN
     730     1491217 : gprec_wtrunc(GEN x, long pr)
     731             : {
     732             :   long lx, i;
     733             :   GEN y;
     734             : 
     735     1491217 :   switch(typ(x))
     736             :   {
     737     1247593 :     case t_REAL:
     738     1247593 :       return (signe(x) && realprec(x) > pr)? rtor(x,pr): x;
     739      166198 :     case t_COMPLEX:
     740      166198 :       y = cgetg(3, t_COMPLEX);
     741      166198 :       gel(y,1) = gprec_wtrunc(gel(x,1),pr);
     742      166198 :       gel(y,2) = gprec_wtrunc(gel(x,2),pr);
     743      166198 :       break;
     744        4179 :     case t_POL:
     745             :     case t_SER:
     746        4179 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     747       25774 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     748        4179 :       break;
     749       22060 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     750       22060 :       y = cgetg_copy(x, &lx);
     751      195294 :       for (i=1; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     752       22060 :       break;
     753       51187 :     default: return x;
     754             :   }
     755      192437 :   return y;
     756             : }
     757             : 
     758             : /********************************************************************/
     759             : /**                                                                **/
     760             : /**                      SERIES TRANSFORMS                         **/
     761             : /**                                                                **/
     762             : /********************************************************************/
     763             : /**                  LAPLACE TRANSFORM (OF A SERIES)               **/
     764             : /********************************************************************/
     765             : static GEN
     766          14 : serlaplace(GEN x)
     767             : {
     768          14 :   long i, l = lg(x), e = valp(x);
     769          14 :   GEN t, y = cgetg(l,t_SER);
     770          14 :   if (e < 0) pari_err_DOMAIN("laplace","valuation","<",gen_0,stoi(e));
     771          14 :   t = mpfact(e); y[1] = x[1];
     772         154 :   for (i=2; i<l; i++)
     773             :   {
     774         140 :     gel(y,i) = gmul(t, gel(x,i));
     775         140 :     e++; t = mului(e,t);
     776             :   }
     777          14 :   return y;
     778             : }
     779             : static GEN
     780          14 : pollaplace(GEN x)
     781             : {
     782          14 :   long i, e = 0, l = lg(x);
     783          14 :   GEN t = gen_1, y = cgetg(l,t_POL);
     784          14 :   y[1] = x[1];
     785          63 :   for (i=2; i<l; i++)
     786             :   {
     787          49 :     gel(y,i) = gmul(t, gel(x,i));
     788          49 :     e++; t = mului(e,t);
     789             :   }
     790          14 :   return y;
     791             : }
     792             : GEN
     793          35 : laplace(GEN x)
     794             : {
     795          35 :   pari_sp av = avma;
     796          35 :   switch(typ(x))
     797             :   {
     798          14 :     case t_POL: x = pollaplace(x); break;
     799          14 :     case t_SER: x = serlaplace(x); break;
     800           7 :     default: if (is_scalar_t(typ(x))) return gcopy(x);
     801           0 :              pari_err_TYPE("laplace",x);
     802             :   }
     803          28 :   return gerepilecopy(av, x);
     804             : }
     805             : 
     806             : /********************************************************************/
     807             : /**              CONVOLUTION PRODUCT (OF TWO SERIES)               **/
     808             : /********************************************************************/
     809             : GEN
     810           7 : convol(GEN x, GEN y)
     811             : {
     812           7 :   long j, lx, ly, ex, ey, vx = varn(x);
     813             :   GEN z;
     814             : 
     815           7 :   if (typ(x) != t_SER) pari_err_TYPE("convol",x);
     816           7 :   if (typ(y) != t_SER) pari_err_TYPE("convol",y);
     817           7 :   if (varn(y) != vx) pari_err_VAR("convol", x,y);
     818           7 :   ex = valp(x);
     819           7 :   ey = valp(y);
     820           7 :   if (ser_isexactzero(x))
     821           0 :     return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), maxss(ex,ey));
     822           7 :   lx = lg(x) + ex; x -= ex;
     823           7 :   ly = lg(y) + ey; y -= ey;
     824             :   /* inputs shifted: x[i] and y[i] now correspond to monomials of same degree */
     825           7 :   if (ly < lx) lx = ly; /* min length */
     826           7 :   if (ex < ey) ex = ey; /* max valuation */
     827           7 :   if (lx - ex < 3) return zeroser(vx, lx-2);
     828             : 
     829           7 :   z = cgetg(lx - ex, t_SER);
     830           7 :   z[1] = evalvalp(ex) | evalvarn(vx);
     831         119 :   for (j = ex+2; j<lx; j++) gel(z,j-ex) = gmul(gel(x,j),gel(y,j));
     832           7 :   return normalize(z);
     833             : }
     834             : 
     835             : /***********************************************************************/
     836             : /*               OPERATIONS ON DIRICHLET SERIES: *, /                  */
     837             : /* (+, -, scalar multiplication are done on the corresponding vectors) */
     838             : /***********************************************************************/
     839             : static long
     840        2016 : dirval(GEN x)
     841             : {
     842        2016 :   long i = 1, lx = lg(x);
     843        2037 :   while (i < lx && gequal0(gel(x,i))) i++;
     844        2016 :   return i;
     845             : }
     846             : 
     847             : GEN
     848         343 : dirmul(GEN x, GEN y)
     849             : {
     850         343 :   pari_sp av = avma, av2;
     851             :   long nx, ny, nz, dx, dy, i, j, k;
     852             :   GEN z;
     853             : 
     854         343 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirmul",x);
     855         343 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirmul",y);
     856         343 :   dx = dirval(x); nx = lg(x)-1;
     857         343 :   dy = dirval(y); ny = lg(y)-1;
     858         343 :   if (ny-dy < nx-dx) { swap(x,y); lswap(nx,ny); lswap(dx,dy); }
     859         343 :   nz = minss(nx*dy,ny*dx);
     860         343 :   y = RgV_kill0(y);
     861         343 :   av2 = avma;
     862         343 :   z = zerovec(nz);
     863       39291 :   for (j=dx; j<=nx; j++)
     864             :   {
     865       38948 :     GEN c = gel(x,j);
     866       38948 :     if (gequal0(c)) continue;
     867       17255 :     if (gequal1(c))
     868             :     {
     869       95291 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     870       89411 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gel(y,k));
     871             :     }
     872       11375 :     else if (gequalm1(c))
     873             :     {
     874        5649 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     875        4298 :         if (gel(y,k)) gel(z,i) = gsub(gel(z,i),gel(y,k));
     876             :     }
     877             :     else
     878             :     {
     879       46431 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     880       36407 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gmul(c,gel(y,k)));
     881             :     }
     882       17255 :     if (gc_needed(av2,3))
     883             :     {
     884           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirmul, %ld/%ld",j,nx);
     885           0 :       z = gerepilecopy(av2,z);
     886             :     }
     887             :   }
     888         343 :   return gerepilecopy(av,z);
     889             : }
     890             : 
     891             : GEN
     892         665 : dirdiv(GEN x, GEN y)
     893             : {
     894         665 :   pari_sp av = avma, av2;
     895             :   long nx,ny,nz, dx,dy, i,j,k;
     896             :   GEN p1;
     897             : 
     898         665 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirdiv",x);
     899         665 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirdiv",y);
     900         665 :   dx = dirval(x); nx = lg(x)-1;
     901         665 :   dy = dirval(y); ny = lg(y)-1;
     902         665 :   if (dy != 1 || !ny) pari_err_INV("dirdiv",y);
     903         665 :   nz = minss(nx,ny*dx);
     904         665 :   p1 = gel(y,1);
     905         665 :   if (gequal1(p1)) p1 = NULL; else y = gdiv(y,p1);
     906         665 :   y = RgV_kill0(y);
     907         665 :   av2 = avma;
     908         665 :   x = p1 ? gdiv(x,p1): leafcopy(x);
     909         672 :   for (j=1; j<dx; j++) gel(x,j) = gen_0;
     910         665 :   setlg(x,nz+1);
     911      384587 :   for (j=dx; j<=nz; j++)
     912             :   {
     913      383922 :     GEN c = gel(x,j);
     914      383922 :     if (gequal0(c)) continue;
     915      122500 :     if (gequal1(c))
     916             :     {
     917      826805 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     918      773101 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gel(y,k));
     919             :     }
     920       68796 :     else if (gequalm1(c))
     921             :     {
     922      667072 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     923      613018 :         if (gel(y,k)) gel(x,i) = gadd(gel(x,i),gel(y,k));
     924             :     }
     925             :     else
     926             :     {
     927       48286 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     928       33544 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gmul(c,gel(y,k)));
     929             :     }
     930      122500 :     if (gc_needed(av2,3))
     931             :     {
     932           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirdiv, %ld/%ld",j,nz);
     933           0 :       x = gerepilecopy(av2,x);
     934             :     }
     935             :   }
     936         665 :   return gerepilecopy(av,x);
     937             : }
     938             : 
     939             : /*******************************************************************/
     940             : /**                                                               **/
     941             : /**                       COMBINATORICS                           **/
     942             : /**                                                               **/
     943             : /*******************************************************************/
     944             : /**                      BINOMIAL COEFFICIENTS                    **/
     945             : /*******************************************************************/
     946             : GEN
     947       78050 : binomialuu(ulong n, ulong k)
     948             : {
     949       78050 :   pari_sp ltop = avma;
     950             :   GEN z;
     951       78050 :   if (k > n) return gen_0;
     952       78043 :   k = minuu(k,n-k);
     953       78043 :   if (!k) return gen_1;
     954       76419 :   if (k == 1) return utoipos(n);
     955       71659 :   z = diviiexact(mulu_interval(n-k+1, n), mulu_interval(2UL, k));
     956       71659 :   return gerepileuptoint(ltop,z);
     957             : }
     958             : 
     959             : GEN
     960      100345 : binomial(GEN n, long k)
     961             : {
     962             :   long i, prec;
     963             :   pari_sp av;
     964             :   GEN y;
     965             : 
     966      100345 :   if (k <= 1)
     967             :   {
     968       59311 :     if (is_noncalc_t(typ(n))) pari_err_TYPE("binomial",n);
     969       59311 :     if (k < 0) return gen_0;
     970       59311 :     if (k == 0) return gen_1;
     971       25963 :     return gcopy(n);
     972             :   }
     973       41034 :   av = avma;
     974       41034 :   if (typ(n) == t_INT)
     975             :   {
     976       40887 :     if (signe(n) > 0)
     977             :     {
     978       40873 :       GEN z = subiu(n,k);
     979       40873 :       if (cmpis(z,k) < 0)
     980             :       {
     981         924 :         k = itos(z); set_avma(av);
     982         924 :         if (k <= 1)
     983             :         {
     984         371 :           if (k < 0) return gen_0;
     985         371 :           if (k == 0) return gen_1;
     986         343 :           return icopy(n);
     987             :         }
     988             :       }
     989             :     }
     990             :     /* k > 1 */
     991       40516 :     if (lgefint(n) == 3 && signe(n) > 0)
     992             :     {
     993       40495 :       y = binomialuu(itou(n),(ulong)k);
     994       40495 :       return gerepileupto(av, y);
     995             :     }
     996             :     else
     997             :     {
     998          21 :       y = cgetg(k+1,t_VEC);
     999        1603 :       for (i=1; i<=k; i++) gel(y,i) = subiu(n,i-1);
    1000          21 :       y = ZV_prod(y);
    1001             :     }
    1002          21 :     y = diviiexact(y, mpfact(k));
    1003          21 :     return gerepileuptoint(av, y);
    1004             :   }
    1005             : 
    1006         147 :   prec = precision(n);
    1007         147 :   if (prec && k > 200 + 0.8*prec2nbits(prec)) {
    1008           7 :     GEN A = mpfactr(k, prec), B = ggamma(gsubgs(n,k-1), prec);
    1009           7 :     return gerepileupto(av, gdiv(ggamma(gaddgs(n,1), prec), gmul(A,B)));
    1010             :   }
    1011             : 
    1012         140 :   y = cgetg(k+1,t_VEC);
    1013       11571 :   for (i=1; i<=k; i++) gel(y,i) = gsubgs(n,i-1);
    1014         140 :   return gerepileupto(av, gdiv(RgV_prod(y), mpfact(k)));
    1015             : }
    1016             : 
    1017             : GEN
    1018         413 : binomial0(GEN x, GEN k)
    1019             : {
    1020         413 :   if (!k)
    1021             :   {
    1022          21 :     if (typ(x) != t_INT || signe(x) < 0) pari_err_TYPE("binomial", x);
    1023           7 :     return vecbinomial(itos(x));
    1024             :   }
    1025         392 :   if (typ(k) != t_INT) pari_err_TYPE("binomial", k);
    1026         385 :   return binomial(x, itos(k));
    1027             : }
    1028             : 
    1029             : /* Assume n >= 0, return bin, bin[k+1] = binomial(n, k) */
    1030             : GEN
    1031       31885 : vecbinomial(long n)
    1032             : {
    1033             :   long d, k;
    1034             :   GEN C;
    1035       31885 :   if (!n) return mkvec(gen_1);
    1036       31528 :   C = cgetg(n+2, t_VEC) + 1; /* C[k] = binomial(n, k) */
    1037       31528 :   gel(C,0) = gen_1;
    1038       31528 :   gel(C,1) = utoipos(n); d = (n + 1) >> 1;
    1039      158466 :   for (k=2; k <= d; k++)
    1040             :   {
    1041      126938 :     pari_sp av = avma;
    1042      126938 :     gel(C,k) = gerepileuptoint(av, diviuexact(mului(n-k+1, gel(C,k-1)), k));
    1043             :   }
    1044      175049 :   for (   ; k <= n; k++) gel(C,k) = gel(C,n-k);
    1045       31528 :   return C - 1;
    1046             : }
    1047             : 
    1048             : /********************************************************************/
    1049             : /**                  STIRLING NUMBERS                              **/
    1050             : /********************************************************************/
    1051             : /* Stirling number of the 2nd kind. The number of ways of partitioning
    1052             :    a set of n elements into m non-empty subsets. */
    1053             : GEN
    1054        1694 : stirling2(ulong n, ulong m)
    1055             : {
    1056        1694 :   pari_sp av = avma;
    1057             :   GEN s, bmk;
    1058             :   ulong k;
    1059        1694 :   if (n==0) return (m == 0)? gen_1: gen_0;
    1060        1694 :   if (m > n || m == 0) return gen_0;
    1061        1694 :   if (m==n) return gen_1;
    1062             :   /* k = 0 */
    1063        1694 :   bmk = gen_1; s  = powuu(m, n);
    1064       20314 :   for (k = 1; k <= ((m-1)>>1); ++k)
    1065             :   { /* bmk = binomial(m, k) */
    1066             :     GEN c, kn, mkn;
    1067       18620 :     bmk = diviuexact(mului(m-k+1, bmk), k);
    1068       18620 :     kn  = powuu(k, n); mkn = powuu(m-k, n);
    1069       18620 :     c = odd(m)? subii(mkn,kn): addii(mkn,kn);
    1070       18620 :     c = mulii(bmk, c);
    1071       18620 :     s = odd(k)? subii(s, c): addii(s, c);
    1072       18620 :     if (gc_needed(av,2))
    1073             :     {
    1074           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"stirling2");
    1075           0 :       gerepileall(av, 2, &s, &bmk);
    1076             :     }
    1077             :   }
    1078             :   /* k = m/2 */
    1079        1694 :   if (!odd(m))
    1080             :   {
    1081             :     GEN c;
    1082         805 :     bmk = diviuexact(mului(k+1, bmk), k);
    1083         805 :     c = mulii(bmk, powuu(k,n));
    1084         805 :     s = odd(k)? subii(s, c): addii(s, c);
    1085             :   }
    1086        1694 :   return gerepileuptoint(av, diviiexact(s, mpfact(m)));
    1087             : }
    1088             : 
    1089             : /* Stirling number of the first kind. Up to the sign, the number of
    1090             :    permutations of n symbols which have exactly m cycles. */
    1091             : GEN
    1092         154 : stirling1(ulong n, ulong m)
    1093             : {
    1094         154 :   pari_sp ltop=avma;
    1095             :   ulong k;
    1096             :   GEN s, t;
    1097         154 :   if (n < m) return gen_0;
    1098         154 :   else if (n==m) return gen_1;
    1099             :   /* t = binomial(n-1+k, m-1) * binomial(2n-m, n-m-k) */
    1100             :   /* k = n-m > 0 */
    1101         154 :   t = binomialuu(2*n-m-1, m-1);
    1102         154 :   s = mulii(t, stirling2(2*(n-m), n-m));
    1103         154 :   if (odd(n-m)) togglesign(s);
    1104        1547 :   for (k = n-m-1; k > 0; --k)
    1105             :   {
    1106             :     GEN c;
    1107        1393 :     t = diviuuexact(muluui(n-m+k+1, n+k+1, t), n+k, n-m-k);
    1108        1393 :     c = mulii(t, stirling2(n-m+k, k));
    1109        1393 :     s = odd(k)? subii(s, c): addii(s, c);
    1110        1393 :     if ((k & 0x1f) == 0) {
    1111          21 :       t = gerepileuptoint(ltop, t);
    1112          21 :       s = gerepileuptoint(avma, s);
    1113             :     }
    1114             :   }
    1115         154 :   return gerepileuptoint(ltop, s);
    1116             : }
    1117             : 
    1118             : GEN
    1119         301 : stirling(long n, long m, long flag)
    1120             : {
    1121         301 :   if (n < 0) pari_err_DOMAIN("stirling", "n", "<", gen_0, stoi(n));
    1122         301 :   if (m < 0) pari_err_DOMAIN("stirling", "m", "<", gen_0, stoi(m));
    1123         301 :   switch (flag)
    1124             :   {
    1125         154 :     case 1: return stirling1((ulong)n,(ulong)m);
    1126         147 :     case 2: return stirling2((ulong)n,(ulong)m);
    1127           0 :     default: pari_err_FLAG("stirling");
    1128             :   }
    1129             :   return NULL; /*LCOV_EXCL_LINE*/
    1130             : }
    1131             : 
    1132             : /*******************************************************************/
    1133             : /**                                                               **/
    1134             : /**                     RECIPROCAL POLYNOMIAL                     **/
    1135             : /**                                                               **/
    1136             : /*******************************************************************/
    1137             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1138             : GEN
    1139         154 : polrecip(GEN x)
    1140             : {
    1141         154 :   long tx = typ(x);
    1142         154 :   if (is_scalar_t(tx)) return gcopy(x);
    1143         147 :   if (tx != t_POL) pari_err_TYPE("polrecip",x);
    1144         147 :   return RgX_recip(x);
    1145             : }
    1146             : 
    1147             : /********************************************************************/
    1148             : /**                                                                **/
    1149             : /**                  POLYNOMIAL INTERPOLATION                      **/
    1150             : /**                                                                **/
    1151             : /********************************************************************/
    1152             : static GEN
    1153         210 : RgV_polint_fast(GEN X, GEN Y, long v)
    1154             : {
    1155             :   GEN p, pol;
    1156             :   long t, pa;
    1157         210 :   if (X) t = RgV_type2(X,Y, &p, &pol, &pa);
    1158          14 :   else   t = Rg_type(Y, &p, &pol, &pa);
    1159         210 :   if (t != t_INTMOD) return NULL;
    1160           7 :   Y = RgC_to_FpC(Y, p);
    1161           7 :   X = X? RgC_to_FpC(X, p): identity_ZV(lg(Y)-1);
    1162           7 :   return FpX_to_mod(FpV_polint(X, Y, p, v), p);
    1163             : }
    1164             : /* allow X = NULL for [1,...,n] */
    1165             : GEN
    1166         210 : RgV_polint(GEN X, GEN Y, long v)
    1167             : {
    1168         210 :   pari_sp av0 = avma, av;
    1169         210 :   GEN Q, P = NULL;
    1170         210 :   long i, l = lg(Y);
    1171         210 :   if ((Q = RgV_polint_fast(X,Y,v))) return Q;
    1172         203 :   if (!X) X = identity_ZV(l-1);
    1173         203 :   Q = roots_to_pol(X, v); av = avma;
    1174         462 :   for (i=1; i<l; i++)
    1175             :   {
    1176             :     GEN inv, T, dP;
    1177         259 :     if (gequal0(gel(Y,i))) continue;
    1178         161 :     T = RgX_div_by_X_x(Q, gel(X,i), NULL);
    1179         161 :     inv = ginv(poleval(T,gel(X,i)));
    1180         161 :     dP = RgX_Rg_mul(T, gmul(gel(Y,i),inv));
    1181         161 :     P = P? RgX_add(P, dP): dP;
    1182         161 :     if (gc_needed(av,2))
    1183             :     {
    1184           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgV_polint i = %ld/%ld", i, l-1);
    1185           0 :       P = gerepileupto(av, P);
    1186             :     }
    1187             :   }
    1188         203 :   if (!P) { set_avma(av); return zeropol(v); }
    1189         133 :   return gerepileupto(av0, P);
    1190             : }
    1191             : static int
    1192       17357 : inC(GEN x)
    1193             : {
    1194       17357 :   switch(typ(x)) {
    1195        1365 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD: return 1;
    1196       15992 :     default: return 0;
    1197             :   }
    1198             : }
    1199             : static long
    1200       16188 : check_dy(GEN X, GEN x, long n)
    1201             : {
    1202       16188 :   GEN D = NULL;
    1203       16188 :   long i, ns = 0;
    1204       16188 :   if (!inC(x)) return -1;
    1205        1176 :   for (i = 0; i < n; i++)
    1206             :   {
    1207         966 :     GEN t = gsub(x, gel(X,i));
    1208         966 :     if (!inC(t)) return -1;
    1209         952 :     t = gabs(t, DEFAULTPREC);
    1210         952 :     if (!D || gcmp(t,D) < 0) { ns = i; D = t; }
    1211             :   }
    1212             :   /* X[ns] is closest to x */
    1213         210 :   return ns;
    1214             : }
    1215             : /* X,Y are "spec" GEN vectors with n > 0 components ( at X[0], ... X[n-1] ) */
    1216             : GEN
    1217       16223 : polintspec(GEN X, GEN Y, GEN x, long n, long *pe)
    1218             : {
    1219             :   long i, m, ns;
    1220       16223 :   pari_sp av = avma, av2;
    1221       16223 :   GEN y, c, d, dy = NULL; /* gcc -Wall */
    1222             : 
    1223       16223 :   if (pe) *pe = -HIGHEXPOBIT;
    1224       16223 :   if (n == 1) return gmul(gel(Y,0), Rg_get_1(x));
    1225       16188 :   if (!X) X = identity_ZV(n) + 1;
    1226       16188 :   av2 = avma;
    1227       16188 :   ns = check_dy(X, x, n); if (ns < 0) { pe = NULL; ns = 0; }
    1228       16188 :   c = cgetg(n+1, t_VEC);
    1229       81031 :   d = cgetg(n+1, t_VEC); for (i=0; i<n; i++) gel(c,i+1) = gel(d,i+1) = gel(Y,i);
    1230       16188 :   y = gel(d,ns+1);
    1231             :   /* divided differences */
    1232       64836 :   for (m = 1; m < n; m++)
    1233             :   {
    1234      146238 :     for (i = 0; i < n-m; i++)
    1235             :     {
    1236       97590 :       GEN ho = gsub(gel(X,i),x), hp = gsub(gel(X,i+m),x), den = gsub(ho,hp);
    1237       97590 :       if (gequal0(den))
    1238             :       {
    1239           7 :         char *x1 = stack_sprintf("X[%ld]", i+1);
    1240           7 :         char *x2 = stack_sprintf("X[%ld]", i+m+1);
    1241           7 :         pari_err_DOMAIN("polinterpolate",x1,"=",strtoGENstr(x2), X);
    1242             :       }
    1243       97583 :       den = gdiv(gsub(gel(c,i+2),gel(d,i+1)), den);
    1244       97583 :       gel(c,i+1) = gmul(ho,den);
    1245       97583 :       gel(d,i+1) = gmul(hp,den);
    1246             :     }
    1247       48648 :     dy = (2*ns < n-m)? gel(c,ns+1): gel(d,ns--);
    1248       48648 :     y = gadd(y,dy);
    1249       48648 :     if (gc_needed(av2,2))
    1250             :     {
    1251           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"polint, %ld/%ld",m,n-1);
    1252           0 :       gerepileall(av2, 4, &y, &c, &d, &dy);
    1253             :     }
    1254             :   }
    1255       16181 :   if (pe && inC(dy)) *pe = gexpo(dy);
    1256       16181 :   return gerepileupto(av, y);
    1257             : }
    1258             : 
    1259             : GEN
    1260         315 : polint_i(GEN X, GEN Y, GEN t, long *pe)
    1261             : {
    1262         315 :   long lx = lg(X), vt;
    1263             : 
    1264         315 :   if (! is_vec_t(typ(X))) pari_err_TYPE("polinterpolate",X);
    1265         315 :   if (Y)
    1266             :   {
    1267         294 :     if (! is_vec_t(typ(Y))) pari_err_TYPE("polinterpolate",Y);
    1268         294 :     if (lx != lg(Y)) pari_err_DIM("polinterpolate");
    1269             :   }
    1270             :   else
    1271             :   {
    1272          21 :     Y = X;
    1273          21 :     X = NULL;
    1274             :   }
    1275         315 :   if (pe) *pe = -HIGHEXPOBIT;
    1276         315 :   vt = t? gvar(t): 0;
    1277         315 :   if (vt != NO_VARIABLE)
    1278             :   { /* formal interpolation */
    1279             :     pari_sp av;
    1280         210 :     long v0, vY = gvar(Y);
    1281             :     GEN P;
    1282         210 :     if (X) vY = varnmax(vY, gvar(X));
    1283             :     /* shortcut */
    1284         210 :     if (varncmp(vY, vt) > 0 && (!t || gequalX(t))) return RgV_polint(X, Y, vt);
    1285          84 :     av = avma;
    1286             :     /* first interpolate in high priority variable, then substitute t */
    1287          84 :     v0 = fetch_var_higher();
    1288          84 :     P = RgV_polint(X, Y, v0);
    1289          84 :     P = gsubst(P, v0, t? t: pol_x(0));
    1290          84 :     (void)delete_var();
    1291          84 :     return gerepileupto(av, P);
    1292             :   }
    1293             :   /* numerical interpolation */
    1294         105 :   if (lx == 1) return Rg_get_0(t);
    1295          91 :   return polintspec(X? X+1: NULL,Y+1,t,lx-1, pe);
    1296             : }
    1297             : GEN
    1298         315 : polint(GEN X, GEN Y, GEN t, GEN *pe)
    1299             : {
    1300             :   long e;
    1301         315 :   GEN p = polint_i(X, Y, t, &e);
    1302         308 :   if (pe) *pe = stoi(e);
    1303         308 :   return p;
    1304             : }
    1305             : 
    1306             : /********************************************************************/
    1307             : /**                                                                **/
    1308             : /**                       MODREVERSE                               **/
    1309             : /**                                                                **/
    1310             : /********************************************************************/
    1311             : static void
    1312           7 : err_reverse(GEN x, GEN T)
    1313             : {
    1314           7 :   pari_err_DOMAIN("modreverse","deg(minpoly(z))", "<", stoi(degpol(T)),
    1315             :                   mkpolmod(x,T));
    1316           0 : }
    1317             : 
    1318             : /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
    1319             : GEN
    1320         182 : RgXQ_reverse(GEN a, GEN T)
    1321             : {
    1322         182 :   pari_sp av = avma;
    1323         182 :   long n = degpol(T);
    1324             :   GEN y;
    1325             : 
    1326         182 :   if (n <= 1) {
    1327           7 :     if (n <= 0) return gcopy(a);
    1328           7 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1329             :   }
    1330         175 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1331         175 :   y = RgXV_to_RgM(RgXQ_powers(a,n-1,T), n);
    1332         175 :   y = RgM_solve(y, col_ei(n, 2));
    1333         175 :   if (!y) err_reverse(a,T);
    1334         168 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1335             : }
    1336             : GEN
    1337        1498 : QXQ_reverse(GEN a, GEN T)
    1338             : {
    1339        1498 :   pari_sp av = avma;
    1340        1498 :   long n = degpol(T);
    1341             :   GEN y;
    1342             : 
    1343        1498 :   if (n <= 1) {
    1344          14 :     if (n <= 0) return gcopy(a);
    1345          14 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1346             :   }
    1347        1484 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1348        1484 :   if (gequalX(a)) return gcopy(a);
    1349        1463 :   y = RgXV_to_RgM(QXQ_powers(a,n-1,T), n);
    1350        1463 :   y = QM_gauss(y, col_ei(n, 2));
    1351        1463 :   if (!y) err_reverse(a,T);
    1352        1463 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1353             : }
    1354             : 
    1355             : GEN
    1356          28 : modreverse(GEN x)
    1357             : {
    1358             :   long v, n;
    1359             :   GEN T, a;
    1360             : 
    1361          28 :   if (typ(x)!=t_POLMOD) pari_err_TYPE("modreverse",x);
    1362          28 :   T = gel(x,1); n = degpol(T); if (n <= 0) return gcopy(x);
    1363          21 :   a = gel(x,2);
    1364          21 :   v = varn(T);
    1365          21 :   retmkpolmod(RgXQ_reverse(a, T),
    1366             :               (n==1)? gsub(pol_x(v), a): RgXQ_charpoly(a, T, v));
    1367             : }
    1368             : 
    1369             : /********************************************************************/
    1370             : /**                                                                **/
    1371             : /**                          MERGESORT                             **/
    1372             : /**                                                                **/
    1373             : /********************************************************************/
    1374             : static int
    1375          77 : cmp_small(GEN x, GEN y) {
    1376          77 :   long a = (long)x, b = (long)y;
    1377          77 :   return a>b? 1: (a<b? -1: 0);
    1378             : }
    1379             : 
    1380             : static int
    1381      294994 : veccmp(void *data, GEN x, GEN y)
    1382             : {
    1383      294994 :   GEN k = (GEN)data;
    1384      294994 :   long i, s, lk = lg(k), lx = minss(lg(x), lg(y));
    1385             : 
    1386      294994 :   if (!is_vec_t(typ(x))) pari_err_TYPE("lexicographic vecsort",x);
    1387      294994 :   if (!is_vec_t(typ(y))) pari_err_TYPE("lexicographic vecsort",y);
    1388      306663 :   for (i=1; i<lk; i++)
    1389             :   {
    1390      295022 :     long c = k[i];
    1391      295022 :     if (c >= lx)
    1392          14 :       pari_err_TYPE("lexicographic vecsort, index too large", stoi(c));
    1393      295008 :     s = lexcmp(gel(x,c), gel(y,c));
    1394      295008 :     if (s) return s;
    1395             :   }
    1396       11641 :   return 0;
    1397             : }
    1398             : 
    1399             : /* return permutation sorting v[1..n], removing duplicates. Assume n > 0 */
    1400             : static GEN
    1401     1264075 : gen_sortspec_uniq(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1402             : {
    1403             :   pari_sp av;
    1404             :   long NX, nx, ny, m, ix, iy, i;
    1405             :   GEN x, y, w, W;
    1406             :   int s;
    1407     1264075 :   switch(n)
    1408             :   {
    1409       40396 :     case 1: return mkvecsmall(1);
    1410      502056 :     case 2:
    1411      502056 :       s = cmp(E,gel(v,1),gel(v,2));
    1412      502056 :       if      (s < 0) return mkvecsmall2(1,2);
    1413      257859 :       else if (s > 0) return mkvecsmall2(2,1);
    1414        3325 :       return mkvecsmall(1);
    1415      165255 :     case 3:
    1416      165255 :       s = cmp(E,gel(v,1),gel(v,2));
    1417      165255 :       if (s < 0) {
    1418       82228 :         s = cmp(E,gel(v,2),gel(v,3));
    1419       82228 :         if (s < 0) return mkvecsmall3(1,2,3);
    1420       55005 :         else if (s == 0) return mkvecsmall2(1,2);
    1421       54592 :         s = cmp(E,gel(v,1),gel(v,3));
    1422       54592 :         if      (s < 0) return mkvecsmall3(1,3,2);
    1423       27173 :         else if (s > 0) return mkvecsmall3(3,1,2);
    1424        1882 :         return mkvecsmall2(1,2);
    1425       83027 :       } else if (s > 0) {
    1426       81263 :         s = cmp(E,gel(v,1),gel(v,3));
    1427       81263 :         if (s < 0) return mkvecsmall3(2,1,3);
    1428       54600 :         else if (s == 0) return mkvecsmall2(2,1);
    1429       53284 :         s = cmp(E,gel(v,2),gel(v,3));
    1430       53284 :         if (s < 0) return mkvecsmall3(2,3,1);
    1431       27076 :         else if (s > 0) return mkvecsmall3(3,2,1);
    1432         742 :         return mkvecsmall2(2,1);
    1433             :       } else {
    1434        1764 :         s = cmp(E,gel(v,1),gel(v,3));
    1435        1764 :         if (s < 0) return mkvecsmall2(1,3);
    1436        1043 :         else if (s == 0) return mkvecsmall(1);
    1437         798 :         return mkvecsmall2(3,1);
    1438             :       }
    1439             :   }
    1440      556368 :   NX = nx = n>>1; ny = n-nx;
    1441      556368 :   av = avma;
    1442      556368 :   x = gen_sortspec_uniq(v,   nx,E,cmp); nx = lg(x)-1;
    1443      556368 :   y = gen_sortspec_uniq(v+NX,ny,E,cmp); ny = lg(y)-1;
    1444      556368 :   w = cgetg(n+1, t_VECSMALL);
    1445      556368 :   m = ix = iy = 1;
    1446     7874716 :   while (ix<=nx && iy<=ny)
    1447             :   {
    1448     7318348 :     s = cmp(E, gel(v,x[ix]), gel(v,y[iy]+NX));
    1449     7318348 :     if (s < 0)
    1450     3178358 :       w[m++] = x[ix++];
    1451     4139990 :     else if (s > 0)
    1452     3274271 :       w[m++] = y[iy++]+NX;
    1453             :     else {
    1454      865719 :       w[m++] = x[ix++];
    1455      865719 :       iy++;
    1456             :     }
    1457             :   }
    1458      984418 :   while (ix<=nx) w[m++] = x[ix++];
    1459     1055693 :   while (iy<=ny) w[m++] = y[iy++]+NX;
    1460      556368 :   set_avma(av);
    1461      556368 :   W = cgetg(m, t_VECSMALL);
    1462     8802091 :   for (i = 1; i < m; i++) W[i] = w[i];
    1463      556368 :   return W;
    1464             : }
    1465             : 
    1466             : /* return permutation sorting v[1..n]. Assume n > 0 */
    1467             : static GEN
    1468   117033432 : gen_sortspec(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1469             : {
    1470             :   long nx, ny, m, ix, iy;
    1471             :   GEN x, y, w;
    1472   117033432 :   switch(n)
    1473             :   {
    1474     2798121 :     case 1:
    1475     2798121 :       (void)cmp(E,gel(v,1),gel(v,1)); /* check for type error */
    1476     2798125 :       return mkvecsmall(1);
    1477    48083251 :     case 2:
    1478    82481828 :       return cmp(E,gel(v,1),gel(v,2)) <= 0? mkvecsmall2(1,2)
    1479    82481821 :                                           : mkvecsmall2(2,1);
    1480    22608418 :     case 3:
    1481    22608418 :       if (cmp(E,gel(v,1),gel(v,2)) <= 0) {
    1482    16223986 :         if (cmp(E,gel(v,2),gel(v,3)) <= 0) return mkvecsmall3(1,2,3);
    1483     7417357 :         return (cmp(E,gel(v,1),gel(v,3)) <= 0)? mkvecsmall3(1,3,2)
    1484     7417357 :                                               : mkvecsmall3(3,1,2);
    1485             :       } else {
    1486     6384432 :         if (cmp(E,gel(v,1),gel(v,3)) <= 0) return mkvecsmall3(2,1,3);
    1487     6528391 :         return (cmp(E,gel(v,2),gel(v,3)) <= 0)? mkvecsmall3(2,3,1)
    1488     6528392 :                                               : mkvecsmall3(3,2,1);
    1489             :       }
    1490             :   }
    1491    43543642 :   nx = n>>1; ny = n-nx;
    1492    43543642 :   w = cgetg(n+1,t_VECSMALL);
    1493    43543657 :   x = gen_sortspec(v,   nx,E,cmp);
    1494    43543646 :   y = gen_sortspec(v+nx,ny,E,cmp);
    1495    43543648 :   m = ix = iy = 1;
    1496   300097878 :   while (ix<=nx && iy<=ny)
    1497   256554231 :     if (cmp(E, gel(v,x[ix]), gel(v,y[iy]+nx))<=0)
    1498   142130479 :       w[m++] = x[ix++];
    1499             :     else
    1500   114423751 :       w[m++] = y[iy++]+nx;
    1501    66134558 :   while (ix<=nx) w[m++] = x[ix++];
    1502   111027957 :   while (iy<=ny) w[m++] = y[iy++]+nx;
    1503    43543647 :   set_avma((pari_sp)w); return w;
    1504             : }
    1505             : 
    1506             : static void
    1507    27560546 : init_sort(GEN *x, long *tx, long *lx)
    1508             : {
    1509    27560546 :   *tx = typ(*x);
    1510    27560546 :   if (*tx == t_LIST)
    1511             :   {
    1512          35 :     if (list_typ(*x)!=t_LIST_RAW) pari_err_TYPE("sort",*x);
    1513          35 :     *x = list_data(*x);
    1514          35 :     *lx = *x? lg(*x): 1;
    1515             :   } else {
    1516    27560511 :     if (!is_matvec_t(*tx) && *tx != t_VECSMALL) pari_err_TYPE("gen_sort",*x);
    1517    27560511 :     *lx = lg(*x);
    1518             :   }
    1519    27560546 : }
    1520             : 
    1521             : /* (x o y)[1..lx-1], destroy y */
    1522             : INLINE GEN
    1523     1611942 : sort_extract(GEN x, GEN y, long tx, long lx)
    1524             : {
    1525             :   long i;
    1526     1611942 :   switch(tx)
    1527             :   {
    1528           7 :     case t_VECSMALL:
    1529          35 :       for (i=1; i<lx; i++) y[i] = x[y[i]];
    1530           7 :       break;
    1531           7 :     case t_LIST:
    1532           7 :       settyp(y,t_VEC);
    1533          35 :       for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1534           7 :       return gtolist(y);
    1535     1611928 :     default:
    1536     1611928 :       settyp(y,tx);
    1537     5410422 :       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,y[i]));
    1538             :   }
    1539     1611957 :   return y;
    1540             : }
    1541             : 
    1542             : static GEN
    1543      794950 : triv_sort(long tx) { return tx == t_LIST? mklist(): cgetg(1, tx); }
    1544             : /* Sort the vector x, using cmp to compare entries. */
    1545             : GEN
    1546      151844 : gen_sort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1547             : {
    1548             :   long tx, lx;
    1549             :   GEN y;
    1550             : 
    1551      151844 :   init_sort(&x, &tx, &lx);
    1552      151844 :   if (lx==1) return triv_sort(tx);
    1553      151164 :   y = gen_sortspec_uniq(x,lx-1,E,cmp);
    1554      151164 :   return sort_extract(x, y, tx, lg(y)); /* lg(y) <= lx */
    1555             : }
    1556             : /* Sort the vector x, using cmp to compare entries. */
    1557             : GEN
    1558     2255048 : gen_sort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1559             : {
    1560             :   long tx, lx;
    1561             :   GEN y;
    1562             : 
    1563     2255048 :   init_sort(&x, &tx, &lx);
    1564     2255048 :   if (lx==1) return triv_sort(tx);
    1565     1460778 :   y = gen_sortspec(x,lx-1,E,cmp);
    1566     1460772 :   return sort_extract(x, y, tx, lx);
    1567             : }
    1568             : /* indirect sort: return the permutation that would sort x */
    1569             : GEN
    1570         175 : gen_indexsort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1571             : {
    1572             :   long tx, lx;
    1573         175 :   init_sort(&x, &tx, &lx);
    1574         175 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1575         175 :   return gen_sortspec_uniq(x,lx-1,E,cmp);
    1576             : }
    1577             : /* indirect sort: return the permutation that would sort x */
    1578             : GEN
    1579      186716 : gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1580             : {
    1581             :   long tx, lx;
    1582      186716 :   init_sort(&x, &tx, &lx);
    1583      186716 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1584      186415 :   return gen_sortspec(x,lx-1,E,cmp);
    1585             : }
    1586             : 
    1587             : /* Sort the vector x in place, using cmp to compare entries */
    1588             : void
    1589    24835373 : gen_sort_inplace(GEN x, void *E, int (*cmp)(void*,GEN,GEN), GEN *perm)
    1590             : {
    1591             :   long tx, lx, i;
    1592    24835373 :   pari_sp av = avma;
    1593             :   GEN y;
    1594             : 
    1595    24835373 :   init_sort(&x, &tx, &lx);
    1596    24835373 :   if (lx<=2)
    1597             :   {
    1598      255523 :     if (perm) *perm = lx == 1? cgetg(1, t_VECSMALL): mkvecsmall(1);
    1599      255523 :     return;
    1600             :   }
    1601    24579850 :   y = gen_sortspec(x,lx-1, E, cmp);
    1602    24579850 :   if (perm)
    1603             :   {
    1604       10633 :     GEN z = new_chunk(lx);
    1605      106274 :     for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1606      106274 :     for (i=1; i<lx; i++) gel(x,i) = gel(z,i);
    1607       10633 :     *perm = y;
    1608       10633 :     set_avma((pari_sp)y);
    1609             :   } else {
    1610   180197798 :     for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1611   180197798 :     for (i=1; i<lx; i++) gel(x,i) = gel(y,i);
    1612    24569217 :     set_avma(av);
    1613             :   }
    1614             : }
    1615             : GEN
    1616      131362 : gen_sort_shallow(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1617             : {
    1618             :   long tx, lx, i;
    1619             :   pari_sp av;
    1620             :   GEN y, z;
    1621             : 
    1622      131362 :   init_sort(&x, &tx, &lx);
    1623      131362 :   if (lx<=2) return x;
    1624       63805 :   z = cgetg(lx, tx); av = avma;
    1625       63805 :   y = gen_sortspec(x,lx-1, E, cmp);
    1626      408660 :   for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1627       63805 :   return gc_const(av, z);
    1628             : }
    1629             : 
    1630             : static int
    1631        7889 : closurecmp(void *data, GEN x, GEN y)
    1632             : {
    1633        7889 :   pari_sp av = avma;
    1634        7889 :   long s = gsigne(closure_callgen2((GEN)data, x,y));
    1635        7889 :   set_avma(av); return s;
    1636             : }
    1637             : static void
    1638         126 : check_positive_entries(GEN k)
    1639             : {
    1640         126 :   long i, l = lg(k);
    1641         287 :   for (i=1; i<l; i++)
    1642         161 :     if (k[i] <= 0) pari_err_DOMAIN("sort_function", "index", "<", gen_0, stoi(k[i]));
    1643         126 : }
    1644             : 
    1645             : typedef int (*CMP_FUN)(void*,GEN,GEN);
    1646             : /* return NULL if t_CLOSURE k is a "key" (arity 1) and not a sorting func */
    1647             : static CMP_FUN
    1648      126693 : sort_function(void **E, GEN x, GEN k)
    1649             : {
    1650      126693 :   int (*cmp)(GEN,GEN) = &lexcmp;
    1651      126693 :   long tx = typ(x);
    1652      126693 :   if (!k)
    1653             :   {
    1654      126021 :     *E = (void*)((typ(x) == t_VECSMALL)? cmp_small: cmp);
    1655      126021 :     return &cmp_nodata;
    1656             :   }
    1657         672 :   if (tx == t_VECSMALL) pari_err_TYPE("sort_function", x);
    1658         658 :   switch(typ(k))
    1659             :   {
    1660          91 :     case t_INT: k = mkvecsmall(itos(k));  break;
    1661          35 :     case t_VEC: case t_COL: k = ZV_to_zv(k); break;
    1662           0 :     case t_VECSMALL: break;
    1663         532 :     case t_CLOSURE:
    1664         532 :      if (closure_is_variadic(k))
    1665           0 :        pari_err_TYPE("sort_function, variadic cmpf",k);
    1666         532 :      *E = (void*)k;
    1667         532 :      switch(closure_arity(k))
    1668             :      {
    1669          35 :        case 1: return NULL; /* wrt key */
    1670         497 :        case 2: return &closurecmp;
    1671           0 :        default: pari_err_TYPE("sort_function, cmpf arity != 1, 2",k);
    1672             :      }
    1673           0 :     default: pari_err_TYPE("sort_function",k);
    1674             :   }
    1675         126 :   check_positive_entries(k);
    1676         126 :   *E = (void*)k; return &veccmp;
    1677             : }
    1678             : 
    1679             : #define cmp_IND 1
    1680             : #define cmp_LEX 2 /* FIXME: backward compatibility, ignored */
    1681             : #define cmp_REV 4
    1682             : #define cmp_UNIQ 8
    1683             : GEN
    1684         714 : vecsort0(GEN x, GEN k, long flag)
    1685             : {
    1686             :   void *E;
    1687         714 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
    1688             : 
    1689         707 :   if (flag < 0 || flag > (cmp_REV|cmp_LEX|cmp_IND|cmp_UNIQ))
    1690           0 :     pari_err_FLAG("vecsort");
    1691         707 :   if (!CMP)
    1692             :   { /* wrt key: precompute all values, O(n) calls instead of O(n log n) */
    1693          28 :     pari_sp av = avma;
    1694             :     GEN v, y;
    1695             :     long i, tx, lx;
    1696          28 :     init_sort(&x, &tx, &lx);
    1697          28 :     if (lx == 1) return flag&cmp_IND? cgetg(1,t_VECSMALL): triv_sort(tx);
    1698          28 :     v = cgetg(lx, t_VEC);
    1699         140 :     for (i = 1; i < lx; i++) gel(v,i) = closure_callgen1(k, gel(x,i));
    1700          28 :     y = vecsort0(v, NULL, flag | cmp_IND);
    1701          28 :     y = flag&cmp_IND? y: sort_extract(x, y, tx, lg(y));
    1702          28 :     return gerepileupto(av, y);
    1703             :   }
    1704         679 :   if (flag&cmp_UNIQ)
    1705          35 :     x = flag&cmp_IND? gen_indexsort_uniq(x, E, CMP): gen_sort_uniq(x, E, CMP);
    1706             :   else
    1707         644 :     x = flag&cmp_IND? gen_indexsort(x, E, CMP): gen_sort(x, E, CMP);
    1708         665 :   if (flag & cmp_REV)
    1709             :   { /* reverse order */
    1710          35 :     GEN y = x;
    1711          35 :     if (typ(x)==t_LIST) { y = list_data(x); if (!y) return x; }
    1712          28 :     vecreverse_inplace(y);
    1713             :   }
    1714         658 :   return x;
    1715             : }
    1716             : 
    1717             : GEN
    1718       17233 : indexsort(GEN x) { return gen_indexsort(x, (void*)&gcmp, cmp_nodata); }
    1719             : GEN
    1720           0 : indexlexsort(GEN x) { return gen_indexsort(x, (void*)&lexcmp, cmp_nodata); }
    1721             : GEN
    1722          42 : indexvecsort(GEN x, GEN k)
    1723             : {
    1724          42 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1725          42 :   return gen_indexsort(x, (void*)k, &veccmp);
    1726             : }
    1727             : 
    1728             : GEN
    1729      725653 : sort(GEN x) { return gen_sort(x, (void*)gcmp, cmp_nodata); }
    1730             : GEN
    1731           0 : lexsort(GEN x) { return gen_sort(x, (void*)lexcmp, cmp_nodata); }
    1732             : GEN
    1733        2954 : vecsort(GEN x, GEN k)
    1734             : {
    1735        2954 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1736        2954 :   return gen_sort(x, (void*)k, &veccmp);
    1737             : }
    1738             : /* adapted from gen_search; don't export: keys of T[i] should be precomputed */
    1739             : static long
    1740           7 : key_search(GEN T, GEN x, GEN code)
    1741             : {
    1742           7 :   long u = lg(T)-1, i, l, s;
    1743             : 
    1744           7 :   if (!u) return 0;
    1745           7 :   l = 1; x = closure_callgen1(code, x);
    1746             :   do
    1747             :   {
    1748          14 :     i = (l+u)>>1; s = lexcmp(x, closure_callgen1(code, gel(T,i)));
    1749          14 :     if (!s) return i;
    1750           7 :     if (s<0) u=i-1; else l=i+1;
    1751           7 :   } while (u>=l);
    1752           0 :   return 0;
    1753             : }
    1754             : long
    1755      125979 : vecsearch(GEN v, GEN x, GEN k)
    1756             : {
    1757      125979 :   pari_sp av = avma;
    1758             :   void *E;
    1759      125979 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, v, k);
    1760      125972 :   long r, tv = typ(v);
    1761      125972 :   if (tv == t_VECSMALL)
    1762          21 :     x = (GEN)itos(x);
    1763      125951 :   else if (!is_matvec_t(tv)) pari_err_TYPE("vecsearch", v);
    1764      125972 :   r = CMP? gen_search(v, x, 0, E, CMP): key_search(v, x, k);
    1765      125972 :   return gc_long(av, r);
    1766             : }
    1767             : 
    1768             : GEN
    1769         914 : ZV_indexsort(GEN L) { return gen_indexsort(L, (void*)&cmpii, &cmp_nodata); }
    1770             : GEN
    1771      783006 : ZV_sort(GEN L) { return gen_sort(L, (void*)&cmpii, &cmp_nodata); }
    1772             : GEN
    1773        8316 : ZV_sort_uniq(GEN L) { return gen_sort_uniq(L, (void*)&cmpii, &cmp_nodata); }
    1774             : void
    1775      363623 : ZV_sort_inplace(GEN L) { gen_sort_inplace(L, (void*)&cmpii, &cmp_nodata,NULL); }
    1776             : 
    1777             : GEN
    1778         539 : vec_equiv(GEN F)
    1779             : {
    1780         539 :   pari_sp av = avma;
    1781         539 :   long j, k, L = lg(F);
    1782         539 :   GEN w = cgetg(L, t_VEC);
    1783         539 :   GEN perm = gen_indexsort(F, (void*)&cmp_universal, cmp_nodata);
    1784        1407 :   for (j = k = 1; j < L;)
    1785             :   {
    1786         868 :     GEN v = cgetg(L, t_VECSMALL);
    1787         868 :     long l = 1, o = perm[j];
    1788         868 :     v[l++] = o;
    1789        3269 :     for (j++; j < L; v[l++] = perm[j++])
    1790        2730 :       if (!gequal(gel(F,o), gel(F, perm[j]))) break;
    1791         868 :     setlg(v, l); gel(w, k++) = v;
    1792             :   }
    1793         539 :   setlg(w, k); return gerepilecopy(av,w);
    1794             : }
    1795             : 
    1796             : GEN
    1797       14861 : vec_reduce(GEN v, GEN *pE)
    1798             : {
    1799       14861 :   GEN E, F, P = gen_indexsort(v, (void*)cmp_universal, cmp_nodata);
    1800             :   long i, m, l;
    1801       14861 :   F = cgetg_copy(v, &l);
    1802       14861 :   *pE = E = cgetg(l, t_VECSMALL);
    1803       36694 :   for (i = m = 1; i < l;)
    1804             :   {
    1805       21833 :     GEN u = gel(v, P[i]);
    1806             :     long k;
    1807       26362 :     for(k = i + 1; k < l; k++)
    1808       11508 :       if (cmp_universal(gel(v, P[k]), u)) break;
    1809       21833 :     E[m] = k - i; gel(F, m) = u; i = k; m++;
    1810             :   }
    1811       14861 :   setlg(F, m);
    1812       14861 :   setlg(E, m); return F;
    1813             : }
    1814             : 
    1815             : /********************************************************************/
    1816             : /**                      SEARCH IN SORTED VECTOR                   **/
    1817             : /********************************************************************/
    1818             : /* index of x in table T, 0 otherwise */
    1819             : long
    1820      916825 : tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
    1821             : {
    1822      916825 :   long l = 1, u = lg(T)-1, i, s;
    1823             : 
    1824     6791204 :   while (u>=l)
    1825             :   {
    1826     6782818 :     i = (l+u)>>1; s = cmp(x, gel(T,i));
    1827     6782818 :     if (!s) return i;
    1828     5874379 :     if (s<0) u=i-1; else l=i+1;
    1829             :   }
    1830        8386 :   return 0;
    1831             : }
    1832             : 
    1833             : /* looks if x belongs to the set T and returns the index if yes, 0 if no */
    1834             : long
    1835    20296829 : gen_search(GEN T, GEN x, long flag, void *data, int (*cmp)(void*,GEN,GEN))
    1836             : {
    1837    20296829 :   long u = lg(T)-1, i, l, s;
    1838             : 
    1839    20296829 :   if (!u) return flag? 1: 0;
    1840    20296829 :   l = 1;
    1841             :   do
    1842             :   {
    1843    87221302 :     i = (l+u)>>1; s = cmp(data, x, gel(T,i));
    1844    87221302 :     if (!s) return flag? 0: i;
    1845    67026162 :     if (s<0) u=i-1; else l=i+1;
    1846    67026162 :   } while (u>=l);
    1847      101689 :   if (!flag) return 0;
    1848       52724 :   return (s<0)? i: i+1;
    1849             : }
    1850             : 
    1851             : long
    1852      906850 : ZV_search(GEN x, GEN y) { return tablesearch(x, y, cmpii); }
    1853             : 
    1854             : long
    1855    18382066 : zv_search(GEN T, long x)
    1856             : {
    1857    18382066 :   long l = 1, u = lg(T)-1;
    1858    74046811 :   while (u>=l)
    1859             :   {
    1860    58957483 :     long i = (l+u)>>1;
    1861    58957483 :     if (x < T[i]) u = i-1;
    1862    33769638 :     else if (x > T[i]) l = i+1;
    1863     3292738 :     else return i;
    1864             :   }
    1865    15089328 :   return 0;
    1866             : }
    1867             : 
    1868             : /********************************************************************/
    1869             : /**                   COMPARISON FUNCTIONS                         **/
    1870             : /********************************************************************/
    1871             : int
    1872   422178788 : cmp_nodata(void *data, GEN x, GEN y)
    1873             : {
    1874   422178788 :   int (*cmp)(GEN,GEN)=(int (*)(GEN,GEN)) data;
    1875   422178788 :   return cmp(x,y);
    1876             : }
    1877             : 
    1878             : /* assume x and y come from the same idealprimedec call (uniformizer unique) */
    1879             : int
    1880      895194 : cmp_prime_over_p(GEN x, GEN y)
    1881             : {
    1882      895194 :   long k = pr_get_f(x) - pr_get_f(y); /* diff. between residue degree */
    1883       25013 :   return k? ((k > 0)? 1: -1)
    1884      920210 :           : ZV_cmp(pr_get_gen(x), pr_get_gen(y));
    1885             : }
    1886             : 
    1887             : int
    1888       76938 : cmp_prime_ideal(GEN x, GEN y)
    1889             : {
    1890       76938 :   int k = cmpii(pr_get_p(x), pr_get_p(y));
    1891       76938 :   return k? k: cmp_prime_over_p(x,y);
    1892             : }
    1893             : 
    1894             : /* assume x and y are t_POL in the same variable whose coeffs can be
    1895             :  * compared (used to sort polynomial factorizations) */
    1896             : int
    1897     2115607 : gen_cmp_RgX(void *data, GEN x, GEN y)
    1898             : {
    1899     2115607 :   int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
    1900     2115607 :   long i, lx = lg(x), ly = lg(y);
    1901             :   int fl;
    1902     2115607 :   if (lx > ly) return  1;
    1903     2095457 :   if (lx < ly) return -1;
    1904     5015242 :   for (i=lx-1; i>1; i--)
    1905     4773625 :     if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
    1906      241617 :   return 0;
    1907             : }
    1908             : 
    1909             : static int
    1910         511 : cmp_RgX_Rg(GEN x, GEN y)
    1911             : {
    1912         511 :   long lx = lgpol(x), ly;
    1913         511 :   if (lx > 1) return  1;
    1914           0 :   ly = gequal0(y) ? 0:1;
    1915           0 :   if (lx > ly) return  1;
    1916           0 :   if (lx < ly) return -1;
    1917           0 :   if (lx==0) return 0;
    1918           0 :   return gcmp(gel(x,2), y);
    1919             : }
    1920             : int
    1921       50584 : cmp_RgX(GEN x, GEN y)
    1922             : {
    1923       50584 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    1924       50584 :   if (typ(y) == t_POLMOD) y = gel(y,2);
    1925       50584 :   if (typ(x) == t_POL) {
    1926       32717 :     if (typ(y) != t_POL) return cmp_RgX_Rg(x, y);
    1927             :   } else {
    1928       17867 :     if (typ(y) != t_POL) return gcmp(x,y);
    1929         427 :     return - cmp_RgX_Rg(y,x);
    1930             :   }
    1931       32633 :   return gen_cmp_RgX((void*)&gcmp,x,y);
    1932             : }
    1933             : 
    1934             : int
    1935      229016 : cmp_Flx(GEN x, GEN y)
    1936             : {
    1937      229016 :   long i, lx = lg(x), ly = lg(y);
    1938      229016 :   if (lx > ly) return  1;
    1939      213376 :   if (lx < ly) return -1;
    1940      410305 :   for (i=lx-1; i>1; i--)
    1941      381925 :     if (uel(x,i) != uel(y,i)) return uel(x,i)<uel(y,i)? -1: 1;
    1942       28380 :   return 0;
    1943             : }
    1944             : /********************************************************************/
    1945             : /**                   MERGE & SORT FACTORIZATIONS                  **/
    1946             : /********************************************************************/
    1947             : /* merge fx, fy two factorizations, whose 1st column is sorted in strictly
    1948             :  * increasing order wrt cmp */
    1949             : GEN
    1950     1031747 : merge_factor(GEN fx, GEN fy, void *data, int (*cmp)(void *,GEN,GEN))
    1951             : {
    1952     1031747 :   GEN x = gel(fx,1), e = gel(fx,2), M, E;
    1953     1031747 :   GEN y = gel(fy,1), f = gel(fy,2);
    1954     1031747 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    1955             : 
    1956     1031747 :   M = cgetg(l, t_COL);
    1957     1031747 :   E = cgetg(l, t_COL);
    1958             : 
    1959     1031747 :   m = ix = iy = 1;
    1960    10761235 :   while (ix<lx && iy<ly)
    1961             :   {
    1962     9729488 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    1963     9729488 :     if (s < 0)
    1964     8760407 :     { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; }
    1965      969081 :     else if (s == 0)
    1966             :     {
    1967      189906 :       GEN z = gel(x,ix), g = addii(gel(e,ix), gel(f,iy));
    1968      189906 :       iy++; ix++; if (!signe(g)) continue;
    1969      105969 :       gel(M,m) = z; gel(E,m) = g;
    1970             :     }
    1971             :     else
    1972      779175 :     { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; }
    1973     9645551 :     m++;
    1974             :   }
    1975     5482965 :   while (ix<lx) { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; m++; }
    1976     1293894 :   while (iy<ly) { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; m++; }
    1977     1031747 :   setlg(M, m);
    1978     1031747 :   setlg(E, m); return mkmat2(M, E);
    1979             : }
    1980             : /* merge two sorted vectors, removing duplicates. Shallow */
    1981             : GEN
    1982       43063 : merge_sort_uniq(GEN x, GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    1983             : {
    1984       43063 :   long i, j, k, lx = lg(x), ly = lg(y);
    1985       43063 :   GEN z = cgetg(lx + ly - 1, typ(x));
    1986       43063 :   i = j = k = 1;
    1987       70964 :   while (i<lx && j<ly)
    1988             :   {
    1989       27901 :     int s = cmp(data, gel(x,i), gel(y,j));
    1990       27901 :     if (s < 0)
    1991       17873 :       gel(z,k++) = gel(x,i++);
    1992       10028 :     else if (s > 0)
    1993        9965 :       gel(z,k++) = gel(y,j++);
    1994             :     else
    1995          63 :     { gel(z,k++) = gel(x,i++); j++; }
    1996             :   }
    1997       78765 :   while (i<lx) gel(z,k++) = gel(x,i++);
    1998       68794 :   while (j<ly) gel(z,k++) = gel(y,j++);
    1999       43063 :   setlg(z, k); return z;
    2000             : }
    2001             : /* in case of equal keys in x,y, take the key from x */
    2002             : static GEN
    2003      364129 : ZV_union_shallow_t(GEN x, GEN y, long t)
    2004             : {
    2005      364129 :   long i, j, k, lx = lg(x), ly = lg(y);
    2006      364129 :   GEN z = cgetg(lx + ly - 1, t);
    2007      364129 :   i = j = k = 1;
    2008     1023353 :   while (i<lx && j<ly)
    2009             :   {
    2010      659224 :     int s = cmpii(gel(x,i), gel(y,j));
    2011      659224 :     if (s < 0)
    2012        9980 :       gel(z,k++) = gel(x,i++);
    2013      649244 :     else if (s > 0)
    2014      242281 :       gel(z,k++) = gel(y,j++);
    2015             :     else
    2016      406963 :     { gel(z,k++) = gel(x,i++); j++; }
    2017             :   }
    2018      390310 :   while (i < lx) gel(z,k++) = gel(x,i++);
    2019      391641 :   while (j < ly) gel(z,k++) = gel(y,j++);
    2020      364129 :   setlg(z, k); return z;
    2021             : }
    2022             : GEN
    2023          49 : ZV_union_shallow(GEN x, GEN y)
    2024          49 : { return ZV_union_shallow_t(x, y, t_VEC); }
    2025             : GEN
    2026      364080 : ZC_union_shallow(GEN x, GEN y)
    2027      364080 : { return ZV_union_shallow_t(x, y, t_COL); }
    2028             : 
    2029             : /* sort generic factorization, in place */
    2030             : GEN
    2031     3664046 : sort_factor(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    2032             : {
    2033             :   GEN a, b, A, B, w;
    2034             :   pari_sp av;
    2035             :   long n, i;
    2036             : 
    2037     3664046 :   a = gel(y,1); n = lg(a); if (n == 1) return y;
    2038     3655282 :   b = gel(y,2); av = avma;
    2039     3655282 :   A = new_chunk(n);
    2040     3655283 :   B = new_chunk(n);
    2041     3655284 :   w = gen_sortspec(a, n-1, data, cmp);
    2042    10327920 :   for (i=1; i<n; i++) { long k=w[i]; gel(A,i) = gel(a,k); gel(B,i) = gel(b,k); }
    2043    10327920 :   for (i=1; i<n; i++) { gel(a,i) = gel(A,i); gel(b,i) = gel(B,i); }
    2044     3655292 :   set_avma(av); return y;
    2045             : }
    2046             : /* sort polynomial factorization, in place */
    2047             : GEN
    2048      653402 : sort_factor_pol(GEN y,int (*cmp)(GEN,GEN))
    2049             : {
    2050      653402 :   (void)sort_factor(y,(void*)cmp, &gen_cmp_RgX);
    2051      653413 :   return y;
    2052             : }
    2053             : 
    2054             : /***********************************************************************/
    2055             : /*                                                                     */
    2056             : /*                          SET OPERATIONS                             */
    2057             : /*                                                                     */
    2058             : /***********************************************************************/
    2059             : GEN
    2060      140168 : gtoset(GEN x)
    2061             : {
    2062             :   long lx;
    2063      140168 :   if (!x) return cgetg(1, t_VEC);
    2064      140168 :   switch(typ(x))
    2065             :   {
    2066      140140 :     case t_VEC:
    2067      140140 :     case t_COL: lx = lg(x); break;
    2068          14 :     case t_LIST:
    2069          14 :       if (list_typ(x)==t_LIST_MAP) return mapdomain(x);
    2070          14 :       x = list_data(x); lx = x? lg(x): 1; break;
    2071           7 :     case t_VECSMALL: lx = lg(x); x = zv_to_ZV(x); break;
    2072           7 :     default: return mkveccopy(x);
    2073             :   }
    2074      140161 :   if (lx==1) return cgetg(1,t_VEC);
    2075      140154 :   x = gen_sort_uniq(x, (void*)&cmp_universal, cmp_nodata);
    2076      140154 :   settyp(x, t_VEC); /* it may be t_COL */
    2077      140154 :   return x;
    2078             : }
    2079             : 
    2080             : long
    2081          14 : setisset(GEN x)
    2082             : {
    2083          14 :   long i, lx = lg(x);
    2084             : 
    2085          14 :   if (typ(x) != t_VEC) return 0;
    2086          14 :   if (lx == 1) return 1;
    2087          70 :   for (i=1; i<lx-1; i++)
    2088          63 :     if (cmp_universal(gel(x,i+1), gel(x,i)) <= 0) return 0;
    2089           7 :   return 1;
    2090             : }
    2091             : 
    2092             : long
    2093         119 : setsearch(GEN T, GEN y, long flag)
    2094             : {
    2095             :   long lx;
    2096         119 :   switch(typ(T))
    2097             :   {
    2098         105 :     case t_VEC: lx = lg(T); break;
    2099           7 :     case t_LIST:
    2100           7 :     if (list_typ(T) != t_LIST_RAW) pari_err_TYPE("setsearch",T);
    2101           7 :     T = list_data(T); lx = T? lg(T): 1; break;
    2102           7 :     default: pari_err_TYPE("setsearch",T);
    2103             :       return 0; /*LCOV_EXCL_LINE*/
    2104             :   }
    2105         112 :   if (lx==1) return flag? 1: 0;
    2106         112 :   return gen_search(T,y,flag,(void*)cmp_universal,cmp_nodata);
    2107             : }
    2108             : 
    2109             : GEN
    2110          63 : setunion_i(GEN x, GEN y)
    2111          63 : { return merge_sort_uniq(x,y, (void*)cmp_universal, cmp_nodata); }
    2112             : 
    2113             : GEN
    2114           7 : setunion(GEN x, GEN y)
    2115             : {
    2116           7 :   pari_sp av = avma;
    2117           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setunion",x);
    2118           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setunion",y);
    2119           7 :   return gerepilecopy(av, setunion_i(x, y));
    2120             : }
    2121             : 
    2122             : GEN
    2123           7 : setintersect(GEN x, GEN y)
    2124             : {
    2125           7 :   long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
    2126           7 :   pari_sp av = avma;
    2127           7 :   GEN z = cgetg(lx,t_VEC);
    2128           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setintersect",x);
    2129           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setintersect",y);
    2130          70 :   while (ix < lx && iy < ly)
    2131             :   {
    2132          63 :     int c = cmp_universal(gel(x,ix), gel(y,iy));
    2133          63 :     if      (c < 0) ix++;
    2134          35 :     else if (c > 0) iy++;
    2135          21 :     else { gel(z, iz++) = gel(x,ix); ix++; iy++; }
    2136             :   }
    2137           7 :   setlg(z,iz); return gerepilecopy(av,z);
    2138             : }
    2139             : 
    2140             : GEN
    2141           7 : gen_setminus(GEN A, GEN B, int (*cmp)(GEN,GEN))
    2142             : {
    2143           7 :   pari_sp ltop = avma;
    2144           7 :   long i = 1, j = 1, k = 1, lx = lg(A), ly = lg(B);
    2145           7 :   GEN  diff = cgetg(lx,t_VEC);
    2146          84 :   while (i < lx && j < ly)
    2147          77 :     switch ( cmp(gel(A,i),gel(B,j)) )
    2148             :     {
    2149          28 :       case -1: gel(diff,k++) = gel(A,i++); break;
    2150          28 :       case 1: j++; break;
    2151          21 :       case 0: i++; break;
    2152             :     }
    2153          84 :   while (i < lx) gel(diff,k++) = gel(A,i++);
    2154           7 :   setlg(diff,k);
    2155           7 :   return gerepilecopy(ltop,diff);
    2156             : }
    2157             : 
    2158             : GEN
    2159           7 : setminus(GEN x, GEN y)
    2160             : {
    2161           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setminus",x);
    2162           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setminus",y);
    2163           7 :   return gen_setminus(x,y,cmp_universal);
    2164             : }
    2165             : 
    2166             : GEN
    2167          21 : setbinop(GEN f, GEN x, GEN y)
    2168             : {
    2169          21 :   pari_sp av = avma;
    2170          21 :   long i, j, lx, ly, k = 1;
    2171             :   GEN z;
    2172          21 :   if (typ(f) != t_CLOSURE || closure_arity(f) != 2 || closure_is_variadic(f))
    2173           7 :     pari_err_TYPE("setbinop [function needs exactly 2 arguments]",f);
    2174          14 :   lx = lg(x);
    2175          14 :   if (typ(x) != t_VEC) pari_err_TYPE("setbinop", x);
    2176          14 :   if (y == NULL) { /* assume x = y and f symmetric */
    2177           7 :     z = cgetg((((lx-1)*lx) >> 1) + 1, t_VEC);
    2178          28 :     for (i = 1; i < lx; i++)
    2179          63 :       for (j = i; j < lx; j++)
    2180          42 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(x,j));
    2181             :   } else {
    2182           7 :     ly = lg(y);
    2183           7 :     if (typ(y) != t_VEC) pari_err_TYPE("setbinop", y);
    2184           7 :     z = cgetg((lx-1)*(ly-1) + 1, t_VEC);
    2185          28 :     for (i = 1; i < lx; i++)
    2186          84 :       for (j = 1; j < ly; j++)
    2187          63 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(y,j));
    2188             :   }
    2189          14 :   return gerepileupto(av, gtoset(z));
    2190             : }

Generated by: LCOV version 1.13