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

Generated by: LCOV version 1.13