Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - bibli2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19044-129ab8a) Lines: 1042 1076 96.8 %
Date: 2016-06-27 Functions: 90 93 96.8 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 794 957 83.0 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : #include "pari.h"
      15                 :            : #include "paripriv.h"
      16                 :            : 
      17                 :            : /*******************************************************************/
      18                 :            : /**                                                               **/
      19                 :            : /**                      SPECIAL POLYNOMIALS                      **/
      20                 :            : /**                                                               **/
      21                 :            : /*******************************************************************/
      22                 :            : /* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
      23                 :            :  * T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
      24                 :            :  *   where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
      25                 :            :  *   and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
      26                 :            : GEN
      27                 :       2156 : polchebyshev1(long n, long v) /* Assume 4*n < LONG_MAX */
      28                 :            : {
      29                 :            :   long k, l;
      30                 :            :   pari_sp av;
      31                 :            :   GEN q,a,r;
      32                 :            : 
      33         [ +  + ]:       2156 :   if (v<0) v = 0;
      34                 :            :   /* polchebyshev(-n,1) = polchebyshev(n,1) */
      35         [ +  + ]:       2156 :   if (n < 0) n = -n;
      36         [ +  + ]:       2156 :   if (n==0) return pol_1(v);
      37         [ +  + ]:       2135 :   if (n==1) return pol_x(v);
      38                 :            : 
      39                 :       2093 :   q = cgetg(n+3, t_POL); r = q + n+2;
      40                 :       2093 :   a = int2n(n-1);
      41                 :       2093 :   gel(r--,0) = a;
      42                 :       2093 :   gel(r--,0) = gen_0;
      43         [ +  + ]:      31955 :   for (k=1,l=n; l>1; k++,l-=2)
      44                 :            :   {
      45                 :      29862 :     av = avma;
      46                 :      29862 :     a = diviuuexact(muluui(l, l-1, a), 4*k, n-k);
      47                 :      29862 :     togglesign(a); a = gerepileuptoint(av, a);
      48                 :      29862 :     gel(r--,0) = a;
      49                 :      29862 :     gel(r--,0) = gen_0;
      50                 :            :   }
      51                 :       2093 :   q[1] = evalsigne(1) | evalvarn(v);
      52                 :       2156 :   return q;
      53                 :            : }
      54                 :            : static void
      55                 :         70 : polchebyshev1_eval_aux(long n, GEN x, GEN *pt1, GEN *pt2)
      56                 :            : {
      57                 :            :   GEN t1, t2, b;
      58         [ +  + ]:         70 :   if (n == 1) { *pt1 = gen_1; *pt2 = x; return; }
      59         [ -  + ]:         56 :   if (n == 0) { *pt1 = x; *pt2 = gen_1; return; }
      60                 :         56 :   polchebyshev1_eval_aux((n+1) >> 1, x, &t1, &t2);
      61                 :         56 :   b = gsub(gmul(gmul2n(t1,1), t2), x);
      62         [ +  + ]:         56 :   if (odd(n)) { *pt1 = gadd(gmul2n(gsqr(t1), 1), gen_m1); *pt2 = b; }
      63                 :         70 :   else        { *pt1 = b; *pt2 = gadd(gmul2n(gsqr(t2), 1), gen_m1); }
      64                 :            : }
      65                 :            : static GEN
      66                 :         14 : polchebyshev1_eval(long n, GEN x)
      67                 :            : {
      68                 :            :   GEN t1, t2;
      69                 :            :   long i, v;
      70                 :            :   pari_sp av;
      71                 :            : 
      72         [ -  + ]:         14 :   if (n < 0) n = -n;
      73         [ -  + ]:         14 :   if (n==0) return gen_1;
      74         [ -  + ]:         14 :   if (n==1) return gcopy(x);
      75                 :         14 :   av = avma;
      76                 :         14 :   v = u_lvalrem(n, 2, (ulong*)&n);
      77                 :         14 :   polchebyshev1_eval_aux((n+1)>>1, x, &t1, &t2);
      78         [ +  - ]:         14 :   if (n != 1) t2 = gsub(gmul(gmul2n(t1,1), t2), x);
      79         [ +  + ]:         35 :   for (i = 1; i <= v; i++) t2 = gadd(gmul2n(gsqr(t2), 1), gen_m1);
      80                 :         14 :   return gerepileupto(av, t2);
      81                 :            : }
      82                 :            : 
      83                 :            : /* Chebychev  polynomial of the second kind U(n,x): the coefficient in front of
      84                 :            :  * x^(n-2*m) is (-1)^m * 2^(n-2m)*(n-m)!/m!/(n-2m)!  for m=0,1,...,n/2 */
      85                 :            : GEN
      86                 :       2135 : polchebyshev2(long n, long v)
      87                 :            : {
      88                 :            :   pari_sp av;
      89                 :            :   GEN q, a, r;
      90                 :            :   long m;
      91                 :       2135 :   int neg = 0;
      92                 :            : 
      93         [ -  + ]:       2135 :   if (v<0) v = 0;
      94                 :            :   /* polchebyshev(-n,2) = -polchebyshev(n-2,2) */
      95         [ +  + ]:       2135 :   if (n < 0) {
      96         [ +  + ]:       1050 :     if (n == -1) return zeropol(v);
      97                 :       1029 :     neg = 1; n = -n-2;
      98                 :            :   }
      99 [ +  + ][ +  + ]:       2114 :   if (n==0) return neg ? scalar_ZX_shallow(gen_m1, v): pol_1(v);
     100                 :            : 
     101                 :       2072 :   q = cgetg(n+3, t_POL); r = q + n+2;
     102                 :       2072 :   a = int2n(n);
     103         [ +  + ]:       2072 :   if (neg) togglesign(a);
     104                 :       2072 :   gel(r--,0) = a;
     105                 :       2072 :   gel(r--,0) = gen_0;
     106         [ +  + ]:      30807 :   for (m=1; 2*m<= n; m++)
     107                 :            :   {
     108                 :      28735 :     av = avma;
     109                 :      28735 :     a = diviuuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m, n-m+1);
     110                 :      28735 :     togglesign(a); a = gerepileuptoint(av, a);
     111                 :      28735 :     gel(r--,0) = a;
     112                 :      28735 :     gel(r--,0) = gen_0;
     113                 :            :   }
     114                 :       2072 :   q[1] = evalsigne(1) | evalvarn(v);
     115                 :       2135 :   return q;
     116                 :            : }
     117                 :            : static void
     118                 :         91 : polchebyshev2_eval_aux(long n, GEN x, GEN *pu1, GEN *pu2)
     119                 :            : {
     120                 :            :   GEN u1, u2, u, mu1;
     121         [ +  + ]:         91 :   if (n == 1) { *pu1 = gen_1; *pu2 = gmul2n(x,1); return; }
     122         [ -  + ]:         70 :   if (n == 0) { *pu1 = gen_0; *pu2 = gen_1; return; }
     123                 :         70 :   polchebyshev2_eval_aux(n >> 1, x, &u1, &u2);
     124                 :         70 :   mu1 = gneg(u1);
     125                 :         70 :   u = gmul(gadd(u2,u1), gadd(u2,mu1));
     126         [ +  + ]:         70 :   if (odd(n)) { *pu1 = u; *pu2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1)); }
     127                 :         91 :   else        { *pu2 = u; *pu1 = gmul(gmul2n(u1,1), gadd(u2, gmul(x,mu1))); }
     128                 :            : }
     129                 :            : static GEN
     130                 :         35 : polchebyshev2_eval(long n, GEN x)
     131                 :            : {
     132                 :            :   GEN u1, u2, mu1;
     133                 :         35 :   long neg = 0;
     134                 :            :   pari_sp av;
     135                 :            : 
     136         [ +  + ]:         35 :   if (n < 0) {
     137         [ +  + ]:         14 :     if (n == -1) return gen_0;
     138                 :          7 :     neg = 1; n = -n-2;
     139                 :            :   }
     140 [ +  + ][ +  - ]:         28 :   if (n==0) return neg ? gen_m1: gen_1;
     141                 :         21 :   av = avma;
     142                 :         21 :   polchebyshev2_eval_aux(n>>1, x, &u1, &u2);
     143                 :         21 :   mu1 = gneg(u1);
     144         [ +  + ]:         21 :   if (odd(n)) u2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1));
     145                 :         14 :   else        u2 = gmul(gadd(u2,u1), gadd(u2,mu1));
     146         [ -  + ]:         21 :   if (neg) u2 = gneg(u2);
     147                 :         35 :   return gerepileupto(av, u2);
     148                 :            : }
     149                 :            : 
     150                 :            : GEN
     151                 :       4284 : polchebyshev(long n, long kind, long v)
     152                 :            : {
     153      [ +  +  - ]:       4284 :   switch (kind)
     154                 :            :   {
     155                 :       2149 :     case 1: return polchebyshev1(n, v);
     156                 :       2135 :     case 2: return polchebyshev2(n, v);
     157                 :          0 :     default: pari_err_FLAG("polchebyshev");
     158                 :            :   }
     159                 :       4284 :   return NULL; /* not reached */
     160                 :            : }
     161                 :            : GEN
     162                 :       4333 : polchebyshev_eval(long n, long kind, GEN x)
     163                 :            : {
     164         [ +  + ]:       4333 :   if (!x) return polchebyshev(n, kind, 0);
     165         [ +  + ]:         63 :   if (gequalX(x)) return polchebyshev(n, kind, varn(x));
     166      [ +  +  - ]:         49 :   switch (kind)
     167                 :            :   {
     168                 :         14 :     case 1: return polchebyshev1_eval(n, x);
     169                 :         35 :     case 2: return polchebyshev2_eval(n, x);
     170                 :          0 :     default: pari_err_FLAG("polchebyshev");
     171                 :            :   }
     172                 :       4333 :   return NULL; /* not reached */
     173                 :            : }
     174                 :            : 
     175                 :            : /* Hermite polynomial H(n,x):  H(n+1) = 2x H(n) - 2n H(n-1)
     176                 :            :  * The coefficient in front of x^(n-2*m) is
     177                 :            :  * (-1)^m * n! * 2^(n-2m)/m!/(n-2m)!  for m=0,1,...,n/2.. */
     178                 :            : GEN
     179                 :       1428 : polhermite(long n, long v)
     180                 :            : {
     181                 :            :   long m;
     182                 :            :   pari_sp av;
     183                 :            :   GEN q,a,r;
     184                 :            : 
     185         [ -  + ]:       1428 :   if (v<0) v = 0;
     186         [ -  + ]:       1428 :   if (n < 0) pari_err_DOMAIN("polhermite", "degree", "<", gen_0, stoi(n));
     187         [ +  + ]:       1428 :   if (n==0) return pol_1(v);
     188                 :            : 
     189                 :       1421 :   q = cgetg(n+3, t_POL); r = q + n+2;
     190                 :       1421 :   a = int2n(n);
     191                 :       1421 :   gel(r--,0) = a;
     192                 :       1421 :   gel(r--,0) = gen_0;
     193         [ +  + ]:      40285 :   for (m=1; 2*m<= n; m++)
     194                 :            :   {
     195                 :      38864 :     av = avma;
     196                 :      38864 :     a = diviuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m);
     197                 :      38864 :     togglesign(a);
     198                 :      38864 :     gel(r--,0) = a = gerepileuptoint(av, a);
     199                 :      38864 :     gel(r--,0) = gen_0;
     200                 :            :   }
     201                 :       1421 :   q[1] = evalsigne(1) | evalvarn(v);
     202                 :       1428 :   return q;
     203                 :            : }
     204                 :            : GEN
     205                 :       1442 : polhermite_eval(long n, GEN x)
     206                 :            : {
     207                 :            :   long i;
     208                 :            :   pari_sp av, av2;
     209                 :            :   GEN x2, u, v;
     210                 :            : 
     211         [ +  + ]:       1442 :   if (!x) return polhermite(n, 0);
     212         [ -  + ]:         14 :   if (gequalX(x)) return polhermite(n, varn(x));
     213         [ -  + ]:         14 :   if (n==0) return gen_1;
     214         [ -  + ]:         14 :   if (n==1) return gmul2n(x,1);
     215                 :         14 :   av = avma; x2 = gmul2n(x,1); v = gen_1; u = x2;
     216                 :         14 :   av2= avma;
     217         [ +  + ]:       7035 :   for (i=1; i<n; i++)
     218                 :            :   { /* u = H_i(x), v = H_{i-1}(x), compute t = H_{i+1}(x) */
     219                 :            :     GEN t;
     220         [ +  + ]:       7021 :     if ((i & 0xff) == 0) gerepileall(av2,2,&u, &v);
     221                 :       7021 :     t = gsub(gmul(x2, u), gmulsg(2*i,v));
     222                 :       7021 :     v = u; u = t;
     223                 :            :   }
     224                 :       1442 :   return gerepileupto(av, u);
     225                 :            : }
     226                 :            : 
     227                 :            : /* Legendre polynomial
     228                 :            :  * L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1)
     229                 :            :  * L(n) = 2^-n sum_{k=0}^{n/2} a_k x^(n-2k)
     230                 :            :  *   where a_k = (-1)^k (2n-2k)! / k! (n-k)! (n-2k)! is an integer
     231                 :            :  *   and a_0 = binom(2n,n), a_k / a_{k-1} = - (n-2k+1)(n-2k+2) / 2k (2n-2k+1) */
     232                 :            : GEN
     233                 :       2184 : pollegendre(long n, long v)
     234                 :            : {
     235                 :            :   long k, l;
     236                 :            :   pari_sp av;
     237                 :            :   GEN a, r, q;
     238                 :            : 
     239         [ -  + ]:       2184 :   if (v<0) v = 0;
     240                 :            :   /* pollegendre(-n) = pollegendre(n-1) */
     241         [ +  + ]:       2184 :   if (n < 0) n = -n-1;
     242         [ +  + ]:       2184 :   if (n==0) return pol_1(v);
     243         [ +  + ]:       2142 :   if (n==1) return pol_x(v);
     244                 :            : 
     245                 :       2100 :   av = avma;
     246                 :       2100 :   q = cgetg(n+3, t_POL); r = q + n+2;
     247                 :       2100 :   gel(r--,0) = a = binomialuu(n<<1,n);
     248                 :       2100 :   gel(r--,0) = gen_0;
     249         [ +  + ]:      32354 :   for (k=1,l=n; l>1; k++,l-=2)
     250                 :            :   { /* l = n-2*k+2 */
     251                 :      30254 :     av = avma;
     252                 :      30254 :     a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
     253                 :      30254 :     togglesign(a); a = gerepileuptoint(av, a);
     254                 :      30254 :     gel(r--,0) = a;
     255                 :      30254 :     gel(r--,0) = gen_0;
     256                 :            :   }
     257                 :       2100 :   q[1] = evalsigne(1) | evalvarn(v);
     258                 :       2184 :   return gerepileupto(av, gmul2n(q,-n));
     259                 :            : }
     260                 :            : 
     261                 :            : GEN
     262                 :       2163 : pollegendre_eval(long n, GEN x)
     263                 :            : {
     264                 :            :   long i;
     265                 :            :   pari_sp av;
     266                 :            :   GEN u, v;
     267                 :            : 
     268         [ +  + ]:       2163 :   if (!x) return pollegendre(n, 0);
     269         [ -  + ]:         14 :   if (gequalX(x)) return pollegendre(n, varn(x));
     270                 :            :   /* pollegendre(-n) = pollegendre(n-1) */
     271         [ -  + ]:         14 :   if (n < 0) n = -n-1;
     272         [ -  + ]:         14 :   if (n==0) return gen_1;
     273         [ -  + ]:         14 :   if (n==1) return gcopy(x);
     274                 :         14 :   av = avma; v = gen_1; u = x;
     275         [ +  + ]:       7035 :   for (i=1; i<n; i++)
     276                 :            :   { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
     277                 :            :     GEN t;
     278         [ +  + ]:       7021 :     if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
     279                 :       7021 :     t = gdivgs(gsub(gmul(gmulsg(2*i+1,x), u), gmulsg(i,v)), i+1);
     280                 :       7021 :     v = u; u = t;
     281                 :            :   }
     282                 :       2163 :   return gerepileupto(av, u);
     283                 :            : }
     284                 :            : 
     285                 :            : /* polcyclo(p) = X^(p-1) + ... + 1 */
     286                 :            : static GEN
     287                 :      35170 : polcyclo_prime(long p, long v)
     288                 :            : {
     289                 :      35170 :   GEN T = cgetg(p+2, t_POL);
     290                 :            :   long i;
     291                 :      35170 :   T[1] = evalsigne(1) | evalvarn(v);
     292         [ +  + ]:     217248 :   for (i = 2; i < p+2; i++) gel(T,i) = gen_1;
     293                 :      35170 :   return T;
     294                 :            : }
     295                 :            : 
     296                 :            : /* cyclotomic polynomial */
     297                 :            : GEN
     298                 :      39918 : polcyclo(long n, long v)
     299                 :            : {
     300                 :            :   long s, q, i, l;
     301                 :      39918 :   pari_sp av=avma;
     302                 :            :   GEN T, P;
     303                 :            : 
     304         [ +  + ]:      39918 :   if (v<0) v = 0;
     305         [ +  + ]:      39918 :   if (n < 3)
     306      [ +  +  - ]:       4747 :     switch(n)
     307                 :            :     {
     308                 :        868 :       case 1: return deg1pol_shallow(gen_1, gen_m1, v);
     309                 :       3879 :       case 2: return deg1pol_shallow(gen_1, gen_1, v);
     310                 :          0 :       default: pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     311                 :            :     }
     312                 :      35171 :   P = gel(factoru(n), 1); l = lg(P);
     313                 :      35170 :   s = P[1]; T = polcyclo_prime(s, v);
     314         [ +  + ]:      38033 :   for (i = 2; i < l; i++)
     315                 :            :   { /* Phi_{np}(X) = Phi_n(X^p) / Phi_n(X) */
     316                 :       2863 :     s *= P[i];
     317                 :       2863 :     T = RgX_div(RgX_inflate(T, P[i]), T);
     318                 :            :   }
     319                 :            :   /* s = squarefree part of n */
     320                 :      35170 :   q = n / s;
     321         [ +  + ]:      35170 :   if (q == 1) return gerepileupto(av, T);
     322                 :      39917 :   return gerepilecopy(av, RgX_inflate(T,q));
     323                 :            : }
     324                 :            : 
     325                 :            : /* cyclotomic polynomial */
     326                 :            : GEN
     327                 :      14385 : polcyclo_eval(long n, GEN x)
     328                 :            : {
     329                 :      14385 :   pari_sp av= avma;
     330                 :            :   GEN P, md, xd, yneg, ypos;
     331                 :            :   long l, s, i, j, q, tx;
     332                 :      14385 :   long root_of_1 = 0;
     333                 :            : 
     334         [ +  + ]:      14385 :   if (!x) return polcyclo(n, 0);
     335                 :      12943 :   tx = typ(x);
     336         [ +  + ]:      12943 :   if (gequalX(x)) return polcyclo(n, varn(x));
     337         [ -  + ]:      12719 :   if (n <= 0) pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
     338         [ -  + ]:      12719 :   if (n == 1) return gsubgs(x, 1);
     339 [ +  + ][ -  + ]:      12719 :   if (tx == t_INT && !signe(x)) return gen_1;
     340         [ +  + ]:      18109 :   while ((n & 3) == 0) { n >>= 1; x = gsqr(x); } /* Phi_4n(x) = Phi_2n(x^2) */
     341                 :            :   /* n not divisible by 4 */
     342         [ +  + ]:      12719 :   if (n == 2) return gerepileupto(av, gaddgs(x,1));
     343         [ +  + ]:        819 :   if (!odd(n)) { n >>= 1; x = gneg(x); } /* Phi_2n(x) = Phi_n(-x) for n>1 odd */
     344                 :            :   /* n odd > 2.  s largest squarefree divisor of n */
     345                 :        819 :   P = gel(factoru(n), 1); s = zv_prod(P);
     346                 :            :   /* replace n by largest squarefree divisor */
     347         [ +  + ]:        819 :   q = n/s; if (q != 1) { x = gpowgs(x, q); n = s; }
     348                 :        819 :   l = lg(P)-1;
     349                 :            :   /* n squarefree odd > 2, l distinct prime divisors. Now handle x = 1 or -1 */
     350         [ +  + ]:        819 :   if (tx == t_INT) { /* shortcut */
     351         [ +  + ]:        749 :     if (is_pm1(x))
     352                 :            :     {
     353                 :         56 :       avma = av;
     354 [ +  + ][ +  - ]:         56 :       if (signe(x) > 0 && l == 1) return utoipos(P[1]);
     355                 :         35 :       return gen_1;
     356                 :            :     }
     357                 :            :   } else {
     358         [ +  + ]:         70 :     if (gequal1(x))
     359                 :            :     { /* n is prime, return n; multiply by x to keep the type */
     360         [ +  + ]:         14 :       if (l == 1) return gerepileupto(av, gmulgs(x,n));
     361                 :          7 :       return gerepilecopy(av, x); /* else 1 */
     362                 :            :     }
     363         [ +  + ]:         56 :     if (gequalm1(x)) return gerepileupto(av, gneg(x)); /* -1 */
     364                 :            :   }
     365                 :            :   /* Heuristic: evaluation will probably not improve things */
     366 [ +  + ][ +  - ]:        742 :   if (tx == t_POL || tx == t_MAT || lg(x) > n)
                 [ +  + ]
     367                 :         17 :     return gerepileupto(av, poleval(polcyclo(n,0), x));
     368                 :            : 
     369                 :        725 :   xd = cgetg((1L<<l) + 1, t_VEC); /* the x^d, where d | n */
     370                 :        725 :   md = cgetg((1L<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
     371                 :        725 :   gel(xd, 1) = x;
     372                 :        725 :   md[1] = 1;
     373                 :            :   /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d).
     374                 :            :    * If x has exact order D, n = Dq, then the result is 0 if q = 1. Otherwise
     375                 :            :    * the factors with x^d-1, D|d are omitted and we multiply at the end by
     376                 :            :    *   prod_{d | q} d^mu(q/d) = q if prime, 1 otherwise */
     377                 :            :   /* We store the factors with mu(d)= 1 (resp.-1) in ypos (resp yneg).
     378                 :            :    * At the end we return ypos/yneg if mu(n)=1 and yneg/ypos if mu(n)=-1 */
     379                 :        725 :   ypos = gsubgs(x,1);
     380                 :        725 :   yneg = gen_1;
     381         [ +  + ]:       1520 :   for (i = 1; i <= l; i++)
     382                 :            :   {
     383                 :        795 :     long ti = 1L<<(i-1), p = P[i];
     384         [ +  + ]:       1674 :     for (j = 1; j <= ti; j++) {
     385                 :        879 :       GEN X = gpowgs(gel(xd,j), p), t = gsubgs(X,1);
     386                 :        879 :       gel(xd,ti+j) = X;
     387                 :        879 :       md[ti+j] = -md[j];
     388         [ +  + ]:        879 :       if (gequal0(t))
     389                 :            :       { /* x^d = 1; root_of_1 := the smallest index ti+j such that X == 1
     390                 :            :         * (whose bits code d: bit i-1 is set iff P[i] | d). If no such index
     391                 :            :         * exists, then root_of_1 remains 0. Do not multiply with X-1 if X = 1,
     392                 :            :         * we handle these factors at the end */
     393         [ +  + ]:         28 :         if (!root_of_1) root_of_1 = ti+j;
     394                 :            :       }
     395                 :            :       else
     396                 :            :       {
     397         [ +  + ]:        851 :         if (md[ti+j] == 1) ypos = gmul(ypos, t);
     398                 :        788 :         else               yneg = gmul(yneg, t);
     399                 :            :       }
     400                 :            :     }
     401                 :            :   }
     402         [ +  + ]:        725 :   ypos = odd(l)? gdiv(yneg,ypos): gdiv(ypos,yneg);
     403         [ +  + ]:        725 :   if (root_of_1)
     404                 :            :   {
     405                 :         21 :     GEN X = gel(xd,(1<<l)); /* = x^n = 1 */
     406                 :         21 :     long bitmask_q = (1<<l) - root_of_1;
     407                 :            :     /* bitmask_q encodes q = n/d: bit (i-1) is 1 iff P[i] | q */
     408                 :            : 
     409                 :            :     /* x is a root of unity.  If bitmask_q = 0, then x was a primitive n-th
     410                 :            :      * root of 1 and the result is zero. Return X - 1 to preserve type. */
     411         [ +  + ]:         21 :     if (!bitmask_q) return gerepileupto(av, gsubgs(X, 1));
     412                 :            :     /* x is a primitive d-th root of unity, where d|n and d<n: we
     413                 :            :      * must multiply ypos by if(isprime(n/d), n/d, 1) */
     414                 :          7 :     ypos = gmul(ypos, X); /* multiply by X = 1 to preserve type */
     415                 :            :     /* If bitmask_q = 1<<(i-1) for some i <= l, then q == P[i] and we multiply
     416                 :            :      * by P[i]; otherwise q is composite and nothing more needs to be done */
     417         [ +  - ]:          7 :     if (!(bitmask_q & (bitmask_q-1))) /* detects power of 2, since bitmask!=0 */
     418                 :            :     {
     419                 :          7 :       i = vals(bitmask_q)+1; /* q = P[i] */
     420                 :          7 :       ypos = gmulgs(ypos, P[i]);
     421                 :            :     }
     422                 :            :   }
     423                 :      14385 :   return gerepileupto(av, ypos);
     424                 :            : }
     425                 :            : /********************************************************************/
     426                 :            : /**                                                                **/
     427                 :            : /**                  HILBERT & PASCAL MATRICES                     **/
     428                 :            : /**                                                                **/
     429                 :            : /********************************************************************/
     430                 :            : GEN
     431                 :        133 : mathilbert(long n) /* Hilbert matrix of order n */
     432                 :            : {
     433                 :            :   long i,j;
     434                 :            :   GEN p;
     435                 :            : 
     436         [ -  + ]:        133 :   if (n < 0) pari_err_DOMAIN("mathilbert", "dimension", "<", gen_0, stoi(n));
     437                 :        133 :   p = cgetg(n+1,t_MAT);
     438         [ +  + ]:       1120 :   for (j=1; j<=n; j++)
     439                 :            :   {
     440                 :        987 :     gel(p,j) = cgetg(n+1,t_COL);
     441 [ +  + ][ +  + ]:      16583 :     for (i=1+(j==1); i<=n; i++)
     442                 :      15596 :       gcoeff(p,i,j) = mkfrac(gen_1, utoipos(i+j-1));
     443                 :            :   }
     444         [ +  - ]:        133 :   if (n) gcoeff(p,1,1) = gen_1;
     445                 :        133 :   return p;
     446                 :            : }
     447                 :            : 
     448                 :            : /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
     449                 :            : GEN
     450                 :       1057 : matqpascal(long n, GEN q)
     451                 :            : {
     452                 :            :   long i, j, I;
     453                 :       1057 :   pari_sp av = avma;
     454                 :       1057 :   GEN m, qpow = NULL; /* gcc -Wall */
     455                 :            : 
     456         [ -  + ]:       1057 :   if (n < -1)  pari_err_DOMAIN("matpascal", "n", "<", gen_m1, stoi(n));
     457                 :       1057 :   n++; m = cgetg(n+1,t_MAT);
     458         [ +  + ]:      22120 :   for (j=1; j<=n; j++) gel(m,j) = cgetg(n+1,t_COL);
     459         [ +  + ]:       1057 :   if (q)
     460                 :            :   {
     461                 :         42 :     I = (n+1)/2;
     462         [ +  - ]:         42 :     if (I > 1) { qpow = new_chunk(I+1); gel(qpow,2)=q; }
     463         [ +  + ]:         84 :     for (j=3; j<=I; j++) gel(qpow,j) = gmul(q, gel(qpow,j-1));
     464                 :            :   }
     465         [ +  + ]:      22120 :   for (i=1; i<=n; i++)
     466                 :            :   {
     467                 :      21063 :     I = (i+1)/2; gcoeff(m,i,1)= gen_1;
     468         [ +  + ]:      21063 :     if (q)
     469                 :            :     {
     470         [ +  + ]:        483 :       for (j=2; j<=I; j++)
     471                 :        238 :         gcoeff(m,i,j) = gadd(gmul(gel(qpow,j),gcoeff(m,i-1,j)),
     472                 :        238 :                              gcoeff(m,i-1,j-1));
     473                 :            :     }
     474                 :            :     else
     475                 :            :     {
     476         [ +  + ]:    1023463 :       for (j=2; j<=I; j++)
     477                 :    1002645 :         gcoeff(m,i,j) = addii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
     478                 :            :     }
     479         [ +  + ]:    1034306 :     for (   ; j<=i; j++) gcoeff(m,i,j) = gcoeff(m,i,i+1-j);
     480         [ +  + ]:    2037189 :     for (   ; j<=n; j++) gcoeff(m,i,j) = gen_0;
     481                 :            :   }
     482                 :       1057 :   return gerepilecopy(av, m);
     483                 :            : }
     484                 :            : 
     485                 :            : /******************************************************************/
     486                 :            : /**                                                              **/
     487                 :            : /**                       PRECISION CHANGES                      **/
     488                 :            : /**                                                              **/
     489                 :            : /******************************************************************/
     490                 :            : 
     491                 :            : GEN
     492                 :        497 : gprec(GEN x, long l)
     493                 :            : {
     494                 :            :   long lx, i;
     495                 :            :   GEN y;
     496                 :            : 
     497         [ -  + ]:        497 :   if (l <= 0) pari_err_DOMAIN("gprec", "precision", "<=", gen_0, stoi(l));
     498   [ +  +  +  +  :        497 :   switch(typ(x))
                +  +  + ]
     499                 :            :   {
     500                 :            :     case t_REAL:
     501                 :         21 :       return rtor(x, ndec2prec(l));
     502                 :            :     case t_COMPLEX:
     503                 :         14 :       y = cgetg(3, t_COMPLEX);
     504                 :         14 :       gel(y,1) = gprec(gel(x,1),l);
     505                 :         14 :       gel(y,2) = gprec(gel(x,2),l);
     506                 :         14 :       break;
     507                 :            :     case t_PADIC:
     508         [ +  + ]:        343 :       if (!signe(gel(x,4))) return zeropadic(gel(x,2), l+precp(x));
     509                 :        336 :       y=cgetg(5,t_PADIC);
     510                 :        336 :       y[1]=x[1]; setprecp(y,l);
     511                 :        336 :       gel(y,2) = icopy(gel(x,2));
     512                 :        336 :       gel(y,3) = powiu(gel(x,2),l);
     513                 :        336 :       gel(y,4) = modii(gel(x,4), gel(y,3));
     514                 :        336 :       break;
     515                 :            : 
     516                 :            :     case t_SER:
     517         [ -  + ]:         28 :       if (lg(x) == 2) return zeroser(varn(x), l);
     518                 :         28 :       y=cgetg(l+2,t_SER); y[1]=x[1]; l++; i=l;
     519                 :         28 :       lx = lg(x);
     520         [ +  - ]:         28 :       if (l>=lx)
     521         [ +  + ]:        350 :         for ( ; i>=lx; i--) gel(y,i) = gen_0;
     522         [ +  + ]:        196 :       for ( ; i>=2; i--) gel(y,i) = gcopy(gel(x,i));
     523                 :         28 :       break;
     524                 :            :    case t_POL:
     525                 :         21 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     526         [ +  + ]:         63 :       for (i=2; i<lx; i++) gel(y,i) = gprec(gel(x,i),l);
     527                 :         21 :       break;
     528                 :            :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     529                 :         21 :       y = cgetg_copy(x, &lx);
     530         [ +  + ]:         63 :       for (i=1; i<lx; i++) gel(y,i) = gprec(gel(x,i),l);
     531                 :         21 :       break;
     532                 :         49 :     default: y = gcopy(x);
     533                 :            :   }
     534                 :        497 :   return y;
     535                 :            : }
     536                 :            : 
     537                 :            : /* internal: precision given in word length (including codewords) */
     538                 :            : GEN
     539                 :    2076267 : gprec_w(GEN x, long pr)
     540                 :            : {
     541                 :            :   long lx, i;
     542                 :            :   GEN y;
     543                 :            : 
     544   [ +  +  +  +  :    2076267 :   switch(typ(x))
                      + ]
     545                 :            :   {
     546                 :            :     case t_REAL:
     547         [ +  + ]:    1742581 :       if (signe(x)) return rtor(x,pr);
     548                 :      11909 :       i = -prec2nbits(pr);
     549         [ +  + ]:      11909 :       if (i < expo(x)) return real_0_bit(i);
     550                 :      11423 :       y = cgetr(2); y[1] = x[1]; return y;
     551                 :            :     case t_COMPLEX:
     552                 :     165023 :       y = cgetg(3, t_COMPLEX);
     553                 :     165023 :       gel(y,1) = gprec_w(gel(x,1),pr);
     554                 :     165023 :       gel(y,2) = gprec_w(gel(x,2),pr);
     555                 :     165023 :       break;
     556                 :            :    case t_POL: case t_SER:
     557                 :      45363 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     558         [ +  + ]:     354277 :       for (i=2; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     559                 :      45363 :       break;
     560                 :            :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     561                 :      26118 :       y = cgetg_copy(x, &lx);
     562         [ +  + ]:     151340 :       for (i=1; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
     563                 :      26118 :       break;
     564                 :      97182 :     default: return x;
     565                 :            :   }
     566                 :    2076267 :   return y;
     567                 :            : }
     568                 :            : 
     569                 :            : /* internal: precision given in word length (including codewords), truncate
     570                 :            :  * mantissa to precision 'pr' but never _increase_ it */
     571                 :            : GEN
     572                 :      22302 : gprec_wtrunc(GEN x, long pr)
     573                 :            : {
     574                 :            :   long lx, i;
     575                 :            :   GEN y;
     576                 :            : 
     577   [ +  +  -  +  :      22302 :   switch(typ(x))
                      + ]
     578                 :            :   {
     579                 :            :     case t_REAL:
     580 [ +  + ][ +  + ]:      17073 :       return (signe(x) && realprec(x) > pr)? rtor(x,pr): x;
     581                 :            :     case t_COMPLEX:
     582                 :       5208 :       y = cgetg(3, t_COMPLEX);
     583                 :       5208 :       gel(y,1) = gprec_wtrunc(gel(x,1),pr);
     584                 :       5208 :       gel(y,2) = gprec_wtrunc(gel(x,2),pr);
     585                 :       5208 :       break;
     586                 :            :     case t_POL:
     587                 :            :     case t_SER:
     588                 :          0 :       y = cgetg_copy(x, &lx); y[1] = x[1];
     589         [ #  # ]:          0 :       for (i=2; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     590                 :          0 :       break;
     591                 :            :     case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
     592                 :          7 :       y = cgetg_copy(x, &lx);
     593         [ +  + ]:        385 :       for (i=1; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
     594                 :          7 :       break;
     595                 :         14 :     default: return x;
     596                 :            :   }
     597                 :      22302 :   return y;
     598                 :            : }
     599                 :            : 
     600                 :            : /********************************************************************/
     601                 :            : /**                                                                **/
     602                 :            : /**                      SERIES TRANSFORMS                         **/
     603                 :            : /**                                                                **/
     604                 :            : /********************************************************************/
     605                 :            : /**                  LAPLACE TRANSFORM (OF A SERIES)               **/
     606                 :            : /********************************************************************/
     607                 :            : static GEN
     608                 :          7 : serlaplace(GEN x)
     609                 :            : {
     610                 :          7 :   long i, l = lg(x), e = valp(x);
     611                 :          7 :   GEN t, y = cgetg(l,t_SER);
     612         [ -  + ]:          7 :   if (e < 0) pari_err_DOMAIN("laplace","valuation","<",gen_0,stoi(e));
     613                 :          7 :   t = mpfact(e); y[1] = x[1];
     614         [ +  + ]:        119 :   for (i=2; i<l; i++)
     615                 :            :   {
     616                 :        112 :     gel(y,i) = gmul(t, gel(x,i));
     617                 :        112 :     e++; t = mului(e,t);
     618                 :            :   }
     619                 :          7 :   return y;
     620                 :            : }
     621                 :            : static GEN
     622                 :         14 : pollaplace(GEN x)
     623                 :            : {
     624                 :         14 :   long i, e = 0, l = lg(x);
     625                 :         14 :   GEN t = gen_1, y = cgetg(l,t_POL);
     626                 :         14 :   y[1] = x[1];
     627         [ +  + ]:         63 :   for (i=2; i<l; i++)
     628                 :            :   {
     629                 :         49 :     gel(y,i) = gmul(t, gel(x,i));
     630                 :         49 :     e++; t = mului(e,t);
     631                 :            :   }
     632                 :         14 :   return y;
     633                 :            : }
     634                 :            : GEN
     635                 :         21 : laplace(GEN x)
     636                 :            : {
     637                 :         21 :   pari_sp av = avma;
     638      [ +  +  - ]:         21 :   switch(typ(x))
     639                 :            :   {
     640                 :         14 :     case t_POL: x = pollaplace(x); break;
     641                 :          7 :     case t_SER: x = serlaplace(x); break;
     642                 :          0 :     default: pari_err_TYPE("laplace",x);
     643                 :            :   }
     644                 :         21 :   return gerepilecopy(av, x);
     645                 :            : }
     646                 :            : 
     647                 :            : /********************************************************************/
     648                 :            : /**              CONVOLUTION PRODUCT (OF TWO SERIES)               **/
     649                 :            : /********************************************************************/
     650                 :            : GEN
     651                 :          7 : convol(GEN x, GEN y)
     652                 :            : {
     653                 :          7 :   long j, lx, ly, ex, ey, vx = varn(x);
     654                 :            :   GEN z;
     655                 :            : 
     656         [ -  + ]:          7 :   if (typ(x) != t_SER) pari_err_TYPE("convol",x);
     657         [ -  + ]:          7 :   if (typ(y) != t_SER) pari_err_TYPE("convol",y);
     658         [ -  + ]:          7 :   if (varn(y) != vx) pari_err_VAR("convol", x,y);
     659                 :          7 :   ex = valp(x);
     660                 :          7 :   ey = valp(y);
     661         [ -  + ]:          7 :   if (ser_isexactzero(x))
     662                 :          0 :     return scalarser(gadd(RgX_get_0(x), RgX_get_0(y)), varn(x), maxss(ex,ey));
     663                 :          7 :   lx = lg(x) + ex; x -= ex;
     664                 :          7 :   ly = lg(y) + ey; y -= ey;
     665                 :            :   /* inputs shifted: x[i] and y[i] now correspond to monomials of same degree */
     666         [ -  + ]:          7 :   if (ly < lx) lx = ly; /* min length */
     667         [ -  + ]:          7 :   if (ex < ey) ex = ey; /* max valuation */
     668         [ -  + ]:          7 :   if (lx - ex < 3) return zeroser(vx, lx-2);
     669                 :            : 
     670                 :          7 :   z = cgetg(lx - ex, t_SER);
     671                 :          7 :   z[1] = evalvalp(ex) | evalvarn(vx);
     672         [ +  + ]:        119 :   for (j = ex+2; j<lx; j++) gel(z,j-ex) = gmul(gel(x,j),gel(y,j));
     673                 :          7 :   return normalize(z);
     674                 :            : }
     675                 :            : 
     676                 :            : /***********************************************************************/
     677                 :            : /*               OPERATIONS ON DIRICHLET SERIES: *, /                  */
     678                 :            : /* (+, -, scalar multiplication are done on the corresponding vectors) */
     679                 :            : /***********************************************************************/
     680                 :            : static long
     681                 :       1302 : dirval(GEN x)
     682                 :            : {
     683                 :       1302 :   long i = 1, lx = lg(x);
     684 [ +  - ][ +  + ]:       1323 :   while (i < lx && gequal0(gel(x,i))) i++;
     685                 :       1302 :   return i;
     686                 :            : }
     687                 :            : 
     688                 :            : GEN
     689                 :        133 : dirmul(GEN x, GEN y)
     690                 :            : {
     691                 :        133 :   pari_sp av = avma, av2;
     692                 :            :   long nx, ny, nz, dx, dy, i, j, k;
     693                 :            :   GEN z;
     694                 :            : 
     695         [ -  + ]:        133 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirmul",x);
     696         [ -  + ]:        133 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirmul",y);
     697                 :        133 :   dx = dirval(x); nx = lg(x)-1;
     698                 :        133 :   dy = dirval(y); ny = lg(y)-1;
     699         [ -  + ]:        133 :   if (ny-dy < nx-dx) { swap(x,y); lswap(nx,ny); lswap(dx,dy); }
     700                 :        133 :   nz = minss(nx*dy,ny*dx);
     701                 :        133 :   y = RgV_kill0(y);
     702                 :        133 :   av2 = avma;
     703                 :        133 :   z = zerovec(nz);
     704         [ +  + ]:       6454 :   for (j=dx; j<=nx; j++)
     705                 :            :   {
     706                 :       6321 :     GEN c = gel(x,j);
     707         [ +  + ]:       6321 :     if (gequal0(c)) continue;
     708         [ +  + ]:       3101 :     if (gequal1(c))
     709                 :            :     {
     710         [ +  + ]:      13027 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     711         [ +  + ]:      11522 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gel(y,k));
     712                 :            :     }
     713         [ +  + ]:       1596 :     else if (gequalm1(c))
     714                 :            :     {
     715         [ +  + ]:       4669 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     716         [ +  + ]:       3563 :         if (gel(y,k)) gel(z,i) = gsub(gel(z,i),gel(y,k));
     717                 :            :     }
     718                 :            :     else
     719                 :            :     {
     720         [ +  + ]:       2233 :       for (k=dy,i=j*dy; i<=nz; i+=j,k++)
     721         [ +  + ]:       1743 :         if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gmul(c,gel(y,k)));
     722                 :            :     }
     723         [ -  + ]:       3101 :     if (gc_needed(av2,3))
     724                 :            :     {
     725         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirmul, %ld/%ld",j,nx);
     726                 :          0 :       z = gerepilecopy(av2,z);
     727                 :            :     }
     728                 :            :   }
     729                 :        133 :   return gerepilecopy(av,z);
     730                 :            : }
     731                 :            : 
     732                 :            : GEN
     733                 :        518 : dirdiv(GEN x, GEN y)
     734                 :            : {
     735                 :        518 :   pari_sp av = avma, av2;
     736                 :            :   long nx,ny,nz, dx,dy, i,j,k;
     737                 :            :   GEN p1;
     738                 :            : 
     739         [ -  + ]:        518 :   if (typ(x)!=t_VEC) pari_err_TYPE("dirdiv",x);
     740         [ -  + ]:        518 :   if (typ(y)!=t_VEC) pari_err_TYPE("dirdiv",y);
     741                 :        518 :   dx = dirval(x); nx = lg(x)-1;
     742                 :        518 :   dy = dirval(y); ny = lg(y)-1;
     743 [ +  - ][ -  + ]:        518 :   if (dy != 1 || !ny) pari_err_INV("dirdiv",y);
     744                 :        518 :   nz = minss(nx,ny*dx);
     745                 :        518 :   p1 = gel(y,1);
     746         [ +  + ]:        518 :   if (gequal1(p1)) p1 = NULL; else y = gdiv(y,p1);
     747                 :        518 :   y = RgV_kill0(y);
     748                 :        518 :   av2 = avma;
     749         [ +  + ]:        518 :   x = p1 ? gdiv(x,p1): leafcopy(x);
     750         [ +  + ]:        525 :   for (j=1; j<dx; j++) gel(x,j) = gen_0;
     751                 :        518 :   setlg(x,nz+1);
     752         [ +  + ]:     514486 :   for (j=dx; j<=nz; j++)
     753                 :            :   {
     754                 :     513968 :     GEN c = gel(x,j);
     755         [ +  + ]:     513968 :     if (gequal0(c)) continue;
     756         [ +  + ]:     237489 :     if (gequal1(c))
     757                 :            :     {
     758         [ +  + ]:    1446046 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     759         [ +  + ]:    1344371 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gel(y,k));
     760                 :            :     }
     761         [ +  + ]:     135814 :     else if (gequalm1(c))
     762                 :            :     {
     763         [ +  + ]:    1321964 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     764         [ +  + ]:    1217314 :         if (gel(y,k)) gel(x,i) = gadd(gel(x,i),gel(y,k));
     765                 :            :     }
     766                 :            :     else
     767                 :            :     {
     768         [ +  + ]:     132888 :       for (i=j+j,k=2; i<=nz; i+=j,k++)
     769         [ +  + ]:     101724 :         if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gmul(c,gel(y,k)));
     770                 :            :     }
     771         [ -  + ]:     237489 :     if (gc_needed(av2,3))
     772                 :            :     {
     773         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"dirdiv, %ld/%ld",j,nz);
     774                 :          0 :       x = gerepilecopy(av2,x);
     775                 :            :     }
     776                 :            :   }
     777                 :        518 :   return gerepilecopy(av,x);
     778                 :            : }
     779                 :            : 
     780                 :            : /*******************************************************************/
     781                 :            : /**                                                               **/
     782                 :            : /**                       COMBINATORICS                           **/
     783                 :            : /**                                                               **/
     784                 :            : /*******************************************************************/
     785                 :            : /**                      BINOMIAL COEFFICIENTS                    **/
     786                 :            : /*******************************************************************/
     787                 :            : GEN
     788                 :      73038 : binomialuu(ulong n, ulong k)
     789                 :            : {
     790                 :      73038 :   pari_sp ltop = avma;
     791                 :            :   GEN z;
     792         [ +  + ]:      73038 :   if (k > n) return gen_0;
     793                 :      73031 :   k = minuu(k,n-k);
     794         [ +  + ]:      73031 :   if (!k) return gen_1;
     795         [ +  + ]:      64428 :   if (k == 1) return utoipos(n);
     796                 :      58709 :   z = diviiexact(mulu_interval(n-k+1, n), mulu_interval(2UL, k));
     797                 :      73038 :   return gerepileuptoint(ltop,z);
     798                 :            : }
     799                 :            : 
     800                 :            : GEN
     801                 :      99967 : binomial(GEN n, long k)
     802                 :            : {
     803                 :            :   long i, prec;
     804                 :            :   pari_sp av;
     805                 :            :   GEN y;
     806                 :            : 
     807         [ +  + ]:      99967 :   if (k <= 1)
     808                 :            :   {
     809         [ -  + ]:      59339 :     if (is_noncalc_t(typ(n))) pari_err_TYPE("binomial",n);
     810         [ -  + ]:      59339 :     if (k < 0) return gen_0;
     811         [ +  + ]:      59339 :     if (k == 0) return gen_1;
     812                 :      26019 :     return gcopy(n);
     813                 :            :   }
     814                 :      40628 :   av = avma;
     815         [ +  + ]:      40628 :   if (typ(n) == t_INT)
     816                 :            :   {
     817         [ +  - ]:      40614 :     if (signe(n) > 0)
     818                 :            :     {
     819                 :      40614 :       GEN z = subis(n,k);
     820         [ +  + ]:      40614 :       if (cmpis(z,k) < 0)
     821                 :            :       {
     822                 :        889 :         k = itos(z); avma = av;
     823         [ +  + ]:        889 :         if (k <= 1)
     824                 :            :         {
     825         [ -  + ]:        336 :           if (k < 0) return gen_0;
     826         [ -  + ]:        336 :           if (k == 0) return gen_1;
     827                 :        336 :           return icopy(n);
     828                 :            :         }
     829                 :            :       }
     830                 :            :     }
     831                 :            :     /* k > 1 */
     832 [ +  + ][ +  - ]:      40278 :     if (lgefint(n) == 3 && signe(n) > 0)
     833                 :            :     {
     834                 :      40271 :       y = binomialuu(itou(n),(ulong)k);
     835                 :      40271 :       return gerepileupto(av, y);
     836                 :            :     }
     837                 :            :     else
     838                 :            :     {
     839                 :          7 :       y = cgetg(k+1,t_VEC);
     840         [ +  + ]:         21 :       for (i=1; i<=k; i++) gel(y,i) = subis(n,i-1);
     841                 :          7 :       y = ZV_prod(y);
     842                 :            :     }
     843                 :          7 :     y = diviiexact(y, mpfact(k));
     844                 :          7 :     return gerepileuptoint(av, y);
     845                 :            :   }
     846                 :            : 
     847                 :         14 :   prec = precision(n);
     848 [ +  - ][ +  + ]:         14 :   if (prec && k > 200 + 0.8*prec2nbits(prec)) {
     849                 :          7 :     GEN A = mpfactr(k, prec), B = ggamma(gsubgs(n,k-1), prec);
     850                 :          7 :     return gerepileupto(av, gdiv(ggamma(gaddgs(n,1), prec), gmul(A,B)));
     851                 :            :   }
     852                 :            : 
     853                 :          7 :   y = cgetg(k+1,t_VEC);
     854         [ +  + ]:         42 :   for (i=1; i<=k; i++) gel(y,i) = gsubgs(n,i-1);
     855                 :      99967 :   return gerepileupto(av, gdiv(RgV_prod(y), mpfact(k)));
     856                 :            : }
     857                 :            : 
     858                 :            : /* Assume n >= 0, return bin, bin[k+1] = binomial(n, k) */
     859                 :            : GEN
     860                 :      13993 : vecbinome(long n)
     861                 :            : {
     862                 :            :   long d, k;
     863                 :            :   GEN C;
     864         [ +  + ]:      13993 :   if (!n) return mkvec(gen_1);
     865                 :      13706 :   C = cgetg(n+2, t_VEC) + 1; /* C[k] = binomial(n, k) */
     866                 :      13706 :   gel(C,0) = gen_1;
     867                 :      13706 :   gel(C,1) = utoipos(n); d = (n + 1) >> 1;
     868         [ +  + ]:     103082 :   for (k=2; k <= d; k++)
     869                 :            :   {
     870                 :      89376 :     pari_sp av = avma;
     871                 :      89376 :     gel(C,k) = gerepileuptoint(av, diviuexact(mului(n-k+1, gel(C,k-1)), k));
     872                 :            :   }
     873         [ +  + ]:     110019 :   for (   ; k <= n; k++) gel(C,k) = gel(C,n-k);
     874                 :      13993 :   return C - 1;
     875                 :            : }
     876                 :            : 
     877                 :            : /********************************************************************/
     878                 :            : /**                  STIRLING NUMBERS                              **/
     879                 :            : /********************************************************************/
     880                 :            : /* Stirling number of the 2nd kind. The number of ways of partitioning
     881                 :            :    a set of n elements into m non-empty subsets. */
     882                 :            : GEN
     883                 :       1694 : stirling2(ulong n, ulong m)
     884                 :            : {
     885                 :       1694 :   pari_sp av = avma;
     886                 :            :   GEN s, bmk;
     887                 :            :   ulong k;
     888 [ -  + ][ #  # ]:       1694 :   if (n==0) return (m == 0)? gen_1: gen_0;
     889 [ +  - ][ -  + ]:       1694 :   if (m > n || m == 0) return gen_0;
     890         [ -  + ]:       1694 :   if (m==n) return gen_1;
     891                 :            :   /* k = 0 */
     892                 :       1694 :   bmk = gen_1; s  = powuu(m, n);
     893         [ +  + ]:      20314 :   for (k = 1; k <= ((m-1)>>1); ++k)
     894                 :            :   { /* bmk = binomial(m, k) */
     895                 :            :     GEN c, kn, mkn;
     896                 :      18620 :     bmk = diviuexact(mului(m-k+1, bmk), k);
     897                 :      18620 :     kn  = powuu(k, n); mkn = powuu(m-k, n);
     898         [ +  + ]:      18620 :     c = odd(m)? subii(mkn,kn): addii(mkn,kn);
     899                 :      18620 :     c = mulii(bmk, c);
     900         [ +  + ]:      18620 :     s = odd(k)? subii(s, c): addii(s, c);
     901         [ -  + ]:      18620 :     if (gc_needed(av,2))
     902                 :            :     {
     903         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"stirling2");
     904                 :          0 :       gerepileall(av, 2, &s, &bmk);
     905                 :            :     }
     906                 :            :   }
     907                 :            :   /* k = m/2 */
     908         [ +  + ]:       1694 :   if (!odd(m))
     909                 :            :   {
     910                 :            :     GEN c;
     911                 :        805 :     bmk = diviuexact(mului(k+1, bmk), k);
     912                 :        805 :     c = mulii(bmk, powuu(k,n));
     913         [ +  + ]:        805 :     s = odd(k)? subii(s, c): addii(s, c);
     914                 :            :   }
     915                 :       1694 :   return gerepileuptoint(av, diviiexact(s, mpfact(m)));
     916                 :            : }
     917                 :            : 
     918                 :            : /* Stirling number of the first kind. Up to the sign, the number of
     919                 :            :    permutations of n symbols which have exactly m cycles. */
     920                 :            : GEN
     921                 :        154 : stirling1(ulong n, ulong m)
     922                 :            : {
     923                 :        154 :   pari_sp ltop=avma;
     924                 :            :   ulong k;
     925                 :            :   GEN s, t;
     926         [ -  + ]:        154 :   if (n < m) return gen_0;
     927         [ -  + ]:        154 :   else if (n==m) return gen_1;
     928                 :            :   /* t = binomial(n-1+k, m-1) * binomial(2n-m, n-m-k) */
     929                 :            :   /* k = n-m > 0 */
     930                 :        154 :   t = binomialuu(2*n-m-1, m-1);
     931                 :        154 :   s = mulii(t, stirling2(2*(n-m), n-m));
     932         [ +  + ]:        154 :   if (odd(n-m)) togglesign(s);
     933         [ +  + ]:       1547 :   for (k = n-m-1; k > 0; --k)
     934                 :            :   {
     935                 :            :     GEN c;
     936                 :       1393 :     t = diviuuexact(muluui(n-m+k+1, n+k+1, t), n+k, n-m-k);
     937                 :       1393 :     c = mulii(t, stirling2(n-m+k, k));
     938         [ +  + ]:       1393 :     s = odd(k)? subii(s, c): addii(s, c);
     939         [ +  + ]:       1393 :     if ((k & 0x1f) == 0) {
     940                 :         21 :       t = gerepileuptoint(ltop, t);
     941                 :         21 :       s = gerepileuptoint(avma, s);
     942                 :            :     }
     943                 :            :   }
     944                 :        154 :   return gerepileuptoint(ltop, s);
     945                 :            : }
     946                 :            : 
     947                 :            : GEN
     948                 :        301 : stirling(long n, long m, long flag)
     949                 :            : {
     950         [ -  + ]:        301 :   if (n < 0) pari_err_DOMAIN("stirling", "n", "<", gen_0, stoi(n));
     951         [ -  + ]:        301 :   if (m < 0) pari_err_DOMAIN("stirling", "m", "<", gen_0, stoi(m));
     952      [ +  +  - ]:        301 :   switch (flag)
     953                 :            :   {
     954                 :        154 :     case 1: return stirling1((ulong)n,(ulong)m);
     955                 :        147 :     case 2: return stirling2((ulong)n,(ulong)m);
     956                 :          0 :     default: pari_err_FLAG("stirling");
     957                 :            :   }
     958                 :        301 :   return NULL; /*NOT REACHED*/
     959                 :            : }
     960                 :            : 
     961                 :            : /***********************************************************************/
     962                 :            : /**                          PERMUTATIONS                             **/
     963                 :            : /***********************************************************************/
     964                 :            : GEN
     965                 :       9485 : numtoperm(long n, GEN x)
     966                 :            : {
     967                 :            :   pari_sp av;
     968                 :            :   ulong i, r;
     969                 :            :   GEN v;
     970                 :            : 
     971         [ -  + ]:       9485 :   if (n < 0) pari_err_DOMAIN("numtoperm", "n", "<", gen_0, stoi(n));
     972         [ -  + ]:       9485 :   if (typ(x) != t_INT) pari_err_TYPE("numtoperm",x);
     973         [ -  + ]:       9485 :   v = cgetg(n+1, t_VEC); if (n==0) return v;
     974                 :       9485 :   uel(v,n) = 1; av = avma;
     975         [ +  + ]:       9485 :   if (signe(x) <= 0) x = modii(x, mpfact(n));
     976         [ +  + ]:      38143 :   for (r=n-1; r>=1; r--)
     977                 :            :   {
     978                 :            :     ulong a;
     979                 :      28658 :     x = diviu_rem(x, n+1-r,&a);
     980         [ +  + ]:      87157 :     for (i=r+1; i<=(ulong)n; i++)
     981         [ +  + ]:      58499 :       if (uel(v,i) > a) uel(v,i)++;
     982                 :      28658 :     uel(v,r) = a+1;
     983         [ -  + ]:      28658 :     if (gc_needed(av,2))
     984                 :          0 :       x = gerepileuptoint(av, x);
     985                 :            :   }
     986                 :       9485 :   avma = av;
     987         [ +  + ]:      47628 :   for (i=1; i<=(ulong)n; i++) gel(v,i) = utoipos(uel(v,i));
     988                 :       9485 :   return v;
     989                 :            : }
     990                 :            : 
     991                 :            : GEN
     992                 :       5166 : permtonum(GEN p)
     993                 :            : {
     994                 :       5166 :   long n = lg(p)-1, i, r;
     995                 :       5166 :   pari_sp av = avma, av2;
     996                 :            :   GEN v, x;
     997                 :            : 
     998         [ -  + ]:       5166 :   if (!is_vec_t(typ(p))) pari_err_TYPE("permtonum",p);
     999                 :       5166 :   v = cgetg(n+1,t_VECSMALL);
    1000         [ +  + ]:      26362 :   for (i=1; i<=n; i++)
    1001                 :            :   {
    1002                 :      21196 :     GEN pi = gel(p, i);
    1003         [ -  + ]:      21196 :     if (typ(pi) != t_INT) pari_err_TYPE("permtonum",pi);
    1004                 :      21196 :     v[i] = itos(pi);
    1005                 :            :   }
    1006                 :       5166 :   x = gen_0; av2 = avma;
    1007         [ +  + ]:      26362 :   for (i=1; i<=n; i++)
    1008                 :            :   {
    1009                 :      21196 :     long vi = v[i];
    1010         [ +  + ]:      21196 :     x = i==1 ? stoi(v[1]-1): addiu(mulis(x,n+1-i),vi-1);
    1011         [ +  + ]:      54719 :     for (r=i+1; r<=n; r++)
    1012         [ +  + ]:      33523 :       if (v[r]>vi) v[r]--;
    1013         [ -  + ]:      21196 :     if (gc_needed(av,2))
    1014                 :          0 :       x = gerepileuptoint(av2, x);
    1015                 :            :   }
    1016                 :       5166 :   return gerepileuptoint(av, x);
    1017                 :            : }
    1018                 :            : 
    1019                 :            : /*******************************************************************/
    1020                 :            : /**                                                               **/
    1021                 :            : /**                     RECIPROCAL POLYNOMIAL                     **/
    1022                 :            : /**                                                               **/
    1023                 :            : /*******************************************************************/
    1024                 :            : /* return coefficients s.t x = x_0 X^n + ... + x_n */
    1025                 :            : GEN
    1026                 :      62118 : polrecip(GEN x)
    1027                 :            : {
    1028         [ -  + ]:      62118 :   if (typ(x) != t_POL) pari_err_TYPE("polrecip",x);
    1029                 :      62118 :   return RgX_recip(x);
    1030                 :            : }
    1031                 :            : 
    1032                 :            : /********************************************************************/
    1033                 :            : /**                                                                **/
    1034                 :            : /**                  POLYNOMIAL INTERPOLATION                      **/
    1035                 :            : /**                                                                **/
    1036                 :            : /********************************************************************/
    1037                 :            : /* allow X = NULL for [1,...,n] */
    1038                 :            : GEN
    1039                 :         28 : RgV_polint(GEN X, GEN Y, long v)
    1040                 :            : {
    1041                 :         28 :   pari_sp av0 = avma, av;
    1042                 :         28 :   GEN Q, P = NULL;
    1043                 :         28 :   long i, l = lg(Y);
    1044         [ +  + ]:         28 :   if (!X)
    1045                 :            :   {
    1046                 :         14 :     X = cgetg(l, t_VEC);
    1047         [ +  + ]:         56 :     for (i=1; i<l; i++) gel(X,i) = utoipos(i);
    1048                 :            :   }
    1049                 :         28 :   Q = roots_to_pol(X, v); av = avma;
    1050         [ +  + ]:        105 :   for (i=1; i<l; i++)
    1051                 :            :   {
    1052                 :            :     GEN inv, T, dP;
    1053         [ +  + ]:         77 :     if (gequal0(gel(Y,i))) continue;
    1054                 :         63 :     T = RgX_div_by_X_x(Q, gel(X,i), NULL);
    1055                 :         63 :     inv = ginv(poleval(T,gel(X,i)));
    1056                 :         63 :     dP = RgX_Rg_mul(T, gmul(gel(Y,i),inv));
    1057         [ +  + ]:         63 :     P = P? RgX_add(P, dP): dP;
    1058         [ -  + ]:         63 :     if (gc_needed(av,2))
    1059                 :            :     {
    1060         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpV_polint");
    1061                 :          0 :       P = gerepileupto(av, P);
    1062                 :            :     }
    1063                 :            :   }
    1064         [ +  + ]:         28 :   if (!P) { avma = av; return zeropol(v); }
    1065                 :         28 :   return gerepileupto(av0, P);
    1066                 :            : }
    1067                 :            : /* X,Y are "spec" GEN vectors with n > 1 components ( at X[0], ... X[n-1] ) */
    1068                 :            : GEN
    1069                 :        417 : polint_i(GEN X, GEN Y, GEN x, long n, GEN *ptdy)
    1070                 :            : {
    1071                 :        417 :   long i, m, ns = 0;
    1072                 :        417 :   pari_sp av = avma;
    1073                 :        417 :   GEN y, c, d, dy = NULL; /* gcc -Wall */
    1074                 :            : 
    1075         [ +  + ]:        417 :   if (!X)
    1076                 :            :   {
    1077                 :          7 :     X = cgetg(n+1, t_VEC);
    1078         [ +  + ]:         28 :     for (i=1; i<=n; i++) gel(X,i) = utoipos(i);
    1079                 :          7 :     X++;
    1080                 :            :   }
    1081         [ +  - ]:        417 :   switch(typ(x)) {
    1082                 :            :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1083                 :            :     {
    1084                 :        417 :       GEN D = NULL;
    1085         [ +  + ]:      16888 :       for (i=0; i<n; i++)
    1086                 :            :       {
    1087                 :      16478 :         GEN t = gsub(x,gel(X,i));
    1088         [ +  + ]:      16478 :         switch(typ(t))
    1089                 :            :         {
    1090                 :            :           case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1091                 :      16471 :             t = gabs(t, DEFAULTPREC);
    1092 [ +  + ][ +  + ]:      16471 :             if (!D || gcmp(t,D) < 0) { ns = i; D = t; }
    1093                 :      16471 :             break;
    1094                 :            :           default:
    1095                 :          7 :             goto NODY;
    1096                 :            :         }
    1097                 :            :       }
    1098                 :        410 :       break;
    1099                 :            :       /* X[ns] is closest to x */
    1100                 :            :     }
    1101                 :            : NODY:
    1102                 :            :     default:
    1103         [ -  + ]:          7 :       if (ptdy) { *ptdy = gen_0; ptdy = NULL; }
    1104                 :            :   }
    1105                 :        417 :   c = new_chunk(n);
    1106         [ +  + ]:      16909 :   d = new_chunk(n); for (i=0; i<n; i++) gel(c,i) = gel(d,i) = gel(Y,i);
    1107                 :        417 :   y = gel(d,ns--);
    1108                 :            :   /* divided differences */
    1109         [ +  + ]:      16485 :   for (m=1; m<n; m++)
    1110                 :            :   {
    1111         [ +  + ]:     455241 :     for (i=0; i<n-m; i++)
    1112                 :            :     {
    1113                 :     439173 :       GEN ho = gsub(gel(X,i),x), hp = gsub(gel(X,i+m),x), den = gsub(ho,hp);
    1114         [ +  + ]:     439173 :       if (gequal0(den))
    1115                 :            :       {
    1116                 :          7 :         char *x1 = stack_sprintf("X[%ld]", i+1);
    1117                 :          7 :         char *x2 = stack_sprintf("X[%ld]", i+m+1);
    1118                 :          7 :         pari_err_DOMAIN("polinterpolate",x1,"=",strtoGENstr(x2), X);
    1119                 :            :       }
    1120                 :     439166 :       den = gdiv(gsub(gel(c,i+1),gel(d,i)), den);
    1121                 :     439166 :       gel(c,i) = gmul(ho,den);
    1122                 :     439166 :       gel(d,i) = gmul(hp,den);
    1123                 :            :     }
    1124         [ +  + ]:      16068 :     dy = (2*(ns+1) < n-m)? gel(c,ns+1): gel(d,ns--);
    1125                 :      16068 :     y = gadd(y,dy);
    1126                 :            :   }
    1127         [ +  + ]:        410 :   if (!ptdy) return gerepileupto(av, y);
    1128                 :         98 :   *ptdy = dy;
    1129                 :         98 :   gerepileall(av, 2, &y, ptdy);
    1130                 :        410 :   return y;
    1131                 :            : }
    1132                 :            : 
    1133                 :            : GEN
    1134                 :        368 : polint(GEN X, GEN Y, GEN x, GEN *ptdy)
    1135                 :            : {
    1136                 :        368 :   long lx = lg(X);
    1137                 :            : 
    1138         [ -  + ]:        368 :   if (! is_vec_t(typ(X))) pari_err_TYPE("polinterpolate",X);
    1139         [ +  + ]:        368 :   if (Y)
    1140                 :            :   {
    1141         [ -  + ]:        347 :     if (! is_vec_t(typ(Y))) pari_err_TYPE("polinterpolate",Y);
    1142         [ -  + ]:        347 :     if (lx != lg(Y)) pari_err_DIM("polinterpolate");
    1143                 :            :   }
    1144                 :            :   else
    1145                 :            :   {
    1146                 :         21 :     Y = X;
    1147                 :         21 :     X = NULL;
    1148                 :            :   }
    1149                 :            : 
    1150         [ +  + ]:        368 :   if (lx <= 2)
    1151                 :            :   {
    1152         [ -  + ]:         14 :     if (ptdy) *ptdy = gen_0;
    1153         [ +  + ]:         14 :     if (lx == 1) return zeropol(0);
    1154                 :          7 :     Y = gel(Y,1);
    1155         [ -  + ]:          7 :     if (gvar(Y) == 0) pari_err_PRIORITY("polinterpolate", Y, "=", 0);
    1156                 :          7 :     return scalarpol(Y, 0);
    1157                 :            :   }
    1158         [ +  + ]:        354 :   if (!x) return RgV_polint(X, Y, 0);
    1159         [ -  + ]:        326 :   if (gequalX(x)) return RgV_polint(X, Y, varn(x));
    1160         [ +  + ]:        368 :   return polint_i(X? X+1: NULL,Y+1,x,lx-1,ptdy);
    1161                 :            : }
    1162                 :            : 
    1163                 :            : /********************************************************************/
    1164                 :            : /**                                                                **/
    1165                 :            : /**                       MODREVERSE                               **/
    1166                 :            : /**                                                                **/
    1167                 :            : /********************************************************************/
    1168                 :            : static void
    1169                 :          7 : err_reverse(GEN x, GEN T)
    1170                 :            : {
    1171                 :          7 :   pari_err_DOMAIN("modreverse","deg(minpoly(z))", "<", stoi(degpol(T)),
    1172                 :            :                   mkpolmod(x,T));
    1173                 :          0 : }
    1174                 :            : 
    1175                 :            : /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
    1176                 :            : GEN
    1177                 :         28 : RgXQ_reverse(GEN a, GEN T)
    1178                 :            : {
    1179                 :         28 :   pari_sp av = avma;
    1180                 :         28 :   long n = degpol(T);
    1181                 :            :   GEN y;
    1182                 :            : 
    1183         [ +  + ]:         28 :   if (n <= 1) {
    1184         [ -  + ]:          7 :     if (n <= 0) return gcopy(a);
    1185                 :          7 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1186                 :            :   }
    1187 [ +  - ][ -  + ]:         21 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1188                 :         21 :   y = RgXV_to_RgM(RgXQ_powers(a,n-1,T), n);
    1189                 :         21 :   y = RgM_solve(y, col_ei(n, 2));
    1190         [ +  + ]:         21 :   if (!y) err_reverse(a,T);
    1191                 :         21 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1192                 :            : }
    1193                 :            : GEN
    1194                 :        560 : QXQ_reverse(GEN a, GEN T)
    1195                 :            : {
    1196                 :        560 :   pari_sp av = avma;
    1197                 :        560 :   long n = degpol(T);
    1198                 :            :   GEN y;
    1199                 :            : 
    1200         [ +  + ]:        560 :   if (n <= 1) {
    1201         [ -  + ]:         14 :     if (n <= 0) return gcopy(a);
    1202                 :         14 :     return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
    1203                 :            :   }
    1204 [ +  - ][ -  + ]:        546 :   if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
    1205         [ -  + ]:        546 :   if (gequalX(a)) return gcopy(a);
    1206                 :        546 :   y = RgXV_to_RgM(QXQ_powers(a,n-1,T), n);
    1207                 :        546 :   y = RgM_solve(y, col_ei(n, 2));
    1208         [ -  + ]:        546 :   if (!y) err_reverse(a,T);
    1209                 :        560 :   return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
    1210                 :            : }
    1211                 :            : 
    1212                 :            : GEN
    1213                 :         28 : modreverse(GEN x)
    1214                 :            : {
    1215                 :            :   long v, n;
    1216                 :            :   GEN T, a, y;
    1217                 :            : 
    1218         [ -  + ]:         28 :   if (typ(x)!=t_POLMOD) pari_err_TYPE("modreverse",x);
    1219         [ +  + ]:         28 :   T = gel(x,1); n = degpol(T); if (n <= 0) return gcopy(x);
    1220                 :         21 :   a = gel(x,2);
    1221                 :         21 :   v = varn(T);
    1222                 :         21 :   y = cgetg(3,t_POLMOD);
    1223         [ +  + ]:         21 :   gel(y,1) = (n==1)? gsub(pol_x(v), a): RgXQ_charpoly(a, T, v);
    1224                 :         28 :   gel(y,2) = RgXQ_reverse(a, T); return y;
    1225                 :            : }
    1226                 :            : 
    1227                 :            : /********************************************************************/
    1228                 :            : /**                                                                **/
    1229                 :            : /**                          MERGESORT                             **/
    1230                 :            : /**                                                                **/
    1231                 :            : /********************************************************************/
    1232                 :            : static int
    1233                 :   28204637 : cmp_small(GEN x, GEN y) {
    1234                 :   28204637 :   long a = (long)x, b = (long)y;
    1235 [ +  + ][ +  + ]:   28204637 :   return a>b? 1: (a<b? -1: 0);
    1236                 :            : }
    1237                 :            : 
    1238                 :            : static int
    1239                 :       1120 : veccmp(void *data, GEN x, GEN y)
    1240                 :            : {
    1241                 :       1120 :   GEN k = (GEN)data;
    1242                 :       1120 :   long i, s, lk = lg(k), lx = minss(lg(x), lg(y));
    1243                 :            : 
    1244         [ -  + ]:       1120 :   if (!is_vec_t(typ(x))) pari_err_TYPE("lexicographic vecsort",x);
    1245         [ -  + ]:       1120 :   if (!is_vec_t(typ(y))) pari_err_TYPE("lexicographic vecsort",y);
    1246         [ +  + ]:       1197 :   for (i=1; i<lk; i++)
    1247                 :            :   {
    1248                 :       1148 :     long c = k[i];
    1249         [ +  + ]:       1148 :     if (c >= lx)
    1250                 :         14 :       pari_err_TYPE("lexicographic vecsort, index too large", stoi(c));
    1251                 :       1134 :     s = lexcmp(gel(x,c), gel(y,c));
    1252         [ +  + ]:       1134 :     if (s) return s;
    1253                 :            :   }
    1254                 :       1106 :   return 0;
    1255                 :            : }
    1256                 :            : 
    1257                 :            : /* return permutation sorting v[1..n], removing duplicates. Assume n > 0 */
    1258                 :            : static GEN
    1259                 :     973217 : gen_sortspec_uniq(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1260                 :            : {
    1261                 :            :   pari_sp av;
    1262                 :            :   long NX, nx, ny, m, ix, iy, i;
    1263                 :            :   GEN x, y, w, W;
    1264                 :            :   int s;
    1265   [ +  +  +  + ]:     973217 :   switch(n)
    1266                 :            :   {
    1267                 :       1043 :     case 1: return mkvecsmall(1);
    1268                 :            :     case 2:
    1269                 :     323918 :       s = cmp(E,gel(v,1),gel(v,2));
    1270         [ +  + ]:     323918 :       if      (s < 0) return mkvecsmall2(1,2);
    1271         [ +  + ]:     166754 :       else if (s > 0) return mkvecsmall2(2,1);
    1272                 :       2100 :       return mkvecsmall(1);
    1273                 :            :     case 3:
    1274                 :     164080 :       s = cmp(E,gel(v,1),gel(v,2));
    1275         [ +  + ]:     164080 :       if (s < 0) {
    1276                 :      82495 :         s = cmp(E,gel(v,2),gel(v,3));
    1277         [ +  + ]:      82495 :         if (s < 0) return mkvecsmall3(1,2,3);
    1278         [ +  + ]:      53074 :         else if (s == 0) return mkvecsmall2(1,2);
    1279                 :      52682 :         s = cmp(E,gel(v,1),gel(v,3));
    1280         [ +  + ]:      52682 :         if      (s < 0) return mkvecsmall3(1,3,2);
    1281         [ +  + ]:      26712 :         else if (s > 0) return mkvecsmall3(3,1,2);
    1282                 :       1715 :         return mkvecsmall2(1,2);
    1283         [ +  + ]:      81585 :       } else if (s > 0) {
    1284                 :      79982 :         s = cmp(E,gel(v,1),gel(v,3));
    1285         [ +  + ]:      79982 :         if (s < 0) return mkvecsmall3(2,1,3);
    1286         [ +  + ]:      54859 :         else if (s == 0) return mkvecsmall2(2,1);
    1287                 :      53781 :         s = cmp(E,gel(v,2),gel(v,3));
    1288         [ +  + ]:      53781 :         if (s < 0) return mkvecsmall3(2,3,1);
    1289         [ +  + ]:      28105 :         else if (s > 0) return mkvecsmall3(3,2,1);
    1290                 :        630 :         return mkvecsmall2(2,1);
    1291                 :            :       } else {
    1292                 :       1603 :         s = cmp(E,gel(v,1),gel(v,3));
    1293         [ +  + ]:       1603 :         if (s < 0) return mkvecsmall2(1,3);
    1294         [ +  + ]:       1015 :         else if (s == 0) return mkvecsmall(1);
    1295                 :        798 :         return mkvecsmall2(3,1);
    1296                 :            :       }
    1297                 :            :   }
    1298                 :     484176 :   NX = nx = n>>1; ny = n-nx;
    1299                 :     484176 :   av = avma;
    1300                 :     484176 :   x = gen_sortspec_uniq(v,   nx,E,cmp); nx = lg(x)-1;
    1301                 :     484176 :   y = gen_sortspec_uniq(v+NX,ny,E,cmp); ny = lg(y)-1;
    1302                 :     484176 :   w = cgetg(n+1, t_VECSMALL);
    1303                 :     484176 :   m = ix = iy = 1;
    1304 [ +  + ][ +  + ]:    7226415 :   while (ix<=nx && iy<=ny)
    1305                 :            :   {
    1306                 :    6742239 :     s = cmp(E, gel(v,x[ix]), gel(v,y[iy]+NX));
    1307         [ +  + ]:    6742239 :     if (s < 0)
    1308                 :    2874354 :       w[m++] = x[ix++];
    1309         [ +  + ]:    3867885 :     else if (s > 0)
    1310                 :    2990631 :       w[m++] = y[iy++]+NX;
    1311                 :            :     else {
    1312                 :     877254 :       w[m++] = x[ix++];
    1313                 :     877254 :       iy++;
    1314                 :            :     }
    1315                 :            :   }
    1316         [ +  + ]:     867489 :   while (ix<=nx) w[m++] = x[ix++];
    1317         [ +  + ]:     923895 :   while (iy<=ny) w[m++] = y[iy++]+NX;
    1318                 :     484176 :   avma = av;
    1319                 :     484176 :   W = cgetg(m, t_VECSMALL);
    1320         [ +  + ]:    8049447 :   for (i = 1; i < m; i++) W[i] = w[i];
    1321                 :     973217 :   return W;
    1322                 :            : }
    1323                 :            : 
    1324                 :            : /* return permutation sorting v[1..n]. Assume n > 0 */
    1325                 :            : static GEN
    1326                 :  108317962 : gen_sortspec(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
    1327                 :            : {
    1328                 :            :   long nx, ny, m, ix, iy;
    1329                 :            :   GEN x, y, w;
    1330   [ +  +  +  + ]:  108317962 :   switch(n)
    1331                 :            :   {
    1332                 :            :     case 1:
    1333                 :    2099328 :       (void)cmp(E,gel(v,1),gel(v,1)); /* check for type error */
    1334                 :    2099328 :       return mkvecsmall(1);
    1335                 :            :     case 2:
    1336                 :   89633804 :       return cmp(E,gel(v,1),gel(v,2)) <= 0? mkvecsmall2(1,2)
    1337         [ +  + ]:   44816895 :                                           : mkvecsmall2(2,1);
    1338                 :            :     case 3:
    1339         [ +  + ]:   20852550 :       if (cmp(E,gel(v,1),gel(v,2)) <= 0) {
    1340         [ +  + ]:   15029770 :         if (cmp(E,gel(v,2),gel(v,3)) <= 0) return mkvecsmall3(1,2,3);
    1341                 :    8969942 :         return (cmp(E,gel(v,1),gel(v,3)) <= 0)? mkvecsmall3(1,3,2)
    1342         [ +  + ]:    4484971 :                                               : mkvecsmall3(3,1,2);
    1343                 :            :       } else {
    1344         [ +  + ]:    5822780 :         if (cmp(E,gel(v,1),gel(v,3)) <= 0) return mkvecsmall3(2,1,3);
    1345                 :    7244674 :         return (cmp(E,gel(v,2),gel(v,3)) <= 0)? mkvecsmall3(2,3,1)
    1346         [ +  + ]:    3622337 :                                               : mkvecsmall3(3,2,1);
    1347                 :            :       }
    1348                 :            :   }
    1349                 :   40549175 :   nx = n>>1; ny = n-nx;
    1350                 :   40549175 :   w = cgetg(n+1,t_VECSMALL);
    1351                 :   40549175 :   x = gen_sortspec(v,   nx,E,cmp);
    1352                 :   40549161 :   y = gen_sortspec(v+nx,ny,E,cmp);
    1353                 :   40549161 :   m = ix = iy = 1;
    1354 [ +  + ][ +  + ]:  256367107 :   while (ix<=nx && iy<=ny)
    1355         [ +  + ]:  215817946 :     if (cmp(E, gel(v,x[ix]), gel(v,y[iy]+nx))<=0)
    1356                 :  119810223 :       w[m++] = x[ix++];
    1357                 :            :     else
    1358                 :   96007723 :       w[m++] = y[iy++]+nx;
    1359         [ +  + ]:   62592116 :   while (ix<=nx) w[m++] = x[ix++];
    1360         [ +  + ]:  101891318 :   while (iy<=ny) w[m++] = y[iy++]+nx;
    1361                 :  108317934 :   avma = (pari_sp)w; return w;
    1362                 :            : }
    1363                 :            : 
    1364                 :            : static void
    1365                 :   25457593 : init_sort(GEN *x, long *tx, long *lx)
    1366                 :            : {
    1367                 :   25457593 :   *tx = typ(*x);
    1368         [ +  + ]:   25457593 :   if (*tx == t_LIST)
    1369                 :            :   {
    1370         [ -  + ]:         35 :     if (list_typ(*x)!=t_LIST_RAW) pari_err_TYPE("sort",*x);
    1371                 :         35 :     *x = list_data(*x);
    1372         [ +  + ]:         35 :     *lx = *x? lg(*x): 1;
    1373                 :            :   } else {
    1374 [ +  + ][ -  + ]:   25457558 :     if (!is_matvec_t(*tx) && *tx != t_VECSMALL) pari_err_TYPE("gen_sort",*x);
    1375                 :   25457558 :     *lx = lg(*x);
    1376                 :            :   }
    1377                 :   25457593 : }
    1378                 :            : 
    1379                 :            : /* (x o y)[1..lx-1], destroy y */
    1380                 :            : INLINE GEN
    1381                 :    1261910 : sort_extract(GEN x, GEN y, long tx, long lx)
    1382                 :            : {
    1383                 :            :   long i;
    1384      [ +  +  + ]:    1261910 :   switch(tx)
    1385                 :            :   {
    1386                 :            :     case t_VECSMALL:
    1387         [ +  + ]:         56 :       for (i=1; i<lx; i++) y[i] = x[y[i]];
    1388                 :         14 :       break;
    1389                 :            :     case t_LIST:
    1390                 :          7 :       settyp(y,t_VEC);
    1391         [ +  + ]:         35 :       for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1392                 :          7 :       return gtolist(y);
    1393                 :            :     default:
    1394                 :    1261889 :       settyp(y,tx);
    1395         [ +  + ]:    8390750 :       for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,y[i]));
    1396                 :            :   }
    1397                 :    1261910 :   return y;
    1398                 :            : }
    1399                 :            : 
    1400                 :            : /* Sort the vector x, using cmp to compare entries. */
    1401                 :            : GEN
    1402                 :       4907 : gen_sort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1403                 :            : {
    1404                 :            :   long tx, lx;
    1405                 :            :   GEN y;
    1406                 :            : 
    1407                 :       4907 :   init_sort(&x, &tx, &lx);
    1408 [ +  + ][ -  + ]:       4907 :   if (lx==1) return tx == t_LIST? listcreate(): cgetg(1, tx);
    1409                 :       4739 :   y = gen_sortspec_uniq(x,lx-1,E,cmp);
    1410                 :       4907 :   return sort_extract(x, y, tx, lg(y)); /* lg(y) <= lx */
    1411                 :            : }
    1412                 :            : /* Sort the vector x, using cmp to compare entries. */
    1413                 :            : GEN
    1414                 :    1920541 : gen_sort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1415                 :            : {
    1416                 :            :   long tx, lx;
    1417                 :            :   GEN y;
    1418                 :            : 
    1419                 :    1920541 :   init_sort(&x, &tx, &lx);
    1420 [ +  + ][ +  + ]:    1920541 :   if (lx==1) return tx == t_LIST? listcreate(): cgetg(1, tx);
    1421                 :    1257185 :   y = gen_sortspec(x,lx-1,E,cmp);
    1422                 :    1920527 :   return sort_extract(x, y, tx, lx);
    1423                 :            : }
    1424                 :            : /* indirect sort: return the permutation that would sort x */
    1425                 :            : GEN
    1426                 :        126 : gen_indexsort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1427                 :            : {
    1428                 :            :   long tx, lx;
    1429                 :        126 :   init_sort(&x, &tx, &lx);
    1430         [ -  + ]:        126 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1431                 :        126 :   return gen_sortspec_uniq(x,lx-1,E,cmp);
    1432                 :            : }
    1433                 :            : /* indirect sort: return the permutation that would sort x */
    1434                 :            : GEN
    1435                 :    1295415 : gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
    1436                 :            : {
    1437                 :            :   long tx, lx;
    1438                 :    1295415 :   init_sort(&x, &tx, &lx);
    1439         [ +  + ]:    1295415 :   if (lx==1) return cgetg(1, t_VECSMALL);
    1440                 :    1295415 :   return gen_sortspec(x,lx-1,E,cmp);
    1441                 :            : }
    1442                 :            : 
    1443                 :            : /* Sort the vector x in place, using cmp to compare entries */
    1444                 :            : void
    1445                 :   22236604 : gen_sort_inplace(GEN x, void *E, int (*cmp)(void*,GEN,GEN), GEN *perm)
    1446                 :            : {
    1447                 :            :   long tx, lx, i;
    1448                 :   22236604 :   pari_sp av = avma;
    1449                 :            :   GEN y;
    1450                 :            : 
    1451                 :   22236604 :   init_sort(&x, &tx, &lx);
    1452         [ +  + ]:   22236604 :   if (lx<=2)
    1453                 :            :   {
    1454 [ +  + ][ -  + ]:     118335 :     if (perm) *perm = lx == 1? cgetg(1, t_VECSMALL): mkvecsmall(1);
    1455                 :   22236604 :     return;
    1456                 :            :   }
    1457                 :   22118269 :   y = gen_sortspec(x,lx-1, E, cmp);
    1458         [ +  + ]:   22118269 :   if (perm)
    1459                 :            :   {
    1460                 :       9212 :     GEN z = new_chunk(lx);
    1461         [ +  + ]:     129913 :     for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
    1462         [ +  + ]:     129913 :     for (i=1; i<lx; i++) gel(x,i) = gel(z,i);
    1463                 :       9212 :     *perm = y;
    1464                 :       9212 :     avma = (pari_sp)y;
    1465                 :            :   } else {
    1466         [ +  + ]:  125505889 :     for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
    1467         [ +  + ]:  125505889 :     for (i=1; i<lx; i++) gel(x,i) = gel(y,i);
    1468                 :   22236604 :     avma = av;
    1469                 :            :   }
    1470                 :            : }
    1471                 :            : 
    1472                 :            : static int
    1473                 :       7889 : closurecmp(void *data, GEN x, GEN y)
    1474                 :            : {
    1475                 :       7889 :   pari_sp av = avma;
    1476                 :       7889 :   GEN z = closure_callgen2((GEN)data, x,y);
    1477         [ -  + ]:       7889 :   if (typ(z) != t_INT)
    1478                 :          0 :     pari_err_TYPE("closurecmp, cmp. fun. must return an integer", z);
    1479                 :       7889 :   avma = av; return signe(z);
    1480                 :            : }
    1481                 :            : 
    1482                 :            : static void
    1483                 :        126 : check_positive_entries(GEN k)
    1484                 :            : {
    1485                 :        126 :   long i, l = lg(k);
    1486         [ +  + ]:        287 :   for (i=1; i<l; i++)
    1487         [ -  + ]:        161 :     if (k[i] <= 0) pari_err_DOMAIN("sort_function", "index", "<", gen_0, stoi(k[i]));
    1488                 :        126 : }
    1489                 :            : 
    1490                 :            : typedef int (*CMP_FUN)(void*,GEN,GEN);
    1491                 :            : static CMP_FUN
    1492                 :        735 : sort_function(void **E, GEN x, GEN k)
    1493                 :            : {
    1494                 :        735 :   int (*cmp)(GEN,GEN) = &lexcmp;
    1495         [ +  + ]:        735 :   if (!k)
    1496                 :            :   {
    1497         [ +  + ]:        112 :     *E = (void*)((typ(x) == t_VECSMALL)? cmp_small: cmp);
    1498                 :        112 :     return &cmp_nodata;
    1499                 :            :   }
    1500         [ -  + ]:        623 :   if (typ(x) == t_VECSMALL) pari_err_TYPE("sort_function", x);
    1501   [ +  +  -  +  :        623 :   switch(typ(k))
                      - ]
    1502                 :            :   {
    1503                 :         91 :     case t_INT: k = mkvecsmall(itos(k));  break;
    1504                 :         35 :     case t_VEC: case t_COL: k = ZV_to_zv(k); break;
    1505                 :          0 :     case t_VECSMALL: break;
    1506                 :            :     case t_CLOSURE:
    1507 [ +  - ][ -  + ]:        497 :      if (closure_arity(k) != 2 || closure_is_variadic(k))
    1508                 :          0 :        pari_err_TYPE("sort_function, cmp. fun. needs exactly 2 arguments",k);
    1509                 :        497 :      *E = (void*)k;
    1510                 :        497 :      return &closurecmp;
    1511                 :          0 :     default: pari_err_TYPE("sort_function",k);
    1512                 :            :   }
    1513                 :        126 :   check_positive_entries(k);
    1514                 :        126 :   *E = (void*)k;
    1515                 :        735 :   return &veccmp;
    1516                 :            : }
    1517                 :            : 
    1518                 :            : #define cmp_IND 1
    1519                 :            : #define cmp_LEX 2 /* FIXME: backward compatibility, ignored */
    1520                 :            : #define cmp_REV 4
    1521                 :            : #define cmp_UNIQ 8
    1522                 :            : GEN
    1523                 :        665 : vecsort0(GEN x, GEN k, long flag)
    1524                 :            : {
    1525                 :            :   void *E;
    1526                 :        665 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
    1527                 :            : 
    1528 [ +  - ][ -  + ]:        665 :   if (flag < 0 || flag > (cmp_REV|cmp_LEX|cmp_IND|cmp_UNIQ))
    1529                 :          0 :     pari_err_FLAG("vecsort");
    1530         [ +  + ]:        665 :   if (flag & cmp_UNIQ)
    1531         [ +  + ]:         28 :     x = flag & cmp_IND? gen_indexsort_uniq(x, E, CMP): gen_sort_uniq(x, E, CMP);
    1532                 :            :   else
    1533         [ +  + ]:        637 :     x = flag & cmp_IND? gen_indexsort(x, E, CMP): gen_sort(x, E, CMP);
    1534         [ +  + ]:        651 :   if (flag & cmp_REV) { /* reverse order */
    1535                 :            :     long j, lx;
    1536                 :            :     GEN y;
    1537         [ +  + ]:         21 :     if (typ(x)==t_LIST)
    1538                 :            :     {
    1539                 :          7 :       y = list_data(x);
    1540         [ +  - ]:          7 :       if (!y) return x;
    1541                 :            :     }
    1542                 :            :     else
    1543                 :         14 :       y = x;
    1544                 :         14 :     lx = lg(y);
    1545         [ +  + ]:         42 :     for (j=1; j<=(lx-1)>>1; j++) swap(gel(y,j), gel(y,lx-j));
    1546                 :            :   }
    1547                 :        651 :   return x;
    1548                 :            : }
    1549                 :            : 
    1550                 :            : GEN
    1551                 :       5067 : indexsort(GEN x) { return gen_indexsort(x, (void*)&gcmp, cmp_nodata); }
    1552                 :            : GEN
    1553                 :          0 : indexlexsort(GEN x) { return gen_indexsort(x, (void*)&lexcmp, cmp_nodata); }
    1554                 :            : GEN
    1555                 :          0 : indexvecsort(GEN x, GEN k)
    1556                 :            : {
    1557         [ #  # ]:          0 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1558                 :          0 :   return gen_indexsort(x, (void*)k, &veccmp);
    1559                 :            : }
    1560                 :            : 
    1561                 :            : GEN
    1562                 :     589763 : sort(GEN x) { return gen_sort(x, (void*)gcmp, cmp_nodata); }
    1563                 :            : GEN
    1564                 :          7 : lexsort(GEN x) { return gen_sort(x, (void*)lexcmp, cmp_nodata); }
    1565                 :            : GEN
    1566                 :          0 : vecsort(GEN x, GEN k)
    1567                 :            : {
    1568         [ #  # ]:          0 :   if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
    1569                 :          0 :   return gen_sort(x, (void*)k, &veccmp);
    1570                 :            : }
    1571                 :            : long
    1572                 :         70 : vecsearch(GEN v, GEN x, GEN k)
    1573                 :            : {
    1574                 :         70 :   pari_sp av = avma;
    1575                 :            :   void *E;
    1576                 :         70 :   int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
    1577                 :            :   long r;
    1578         [ -  + ]:         70 :   if (!is_matvec_t(typ(v))) pari_err_TYPE("vecsearch", v);
    1579                 :         70 :   r = gen_search(v, x, 0, E, CMP);
    1580                 :         70 :   avma = av; return r;
    1581                 :            : }
    1582                 :            : 
    1583                 :            : GEN
    1584                 :        448 : ZV_indexsort(GEN L) { return gen_indexsort(L, (void*)&cmpii, &cmp_nodata); }
    1585                 :            : GEN
    1586                 :     797671 : ZV_sort(GEN L) { return gen_sort(L, (void*)&cmpii, &cmp_nodata); }
    1587                 :            : GEN
    1588                 :       2177 : ZV_sort_uniq(GEN L) { return gen_sort_uniq(L, (void*)&cmpii, &cmp_nodata); }
    1589                 :            : 
    1590                 :            : /********************************************************************/
    1591                 :            : /**                      SEARCH IN SORTED VECTOR                   **/
    1592                 :            : /********************************************************************/
    1593                 :            : /* index of x in table T, 0 otherwise */
    1594                 :            : long
    1595                 :    7416539 : tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
    1596                 :            : {
    1597                 :    7416539 :   long l = 1, u = lg(T)-1, i, s;
    1598                 :            : 
    1599         [ +  + ]:   40361175 :   while (u>=l)
    1600                 :            :   {
    1601                 :   34944797 :     i = (l+u)>>1; s = cmp(x, gel(T,i));
    1602         [ +  + ]:   34944797 :     if (!s) return i;
    1603         [ +  + ]:   32944636 :     if (s<0) u=i-1; else l=i+1;
    1604                 :            :   }
    1605                 :    7416539 :   return 0;
    1606                 :            : }
    1607                 :            : 
    1608                 :            : /* looks if x belongs to the set T and returns the index if yes, 0 if no */
    1609                 :            : long
    1610                 :   14044674 : gen_search(GEN T, GEN x, long flag, void *data, int (*cmp)(void*,GEN,GEN))
    1611                 :            : {
    1612                 :   14044674 :   long lx = lg(T), i, l, u, s;
    1613                 :            : 
    1614         [ -  + ]:   14044674 :   if (lx==1) return flag? 1: 0;
    1615                 :   14044674 :   l = 1; u = lx-1;
    1616                 :            :   do
    1617                 :            :   {
    1618                 :   47536398 :     i = (l+u)>>1; s = cmp(data, x, gel(T,i));
    1619 [ +  + ][ +  - ]:   47536398 :     if (!s) return flag? 0: i;
    1620         [ +  + ]:   33570936 :     if (s<0) u=i-1; else l=i+1;
    1621         [ +  + ]:   33570936 :   } while (u>=l);
    1622         [ +  + ]:      79212 :   if (!flag) return 0;
    1623         [ +  + ]:   14044674 :   return (s<0)? i: i+1;
    1624                 :            : }
    1625                 :            : 
    1626                 :            : long
    1627                 :     894145 : ZV_search(GEN x, GEN y) { return tablesearch(x, y, cmpii); }
    1628                 :            : 
    1629                 :            : long
    1630                 :    6519986 : zv_search(GEN x, long y) { return tablesearch(x, (GEN)y, cmp_small); }
    1631                 :            : 
    1632                 :            : /********************************************************************/
    1633                 :            : /**                   COMPARISON FUNCTIONS                         **/
    1634                 :            : /********************************************************************/
    1635                 :            : int
    1636                 :  364039745 : cmp_nodata(void *data, GEN x, GEN y)
    1637                 :            : {
    1638                 :  364039745 :   int (*cmp)(GEN,GEN)=(int (*)(GEN,GEN)) data;
    1639                 :  364039745 :   return cmp(x,y);
    1640                 :            : }
    1641                 :            : 
    1642                 :            : /* assume x and y come from the same idealprimedec call (uniformizer unique) */
    1643                 :            : int
    1644                 :     520769 : cmp_prime_over_p(GEN x, GEN y)
    1645                 :            : {
    1646                 :     520769 :   long k = pr_get_f(x) - pr_get_f(y); /* diff. between residue degree */
    1647         [ +  + ]:     543249 :   return k? ((k > 0)? 1: -1)
    1648         [ +  + ]:     543249 :           : ZV_cmp(pr_get_gen(x), pr_get_gen(y));
    1649                 :            : }
    1650                 :            : 
    1651                 :            : int
    1652                 :      55153 : cmp_prime_ideal(GEN x, GEN y)
    1653                 :            : {
    1654                 :      55153 :   int k = cmpii(pr_get_p(x), pr_get_p(y));
    1655         [ +  + ]:      55153 :   return k? k: cmp_prime_over_p(x,y);
    1656                 :            : }
    1657                 :            : 
    1658                 :            : /* assume x and y are t_POL in the same variable whose coeffs can be
    1659                 :            :  * compared (used to sort polynomial factorizations) */
    1660                 :            : int
    1661                 :     821125 : gen_cmp_RgX(void *data, GEN x, GEN y)
    1662                 :            : {
    1663                 :     821125 :   int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
    1664                 :     821125 :   long i, lx = lg(x), ly = lg(y);
    1665                 :            :   int fl;
    1666         [ +  + ]:     821125 :   if (lx > ly) return  1;
    1667         [ +  + ]:     766995 :   if (lx < ly) return -1;
    1668         [ +  + ]:    1853195 :   for (i=lx-1; i>1; i--)
    1669         [ +  + ]:    1720640 :     if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
    1670                 :     821125 :   return 0;
    1671                 :            : }
    1672                 :            : 
    1673                 :            : static int
    1674                 :       2576 : cmp_RgX_Rg(GEN x, GEN y)
    1675                 :            : {
    1676                 :       2576 :   long lx = lgpol(x), ly;
    1677         [ +  + ]:       2576 :   if (lx > 1) return  1;
    1678                 :       2177 :   ly = gequal0(y) ? 0:1;
    1679         [ -  + ]:       2177 :   if (lx > ly) return  1;
    1680         [ -  + ]:       2177 :   if (lx < ly) return -1;
    1681         [ +  + ]:       2177 :   if (lx==0) return 0;
    1682                 :       2576 :   return gcmp(gel(x,2), y);
    1683                 :            : }
    1684                 :            : int
    1685                 :      90403 : cmp_RgX(GEN x, GEN y)
    1686                 :            : {
    1687         [ +  + ]:      90403 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    1688         [ +  + ]:      90403 :   if (typ(y) == t_POLMOD) y = gel(y,2);
    1689         [ +  + ]:      90403 :   if (typ(x) == t_POL) {
    1690         [ +  + ]:      43263 :     if (typ(y) != t_POL) return cmp_RgX_Rg(x, y);
    1691                 :            :   } else {
    1692         [ +  + ]:      47140 :     if (typ(y) != t_POL) return gcmp(x,y);
    1693                 :       1225 :     return - cmp_RgX_Rg(y,x);
    1694                 :            :   }
    1695                 :      90403 :   return gen_cmp_RgX((void*)&gcmp,x,y);
    1696                 :            : }
    1697                 :            : 
    1698                 :            : int
    1699                 :     125193 : cmp_Flx(GEN x, GEN y)
    1700                 :            : {
    1701                 :     125193 :   long i, lx = lg(x), ly = lg(y);
    1702         [ +  + ]:     125193 :   if (lx > ly) return  1;
    1703         [ +  + ]:     114021 :   if (lx < ly) return -1;
    1704         [ +  + ]:     243821 :   for (i=lx-1; i>1; i--)
    1705 [ +  + ][ +  + ]:     242075 :     if (uel(x,i) != uel(y,i)) return uel(x,i)<uel(y,i)? -1: 1;
    1706                 :     125193 :   return 0;
    1707                 :            : }
    1708                 :            : /********************************************************************/
    1709                 :            : /**                   MERGE & SORT FACTORIZATIONS                  **/
    1710                 :            : /********************************************************************/
    1711                 :            : /* merge fx, fy two factorizations, whose 1st column is sorted in strictly
    1712                 :            :  * increasing order wrt cmp */
    1713                 :            : GEN
    1714                 :     690863 : merge_factor(GEN fx, GEN fy, void *data, int (*cmp)(void *,GEN,GEN))
    1715                 :            : {
    1716                 :     690863 :   GEN x = gel(fx,1), e = gel(fx,2), M, E;
    1717                 :     690863 :   GEN y = gel(fy,1), f = gel(fy,2);
    1718                 :     690863 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    1719                 :            : 
    1720                 :     690863 :   M = cgetg(l, t_COL);
    1721                 :     690863 :   E = cgetg(l, t_COL);
    1722                 :            : 
    1723                 :     690863 :   m = ix = iy = 1;
    1724 [ +  + ][ +  + ]:    9928625 :   while (ix<lx && iy<ly)
    1725                 :            :   {
    1726                 :    9237762 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    1727         [ +  + ]:    9237762 :     if (s < 0)
    1728                 :    8583302 :     { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; }
    1729         [ +  + ]:     654460 :     else if (s == 0)
    1730                 :            :     {
    1731                 :     103768 :       GEN z = gel(x,ix), g = addii(gel(e,ix), gel(f,iy));
    1732         [ +  + ]:     103768 :       iy++; ix++; if (!signe(g)) continue;
    1733                 :      26740 :       gel(M,m) = z; gel(E,m) = g;
    1734                 :            :     }
    1735                 :            :     else
    1736                 :     550692 :     { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; }
    1737                 :    9160734 :     m++;
    1738                 :            :   }
    1739         [ +  + ]:    4813717 :   while (ix<lx) { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; m++; }
    1740         [ +  + ]:     910643 :   while (iy<ly) { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; m++; }
    1741                 :     690863 :   setlg(M, m);
    1742                 :     690863 :   setlg(E, m); return mkmat2(M, E);
    1743                 :            : }
    1744                 :            : /* merge two sorted vectors, removing duplicates. Shallow */
    1745                 :            : GEN
    1746                 :      32431 : merge_sort_uniq(GEN x, GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    1747                 :            : {
    1748                 :      32431 :   long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
    1749                 :            :   GEN M;
    1750                 :            : 
    1751                 :      32431 :   M = cgetg(l, t_COL);
    1752                 :      32431 :   m = ix = iy = 1;
    1753 [ +  + ][ +  + ]:      53284 :   while (ix<lx && iy<ly)
    1754                 :            :   {
    1755                 :      20853 :     int s = cmp(data, gel(x,ix), gel(y,iy));
    1756         [ +  + ]:      20853 :     if (s < 0)
    1757                 :      13321 :     { gel(M,m) = gel(x,ix); ix++; }
    1758         [ -  + ]:       7532 :     else if (s == 0)
    1759                 :          0 :     { gel(M,m) = gel(x,ix); iy++; ix++; }
    1760                 :            :     else
    1761                 :       7532 :     { gel(M,m) = gel(y,iy); iy++; }
    1762                 :      20853 :     m++;
    1763                 :            :   }
    1764         [ +  + ]:      58940 :   while (ix<lx) { gel(M,m) = gel(x,ix); ix++; m++; }
    1765         [ +  + ]:      52269 :   while (iy<ly) { gel(M,m) = gel(y,iy); iy++; m++; }
    1766                 :      32431 :   setlg(M, m); return M;
    1767                 :            : }
    1768                 :            : 
    1769                 :            : /* sort generic factorization, in place */
    1770                 :            : GEN
    1771                 :    2634659 : sort_factor(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
    1772                 :            : {
    1773                 :            :   GEN a, b, A, B, w;
    1774                 :            :   pari_sp av;
    1775                 :            :   long n, i;
    1776                 :            : 
    1777         [ +  + ]:    2634659 :   a = gel(y,1); n = lg(a); if (n == 1) return y;
    1778                 :    2632447 :   b = gel(y,2); av = avma;
    1779                 :    2632447 :   A = new_chunk(n);
    1780                 :    2632447 :   B = new_chunk(n);
    1781                 :    2632447 :   w = gen_sortspec(a, n-1, data, cmp);
    1782         [ +  + ]:    6717527 :   for (i=1; i<n; i++) { long k = w[i]; A[i] = a[k]; B[i] = b[k]; }
    1783         [ +  + ]:    6717527 :   for (i=1; i<n; i++) { a[i] = A[i]; b[i] = B[i]; }
    1784                 :    2634659 :   avma = av; return y;
    1785                 :            : }
    1786                 :            : /* sort polynomial factorization, in place */
    1787                 :            : GEN
    1788                 :     316757 : sort_factor_pol(GEN y,int (*cmp)(GEN,GEN))
    1789                 :            : {
    1790                 :     316757 :   (void)sort_factor(y,(void*)cmp, &gen_cmp_RgX);
    1791                 :     316757 :   return y;
    1792                 :            : }
    1793                 :            : 
    1794                 :            : /* assume f and g coprime integer factorizations */
    1795                 :            : GEN
    1796                 :        350 : merge_factor_i(GEN f, GEN g)
    1797                 :            : {
    1798         [ -  + ]:        350 :   if (lg(f) == 1) return g;
    1799         [ -  + ]:        350 :   if (lg(g) == 1) return f;
    1800                 :        350 :   return sort_factor(famat_mul_shallow(f,g), (void*)&cmpii, &cmp_nodata);
    1801                 :            : }
    1802                 :            : 
    1803                 :            : /***********************************************************************/
    1804                 :            : /*                                                                     */
    1805                 :            : /*                          SET OPERATIONS                             */
    1806                 :            : /*                                                                     */
    1807                 :            : /***********************************************************************/
    1808                 :            : GEN
    1809                 :        896 : gtoset(GEN x)
    1810                 :            : {
    1811                 :            :   long lx;
    1812         [ -  + ]:        896 :   if (!x) return cgetg(1, t_VEC);
    1813   [ +  +  +  + ]:        896 :   switch(typ(x))
    1814                 :            :   {
    1815                 :            :     case t_VEC:
    1816                 :        868 :     case t_COL: lx = lg(x); break;
    1817                 :            :     case t_LIST:
    1818         [ -  + ]:         14 :       if (list_typ(x)==t_LIST_MAP) return mapdomain(x);
    1819         [ +  + ]:         14 :       x = list_data(x); lx = x? lg(x): 1; break;
    1820                 :          7 :     case t_VECSMALL: lx = lg(x); x = zv_to_ZV(x); break;
    1821                 :          7 :     default: return mkveccopy(x);
    1822                 :            :   }
    1823         [ +  + ]:        889 :   if (lx==1) return cgetg(1,t_VEC);
    1824                 :        882 :   x = gen_sort_uniq(x, (void*)&cmp_universal, cmp_nodata);
    1825                 :        882 :   settyp(x, t_VEC); /* it may be t_COL */
    1826                 :        896 :   return x;
    1827                 :            : }
    1828                 :            : 
    1829                 :            : long
    1830                 :         14 : setisset(GEN x)
    1831                 :            : {
    1832                 :         14 :   long i, lx = lg(x);
    1833                 :            : 
    1834         [ -  + ]:         14 :   if (typ(x) != t_VEC) return 0;
    1835         [ -  + ]:         14 :   if (lx == 1) return 1;
    1836         [ +  + ]:         70 :   for (i=1; i<lx-1; i++)
    1837         [ +  + ]:         63 :     if (cmp_universal(gel(x,i+1), gel(x,i)) <= 0) return 0;
    1838                 :         14 :   return 1;
    1839                 :            : }
    1840                 :            : 
    1841                 :            : long
    1842                 :     105903 : setsearch(GEN T, GEN y, long flag)
    1843                 :            : {
    1844                 :            :   long lx;
    1845      [ +  +  + ]:     105903 :   switch(typ(T))
    1846                 :            :   {
    1847                 :     105889 :     case t_VEC: lx = lg(T); break;
    1848                 :            :     case t_LIST:
    1849         [ -  + ]:          7 :     if (list_typ(T) != t_LIST_RAW) pari_err_TYPE("setsearch",T);
    1850         [ +  - ]:          7 :     T = list_data(T); lx = T? lg(T): 1; break;
    1851                 :          7 :     default: pari_err_TYPE("setsearch",T);
    1852                 :          0 :       return 0; /*not reached*/
    1853                 :            :   }
    1854         [ -  + ]:     105896 :   if (lx==1) return flag? 1: 0;
    1855                 :     105896 :   return gen_search(T,y,flag,(void*)cmp_universal,cmp_nodata);
    1856                 :            : }
    1857                 :            : 
    1858                 :            : GEN
    1859                 :         21 : setunion(GEN x, GEN y)
    1860                 :            : {
    1861                 :         21 :   pari_sp av = avma;
    1862                 :         21 :   long i, j, k, lx = lg(x), ly = lg(y);
    1863                 :         21 :   GEN z = cgetg(lx + ly - 1, t_VEC);
    1864         [ -  + ]:         21 :   if (typ(x) != t_VEC) pari_err_TYPE("setunion",x);
    1865         [ -  + ]:         21 :   if (typ(y) != t_VEC) pari_err_TYPE("setunion",y);
    1866                 :         21 :   i = j = k = 1;
    1867 [ +  + ][ +  - ]:         98 :   while (i<lx && j<ly)
    1868                 :            :   {
    1869                 :         77 :     int s = cmp_universal(gel(x,i), gel(y,j));
    1870         [ +  + ]:         77 :     if (s < 0)
    1871                 :         28 :       z[k++] = x[i++];
    1872         [ +  + ]:         49 :     else if (s > 0)
    1873                 :         14 :       z[k++] = y[j++];
    1874                 :            :     else {
    1875                 :         35 :       z[k++] = x[i++];
    1876                 :         35 :       j++;
    1877                 :            :     }
    1878                 :            :   }
    1879         [ -  + ]:         21 :   while (i<lx) z[k++] = x[i++];
    1880         [ -  + ]:         21 :   while (j<ly) z[k++] = y[j++];
    1881                 :         21 :   setlg(z, k);
    1882                 :         21 :   return gerepilecopy(av, z);
    1883                 :            : }
    1884                 :            : /* in case of equal keys in x,y, take the key from x */
    1885                 :            : GEN
    1886                 :      74975 : ZV_union_shallow(GEN x, GEN y)
    1887                 :            : {
    1888                 :      74975 :   long i, j, k, lx = lg(x), ly = lg(y);
    1889                 :      74975 :   GEN z = cgetg(lx + ly - 1, t_VEC);
    1890                 :      74975 :   i = j = k = 1;
    1891 [ +  + ][ +  + ]:     212091 :   while (i<lx && j<ly)
    1892                 :            :   {
    1893                 :     137116 :     int s = cmpii(gel(x,i), gel(y,j));
    1894         [ +  + ]:     137116 :     if (s < 0)
    1895                 :       2139 :       z[k++] = x[i++];
    1896         [ +  + ]:     134977 :     else if (s > 0)
    1897                 :      50185 :       z[k++] = y[j++];
    1898                 :            :     else {
    1899                 :      84792 :       z[k++] = x[i++];
    1900                 :      84792 :       j++;
    1901                 :            :     }
    1902                 :            :   }
    1903         [ +  + ]:      79962 :   while (i<lx) z[k++] = x[i++];
    1904         [ +  + ]:      80513 :   while (j<ly) z[k++] = y[j++];
    1905                 :      74975 :   setlg(z, k); return z;
    1906                 :            : }
    1907                 :            : 
    1908                 :            : 
    1909                 :            : GEN
    1910                 :          7 : setintersect(GEN x, GEN y)
    1911                 :            : {
    1912                 :          7 :   long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
    1913                 :          7 :   pari_sp av = avma;
    1914                 :          7 :   GEN z = cgetg(lx,t_VEC);
    1915         [ -  + ]:          7 :   if (typ(x) != t_VEC) pari_err_TYPE("setintersect",x);
    1916         [ -  + ]:          7 :   if (typ(y) != t_VEC) pari_err_TYPE("setintersect",y);
    1917 [ +  + ][ +  - ]:         70 :   while (ix < lx && iy < ly)
    1918                 :            :   {
    1919                 :         63 :     int c = cmp_universal(gel(x,ix), gel(y,iy));
    1920         [ +  + ]:         63 :     if      (c < 0) ix++;
    1921         [ +  + ]:         35 :     else if (c > 0) iy++;
    1922                 :         21 :     else { gel(z, iz++) = gel(x,ix); ix++; iy++; }
    1923                 :            :   }
    1924                 :          7 :   setlg(z,iz); return gerepilecopy(av,z);
    1925                 :            : }
    1926                 :            : 
    1927                 :            : GEN
    1928                 :          7 : gen_setminus(GEN A, GEN B, int (*cmp)(GEN,GEN))
    1929                 :            : {
    1930                 :          7 :   pari_sp ltop = avma;
    1931                 :          7 :   long i = 1, j = 1, k = 1, lx = lg(A), ly = lg(B);
    1932                 :          7 :   GEN  diff = cgetg(lx,t_VEC);
    1933 [ +  + ][ +  - ]:         84 :   while (i < lx && j < ly)
    1934   [ +  +  +  - ]:         77 :     switch ( cmp(gel(A,i),gel(B,j)) )
    1935                 :            :     {
    1936                 :         28 :       case -1: gel(diff,k++) = gel(A,i++); break;
    1937                 :         28 :       case 1: j++; break;
    1938                 :         21 :       case 0: i++; break;
    1939                 :            :     }
    1940         [ -  + ]:          7 :   while (i < lx) gel(diff,k++) = gel(A,i++);
    1941                 :          7 :   setlg(diff,k);
    1942                 :          7 :   return gerepilecopy(ltop,diff);
    1943                 :            : }
    1944                 :            : 
    1945                 :            : GEN
    1946                 :          7 : setminus(GEN x, GEN y)
    1947                 :            : {
    1948         [ -  + ]:          7 :   if (typ(x) != t_VEC) pari_err_TYPE("setminus",x);
    1949         [ -  + ]:          7 :   if (typ(y) != t_VEC) pari_err_TYPE("setminus",y);
    1950                 :          7 :   return gen_setminus(x,y,cmp_universal);
    1951                 :            : }
    1952                 :            : 
    1953                 :            : GEN
    1954                 :         21 : setbinop(GEN f, GEN x, GEN y)
    1955                 :            : {
    1956                 :         21 :   pari_sp av = avma;
    1957                 :         21 :   long i, j, lx, ly, k = 1;
    1958                 :            :   GEN z;
    1959 [ +  - ][ +  + ]:         21 :   if (typ(f) != t_CLOSURE || closure_arity(f) != 2 || closure_is_variadic(f))
                 [ -  + ]
    1960                 :          7 :     pari_err_TYPE("setbinop [function needs exactly 2 arguments]",f);
    1961                 :         14 :   lx = lg(x);
    1962         [ -  + ]:         14 :   if (typ(x) != t_VEC) pari_err_TYPE("setbinop", x);
    1963         [ +  + ]:         14 :   if (y == NULL) { /* assume x = y and f symmetric */
    1964                 :          7 :     z = cgetg((((lx-1)*lx) >> 1) + 1, t_VEC);
    1965         [ +  + ]:         28 :     for (i = 1; i < lx; i++)
    1966         [ +  + ]:         63 :       for (j = i; j < lx; j++)
    1967                 :         42 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(x,j));
    1968                 :            :   } else {
    1969                 :          7 :     ly = lg(y);
    1970         [ -  + ]:          7 :     if (typ(y) != t_VEC) pari_err_TYPE("setbinop", y);
    1971                 :          7 :     z = cgetg((lx-1)*(ly-1) + 1, t_VEC);
    1972         [ +  + ]:         28 :     for (i = 1; i < lx; i++)
    1973         [ +  + ]:         84 :       for (j = 1; j < ly; j++)
    1974                 :         63 :         gel(z, k++) = closure_callgen2(f, gel(x,i),gel(y,j));
    1975                 :            :   }
    1976                 :         14 :   return gerepileupto(av, gtoset(z));
    1977                 :            : }

Generated by: LCOV version 1.9