Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - bibli2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28695-49bb1ac00f) Lines: 1251 1311 95.4 %
Date: 2023-09-24 07:47:42 Functions: 114 120 95.0 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : /*******************************************************************/
      19             : /**                                                               **/
      20             : /**                      SPECIAL POLYNOMIALS                      **/
      21             : /**                                                               **/
      22             : /*******************************************************************/
      23             : /* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
      24             :  * T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
      25             :  *   where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
      26             :  *   and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
      27             : GEN
      28        2156 : polchebyshev1(long n, long v) /* Assume 4*n < LONG_MAX */
      29             : {
      30             :   long k, l;
      31             :   pari_sp av;
      32             :   GEN q,a,r;
      33             : 
      34        2156 :   if (v<0) v = 0;
      35             :   /* polchebyshev(-n,1) = polchebyshev(n,1) */
      36        2156 :   if (n < 0) n = -n;
      37        2156 :   if (n==0) return pol_1(v);
      38        2135 :   if (n==1) return pol_x(v);
      39             : 
      40        2093 :   q = cgetg(n+3, t_POL); r = q + n+2;
      41        2093 :   a = int2n(n-1);
      42        2093 :   gel(r--,0) = a;
      43        2093 :   gel(r--,0) = gen_0;
      44       31955 :   for (k=1,l=n; l>1; k++,l-=2)
      45             :   {
      46       29862 :     av = avma;
      47       29862 :     a = diviuuexact(muluui(l, l-1, a), 4*k, n-k);
      48       29862 :     togglesign(a); a = gerepileuptoint(av, a);
      49       29862 :     gel(r--,0) = a;
      50       29862 :     gel(r--,0) = gen_0;
      51             :   }
      52        2093 :   q[1] = evalsigne(1) | evalvarn(v);
      53        2093 :   return q;
      54             : }
      55             : static void
      56          70 : polchebyshev1_eval_aux(long n, GEN x, GEN *pt1, GEN *pt2)
      57             : {
      58             :   GEN t1, t2, b;
      59          70 :   if (n == 1) { *pt1 = gen_1; *pt2 = x; return; }
      60          56 :   if (n == 0) { *pt1 = x; *pt2 = gen_1; return; }
      61          56 :   polchebyshev1_eval_aux((n+1) >> 1, x, &t1, &t2);
      62          56 :   b = gsub(gmul(gmul2n(t1,1), t2), x);
      63          56 :   if (odd(n)) { *pt1 = gadd(gmul2n(gsqr(t1), 1), gen_m1); *pt2 = b; }
      64          42 :   else        { *pt1 = b; *pt2 = gadd(gmul2n(gsqr(t2), 1), gen_m1); }
      65             : }
      66             : static GEN
      67          14 : polchebyshev1_eval(long n, GEN x)
      68             : {
      69             :   GEN t1, t2;
      70             :   long i, v;
      71             :   pari_sp av;
      72             : 
      73          14 :   if (n < 0) n = -n;
      74          14 :   if (n==0) return gen_1;
      75          14 :   if (n==1) return gcopy(x);
      76          14 :   av = avma;
      77          14 :   v = u_lvalrem(n, 2, (ulong*)&n);
      78          14 :   polchebyshev1_eval_aux((n+1)>>1, x, &t1, &t2);
      79          14 :   if (n != 1) t2 = gsub(gmul(gmul2n(t1,1), t2), x);
      80          35 :   for (i = 1; i <= v; i++) t2 = gadd(gmul2n(gsqr(t2), 1), gen_m1);
      81          14 :   return gerepileupto(av, t2);
      82             : }
      83             : 
      84             : /* Chebychev  polynomial of the second kind U(n,x): the coefficient in front of
      85             :  * x^(n-2*m) is (-1)^m * 2^(n-2m)*(n-m)!/m!/(n-2m)!  for m=0,1,...,n/2 */
      86             : GEN
      87        2135 : polchebyshev2(long n, long v)
      88             : {
      89             :   pari_sp av;
      90             :   GEN q, a, r;
      91             :   long m;
      92        2135 :   int neg = 0;
      93             : 
      94        2135 :   if (v<0) v = 0;
      95             :   /* polchebyshev(-n,2) = -polchebyshev(n-2,2) */
      96        2135 :   if (n < 0) {
      97        1050 :     if (n == -1) return zeropol(v);
      98        1029 :     neg = 1; n = -n-2;
      99             :   }
     100        2114 :   if (n==0) return neg ? scalar_ZX_shallow(gen_m1, v): pol_1(v);
     101             : 
     102        2072 :   q = cgetg(n+3, t_POL); r = q + n+2;
     103        2072 :   a = int2n(n);
     104        2072 :   if (neg) togglesign(a);
     105        2072 :   gel(r--,0) = a;
     106        2072 :   gel(r--,0) = gen_0;
     107       30807 :   for (m=1; 2*m<= n; m++)
     108             :   {
     109       28735 :     av = avma;
     110       28735 :     a = diviuuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m, n-m+1);
     111       28735 :     togglesign(a); a = gerepileuptoint(av, a);
     112       28735 :     gel(r--,0) = a;
     113       28735 :     gel(r--,0) = gen_0;
     114             :   }
     115        2072 :   q[1] = evalsigne(1) | evalvarn(v);
     116        2072 :   return q;
     117             : }
     118             : static void
     119          91 : polchebyshev2_eval_aux(long n, GEN x, GEN *pu1, GEN *pu2)
     120             : {
     121             :   GEN u1, u2, u, mu1;
     122          91 :   if (n == 1) { *pu1 = gen_1; *pu2 = gmul2n(x,1); return; }
     123          70 :   if (n == 0) { *pu1 = gen_0; *pu2 = gen_1; return; }
     124          70 :   polchebyshev2_eval_aux(n >> 1, x, &u1, &u2);
     125          70 :   mu1 = gneg(u1);
     126          70 :   u = gmul(gadd(u2,u1), gadd(u2,mu1));
     127          70 :   if (odd(n)) { *pu1 = u; *pu2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1)); }
     128          35 :   else        { *pu2 = u; *pu1 = gmul(gmul2n(u1,1), gadd(u2, gmul(x,mu1))); }
     129             : }
     130             : static GEN
     131          35 : polchebyshev2_eval(long n, GEN x)
     132             : {
     133             :   GEN u1, u2, mu1;
     134          35 :   long neg = 0;
     135             :   pari_sp av;
     136             : 
     137          35 :   if (n < 0) {
     138          14 :     if (n == -1) return gen_0;
     139           7 :     neg = 1; n = -n-2;
     140             :   }
     141          28 :   if (n==0) return neg ? gen_m1: gen_1;
     142          21 :   av = avma;
     143          21 :   polchebyshev2_eval_aux(n>>1, x, &u1, &u2);
     144          21 :   mu1 = gneg(u1);
     145          21 :   if (odd(n)) u2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1));
     146          14 :   else        u2 = gmul(gadd(u2,u1), gadd(u2,mu1));
     147          21 :   if (neg) u2 = gneg(u2);
     148          21 :   return gerepileupto(av, u2);
     149             : }
     150             : 
     151             : GEN
     152        4284 : polchebyshev(long n, long kind, long v)
     153             : {
     154        4284 :   switch (kind)
     155             :   {
     156        2149 :     case 1: return polchebyshev1(n, v);
     157        2135 :     case 2: return polchebyshev2(n, v);
     158           0 :     default: pari_err_FLAG("polchebyshev");
     159             :   }
     160             :   return NULL; /* LCOV_EXCL_LINE */
     161             : }
     162             : GEN
     163        4333 : polchebyshev_eval(long n, long kind, GEN x)
     164             : {
     165        4333 :   if (!x) return polchebyshev(n, kind, 0);
     166          63 :   if (gequalX(x)) return polchebyshev(n, kind, varn(x));
     167          49 :   switch (kind)
     168             :   {
     169          14 :     case 1: return polchebyshev1_eval(n, x);
     170          35 :     case 2: return polchebyshev2_eval(n, x);
     171           0 :     default: pari_err_FLAG("polchebyshev");
     172             :   }
     173             :   return NULL; /* LCOV_EXCL_LINE */
     174             : }
     175             : 
     176             : /* Hermite polynomial H(n,x):  H(n+1) = 2x H(n) - 2n H(n-1)
     177             :  * The coefficient in front of x^(n-2*m) is
     178             :  * (-1)^m * n! * 2^(n-2m)/m!/(n-2m)!  for m=0,1,...,n/2.. */
     179             : GEN
     180        1442 : polhermite(long n, long v)
     181             : {
     182             :   long m;
     183             :   pari_sp av;
     184             :   GEN q,a,r;
     185             : 
     186        1442 :   if (v<0) v = 0;
     187        1442 :   if (n==0) return pol_1(v);
     188             : 
     189        1435 :   q = cgetg(n+3, t_POL); r = q + n+2;
     190        1435 :   a = int2n(n);
     191        1435 :   gel(r--,0) = a;
     192        1435 :   gel(r--,0) = gen_0;
     193       40327 :   for (m=1; 2*m<= n; m++)
     194             :   {
     195       38892 :     av = avma;
     196       38892 :     a = diviuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m);
     197       38892 :     togglesign(a);
     198       38892 :     gel(r--,0) = a = gerepileuptoint(av, a);
     199       38892 :     gel(r--,0) = gen_0;
     200             :   }
     201        1435 :   q[1] = evalsigne(1) | evalvarn(v);
     202        1435 :   return q;
     203             : }
     204             : static void
     205          21 : err_hermite(long n)
     206          21 : { pari_err_DOMAIN("polhermite", "degree", "<", gen_0, stoi(n)); }
     207             : GEN
     208        1477 : polhermite_eval0(long n, GEN x, long flag)
     209             : {
     210             :   long i;
     211             :   pari_sp av, av2;
     212             :   GEN x2, u, v;
     213             : 
     214        1477 :   if (n < 0) err_hermite(n);
     215        1470 :   if (!x || gequalX(x))
     216             :   {
     217        1442 :     long v = x? varn(x): 0;
     218        1442 :     if (flag)
     219             :     {
     220          14 :       if (!n) err_hermite(-1);
     221           7 :       retmkvec2(polhermite(n-1,v),polhermite(n,v));
     222             :     }
     223        1428 :     return polhermite(n, v);
     224             :   }
     225          28 :   if (n==0)
     226             :   {
     227           7 :     if (flag) err_hermite(-1);
     228           0 :     return gen_1;
     229             :   }
     230          21 :   if (n==1)
     231             :   {
     232           0 :     if (flag) retmkvec2(gen_1, gmul2n(x,1));
     233           0 :     return gmul2n(x,1);
     234             :   }
     235          21 :   av = avma; x2 = gmul2n(x,1); v = gen_1; u = x2;
     236          21 :   av2= avma;
     237        7070 :   for (i=1; i<n; i++)
     238             :   { /* u = H_i(x), v = H_{i-1}(x), compute t = H_{i+1}(x) */
     239             :     GEN t;
     240        7049 :     if ((i & 0xff) == 0) gerepileall(av2,2,&u, &v);
     241        7049 :     t = gsub(gmul(x2, u), gmulsg(2*i,v));
     242        7049 :     v = u; u = t;
     243             :   }
     244          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     245          14 :   return gerepileupto(av, u);
     246             : }
     247             : GEN
     248           0 : polhermite_eval(long n, GEN x) { return polhermite_eval0(n, x, 0); }
     249             : 
     250             : /* Legendre polynomial
     251             :  * L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1)
     252             :  * L(n) = 2^-n sum_{k=0}^{n/2} a_k x^(n-2k)
     253             :  *   where a_k = (-1)^k (2n-2k)! / k! (n-k)! (n-2k)! is an integer
     254             :  *   and a_0 = binom(2n,n), a_k / a_{k-1} = - (n-2k+1)(n-2k+2) / 2k (2n-2k+1) */
     255             : GEN
     256        2163 : pollegendre(long n, long v)
     257             : {
     258             :   long k, l;
     259             :   pari_sp av;
     260             :   GEN a, r, q;
     261             : 
     262        2163 :   if (v<0) v = 0;
     263             :   /* pollegendre(-n) = pollegendre(n-1) */
     264        2163 :   if (n < 0) n = -n-1;
     265        2163 :   if (n==0) return pol_1(v);
     266        2121 :   if (n==1) return pol_x(v);
     267             : 
     268        2079 :   av = avma;
     269        2079 :   q = cgetg(n+3, t_POL); r = q + n+2;
     270        2079 :   gel(r--,0) = a = binomialuu(n<<1,n);
     271        2079 :   gel(r--,0) = gen_0;
     272       31423 :   for (k=1,l=n; l>1; k++,l-=2)
     273             :   { /* l = n-2*k+2 */
     274       29344 :     av = avma;
     275       29344 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     276       29344 :     togglesign(a); a = gerepileuptoint(av, a);
     277       29344 :     gel(r--,0) = a;
     278       29344 :     gel(r--,0) = gen_0;
     279             :   }
     280        2079 :   q[1] = evalsigne(1) | evalvarn(v);
     281        2079 :   return gerepileupto(av, gmul2n(q,-n));
     282             : }
     283             : /* q such that Ln * 2^n = q(x^2) [n even] or x q(x^2) [n odd] */
     284             : GEN
     285           0 : pollegendre_reduced(long n, long v)
     286             : {
     287             :   long k, l, N;
     288             :   pari_sp av;
     289             :   GEN a, r, q;
     290             : 
     291           0 :   if (v<0) v = 0;
     292             :   /* pollegendre(-n) = pollegendre(n-1) */
     293           0 :   if (n < 0) n = -n-1;
     294           0 :   if (n<=1) return n? scalarpol_shallow(gen_2,v): pol_1(v);
     295             : 
     296           0 :   N = n >> 1;
     297           0 :   q = cgetg(N+3, t_POL); r = q + N+2;
     298           0 :   gel(r--,0) = a = binomialuu(n<<1,n);
     299           0 :   for (k=1,l=n; l>1; k++,l-=2)
     300             :   { /* l = n-2*k+2 */
     301           0 :     av = avma;
     302           0 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     303           0 :     togglesign(a);
     304           0 :     gel(r--,0) = a = gerepileuptoint(av, a);
     305             :   }
     306           0 :   q[1] = evalsigne(1) | evalvarn(v);
     307           0 :   return q;
     308             : }
     309             : 
     310             : GEN
     311        2177 : pollegendre_eval0(long n, GEN x, long flag)
     312             : {
     313             :   pari_sp av;
     314             :   GEN u, v;
     315             :   long i;
     316             : 
     317        2177 :   if (n < 0) n = -n-1; /* L(-n) = L(n-1) */
     318             :   /* n >= 0 */
     319        2177 :   if (flag && flag != 1) pari_err_FLAG("pollegendre");
     320        2177 :   if (!x || gequalX(x))
     321             :   {
     322        2156 :     long v = x? varn(x): 0;
     323        2156 :     if (flag) retmkvec2(pollegendre(n-1,v), pollegendre(n,v));
     324        2149 :     return pollegendre(n, v);
     325             :   }
     326          21 :   if (n==0)
     327             :   {
     328           0 :     if (flag) retmkvec2(gen_1, gcopy(x));
     329           0 :     return gen_1;
     330             :   }
     331          21 :   if (n==1)
     332             :   {
     333           0 :     if (flag) retmkvec2(gcopy(x), gen_1);
     334           0 :     return gcopy(x);
     335             :   }
     336          21 :   av = avma; v = gen_1; u = x;
     337        7070 :   for (i=1; i<n; i++)
     338             :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     339             :     GEN t;
     340        7049 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     341        7049 :     t = gdivgu(gsub(gmul(gmulsg(2*i+1,x), u), gmulsg(i,v)), i+1);
     342        7049 :     v = u; u = t;
     343             :   }
     344          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     345          14 :   return gerepileupto(av, u);
     346             : }
     347             : GEN
     348           0 : pollegendre_eval(long n, GEN x) { return pollegendre_eval0(n, x, 0); }
     349             : 
     350             : /* Laguerre polynomial
     351             :  * L0^a = 1; L1^a = -X+a+1;
     352             :  * (n+1)*L^a(n+1) = (-X+(2*n+a+1))*L^a(n) - (n+a)*L^a(n-1)
     353             :  * L^a(n) = sum_{k=0}^n (-1)^k * binom(n+a,n-k) * x^k/k! */
     354             : GEN
     355        2128 : pollaguerre(long n, GEN a, long v)
     356             : {
     357        2128 :   pari_sp av = avma;
     358        2128 :   GEN L = cgetg(n+3, t_POL), c1 = gen_1, c2 = mpfact(n);
     359             :   long i;
     360             : 
     361        2128 :   L[1] = evalsigne(1) | evalvarn(v);
     362        2128 :   if (odd(n)) togglesign_safe(&c2);
     363      117404 :   for (i = n; i >= 0; i--)
     364             :   {
     365      115276 :     gel(L, i+2) = gdiv(c1, c2);
     366      115276 :     if (i)
     367             :     {
     368      113148 :       c2 = divis(c2,-i);
     369      113148 :       c1 = gdivgu(gmul(c1, gaddsg(i,a)), n+1-i);
     370             :     }
     371             :   }
     372        2128 :   return gerepilecopy(av, L);
     373             : }
     374             : static void
     375          21 : err_lag(long n)
     376          21 : { pari_err_DOMAIN("pollaguerre", "degree", "<", gen_0, stoi(n)); }
     377             : GEN
     378        2163 : pollaguerre_eval0(long n, GEN a, GEN x, long flag)
     379             : {
     380        2163 :   pari_sp av = avma;
     381             :   long i;
     382             :   GEN v, u;
     383             : 
     384        2163 :   if (n < 0) err_lag(n);
     385        2156 :   if (flag && flag != 1) pari_err_FLAG("pollaguerre");
     386        2156 :   if (!a) a = gen_0;
     387        2156 :   if (!x || gequalX(x))
     388             :   {
     389        2128 :     long v = x? varn(x): 0;
     390        2128 :     if (flag)
     391             :     {
     392          14 :       if (!n) err_lag(-1);
     393           7 :       retmkvec2(pollaguerre(n-1,a,v), pollaguerre(n,a,v));
     394             :     }
     395        2114 :     return pollaguerre(n,a,v);
     396             :   }
     397          28 :   if (n==0)
     398             :   {
     399           7 :     if (flag) err_lag(-1);
     400           0 :     return gen_1;
     401             :   }
     402          21 :   if (n==1)
     403             :   {
     404           0 :     if (flag) retmkvec2(gsub(gaddgs(a,1),x), gen_1);
     405           0 :     return gsub(gaddgs(a,1),x);
     406             :   }
     407          21 :   av = avma; v = gen_1; u = gsub(gaddgs(a,1),x);
     408        7070 :   for (i=1; i<n; i++)
     409             :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     410             :     GEN t;
     411        7049 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     412        7049 :     t = gdivgu(gsub(gmul(gsub(gaddsg(2*i+1,a),x), u), gmul(gaddsg(i,a),v)), i+1);
     413        7049 :     v = u; u = t;
     414             :   }
     415          21 :   if (flag) return gerepilecopy(av, mkvec2(v, u));
     416          14 :   return gerepileupto(av, u);
     417             : }
     418             : GEN
     419           0 : pollaguerre_eval(long n, GEN x, GEN a) { return pollaguerre_eval0(n, x, a, 0); }
     420             : 
     421             : /* polcyclo(p) = X^(p-1) + ... + 1 */
     422             : static GEN
     423      505820 : polcyclo_prime(long p, long v)
     424             : {
     425      505820 :   GEN T = cgetg(p+2, t_POL);
     426             :   long i;
     427      505820 :   T[1] = evalsigne(1) | evalvarn(v);
     428     3438405 :   for (i = 2; i < p+2; i++) gel(T,i) = gen_1;
     429      505820 :   return T;
     430             : }
     431             : 
     432             : /* cyclotomic polynomial */
     433             : GEN
     434      628228 : polcyclo(long n, long v)
     435             : {
     436             :   long s, q, i, l;
     437      628228 :   pari_sp av=avma;
     438             :   GEN T, P;
     439             : 
     440      628228 :   if (v<0) v = 0;
     441      628228 :   if (n < 3)
     442      122408 :     switch(n)
     443             :     {
     444       30821 :       case 1: return deg1pol_shallow(gen_1, gen_m1, v);
     445       91587 :       case 2: return deg1pol_shallow(gen_1, gen_1, v);
     446           0 :       default: pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     447             :     }
     448      505820 :   P = gel(factoru(n), 1); l = lg(P);
     449      505820 :   s = P[1]; T = polcyclo_prime(s, v);
     450      794220 :   for (i = 2; i < l; i++)
     451             :   { /* Phi_{np}(X) = Phi_n(X^p) / Phi_n(X) */
     452      288400 :     s *= P[i];
     453      288400 :     T = RgX_div(RgX_inflate(T, P[i]), T);
     454             :   }
     455             :   /* s = squarefree part of n */
     456      505820 :   q = n / s;
     457      505820 :   if (q == 1) return gerepileupto(av, T);
     458      243378 :   return gerepilecopy(av, RgX_inflate(T,q));
     459             : }
     460             : 
     461             : /* cyclotomic polynomial */
     462             : GEN
     463      100186 : polcyclo_eval(long n, GEN x)
     464             : {
     465      100186 :   pari_sp av= avma;
     466             :   GEN P, md, xd, yneg, ypos;
     467      100186 :   long vpx, l, s, i, j, q, tx, root_of_1 = 0;
     468             : 
     469      100186 :   if (!x) return polcyclo(n, 0);
     470       15101 :   tx = typ(x);
     471       15101 :   if (gequalX(x)) return polcyclo(n, varn(x));
     472       14506 :   if (n <= 0) pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     473       14506 :   if (n == 1) return gsubgs(x, 1);
     474       14506 :   if (tx == t_INT && !signe(x)) return gen_1;
     475       15682 :   while ((n & 3) == 0) { n >>= 1; x = gsqr(x); } /* Phi_4n(x) = Phi_2n(x^2) */
     476             :   /* n not divisible by 4 */
     477       14506 :   if (n == 2) return gerepileupto(av, gaddgs(x,1));
     478        6344 :   if (!odd(n)) { n >>= 1; x = gneg(x); } /* Phi_2n(x) = Phi_n(-x) for n>1 odd */
     479             :   /* n odd > 2.  s largest squarefree divisor of n */
     480        6344 :   P = gel(factoru(n), 1); s = zv_prod(P);
     481             :   /* replace n by largest squarefree divisor */
     482        6344 :   q = n/s; if (q != 1) { x = gpowgs(x, q); n = s; }
     483        6344 :   l = lg(P)-1;
     484             :   /* n squarefree odd > 2, l distinct prime divisors. Now handle x = 1 or -1 */
     485        6344 :   if (tx == t_INT) { /* shortcut */
     486        1708 :     if (is_pm1(x))
     487             :     {
     488          56 :       set_avma(av);
     489          56 :       if (signe(x) > 0 && l == 1) return utoipos(P[1]);
     490          35 :       return gen_1;
     491             :     }
     492             :   } else {
     493        4636 :     if (gequal1(x))
     494             :     { /* n is prime, return n; multiply by x to keep the type */
     495          14 :       if (l == 1) return gerepileupto(av, gmulgu(x,n));
     496           7 :       return gerepilecopy(av, x); /* else 1 */
     497             :     }
     498        4622 :     if (gequalm1(x)) return gerepileupto(av, gneg(x)); /* -1 */
     499             :   }
     500             :   /* Heuristic: evaluation will probably not improve things */
     501        6267 :   if (tx == t_POL || tx == t_MAT || lg(x) > n)
     502          24 :     return gerepileupto(av, poleval(polcyclo(n,0), x));
     503             : 
     504        6243 :   xd = cgetg((1L<<l) + 1, t_VEC); /* the x^d, where d | n */
     505        6243 :   md = cgetg((1L<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
     506        6243 :   gel(xd, 1) = x;
     507        6243 :   md[1] = 1;
     508             :   /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d).
     509             :    * If x has exact order D, n = Dq, then the result is 0 if q = 1. Otherwise
     510             :    * the factors with x^d-1, D|d are omitted and we multiply at the end by
     511             :    *   prod_{d | q} d^mu(q/d) = q if prime, 1 otherwise */
     512             :   /* We store the factors with mu(d)= 1 (resp.-1) in ypos (resp yneg).
     513             :    * At the end we return ypos/yneg if mu(n)=1 and yneg/ypos if mu(n)=-1 */
     514        6243 :   ypos = gsubgs(x,1);
     515        6243 :   yneg = gen_1;
     516        6243 :   vpx = (typ(x) == t_PADIC)? valp(x): 0;
     517       13809 :   for (i = 1; i <= l; i++)
     518             :   {
     519        7566 :     long ti = 1L<<(i-1), p = P[i];
     520       16553 :     for (j = 1; j <= ti; j++) {
     521        8987 :       GEN X = gel(xd,j), t;
     522        8987 :       if (vpx > 0)
     523             :       { /* ypos, X t_PADIC */
     524          98 :         ulong a = umuluu_or_0(p, valp(X)), b = precp(ypos) - 1;
     525          98 :         long e = (a && a < b) ? b - a : 0;
     526          98 :         if (precp(X) > e) X = cvtop(X, gel(ypos,2), e);
     527          98 :         if (e > 0) X = gpowgs(X, p); /* avoid valp overflow of p-adic 0*/
     528             :       }
     529             :       else
     530        8889 :         X = gpowgs(X, p);
     531        8987 :       md[ti+j] = -md[j];
     532        8987 :       gel(xd,ti+j) = X;
     533             :       /* avoid precp overflow */
     534        8987 :       t = (vpx > 0 && gequal0(X))? gen_m1: gsubgs(X,1);
     535        8987 :       if (gequal0(t))
     536             :       { /* x^d = 1; root_of_1 := the smallest index ti+j such that X == 1
     537             :         * (whose bits code d: bit i-1 is set iff P[i] | d). If no such index
     538             :         * exists, then root_of_1 remains 0. Do not multiply with X-1 if X = 1,
     539             :         * we handle these factors at the end */
     540          28 :         if (!root_of_1) root_of_1 = ti+j;
     541             :       }
     542             :       else
     543             :       {
     544        8959 :         if (md[ti+j] == 1) ypos = gmul(ypos, t);
     545        7601 :         else               yneg = gmul(yneg, t);
     546             :       }
     547             :     }
     548             :   }
     549        6243 :   ypos = odd(l)? gdiv(yneg,ypos): gdiv(ypos,yneg);
     550        6243 :   if (root_of_1)
     551             :   {
     552          21 :     GEN X = gel(xd,(1<<l)); /* = x^n = 1 */
     553          21 :     long bitmask_q = (1<<l) - root_of_1;
     554             :     /* bitmask_q encodes q = n/d: bit (i-1) is 1 iff P[i] | q */
     555             : 
     556             :     /* x is a root of unity.  If bitmask_q = 0, then x was a primitive n-th
     557             :      * root of 1 and the result is zero. Return X - 1 to preserve type. */
     558          21 :     if (!bitmask_q) return gerepileupto(av, gsubgs(X, 1));
     559             :     /* x is a primitive d-th root of unity, where d|n and d<n: we
     560             :      * must multiply ypos by if(isprime(n/d), n/d, 1) */
     561           7 :     ypos = gmul(ypos, X); /* multiply by X = 1 to preserve type */
     562             :     /* If bitmask_q = 1<<(i-1) for some i <= l, then q == P[i] and we multiply
     563             :      * by P[i]; otherwise q is composite and nothing more needs to be done */
     564           7 :     if (!(bitmask_q & (bitmask_q-1))) /* detects power of 2, since bitmask!=0 */
     565             :     {
     566           7 :       i = vals(bitmask_q)+1; /* q = P[i] */
     567           7 :       ypos = gmulgu(ypos, P[i]);
     568             :     }
     569             :   }
     570        6229 :   return gerepileupto(av, ypos);
     571             : }
     572             : /********************************************************************/
     573             : /**                                                                **/
     574             : /**                  HILBERT & PASCAL MATRICES                     **/
     575             : /**                                                                **/
     576             : /********************************************************************/
     577             : GEN
     578         133 : mathilbert(long n) /* Hilbert matrix of order n */
     579             : {
     580             :   long i,j;
     581             :   GEN p;
     582             : 
     583         133 :   if (n < 0) pari_err_DOMAIN("mathilbert", "dimension", "<", gen_0, stoi(n));
     584         133 :   p = cgetg(n+1,t_MAT);
     585        1120 :   for (j=1; j<=n; j++)
     586             :   {
     587         987 :     gel(p,j) = cgetg(n+1,t_COL);
     588       16583 :     for (i=1+(j==1); i<=n; i++)
     589       15596 :       gcoeff(p,i,j) = mkfrac(gen_1, utoipos(i+j-1));
     590             :   }
     591         133 :   if (n) gcoeff(p,1,1) = gen_1;
     592         133 :   return p;
     593             : }
     594             : 
     595             : /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
     596             : GEN
     597       34283 : matqpascal(long n, GEN q)
     598             : {
     599             :   long i, j, I;
     600       34283 :   pari_sp av = avma;
     601       34283 :   GEN m, qpow = NULL; /* gcc -Wall */
     602             : 
     603       34283 :   if (n < -1)  pari_err_DOMAIN("matpascal", "n", "<", gen_m1, stoi(n));
     604       34283 :   n++; m = cgetg(n+1,t_MAT);
     605      143951 :   for (j=1; j<=n; j++) gel(m,j) = cgetg(n+1,t_COL);
     606       34285 :   if (q)
     607             :   {
     608          42 :     I = (n+1)/2;
     609          42 :     if (I > 1) { qpow = new_chunk(I+1); gel(qpow,2)=q; }
     610          84 :     for (j=3; j<=I; j++) gel(qpow,j) = gmul(q, gel(qpow,j-1));
     611             :   }
     612      143957 :   for (i=1; i<=n; i++)
     613             :   {
     614      109672 :     I = (i+1)/2; gcoeff(m,i,1)= gen_1;
     615      109672 :     if (q)
     616             :     {
     617         483 :       for (j=2; j<=I; j++)
     618         238 :         gcoeff(m,i,j) = gadd(gmul(gel(qpow,j),gcoeff(m,i-1,j)),
     619         238 :                              gcoeff(m,i-1,j-1));
     620             :     }
     621             :     else
     622             :     {
     623     1092349 :       for (j=2; j<=I; j++)
     624      982922 :         gcoeff(m,i,j) = addii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
     625             :     }
     626     1146927 :     for (   ; j<=i; j++) gcoeff(m,i,j) = gcoeff(m,i,i+1-j);
     627     2130084 :     for (   ; j<=n; j++) gcoeff(m,i,j) = gen_0;
     628             :   }
     629       34285 :   return gerepilecopy(av, m);
     630             : }
     631             : 
     632             : GEN
     633          77 : eulerianpol(long N, long v)
     634             : {
     635          77 :   pari_sp av = avma;
     636          77 :   long n, n2, k = 0;
     637             :   GEN A;
     638          77 :   if (v < 0) v = 0;
     639          77 :   if (N < 0) pari_err_DOMAIN("eulerianpol", "index", "<", gen_0, stoi(N));
     640          70 :   if (N <= 1) return pol_1(v);
     641          42 :   if (N == 2) return deg1pol_shallow(gen_1, gen_1, v);
     642          35 :   A = cgetg(N+1, t_VEC);
     643          35 :   gel(A,1) = gen_1; gel(A,2) = gen_1; /* A_2 = x+1 */
     644         567 :   for (n = 3; n <= N; n++)
     645             :   { /* A(n,k) = (n-k)A(n-1,k-1) + (k+1)A(n-1,k) */
     646         532 :     n2 = n >> 1;
     647         532 :     if (odd(n)) gel(A,n2+1) = mului(n+1, gel(A,n2));
     648        8652 :     for (k = n2-1; k; k--)
     649        8120 :       gel(A,k+1) = addii(mului(n-k, gel(A,k)), mului(k+1, gel(A,k+1)));
     650         532 :     if (gc_needed(av,1))
     651             :     {
     652           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"eulerianpol, %ld/%ld",n,N);
     653           0 :       for (k = odd(n)? n2+1: n2; k < N; k++) gel(A,k+1) = gen_0;
     654           0 :       A = gerepilecopy(av, A);
     655             :     }
     656             :   }
     657          35 :   k = N >> 1; if (odd(N)) k++;
     658         329 :   for (; k < N; k++) gel(A,k+1) = gel(A, N-k);
     659          35 :   return gerepilecopy(av, RgV_to_RgX(A, v));
     660             : }
     661             : 
     662             : /******************************************************************/
     663             : /**                                                              **/
     664             : /**                       PRECISION CHANGES                      **/
     665             : /**                                                              **/
     666             : /******************************************************************/
     667             : 
     668             : GEN
     669          98 : gprec(GEN x, long d)
     670             : {
     671          98 :   pari_sp av = avma;
     672          98 :   if (d <= 0) pari_err_DOMAIN("gprec", "precision", "<=", gen_0, stoi(d));
     673          98 :   return gerepilecopy(av, gprec_w(x, ndec2prec(d)));
     674             : }
     675             : 
     676             : /* not GC-safe; precision given in word length (including codewords) */
     677             : GEN
     678    10183146 : gprec_w(GEN x, long pr)
     679             : {
     680             :   long lx, i;
     681             :   GEN y;
     682             : 
     683    10183146 :   switch(typ(x))
     684             :   {
     685     6741920 :     case t_REAL:
     686     6741920 :       if (signe(x)) return realprec(x) != pr? rtor(x,pr): x;
     687       62988 :       i = -prec2nbits(pr);
     688       62987 :       if (i < expo(x)) return real_0_bit(i);
     689       60596 :       y = cgetr(2); y[1] = x[1]; return y;
     690     1868419 :     case t_COMPLEX:
     691     1868419 :       y = cgetg(3, t_COMPLEX);
     692     1868419 :       gel(y,1) = gprec_w(gel(x,1),pr);
     693     1868419 :       gel(y,2) = gprec_w(gel(x,2),pr);
     694     1868420 :       break;
     695      401300 :    case t_POL: case t_SER:
     696      401300 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     697     2854346 :       for (i=2; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     698      401302 :       break;
     699      467955 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     700     2111711 :       pari_APPLY_same(gprec_w(gel(x,i), pr));
     701      703552 :     default: return x;
     702             :   }
     703     2269722 :   return y;
     704             : }
     705             : /* not GC-safe; precision given in word length (including codewords) */
     706             : GEN
     707     5942245 : gprec_wensure(GEN x, long pr)
     708             : {
     709             :   long lx, i;
     710             :   GEN y;
     711             : 
     712     5942245 :   switch(typ(x))
     713             :   {
     714     5226254 :     case t_REAL:
     715     5226254 :       if (signe(x)) return realprec(x) < pr? rtor(x,pr): x;
     716       10340 :       i = -prec2nbits(pr);
     717       10340 :       if (i < expo(x)) return real_0_bit(i);
     718        7920 :       y = cgetr(2); y[1] = x[1]; return y;
     719      201787 :     case t_COMPLEX:
     720      201787 :       y = cgetg(3, t_COMPLEX);
     721      201787 :       gel(y,1) = gprec_wensure(gel(x,1),pr);
     722      201787 :       gel(y,2) = gprec_wensure(gel(x,2),pr);
     723      201787 :       break;
     724       49784 :    case t_POL: case t_SER:
     725       49784 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     726      868336 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wensure(gel(x,i),pr);
     727       49784 :       break;
     728       83522 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     729     1300842 :       pari_APPLY_same(gprec_wensure(gel(x,i), pr));
     730      380898 :     default: return x;
     731             :   }
     732      251571 :   return y;
     733             : }
     734             : 
     735             : /* not GC-safe; precision given in word length (including codewords),
     736             :  * truncate mantissa to precision 'pr' but never increase it */
     737             : GEN
     738     3531897 : gprec_wtrunc(GEN x, long pr)
     739             : {
     740             :   long lx, i;
     741             :   GEN y;
     742             : 
     743     3531897 :   switch(typ(x))
     744             :   {
     745     3088093 :     case t_REAL:
     746     3088093 :       return (signe(x) && realprec(x) > pr)? rtor(x,pr): x;
     747      279679 :     case t_COMPLEX:
     748      279679 :       y = cgetg(3, t_COMPLEX);
     749      279680 :       gel(y,1) = gprec_wtrunc(gel(x,1),pr);
     750      279679 :       gel(y,2) = gprec_wtrunc(gel(x,2),pr);
     751      279679 :       break;
     752        4200 :     case t_POL:
     753             :     case t_SER:
     754        4200 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     755       25970 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     756        4200 :       break;
     757       86399 :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     758      397092 :       pari_APPLY_same(gprec_wtrunc(gel(x,i), pr));
     759       73526 :     default: return x;
     760             :   }
     761      283879 :   return y;
     762             : }
     763             : 
     764             : /********************************************************************/
     765             : /**                                                                **/
     766             : /**                      SERIES TRANSFORMS                         **/
     767             : /**                                                                **/
     768             : /********************************************************************/
     769             : /**                  LAPLACE TRANSFORM (OF A SERIES)               **/
     770             : /********************************************************************/
     771             : static GEN
     772          14 : serlaplace(GEN x)
     773             : {
     774          14 :   long i, l = lg(x), e = valser(x);
     775          14 :   GEN t, y = cgetg(l,t_SER);
     776          14 :   if (e < 0) pari_err_DOMAIN("laplace","valuation","<",gen_0,stoi(e));
     777          14 :   t = mpfact(e); y[1] = x[1];
     778         154 :   for (i=2; i<l; i++)
     779             :   {
     780         140 :     gel(y,i) = gmul(t, gel(x,i));
     781         140 :     e++; t = mului(e,t);
     782             :   }
     783          14 :   return y;
     784             : }
     785             : static GEN
     786          14 : pollaplace(GEN x)
     787             : {
     788          14 :   long i, e = 0, l = lg(x);
     789          14 :   GEN t = gen_1, y = cgetg(l,t_POL);
     790          14 :   y[1] = x[1];
     791          63 :   for (i=2; i<l; i++)
     792             :   {
     793          49 :     gel(y,i) = gmul(t, gel(x,i));
     794          49 :     e++; t = mului(e,t);
     795             :   }
     796          14 :   return y;
     797             : }
     798             : GEN
     799          35 : laplace(GEN x)
     800             : {
     801          35 :   pari_sp av = avma;
     802          35 :   switch(typ(x))
     803             :   {
     804          14 :     case t_POL: x = pollaplace(x); break;
     805          14 :     case t_SER: x = serlaplace(x); break;
     806           7 :     default: if (is_scalar_t(typ(x))) return gcopy(x);
     807           0 :              pari_err_TYPE("laplace",x);
     808             :   }
     809          28 :   return gerepilecopy(av, x);
     810             : }
     811             : 
     812             : /********************************************************************/
     813             : /**              CONVOLUTION PRODUCT (OF TWO SERIES)               **/
     814             : /********************************************************************/
     815             : GEN
     816          14 : convol(GEN x, GEN y)
     817             : {
     818          14 :   long j, lx, ly, ex, ey, vx = varn(x);
     819             :   GEN z;
     820             : 
     821          14 :   if (typ(x) != t_SER) pari_err_TYPE("convol",x);
     822          14 :   if (typ(y) != t_SER) pari_err_TYPE("convol",y);
     823          14 :   if (varn(y) != vx) pari_err_VAR("convol", x,y);
     824          14 :   ex = valser(x);
     825          14 :   ey = valser(y);
     826          14 :   if (ser_isexactzero(x))
     827             :   {
     828           7 :     z = scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), 1);
     829           7 :     setvalser(z, maxss(ex,ey)); return z;
     830             :   }
     831           7 :   lx = lg(x) + ex; x -= ex;
     832           7 :   ly = lg(y) + ey; y -= ey;
     833             :   /* inputs shifted: x[i] and y[i] now correspond to monomials of same degree */
     834           7 :   if (ly < lx) lx = ly; /* min length */
     835           7 :   if (ex < ey) ex = ey; /* max valuation */
     836           7 :   if (lx - ex < 3) return zeroser(vx, lx-2);
     837             : 
     838           7 :   z = cgetg(lx - ex, t_SER);
     839           7 :   z[1] = evalvalser(ex) | evalvarn(vx);
     840         119 :   for (j = ex+2; j<lx; j++) gel(z,j-ex) = gmul(gel(x,j),gel(y,j));
     841           7 :   return normalizeser(z);
     842             : }
     843             : 
     844             : /***********************************************************************/
     845             : /*               OPERATIONS ON DIRICHLET SERIES: *, /                  */
     846             : /* (+, -, scalar multiplication are done on the corresponding vectors) */
     847             : /***********************************************************************/
     848             : static long
     849      869288 : dirval(GEN x)
     850             : {
     851      869288 :   long i = 1, lx = lg(x);
     852      869309 :   while (i < lx && gequal0(gel(x,i))) i++;
     853      869288 :   return i;
     854             : }
     855             : 
     856             : GEN
     857         336 : dirmul(GEN x, GEN y)
     858             : {
     859         336 :   pari_sp av = avma, av2;
     860             :   long nx, ny, nz, dx, dy, i, j, k;
     861             :   GEN z;
     862             : 
     863         336 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirmul",x);
     864         336 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirmul",y);
     865         336 :   dx = dirval(x); nx = lg(x)-1;
     866         336 :   dy = dirval(y); ny = lg(y)-1;
     867         336 :   if (ny-dy < nx-dx) { swap(x,y); lswap(nx,ny); lswap(dx,dy); }
     868         336 :   nz = minss(nx*dy,ny*dx);
     869         336 :   y = RgV_kill0(y);
     870         336 :   av2 = avma;
     871         336 :   z = zerovec(nz);
     872       39095 :   for (j=dx; j<=nx; j++)
     873             :   {
     874       38759 :     GEN c = gel(x,j);
     875       38759 :     if (gequal0(c)) continue;
     876       17031 :     if (gequal1(c))
     877             :     {
     878       94199 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     879       88550 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gel(y,k));
     880             :     }
     881       11382 :     else if (gequalm1(c))
     882             :     {
     883        5649 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     884        4298 :         if (gel(y,k)) gel(z,i) = gsub(gel(z,i),gel(y,k));
     885             :     }
     886             :     else
     887             :     {
     888       46508 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     889       36477 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gmul(c,gel(y,k)));
     890             :     }
     891       17031 :     if (gc_needed(av2,3))
     892             :     {
     893           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirmul, %ld/%ld",j,nx);
     894           0 :       z = gerepilecopy(av2,z);
     895             :     }
     896             :   }
     897         336 :   return gerepilecopy(av,z);
     898             : }
     899             : 
     900             : GEN
     901      434308 : dirdiv(GEN x, GEN y)
     902             : {
     903      434308 :   pari_sp av = avma, av2;
     904             :   long nx,ny,nz, dx,dy, i,j,k;
     905             :   GEN p1;
     906             : 
     907      434308 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirdiv",x);
     908      434308 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirdiv",y);
     909      434308 :   dx = dirval(x); nx = lg(x)-1;
     910      434308 :   dy = dirval(y); ny = lg(y)-1;
     911      434308 :   if (dy != 1 || !ny) pari_err_INV("dirdiv",y);
     912      434308 :   nz = minss(nx,ny*dx);
     913      434308 :   p1 = gel(y,1);
     914      434308 :   if (gequal1(p1)) p1 = NULL; else y = gdiv(y,p1);
     915      434308 :   y = RgV_kill0(y);
     916      434308 :   av2 = avma;
     917      434308 :   x = p1 ? gdiv(x,p1): leafcopy(x);
     918      434315 :   for (j=1; j<dx; j++) gel(x,j) = gen_0;
     919      434308 :   setlg(x,nz+1);
     920   109756234 :   for (j=dx; j<=nz; j++)
     921             :   {
     922   109321926 :     GEN c = gel(x,j);
     923   109321926 :     if (gequal0(c)) continue;
     924    75811155 :     if (gequal1(c))
     925             :     {
     926   133665742 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     927   131901014 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gel(y,k));
     928             :     }
     929    74046427 :     else if (gequalm1(c))
     930             :     {
     931    28792855 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     932    27244343 :         if (gel(y,k)) gel(x,i) = gadd(gel(x,i),gel(y,k));
     933             :     }
     934             :     else
     935             :     {
     936   331181123 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     937   258683208 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gmul(c,gel(y,k)));
     938             :     }
     939    75811155 :     if (gc_needed(av2,3))
     940             :     {
     941           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirdiv, %ld/%ld",j,nz);
     942           0 :       x = gerepilecopy(av2,x);
     943             :     }
     944             :   }
     945      434308 :   return gerepilecopy(av,x);
     946             : }
     947             : 
     948             : /*******************************************************************/
     949             : /**                                                               **/
     950             : /**                       COMBINATORICS                           **/
     951             : /**                                                               **/
     952             : /*******************************************************************/
     953             : /**                      BINOMIAL COEFFICIENTS                    **/
     954             : /*******************************************************************/
     955             : /* Lucas's formula for v_p(\binom{n}{k}), used in the tough case p <= sqrt(n) */
     956             : static long
     957        3206 : binomial_lval(ulong n, ulong k, ulong p)
     958             : {
     959        3206 :   ulong r = 0, e = 0;
     960             :   do
     961             :   {
     962       10290 :     ulong a = n % p, b = k % p + r;
     963       10290 :     n /= p; k /= p;
     964       10290 :     if (a < b) { e++; r = 1; } else r = 0;
     965       10290 :   } while (n);
     966        3206 :   return e;
     967             : }
     968             : GEN
     969       79995 : binomialuu(ulong n, ulong k)
     970             : {
     971       79995 :   pari_sp av = avma;
     972             :   ulong p, nk, sn;
     973             :   long c, l;
     974             :   forprime_t T;
     975             :   GEN v, z;
     976       79995 :   if (k > n) return gen_0;
     977       79988 :   nk = n-k; if (k > nk) lswap(nk, k);
     978       79988 :   if (!k) return gen_1;
     979       78364 :   if (k == 1) return utoipos(n);
     980       73016 :   if (k == 2) return muluu(odd(n)? n: n-1, n>>1);
     981       53045 :   if (k < 1000 || ((double)k/ n) * log((double)n) < 0.5)
     982             :   { /* k "small" */
     983       53031 :     z = diviiexact(mulu_interval(n-k+1, n), mulu_interval(2UL, k));
     984       53031 :     return gerepileuptoint(av, z);
     985             :   }
     986          14 :   sn = usqrt(n);
     987             :   /* use Lucas's formula, k <= n/2 */
     988          14 :   l = minuu(1UL << 20, n); v = cgetg(l+1, t_VECSMALL); c = 1;
     989          14 :   u_forprime_init(&T, nk+1, n);
     990     1553958 :   while ((p = u_forprime_next(&T))) /* all primes n-k < p <= n occur, v_p = 1 */
     991             :   {
     992     1553944 :     if (c == l) { ulong L = l << 1; v = vecsmall_lengthen(v, L); l = L; }
     993     1553944 :     v[c++] = p;
     994             :   }
     995          14 :   u_forprime_init(&T, sn+1, n >> 1);
     996     2437785 :   while ((p = u_forprime_next(&T))) /* p^2 > n, v_p <= 1 */
     997     2437771 :     if (n % p < k % p)
     998             :     {
     999     1428679 :       if (c == l) { ulong L = l << 1; v = vecsmall_lengthen(v, L); l = L; }
    1000     1428679 :       v[c++] = p;
    1001             :     }
    1002          14 :   setlg(v, c); z = zv_prod_Z(v);
    1003          14 :   u_forprime_init(&T, 3, sn);
    1004          14 :   l = minuu(1UL << 20, sn); v = cgetg(l + 1, t_VEC); c = 1;
    1005        3220 :   while ((p = u_forprime_next(&T))) /* p <= sqrt(n) */
    1006             :   {
    1007        3206 :     ulong e = binomial_lval(n, k, p);
    1008        3206 :     if (e)
    1009             :     {
    1010        2541 :       if (c == l) { ulong L = l << 1; v = vec_lengthen(v, L); l = L; }
    1011        2541 :       gel(v, c++) = powuu(p, e);
    1012             :     }
    1013             :   }
    1014          14 :   setlg(v, c); z = mulii(z, ZV_prod(v));
    1015             :   { /* p = 2 */
    1016          14 :     ulong e = hammingl(k);
    1017          14 :     e += (k == nk)? e: hammingl(nk);
    1018          14 :     e -= hammingl(n); if (e) z = shifti(z, e);
    1019             :   }
    1020          14 :   return gerepileuptoint(av, z);
    1021             : }
    1022             : 
    1023             : GEN
    1024      100989 : binomial(GEN n, long k)
    1025             : {
    1026      100989 :   long i, prec, tn = typ(n);
    1027             :   pari_sp av;
    1028             :   GEN y;
    1029             : 
    1030      100989 :   av = avma;
    1031      100989 :   if (tn == t_INT)
    1032             :   {
    1033             :     long sn;
    1034             :     GEN z;
    1035      100814 :     if (k == 0) return gen_1;
    1036       67431 :     sn = signe(n);
    1037       67431 :     if (sn == 0) return gen_0; /* k != 0 */
    1038       67431 :     if (sn > 0)
    1039             :     { /* n > 0 */
    1040       67340 :       if (k < 0) return gen_0;
    1041       67340 :       if (k == 1) return icopy(n);
    1042       41349 :       z = subiu(n, k);
    1043       41349 :       if (cmpiu(z, k) < 0)
    1044             :       {
    1045        1043 :         switch(signe(z))
    1046             :         {
    1047           7 :           case -1: return gc_const(av, gen_0);
    1048          63 :           case 0: return gc_const(av, gen_1);
    1049             :         }
    1050         973 :         k = z[2];
    1051         973 :         if (k == 1) { set_avma(av); return icopy(n); }
    1052             :       }
    1053       40880 :       set_avma(av);
    1054       40880 :       if (lgefint(n) == 3) return binomialuu(n[2],(ulong)k);
    1055             :     }
    1056             :     else
    1057             :     { /* n < 0, k != 0; use Kronenburg's definition */
    1058          91 :       if (k > 0)
    1059          70 :         z = binomial(subsi(k - 1, n), k);
    1060             :       else
    1061             :       {
    1062          21 :         z = subis(n, k); if (signe(z) < 0) return gen_0;
    1063          14 :         n = stoi(-k-1); k = itos(z);
    1064          14 :         z = binomial(n, k);
    1065             :       }
    1066          84 :       if (odd(k)) togglesign_safe(&z);
    1067          84 :       return gerepileuptoint(av, z);
    1068             :     }
    1069             :     /* n >= 0 and huge, k != 0 */
    1070           8 :     if (k < 0) return gen_0;
    1071           8 :     if (k == 1) return icopy(n);
    1072             :     /* k > 1 */
    1073           8 :     y = cgetg(k+1,t_VEC); gel(y,1) = n;
    1074          18 :     for (i = 2; i <= k; i++) gel(y,i) = subiu(n,i-1);
    1075           8 :     y = diviiexact(ZV_prod(y), mpfact(k));
    1076           8 :     return gerepileuptoint(av, y);
    1077             :   }
    1078         175 :   if (is_noncalc_t(tn)) pari_err_TYPE("binomial",n);
    1079         175 :   if (k <= 1)
    1080             :   {
    1081          14 :     if (k < 0) return Rg_get_0(n);
    1082           7 :     if (k == 0) return Rg_get_1(n);
    1083           0 :     return gcopy(n);
    1084             :   }
    1085         161 :   prec = precision(n);
    1086         161 :   if (prec && k > 200 + 0.8*prec2nbits(prec)) {
    1087           7 :     GEN A = mpfactr(k, prec), B = ggamma(gsubgs(n,k-1), prec);
    1088           7 :     return gerepileupto(av, gdiv(ggamma(gaddgs(n,1), prec), gmul(A,B)));
    1089             :   }
    1090             : 
    1091         154 :   y = cgetg(k+1,t_VEC);
    1092       12236 :   for (i=1; i<=k; i++) gel(y,i) = gsubgs(n,i-1);
    1093         154 :   return gerepileupto(av, gdiv(RgV_prod(y), mpfact(k)));
    1094             : }
    1095             : 
    1096             : GEN
    1097         924 : binomial0(GEN x, GEN k)
    1098             : {
    1099         924 :   if (!k)
    1100             :   {
    1101          21 :     if (typ(x) != t_INT || signe(x) < 0) pari_err_TYPE("binomial", x);
    1102           7 :     return vecbinomial(itos(x));
    1103             :   }
    1104         903 :   if (typ(k) != t_INT) pari_err_TYPE("binomial", k);
    1105         896 :   return binomial(x, itos(k));
    1106             : }
    1107             : 
    1108             : /* Assume n >= 0, return bin, bin[k+1] = binomial(n, k) */
    1109             : GEN
    1110      147354 : vecbinomial(long n)
    1111             : {
    1112             :   long d, k;
    1113             :   GEN C;
    1114      147354 :   if (!n) return mkvec(gen_1);
    1115      146997 :   C = cgetg(n+2, t_VEC) + 1; /* C[k] = binomial(n, k) */
    1116      146997 :   gel(C,0) = gen_1;
    1117      146997 :   gel(C,1) = utoipos(n); d = (n + 1) >> 1;
    1118      637391 :   for (k=2; k <= d; k++)
    1119             :   {
    1120      490393 :     pari_sp av = avma;
    1121      490393 :     gel(C,k) = gerepileuptoint(av, diviuexact(mului(n-k+1, gel(C,k-1)), k));
    1122             :   }
    1123      711654 :   for (   ; k <= n; k++) gel(C,k) = gel(C,n-k);
    1124      146998 :   return C - 1;
    1125             : }
    1126             : 
    1127             : /********************************************************************/
    1128             : /**                  STIRLING NUMBERS                              **/
    1129             : /********************************************************************/
    1130             : /* Stirling number of the 2nd kind. The number of ways of partitioning
    1131             :    a set of n elements into m nonempty subsets. */
    1132             : GEN
    1133        1694 : stirling2(ulong n, ulong m)
    1134             : {
    1135        1694 :   pari_sp av = avma;
    1136             :   GEN s, bmk;
    1137             :   ulong k;
    1138        1694 :   if (n==0) return (m == 0)? gen_1: gen_0;
    1139        1694 :   if (m > n || m == 0) return gen_0;
    1140        1694 :   if (m==n) return gen_1;
    1141             :   /* k = 0 */
    1142        1694 :   bmk = gen_1; s  = powuu(m, n);
    1143       20314 :   for (k = 1; k <= ((m-1)>>1); ++k)
    1144             :   { /* bmk = binomial(m, k) */
    1145             :     GEN c, kn, mkn;
    1146       18620 :     bmk = diviuexact(mului(m-k+1, bmk), k);
    1147       18620 :     kn  = powuu(k, n); mkn = powuu(m-k, n);
    1148       18620 :     c = odd(m)? subii(mkn,kn): addii(mkn,kn);
    1149       18620 :     c = mulii(bmk, c);
    1150       18620 :     s = odd(k)? subii(s, c): addii(s, c);
    1151       18620 :     if (gc_needed(av,2))
    1152             :     {
    1153           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"stirling2");
    1154           0 :       gerepileall(av, 2, &s, &bmk);
    1155             :     }
    1156             :   }
    1157             :   /* k = m/2 */
    1158        1694 :   if (!odd(m))
    1159             :   {
    1160             :     GEN c;
    1161         805 :     bmk = diviuexact(mului(k+1, bmk), k);
    1162         805 :     c = mulii(bmk, powuu(k,n));
    1163         805 :     s = odd(k)? subii(s, c): addii(s, c);
    1164             :   }
    1165        1694 :   return gerepileuptoint(av, diviiexact(s, mpfact(m)));
    1166             : }
    1167             : 
    1168             : /* Stirling number of the first kind. Up to the sign, the number of
    1169             :    permutations of n symbols which have exactly m cycles. */
    1170             : GEN
    1171         154 : stirling1(ulong n, ulong m)
    1172             : {
    1173         154 :   pari_sp ltop=avma;
    1174             :   ulong k;
    1175             :   GEN s, t;
    1176         154 :   if (n < m) return gen_0;
    1177         154 :   else if (n==m) return gen_1;
    1178             :   /* t = binomial(n-1+k, m-1) * binomial(2n-m, n-m-k) */
    1179             :   /* k = n-m > 0 */
    1180         154 :   t = binomialuu(2*n-m-1, m-1);
    1181         154 :   s = mulii(t, stirling2(2*(n-m), n-m));
    1182         154 :   if (odd(n-m)) togglesign(s);
    1183        1547 :   for (k = n-m-1; k > 0; --k)
    1184             :   {
    1185             :     GEN c;
    1186        1393 :     t = diviuuexact(muluui(n-m+k+1, n+k+1, t), n+k, n-m-k);
    1187        1393 :     c = mulii(t, stirling2(n-m+k, k));
    1188        1393 :     s = odd(k)? subii(s, c): addii(s, c);
    1189        1393 :     if ((k & 0x1f) == 0) {
    1190          21 :       t = gerepileuptoint(ltop, t);
    1191          21 :       s = gerepileuptoint(avma, s);
    1192             :     }
    1193             :   }
    1194         154 :   return gerepileuptoint(ltop, s);
    1195             : }
    1196             : 
    1197             : GEN
    1198         301 : stirling(long n, long m, long flag)
    1199             : {
    1200         301 :   if (n < 0) pari_err_DOMAIN("stirling", "n", "<", gen_0, stoi(n));
    1201         301 :   if (m < 0) pari_err_DOMAIN("stirling", "m", "<", gen_0, stoi(m));
    1202         301 :   switch (flag)
    1203             :   {
    1204         154 :     case 1: return stirling1((ulong)n,(ulong)m);
    1205         147 :     case 2: return stirling2((ulong)n,(ulong)m);
    1206           0 :     default: pari_err_FLAG("stirling");
    1207             :   }
    1208             :   return NULL; /*LCOV_EXCL_LINE*/
    1209             : }
    1210             : 
    1211             : /*******************************************************************/
    1212             : /**                                                               **/
    1213             : /**                     RECIPROCAL POLYNOMIAL                     **/
    1214             : /**                                                               **/
    1215             : /*******************************************************************/
    1216             : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1217             : GEN
    1218         161 : polrecip(GEN x)
    1219             : {
    1220         161 :   long tx = typ(x);
    1221         161 :   if (is_scalar_t(tx)) return gcopy(x);
    1222         154 :   if (tx != t_POL) pari_err_TYPE("polrecip",x);
    1223         154 :   return RgX_recip(x);
    1224             : }
    1225             : 
    1226             : /********************************************************************/
    1227             : /**                                                                **/
    1228             : /**                  POLYNOMIAL INTERPOLATION                      **/
    1229             : /**                                                                **/
    1230             : /********************************************************************/
    1231             : /* given complex roots L[i], i <= n of some monic T in C[X], return
    1232             :  * the T'(L[i]), computed stably via products of differences */
    1233             : GEN
    1234       83788 : vandermondeinverseinit(GEN L)
    1235             : {
    1236       83788 :   long i, j, l = lg(L);
    1237       83788 :   GEN V = cgetg(l, t_VEC);
    1238      470805 :   for (i = 1; i < l; i++)
    1239             :   {
    1240      387017 :     pari_sp av = avma;
    1241      387017 :     GEN W = cgetg(l-1,t_VEC);
    1242      387015 :     long k = 1;
    1243     4229724 :     for (j = 1; j < l; j++)
    1244     3842870 :       if (i != j) gel(W, k++) = gsub(gel(L,i), gel(L,j));
    1245      386854 :     gel(V,i) = gerepileupto(av, RgV_prod(W));
    1246             :   }
    1247       83788 :   return V;
    1248             : }
    1249             : 
    1250             : /* Compute the inverse of the van der Monde matrix of T multiplied by den */
    1251             : GEN
    1252       52911 : vandermondeinverse(GEN L, GEN T, GEN den, GEN V)
    1253             : {
    1254       52911 :   pari_sp av = avma;
    1255       52911 :   long i, n = lg(L)-1;
    1256       52911 :   GEN M = cgetg(n+1, t_MAT);
    1257             : 
    1258       52911 :   if (!V) V = vandermondeinverseinit(L);
    1259       52911 :   if (den && equali1(den)) den = NULL;
    1260      290028 :   for (i = 1; i <= n; i++)
    1261             :   {
    1262      474226 :     GEN d = gel(V,i), P = RgX_Rg_mul(RgX_div_by_X_x(T, gel(L,i), NULL),
    1263      237115 :                                      den? gdiv(den,d): ginv(d));
    1264      237112 :     gel(M,i) = RgX_to_RgC(P, n);
    1265             :   }
    1266       52913 :   return gerepilecopy(av, M);
    1267             : }
    1268             : 
    1269             : static GEN
    1270         224 : RgV_polint_fast(GEN X, GEN Y, long v)
    1271             : {
    1272             :   GEN p, pol;
    1273             :   long t, pa;
    1274         224 :   if (X) t = RgV_type2(X,Y, &p, &pol, &pa);
    1275          21 :   else   t = Rg_type(Y, &p, &pol, &pa);
    1276         224 :   if (t != t_INTMOD) return NULL;
    1277           7 :   Y = RgC_to_FpC(Y, p);
    1278           7 :   X = X? RgC_to_FpC(X, p): identity_ZV(lg(Y)-1);
    1279           7 :   return FpX_to_mod(FpV_polint(X, Y, p, v), p);
    1280             : }
    1281             : /* allow X = NULL for [1,...,n] */
    1282             : GEN
    1283         224 : RgV_polint(GEN X, GEN Y, long v)
    1284             : {
    1285         224 :   pari_sp av0 = avma, av;
    1286         224 :   GEN Q, L, P = NULL;
    1287         224 :   long i, l = lg(Y);
    1288         224 :   if ((Q = RgV_polint_fast(X,Y,v))) return Q;
    1289         217 :   if (!X) X = identity_ZV(l-1);
    1290         217 :   L = vandermondeinverseinit(X);
    1291         217 :   Q = roots_to_pol(X, v); av = avma;
    1292         553 :   for (i=1; i<l; i++)
    1293             :   {
    1294             :     GEN T, dP;
    1295         336 :     if (gequal0(gel(Y,i))) continue;
    1296         238 :     T = RgX_div_by_X_x(Q, gel(X,i), NULL);
    1297         238 :     dP = RgX_Rg_mul(T, gdiv(gel(Y,i), gel(L,i)));
    1298         238 :     P = P? RgX_add(P, dP): dP;
    1299         238 :     if (gc_needed(av,2))
    1300             :     {
    1301           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"RgV_polint i = %ld/%ld", i, l-1);
    1302           0 :       P = gerepileupto(av, P);
    1303             :     }
    1304             :   }
    1305         217 :   if (!P) { set_avma(av); return zeropol(v); }
    1306         147 :   return gerepileupto(av0, P);
    1307             : }
    1308             : static int
    1309       17357 : inC(GEN x)
    1310             : {
    1311       17357 :   switch(typ(x)) {
    1312        1365 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD: return 1;
    1313       15992 :     default: return 0;
    1314             :   }
    1315             : }
    1316             : static long
    1317       16188 : check_dy(GEN X, GEN x, long n)
    1318             : {
    1319       16188 :   GEN D = NULL;
    1320       16188 :   long i, ns = 0;
    1321       16188 :   if (!inC(x)) return -1;
    1322        1176 :   for (i = 0; i < n; i++)
    1323             :   {
    1324         966 :     GEN t = gsub(x, gel(X,i));
    1325         966 :     if (!inC(t)) return -1;
    1326         952 :     t = gabs(t, DEFAULTPREC);
    1327         952 :     if (!D || gcmp(t,D) < 0) { ns = i; D = t; }
    1328             :   }
    1329             :   /* X[ns] is closest to x */
    1330         210 :   return ns;
    1331             : }
    1332             : /* X,Y are "spec" GEN vectors with n > 0 components ( at X[0], ... X[n-1] ) */
    1333             : GEN
    1334       16223 : polintspec(GEN X, GEN Y, GEN x, long n, long *pe)
    1335             : {
    1336             :   long i, m, ns;
    1337       16223 :   pari_sp av = avma, av2;
    1338       16223 :   GEN y, c, d, dy = NULL; /* gcc -Wall */
    1339             : 
    1340       16223 :   if (pe) *pe = -HIGHEXPOBIT;
    1341       16223 :   if (n == 1) return gmul(gel(Y,0), Rg_get_1(x));
    1342       16188 :   if (!X) X = identity_ZV(n) + 1;
    1343       16188 :   av2 = avma;
    1344       16188 :   ns = check_dy(X, x, n); if (ns < 0) { pe = NULL; ns = 0; }
    1345       16188 :   c = cgetg(n+1, t_VEC);
    1346       81031 :   d = cgetg(n+1, t_VEC); for (i=0; i<n; i++) gel(c,i+1) = gel(d,i+1) = gel(Y,i);
    1347       16188 :   y = gel(d,ns+1);
    1348             :   /* divided differences */
    1349       64836 :   for (m = 1; m < n; m++)
    1350             :   {
    1351      146238 :     for (i = 0; i < n-m; i++)
    1352             :     {
    1353       97590 :       GEN ho = gsub(gel(X,i),x), hp = gsub(gel(X,i+m),x), den = gsub(ho,hp);
    1354       97590 :       if (gequal0(den))
    1355             :       {
    1356           7 :         char *x1 = stack_sprintf("X[%ld]", i+1);
    1357           7 :         char *x2 = stack_sprintf("X[%ld]", i+m+1);
    1358           7 :         pari_err_DOMAIN("polinterpolate",x1,"=",strtoGENstr(x2), X);
    1359             :       }
    1360       97583 :       den = gdiv(gsub(gel(c,i+2),gel(d,i+1)), den);
    1361       97583 :       gel(c,i+1) = gmul(ho,den);
    1362       97583 :       gel(d,i+1) = gmul(hp,den);
    1363             :     }
    1364       48648 :     dy = (2*ns < n-m)? gel(c,ns+1): gel(d,ns--);
    1365       48648 :     y = gadd(y,dy);
    1366       48648 :     if (gc_needed(av2,2))
    1367             :     {
    1368           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"polint, %ld/%ld",m,n-1);
    1369           0 :       gerepileall(av2, 4, &y, &c, &d, &dy);
    1370             :     }
    1371             :   }
    1372       16181 :   if (pe && inC(dy)) *pe = gexpo(dy);
    1373       16181 :   return gerepileupto(av, y);
    1374             : }
    1375             : 
    1376             : GEN
    1377         329 : polint_i(GEN X, GEN Y, GEN t, long *pe)
    1378             : {
    1379         329 :   long lx = lg(X), vt;
    1380             : 
    1381         329 :   if (! is_vec_t(typ(X))) pari_err_TYPE("polinterpolate",X);
    1382         329 :   if (Y)
    1383             :   {
    1384         301 :     if (! is_vec_t(typ(Y))) pari_err_TYPE("polinterpolate",Y);
    1385         301 :     if (lx != lg(Y)) pari_err_DIM("polinterpolate");
    1386             :   }
    1387             :   else
    1388             :   {
    1389          28 :     Y = X;
    1390          28 :     X = NULL;
    1391             :   }
    1392         329 :   if (pe) *pe = -HIGHEXPOBIT;
    1393         329 :   vt = t? gvar(t): 0;
    1394         329 :   if (vt != NO_VARIABLE)
    1395             :   { /* formal interpolation */
    1396             :     pari_sp av;
    1397         224 :     long v0, vY = gvar(Y);
    1398             :     GEN P;
    1399         224 :     if (X) vY = varnmax(vY, gvar(X));
    1400             :     /* shortcut */
    1401         224 :     if (varncmp(vY, vt) > 0 && (!t || gequalX(t))) return RgV_polint(X, Y, vt);
    1402          84 :     av = avma;
    1403             :     /* first interpolate in high priority variable, then substitute t */
    1404          84 :     v0 = fetch_var_higher();
    1405          84 :     P = RgV_polint(X, Y, v0);
    1406          84 :     P = gsubst(P, v0, t? t: pol_x(0));
    1407          84 :     (void)delete_var();
    1408          84 :     return gerepileupto(av, P);
    1409             :   }
    1410             :   /* numerical interpolation */
    1411         105 :   if (lx == 1) return Rg_get_0(t);
    1412          91 :   return polintspec(X? X+1: NULL,Y+1,t,lx-1, pe);
    1413             : }
    1414             : GEN
    1415         329 : polint(GEN X, GEN Y, GEN t, GEN *pe)
    1416             : {
    1417             :   long e;
    1418         329 :   GEN p = polint_i(X, Y, t, &e);
    1419         322 :   if (pe) *pe = stoi(e);
    1420         322 :   return p;
    1421             : }
    1422             : 
    1423             : /********************************************************************/
    1424             : /**                                                                **/
    1425             : /**                       MODREVERSE                               **/
    1426             : /**                                                                **/
    1427             : /********************************************************************/
    1428             : static void
    1429           7 : err_reverse(GEN x, GEN T)
    1430             : {
    1431           7 :   pari_err_DOMAIN("modreverse","deg(minpoly(z))", "<", stoi(degpol(T)),
    1432             :                   mkpolmod(x,T));
    1433           0 : }
    1434             : 
    1435             : /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
    1436             : GEN
    1437         175 : RgXQ_reverse(GEN a, GEN T)
    1438             : {
    1439         175 :   pari_sp av = avma;
    1440         175 :   long n = degpol(T);
    1441             :   GEN y;
    1442             : 
    1443         175 :   if (n <= 1) {
    1444           7 :     if (n <= 0) return gcopy(a);
    1445           7 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1446             :   }
    1447         168 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1448         168 :   y = RgXV_to_RgM(RgXQ_powers(a,n-1,T), n);
    1449         168 :   y = RgM_solve(y, col_ei(n, 2));
    1450         168 :   if (!y) err_reverse(a,T);
    1451         161 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1452             : }
    1453             : GEN
    1454        5931 : QXQ_reverse(GEN a, GEN T)
    1455             : {
    1456        5931 :   pari_sp av = avma;
    1457        5931 :   long n = degpol(T);
    1458             :   GEN y;
    1459             : 
    1460        5931 :   if (n <= 1) {
    1461          14 :     if (n <= 0) return gcopy(a);
    1462          14 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1463             :   }
    1464        5917 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1465        5917 :   if (gequalX(a)) return gcopy(a);
    1466        5695 :   y = RgXV_to_RgM(QXQ_powers(a,n-1,T), n);
    1467        5695 :   y = QM_gauss(y, col_ei(n, 2));
    1468        5695 :   if (!y) err_reverse(a,T);
    1469        5695 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1470             : }
    1471             : 
    1472             : GEN
    1473          28 : modreverse(GEN x)
    1474             : {
    1475             :   long v, n;
    1476             :   GEN T, a;
    1477             : 
    1478          28 :   if (typ(x)!=t_POLMOD) pari_err_TYPE("modreverse",x);
    1479          28 :   T = gel(x,1); n = degpol(T); if (n <= 0) return gcopy(x);
    1480          21 :   a = gel(x,2);
    1481          21 :   v = varn(T);
    1482          21 :   retmkpolmod(RgXQ_reverse(a, T),
    1483             :               (n==1)? gsub(pol_x(v), a): RgXQ_charpoly(a, T, v));
    1484             : }
    1485             : 
    1486             : /********************************************************************/
    1487             : /**                                                                **/
    1488             : /**                          MERGESORT                             **/
    1489             : /**                                                                **/
    1490             : /********************************************************************/
    1491             : static int
    1492          77 : cmp_small(GEN x, GEN y) {
    1493          77 :   long a = (long)x, b = (long)y;
    1494          77 :   return a>b? 1: (a<b? -1: 0);
    1495             : }
    1496             : 
    1497             : static int
    1498      295015 : veccmp(void *data, GEN x, GEN y)
    1499             : {
    1500      295015 :   GEN k = (GEN)data;
    1501      295015 :   long i, s, lk = lg(k), lx = minss(lg(x), lg(y));
    1502             : 
    1503      295015 :   if (!is_vec_t(typ(x))) pari_err_TYPE("lexicographic vecsort",x);
    1504      295015 :   if (!is_vec_t(typ(y))) pari_err_TYPE("lexicographic vecsort",y);
    1505      306684 :   for (i=1; i<lk; i++)
    1506             :   {
    1507      295043 :     long c = k[i];
    1508      295043 :     if (c >= lx)
    1509          14 :       pari_err_TYPE("lexicographic vecsort, index too large", stoi(c));
    1510      295029 :     s = lexcmp(gel(x,c), gel(y,c));
    1511      295029 :     if (s) return s;
    1512             :   }
    1513       11641 :   return 0;
    1514             : }
    1515             : 
    1516             : /* return permutation sorting v[1..n], removing duplicates. Assume n > 0 */
    1517             : static GEN
    1518     2011356 : gen_sortspec_uniq(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1519             : {
    1520             :   pari_sp av;
    1521             :   long NX, nx, ny, m, ix, iy, i;
    1522             :   GEN x, y, w, W;
    1523             :   int s;
    1524     2011356 :   switch(n)
    1525             :   {
    1526       84006 :     case 1: return mkvecsmall(1);
    1527      816026 :     case 2:
    1528      816026 :       s = cmp(E,gel(v,1),gel(v,2));
    1529      816031 :       if      (s < 0) return mkvecsmall2(1,2);
    1530      375907 :       else if (s > 0) return mkvecsmall2(2,1);
    1531        6244 :       return mkvecsmall(1);
    1532      268175 :     case 3:
    1533      268175 :       s = cmp(E,gel(v,1),gel(v,2));
    1534      268175 :       if (s < 0) {
    1535      166673 :         s = cmp(E,gel(v,2),gel(v,3));
    1536      166673 :         if (s < 0) return mkvecsmall3(1,2,3);
    1537       65688 :         else if (s == 0) return mkvecsmall2(1,2);
    1538       65057 :         s = cmp(E,gel(v,1),gel(v,3));
    1539       65057 :         if      (s < 0) return mkvecsmall3(1,3,2);
    1540       33011 :         else if (s > 0) return mkvecsmall3(3,1,2);
    1541        2212 :         return mkvecsmall2(1,2);
    1542      101502 :       } else if (s > 0) {
    1543       98394 :         s = cmp(E,gel(v,1),gel(v,3));
    1544       98394 :         if (s < 0) return mkvecsmall3(2,1,3);
    1545       66612 :         else if (s == 0) return mkvecsmall2(2,1);
    1546       65227 :         s = cmp(E,gel(v,2),gel(v,3));
    1547       65227 :         if (s < 0) return mkvecsmall3(2,3,1);
    1548       32067 :         else if (s > 0) return mkvecsmall3(3,2,1);
    1549         721 :         return mkvecsmall2(2,1);
    1550             :       } else {
    1551        3108 :         s = cmp(E,gel(v,1),gel(v,3));
    1552        3108 :         if (s < 0) return mkvecsmall2(1,3);
    1553        1603 :         else if (s == 0) return mkvecsmall(1);
    1554         763 :         return mkvecsmall2(3,1);
    1555             :       }
    1556             :   }
    1557      843149 :   NX = nx = n>>1; ny = n-nx;
    1558      843149 :   av = avma;
    1559      843149 :   x = gen_sortspec_uniq(v,   nx,E,cmp); nx = lg(x)-1;
    1560      843195 :   y = gen_sortspec_uniq(v+NX,ny,E,cmp); ny = lg(y)-1;
    1561      843200 :   w = cgetg(n+1, t_VECSMALL);
    1562      843195 :   m = ix = iy = 1;
    1563     9992232 :   while (ix<=nx && iy<=ny)
    1564             :   {
    1565     9149037 :     s = cmp(E, gel(v,x[ix]), gel(v,y[iy]+NX));
    1566     9149037 :     if (s < 0)
    1567     4226929 :       w[m++] = x[ix++];
    1568     4922108 :     else if (s > 0)
    1569     3841753 :       w[m++] = y[iy++]+NX;
    1570             :     else {
    1571     1080355 :       w[m++] = x[ix++];
    1572     1080355 :       iy++;
    1573             :     }
    1574             :   }
    1575     1367320 :   while (ix<=nx) w[m++] = x[ix++];
    1576     2000305 :   while (iy<=ny) w[m++] = y[iy++]+NX;
    1577      843195 :   set_avma(av);
    1578      843196 :   W = cgetg(m, t_VECSMALL);
    1579    11673477 :   for (i = 1; i < m; i++) W[i] = w[i];
    1580      843195 :   return W;
    1581             : }
    1582             : 
    1583             : /* return permutation sorting v[1..n]. Assume n > 0 */
    1584             : static GEN
    1585   189967034 : gen_sortspec(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1586             : {
    1587             :   long nx, ny, m, ix, iy;
    1588             :   GEN x, y, w;
    1589   189967034 :   switch(n)
    1590             :   {
    1591     5450231 :     case 1:
    1592     5450231 :       (void)cmp(E,gel(v,1),gel(v,1)); /* check for type error */
    1593     5450316 :       return mkvecsmall(1);
    1594    78679996 :     case 2:
    1595   135122580 :       return cmp(E,gel(v,1),gel(v,2)) <= 0? mkvecsmall2(1,2)
    1596   135122447 :                                           : mkvecsmall2(2,1);
    1597    37089727 :     case 3:
    1598    37089727 :       if (cmp(E,gel(v,1),gel(v,2)) <= 0) {
    1599    27034977 :         if (cmp(E,gel(v,2),gel(v,3)) <= 0) return mkvecsmall3(1,2,3);
    1600    11902040 :         return (cmp(E,gel(v,1),gel(v,3)) <= 0)? mkvecsmall3(1,3,2)
    1601    11902047 :                                               : mkvecsmall3(3,1,2);
    1602             :       } else {
    1603    10054762 :         if (cmp(E,gel(v,1),gel(v,3)) <= 0) return mkvecsmall3(2,1,3);
    1604    10545390 :         return (cmp(E,gel(v,2),gel(v,3)) <= 0)? mkvecsmall3(2,3,1)
    1605    10545391 :                                               : mkvecsmall3(3,2,1);
    1606             :       }
    1607             :   }
    1608    68747080 :   nx = n>>1; ny = n-nx;
    1609    68747080 :   w = cgetg(n+1,t_VECSMALL);
    1610    68749645 :   x = gen_sortspec(v,   nx,E,cmp);
    1611    68749621 :   y = gen_sortspec(v+nx,ny,E,cmp);
    1612    68749665 :   m = ix = iy = 1;
    1613   459873796 :   while (ix<=nx && iy<=ny)
    1614   391124170 :     if (cmp(E, gel(v,x[ix]), gel(v,y[iy]+nx))<=0)
    1615   217058228 :       w[m++] = x[ix++];
    1616             :     else
    1617   174065903 :       w[m++] = y[iy++]+nx;
    1618   104646492 :   while (ix<=nx) w[m++] = x[ix++];
    1619   174980842 :   while (iy<=ny) w[m++] = y[iy++]+nx;
    1620    68749626 :   set_avma((pari_sp)w); return w;
    1621             : }
    1622             : 
    1623             : static void
    1624    44823002 : init_sort(GEN *x, long *tx, long *lx)
    1625             : {
    1626    44823002 :   *tx = typ(*x);
    1627    44823002 :   if (*tx == t_LIST)
    1628             :   {
    1629          35 :     if (list_typ(*x)!=t_LIST_RAW) pari_err_TYPE("sort",*x);
    1630          35 :     *x = list_data(*x);
    1631          35 :     *lx = *x? lg(*x): 1;
    1632             :   } else {
    1633    44822967 :     if (!is_matvec_t(*tx) && *tx != t_VECSMALL) pari_err_TYPE("gen_sort",*x);
    1634    44822967 :     *lx = lg(*x);
    1635             :   }
    1636    44823002 : }
    1637             : 
    1638             : /* (x o y)[1..lx-1], destroy y */
    1639             : INLINE GEN
    1640     2778404 : sort_extract(GEN x, GEN y, long tx, long lx)
    1641             : {
    1642             :   long i;
    1643     2778404 :   switch(tx)
    1644             :   {
    1645           7 :     case t_VECSMALL:
    1646          35 :       for (i=1; i<lx; i++) y[i] = x[y[i]];
    1647           7 :       break;
    1648           7 :     case t_LIST:
    1649           7 :       settyp(y,t_VEC);
    1650          35 :       for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1651           7 :       return gtolist(y);
    1652     2778390 :     default:
    1653     2778390 :       settyp(y,tx);
    1654     8552066 :       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,y[i]));
    1655             :   }
    1656     2778432 :   return y;
    1657             : }
    1658             : 
    1659             : static GEN
    1660     1058588 : triv_sort(long tx) { return tx == t_LIST? mklist(): cgetg(1, tx); }
    1661             : /* Sort the vector x, using cmp to compare entries. */
    1662             : GEN
    1663      276388 : gen_sort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1664             : {
    1665             :   long tx, lx;
    1666             :   GEN y;
    1667             : 
    1668      276388 :   init_sort(&x, &tx, &lx);
    1669      276388 :   if (lx==1) return triv_sort(tx);
    1670      273924 :   y = gen_sortspec_uniq(x,lx-1,E,cmp);
    1671      273924 :   return sort_extract(x, y, tx, lg(y)); /* lg(y) <= lx */
    1672             : }
    1673             : /* Sort the vector x, using cmp to compare entries. */
    1674             : GEN
    1675     3560622 : gen_sort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1676             : {
    1677             :   long tx, lx;
    1678             :   GEN y;
    1679             : 
    1680     3560622 :   init_sort(&x, &tx, &lx);
    1681     3560621 :   if (lx==1) return triv_sort(tx);
    1682     2504497 :   y = gen_sortspec(x,lx-1,E,cmp);
    1683     2504486 :   return sort_extract(x, y, tx, lx);
    1684             : }
    1685             : /* indirect sort: return the permutation that would sort x */
    1686             : GEN
    1687       58905 : gen_indexsort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1688             : {
    1689             :   long tx, lx;
    1690       58905 :   init_sort(&x, &tx, &lx);
    1691       58905 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1692       51106 :   return gen_sortspec_uniq(x,lx-1,E,cmp);
    1693             : }
    1694             : /* indirect sort: return the permutation that would sort x */
    1695             : GEN
    1696      813050 : gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1697             : {
    1698             :   long tx, lx;
    1699      813050 :   init_sort(&x, &tx, &lx);
    1700      813050 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1701      812735 :   return gen_sortspec(x,lx-1,E,cmp);
    1702             : }
    1703             : 
    1704             : /* Sort the vector x in place, using cmp to compare entries */
    1705             : void
    1706    39716559 : gen_sort_inplace(GEN x, void *E, int (*cmp)(void*,GEN,GEN), GEN *perm)
    1707             : {
    1708             :   long tx, lx, i;
    1709    39716559 :   pari_sp av = avma;
    1710             :   GEN y;
    1711             : 
    1712    39716559 :   init_sort(&x, &tx, &lx);
    1713    39716560 :   if (lx<=2)
    1714             :   {
    1715      552112 :     if (perm) *perm = lx == 1? cgetg(1, t_VECSMALL): mkvecsmall(1);
    1716      552112 :     return;
    1717             :   }
    1718    39164448 :   y = gen_sortspec(x,lx-1, E, cmp);
    1719    39164442 :   if (perm)
    1720             :   {
    1721       10997 :     GEN z = new_chunk(lx);
    1722      108262 :     for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1723      108262 :     for (i=1; i<lx; i++) gel(x,i) = gel(z,i);
    1724       10997 :     *perm = y;
    1725       10997 :     set_avma((pari_sp)y);
    1726             :   } else {
    1727   281998340 :     for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1728   281998343 :     for (i=1; i<lx; i++) gel(x,i) = gel(y,i);
    1729    39153445 :     set_avma(av);
    1730             :   }
    1731             : }
    1732             : GEN
    1733      397516 : gen_sort_shallow(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1734             : {
    1735             :   long tx, lx, i;
    1736             :   pari_sp av;
    1737             :   GEN y, z;
    1738             : 
    1739      397516 :   init_sort(&x, &tx, &lx);
    1740      397516 :   if (lx<=2) return x;
    1741      236166 :   z = cgetg(lx, tx); av = avma;
    1742      236166 :   y = gen_sortspec(x,lx-1, E, cmp);
    1743     1226582 :   for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1744      236166 :   return gc_const(av, z);
    1745             : }
    1746             : 
    1747             : static int
    1748        7889 : closurecmp(void *data, GEN x, GEN y)
    1749             : {
    1750        7889 :   pari_sp av = avma;
    1751        7889 :   long s = gsigne(closure_callgen2((GEN)data, x,y));
    1752        7889 :   set_avma(av); return s;
    1753             : }
    1754             : static void
    1755         133 : check_positive_entries(GEN k)
    1756             : {
    1757         133 :   long i, l = lg(k);
    1758         301 :   for (i=1; i<l; i++)
    1759         168 :     if (k[i] <= 0) pari_err_DOMAIN("sort_function", "index", "<", gen_0, stoi(k[i]));
    1760         133 : }
    1761             : 
    1762             : typedef int (*CMP_FUN)(void*,GEN,GEN);
    1763             : /* return NULL if t_CLOSURE k is a "key" (arity 1) and not a sorting func */
    1764             : static CMP_FUN
    1765      126833 : sort_function(void **E, GEN x, GEN k)
    1766             : {
    1767      126833 :   int (*cmp)(GEN,GEN) = &lexcmp;
    1768      126833 :   long tx = typ(x);
    1769      126833 :   if (!k)
    1770             :   {
    1771      126154 :     *E = (void*)((typ(x) == t_VECSMALL)? cmp_small: cmp);
    1772      126154 :     return &cmp_nodata;
    1773             :   }
    1774         679 :   if (tx == t_VECSMALL) pari_err_TYPE("sort_function", x);
    1775         665 :   switch(typ(k))
    1776             :   {
    1777          98 :     case t_INT: k = mkvecsmall(itos(k));  break;
    1778          35 :     case t_VEC: case t_COL: k = ZV_to_zv(k); break;
    1779           0 :     case t_VECSMALL: break;
    1780         532 :     case t_CLOSURE:
    1781         532 :      if (closure_is_variadic(k))
    1782           0 :        pari_err_TYPE("sort_function, variadic cmpf",k);
    1783         532 :      *E = (void*)k;
    1784         532 :      switch(closure_arity(k))
    1785             :      {
    1786          35 :        case 1: return NULL; /* wrt key */
    1787         497 :        case 2: return &closurecmp;
    1788           0 :        default: pari_err_TYPE("sort_function, cmpf arity != 1, 2",k);
    1789             :      }
    1790           0 :     default: pari_err_TYPE("sort_function",k);
    1791             :   }
    1792         133 :   check_positive_entries(k);
    1793         133 :   *E = (void*)k; return &veccmp;
    1794             : }
    1795             : 
    1796             : #define cmp_IND 1
    1797             : #define cmp_LEX 2 /* FIXME: backward compatibility, ignored */
    1798             : #define cmp_REV 4
    1799             : #define cmp_UNIQ 8
    1800             : GEN
    1801         728 : vecsort0(GEN x, GEN k, long flag)
    1802             : {
    1803             :   void *E;
    1804         728 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
    1805             : 
    1806         721 :   if (flag < 0 || flag > (cmp_REV|cmp_LEX|cmp_IND|cmp_UNIQ))
    1807           0 :     pari_err_FLAG("vecsort");
    1808         721 :   if (!CMP)
    1809             :   { /* wrt key: precompute all values, O(n) calls instead of O(n log n) */
    1810          28 :     pari_sp av = avma;
    1811             :     GEN v, y;
    1812             :     long i, tx, lx;
    1813          28 :     init_sort(&x, &tx, &lx);
    1814          28 :     if (lx == 1) return flag&cmp_IND? cgetg(1,t_VECSMALL): triv_sort(tx);
    1815          28 :     v = cgetg(lx, t_VEC);
    1816         140 :     for (i = 1; i < lx; i++) gel(v,i) = closure_callgen1(k, gel(x,i));
    1817          28 :     y = vecsort0(v, NULL, flag | cmp_IND);
    1818          28 :     y = flag&cmp_IND? y: sort_extract(x, y, tx, lg(y));
    1819          28 :     return gerepileupto(av, y);
    1820             :   }
    1821         693 :   if (flag&cmp_UNIQ)
    1822          35 :     x = flag&cmp_IND? gen_indexsort_uniq(x, E, CMP): gen_sort_uniq(x, E, CMP);
    1823             :   else
    1824         658 :     x = flag&cmp_IND? gen_indexsort(x, E, CMP): gen_sort(x, E, CMP);
    1825         679 :   if (flag & cmp_REV)
    1826             :   { /* reverse order */
    1827          35 :     GEN y = x;
    1828          35 :     if (typ(x)==t_LIST) { y = list_data(x); if (!y) return x; }
    1829          28 :     vecreverse_inplace(y);
    1830             :   }
    1831         672 :   return x;
    1832             : }
    1833             : 
    1834             : GEN
    1835      202176 : indexsort(GEN x) { return gen_indexsort(x, (void*)&gcmp, cmp_nodata); }
    1836             : GEN
    1837           0 : indexlexsort(GEN x) { return gen_indexsort(x, (void*)&lexcmp, cmp_nodata); }
    1838             : GEN
    1839          42 : indexvecsort(GEN x, GEN k)
    1840             : {
    1841          42 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1842          42 :   return gen_indexsort(x, (void*)k, &veccmp);
    1843             : }
    1844             : 
    1845             : GEN
    1846      912652 : sort(GEN x) { return gen_sort(x, (void*)gcmp, cmp_nodata); }
    1847             : GEN
    1848           0 : lexsort(GEN x) { return gen_sort(x, (void*)lexcmp, cmp_nodata); }
    1849             : GEN
    1850        2954 : vecsort(GEN x, GEN k)
    1851             : {
    1852        2954 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1853        2954 :   return gen_sort(x, (void*)k, &veccmp);
    1854             : }
    1855             : /* adapted from gen_search; don't export: keys of T[i] should be precomputed */
    1856             : static long
    1857           7 : key_search(GEN T, GEN x, GEN code)
    1858             : {
    1859           7 :   long u = lg(T)-1, i, l, s;
    1860             : 
    1861           7 :   if (!u) return 0;
    1862           7 :   l = 1; x = closure_callgen1(code, x);
    1863             :   do
    1864             :   {
    1865          14 :     i = (l+u)>>1; s = lexcmp(x, closure_callgen1(code, gel(T,i)));
    1866          14 :     if (!s) return i;
    1867           7 :     if (s<0) u=i-1; else l=i+1;
    1868           7 :   } while (u>=l);
    1869           0 :   return 0;
    1870             : }
    1871             : long
    1872      126105 : vecsearch(GEN v, GEN x, GEN k)
    1873             : {
    1874      126105 :   pari_sp av = avma;
    1875             :   void *E;
    1876      126105 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, v, k);
    1877      126098 :   long r, tv = typ(v);
    1878      126098 :   if (tv == t_VECSMALL)
    1879          21 :     x = (GEN)itos(x);
    1880      126077 :   else if (!is_matvec_t(tv)) pari_err_TYPE("vecsearch", v);
    1881      126098 :   r = CMP? gen_search(v, x, E, CMP): key_search(v, x, k);
    1882      126098 :   return gc_long(av, r < 0? 0: r);
    1883             : }
    1884             : 
    1885             : GEN
    1886        1661 : ZV_indexsort(GEN L) { return gen_indexsort(L, (void*)&cmpii, &cmp_nodata); }
    1887             : GEN
    1888          63 : ZV_sort(GEN L) { return gen_sort(L, (void*)&cmpii, &cmp_nodata); }
    1889             : GEN
    1890       40313 : ZV_sort_uniq(GEN L) { return gen_sort_uniq(L, (void*)&cmpii, &cmp_nodata); }
    1891             : void
    1892     1153359 : ZV_sort_inplace(GEN L) { gen_sort_inplace(L, (void*)&cmpii, &cmp_nodata,NULL); }
    1893             : GEN
    1894       21511 : ZV_sort_uniq_shallow(GEN L)
    1895             : {
    1896       21511 :   GEN v = gen_indexsort_uniq(L, (void*)&cmpii, &cmp_nodata);
    1897       21511 :   return vecpermute(L, v);
    1898             : }
    1899             : GEN
    1900        1372 : ZV_sort_shallow(GEN L)
    1901             : {
    1902        1372 :   GEN v = gen_indexsort(L, (void*)&cmpii, &cmp_nodata);
    1903        1372 :   return vecpermute(L, v);
    1904             : }
    1905             : 
    1906             : GEN
    1907        1085 : vec_equiv(GEN F)
    1908             : {
    1909        1085 :   pari_sp av = avma;
    1910        1085 :   long j, k, L = lg(F);
    1911        1085 :   GEN w = cgetg(L, t_VEC);
    1912        1085 :   GEN perm = gen_indexsort(F, (void*)&cmp_universal, cmp_nodata);
    1913        3094 :   for (j = k = 1; j < L;)
    1914             :   {
    1915        2009 :     GEN v = cgetg(L, t_VECSMALL);
    1916        2009 :     long l = 1, o = perm[j];
    1917        2009 :     v[l++] = o;
    1918        4767 :     for (j++; j < L; v[l++] = perm[j++])
    1919        3682 :       if (!gequal(gel(F,o), gel(F, perm[j]))) break;
    1920        2009 :     setlg(v, l); gel(w, k++) = v;
    1921             :   }
    1922        1085 :   setlg(w, k); return gerepilecopy(av,w);
    1923             : }
    1924             : 
    1925             : GEN
    1926       15246 : vec_reduce(GEN v, GEN *pE)
    1927             : {
    1928       15246 :   GEN E, F, P = gen_indexsort(v, (void*)cmp_universal, cmp_nodata);
    1929             :   long i, m, l;
    1930       15246 :   F = cgetg_copy(v, &l);
    1931       15246 :   *pE = E = cgetg(l, t_VECSMALL);
    1932       37821 :   for (i = m = 1; i < l;)
    1933             :   {
    1934       22575 :     GEN u = gel(v, P[i]);
    1935             :     long k;
    1936       27832 :     for(k = i + 1; k < l; k++)
    1937       12593 :       if (cmp_universal(gel(v, P[k]), u)) break;
    1938       22575 :     E[m] = k - i; gel(F, m) = u; i = k; m++;
    1939             :   }
    1940       15246 :   setlg(F, m);
    1941       15246 :   setlg(E, m); return F;
    1942             : }
    1943             : 
    1944             : /********************************************************************/
    1945             : /**                      SEARCH IN SORTED VECTOR                   **/
    1946             : /********************************************************************/
    1947             : /* index of x in table T, 0 otherwise */
    1948             : long
    1949     1124417 : tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
    1950             : {
    1951     1124417 :   long l = 1, u = lg(T)-1, i, s;
    1952             : 
    1953     8163421 :   while (u>=l)
    1954             :   {
    1955     8119617 :     i = (l+u)>>1; s = cmp(x, gel(T,i));
    1956     8119615 :     if (!s) return i;
    1957     7039004 :     if (s<0) u=i-1; else l=i+1;
    1958             :   }
    1959       43804 :   return 0;
    1960             : }
    1961             : 
    1962             : /* looks if x belongs to the set T and returns the index if yes, 0 if no */
    1963             : long
    1964    23574705 : gen_search(GEN T, GEN x, void *data, int (*cmp)(void*,GEN,GEN))
    1965             : {
    1966    23574705 :   long u = lg(T)-1, i, l, s;
    1967             : 
    1968    23574705 :   if (!u) return -1;
    1969    23574684 :   l = 1;
    1970             :   do
    1971             :   {
    1972   110212893 :     i = (l+u) >> 1; s = cmp(data, x, gel(T,i));
    1973   110212893 :     if (!s) return i;
    1974    86691717 :     if (s < 0) u = i-1; else l = i+1;
    1975    86691717 :   } while (u >= l);
    1976       53508 :   return -((s < 0)? i: i+1);
    1977             : }
    1978             : 
    1979             : long
    1980     1073674 : ZV_search(GEN x, GEN y) { return tablesearch(x, y, cmpii); }
    1981             : 
    1982             : long
    1983     9338308 : zv_search(GEN T, long x)
    1984             : {
    1985     9338308 :   long l = 1, u = lg(T)-1;
    1986    42643550 :   while (u>=l)
    1987             :   {
    1988    35642360 :     long i = (l+u)>>1;
    1989    35642360 :     if (x < T[i]) u = i-1;
    1990    19201470 :     else if (x > T[i]) l = i+1;
    1991     2337118 :     else return i;
    1992             :   }
    1993     7001190 :   return 0;
    1994             : }
    1995             : 
    1996             : /********************************************************************/
    1997             : /**                   COMPARISON FUNCTIONS                         **/
    1998             : /********************************************************************/
    1999             : int
    2000   644502689 : cmp_nodata(void *data, GEN x, GEN y)
    2001             : {
    2002   644502689 :   int (*cmp)(GEN,GEN)=(int (*)(GEN,GEN)) data;
    2003   644502689 :   return cmp(x,y);
    2004             : }
    2005             : 
    2006             : /* assume x and y come from the same idealprimedec call (uniformizer unique) */
    2007             : int
    2008     3100416 : cmp_prime_over_p(GEN x, GEN y)
    2009             : {
    2010     3100416 :   long k = pr_get_f(x) - pr_get_f(y); /* diff. between residue degree */
    2011      137143 :   return k? ((k > 0)? 1: -1)
    2012     3237539 :           : ZV_cmp(pr_get_gen(x), pr_get_gen(y));
    2013             : }
    2014             : 
    2015             : int
    2016      347312 : cmp_prime_ideal(GEN x, GEN y)
    2017             : {
    2018      347312 :   int k = cmpii(pr_get_p(x), pr_get_p(y));
    2019      347312 :   return k? k: cmp_prime_over_p(x,y);
    2020             : }
    2021             : 
    2022             : /* assume x and y are t_POL in the same variable whose coeffs can be
    2023             :  * compared (used to sort polynomial factorizations) */
    2024             : int
    2025     5358817 : gen_cmp_RgX(void *data, GEN x, GEN y)
    2026             : {
    2027     5358817 :   int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
    2028     5358817 :   long i, lx = lg(x), ly = lg(y);
    2029             :   int fl;
    2030     5358817 :   if (lx > ly) return  1;
    2031     5320743 :   if (lx < ly) return -1;
    2032    12086841 :   for (i=lx-1; i>1; i--)
    2033    11449044 :     if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
    2034      637797 :   return 0;
    2035             : }
    2036             : 
    2037             : static int
    2038        3613 : cmp_RgX_Rg(GEN x, GEN y)
    2039             : {
    2040        3613 :   long lx = lgpol(x), ly;
    2041        3613 :   if (lx > 1) return  1;
    2042           0 :   ly = gequal0(y) ? 0:1;
    2043           0 :   if (lx > ly) return  1;
    2044           0 :   if (lx < ly) return -1;
    2045           0 :   if (lx==0) return 0;
    2046           0 :   return gcmp(gel(x,2), y);
    2047             : }
    2048             : int
    2049      110606 : cmp_RgX(GEN x, GEN y)
    2050             : {
    2051      110606 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    2052      110606 :   if (typ(y) == t_POLMOD) y = gel(y,2);
    2053      110606 :   if (typ(x) == t_POL) {
    2054       55695 :     if (typ(y) != t_POL) return cmp_RgX_Rg(x, y);
    2055             :   } else {
    2056       54911 :     if (typ(y) != t_POL) return gcmp(x,y);
    2057        3361 :     return - cmp_RgX_Rg(y,x);
    2058             :   }
    2059       55443 :   return gen_cmp_RgX((void*)&gcmp,x,y);
    2060             : }
    2061             : 
    2062             : int
    2063      322012 : cmp_Flx(GEN x, GEN y)
    2064             : {
    2065      322012 :   long i, lx = lg(x), ly = lg(y);
    2066      322012 :   if (lx > ly) return  1;
    2067      306000 :   if (lx < ly) return -1;
    2068      547002 :   for (i=lx-1; i>1; i--)
    2069      462135 :     if (uel(x,i) != uel(y,i)) return uel(x,i)<uel(y,i)? -1: 1;
    2070       84867 :   return 0;
    2071             : }
    2072             : /********************************************************************/
    2073             : /**                   MERGE & SORT FACTORIZATIONS                  **/
    2074             : /********************************************************************/
    2075             : /* merge fx, fy two factorizations, whose 1st column is sorted in strictly
    2076             :  * increasing order wrt cmp */
    2077             : GEN
    2078      690170 : merge_factor(GEN fx, GEN fy, void *data, int (*cmp)(void *,GEN,GEN))
    2079             : {
    2080      690170 :   GEN x = gel(fx,1), e = gel(fx,2), M, E;
    2081      690170 :   GEN y = gel(fy,1), f = gel(fy,2);
    2082      690170 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    2083             : 
    2084      690170 :   M = cgetg(l, t_COL);
    2085      690170 :   E = cgetg(l, t_COL);
    2086             : 
    2087      690170 :   m = ix = iy = 1;
    2088    10039932 :   while (ix<lx && iy<ly)
    2089             :   {
    2090     9349762 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    2091     9349762 :     if (s < 0)
    2092     8715030 :     { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; }
    2093      634732 :     else if (s == 0)
    2094             :     {
    2095       94983 :       GEN z = gel(x,ix), g = addii(gel(e,ix), gel(f,iy));
    2096       94983 :       iy++; ix++; if (!signe(g)) continue;
    2097       11046 :       gel(M,m) = z; gel(E,m) = g;
    2098             :     }
    2099             :     else
    2100      539749 :     { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; }
    2101     9265825 :     m++;
    2102             :   }
    2103     4858237 :   while (ix<lx) { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; m++; }
    2104      932333 :   while (iy<ly) { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; m++; }
    2105      690170 :   setlg(M, m);
    2106      690170 :   setlg(E, m); return mkmat2(M, E);
    2107             : }
    2108             : /* merge two sorted vectors, removing duplicates. Shallow */
    2109             : GEN
    2110      454532 : merge_sort_uniq(GEN x, GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    2111             : {
    2112      454532 :   long i, j, k, lx = lg(x), ly = lg(y);
    2113      454532 :   GEN z = cgetg(lx + ly - 1, typ(x));
    2114      454531 :   i = j = k = 1;
    2115      594686 :   while (i<lx && j<ly)
    2116             :   {
    2117      140153 :     int s = cmp(data, gel(x,i), gel(y,j));
    2118      140155 :     if (s < 0)
    2119      118962 :       gel(z,k++) = gel(x,i++);
    2120       21193 :     else if (s > 0)
    2121       21172 :       gel(z,k++) = gel(y,j++);
    2122             :     else
    2123          21 :     { gel(z,k++) = gel(x,i++); j++; }
    2124             :   }
    2125      809215 :   while (i<lx) gel(z,k++) = gel(x,i++);
    2126      582977 :   while (j<ly) gel(z,k++) = gel(y,j++);
    2127      454533 :   setlg(z, k); return z;
    2128             : }
    2129             : /* in case of equal keys in x,y, take the key from x */
    2130             : static GEN
    2131       34615 : ZV_union_shallow_t(GEN x, GEN y, long t)
    2132             : {
    2133       34615 :   long i, j, k, lx = lg(x), ly = lg(y);
    2134       34615 :   GEN z = cgetg(lx + ly - 1, t);
    2135       34615 :   i = j = k = 1;
    2136       78001 :   while (i<lx && j<ly)
    2137             :   {
    2138       43386 :     int s = cmpii(gel(x,i), gel(y,j));
    2139       43386 :     if (s < 0)
    2140       23576 :       gel(z,k++) = gel(x,i++);
    2141       19810 :     else if (s > 0)
    2142       10332 :       gel(z,k++) = gel(y,j++);
    2143             :     else
    2144        9478 :     { gel(z,k++) = gel(x,i++); j++; }
    2145             :   }
    2146       41734 :   while (i < lx) gel(z,k++) = gel(x,i++);
    2147       70098 :   while (j < ly) gel(z,k++) = gel(y,j++);
    2148       34615 :   setlg(z, k); return z;
    2149             : }
    2150             : GEN
    2151       34433 : ZV_union_shallow(GEN x, GEN y)
    2152       34433 : { return ZV_union_shallow_t(x, y, t_VEC); }
    2153             : GEN
    2154         182 : ZC_union_shallow(GEN x, GEN y)
    2155         182 : { return ZV_union_shallow_t(x, y, t_COL); }
    2156             : 
    2157             : /* sort generic factorization, in place */
    2158             : GEN
    2159     9769377 : sort_factor(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    2160             : {
    2161             :   GEN a, b, A, B, w;
    2162             :   pari_sp av;
    2163             :   long n, i;
    2164             : 
    2165     9769377 :   a = gel(y,1); n = lg(a); if (n == 1) return y;
    2166     9749876 :   b = gel(y,2); av = avma;
    2167     9749876 :   A = new_chunk(n);
    2168     9750068 :   B = new_chunk(n);
    2169     9750696 :   w = gen_sortspec(a, n-1, data, cmp);
    2170    29780669 :   for (i=1; i<n; i++) { long k=w[i]; gel(A,i) = gel(a,k); gel(B,i) = gel(b,k); }
    2171    29781781 :   for (i=1; i<n; i++) { gel(a,i) = gel(A,i); gel(b,i) = gel(B,i); }
    2172     9751004 :   set_avma(av); return y;
    2173             : }
    2174             : /* sort polynomial factorization, in place */
    2175             : GEN
    2176     1791951 : sort_factor_pol(GEN y,int (*cmp)(GEN,GEN))
    2177             : {
    2178     1791951 :   (void)sort_factor(y,(void*)cmp, &gen_cmp_RgX);
    2179     1791960 :   return y;
    2180             : }
    2181             : 
    2182             : /***********************************************************************/
    2183             : /*                                                                     */
    2184             : /*                          SET OPERATIONS                             */
    2185             : /*                                                                     */
    2186             : /***********************************************************************/
    2187             : GEN
    2188      227045 : gtoset(GEN x)
    2189             : {
    2190             :   long lx;
    2191      227045 :   if (!x) return cgetg(1, t_VEC);
    2192      227045 :   switch(typ(x))
    2193             :   {
    2194      227017 :     case t_VEC:
    2195      227017 :     case t_COL: lx = lg(x); break;
    2196          14 :     case t_LIST:
    2197          14 :       if (list_typ(x)==t_LIST_MAP) return mapdomain(x);
    2198          14 :       x = list_data(x); lx = x? lg(x): 1; break;
    2199           7 :     case t_VECSMALL: lx = lg(x); x = zv_to_ZV(x); break;
    2200           7 :     default: return mkveccopy(x);
    2201             :   }
    2202      227038 :   if (lx==1) return cgetg(1,t_VEC);
    2203      226863 :   x = gen_sort_uniq(x, (void*)&cmp_universal, cmp_nodata);
    2204      226863 :   settyp(x, t_VEC); /* it may be t_COL */
    2205      226863 :   return x;
    2206             : }
    2207             : 
    2208             : long
    2209          14 : setisset(GEN x)
    2210             : {
    2211          14 :   long i, lx = lg(x);
    2212             : 
    2213          14 :   if (typ(x) != t_VEC) return 0;
    2214          14 :   if (lx == 1) return 1;
    2215          70 :   for (i=1; i<lx-1; i++)
    2216          63 :     if (cmp_universal(gel(x,i+1), gel(x,i)) <= 0) return 0;
    2217           7 :   return 1;
    2218             : }
    2219             : 
    2220             : long
    2221       83622 : setsearch(GEN T, GEN y, long flag)
    2222             : {
    2223             :   long i, lx;
    2224       83622 :   switch(typ(T))
    2225             :   {
    2226       83608 :     case t_VEC: lx = lg(T); break;
    2227           7 :     case t_LIST:
    2228           7 :     if (list_typ(T) != t_LIST_RAW) pari_err_TYPE("setsearch",T);
    2229           7 :     T = list_data(T); lx = T? lg(T): 1; break;
    2230           7 :     default: pari_err_TYPE("setsearch",T);
    2231             :       return 0; /*LCOV_EXCL_LINE*/
    2232             :   }
    2233       83615 :   if (lx==1) return flag? 1: 0;
    2234       83615 :   i = gen_search(T,y,(void*)cmp_universal,cmp_nodata);
    2235       83615 :   if (i > 0) return flag? 0: i;
    2236          56 :   return flag ? -i: 0;
    2237             : }
    2238             : 
    2239             : GEN
    2240           7 : setunion_i(GEN x, GEN y)
    2241           7 : { return merge_sort_uniq(x,y, (void*)cmp_universal, cmp_nodata); }
    2242             : 
    2243             : GEN
    2244           7 : setunion(GEN x, GEN y)
    2245             : {
    2246           7 :   pari_sp av = avma;
    2247           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setunion",x);
    2248           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setunion",y);
    2249           7 :   return gerepilecopy(av, setunion_i(x, y));
    2250             : }
    2251             : 
    2252             : GEN
    2253          14 : setdelta(GEN x, GEN y)
    2254             : {
    2255          14 :   long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
    2256          14 :   pari_sp av = avma;
    2257          14 :   GEN z = cgetg(lx + ly - 1,t_VEC);
    2258          14 :   if (typ(x) != t_VEC) pari_err_TYPE("setdelta",x);
    2259          14 :   if (typ(y) != t_VEC) pari_err_TYPE("setdelta",y);
    2260          84 :   while (ix < lx && iy < ly)
    2261             :   {
    2262          70 :     int c = cmp_universal(gel(x,ix), gel(y,iy));
    2263          70 :     if      (c < 0) gel(z, iz++) = gel(x,ix++);
    2264          42 :     else if (c > 0) gel(z, iz++) = gel(y,iy++);
    2265          28 :     else { ix++; iy++; }
    2266             :   }
    2267          21 :   while (ix<lx) gel(z,iz++) = gel(x,ix++);
    2268          14 :   while (iy<ly) gel(z,iz++) = gel(y,iy++);
    2269          14 :   setlg(z,iz); return gerepilecopy(av,z);
    2270             : }
    2271             : 
    2272             : GEN
    2273           7 : setintersect(GEN x, GEN y)
    2274             : {
    2275           7 :   long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
    2276           7 :   pari_sp av = avma;
    2277           7 :   GEN z = cgetg(lx,t_VEC);
    2278           7 :   if (typ(x) != t_VEC) pari_err_TYPE("setintersect",x);
    2279           7 :   if (typ(y) != t_VEC) pari_err_TYPE("setintersect",y);
    2280          70 :   while (ix < lx && iy < ly)
    2281             :   {
    2282          63 :     int c = cmp_universal(gel(x,ix), gel(y,iy));
    2283          63 :     if      (c < 0) ix++;
    2284          35 :     else if (c > 0) iy++;
    2285          21 :     else { gel(z, iz++) = gel(x,ix); ix++; iy++; }
    2286             :   }
    2287           7 :   setlg(z,iz); return gerepilecopy(av,z);
    2288             : }
    2289             : 
    2290             : GEN
    2291        1001 : gen_setminus(GEN A, GEN B, int (*cmp)(GEN,GEN))
    2292             : {
    2293        1001 :   pari_sp ltop = avma;
    2294        1001 :   long i = 1, j = 1, k = 1, lx = lg(A), ly = lg(B);
    2295        1001 :   GEN  diff = cgetg(lx,t_VEC);
    2296        7460 :   while (i < lx && j < ly)
    2297        5844 :     switch ( cmp(gel(A,i),gel(B,j)) )
    2298             :     {
    2299        1050 :       case -1: gel(diff,k++) = gel(A,i++); break;
    2300        2246 :       case 1: j++; break;
    2301        2548 :       case 0: i++; break;
    2302             :     }
    2303       11281 :   while (i < lx) gel(diff,k++) = gel(A,i++);
    2304        1001 :   setlg(diff,k);
    2305        1001 :   return gerepilecopy(ltop,diff);
    2306             : }
    2307             : 
    2308             : GEN
    2309        1001 : setminus(GEN x, GEN y)
    2310             : {
    2311        1001 :   if (typ(x) != t_VEC) pari_err_TYPE("setminus",x);
    2312        1001 :   if (typ(y) != t_VEC) pari_err_TYPE("setminus",y);
    2313        1001 :   return gen_setminus(x,y,cmp_universal);
    2314             : }
    2315             : 
    2316             : GEN
    2317          21 : setbinop(GEN f, GEN x, GEN y)
    2318             : {
    2319          21 :   pari_sp av = avma;
    2320          21 :   long i, j, lx, ly, k = 1;
    2321             :   GEN z;
    2322          21 :   if (typ(f) != t_CLOSURE || closure_arity(f) != 2 || closure_is_variadic(f))
    2323           7 :     pari_err_TYPE("setbinop [function needs exactly 2 arguments]",f);
    2324          14 :   lx = lg(x);
    2325          14 :   if (typ(x) != t_VEC) pari_err_TYPE("setbinop", x);
    2326          14 :   if (y == NULL) { /* assume x = y and f symmetric */
    2327           7 :     z = cgetg((((lx-1)*lx) >> 1) + 1, t_VEC);
    2328          28 :     for (i = 1; i < lx; i++)
    2329          63 :       for (j = i; j < lx; j++)
    2330          42 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(x,j));
    2331             :   } else {
    2332           7 :     ly = lg(y);
    2333           7 :     if (typ(y) != t_VEC) pari_err_TYPE("setbinop", y);
    2334           7 :     z = cgetg((lx-1)*(ly-1) + 1, t_VEC);
    2335          28 :     for (i = 1; i < lx; i++)
    2336          84 :       for (j = 1; j < ly; j++)
    2337          63 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(y,j));
    2338             :   }
    2339          14 :   return gerepileupto(av, gtoset(z));
    2340             : }

Generated by: LCOV version 1.14