Code coverage tests

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

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

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - bibli2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21072-998352a) Lines: 1051 1079 97.4 %
Date: 2017-09-26 06:25:23 Functions: 95 97 97.9 %
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          84 :   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          14 :   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         112 :   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        1428 : polhermite(long n, long v)
     180             : {
     181             :   long m;
     182             :   pari_sp av;
     183             :   GEN q,a,r;
     184             : 
     185        1428 :   if (v<0) v = 0;
     186        1428 :   if (n < 0) pari_err_DOMAIN("polhermite", "degree", "<", gen_0, stoi(n));
     187        1428 :   if (n==0) return pol_1(v);
     188             : 
     189        1421 :   q = cgetg(n+3, t_POL); r = q + n+2;
     190        1421 :   a = int2n(n);
     191        1421 :   gel(r--,0) = a;
     192        1421 :   gel(r--,0) = gen_0;
     193       40285 :   for (m=1; 2*m<= n; m++)
     194             :   {
     195       38864 :     av = avma;
     196       38864 :     a = diviuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m);
     197       38864 :     togglesign(a);
     198       38864 :     gel(r--,0) = a = gerepileuptoint(av, a);
     199       38864 :     gel(r--,0) = gen_0;
     200             :   }
     201        1421 :   q[1] = evalsigne(1) | evalvarn(v);
     202        1421 :   return q;
     203             : }
     204             : GEN
     205        1442 : polhermite_eval(long n, GEN x)
     206             : {
     207             :   long i;
     208             :   pari_sp av, av2;
     209             :   GEN x2, u, v;
     210             : 
     211        1442 :   if (!x) return polhermite(n, 0);
     212          14 :   if (gequalX(x)) return polhermite(n, varn(x));
     213          14 :   if (n==0) return gen_1;
     214          14 :   if (n==1) return gmul2n(x,1);
     215          14 :   av = avma; x2 = gmul2n(x,1); v = gen_1; u = x2;
     216          14 :   av2= avma;
     217        7035 :   for (i=1; i<n; i++)
     218             :   { /* u = H_i(x), v = H_{i-1}(x), compute t = H_{i+1}(x) */
     219             :     GEN t;
     220        7021 :     if ((i & 0xff) == 0) gerepileall(av2,2,&u, &v);
     221        7021 :     t = gsub(gmul(x2, u), gmulsg(2*i,v));
     222        7021 :     v = u; u = t;
     223             :   }
     224          14 :   return gerepileupto(av, u);
     225             : }
     226             : 
     227             : /* Legendre polynomial
     228             :  * L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1)
     229             :  * L(n) = 2^-n sum_{k=0}^{n/2} a_k x^(n-2k)
     230             :  *   where a_k = (-1)^k (2n-2k)! / k! (n-k)! (n-2k)! is an integer
     231             :  *   and a_0 = binom(2n,n), a_k / a_{k-1} = - (n-2k+1)(n-2k+2) / 2k (2n-2k+1) */
     232             : GEN
     233        2184 : pollegendre(long n, long v)
     234             : {
     235             :   long k, l;
     236             :   pari_sp av;
     237             :   GEN a, r, q;
     238             : 
     239        2184 :   if (v<0) v = 0;
     240             :   /* pollegendre(-n) = pollegendre(n-1) */
     241        2184 :   if (n < 0) n = -n-1;
     242        2184 :   if (n==0) return pol_1(v);
     243        2142 :   if (n==1) return pol_x(v);
     244             : 
     245        2100 :   av = avma;
     246        2100 :   q = cgetg(n+3, t_POL); r = q + n+2;
     247        2100 :   gel(r--,0) = a = binomialuu(n<<1,n);
     248        2100 :   gel(r--,0) = gen_0;
     249       32354 :   for (k=1,l=n; l>1; k++,l-=2)
     250             :   { /* l = n-2*k+2 */
     251       30254 :     av = avma;
     252       30254 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     253       30254 :     togglesign(a); a = gerepileuptoint(av, a);
     254       30254 :     gel(r--,0) = a;
     255       30254 :     gel(r--,0) = gen_0;
     256             :   }
     257        2100 :   q[1] = evalsigne(1) | evalvarn(v);
     258        2100 :   return gerepileupto(av, gmul2n(q,-n));
     259             : }
     260             : 
     261             : GEN
     262        2163 : pollegendre_eval(long n, GEN x)
     263             : {
     264             :   long i;
     265             :   pari_sp av;
     266             :   GEN u, v;
     267             : 
     268        2163 :   if (!x) return pollegendre(n, 0);
     269          14 :   if (gequalX(x)) return pollegendre(n, varn(x));
     270             :   /* pollegendre(-n) = pollegendre(n-1) */
     271          14 :   if (n < 0) n = -n-1;
     272          14 :   if (n==0) return gen_1;
     273          14 :   if (n==1) return gcopy(x);
     274          14 :   av = avma; v = gen_1; u = x;
     275        7035 :   for (i=1; i<n; i++)
     276             :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     277             :     GEN t;
     278        7021 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     279        7021 :     t = gdivgs(gsub(gmul(gmulsg(2*i+1,x), u), gmulsg(i,v)), i+1);
     280        7021 :     v = u; u = t;
     281             :   }
     282          14 :   return gerepileupto(av, u);
     283             : }
     284             : 
     285             : /* polcyclo(p) = X^(p-1) + ... + 1 */
     286             : static GEN
     287      212628 : polcyclo_prime(long p, long v)
     288             : {
     289      212628 :   GEN T = cgetg(p+2, t_POL);
     290             :   long i;
     291      212628 :   T[1] = evalsigne(1) | evalvarn(v);
     292      212628 :   for (i = 2; i < p+2; i++) gel(T,i) = gen_1;
     293      212628 :   return T;
     294             : }
     295             : 
     296             : /* cyclotomic polynomial */
     297             : GEN
     298      227571 : polcyclo(long n, long v)
     299             : {
     300             :   long s, q, i, l;
     301      227571 :   pari_sp av=avma;
     302             :   GEN T, P;
     303             : 
     304      227571 :   if (v<0) v = 0;
     305      227571 :   if (n < 3)
     306       14943 :     switch(n)
     307             :     {
     308       10724 :       case 1: return deg1pol_shallow(gen_1, gen_m1, v);
     309        4219 :       case 2: return deg1pol_shallow(gen_1, gen_1, v);
     310           0 :       default: pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     311             :     }
     312      212628 :   P = gel(factoru(n), 1); l = lg(P);
     313      212628 :   s = P[1]; T = polcyclo_prime(s, v);
     314      238185 :   for (i = 2; i < l; i++)
     315             :   { /* Phi_{np}(X) = Phi_n(X^p) / Phi_n(X) */
     316       25557 :     s *= P[i];
     317       25557 :     T = RgX_div(RgX_inflate(T, P[i]), T);
     318             :   }
     319             :   /* s = squarefree part of n */
     320      212628 :   q = n / s;
     321      212628 :   if (q == 1) return gerepileupto(av, T);
     322      101899 :   return gerepilecopy(av, RgX_inflate(T,q));
     323             : }
     324             : 
     325             : /* cyclotomic polynomial */
     326             : GEN
     327       17199 : polcyclo_eval(long n, GEN x)
     328             : {
     329       17199 :   pari_sp av= avma;
     330             :   GEN P, md, xd, yneg, ypos;
     331             :   long l, s, i, j, q, tx;
     332       17199 :   long root_of_1 = 0;
     333             : 
     334       17199 :   if (!x) return polcyclo(n, 0);
     335       15736 :   tx = typ(x);
     336       15736 :   if (gequalX(x)) return polcyclo(n, varn(x));
     337       15197 :   if (n <= 0) pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     338       15197 :   if (n == 1) return gsubgs(x, 1);
     339       15197 :   if (tx == t_INT && !signe(x)) return gen_1;
     340       15197 :   while ((n & 3) == 0) { n >>= 1; x = gsqr(x); } /* Phi_4n(x) = Phi_2n(x^2) */
     341             :   /* n not divisible by 4 */
     342       15197 :   if (n == 2) return gerepileupto(av, gaddgs(x,1));
     343        1050 :   if (!odd(n)) { n >>= 1; x = gneg(x); } /* Phi_2n(x) = Phi_n(-x) for n>1 odd */
     344             :   /* n odd > 2.  s largest squarefree divisor of n */
     345        1050 :   P = gel(factoru(n), 1); s = zv_prod(P);
     346             :   /* replace n by largest squarefree divisor */
     347        1050 :   q = n/s; if (q != 1) { x = gpowgs(x, q); n = s; }
     348        1050 :   l = lg(P)-1;
     349             :   /* n squarefree odd > 2, l distinct prime divisors. Now handle x = 1 or -1 */
     350        1050 :   if (tx == t_INT) { /* shortcut */
     351         980 :     if (is_pm1(x))
     352             :     {
     353          56 :       avma = av;
     354          56 :       if (signe(x) > 0 && l == 1) return utoipos(P[1]);
     355          35 :       return gen_1;
     356             :     }
     357             :   } else {
     358          70 :     if (gequal1(x))
     359             :     { /* n is prime, return n; multiply by x to keep the type */
     360          14 :       if (l == 1) return gerepileupto(av, gmulgs(x,n));
     361           7 :       return gerepilecopy(av, x); /* else 1 */
     362             :     }
     363          56 :     if (gequalm1(x)) return gerepileupto(av, gneg(x)); /* -1 */
     364             :   }
     365             :   /* Heuristic: evaluation will probably not improve things */
     366         973 :   if (tx == t_POL || tx == t_MAT || lg(x) > n)
     367          17 :     return gerepileupto(av, poleval(polcyclo(n,0), x));
     368             : 
     369         956 :   xd = cgetg((1L<<l) + 1, t_VEC); /* the x^d, where d | n */
     370         956 :   md = cgetg((1L<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
     371         956 :   gel(xd, 1) = x;
     372         956 :   md[1] = 1;
     373             :   /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d).
     374             :    * If x has exact order D, n = Dq, then the result is 0 if q = 1. Otherwise
     375             :    * the factors with x^d-1, D|d are omitted and we multiply at the end by
     376             :    *   prod_{d | q} d^mu(q/d) = q if prime, 1 otherwise */
     377             :   /* We store the factors with mu(d)= 1 (resp.-1) in ypos (resp yneg).
     378             :    * At the end we return ypos/yneg if mu(n)=1 and yneg/ypos if mu(n)=-1 */
     379         956 :   ypos = gsubgs(x,1);
     380         956 :   yneg = gen_1;
     381        1996 :   for (i = 1; i <= l; i++)
     382             :   {
     383        1040 :     long ti = 1L<<(i-1), p = P[i];
     384        2178 :     for (j = 1; j <= ti; j++) {
     385        1138 :       GEN X = gpowgs(gel(xd,j), p), t = gsubgs(X,1);
     386        1138 :       gel(xd,ti+j) = X;
     387        1138 :       md[ti+j] = -md[j];
     388        1138 :       if (gequal0(t))
     389             :       { /* x^d = 1; root_of_1 := the smallest index ti+j such that X == 1
     390             :         * (whose bits code d: bit i-1 is set iff P[i] | d). If no such index
     391             :         * exists, then root_of_1 remains 0. Do not multiply with X-1 if X = 1,
     392             :         * we handle these factors at the end */
     393          28 :         if (!root_of_1) root_of_1 = ti+j;
     394             :       }
     395             :       else
     396             :       {
     397        1110 :         if (md[ti+j] == 1) ypos = gmul(ypos, t);
     398        1033 :         else               yneg = gmul(yneg, t);
     399             :       }
     400             :     }
     401             :   }
     402         956 :   ypos = odd(l)? gdiv(yneg,ypos): gdiv(ypos,yneg);
     403         956 :   if (root_of_1)
     404             :   {
     405          21 :     GEN X = gel(xd,(1<<l)); /* = x^n = 1 */
     406          21 :     long bitmask_q = (1<<l) - root_of_1;
     407             :     /* bitmask_q encodes q = n/d: bit (i-1) is 1 iff P[i] | q */
     408             : 
     409             :     /* x is a root of unity.  If bitmask_q = 0, then x was a primitive n-th
     410             :      * root of 1 and the result is zero. Return X - 1 to preserve type. */
     411          21 :     if (!bitmask_q) return gerepileupto(av, gsubgs(X, 1));
     412             :     /* x is a primitive d-th root of unity, where d|n and d<n: we
     413             :      * must multiply ypos by if(isprime(n/d), n/d, 1) */
     414           7 :     ypos = gmul(ypos, X); /* multiply by X = 1 to preserve type */
     415             :     /* If bitmask_q = 1<<(i-1) for some i <= l, then q == P[i] and we multiply
     416             :      * by P[i]; otherwise q is composite and nothing more needs to be done */
     417           7 :     if (!(bitmask_q & (bitmask_q-1))) /* detects power of 2, since bitmask!=0 */
     418             :     {
     419           7 :       i = vals(bitmask_q)+1; /* q = P[i] */
     420           7 :       ypos = gmulgs(ypos, P[i]);
     421             :     }
     422             :   }
     423         942 :   return gerepileupto(av, ypos);
     424             : }
     425             : /********************************************************************/
     426             : /**                                                                **/
     427             : /**                  HILBERT & PASCAL MATRICES                     **/
     428             : /**                                                                **/
     429             : /********************************************************************/
     430             : GEN
     431         133 : mathilbert(long n) /* Hilbert matrix of order n */
     432             : {
     433             :   long i,j;
     434             :   GEN p;
     435             : 
     436         133 :   if (n < 0) pari_err_DOMAIN("mathilbert", "dimension", "<", gen_0, stoi(n));
     437         133 :   p = cgetg(n+1,t_MAT);
     438        1120 :   for (j=1; j<=n; j++)
     439             :   {
     440         987 :     gel(p,j) = cgetg(n+1,t_COL);
     441       16583 :     for (i=1+(j==1); i<=n; i++)
     442       15596 :       gcoeff(p,i,j) = mkfrac(gen_1, utoipos(i+j-1));
     443             :   }
     444         133 :   if (n) gcoeff(p,1,1) = gen_1;
     445         133 :   return p;
     446             : }
     447             : 
     448             : /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
     449             : GEN
     450        1148 : matqpascal(long n, GEN q)
     451             : {
     452             :   long i, j, I;
     453        1148 :   pari_sp av = avma;
     454        1148 :   GEN m, qpow = NULL; /* gcc -Wall */
     455             : 
     456        1148 :   if (n < -1)  pari_err_DOMAIN("matpascal", "n", "<", gen_m1, stoi(n));
     457        1148 :   n++; m = cgetg(n+1,t_MAT);
     458        1148 :   for (j=1; j<=n; j++) gel(m,j) = cgetg(n+1,t_COL);
     459        1148 :   if (q)
     460             :   {
     461          42 :     I = (n+1)/2;
     462          42 :     if (I > 1) { qpow = new_chunk(I+1); gel(qpow,2)=q; }
     463          42 :     for (j=3; j<=I; j++) gel(qpow,j) = gmul(q, gel(qpow,j-1));
     464             :   }
     465       24612 :   for (i=1; i<=n; i++)
     466             :   {
     467       23464 :     I = (i+1)/2; gcoeff(m,i,1)= gen_1;
     468       23464 :     if (q)
     469             :     {
     470         483 :       for (j=2; j<=I; j++)
     471         476 :         gcoeff(m,i,j) = gadd(gmul(gel(qpow,j),gcoeff(m,i-1,j)),
     472         238 :                              gcoeff(m,i-1,j-1));
     473             :     }
     474             :     else
     475             :     {
     476      971957 :       for (j=2; j<=I; j++)
     477      948738 :         gcoeff(m,i,j) = addii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
     478             :     }
     479       23464 :     for (   ; j<=i; j++) gcoeff(m,i,j) = gcoeff(m,i,i+1-j);
     480       23464 :     for (   ; j<=n; j++) gcoeff(m,i,j) = gen_0;
     481             :   }
     482        1148 :   return gerepilecopy(av, m);
     483             : }
     484             : 
     485             : /******************************************************************/
     486             : /**                                                              **/
     487             : /**                       PRECISION CHANGES                      **/
     488             : /**                                                              **/
     489             : /******************************************************************/
     490             : 
     491             : GEN
     492         182 : gprec(GEN x, long l)
     493             : {
     494             :   long lx, i;
     495             :   GEN y;
     496             : 
     497         182 :   if (l <= 0) pari_err_DOMAIN("gprec", "precision", "<=", gen_0, stoi(l));
     498         182 :   switch(typ(x))
     499             :   {
     500             :     case t_REAL:
     501          21 :       return rtor(x, ndec2prec(l));
     502             :     case t_COMPLEX:
     503           7 :       y = cgetg(3, t_COMPLEX);
     504           7 :       gel(y,1) = gprec(gel(x,1),l);
     505           7 :       gel(y,2) = gprec(gel(x,2),l);
     506           7 :       break;
     507             :     case t_POL: case t_SER:
     508          28 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     509          28 :       for (i=2; i<lx; i++) gel(y,i) = gprec(gel(x,i),l);
     510          28 :       break;
     511             :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     512          14 :       y = cgetg_copy(x, &lx);
     513          14 :       for (i=1; i<lx; i++) gel(y,i) = gprec(gel(x,i),l);
     514          14 :       break;
     515         112 :     default: y = gcopy(x);
     516             :   }
     517         161 :   return y;
     518             : }
     519             : 
     520             : /* internal: precision given in word length (including codewords) */
     521             : GEN
     522     2407201 : gprec_w(GEN x, long pr)
     523             : {
     524             :   long lx, i;
     525             :   GEN y;
     526             : 
     527     2407201 :   switch(typ(x))
     528             :   {
     529             :     case t_REAL:
     530     1965873 :       if (signe(x)) return rtor(x,pr);
     531       15630 :       i = -prec2nbits(pr);
     532       15630 :       if (i < expo(x)) return real_0_bit(i);
     533       15139 :       y = cgetr(2); y[1] = x[1]; return y;
     534             :     case t_COMPLEX:
     535      218740 :       y = cgetg(3, t_COMPLEX);
     536      218740 :       gel(y,1) = gprec_w(gel(x,1),pr);
     537      218740 :       gel(y,2) = gprec_w(gel(x,2),pr);
     538      218740 :       break;
     539             :    case t_POL: case t_SER:
     540       59220 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     541       59220 :       for (i=2; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     542       59220 :       break;
     543             :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     544       50498 :       y = cgetg_copy(x, &lx);
     545       50498 :       for (i=1; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     546       50498 :       break;
     547      112870 :     default: return x;
     548             :   }
     549      328458 :   return y;
     550             : }
     551             : 
     552             : /* internal: precision given in word length (including codewords), truncate
     553             :  * mantissa to precision 'pr' but never _increase_ it */
     554             : GEN
     555      351882 : gprec_wtrunc(GEN x, long pr)
     556             : {
     557             :   long lx, i;
     558             :   GEN y;
     559             : 
     560      351882 :   switch(typ(x))
     561             :   {
     562             :     case t_REAL:
     563      320956 :       return (signe(x) && realprec(x) > pr)? rtor(x,pr): x;
     564             :     case t_COMPLEX:
     565       30079 :       y = cgetg(3, t_COMPLEX);
     566       30079 :       gel(y,1) = gprec_wtrunc(gel(x,1),pr);
     567       30079 :       gel(y,2) = gprec_wtrunc(gel(x,2),pr);
     568       30079 :       break;
     569             :     case t_POL:
     570             :     case t_SER:
     571         301 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     572         301 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     573         301 :       break;
     574             :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     575         434 :       y = cgetg_copy(x, &lx);
     576         434 :       for (i=1; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     577         434 :       break;
     578         112 :     default: return x;
     579             :   }
     580       30814 :   return y;
     581             : }
     582             : 
     583             : /********************************************************************/
     584             : /**                                                                **/
     585             : /**                      SERIES TRANSFORMS                         **/
     586             : /**                                                                **/
     587             : /********************************************************************/
     588             : /**                  LAPLACE TRANSFORM (OF A SERIES)               **/
     589             : /********************************************************************/
     590             : static GEN
     591           7 : serlaplace(GEN x)
     592             : {
     593           7 :   long i, l = lg(x), e = valp(x);
     594           7 :   GEN t, y = cgetg(l,t_SER);
     595           7 :   if (e < 0) pari_err_DOMAIN("laplace","valuation","<",gen_0,stoi(e));
     596           7 :   t = mpfact(e); y[1] = x[1];
     597         119 :   for (i=2; i<l; i++)
     598             :   {
     599         112 :     gel(y,i) = gmul(t, gel(x,i));
     600         112 :     e++; t = mului(e,t);
     601             :   }
     602           7 :   return y;
     603             : }
     604             : static GEN
     605          14 : pollaplace(GEN x)
     606             : {
     607          14 :   long i, e = 0, l = lg(x);
     608          14 :   GEN t = gen_1, y = cgetg(l,t_POL);
     609          14 :   y[1] = x[1];
     610          63 :   for (i=2; i<l; i++)
     611             :   {
     612          49 :     gel(y,i) = gmul(t, gel(x,i));
     613          49 :     e++; t = mului(e,t);
     614             :   }
     615          14 :   return y;
     616             : }
     617             : GEN
     618          21 : laplace(GEN x)
     619             : {
     620          21 :   pari_sp av = avma;
     621          21 :   switch(typ(x))
     622             :   {
     623          14 :     case t_POL: x = pollaplace(x); break;
     624           7 :     case t_SER: x = serlaplace(x); break;
     625           0 :     default: pari_err_TYPE("laplace",x);
     626             :   }
     627          21 :   return gerepilecopy(av, x);
     628             : }
     629             : 
     630             : /********************************************************************/
     631             : /**              CONVOLUTION PRODUCT (OF TWO SERIES)               **/
     632             : /********************************************************************/
     633             : GEN
     634           7 : convol(GEN x, GEN y)
     635             : {
     636           7 :   long j, lx, ly, ex, ey, vx = varn(x);
     637             :   GEN z;
     638             : 
     639           7 :   if (typ(x) != t_SER) pari_err_TYPE("convol",x);
     640           7 :   if (typ(y) != t_SER) pari_err_TYPE("convol",y);
     641           7 :   if (varn(y) != vx) pari_err_VAR("convol", x,y);
     642           7 :   ex = valp(x);
     643           7 :   ey = valp(y);
     644           7 :   if (ser_isexactzero(x))
     645           0 :     return scalarser(gadd(RgX_get_0(x), RgX_get_0(y)), varn(x), maxss(ex,ey));
     646           7 :   lx = lg(x) + ex; x -= ex;
     647           7 :   ly = lg(y) + ey; y -= ey;
     648             :   /* inputs shifted: x[i] and y[i] now correspond to monomials of same degree */
     649           7 :   if (ly < lx) lx = ly; /* min length */
     650           7 :   if (ex < ey) ex = ey; /* max valuation */
     651           7 :   if (lx - ex < 3) return zeroser(vx, lx-2);
     652             : 
     653           7 :   z = cgetg(lx - ex, t_SER);
     654           7 :   z[1] = evalvalp(ex) | evalvarn(vx);
     655           7 :   for (j = ex+2; j<lx; j++) gel(z,j-ex) = gmul(gel(x,j),gel(y,j));
     656           7 :   return normalize(z);
     657             : }
     658             : 
     659             : /***********************************************************************/
     660             : /*               OPERATIONS ON DIRICHLET SERIES: *, /                  */
     661             : /* (+, -, scalar multiplication are done on the corresponding vectors) */
     662             : /***********************************************************************/
     663             : static long
     664        1652 : dirval(GEN x)
     665             : {
     666        1652 :   long i = 1, lx = lg(x);
     667        1652 :   while (i < lx && gequal0(gel(x,i))) i++;
     668        1652 :   return i;
     669             : }
     670             : 
     671             : GEN
     672         189 : dirmul(GEN x, GEN y)
     673             : {
     674         189 :   pari_sp av = avma, av2;
     675             :   long nx, ny, nz, dx, dy, i, j, k;
     676             :   GEN z;
     677             : 
     678         189 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirmul",x);
     679         189 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirmul",y);
     680         189 :   dx = dirval(x); nx = lg(x)-1;
     681         189 :   dy = dirval(y); ny = lg(y)-1;
     682         189 :   if (ny-dy < nx-dx) { swap(x,y); lswap(nx,ny); lswap(dx,dy); }
     683         189 :   nz = minss(nx*dy,ny*dx);
     684         189 :   y = RgV_kill0(y);
     685         189 :   av2 = avma;
     686         189 :   z = zerovec(nz);
     687        8869 :   for (j=dx; j<=nx; j++)
     688             :   {
     689        8680 :     GEN c = gel(x,j);
     690        8680 :     if (gequal0(c)) continue;
     691        5460 :     if (gequal1(c))
     692             :     {
     693       21518 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     694       18431 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gel(y,k));
     695             :     }
     696        2373 :     else if (gequalm1(c))
     697             :     {
     698        4669 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     699        3563 :         if (gel(y,k)) gel(z,i) = gsub(gel(z,i),gel(y,k));
     700             :     }
     701             :     else
     702             :     {
     703        5565 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     704        4298 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gmul(c,gel(y,k)));
     705             :     }
     706        5460 :     if (gc_needed(av2,3))
     707             :     {
     708           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirmul, %ld/%ld",j,nx);
     709           0 :       z = gerepilecopy(av2,z);
     710             :     }
     711             :   }
     712         189 :   return gerepilecopy(av,z);
     713             : }
     714             : 
     715             : GEN
     716         637 : dirdiv(GEN x, GEN y)
     717             : {
     718         637 :   pari_sp av = avma, av2;
     719             :   long nx,ny,nz, dx,dy, i,j,k;
     720             :   GEN p1;
     721             : 
     722         637 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirdiv",x);
     723         637 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirdiv",y);
     724         637 :   dx = dirval(x); nx = lg(x)-1;
     725         637 :   dy = dirval(y); ny = lg(y)-1;
     726         637 :   if (dy != 1 || !ny) pari_err_INV("dirdiv",y);
     727         637 :   nz = minss(nx,ny*dx);
     728         637 :   p1 = gel(y,1);
     729         637 :   if (gequal1(p1)) p1 = NULL; else y = gdiv(y,p1);
     730         637 :   y = RgV_kill0(y);
     731         637 :   av2 = avma;
     732         637 :   x = p1 ? gdiv(x,p1): leafcopy(x);
     733         637 :   for (j=1; j<dx; j++) gel(x,j) = gen_0;
     734         637 :   setlg(x,nz+1);
     735      305606 :   for (j=dx; j<=nz; j++)
     736             :   {
     737      304969 :     GEN c = gel(x,j);
     738      304969 :     if (gequal0(c)) continue;
     739      100359 :     if (gequal1(c))
     740             :     {
     741      618845 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     742      577079 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gel(y,k));
     743             :     }
     744       58593 :     else if (gequalm1(c))
     745             :     {
     746      506737 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     747      462945 :         if (gel(y,k)) gel(x,i) = gadd(gel(x,i),gel(y,k));
     748             :     }
     749             :     else
     750             :     {
     751       72832 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     752       58031 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gmul(c,gel(y,k)));
     753             :     }
     754      100359 :     if (gc_needed(av2,3))
     755             :     {
     756           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirdiv, %ld/%ld",j,nz);
     757           0 :       x = gerepilecopy(av2,x);
     758             :     }
     759             :   }
     760         637 :   return gerepilecopy(av,x);
     761             : }
     762             : 
     763             : /*******************************************************************/
     764             : /**                                                               **/
     765             : /**                       COMBINATORICS                           **/
     766             : /**                                                               **/
     767             : /*******************************************************************/
     768             : /**                      BINOMIAL COEFFICIENTS                    **/
     769             : /*******************************************************************/
     770             : GEN
     771       75817 : binomialuu(ulong n, ulong k)
     772             : {
     773       75817 :   pari_sp ltop = avma;
     774             :   GEN z;
     775       75817 :   if (k > n) return gen_0;
     776       75810 :   k = minuu(k,n-k);
     777       75810 :   if (!k) return gen_1;
     778       65121 :   if (k == 1) return utoipos(n);
     779       58877 :   z = diviiexact(mulu_interval(n-k+1, n), mulu_interval(2UL, k));
     780       58877 :   return gerepileuptoint(ltop,z);
     781             : }
     782             : 
     783             : GEN
     784       99981 : binomial(GEN n, long k)
     785             : {
     786             :   long i, prec;
     787             :   pari_sp av;
     788             :   GEN y;
     789             : 
     790       99981 :   if (k <= 1)
     791             :   {
     792       59283 :     if (is_noncalc_t(typ(n))) pari_err_TYPE("binomial",n);
     793       59283 :     if (k < 0) return gen_0;
     794       59283 :     if (k == 0) return gen_1;
     795       25942 :     return gcopy(n);
     796             :   }
     797       40698 :   av = avma;
     798       40698 :   if (typ(n) == t_INT)
     799             :   {
     800       40642 :     if (signe(n) > 0)
     801             :     {
     802       40635 :       GEN z = subiu(n,k);
     803       40635 :       if (cmpis(z,k) < 0)
     804             :       {
     805         910 :         k = itos(z); avma = av;
     806         910 :         if (k <= 1)
     807             :         {
     808         371 :           if (k < 0) return gen_0;
     809         371 :           if (k == 0) return gen_1;
     810         343 :           return icopy(n);
     811             :         }
     812             :       }
     813             :     }
     814             :     /* k > 1 */
     815       40271 :     if (lgefint(n) == 3 && signe(n) > 0)
     816             :     {
     817       40257 :       y = binomialuu(itou(n),(ulong)k);
     818       40257 :       return gerepileupto(av, y);
     819             :     }
     820             :     else
     821             :     {
     822          14 :       y = cgetg(k+1,t_VEC);
     823          14 :       for (i=1; i<=k; i++) gel(y,i) = subiu(n,i-1);
     824          14 :       y = ZV_prod(y);
     825             :     }
     826          14 :     y = diviiexact(y, mpfact(k));
     827          14 :     return gerepileuptoint(av, y);
     828             :   }
     829             : 
     830          56 :   prec = precision(n);
     831          56 :   if (prec && k > 200 + 0.8*prec2nbits(prec)) {
     832           7 :     GEN A = mpfactr(k, prec), B = ggamma(gsubgs(n,k-1), prec);
     833           7 :     return gerepileupto(av, gdiv(ggamma(gaddgs(n,1), prec), gmul(A,B)));
     834             :   }
     835             : 
     836          49 :   y = cgetg(k+1,t_VEC);
     837          49 :   for (i=1; i<=k; i++) gel(y,i) = gsubgs(n,i-1);
     838          49 :   return gerepileupto(av, gdiv(RgV_prod(y), mpfact(k)));
     839             : }
     840             : 
     841             : GEN
     842         210 : binomial0(GEN x, GEN k)
     843             : {
     844         210 :   if (!k)
     845             :   {
     846          21 :     if (typ(x) != t_INT || signe(x) < 0) pari_err_TYPE("binomial", x);
     847           7 :     return vecbinomial(itos(x));
     848             :   }
     849         189 :   if (typ(k) != t_INT) pari_err_TYPE("binomial", k);
     850         182 :   return binomial(x, itos(k));
     851             : }
     852             : 
     853             : /* Assume n >= 0, return bin, bin[k+1] = binomial(n, k) */
     854             : GEN
     855       18125 : vecbinomial(long n)
     856             : {
     857             :   long d, k;
     858             :   GEN C;
     859       18125 :   if (!n) return mkvec(gen_1);
     860       17838 :   C = cgetg(n+2, t_VEC) + 1; /* C[k] = binomial(n, k) */
     861       17838 :   gel(C,0) = gen_1;
     862       17838 :   gel(C,1) = utoipos(n); d = (n + 1) >> 1;
     863      105026 :   for (k=2; k <= d; k++)
     864             :   {
     865       87188 :     pari_sp av = avma;
     866       87188 :     gel(C,k) = gerepileuptoint(av, diviuexact(mului(n-k+1, gel(C,k-1)), k));
     867             :   }
     868       17838 :   for (   ; k <= n; k++) gel(C,k) = gel(C,n-k);
     869       17838 :   return C - 1;
     870             : }
     871             : 
     872             : /********************************************************************/
     873             : /**                  STIRLING NUMBERS                              **/
     874             : /********************************************************************/
     875             : /* Stirling number of the 2nd kind. The number of ways of partitioning
     876             :    a set of n elements into m non-empty subsets. */
     877             : GEN
     878        1694 : stirling2(ulong n, ulong m)
     879             : {
     880        1694 :   pari_sp av = avma;
     881             :   GEN s, bmk;
     882             :   ulong k;
     883        1694 :   if (n==0) return (m == 0)? gen_1: gen_0;
     884        1694 :   if (m > n || m == 0) return gen_0;
     885        1694 :   if (m==n) return gen_1;
     886             :   /* k = 0 */
     887        1694 :   bmk = gen_1; s  = powuu(m, n);
     888       20314 :   for (k = 1; k <= ((m-1)>>1); ++k)
     889             :   { /* bmk = binomial(m, k) */
     890             :     GEN c, kn, mkn;
     891       18620 :     bmk = diviuexact(mului(m-k+1, bmk), k);
     892       18620 :     kn  = powuu(k, n); mkn = powuu(m-k, n);
     893       18620 :     c = odd(m)? subii(mkn,kn): addii(mkn,kn);
     894       18620 :     c = mulii(bmk, c);
     895       18620 :     s = odd(k)? subii(s, c): addii(s, c);
     896       18620 :     if (gc_needed(av,2))
     897             :     {
     898           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"stirling2");
     899           0 :       gerepileall(av, 2, &s, &bmk);
     900             :     }
     901             :   }
     902             :   /* k = m/2 */
     903        1694 :   if (!odd(m))
     904             :   {
     905             :     GEN c;
     906         805 :     bmk = diviuexact(mului(k+1, bmk), k);
     907         805 :     c = mulii(bmk, powuu(k,n));
     908         805 :     s = odd(k)? subii(s, c): addii(s, c);
     909             :   }
     910        1694 :   return gerepileuptoint(av, diviiexact(s, mpfact(m)));
     911             : }
     912             : 
     913             : /* Stirling number of the first kind. Up to the sign, the number of
     914             :    permutations of n symbols which have exactly m cycles. */
     915             : GEN
     916         154 : stirling1(ulong n, ulong m)
     917             : {
     918         154 :   pari_sp ltop=avma;
     919             :   ulong k;
     920             :   GEN s, t;
     921         154 :   if (n < m) return gen_0;
     922         154 :   else if (n==m) return gen_1;
     923             :   /* t = binomial(n-1+k, m-1) * binomial(2n-m, n-m-k) */
     924             :   /* k = n-m > 0 */
     925         154 :   t = binomialuu(2*n-m-1, m-1);
     926         154 :   s = mulii(t, stirling2(2*(n-m), n-m));
     927         154 :   if (odd(n-m)) togglesign(s);
     928        1547 :   for (k = n-m-1; k > 0; --k)
     929             :   {
     930             :     GEN c;
     931        1393 :     t = diviuuexact(muluui(n-m+k+1, n+k+1, t), n+k, n-m-k);
     932        1393 :     c = mulii(t, stirling2(n-m+k, k));
     933        1393 :     s = odd(k)? subii(s, c): addii(s, c);
     934        1393 :     if ((k & 0x1f) == 0) {
     935          21 :       t = gerepileuptoint(ltop, t);
     936          21 :       s = gerepileuptoint(avma, s);
     937             :     }
     938             :   }
     939         154 :   return gerepileuptoint(ltop, s);
     940             : }
     941             : 
     942             : GEN
     943         301 : stirling(long n, long m, long flag)
     944             : {
     945         301 :   if (n < 0) pari_err_DOMAIN("stirling", "n", "<", gen_0, stoi(n));
     946         301 :   if (m < 0) pari_err_DOMAIN("stirling", "m", "<", gen_0, stoi(m));
     947         301 :   switch (flag)
     948             :   {
     949         154 :     case 1: return stirling1((ulong)n,(ulong)m);
     950         147 :     case 2: return stirling2((ulong)n,(ulong)m);
     951           0 :     default: pari_err_FLAG("stirling");
     952             :   }
     953             :   return NULL; /*LCOV_EXCL_LINE*/
     954             : }
     955             : 
     956             : /***********************************************************************/
     957             : /**                          PERMUTATIONS                             **/
     958             : /***********************************************************************/
     959             : GEN
     960        5915 : Z_to_perm(long n, GEN x)
     961             : {
     962             :   pari_sp av;
     963             :   ulong i, r;
     964        5915 :   GEN v = cgetg(n+1, t_VECSMALL);
     965        5915 :   if (n==0) return v;
     966        5908 :   uel(v,n) = 1; av = avma;
     967        5908 :   if (signe(x) <= 0) x = modii(x, mpfact(n));
     968       27146 :   for (r=n-1; r>=1; r--)
     969             :   {
     970             :     ulong a;
     971       21238 :     x = diviu_rem(x, n+1-r, &a);
     972       71687 :     for (i=r+1; i<=(ulong)n; i++)
     973       50449 :       if (uel(v,i) > a) uel(v,i)++;
     974       21238 :     uel(v,r) = a+1;
     975             :   }
     976        5908 :   avma = av; return v;
     977             : }
     978             : GEN
     979        5915 : numtoperm(long n, GEN x)
     980             : {
     981             :   long i;
     982             :   GEN v;
     983             : 
     984        5915 :   if (n < 0) pari_err_DOMAIN("numtoperm", "n", "<", gen_0, stoi(n));
     985        5915 :   if (typ(x) != t_INT) pari_err_TYPE("numtoperm",x);
     986        5915 :   v = Z_to_perm(n, x); settyp(v, t_VEC);
     987        5915 :   for (i = 1; i <= n; i++) gel(v,i) = utoipos(uel(v,i));
     988        5915 :   return v;
     989             : }
     990             : 
     991             : /* destroys v */
     992             : static GEN
     993        1701 : perm_to_Z_inplace(GEN v)
     994             : {
     995        1701 :   long l = lg(v), i, r;
     996        1701 :   GEN x = gen_0;
     997       10164 :   for (i = 1; i < l; i++)
     998             :   {
     999        8470 :     long vi = v[i];
    1000        8470 :     if (vi <= 0) return NULL;
    1001        8463 :     x = i==1 ? utoi(vi-1): addiu(muliu(x,l-i), vi-1);
    1002       25431 :     for (r = i+1; r < l; r++)
    1003       16968 :       if (v[r] > vi) v[r]--;
    1004             :   }
    1005        1694 :   return x;
    1006             : }
    1007             : GEN
    1008         840 : perm_to_Z(GEN v)
    1009             : {
    1010         840 :   pari_sp av = avma;
    1011         840 :   GEN x = perm_to_Z_inplace(leafcopy(v));
    1012         840 :   if (!x) pari_err_TYPE("permtonum",v);
    1013         840 :   return gerepileuptoint(av, x);
    1014             : }
    1015             : GEN
    1016        1708 : permtonum(GEN p)
    1017             : {
    1018        1708 :   pari_sp av = avma;
    1019             :   GEN v, x;
    1020        1708 :   switch(typ(p))
    1021             :   {
    1022         840 :     case t_VECSMALL: return perm_to_Z(p);
    1023             :     case t_VEC: case t_COL:
    1024         861 :       if (RgV_is_ZV(p)) { v = ZV_to_zv(p); break; }
    1025           7 :     default: pari_err_TYPE("permtonum",p); return NULL;
    1026             :   }
    1027         861 :   x = perm_to_Z_inplace(v);
    1028         861 :   if (!x) pari_err_TYPE("permtonum",p);
    1029         854 :   return gerepileuptoint(av, x);
    1030             : }
    1031             : 
    1032             : /*******************************************************************/
    1033             : /**                                                               **/
    1034             : /**                     RECIPROCAL POLYNOMIAL                     **/
    1035             : /**                                                               **/
    1036             : /*******************************************************************/
    1037             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1038             : GEN
    1039         112 : polrecip(GEN x)
    1040             : {
    1041         112 :   if (typ(x) != t_POL) pari_err_TYPE("polrecip",x);
    1042         112 :   return RgX_recip(x);
    1043             : }
    1044             : 
    1045             : /********************************************************************/
    1046             : /**                                                                **/
    1047             : /**                  POLYNOMIAL INTERPOLATION                      **/
    1048             : /**                                                                **/
    1049             : /********************************************************************/
    1050             : /* allow X = NULL for [1,...,n] */
    1051             : GEN
    1052         210 : RgV_polint(GEN X, GEN Y, long v)
    1053             : {
    1054         210 :   pari_sp av0 = avma, av;
    1055         210 :   GEN Q, P = NULL;
    1056         210 :   long i, l = lg(Y);
    1057         210 :   if (!X)
    1058             :   {
    1059          14 :     X = cgetg(l, t_VEC);
    1060          14 :     for (i=1; i<l; i++) gel(X,i) = utoipos(i);
    1061             :   }
    1062         210 :   Q = roots_to_pol(X, v); av = avma;
    1063         490 :   for (i=1; i<l; i++)
    1064             :   {
    1065             :     GEN inv, T, dP;
    1066         280 :     if (gequal0(gel(Y,i))) continue;
    1067         182 :     T = RgX_div_by_X_x(Q, gel(X,i), NULL);
    1068         182 :     inv = ginv(poleval(T,gel(X,i)));
    1069         182 :     dP = RgX_Rg_mul(T, gmul(gel(Y,i),inv));
    1070         182 :     P = P? RgX_add(P, dP): dP;
    1071         182 :     if (gc_needed(av,2))
    1072             :     {
    1073           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpV_polint");
    1074           0 :       P = gerepileupto(av, P);
    1075             :     }
    1076             :   }
    1077         210 :   if (!P) { avma = av; return zeropol(v); }
    1078         140 :   return gerepileupto(av0, P);
    1079             : }
    1080             : /* X,Y are "spec" GEN vectors with n > 0 components ( at X[0], ... X[n-1] ) */
    1081             : GEN
    1082       16503 : polint_i(GEN X, GEN Y, GEN x, long n, GEN *ptdy)
    1083             : {
    1084       16503 :   long i, m, ns = 0;
    1085       16503 :   pari_sp av = avma;
    1086       16503 :   GEN y, c, d, dy = NULL; /* gcc -Wall */
    1087             : 
    1088       16503 :   if (n == 1)
    1089             :   {
    1090          35 :     if (ptdy) *ptdy = gen_0;
    1091          35 :     return gmul(gel(Y,0), RgX_get_1(x));
    1092             :   }
    1093       16468 :   if (!X)
    1094             :   {
    1095           7 :     X = cgetg(n+1, t_VEC);
    1096           7 :     for (i=1; i<=n; i++) gel(X,i) = utoipos(i);
    1097           7 :     X++;
    1098             :   }
    1099       16468 :   switch(typ(x)) {
    1100             :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1101             :     {
    1102         504 :       GEN D = NULL;
    1103       20069 :       for (i=0; i<n; i++)
    1104             :       {
    1105       19579 :         GEN t = gsub(x,gel(X,i));
    1106       19579 :         switch(typ(t))
    1107             :         {
    1108             :           case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1109       19565 :             t = gabs(t, DEFAULTPREC);
    1110       19565 :             if (!D || gcmp(t,D) < 0) { ns = i; D = t; }
    1111       19565 :             break;
    1112             :           default:
    1113          14 :             goto NODY;
    1114             :         }
    1115             :       }
    1116         490 :       break;
    1117             :       /* X[ns] is closest to x */
    1118             :     }
    1119             : NODY:
    1120             :     default:
    1121       15978 :       if (ptdy) { *ptdy = gen_0; ptdy = NULL; }
    1122             :   }
    1123       16468 :   c = new_chunk(n);
    1124       16468 :   d = new_chunk(n); for (i=0; i<n; i++) gel(c,i) = gel(d,i) = gel(Y,i);
    1125       16468 :   y = gel(d,ns--);
    1126             :   /* divided differences */
    1127       83449 :   for (m=1; m<n; m++)
    1128             :   {
    1129      677013 :     for (i=0; i<n-m; i++)
    1130             :     {
    1131      610032 :       GEN ho = gsub(gel(X,i),x), hp = gsub(gel(X,i+m),x), den = gsub(ho,hp);
    1132      610032 :       if (gequal0(den))
    1133             :       {
    1134           7 :         char *x1 = stack_sprintf("X[%ld]", i+1);
    1135           7 :         char *x2 = stack_sprintf("X[%ld]", i+m+1);
    1136           7 :         pari_err_DOMAIN("polinterpolate",x1,"=",strtoGENstr(x2), X);
    1137             :       }
    1138      610025 :       den = gdiv(gsub(gel(c,i+1),gel(d,i)), den);
    1139      610025 :       gel(c,i) = gmul(ho,den);
    1140      610025 :       gel(d,i) = gmul(hp,den);
    1141             :     }
    1142       66981 :     dy = (2*(ns+1) < n-m)? gel(c,ns+1): gel(d,ns--);
    1143       66981 :     y = gadd(y,dy);
    1144             :   }
    1145       16461 :   if (!ptdy) return gerepileupto(av, y);
    1146         112 :   *ptdy = dy;
    1147         112 :   gerepileall(av, 2, &y, ptdy);
    1148         112 :   return y;
    1149             : }
    1150             : 
    1151             : GEN
    1152         658 : polint(GEN X, GEN Y, GEN t, GEN *ptdy)
    1153             : {
    1154         658 :   long lx = lg(X), vt;
    1155             : 
    1156         658 :   if (! is_vec_t(typ(X))) pari_err_TYPE("polinterpolate",X);
    1157         658 :   if (Y)
    1158             :   {
    1159         637 :     if (! is_vec_t(typ(Y))) pari_err_TYPE("polinterpolate",Y);
    1160         637 :     if (lx != lg(Y)) pari_err_DIM("polinterpolate");
    1161             :   }
    1162             :   else
    1163             :   {
    1164          21 :     Y = X;
    1165          21 :     X = NULL;
    1166             :   }
    1167         658 :   if (ptdy) *ptdy = gen_0;
    1168         658 :   vt = t? gvar(t): 0;
    1169         658 :   if (vt != NO_VARIABLE)
    1170             :   { /* formal interpolation */
    1171             :     pari_sp av;
    1172         210 :     long v0, vY = gvar(Y);
    1173             :     GEN P;
    1174         210 :     if (X) vY = varnmax(vY, gvar(X));
    1175             :     /* shortcut */
    1176         210 :     if (varncmp(vY, vt) > 0 && (!t || gequalX(t))) return RgV_polint(X, Y, vt);
    1177          84 :     av = avma;
    1178             :     /* first interpolate in high priority variable, then substitute t */
    1179          84 :     v0 = fetch_var_higher();
    1180          84 :     P = RgV_polint(X, Y, v0);
    1181          84 :     P = gsubst(P, v0, t? t: pol_x(0));
    1182          84 :     (void)delete_var();
    1183          84 :     return gerepileupto(av, P);
    1184             :   }
    1185             :   /* numerical interpolation */
    1186         448 :   if (lx == 1) return RgX_get_0(t);
    1187         434 :   return polint_i(X? X+1: NULL,Y+1,t,lx-1,ptdy);
    1188             : }
    1189             : 
    1190             : /********************************************************************/
    1191             : /**                                                                **/
    1192             : /**                       MODREVERSE                               **/
    1193             : /**                                                                **/
    1194             : /********************************************************************/
    1195             : static void
    1196           7 : err_reverse(GEN x, GEN T)
    1197             : {
    1198           7 :   pari_err_DOMAIN("modreverse","deg(minpoly(z))", "<", stoi(degpol(T)),
    1199             :                   mkpolmod(x,T));
    1200           0 : }
    1201             : 
    1202             : /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
    1203             : GEN
    1204          70 : RgXQ_reverse(GEN a, GEN T)
    1205             : {
    1206          70 :   pari_sp av = avma;
    1207          70 :   long n = degpol(T);
    1208             :   GEN y;
    1209             : 
    1210          70 :   if (n <= 1) {
    1211           7 :     if (n <= 0) return gcopy(a);
    1212           7 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1213             :   }
    1214          63 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1215          63 :   y = RgXV_to_RgM(RgXQ_powers(a,n-1,T), n);
    1216          63 :   y = RgM_solve(y, col_ei(n, 2));
    1217          63 :   if (!y) err_reverse(a,T);
    1218          56 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1219             : }
    1220             : GEN
    1221         637 : QXQ_reverse(GEN a, GEN T)
    1222             : {
    1223         637 :   pari_sp av = avma;
    1224         637 :   long n = degpol(T);
    1225             :   GEN y;
    1226             : 
    1227         637 :   if (n <= 1) {
    1228          14 :     if (n <= 0) return gcopy(a);
    1229          14 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1230             :   }
    1231         623 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1232         623 :   if (gequalX(a)) return gcopy(a);
    1233         609 :   y = RgXV_to_RgM(QXQ_powers(a,n-1,T), n);
    1234         609 :   y = RgM_solve(y, col_ei(n, 2));
    1235         609 :   if (!y) err_reverse(a,T);
    1236         609 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1237             : }
    1238             : 
    1239             : GEN
    1240          28 : modreverse(GEN x)
    1241             : {
    1242             :   long v, n;
    1243             :   GEN T, a, y;
    1244             : 
    1245          28 :   if (typ(x)!=t_POLMOD) pari_err_TYPE("modreverse",x);
    1246          28 :   T = gel(x,1); n = degpol(T); if (n <= 0) return gcopy(x);
    1247          21 :   a = gel(x,2);
    1248          21 :   v = varn(T);
    1249          21 :   y = cgetg(3,t_POLMOD);
    1250          21 :   gel(y,1) = (n==1)? gsub(pol_x(v), a): RgXQ_charpoly(a, T, v);
    1251          21 :   gel(y,2) = RgXQ_reverse(a, T); return y;
    1252             : }
    1253             : 
    1254             : /********************************************************************/
    1255             : /**                                                                **/
    1256             : /**                          MERGESORT                             **/
    1257             : /**                                                                **/
    1258             : /********************************************************************/
    1259             : static int
    1260    56225604 : cmp_small(GEN x, GEN y) {
    1261    56225604 :   long a = (long)x, b = (long)y;
    1262    56225604 :   return a>b? 1: (a<b? -1: 0);
    1263             : }
    1264             : 
    1265             : static int
    1266      147153 : veccmp(void *data, GEN x, GEN y)
    1267             : {
    1268      147153 :   GEN k = (GEN)data;
    1269      147153 :   long i, s, lk = lg(k), lx = minss(lg(x), lg(y));
    1270             : 
    1271      147153 :   if (!is_vec_t(typ(x))) pari_err_TYPE("lexicographic vecsort",x);
    1272      147153 :   if (!is_vec_t(typ(y))) pari_err_TYPE("lexicographic vecsort",y);
    1273      152886 :   for (i=1; i<lk; i++)
    1274             :   {
    1275      147181 :     long c = k[i];
    1276      147181 :     if (c >= lx)
    1277          14 :       pari_err_TYPE("lexicographic vecsort, index too large", stoi(c));
    1278      147167 :     s = lexcmp(gel(x,c), gel(y,c));
    1279      147167 :     if (s) return s;
    1280             :   }
    1281        5705 :   return 0;
    1282             : }
    1283             : 
    1284             : /* return permutation sorting v[1..n], removing duplicates. Assume n > 0 */
    1285             : static GEN
    1286      956998 : gen_sortspec_uniq(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1287             : {
    1288             :   pari_sp av;
    1289             :   long NX, nx, ny, m, ix, iy, i;
    1290             :   GEN x, y, w, W;
    1291             :   int s;
    1292      956998 :   switch(n)
    1293             :   {
    1294        1176 :     case 1: return mkvecsmall(1);
    1295             :     case 2:
    1296      318927 :       s = cmp(E,gel(v,1),gel(v,2));
    1297      318927 :       if      (s < 0) return mkvecsmall2(1,2);
    1298      165851 :       else if (s > 0) return mkvecsmall2(2,1);
    1299        2639 :       return mkvecsmall(1);
    1300             :     case 3:
    1301      160853 :       s = cmp(E,gel(v,1),gel(v,2));
    1302      160853 :       if (s < 0) {
    1303       79814 :         s = cmp(E,gel(v,2),gel(v,3));
    1304       79814 :         if (s < 0) return mkvecsmall3(1,2,3);
    1305       53865 :         else if (s == 0) return mkvecsmall2(1,2);
    1306       53473 :         s = cmp(E,gel(v,1),gel(v,3));
    1307       53473 :         if      (s < 0) return mkvecsmall3(1,3,2);
    1308       26488 :         else if (s > 0) return mkvecsmall3(3,1,2);
    1309        1729 :         return mkvecsmall2(1,2);
    1310       81039 :       } else if (s > 0) {
    1311       79429 :         s = cmp(E,gel(v,1),gel(v,3));
    1312       79429 :         if (s < 0) return mkvecsmall3(2,1,3);
    1313       53298 :         else if (s == 0) return mkvecsmall2(2,1);
    1314       52192 :         s = cmp(E,gel(v,2),gel(v,3));
    1315       52192 :         if (s < 0) return mkvecsmall3(2,3,1);
    1316       26572 :         else if (s > 0) return mkvecsmall3(3,2,1);
    1317         630 :         return mkvecsmall2(2,1);
    1318             :       } else {
    1319        1610 :         s = cmp(E,gel(v,1),gel(v,3));
    1320        1610 :         if (s < 0) return mkvecsmall2(1,3);
    1321        1015 :         else if (s == 0) return mkvecsmall(1);
    1322         798 :         return mkvecsmall2(3,1);
    1323             :       }
    1324             :   }
    1325      476042 :   NX = nx = n>>1; ny = n-nx;
    1326      476042 :   av = avma;
    1327      476042 :   x = gen_sortspec_uniq(v,   nx,E,cmp); nx = lg(x)-1;
    1328      476042 :   y = gen_sortspec_uniq(v+NX,ny,E,cmp); ny = lg(y)-1;
    1329      476042 :   w = cgetg(n+1, t_VECSMALL);
    1330      476042 :   m = ix = iy = 1;
    1331     7680995 :   while (ix<=nx && iy<=ny)
    1332             :   {
    1333     6728911 :     s = cmp(E, gel(v,x[ix]), gel(v,y[iy]+NX));
    1334     6728911 :     if (s < 0)
    1335     2889208 :       w[m++] = x[ix++];
    1336     3839703 :     else if (s > 0)
    1337     2983750 :       w[m++] = y[iy++]+NX;
    1338             :     else {
    1339      855953 :       w[m++] = x[ix++];
    1340      855953 :       iy++;
    1341             :     }
    1342             :   }
    1343      476042 :   while (ix<=nx) w[m++] = x[ix++];
    1344      476042 :   while (iy<=ny) w[m++] = y[iy++]+NX;
    1345      476042 :   avma = av;
    1346      476042 :   W = cgetg(m, t_VECSMALL);
    1347      476042 :   for (i = 1; i < m; i++) W[i] = w[i];
    1348      476042 :   return W;
    1349             : }
    1350             : 
    1351             : /* return permutation sorting v[1..n]. Assume n > 0 */
    1352             : static GEN
    1353   110455045 : gen_sortspec(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1354             : {
    1355             :   long nx, ny, m, ix, iy;
    1356             :   GEN x, y, w;
    1357   110455045 :   switch(n)
    1358             :   {
    1359             :     case 1:
    1360     2463346 :       (void)cmp(E,gel(v,1),gel(v,1)); /* check for type error */
    1361     2463346 :       return mkvecsmall(1);
    1362             :     case 2:
    1363    90381382 :       return cmp(E,gel(v,1),gel(v,2)) <= 0? mkvecsmall2(1,2)
    1364    45190684 :                                           : mkvecsmall2(2,1);
    1365             :     case 3:
    1366    21516147 :       if (cmp(E,gel(v,1),gel(v,2)) <= 0) {
    1367    15431459 :         if (cmp(E,gel(v,2),gel(v,3)) <= 0) return mkvecsmall3(1,2,3);
    1368     9363312 :         return (cmp(E,gel(v,1),gel(v,3)) <= 0)? mkvecsmall3(1,3,2)
    1369     4681656 :                                               : mkvecsmall3(3,1,2);
    1370             :       } else {
    1371     6084688 :         if (cmp(E,gel(v,1),gel(v,3)) <= 0) return mkvecsmall3(2,1,3);
    1372     7641906 :         return (cmp(E,gel(v,2),gel(v,3)) <= 0)? mkvecsmall3(2,3,1)
    1373     3820953 :                                               : mkvecsmall3(3,2,1);
    1374             :       }
    1375             :   }
    1376    41284854 :   nx = n>>1; ny = n-nx;
    1377    41284854 :   w = cgetg(n+1,t_VECSMALL);
    1378    41284854 :   x = gen_sortspec(v,   nx,E,cmp);
    1379    41284840 :   y = gen_sortspec(v+nx,ny,E,cmp);
    1380    41284841 :   m = ix = iy = 1;
    1381   332576666 :   while (ix<=nx && iy<=ny)
    1382   250006985 :     if (cmp(E, gel(v,x[ix]), gel(v,y[iy]+nx))<=0)
    1383   137326848 :       w[m++] = x[ix++];
    1384             :     else
    1385   112680136 :       w[m++] = y[iy++]+nx;
    1386    41284840 :   while (ix<=nx) w[m++] = x[ix++];
    1387    41284840 :   while (iy<=ny) w[m++] = y[iy++]+nx;
    1388    41284840 :   avma = (pari_sp)w; return w;
    1389             : }
    1390             : 
    1391             : static void
    1392    25977168 : init_sort(GEN *x, long *tx, long *lx)
    1393             : {
    1394    25977168 :   *tx = typ(*x);
    1395    25977168 :   if (*tx == t_LIST)
    1396             :   {
    1397          35 :     if (list_typ(*x)!=t_LIST_RAW) pari_err_TYPE("sort",*x);
    1398          35 :     *x = list_data(*x);
    1399          35 :     *lx = *x? lg(*x): 1;
    1400             :   } else {
    1401    25977133 :     if (!is_matvec_t(*tx) && *tx != t_VECSMALL) pari_err_TYPE("gen_sort",*x);
    1402    25977133 :     *lx = lg(*x);
    1403             :   }
    1404    25977168 : }
    1405             : 
    1406             : /* (x o y)[1..lx-1], destroy y */
    1407             : INLINE GEN
    1408     1402325 : sort_extract(GEN x, GEN y, long tx, long lx)
    1409             : {
    1410             :   long i;
    1411     1402325 :   switch(tx)
    1412             :   {
    1413             :     case t_VECSMALL:
    1414           7 :       for (i=1; i<lx; i++) y[i] = x[y[i]];
    1415           7 :       break;
    1416             :     case t_LIST:
    1417           7 :       settyp(y,t_VEC);
    1418           7 :       for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1419           7 :       return gtolist(y);
    1420             :     default:
    1421     1402311 :       settyp(y,tx);
    1422     1402311 :       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,y[i]));
    1423             :   }
    1424     1402318 :   return y;
    1425             : }
    1426             : 
    1427             : /* Sort the vector x, using cmp to compare entries. */
    1428             : GEN
    1429        4956 : gen_sort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1430             : {
    1431             :   long tx, lx;
    1432             :   GEN y;
    1433             : 
    1434        4956 :   init_sort(&x, &tx, &lx);
    1435        4956 :   if (lx==1) return tx == t_LIST? mklist(): cgetg(1, tx);
    1436        4788 :   y = gen_sortspec_uniq(x,lx-1,E,cmp);
    1437        4788 :   return sort_extract(x, y, tx, lg(y)); /* lg(y) <= lx */
    1438             : }
    1439             : /* Sort the vector x, using cmp to compare entries. */
    1440             : GEN
    1441     2104059 : gen_sort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1442             : {
    1443             :   long tx, lx;
    1444             :   GEN y;
    1445             : 
    1446     2104059 :   init_sort(&x, &tx, &lx);
    1447     2104059 :   if (lx==1) return tx == t_LIST? mklist(): cgetg(1, tx);
    1448     1397551 :   y = gen_sortspec(x,lx-1,E,cmp);
    1449     1397537 :   return sort_extract(x, y, tx, lx);
    1450             : }
    1451             : /* indirect sort: return the permutation that would sort x */
    1452             : GEN
    1453         126 : gen_indexsort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1454             : {
    1455             :   long tx, lx;
    1456         126 :   init_sort(&x, &tx, &lx);
    1457         126 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1458         126 :   return gen_sortspec_uniq(x,lx-1,E,cmp);
    1459             : }
    1460             : /* indirect sort: return the permutation that would sort x */
    1461             : GEN
    1462     1379803 : gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1463             : {
    1464             :   long tx, lx;
    1465     1379803 :   init_sort(&x, &tx, &lx);
    1466     1379803 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1467     1296134 :   return gen_sortspec(x,lx-1,E,cmp);
    1468             : }
    1469             : 
    1470             : /* Sort the vector x in place, using cmp to compare entries */
    1471             : void
    1472    22488224 : gen_sort_inplace(GEN x, void *E, int (*cmp)(void*,GEN,GEN), GEN *perm)
    1473             : {
    1474             :   long tx, lx, i;
    1475    22488224 :   pari_sp av = avma;
    1476             :   GEN y;
    1477             : 
    1478    22488224 :   init_sort(&x, &tx, &lx);
    1479    22488224 :   if (lx<=2)
    1480             :   {
    1481      161896 :     if (perm) *perm = lx == 1? cgetg(1, t_VECSMALL): mkvecsmall(1);
    1482    22650120 :     return;
    1483             :   }
    1484    22326328 :   y = gen_sortspec(x,lx-1, E, cmp);
    1485    22326328 :   if (perm)
    1486             :   {
    1487       13664 :     GEN z = new_chunk(lx);
    1488       13664 :     for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1489       13664 :     for (i=1; i<lx; i++) gel(x,i) = gel(z,i);
    1490       13664 :     *perm = y;
    1491       13664 :     avma = (pari_sp)y;
    1492             :   } else {
    1493    22312664 :     for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1494    22312664 :     for (i=1; i<lx; i++) gel(x,i) = gel(y,i);
    1495    22312664 :     avma = av;
    1496             :   }
    1497             : }
    1498             : 
    1499             : static int
    1500        7889 : closurecmp(void *data, GEN x, GEN y)
    1501             : {
    1502        7889 :   pari_sp av = avma;
    1503        7889 :   GEN z = closure_callgen2((GEN)data, x,y);
    1504        7889 :   if (typ(z) != t_INT)
    1505           0 :     pari_err_TYPE("closurecmp, cmp. fun. must return an integer", z);
    1506        7889 :   avma = av; return signe(z);
    1507             : }
    1508             : 
    1509             : static void
    1510         126 : check_positive_entries(GEN k)
    1511             : {
    1512         126 :   long i, l = lg(k);
    1513         287 :   for (i=1; i<l; i++)
    1514         161 :     if (k[i] <= 0) pari_err_DOMAIN("sort_function", "index", "<", gen_0, stoi(k[i]));
    1515         126 : }
    1516             : 
    1517             : typedef int (*CMP_FUN)(void*,GEN,GEN);
    1518             : static CMP_FUN
    1519      818090 : sort_function(void **E, GEN x, GEN k)
    1520             : {
    1521      818090 :   int (*cmp)(GEN,GEN) = &lexcmp;
    1522      818090 :   if (!k)
    1523             :   {
    1524      817453 :     *E = (void*)((typ(x) == t_VECSMALL)? cmp_small: cmp);
    1525      817453 :     return &cmp_nodata;
    1526             :   }
    1527         637 :   if (typ(x) == t_VECSMALL) pari_err_TYPE("sort_function", x);
    1528         623 :   switch(typ(k))
    1529             :   {
    1530          91 :     case t_INT: k = mkvecsmall(itos(k));  break;
    1531          35 :     case t_VEC: case t_COL: k = ZV_to_zv(k); break;
    1532           0 :     case t_VECSMALL: break;
    1533             :     case t_CLOSURE:
    1534         497 :      if (closure_arity(k) != 2 || closure_is_variadic(k))
    1535           0 :        pari_err_TYPE("sort_function, cmp. fun. needs exactly 2 arguments",k);
    1536         497 :      *E = (void*)k;
    1537         497 :      return &closurecmp;
    1538           0 :     default: pari_err_TYPE("sort_function",k);
    1539             :   }
    1540         126 :   check_positive_entries(k);
    1541         126 :   *E = (void*)k;
    1542         126 :   return &veccmp;
    1543             : }
    1544             : 
    1545             : #define cmp_IND 1
    1546             : #define cmp_LEX 2 /* FIXME: backward compatibility, ignored */
    1547             : #define cmp_REV 4
    1548             : #define cmp_UNIQ 8
    1549             : GEN
    1550         658 : vecsort0(GEN x, GEN k, long flag)
    1551             : {
    1552             :   void *E;
    1553         658 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
    1554             : 
    1555         651 :   if (flag < 0 || flag > (cmp_REV|cmp_LEX|cmp_IND|cmp_UNIQ))
    1556           0 :     pari_err_FLAG("vecsort");
    1557         651 :   if (flag & cmp_UNIQ)
    1558          28 :     x = flag & cmp_IND? gen_indexsort_uniq(x, E, CMP): gen_sort_uniq(x, E, CMP);
    1559             :   else
    1560         623 :     x = flag & cmp_IND? gen_indexsort(x, E, CMP): gen_sort(x, E, CMP);
    1561         637 :   if (flag & cmp_REV) { /* reverse order */
    1562             :     long j, lx;
    1563             :     GEN y;
    1564          21 :     if (typ(x)==t_LIST)
    1565             :     {
    1566           7 :       y = list_data(x);
    1567           7 :       if (!y) return x;
    1568             :     }
    1569             :     else
    1570          14 :       y = x;
    1571          14 :     lx = lg(y);
    1572          14 :     for (j=1; j<=(lx-1)>>1; j++) swap(gel(y,j), gel(y,lx-j));
    1573             :   }
    1574         630 :   return x;
    1575             : }
    1576             : 
    1577             : GEN
    1578       11637 : indexsort(GEN x) { return gen_indexsort(x, (void*)&gcmp, cmp_nodata); }
    1579             : GEN
    1580           0 : indexlexsort(GEN x) { return gen_indexsort(x, (void*)&lexcmp, cmp_nodata); }
    1581             : GEN
    1582          35 : indexvecsort(GEN x, GEN k)
    1583             : {
    1584          35 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1585          35 :   return gen_indexsort(x, (void*)k, &veccmp);
    1586             : }
    1587             : 
    1588             : GEN
    1589      601544 : sort(GEN x) { return gen_sort(x, (void*)gcmp, cmp_nodata); }
    1590             : GEN
    1591           0 : lexsort(GEN x) { return gen_sort(x, (void*)lexcmp, cmp_nodata); }
    1592             : GEN
    1593          56 : vecsort(GEN x, GEN k)
    1594             : {
    1595          56 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1596          56 :   return gen_sort(x, (void*)k, &veccmp);
    1597             : }
    1598             : long
    1599      817432 : vecsearch(GEN v, GEN x, GEN k)
    1600             : {
    1601      817432 :   pari_sp av = avma;
    1602             :   void *E;
    1603      817432 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, v, k);
    1604      817425 :   long r, tv = typ(v);
    1605      817425 :   if (tv == t_VECSMALL)
    1606          21 :     x = (GEN)itos(x);
    1607      817404 :   else if (!is_matvec_t(tv)) pari_err_TYPE("vecsearch", v);
    1608      817425 :   r = gen_search(v, x, 0, E, CMP);
    1609      817425 :   avma = av; return r;
    1610             : }
    1611             : 
    1612             : GEN
    1613         413 : ZV_indexsort(GEN L) { return gen_indexsort(L, (void*)&cmpii, &cmp_nodata); }
    1614             : GEN
    1615      781970 : ZV_sort(GEN L) { return gen_sort(L, (void*)&cmpii, &cmp_nodata); }
    1616             : GEN
    1617        2933 : ZV_sort_uniq(GEN L) { return gen_sort_uniq(L, (void*)&cmpii, &cmp_nodata); }
    1618             : void
    1619      363020 : ZV_sort_inplace(GEN L) { gen_sort_inplace(L, (void*)&cmpii, &cmp_nodata,NULL); }
    1620             : 
    1621             : /********************************************************************/
    1622             : /**                      SEARCH IN SORTED VECTOR                   **/
    1623             : /********************************************************************/
    1624             : /* index of x in table T, 0 otherwise */
    1625             : long
    1626    18372528 : tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
    1627             : {
    1628    18372528 :   long l = 1, u = lg(T)-1, i, s;
    1629             : 
    1630    95657224 :   while (u>=l)
    1631             :   {
    1632    63000953 :     i = (l+u)>>1; s = cmp(x, gel(T,i));
    1633    63000953 :     if (!s) return i;
    1634    58912168 :     if (s<0) u=i-1; else l=i+1;
    1635             :   }
    1636    14283743 :   return 0;
    1637             : }
    1638             : 
    1639             : /* looks if x belongs to the set T and returns the index if yes, 0 if no */
    1640             : long
    1641    15227478 : gen_search(GEN T, GEN x, long flag, void *data, int (*cmp)(void*,GEN,GEN))
    1642             : {
    1643    15227478 :   long lx = lg(T), i, l, u, s;
    1644             : 
    1645    15227478 :   if (lx==1) return flag? 1: 0;
    1646    15227478 :   l = 1; u = lx-1;
    1647             :   do
    1648             :   {
    1649    52861697 :     i = (l+u)>>1; s = cmp(data, x, gel(T,i));
    1650    52861697 :     if (!s) return flag? 0: i;
    1651    37713445 :     if (s<0) u=i-1; else l=i+1;
    1652    37713445 :   } while (u>=l);
    1653       79226 :   if (!flag) return 0;
    1654       41356 :   return (s<0)? i: i+1;
    1655             : }
    1656             : 
    1657             : long
    1658      920290 : ZV_search(GEN x, GEN y) { return tablesearch(x, y, cmpii); }
    1659             : 
    1660             : long
    1661    17447674 : zv_search(GEN x, long y) { return tablesearch(x, (GEN)y, cmp_small); }
    1662             : 
    1663             : /********************************************************************/
    1664             : /**                   COMPARISON FUNCTIONS                         **/
    1665             : /********************************************************************/
    1666             : int
    1667   367936030 : cmp_nodata(void *data, GEN x, GEN y)
    1668             : {
    1669   367936030 :   int (*cmp)(GEN,GEN)=(int (*)(GEN,GEN)) data;
    1670   367936030 :   return cmp(x,y);
    1671             : }
    1672             : 
    1673             : /* assume x and y come from the same idealprimedec call (uniformizer unique) */
    1674             : int
    1675      682705 : cmp_prime_over_p(GEN x, GEN y)
    1676             : {
    1677      682705 :   long k = pr_get_f(x) - pr_get_f(y); /* diff. between residue degree */
    1678      702714 :   return k? ((k > 0)? 1: -1)
    1679      702714 :           : ZV_cmp(pr_get_gen(x), pr_get_gen(y));
    1680             : }
    1681             : 
    1682             : int
    1683       48936 : cmp_prime_ideal(GEN x, GEN y)
    1684             : {
    1685       48936 :   int k = cmpii(pr_get_p(x), pr_get_p(y));
    1686       48936 :   return k? k: cmp_prime_over_p(x,y);
    1687             : }
    1688             : 
    1689             : /* assume x and y are t_POL in the same variable whose coeffs can be
    1690             :  * compared (used to sort polynomial factorizations) */
    1691             : int
    1692     1334180 : gen_cmp_RgX(void *data, GEN x, GEN y)
    1693             : {
    1694     1334180 :   int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
    1695     1334180 :   long i, lx = lg(x), ly = lg(y);
    1696             :   int fl;
    1697     1334180 :   if (lx > ly) return  1;
    1698     1273729 :   if (lx < ly) return -1;
    1699     3341802 :   for (i=lx-1; i>1; i--)
    1700     3139331 :     if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
    1701      202471 :   return 0;
    1702             : }
    1703             : 
    1704             : static int
    1705         378 : cmp_RgX_Rg(GEN x, GEN y)
    1706             : {
    1707         378 :   long lx = lgpol(x), ly;
    1708         378 :   if (lx > 1) return  1;
    1709           0 :   ly = gequal0(y) ? 0:1;
    1710           0 :   if (lx > ly) return  1;
    1711           0 :   if (lx < ly) return -1;
    1712           0 :   if (lx==0) return 0;
    1713           0 :   return gcmp(gel(x,2), y);
    1714             : }
    1715             : int
    1716       24773 : cmp_RgX(GEN x, GEN y)
    1717             : {
    1718       24773 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    1719       24773 :   if (typ(y) == t_POLMOD) y = gel(y,2);
    1720       24773 :   if (typ(x) == t_POL) {
    1721       17749 :     if (typ(y) != t_POL) return cmp_RgX_Rg(x, y);
    1722             :   } else {
    1723        7024 :     if (typ(y) != t_POL) return gcmp(x,y);
    1724         343 :     return - cmp_RgX_Rg(y,x);
    1725             :   }
    1726       17714 :   return gen_cmp_RgX((void*)&gcmp,x,y);
    1727             : }
    1728             : 
    1729             : int
    1730      176087 : cmp_Flx(GEN x, GEN y)
    1731             : {
    1732      176087 :   long i, lx = lg(x), ly = lg(y);
    1733      176087 :   if (lx > ly) return  1;
    1734      162108 :   if (lx < ly) return -1;
    1735      339134 :   for (i=lx-1; i>1; i--)
    1736      312363 :     if (uel(x,i) != uel(y,i)) return uel(x,i)<uel(y,i)? -1: 1;
    1737       26771 :   return 0;
    1738             : }
    1739             : /********************************************************************/
    1740             : /**                   MERGE & SORT FACTORIZATIONS                  **/
    1741             : /********************************************************************/
    1742             : /* merge fx, fy two factorizations, whose 1st column is sorted in strictly
    1743             :  * increasing order wrt cmp */
    1744             : GEN
    1745      962587 : merge_factor(GEN fx, GEN fy, void *data, int (*cmp)(void *,GEN,GEN))
    1746             : {
    1747      962587 :   GEN x = gel(fx,1), e = gel(fx,2), M, E;
    1748      962587 :   GEN y = gel(fy,1), f = gel(fy,2);
    1749      962587 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    1750             : 
    1751      962587 :   M = cgetg(l, t_COL);
    1752      962587 :   E = cgetg(l, t_COL);
    1753             : 
    1754      962587 :   m = ix = iy = 1;
    1755    11449617 :   while (ix<lx && iy<ly)
    1756             :   {
    1757     9524443 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    1758     9524443 :     if (s < 0)
    1759     8619021 :     { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; }
    1760      905422 :     else if (s == 0)
    1761             :     {
    1762      176732 :       GEN z = gel(x,ix), g = addii(gel(e,ix), gel(f,iy));
    1763      176732 :       iy++; ix++; if (!signe(g)) continue;
    1764       99704 :       gel(M,m) = z; gel(E,m) = g;
    1765             :     }
    1766             :     else
    1767      728690 :     { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; }
    1768     9447415 :     m++;
    1769             :   }
    1770      962587 :   while (ix<lx) { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; m++; }
    1771      962587 :   while (iy<ly) { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; m++; }
    1772      962587 :   setlg(M, m);
    1773      962587 :   setlg(E, m); return mkmat2(M, E);
    1774             : }
    1775             : /* merge two sorted vectors, removing duplicates. Shallow */
    1776             : GEN
    1777       33061 : merge_sort_uniq(GEN x, GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    1778             : {
    1779       33061 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    1780             :   GEN M;
    1781             : 
    1782       33061 :   M = cgetg(l, t_COL);
    1783       33061 :   m = ix = iy = 1;
    1784       87591 :   while (ix<lx && iy<ly)
    1785             :   {
    1786       21469 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    1787       21469 :     if (s < 0)
    1788       13636 :     { gel(M,m) = gel(x,ix); ix++; }
    1789        7833 :     else if (s == 0)
    1790           0 :     { gel(M,m) = gel(x,ix); iy++; ix++; }
    1791             :     else
    1792        7833 :     { gel(M,m) = gel(y,iy); iy++; }
    1793       21469 :     m++;
    1794             :   }
    1795       33061 :   while (ix<lx) { gel(M,m) = gel(x,ix); ix++; m++; }
    1796       33061 :   while (iy<ly) { gel(M,m) = gel(y,iy); iy++; m++; }
    1797       33061 :   setlg(M, m); return M;
    1798             : }
    1799             : 
    1800             : /* sort generic factorization, in place */
    1801             : GEN
    1802     2871771 : sort_factor(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    1803             : {
    1804             :   GEN a, b, A, B, w;
    1805             :   pari_sp av;
    1806             :   long n, i;
    1807             : 
    1808     2871771 :   a = gel(y,1); n = lg(a); if (n == 1) return y;
    1809     2865338 :   b = gel(y,2); av = avma;
    1810     2865338 :   A = new_chunk(n);
    1811     2865338 :   B = new_chunk(n);
    1812     2865338 :   w = gen_sortspec(a, n-1, data, cmp);
    1813     2865338 :   for (i=1; i<n; i++) { long k=w[i]; gel(A,i) = gel(a,k); gel(B,i) = gel(b,k); }
    1814     2865338 :   for (i=1; i<n; i++) { gel(a,i) = gel(A,i); gel(b,i) = gel(B,i); }
    1815     2865338 :   avma = av; return y;
    1816             : }
    1817             : /* sort polynomial factorization, in place */
    1818             : GEN
    1819      475738 : sort_factor_pol(GEN y,int (*cmp)(GEN,GEN))
    1820             : {
    1821      475738 :   (void)sort_factor(y,(void*)cmp, &gen_cmp_RgX);
    1822      475738 :   return y;
    1823             : }
    1824             : 
    1825             : /***********************************************************************/
    1826             : /*                                                                     */
    1827             : /*                          SET OPERATIONS                             */
    1828             : /*                                                                     */
    1829             : /***********************************************************************/
    1830             : GEN
    1831         147 : gtoset(GEN x)
    1832             : {
    1833             :   long lx;
    1834         147 :   if (!x) return cgetg(1, t_VEC);
    1835         147 :   switch(typ(x))
    1836             :   {
    1837             :     case t_VEC:
    1838         119 :     case t_COL: lx = lg(x); break;
    1839             :     case t_LIST:
    1840          14 :       if (list_typ(x)==t_LIST_MAP) return mapdomain(x);
    1841          14 :       x = list_data(x); lx = x? lg(x): 1; break;
    1842           7 :     case t_VECSMALL: lx = lg(x); x = zv_to_ZV(x); break;
    1843           7 :     default: return mkveccopy(x);
    1844             :   }
    1845         140 :   if (lx==1) return cgetg(1,t_VEC);
    1846         133 :   x = gen_sort_uniq(x, (void*)&cmp_universal, cmp_nodata);
    1847         133 :   settyp(x, t_VEC); /* it may be t_COL */
    1848         133 :   return x;
    1849             : }
    1850             : 
    1851             : long
    1852          14 : setisset(GEN x)
    1853             : {
    1854          14 :   long i, lx = lg(x);
    1855             : 
    1856          14 :   if (typ(x) != t_VEC) return 0;
    1857          14 :   if (lx == 1) return 1;
    1858          70 :   for (i=1; i<lx-1; i++)
    1859          63 :     if (cmp_universal(gel(x,i+1), gel(x,i)) <= 0) return 0;
    1860           7 :   return 1;
    1861             : }
    1862             : 
    1863             : long
    1864          35 : setsearch(GEN T, GEN y, long flag)
    1865             : {
    1866             :   long lx;
    1867          35 :   switch(typ(T))
    1868             :   {
    1869          21 :     case t_VEC: lx = lg(T); break;
    1870             :     case t_LIST:
    1871           7 :     if (list_typ(T) != t_LIST_RAW) pari_err_TYPE("setsearch",T);
    1872           7 :     T = list_data(T); lx = T? lg(T): 1; break;
    1873           7 :     default: pari_err_TYPE("setsearch",T);
    1874             :       return 0; /*LCOV_EXCL_LINE*/
    1875             :   }
    1876          28 :   if (lx==1) return flag? 1: 0;
    1877          28 :   return gen_search(T,y,flag,(void*)cmp_universal,cmp_nodata);
    1878             : }
    1879             : 
    1880             : GEN
    1881          35 : setunion(GEN x, GEN y)
    1882             : {
    1883          35 :   pari_sp av = avma;
    1884          35 :   long i, j, k, lx = lg(x), ly = lg(y);
    1885          35 :   GEN z = cgetg(lx + ly - 1, t_VEC);
    1886          35 :   if (typ(x) != t_VEC) pari_err_TYPE("setunion",x);
    1887          35 :   if (typ(y) != t_VEC) pari_err_TYPE("setunion",y);
    1888          35 :   i = j = k = 1;
    1889         161 :   while (i<lx && j<ly)
    1890             :   {
    1891          91 :     int s = cmp_universal(gel(x,i), gel(y,j));
    1892          91 :     if (s < 0)
    1893          28 :       z[k++] = x[i++];
    1894          63 :     else if (s > 0)
    1895          14 :       z[k++] = y[j++];
    1896             :     else {
    1897          49 :       z[k++] = x[i++];
    1898          49 :       j++;
    1899             :     }
    1900             :   }
    1901          35 :   while (i<lx) z[k++] = x[i++];
    1902          35 :   while (j<ly) z[k++] = y[j++];
    1903          35 :   setlg(z, k);
    1904          35 :   return gerepilecopy(av, z);
    1905             : }
    1906             : /* in case of equal keys in x,y, take the key from x */
    1907             : GEN
    1908      344018 : ZV_union_shallow(GEN x, GEN y)
    1909             : {
    1910      344018 :   long i, j, k, lx = lg(x), ly = lg(y);
    1911      344018 :   GEN z = cgetg(lx + ly - 1, t_VEC);
    1912      344018 :   i = j = k = 1;
    1913     1309054 :   while (i<lx && j<ly)
    1914             :   {
    1915      621018 :     int s = cmpii(gel(x,i), gel(y,j));
    1916      621018 :     if (s < 0)
    1917        9420 :       z[k++] = x[i++];
    1918      611598 :     else if (s > 0)
    1919      228127 :       z[k++] = y[j++];
    1920             :     else {
    1921      383471 :       z[k++] = x[i++];
    1922      383471 :       j++;
    1923             :     }
    1924             :   }
    1925      344018 :   while (i<lx) z[k++] = x[i++];
    1926      344018 :   while (j<ly) z[k++] = y[j++];
    1927      344018 :   setlg(z, k); return z;
    1928             : }
    1929             : 
    1930             : 
    1931             : GEN
    1932           7 : setintersect(GEN x, GEN y)
    1933             : {
    1934           7 :   long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
    1935           7 :   pari_sp av = avma;
    1936           7 :   GEN z = cgetg(lx,t_VEC);
    1937           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setintersect",x);
    1938           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setintersect",y);
    1939          77 :   while (ix < lx && iy < ly)
    1940             :   {
    1941          63 :     int c = cmp_universal(gel(x,ix), gel(y,iy));
    1942          63 :     if      (c < 0) ix++;
    1943          35 :     else if (c > 0) iy++;
    1944          21 :     else { gel(z, iz++) = gel(x,ix); ix++; iy++; }
    1945             :   }
    1946           7 :   setlg(z,iz); return gerepilecopy(av,z);
    1947             : }
    1948             : 
    1949             : GEN
    1950           7 : gen_setminus(GEN A, GEN B, int (*cmp)(GEN,GEN))
    1951             : {
    1952           7 :   pari_sp ltop = avma;
    1953           7 :   long i = 1, j = 1, k = 1, lx = lg(A), ly = lg(B);
    1954           7 :   GEN  diff = cgetg(lx,t_VEC);
    1955          91 :   while (i < lx && j < ly)
    1956          77 :     switch ( cmp(gel(A,i),gel(B,j)) )
    1957             :     {
    1958          28 :       case -1: gel(diff,k++) = gel(A,i++); break;
    1959          28 :       case 1: j++; break;
    1960          21 :       case 0: i++; break;
    1961             :     }
    1962           7 :   while (i < lx) gel(diff,k++) = gel(A,i++);
    1963           7 :   setlg(diff,k);
    1964           7 :   return gerepilecopy(ltop,diff);
    1965             : }
    1966             : 
    1967             : GEN
    1968           7 : setminus(GEN x, GEN y)
    1969             : {
    1970           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setminus",x);
    1971           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setminus",y);
    1972           7 :   return gen_setminus(x,y,cmp_universal);
    1973             : }
    1974             : 
    1975             : GEN
    1976          21 : setbinop(GEN f, GEN x, GEN y)
    1977             : {
    1978          21 :   pari_sp av = avma;
    1979          21 :   long i, j, lx, ly, k = 1;
    1980             :   GEN z;
    1981          21 :   if (typ(f) != t_CLOSURE || closure_arity(f) != 2 || closure_is_variadic(f))
    1982           7 :     pari_err_TYPE("setbinop [function needs exactly 2 arguments]",f);
    1983          14 :   lx = lg(x);
    1984          14 :   if (typ(x) != t_VEC) pari_err_TYPE("setbinop", x);
    1985          14 :   if (y == NULL) { /* assume x = y and f symmetric */
    1986           7 :     z = cgetg((((lx-1)*lx) >> 1) + 1, t_VEC);
    1987          28 :     for (i = 1; i < lx; i++)
    1988          63 :       for (j = i; j < lx; j++)
    1989          42 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(x,j));
    1990             :   } else {
    1991           7 :     ly = lg(y);
    1992           7 :     if (typ(y) != t_VEC) pari_err_TYPE("setbinop", y);
    1993           7 :     z = cgetg((lx-1)*(ly-1) + 1, t_VEC);
    1994          28 :     for (i = 1; i < lx; i++)
    1995          84 :       for (j = 1; j < ly; j++)
    1996          63 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(y,j));
    1997             :   }
    1998          14 :   return gerepileupto(av, gtoset(z));
    1999             : }

Generated by: LCOV version 1.11