Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - zetamult.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 544 563 96.6 %
Date: 2020-09-18 06:10:04 Functions: 40 40 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : This file is part of the PARI/GP package.
       3             : 
       4             : PARI/GP is free software; you can redistribute it and/or modify it under the
       5             : terms of the GNU General Public License as published by the Free Software
       6             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       7             : ANY WARRANTY WHATSOEVER.
       8             : 
       9             : Check the License for details. You should have received a copy of it, along
      10             : with the package; see the file 'COPYING'. If not, write to the Free Software
      11             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      12             : 
      13             : /********************************************************************/
      14             : /**                                                                **/
      15             : /**                         MULTIPLE ZETA VALUES                   **/
      16             : /**                                                                **/
      17             : /********************************************************************/
      18             : #include "pari.h"
      19             : #include "paripriv.h"
      20             : 
      21             : /********************************************************************/
      22             : /**                           CONVERSIONS                          **/
      23             : /********************************************************************/
      24             : /* vecsmall to binary */
      25             : static long
      26         315 : fd(GEN evec, long a, long b)
      27             : {
      28         315 :   long i, s = 0;
      29        1050 :   for (i = a; i <= b; i++) s = evec[i] | (s << 1);
      30         315 :   return s;
      31             : }
      32             : /* 2^(b-a+1) + fd(evec) */
      33             : static long
      34        7182 : fd1(GEN evec, long a, long b)
      35             : {
      36        7182 :   long i, s = 1;
      37       64694 :   for (i = a; i <= b; i++) s = evec[i] | (s << 1);
      38        7182 :   return s;
      39             : }
      40             : 
      41             : /* m > 0 */
      42             : static GEN
      43        7196 : mtoevec(GEN m)
      44             : {
      45        7196 :   GEN e = vecsmall_append(binary_zv(m), 1);
      46        7196 :   e[1] = 0; return e;
      47             : }
      48             : static GEN
      49        7182 : etoindex(GEN evec) { return utoipos(fd1(evec, 2, lg(evec)-2)); }
      50             : 
      51             : /* dual of evec[1..l-1] */
      52             : static GEN
      53          28 : revslice(GEN evec, long l)
      54             : {
      55          28 :   GEN res = cgetg(l, t_VECSMALL);
      56             :   long i;
      57         168 :   for (i = 1; i < l; i++) res[i] = 1 - evec[l-i];
      58          28 :   return res;
      59             : }
      60             : 
      61             : /* N.B. evec[ne] = 1 */
      62             : static GEN
      63        7861 : etoa(GEN evec)
      64             : {
      65        7861 :   long c = 0, cold = 0, i = 1, l = lg(evec);
      66        7861 :   GEN avec = cgetg(l, t_VECSMALL);
      67       82810 :   while (++c < l)
      68       74949 :     if (evec[c] == 1) { avec[i++] = c - cold; cold = c; }
      69        7861 :   setlg(avec, i); return avec;
      70             : }
      71             : 
      72             : static GEN
      73        7266 : atoe(GEN avec)
      74             : {
      75        7266 :   long i, j, l = lg(avec), n = zv_sum(avec);
      76        7266 :   GEN evec = zero_zv(n);
      77       43603 :   for (i = 1, j = 0; i < l; i++) { long a = avec[i]; j += a; evec[j] = 1; }
      78        7266 :   return evec;
      79             : }
      80             : 
      81             : 
      82             : /* Conversions: types are evec, avec, m (if evec=0y1, m=(1y)_2).
      83             :    fl is respectively 0, 1, 2. Type of a is autodetected. */
      84             : static GEN
      85       14651 : zetamultconvert_i(GEN a, long fl)
      86             : {
      87             :   long i, l;
      88       14651 :   if (fl < 0 || fl > 2) pari_err_FLAG("zetamultconvert");
      89       14651 :   switch(typ(a))
      90             :   {
      91        7203 :     case t_INT:
      92        7203 :       if (signe(a) <= 0) pari_err_TYPE("zetamultconvert",a);
      93             :       switch (fl)
      94             :       {
      95          14 :         case 0: a = mtoevec(a); break;
      96        7182 :         case 1: a = etoa(mtoevec(a)); break;
      97           7 :         case 2: a = icopy(a); break;
      98             :       }
      99        7203 :       break;
     100        7441 :     case t_VEC: case t_COL: case t_VECSMALL:
     101        7441 :       a = gtovecsmall(a);
     102        7441 :       l = lg(a);
     103        7441 :       if (a[1] == 0)
     104             :       {
     105          49 :         if (!a[l-1]) pari_err_TYPE("zetamultconvert", a);
     106         280 :         for (i = 1; i < l; i++)
     107         252 :           if (a[i] & ~1UL) pari_err_TYPE("zetamultconvert", a);
     108             :         switch (fl)
     109             :         {
     110           7 :           case 1: a = etoa(a); break;
     111           7 :           case 2: a = etoindex(a);
     112             :         }
     113          28 :       }
     114             :       else
     115             :       {
     116        7392 :         if (a[1] < 2) pari_err_TYPE("zetamultconvert", a);
     117       37086 :         for (i = 2; i < l; i++)
     118       29715 :           if (a[i] <= 0) pari_err_TYPE("zetamultconvert", a);
     119             :         switch (fl)
     120             :         {
     121          21 :           case 0: a = atoe(a); break;
     122        7175 :           case 2: a = etoindex(atoe(a));
     123             :         }
     124        7399 :       }
     125        7399 :       break;
     126           7 :     default: pari_err_TYPE("zetamultconvert", a);
     127             :   }
     128       14602 :   return a;
     129             : }
     130             : GEN
     131       14399 : zetamultconvert(GEN a, long fl)
     132             : {
     133       14399 :   pari_sp av = avma;
     134       14399 :   return gerepileuptoleaf(av, zetamultconvert_i(a, fl));
     135             : }
     136             : 
     137             : GEN
     138          28 : zetamultdual(GEN s)
     139             : {
     140          28 :   pari_sp av = avma;
     141          28 :   GEN e = zetamultconvert_i(s, 0);
     142          28 :   return gerepileuptoleaf(av, etoa(revslice(e, lg(e))));
     143             : }
     144             : 
     145             : /********************************************************************/
     146             : /**                  AVEC -> LIST OF STAR AVECS                    **/
     147             : /********************************************************************/
     148             : /* all star avecs corresponding to given t_VECSMALL avec */
     149             : static GEN
     150         644 : allstar(GEN avec)
     151             : {
     152         644 :   long i, la = lg(avec), K = 1 << (la - 2);
     153         644 :   GEN w = cgetg(K + 1, t_VEC);
     154             : 
     155         644 :   gel(w, 1) = avec;
     156        1456 :   for (i = 2; i < la; i++)
     157             :   {
     158         812 :     long L = 1 << (i - 2), j;
     159        2044 :     for (j = 1; j <= L; j++)
     160             :     {
     161        1232 :       GEN u = gel(w,j), v;
     162        1232 :       long k, l = lg(u) - 1, ind = l - la + i;
     163        1232 :       gel(w, L + j) = v = cgetg(l, t_VECSMALL);
     164        1708 :       for (k = 1; k < ind; k++) v[k] = u[k];
     165        1232 :       v[ind] = u[ind] + u[ind + 1];
     166        1652 :       for (k = ind + 1; k < l; k++) v[k] = u[k+1];
     167             :     }
     168             :   }
     169         644 :   return w;
     170             : }
     171             : /* same for multipolylogs */
     172             : static GEN
     173           7 : allstar2(GEN avec, GEN zvec)
     174             : {
     175           7 :   long la = lg(avec), i, K = 1 << (la - 2);
     176           7 :   GEN W = cgetg(K + 1, t_VEC), Z = cgetg(K + 1, t_VEC);
     177             : 
     178           7 :   gel(W, 1) = avec;
     179           7 :   gel(Z, 1) = zvec = gtovec(zvec);
     180          35 :   for (i = 2; i < la; i++)
     181             :   {
     182          28 :     long L = 1 << (i - 2), j;
     183         133 :     for (j = 1; j <= L; j++)
     184             :     {
     185         105 :       GEN u = gel(W, j), w, y = gel(Z, j), z;
     186         105 :       long l = lg(u) - 1, ind = l - la + i, k;
     187         105 :       w = cgetg(l, t_VECSMALL);
     188         105 :       z = cgetg(l, t_VEC);
     189         224 :       for (k = 1; k < ind; k++) { w[k] = u[k]; gel(z, k) = gel(y, k); }
     190         105 :       w[ind] = u[ind] + u[ind + 1];
     191         105 :       gel(z, ind) = gmul(gel(y, ind), gel(y, ind + 1));
     192         182 :       for (k = ind + 1; k < l; k++) { w[k] = u[k+1]; gel(z, k) = gel(y, k+1); }
     193         105 :       gel(W, L + j) = w;
     194         105 :       gel(Z, L + j) = z;
     195             :     }
     196             :   }
     197           7 :   return mkvec2(W, Z);
     198             : }
     199             : 
     200             : /**************************************************************/
     201             : /*              ZAGIER & RADCHENKO'S ALGORITHM                */
     202             : /**************************************************************/
     203             : /* accuracy 2^(-b); s << (b/log b)^2, l << b/sqrt(log b) */
     204             : static void
     205         154 : zparams(long *s, long *l, long b)
     206             : {
     207         154 :   double D = b * LOG10_2, E = 3*D/2 / log(3*D);
     208         154 :   *s = (long)floor(E*E);
     209         154 :   *l = (long)floor(sqrt(*s * log((double)*s)/2.));
     210         154 : }
     211             : 
     212             : static GEN
     213         140 : zetamult_Zagier(GEN avec, long bit, long prec)
     214             : {
     215             :   pari_sp av;
     216         140 :   GEN ze, z = NULL, b;
     217         140 :   long h, r, n, s, l, t = lg(avec) - 1;
     218             : 
     219         140 :   zparams(&s, &l, bit);
     220         140 :   ze= cgetg(s + 1, t_VEC);
     221         140 :   b = cgetg(l + 1, t_VEC);
     222       15512 :   for (r = 1; r <= s; r++) gel(ze,r) = cgetr(prec);
     223        2072 :   for (r = 1; r <= l; r++) { gel(b,r) = cgetr(prec); affur(0,gel(b,r)); }
     224         140 :   affur(1, gel(b,1)); av = avma;
     225        1155 :   for (r = 1, h = -1; r <= t; r++)
     226             :   {
     227        1015 :     long m, j = avec[r];
     228             :     GEN q;
     229             : 
     230        1015 :     h += j - 1; z = gen_0;
     231        1015 :     q = h? invr(itor(powuu(s,h), prec)): real_1(prec);
     232       18900 :     for (m = 1; m <= l; m++)
     233             :     {
     234             :       pari_sp av2;
     235       17885 :       GEN S = gel(b, m), C;
     236       17885 :       q = divru(q, s); av2 = avma;
     237       17885 :       C = binomialuu(h + m, m);
     238      224616 :       for (n = 1; n < m; n++)
     239             :       { /* C = binom(h+m, m-n+1) */
     240      206731 :         S = gsub(S, mulir(C, gel(b, n)));
     241      206731 :         C = diviuexact(muliu(C, m-n+1), h+n);
     242             :       }
     243       17885 :       affrr(divru(S, h+m), gel(b,m)); set_avma(av2);
     244       17885 :       z = gadd(z, gmul(gel(b, m), q));
     245             :     }
     246      162239 :     for (m = s; m >= 1; m--)
     247             :     {
     248      161224 :       GEN z1 = r == 1? ginv(powuu(m,j)): gdiv(gel(ze, m), powuu(m,j));
     249      161224 :       z1 = gadd(z, z1);
     250      161224 :       affrr(z, gel(ze, m)); z = z1;
     251             :     }
     252        1015 :     set_avma(av);
     253             :   }
     254         140 :   return z;
     255             : }
     256             : 
     257             : /* Compute t-mzvs; slower than Zagier's code for t=0 */
     258             : static GEN
     259          14 : zetamult_interpolate2_i(GEN avec, GEN t, long prec)
     260             : {
     261             :   pari_sp av;
     262             :   GEN a, b, ze, _1;
     263             :   long i, j, n, s, l;
     264             : 
     265          14 :   zparams(&s, &l, prec2nbits(prec));
     266          14 :   n = lg(avec) - 1;
     267          14 :   a = zeromatcopy(n + 1, l);
     268          14 :   b = zerovec(l + 1);
     269          84 :   for (i = 1; i <= n; i++)
     270             :   {
     271          70 :     long h = avec[n + 1 - i];
     272          70 :     if (i == 1) gel(b, h) = gen_1;
     273             :     else
     274        1176 :       for (j = l + 1; j >= 1; j--)
     275             :       {
     276        1120 :         if (j <= h) gel(b, j) = gen_0;
     277        1036 :         else gel(b, j) = gadd(gcoeff(a, i, j-h), gmul(t, gel(b, j-h)));
     278             :       }
     279          70 :     gcoeff(a, i+1, 1) = gel(b, 2); /* j = 1 */
     280        1330 :     for (j = 2; j <= l; j++)
     281             :     { /* b[j+1] - sum_{0 <= u < j-1} binom(j,u) a[i+1,u+1]*/
     282        1260 :       pari_sp av = avma;
     283        1260 :       GEN S = gel(b, j + 1);
     284        1260 :       S = gsub(S, gcoeff(a, i+1, 1)); /* u = 0 */
     285        1260 :       if (j > 2) S = gsub(S, gmulgs(gcoeff(a, i+1, 2), j)); /* u = 1 */
     286        1260 :       if (j >= 4)
     287             :       {
     288        1120 :         GEN C = utoipos(j*(j-1) / 2);
     289        1120 :         long u, U = (j-1)/2;
     290        5600 :         for (u = 2; u <= U; u++)
     291             :         { /* C = binom(j, u) = binom(j, j-u) */
     292        4480 :           GEN A = gadd(gcoeff(a, i+1, u+1), gcoeff(a, i+1, j-u+1));
     293        4480 :           S = gsub(S, gmul(C, A));
     294        4480 :           C = diviuexact(muliu(C, j-u), u+1);
     295             :         }
     296        1120 :         if (!odd(j)) S = gsub(S, gmul(C, gcoeff(a,i+1, j/2+1)));
     297             :       }
     298        1260 :       gcoeff(a, i+1, j) = gerepileupto(av, gdivgs(S, j));
     299             :     }
     300             :   }
     301          14 :   _1 = real_1(prec + EXTRAPRECWORD); av = avma;
     302          14 :   ze = cgetg(n+1, t_VEC);
     303          84 :   for (i = 1; i <= n; i++)
     304             :   {
     305          70 :     GEN S = gdivgs(gcoeff(a, n+2-i, 1), s), sj = utoipos(s);
     306        1330 :     for (j = 2; j <= l; j++)
     307             :     {
     308        1260 :       sj = muliu(sj, s); /* = s^j */
     309        1260 :       S = gadd(S, gdiv(gcoeff(a, n+2-i, j), sj));
     310             :     }
     311          70 :     gel(ze, i) = S;
     312             :   }
     313        2086 :   for (i = s; i >= 1; i--)
     314             :   {
     315        2072 :     GEN b0 = divri(_1, powuu(i, avec[n])), z;
     316        2072 :     z = gel(ze,n); gel(ze,n) = gadd(z, b0);
     317       10360 :     for (j = n-1; j >= 1; j--)
     318             :     {
     319        8288 :       b0 = gdiv(gadd(gmul(t, b0), z), powuu(i, avec[j]));
     320        8288 :       z = gel(ze,j); gel(ze,j) = gadd(gel(ze,j), b0);
     321             :     }
     322        2072 :     if (gc_needed(av, 1))
     323             :     {
     324           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"zetamult: i = %ld", i);
     325           0 :       ze = gerepilecopy(av, ze);
     326             :     }
     327             :   }
     328          14 :   return gel(ze, 1);
     329             : }
     330             : 
     331             : /********************************************************************/
     332             : /**                      AKHILESH ALGORITHM                        **/
     333             : /********************************************************************/
     334             : /* a t_VECSMALL, upper bound for -log2(zeta(a)) */
     335             : static long
     336          70 : log2zeta_bound(GEN a)
     337          70 : { return ceil(-dbllog2(zetamult_Zagier(a, 32, LOWDEFAULTPREC))); }
     338             : /* ibin[n+1] = 1 / binom(2n, n) as a t_REAL */
     339             : static void
     340         329 : get_ibin(GEN *pibin, GEN *pibin1, long N, long prec)
     341             : {
     342             :   GEN ibin, ibin1;
     343             :   long n;
     344         329 :   *pibin = ibin = cgetg(N + 2, t_VEC);
     345         329 :   *pibin1= ibin1= cgetg(N + 2, t_VEC);
     346         329 :   gel(ibin,1) = gel(ibin1,1) = gen_0; /* unused */
     347         329 :   gel(ibin,2) = gel(ibin1,2) = real2n(-1,prec);
     348       39851 :   for (n = 2; n <= N; n++)
     349             :   {
     350       39522 :     gel(ibin, n+1) = divru(mulru(gel(ibin, n), n), 4*n-2);
     351       39522 :     gel(ibin1, n+1) = divru(gel(ibin, n+1), n);
     352             :   }
     353         329 : }
     354             : /**************************************************************/
     355             : /*                         ALL MZV's                          */
     356             : /**************************************************************/
     357             : 
     358             : /* Generalization to Multiple Polylogarithms.
     359             : The basic function is polylogmult which takes two arguments
     360             : avec, as above, and zvec, the vector of z_j, and the result
     361             : is $\sum_{n_1>n_2>...>n_r>0}z_1^{n_1}...z_r^{n_r}/(n_1^a_1...n_r^{a_r})$. */
     362             : 
     363             : /* Given admissible evec = xe_2....e_{k-1}y, (k>=2), compute a,b,v such that
     364             : evec = x{1}_{a-1}v{0}_{b-1}y with v empty or admissible.
     365             : Input: vector w=evec
     366             : Output: v=wmid, wini, wfin */
     367             : static void
     368        3108 : findabvgen(GEN evec, long _0, long _1, GEN *pwmid, GEN *pwini, GEN *pwfin,
     369             :            long *pa, long *pb)
     370             : {
     371        3108 :   long s = lg(evec) - 1, m, a, b, j, x = evec[1], y = evec[s];
     372             :   GEN wmid, wini, wfin;
     373        3108 :   if (s == 2)
     374             :   {
     375         217 :     *pwmid = cgetg(1, t_VECSMALL);
     376         217 :     *pwini = mkvecsmall(x);
     377         217 :     *pwfin = mkvecsmall(y);
     378         217 :     *pa = *pb = 1; return;
     379             :   }
     380        2891 :   a = s - 1;
     381        2891 :   for (j = 1; j <= s - 2; j++) if (evec[j + 1] != _1) { a = j; break; }
     382        2891 :   *pa = a;
     383        2891 :   b = s - 1;
     384       12467 :   for (j = s - 2; j >= 1; j--) if (evec[j + 1] != _0) { b = s - 1 - j; break; }
     385        2891 :   *pb = b;
     386             : 
     387        2891 :   *pwmid = wmid = a+b < s? vecslice(evec, a+1, s-b): cgetg(1, t_VECSMALL);
     388        2891 :   m = lg(wmid) - 1;
     389        2891 :   *pwini = wini = cgetg(a + m + 1, t_VECSMALL);
     390        2891 :   wini[1] = x; for (j = 2; j <= a; j++) wini[j] = _1;
     391       11669 :   for (; j <= a + m; j++) wini[j] = wmid[j-a];
     392        2891 :   *pwfin = wfin = cgetg(b + m + 1, t_VECSMALL);
     393       11669 :   for (j = 1; j <= m; j++) wfin[j] = wmid[j];
     394       12467 :   for (; j < b + m; j++) wfin[j] = _0;
     395        2891 :   wfin[j] = y;
     396             : }
     397             : 
     398             : /* y != 0,1 */
     399             : static GEN
     400         196 : filllg1(GEN ibin1, GEN r1, GEN y, long N, long prec)
     401             : {
     402         196 :   GEN v, y1 = gsubgs(gmulsg(2, y), 1), y3 = gmul(y, gsubsg(1, y));
     403             :   long n;
     404             : 
     405         196 :   v = cgetg(N + 2, t_VEC); gel(v, N + 1) = gen_0;
     406         196 :   if (gcmpgs(gnorm(y3),1) > 0)
     407             :   {
     408         175 :     GEN y2 = gdiv(r1, y3);
     409       21378 :     for (n = N; n >= 1; n--)
     410             :     {
     411       21203 :       pari_sp av2 = avma;
     412       21203 :       GEN z = gmul(y2, gsub(gel(v, n+1), gmul(y1, gel(ibin1, n+1))));
     413       21203 :       gel(v, n) = gerepileupto(av2, z);
     414             :     }
     415             :   }
     416             :   else
     417             :   {
     418          21 :     pari_sp av0 = avma;
     419          21 :     gel(v, 1) = gerepileupto(av0, glog(gdiv(y, gsubgs(y, 1)), prec));
     420        2247 :     for (n = 1; n < N; n++)
     421             :     {
     422        2226 :       pari_sp av2 = avma;
     423        2226 :       GEN z = gadd(gmul(y3, gel(v, n)), gmul(y1, gel(ibin1, n+1)));
     424        2226 :       gel(v, n + 1) = gerepileupto(av2, z);
     425             :     }
     426             :   }
     427         196 :   return v;
     428             : }
     429             : static GEN
     430        9478 : fillrec(hashtable *H, GEN evec, long _0, long _1, GEN X, GEN pab, GEN r1,
     431             :         long N, long prec)
     432             : {
     433             :   long n, a, b, s, x0;
     434             :   GEN xy1, x, y, r, wmid, wini, wfin, mid, ini, fin;
     435        9478 :   hashentry *ep = hash_search(H, evec);
     436             : 
     437        9478 :   if (ep) return (GEN)ep->val;
     438        3108 :   findabvgen(evec, _0, _1, &wmid, &wini, &wfin, &a, &b);
     439        3108 :   x = gel(X, evec[1]); s = lg(evec)-1; /* > 1 */
     440        3108 :   y = gel(X, evec[s]);
     441        3108 :   mid = fillrec(H, wmid, _0, _1, X, pab, r1, N, prec);
     442        3108 :   ini = fillrec(H, wini, _0, _1, X, pab, r1, N, prec);
     443        3108 :   fin = fillrec(H, wfin, _0, _1, X, pab, r1, N, prec);
     444        3108 :   if (evec[1] == _0) { x0 = 1; xy1 = gdiv(r1, y); }
     445         525 :   else { x0 = 0; xy1 = gdiv(r1, gmul(gsubsg(1, x), y)); }
     446        3108 :   r = cgetg(N+2, t_VEC); gel(r, N+1) = gen_0;
     447      398006 :   for (n = N; n > 1; n--)
     448             :   {
     449      394898 :     pari_sp av = avma;
     450      394898 :     GEN t = gmul(gel(ini, n+1), gmael(pab, n, a));
     451      394898 :     GEN u = gadd(gmul(gel(fin, n+1), gmael(pab, n, b)), gel(mid, n+1));
     452      394898 :     GEN v = gdiv(x0? gadd(t, u): gsub(t, u), gmael(pab, n, a+b));
     453      394898 :     gel(r, n) = gerepileupto(av, gmul(xy1, gadd(gel(r, n+1), v)));
     454             :   }
     455             :   { /* n = 1 */
     456        3108 :     pari_sp av = avma;
     457        3108 :     GEN t = gel(ini, 2), u = gadd(gel(fin, 2), gel(mid, 2));
     458        3108 :     GEN v = x0? gadd(t, u): gsub(t, u);
     459        3108 :     gel(r,1) = gerepileupto(av, gmul(xy1, gadd(gel(r,2), v)));
     460             :   }
     461        3108 :   hash_insert(H, (void*)evec, (void*)r); return r;
     462             : }
     463             : 
     464             : static GEN
     465         161 : aztoe(GEN avec, GEN zvec, long prec)
     466             : {
     467         161 :   GEN y, E, u = subsr(1, real2n(10-prec2nbits(prec), LOWDEFAULTPREC));
     468         161 :   long i, l = lg(avec);
     469             : 
     470         161 :   E = cgetg(l, t_VEC); if (l == 1) return E;
     471         161 :   y = gen_1;
     472         609 :   for (i = 1; i < l; i++)
     473             :   {
     474         455 :     long a = avec[i];
     475         455 :     GEN e, zi = gel(zvec, i);
     476         455 :     if (a <= 0 || (a == 1 && i == 1 && gequal1(zi)))
     477           7 :       pari_err_TYPE("polylogmult [divergent]", avec);
     478         448 :     if (gequal0(zi)) return NULL;
     479         448 :     gel(E, i) = e = zerovec(a);
     480         448 :     gel(e, a) = y = gdiv(y, zi);
     481         448 :     if (gcmp(gnorm(y), u) < 0) pari_err_TYPE("polylogmult [divergent]", zvec);
     482             :   }
     483         154 :   return shallowconcat1(E);
     484             : }
     485             : 
     486             : /***********************************************************/
     487             : /* Special case of zvec = [1,1,...], i.e., zetamult.       */
     488             : /***********************************************************/
     489             : static void
     490         812 : findabvgens(GEN evec, GEN *pwmid, GEN *pwini, GEN *pwfin, long *pa, long *pb)
     491             : {
     492             :   GEN wmid, wini, wfin;
     493         812 :   long s = lg(evec) - 1, a, b, j, m;
     494         812 :   if (s == 2)
     495             :   {
     496          70 :     *pwmid = cgetg(1, t_VECSMALL);
     497          70 :     *pwini = mkvecsmall(0);
     498          70 :     *pwfin = mkvecsmall(1);
     499          70 :     *pa = *pb = 1; return;
     500             :   }
     501         742 :   a = s - 1;
     502        1330 :   for (j = 1; j <= s - 2; j++) if (!evec[j + 1]) { a = j; break; }
     503         742 :   *pa = a;
     504         742 :   b = s - 1;
     505        1442 :   for (j = s - 2; j >= 1; j--) if (evec[j + 1]) { b = s - 1 - j; break; }
     506         742 :   *pb = b;
     507             : 
     508         742 :   *pwmid = wmid = a+b < s? vecslice(evec, a+1, s-b): cgetg(1, t_VECSMALL);
     509         742 :   m = lg(wmid) - 1;
     510         742 :   *pwini = wini = cgetg(a + m + 1, t_VECSMALL);
     511        1330 :   wini[1] = 0; for (j = 2; j <= a; j++) wini[j] = 1;
     512        4704 :   for (; j <= a + m; j++) wini[j] = wmid[j-a];
     513         742 :   *pwfin = wfin = cgetg(b + m + 1, t_VECSMALL);
     514        4704 :   for (j = 1; j <= m; j++) wfin[j] = wmid[j];
     515        1442 :   for (; j < b + m; j++) wfin[j] = 0;
     516         742 :   wfin[j] = 1;
     517             : }
     518             : static GEN
     519        2506 : fillrecs(hashtable *H, GEN evec, GEN pab, long N, long prec)
     520             : {
     521             :   long n, a, b;
     522             :   GEN r, wmid, wini, wfin, mid, ini, fin;
     523        2506 :   hashentry *ep = hash_search(H, evec);
     524             : 
     525        2506 :   if (ep) return (GEN)ep->val;
     526         812 :   findabvgens(evec, &wmid, &wini, &wfin, &a, &b);
     527         812 :   mid = fillrecs(H, wmid, pab, N, prec);
     528         812 :   ini = fillrecs(H, wini, pab, N, prec);
     529         812 :   fin = fillrecs(H, wfin, pab, N, prec);
     530         812 :   r = cgetg(N + 2, t_VEC); gel(r, N+1) = gen_0;
     531      104748 :   for (n = N; n > 1; n--)
     532             :   {
     533      103936 :     GEN z = cgetr(prec);
     534      103936 :     pari_sp av = avma;
     535      103936 :     GEN t = mpmul(gel(ini, n+1), gmael(pab, n, a));
     536      103936 :     GEN u = mpadd(mpmul(gel(fin, n+1), gmael(pab, n, b)), gel(mid,n+1));
     537      103936 :     GEN v = mpdiv(mpadd(t, u), gmael(pab, n, a+b));
     538      103936 :     mpaff(mpadd(gel(r, n+1), v), z); set_avma(av); gel(r,n) = z;
     539             :   }
     540             :   { /* n = 1 */
     541         812 :     GEN z = cgetr(prec);
     542         812 :     pari_sp av = avma;
     543         812 :     GEN t = gel(ini,2), u = mpadd(gel(fin,2), gel(mid,2)), v = mpadd(t, u);
     544         812 :     mpaff(mpadd(gel(r, 2), v), z); set_avma(av); gel(r,1) = z;
     545             :   }
     546         812 :   hash_insert(H, (void*)evec, (void*)r); return r;
     547             : }
     548             : /* [n, ..., n^k] */
     549             : static GEN
     550       39522 : powersu(ulong n, long k)
     551             : {
     552       39522 :   GEN v, gn = utoipos(n);
     553       39522 :   long i, l = k+1;
     554       39522 :   v = cgetg(l, t_VEC); gel(v,1) = gn;
     555      420280 :   for (i = 2; i < l; i++) gel(v,i) = muliu(gel(v,i-1), n);
     556       39522 :   return v;
     557             : }
     558             : /* n^a = pab[n][a] */
     559             : static GEN
     560         329 : get_pab(long N, long k)
     561             : {
     562         329 :   GEN v = cgetg(N+1, t_VEC);
     563             :   long j;
     564         329 :   gel(v, 1) = gen_0; /* not needed */
     565       39851 :   for (j = 2; j <= N; j++) gel(v, j) = powersu(j, k);
     566         329 :   return v;
     567             : }
     568             : static hashtable *
     569         224 : zetamult_hash(long _0, long _1, GEN ibin, GEN ibin1)
     570             : {
     571         224 :   hashtable *H = hash_create(4096, (ulong(*)(void*))&hash_zv,
     572             :                                    (int(*)(void*,void*))&zv_equal, 1);
     573         224 :   hash_insert(H, (void*)cgetg(1, t_VECSMALL), (void*)ibin);
     574         224 :   hash_insert(H, (void*)mkvecsmall(_0), (void*)ibin1);
     575         224 :   hash_insert(H, (void*)mkvecsmall(_1), (void*)ibin1); return H;
     576             : }
     577             : /* Akhilesh recursive algorithm, #a > 1;
     578             :  * e t_VECSMALL, prec final precision, bit required bitprecision */
     579             : static GEN
     580          70 : zetamult_Akhilesh(GEN e, long bit, long prec)
     581             : {
     582          70 :   long k = lg(e) - 1, N = 1 + bit/2, prec2 = nbits2prec(bit);
     583             :   GEN r, pab, ibin, ibin1;
     584             :   hashtable *H;
     585             : 
     586          70 :   get_ibin(&ibin, &ibin1, N, prec2);
     587          70 :   pab = get_pab(N, k);
     588          70 :   H = zetamult_hash(0, 1, ibin, ibin1);
     589          70 :   r = fillrecs(H, e, pab, lg(pab)-1, prec2);
     590          70 :   if (DEBUGLEVEL) err_printf("polylogmult: k = %ld, %ld nodes\n", k, H->nb);
     591          70 :   return gprec_wtrunc(gel(r,1), prec);
     592             : }
     593             : /* evec t_VEC */
     594             : static GEN
     595         154 : zetamultevec(GEN evec, long prec)
     596             : {
     597         154 :   pari_sp av = avma;
     598         154 :   double *x, *y, z = 0;
     599         154 :   long i, j, l, lH, bitprec, prec2, N, _0 = 0, _1 = 0, k = lg(evec) - 1;
     600             :   GEN r1, r, pab, ibin, ibin1, X, Evec, v;
     601             :   hashtable *H;
     602             : 
     603         154 :   if (k == 0) return gen_1;
     604         154 :   v = vec_equiv(evec); l = lg(v);
     605         154 :   Evec = cgetg(k+1, t_VECSMALL);
     606         154 :   X = cgetg(l + 2, t_VEC); /* our alphabet */
     607         490 :   for (i = lH = 1; i < l; i++)
     608             :   {
     609         336 :     GEN vi = gel(v,i), xi = gel(evec, vi[1]);
     610         336 :     long li = lg(vi);
     611         336 :     gel(X,i) = xi;
     612         336 :     if (!_0 && gequal0(xi)) _0 = i;
     613         196 :     else if (!_1 && gequal1(xi)) _1 = i;
     614        2338 :     for (j = 1; j < li; j++) Evec[vi[j]] = i;
     615             :   }
     616             :   /* add 0,1 if needed */
     617         154 :   if (!_0) { gel(X, i) = gen_0; _0 = i++; }
     618         154 :   if (!_1) { gel(X, i) = gen_1; _1 = i++; }
     619         154 :   l = i; setlg(X, l);
     620         154 :   av = avma;
     621         154 :   x = (double*)stack_malloc_align(l * sizeof(double), sizeof(double));
     622         154 :   y = (double*)stack_malloc_align(l * sizeof(double), sizeof(double));
     623         658 :   for (j = 1; j < l; j++)
     624             :   {
     625         504 :     GEN t = gel(X,j);
     626         504 :     x[j] = (j == _1)? 0: -dbllog2(gsubsg(1, t));
     627         504 :     y[j] = (j == _0)? 0: -dbllog2(t);
     628             :   }
     629         658 :   for (i = 1; i < l; i++)
     630        1120 :     for (j = i+1; j < l; j++) z = maxdd(z, x[i] + y[j]);
     631         154 :   set_avma(av);
     632         154 :   if (z >= 2) pari_err_IMPL("polylogmult in this range");
     633         154 :   bitprec = prec2nbits(prec) + 64*(1 + (k >> 5));
     634         154 :   N = 1 + bitprec/ (2 - z);
     635         154 :   bitprec += z * N;
     636         154 :   prec2 = nbits2prec(bitprec);
     637         154 :   X = gprec_wensure(X, prec2);
     638         154 :   get_ibin(&ibin, &ibin1, N, prec2);
     639         154 :   pab = get_pab(N, k);
     640         154 :   H = zetamult_hash(_0, _1, ibin, ibin1);
     641         154 :   r1 = real_1(prec2);
     642         658 :   for (i = 1; i < l; i++)
     643         504 :     if (i != _0 && i != _1)
     644         196 :       hash_insert(H, mkvecsmall(i), filllg1(ibin1, r1, gel(X,i), N, prec2));
     645         154 :   r = fillrec(H, Evec, _0, _1, X, pab, r1, N, prec2);
     646         154 :   if (DEBUGLEVEL) err_printf("polylogmult: k = %ld, %ld nodes\n", k, H->nb);
     647         154 :   return gprec_wtrunc(gel(r,1), prec);
     648             : }
     649             : 
     650             : /* a t_VECSMALL */
     651             : static GEN
     652         147 : zetamult_i(GEN a, long prec)
     653             : {
     654         147 :   long r = lg(a)-1, k, bit;
     655         147 :   if (r == 0) return gen_1;
     656         147 :   if (r == 1) return szeta(a[1], prec);
     657         140 :   bit = prec2nbits(prec);
     658         140 :   if (bit <= 128)
     659          42 :     return zetamult_Zagier(a, bit, prec + EXTRAPRECWORD);
     660          98 :   k = zv_sum(a);
     661          98 :   if (((double)r) / (k*k) * bit / log((double)10*bit) < 0.5)
     662          28 :     return zetamult_Zagier(a, bit, prec + EXTRAPRECWORD);
     663          70 :   bit += maxss(log2zeta_bound(a), 64);
     664          70 :   return zetamult_Akhilesh(atoe(a), bit, prec);
     665             : }
     666             : GEN
     667         203 : zetamult(GEN s, long prec)
     668             : {
     669         203 :   pari_sp av0 = avma, av;
     670         203 :   GEN z, avec, r = cgetr(prec);
     671             : 
     672         203 :   av = avma; avec = zetamultconvert_i(s,1);
     673         154 :   if (lg(avec) == 1) return gc_const(av0, gen_1);
     674         147 :   z = zetamult_i(avec, prec); affrr(z, r); return gc_const(av, r);
     675             : }
     676             : 
     677             : /* If star = NULL: MZV, otherwise Yamamoto interpolation (MZSV for t=1) */
     678             : GEN
     679         210 : zetamult_interpolate(GEN s, GEN t, long prec)
     680             : {
     681         210 :   pari_sp av = avma, av2;
     682             :   long i, k, l, la, bit;
     683             :   GEN avec, v, V;
     684             : 
     685         210 :   if (!t) return zetamult(s, prec);
     686          14 :   avec = zetamultconvert_i(s, 1); k = zv_sum(avec);
     687          14 :   bit = prec2nbits(prec);
     688          14 :   if (bit <= 128 || k > 20 || (bit >> k) < 4)
     689          14 :     return zetamult_interpolate2_i(vecsmall_reverse(avec), t, prec);
     690           0 :   v = allstar(avec); l = lg(v); la = lg(avec);
     691           0 :   V = cgetg(la, t_VEC);
     692           0 :   for (i = 1; i < la; i++)
     693           0 :   { gel(V,i) = cgetr(prec + EXTRAPRECWORD); affur(0, gel(V,i)); }
     694           0 :   av2 = avma;
     695           0 :   for (i = 1; i < l; i++, set_avma(av2))
     696             :   {
     697           0 :     GEN a = gel(v,i); /* avec */
     698           0 :     long n = lg(a)-1; /* > 0 */
     699           0 :     affrr(addrr(gel(V,n), zetamult_i(a, prec)), gel(V,n));
     700             :   }
     701           0 :   return gerepileupto(av, poleval(vecreverse(V),t));
     702             : }
     703             : 
     704             : 
     705             : GEN
     706          63 : polylogmult(GEN a, GEN z, long prec)
     707             : {
     708          63 :   pari_sp av = avma;
     709          63 :   if (!z) return zetamult(a, prec);
     710          56 :   switch(typ(a))
     711             :   {
     712          56 :     case t_VEC:
     713          56 :     case t_COL: a = gtovecsmall(a); break;
     714           0 :     case t_VECSMALL: break;
     715           0 :     default: pari_err_TYPE("polylogmult", a);
     716             :              return NULL;/*LCOV_EXCL_LINE*/
     717             :   }
     718          56 :   switch (typ(z))
     719             :   {
     720           0 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     721           0 :       z = mkvec(z); break;
     722          56 :     case t_VEC: case t_COL: break;
     723           0 :     case t_VECSMALL: z = zv_to_ZV(z); break;
     724           0 :     default: pari_err_TYPE("polylogmult [z]", z);
     725             :   }
     726          56 :   if (lg(z) != lg(a))
     727           7 :     pari_err_TYPE("polylogmult [#s != #z]", mkvec2(a,z));
     728          49 :   return gerepilecopy(av, zetamultevec(aztoe(a,z,prec), prec));
     729             : }
     730             : 
     731             : GEN
     732          70 : polylogmult_interpolate(GEN s, GEN zvec, GEN t, long prec)
     733             : {
     734          70 :   pari_sp av = avma;
     735             :   GEN V, avec, A, AZ, Z;
     736             :   long i, la, l;
     737             : 
     738          70 :   if (!t) return polylogmult(s, zvec, prec);
     739           7 :   if (!zvec) return zetamult_interpolate(s, t, prec);
     740           7 :   avec = zetamultconvert_i(s, 1); la = lg(avec);
     741           7 :   AZ = allstar2(avec, zvec);
     742           7 :   A = gel(AZ, 1); l = lg(A);
     743           7 :   Z = gel(AZ, 2); V = zerovec(la-1);
     744         119 :   for (i = 1; i < l; i++)
     745             :   {
     746         112 :     pari_sp av2 = avma;
     747         112 :     GEN a = gel(A,i), e = aztoe(a, gel(Z,i), prec);
     748         112 :     long n = lg(a)-1; /* > 0 */
     749         112 :     gel(V,n) = gerepileupto(av2, gadd(gel(V,n), zetamultevec(e, prec)));
     750             :   }
     751           7 :   return gerepileupto(av, poleval(vecreverse(V),t));
     752             : }
     753             : 
     754             : /**************************************************************/
     755             : /*                           ALL MZV's                        */
     756             : /**************************************************************/
     757             : 
     758             : /* Given admissible evec w = 0e_2....e_{k-1}1, compute a,b,v such that
     759             :  * w=0{1}_{b-1}v{0}_{a-1}1 with v empty or admissible.
     760             :  * Input: binary vector evec */
     761             : static void
     762         945 : findabv(GEN w, long *pa, long *pb, long *pminit, long *pmmid, long *pmfin)
     763             : {
     764         945 :   long le = lg(w) - 2;
     765         945 :   if (le == 0)
     766             :   {
     767         105 :     *pa = 1;
     768         105 :     *pb = 1;
     769         105 :     *pminit = 2;
     770         105 :     *pmfin = 2;
     771         105 :     *pmmid = 1;
     772             :   }
     773             :   else
     774             :   {
     775             :     long a, b, j, lv;
     776        1050 :     for (j = 1; j <= le; j++)
     777        1050 :       if (!w[j+1]) break;
     778         840 :     *pb = b = j;
     779        1890 :     for (j = le; j >= 1; j--)
     780        1575 :       if (w[j+1]) break;
     781         840 :     *pa = a = le + 1 - j;
     782         840 :     lv = le + 2 - a - b;
     783         840 :     if (lv > 0)
     784             :     {
     785         315 :       long v = fd(w, b + 1, le - a + 2), u = v + (1 << (lv-1));
     786         315 :       *pminit = (((1 << b) - 1) << (lv - 1)) + (v/2) + 2;
     787         315 :       *pmfin = (u << (a - 1)) + 2;
     788         315 :       *pmmid = (u >> 1) + 2;
     789             :     }
     790             :     else
     791             :     {
     792         525 :       *pminit = (1 << (b - 1)) + 1;
     793         525 :       *pmfin = (a == 1) ? 2 : (1 << (a - 2)) + 2;
     794         525 :       *pmmid = 1;
     795             :     }
     796             :   }
     797         945 : }
     798             : 
     799             : /* Returns L:
     800             : * L[1] contains zeta(emptyset)_{n-1,n-1},
     801             : * L[2] contains zeta({0})_{n-1,n-1}=zeta({1})_{n-1,n-1} for n >= 2,
     802             : * L[m+2][n] : 1 <= m < 2^{k-2}, 1 <= n <= N + 1
     803             : * contains zeta(w)_{n-1,n-1}, w corresponding to m,n
     804             : * L[m+2] : 2^{k-2} <= m < 2^{k-1} contains zeta(w), w corresponding to m
     805             : (code: w=0y1 iff m=1y). */
     806             : static GEN
     807         105 : fillL(long k, long bitprec)
     808             : {
     809         105 :   long N = 1 + bitprec/2, prec = nbits2prec(bitprec);
     810         105 :   long s, j, n, m, K = 1 << (k - 1), K2 = K/2;
     811         105 :   GEN p1, p2, pab = get_pab(N, k), L = cgetg(K + 2, t_VEC);
     812             : 
     813         105 :   get_ibin(&gel(L,1), &gel(L,2), N, prec);
     814         840 :   for (m = 1; m < K2; m++)
     815             :   {
     816         735 :     gel(L, m+2) = p1 = cgetg(N+1, t_VEC);
     817       83055 :     for (n = 1; n < N; n++) gel(p1, n) = cgetr(prec);
     818         735 :     gel(p1, n) = gen_0;
     819             :   }
     820         945 :   for (m = K2; m < K; m++) gel(L, m+2) = utor(0, prec);
     821         525 :   for (s = 2; s <= k; s++)
     822             :   { /* Assume length evec < s filled */
     823             :     /* If evec = 0e_2...e_{s-1}1 then m = (1e_2...e_{s-1})_2 */
     824         420 :     GEN w = cgetg(s, t_VECSMALL);
     825         420 :     long M = 1 << (s - 2);
     826         420 :     pari_sp av = avma;
     827        1995 :     for (m = M; m < 2*M; m++)
     828             :     {
     829             :       GEN pinit, pfin, pmid;
     830             :       long comp, a, b, mbar, minit, mfin, mmid, mc;
     831        1575 :       p1 = gel(L, m + 2);
     832        5145 :       for (j = s - 1, mc = m, mbar = 1; j >= 2; j--, mc >>= 1)
     833             :       {
     834        3570 :         w[j] = mc & 1;
     835        3570 :         mbar = (1 - w[j]) | (mbar << 1);
     836             :       }
     837             :       /* m, mbar are dual; handle smallest, copy the other */
     838        1575 :       comp = mbar - m; if (comp < 0) continue; /* m > mbar */
     839         945 :       if (comp)
     840             :       {
     841         630 :         p2 = gel(L, mbar + 2);
     842         630 :         setisclone(p2); /* flag as dual */
     843             :       }
     844             :       else
     845         315 :         p2 = NULL; /* no copy needed if m = mbar */
     846         945 :       findabv(w, &a,&b,&minit,&mmid,&mfin);
     847         945 :       pinit= gel(L, minit);
     848         945 :       pfin = gel(L, mfin);
     849         945 :       pmid = gel(L, mmid);
     850      105840 :       for (n = N-1; n > 1; n--, set_avma(av))
     851             :       {
     852      104895 :         GEN t = mpmul(gel(pinit,n+1), gmael(pab, n, b));
     853      104895 :         GEN u = mpmul(gel(pfin, n+1), gmael(pab, n, a));
     854      104895 :         GEN v = gel(pmid, n+1), S = s < k ? gel(p1, n+1): p1;
     855      104895 :         S = mpadd(S, mpdiv(mpadd(mpadd(t, u), v), gmael(pab, n, a+b)));
     856      104895 :         mpaff(S, s < k ? gel(p1, n) : p1);
     857      104895 :         if (p2 && s < k) mpaff(S, gel(p2, n));
     858             :       }
     859             :       { /* n = 1: same formula simplifies */
     860         945 :         GEN t = gel(pinit,2), u = gel(pfin,2), v = gel(pmid,2);
     861         945 :         GEN S = s < k ? gel(p1,2): p1;
     862         945 :         S = mpadd(S, mpadd(mpadd(t, u), v));
     863         945 :         mpaff(S, s < k ? gel(p1,1) : p1);
     864         945 :         if (p2 && s < k) mpaff(S, gel(p2, 1));
     865         945 :         set_avma(av);
     866             :       }
     867         945 :       if (p2 && s == k) mpaff(p1, p2);
     868             :     }
     869             :   }
     870         105 :   return L;
     871             : }
     872             : 
     873             : /* bit 1 of flag unset: full, otherwise up to duality (~ half)
     874             :  * bit 2 of flag unset: all <= k, otherwise only k
     875             :  * half: 2^(k-3)+ delta_{k even} * 2^(k/2-2), sum = 2^(k-2)+2^(floor(k/2)-1)-1
     876             :  * full: 2^(k-2); sum = 2^(k-1)-1 */
     877             : static GEN
     878         105 : zetamultall_i(long k, long flag, long prec)
     879             : {
     880         105 :   GEN res, ind, L = fillL(k, prec2nbits(prec) + 32);
     881         105 :   long m, K2 = 1 << (k-2), n = lg(L) - 1, m0 = (flag & 4L) ? K2 : 1;
     882             : 
     883         105 :   if (!(flag & 2L))
     884             :   {
     885          77 :     res = cgetg(n - m0, t_VEC);
     886          77 :     ind = cgetg(n - m0, t_VECSMALL);
     887         938 :     for (m = m0; m < n - 1; m++)
     888             :     {
     889         861 :       gel(res, m - m0 + 1) = m < K2 ? gmael(L, m + 2, 1) : gel(L, m + 2);
     890         861 :       ind[m - m0 + 1] = m;
     891             :     }
     892             :   }
     893             :   else
     894             :   { /* up to duality */
     895             :     long nres, c;
     896          28 :     if (k == 2) nres = 1;
     897          28 :     else if (!(flag & 2L))
     898           0 :       nres = (1 << (k - 2)) + (1 << ((k/2) - 1)) - 1;
     899             :     else
     900          28 :       nres = (1 << (k - 1));
     901          28 :     res = cgetg(nres + 1, t_VEC);
     902          28 :     ind = cgetg(nres + 1, t_VECSMALL);
     903         350 :     for (m = m0, c = 1; m < n - 1; m++)
     904             :     {
     905         322 :       GEN z = gel(L,m+2);
     906         322 :       if (isclone(z)) continue; /* dual */
     907         182 :       if (m < K2) z = gel(z,1);
     908         182 :       gel(res, c) = z;
     909         182 :       ind[c] = m; c++;
     910             :     }
     911          28 :     setlg(res, c);
     912          28 :     setlg(ind, c);
     913             :   }
     914         105 :   return mkvec2(res, ind);
     915             : }
     916             : 
     917             : /* fd(e, 2, lg(e)-2), e = atoe(avec) */
     918             : static long
     919        1876 : atom(GEN avec)
     920             : {
     921        1876 :   long i, m, l = lg(avec);
     922        1876 :   if (l < 3) return 0;
     923        1232 :   m = 1; /* avec[1] != 0 */
     924        1708 :   for (i = 2; i < l-1; i++) m = (m << avec[i]) + 1;
     925        1232 :   return m << (avec[i]-1);
     926             : }
     927             : static long
     928        1876 : atoind(GEN avec, long flag)
     929        1876 : { return atom(avec) + (flag? 1: (1 << (zv_sum(avec) - 2))); }
     930             : /* If flag is unset, L has all k1 <= k, otherwise only k */
     931             : static GEN
     932         644 : zetamultstar_i(GEN L, GEN avec, long flag)
     933             : {
     934         644 :   GEN s = allstar(avec), S = gen_0;
     935         644 :   long i, l = lg(s);
     936        2520 :   for (i = 1; i < l; i++) S = gadd(S, gel(L, atoind(gel(s,i), flag)));
     937         644 :   return S;
     938             : }
     939             : 
     940             : /* bit 0: notstar/star
     941             :  * bit 1: full/half (ignored if star, always full)
     942             :  * bit 2: all <= k / only k
     943             :  * bit 3: without / with index */
     944             : GEN
     945         119 : zetamultall(long k, long flag, long prec)
     946             : {
     947         119 :   pari_sp av = avma;
     948             :   GEN Lind, L, res;
     949             :   long K, k1, ct, fl;
     950             : 
     951         119 :   if (flag < 0 || flag > 15) pari_err_FLAG("zetamultall");
     952         119 :   if (k < 1) pari_err_DOMAIN("zetamultall", "k", "<", gen_1, stoi(k));
     953         112 :   if (k >= 64) pari_err_OVERFLOW("zetamultall");
     954         105 :   if (!(flag & 1L))
     955             :   { /* not star */
     956          49 :     Lind = zetamultall_i(k, flag, prec);
     957          49 :     res = (flag & 8L)? Lind : gel(Lind, 1);
     958          49 :     return gerepilecopy(av, res);
     959             :   }
     960             :   /* star */
     961          56 :   fl = flag & 4L; /* 4 if k, else 0 (all <= k) */
     962          56 :   Lind = gerepilecopy(av, zetamultall_i(k, fl, prec)); /* full */
     963          56 :   L = gel(Lind, 1);
     964          56 :   K = 1 << (k - 2);
     965          56 :   res = cgetg(fl? K+1: 2*K, t_VEC);
     966         196 :   for (ct = 1, k1 = fl? k: 2; k1 <= k; k1++)
     967             :   {
     968         140 :     GEN w = cgetg(k1 + 1, t_VECSMALL);
     969         140 :     long M = 1 << (k1 - 1), m;
     970         784 :     for (m = 1; m <= M; m += 2)
     971             :     {
     972         644 :       pari_sp av = avma;
     973         644 :       long j, mc = m;
     974        3556 :       for (j = k1; j >= 1; j--) { w[j] = mc & 1; mc >>= 1; }
     975         644 :       gel(res, ct++) = gerepileupto(av, zetamultstar_i(L, etoa(w), fl));
     976             :     }
     977             :   }
     978          56 :   if (flag & 8L) res = mkvec2(res, gel(Lind, 2));
     979          56 :   return gerepilecopy(av, res);
     980             : }

Generated by: LCOV version 1.13