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 : }
|
