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 14693 : zetamultconvert_i(GEN a, long fl)
89 : {
90 : long i, l;
91 14693 : if (fl < 0 || fl > 2) pari_err_FLAG("zetamultconvert");
92 14693 : 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 7483 : case t_VEC: case t_COL: case t_VECSMALL:
104 7483 : a = gtovecsmall(a);
105 7483 : l = lg(a);
106 7483 : 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 : }
117 : else
118 : {
119 7434 : if (a[1] < 2) pari_err_TYPE("zetamultconvert", a);
120 37205 : for (i = 2; i < l; i++)
121 29792 : 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 : }
128 7441 : break;
129 7 : default: pari_err_TYPE("zetamultconvert", a);
130 : }
131 14644 : 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 28 : allstar2(GEN avec, GEN zvec)
177 : {
178 28 : long la = lg(avec), i, K = 1 << (la - 2);
179 28 : GEN W = cgetg(K + 1, t_VEC), Z = cgetg(K + 1, t_VEC);
180 :
181 28 : gel(W, 1) = avec;
182 28 : gel(Z, 1) = zvec = gtovec(zvec);
183 77 : for (i = 2; i < la; i++)
184 : {
185 49 : long L = 1 << (i - 2), j;
186 175 : for (j = 1; j <= L; j++)
187 : {
188 126 : GEN u = gel(W, j), w, y = gel(Z, j), z;
189 126 : long l = lg(u) - 1, ind = l - la + i, k;
190 126 : w = cgetg(l, t_VECSMALL);
191 126 : z = cgetg(l, t_VEC);
192 245 : for (k = 1; k < ind; k++) { w[k] = u[k]; gel(z, k) = gel(y, k); }
193 126 : w[ind] = u[ind] + u[ind + 1];
194 126 : gel(z, ind) = gmul(gel(y, ind), gel(y, ind + 1));
195 203 : for (k = ind + 1; k < l; k++) { w[k] = u[k+1]; gel(z, k) = gel(y, k+1); }
196 126 : gel(W, L + j) = w;
197 126 : gel(Z, L + j) = z;
198 : }
199 : }
200 28 : 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 182 : zparams(long *s, long *l, long b)
209 : {
210 182 : double D = b * LOG10_2, E = 3*D/2 / log(3*D);
211 182 : *s = (long)floor(E*E);
212 182 : *l = (long)floor(sqrt(*s * log((double)*s)/2.));
213 182 : }
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) = utor(0, prec);
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 42 : 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 42 : zparams(&s, &l, prec2nbits(prec));
269 42 : n = lg(avec) - 1;
270 42 : a = zeromatcopy(n + 1, l);
271 42 : b = zerovec(l + 1);
272 196 : for (i = 1; i <= n; i++)
273 : {
274 154 : long h = avec[n + 1 - i];
275 154 : if (i == 1) gel(b, h) = gen_1;
276 : else
277 2352 : for (j = l + 1; j >= 1; j--)
278 : {
279 2240 : if (j <= h) gel(b, j) = gen_0;
280 1960 : else gel(b, j) = gadd(gcoeff(a, i, j-h), gmul(t, gel(b, j-h)));
281 : }
282 154 : gcoeff(a, i+1, 1) = gel(b, 2); /* j = 1 */
283 2926 : for (j = 2; j <= l; j++)
284 : { /* b[j+1] - sum_{0 <= u < j-1} binom(j,u) a[i+1,u+1]*/
285 2772 : pari_sp av = avma;
286 2772 : GEN S = gel(b, j + 1);
287 2772 : S = gsub(S, gcoeff(a, i+1, 1)); /* u = 0 */
288 2772 : if (j > 2) S = gsub(S, gmulgu(gcoeff(a, i+1, 2), j)); /* u = 1 */
289 2772 : if (j >= 4)
290 : {
291 2464 : GEN C = utoipos(j*(j-1) / 2);
292 2464 : long u, U = (j-1)/2;
293 12320 : for (u = 2; u <= U; u++)
294 : { /* C = binom(j, u) = binom(j, j-u) */
295 9856 : GEN A = gadd(gcoeff(a, i+1, u+1), gcoeff(a, i+1, j-u+1));
296 9856 : S = gsub(S, gmul(C, A));
297 9856 : C = diviuexact(muliu(C, j-u), u+1);
298 : }
299 2464 : if (!odd(j)) S = gsub(S, gmul(C, gcoeff(a,i+1, j/2+1)));
300 : }
301 2772 : gcoeff(a, i+1, j) = gerepileupto(av, gdivgu(S, j));
302 : }
303 : }
304 42 : _1 = real_1(prec + EXTRAPREC64); av = avma;
305 42 : ze = cgetg(n+1, t_VEC);
306 196 : for (i = 1; i <= n; i++)
307 : {
308 154 : GEN S = gdivgu(gcoeff(a, n+2-i, 1), s), sj = utoipos(s);
309 2926 : for (j = 2; j <= l; j++)
310 : {
311 2772 : sj = muliu(sj, s); /* = s^j */
312 2772 : S = gadd(S, gdiv(gcoeff(a, n+2-i, j), sj));
313 : }
314 154 : gel(ze, i) = S;
315 : }
316 6258 : for (i = s; i >= 1; i--)
317 : {
318 6216 : GEN b0 = divri(_1, powuu(i, avec[n])), z;
319 6216 : z = gel(ze,n); gel(ze,n) = gadd(z, b0);
320 22792 : for (j = n-1; j >= 1; j--)
321 : {
322 16576 : b0 = gdiv(gadd(gmul(t, b0), z), powuu(i, avec[j]));
323 16576 : z = gel(ze,j); gel(ze,j) = gadd(gel(ze,j), b0);
324 : }
325 6216 : 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 42 : 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 378 : get_ibin(GEN *pibin, GEN *pibin1, long N, long prec)
344 : {
345 : GEN ibin, ibin1;
346 : long n;
347 378 : *pibin = ibin = cgetg(N + 2, t_VEC);
348 378 : *pibin1= ibin1= cgetg(N + 2, t_VEC);
349 378 : gel(ibin,1) = gel(ibin1,1) = gen_0; /* unused */
350 378 : gel(ibin,2) = gel(ibin1,2) = real2n(-1,prec);
351 50792 : for (n = 2; n <= N; n++)
352 : {
353 50414 : gel(ibin, n+1) = divru(mulru(gel(ibin, n), n), 4*n-2);
354 50414 : gel(ibin1, n+1) = divru(gel(ibin, n+1), n);
355 : }
356 378 : }
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 3234 : findabvgen(GEN evec, long _0, long _1, GEN *pwmid, GEN *pwini, GEN *pwfin,
372 : long *pa, long *pb)
373 : {
374 3234 : long s = lg(evec) - 1, m, a, b, j, x = evec[1], y = evec[s];
375 : GEN wmid, wini, wfin;
376 3234 : if (s == 2)
377 : {
378 287 : *pwmid = cgetg(1, t_VECSMALL);
379 287 : *pwini = mkvecsmall(x);
380 287 : *pwfin = mkvecsmall(y);
381 287 : *pa = *pb = 1; return;
382 : }
383 2947 : a = s - 1;
384 2968 : for (j = 1; j <= s - 2; j++) if (evec[j + 1] != _1) { a = j; break; }
385 2947 : *pa = a;
386 2947 : b = s - 1;
387 12544 : for (j = s - 2; j >= 1; j--) if (evec[j + 1] != _0) { b = s - 1 - j; break; }
388 2947 : *pb = b;
389 :
390 2947 : *pwmid = wmid = a+b < s? vecslice(evec, a+1, s-b): cgetg(1, t_VECSMALL);
391 2947 : m = lg(wmid) - 1;
392 2947 : *pwini = wini = cgetg(a + m + 1, t_VECSMALL);
393 2968 : wini[1] = x; for (j = 2; j <= a; j++) wini[j] = _1;
394 11746 : for (; j <= a + m; j++) wini[j] = wmid[j-a];
395 2947 : *pwfin = wfin = cgetg(b + m + 1, t_VECSMALL);
396 11746 : for (j = 1; j <= m; j++) wfin[j] = wmid[j];
397 12544 : for (; j < b + m; j++) wfin[j] = _0;
398 2947 : wfin[j] = y;
399 : }
400 :
401 : /* y != 0,1 */
402 : static GEN
403 231 : filllg1(GEN ibin1, GEN r1, GEN y, long N, long prec)
404 : {
405 231 : GEN v, y1 = gsubgs(gmulsg(2, y), 1), y3 = gmul(y, gsubsg(1, y));
406 : long n;
407 :
408 231 : v = cgetg(N + 2, t_VEC); gel(v, N + 1) = gen_0;
409 231 : if (gcmpgs(gnorm(y3),1) > 0)
410 : {
411 196 : GEN y2 = gdiv(r1, y3);
412 23436 : for (n = N; n >= 1; n--)
413 : {
414 23240 : pari_sp av2 = avma;
415 23240 : GEN z = gmul(y2, gsub(gel(v, n+1), gmul(y1, gel(ibin1, n+1))));
416 23240 : 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 231 : return v;
431 : }
432 : static GEN
433 9905 : 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 9905 : hashentry *ep = hash_search(H, evec);
439 :
440 9905 : if (ep) return (GEN)ep->val;
441 3234 : findabvgen(evec, _0, _1, &wmid, &wini, &wfin, &a, &b);
442 3234 : x = gel(X, evec[1]); s = lg(evec)-1; /* > 1 */
443 3234 : y = gel(X, evec[s]);
444 3234 : mid = fillrec(H, wmid, _0, _1, X, pab, r1, N);
445 3234 : ini = fillrec(H, wini, _0, _1, X, pab, r1, N);
446 3234 : fin = fillrec(H, wfin, _0, _1, X, pab, r1, N);
447 3234 : if (evec[1] == _0) { x0 = 1; xy1 = gdiv(r1, y); }
448 567 : else { x0 = 0; xy1 = gdiv(r1, gmul(gsubsg(1, x), y)); }
449 3234 : r = cgetg(N+2, t_VEC); gel(r, N+1) = gen_0;
450 428792 : for (n = N; n > 1; n--)
451 : {
452 425558 : pari_sp av = avma;
453 425558 : GEN t = gmul(gel(ini, n+1), gmael(pab, n, a));
454 425558 : GEN u = gadd(gmul(gel(fin, n+1), gmael(pab, n, b)), gel(mid, n+1));
455 425558 : GEN v = gdiv(x0? gadd(t, u): gsub(t, u), gmael(pab, n, a+b));
456 425558 : gel(r, n) = gerepileupto(av, gmul(xy1, gadd(gel(r, n+1), v)));
457 : }
458 : { /* n = 1 */
459 3234 : pari_sp av = avma;
460 3234 : GEN t = gel(ini, 2), u = gadd(gel(fin, 2), gel(mid, 2));
461 3234 : GEN v = x0? gadd(t, u): gsub(t, u);
462 3234 : gel(r,1) = gerepileupto(av, gmul(xy1, gadd(gel(r,2), v)));
463 : }
464 3234 : hash_insert(H, (void*)evec, (void*)r); return r;
465 : }
466 :
467 : static GEN
468 210 : aztoe(GEN avec, GEN zvec, long prec)
469 : {
470 210 : GEN y, E, u = subsr(1, real2n(10-prec2nbits(prec), LOWDEFAULTPREC));
471 210 : long i, l = lg(avec);
472 :
473 210 : E = cgetg(l, t_VEC); if (l == 1) return E;
474 210 : y = gen_1;
475 749 : for (i = 1; i < l; i++)
476 : {
477 546 : long a = avec[i];
478 546 : GEN e, zi = gel(zvec, i);
479 546 : if (a <= 0 || (a == 1 && i == 1 && gequal1(zi)))
480 7 : pari_err_TYPE("polylogmult [divergent]", avec);
481 539 : if (gequal0(zi)) return NULL;
482 539 : gel(E, i) = e = zerovec(a);
483 539 : gel(e, a) = y = gdiv(y, zi);
484 539 : if (gcmp(gnorm(y), u) < 0) pari_err_TYPE("polylogmult [divergent]", zvec);
485 : }
486 203 : 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 50414 : powersu(ulong n, long k)
554 : {
555 50414 : GEN v, gn = utoipos(n);
556 50414 : long i, l = k+1;
557 50414 : v = cgetg(l, t_VEC); gel(v,1) = gn;
558 459816 : for (i = 2; i < l; i++) gel(v,i) = muliu(gel(v,i-1), n);
559 50414 : return v;
560 : }
561 : /* n^a = pab[n][a] */
562 : static GEN
563 378 : get_pab(long N, long k)
564 : {
565 378 : GEN v = cgetg(N+1, t_VEC);
566 : long j;
567 378 : gel(v, 1) = gen_0; /* not needed */
568 50792 : for (j = 2; j <= N; j++) gel(v, j) = powersu(j, k);
569 378 : return v;
570 : }
571 : static hashtable *
572 273 : zetamult_hash(long _0, long _1, GEN ibin, GEN ibin1)
573 : {
574 273 : hashtable *H = hash_create(4096, (ulong(*)(void*))&hash_zv,
575 : (int(*)(void*,void*))&zv_equal, 1);
576 273 : hash_insert(H, (void*)cgetg(1, t_VECSMALL), (void*)ibin);
577 273 : hash_insert(H, (void*)mkvecsmall(_0), (void*)ibin1);
578 273 : 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 203 : vec_round(GEN V, long b)
600 : {
601 203 : GEN v = shallowcopy(V), w = shallowcopy(v);
602 203 : long i, j, l = lg(v);
603 2359 : for (i = 1; i < l; i++)
604 : {
605 : long e;
606 2156 : if (!gel(v,i)) continue;
607 1897 : if (gexpo(gsubgs(gel(v,i), 1)) < b) gel(w,i) = gel(v,i) = gen_1;
608 1897 : e = gexpo(gel(v,i));
609 15589 : for (j = i+1; j < l; j++)
610 13692 : 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 203 : return w;
617 : }
618 :
619 : /* evec t_VEC */
620 : static GEN
621 203 : zetamultevec(GEN evec, long prec)
622 : {
623 203 : pari_sp av = avma;
624 203 : double *x, *y, z = 0;
625 203 : 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 203 : if (k == 0) return gen_1;
630 203 : evec = vec_round(evec, 3 - prec2nbits(prec));
631 203 : v = vec_equiv(evec); l = lg(v);
632 203 : Evec = cgetg(k+1, t_VECSMALL);
633 203 : X = cgetg(l + 2, t_VEC); /* our alphabet */
634 665 : for (i = lH = 1; i < l; i++)
635 : {
636 462 : GEN vi = gel(v,i), xi = gel(evec, vi[1]);
637 462 : long li = lg(vi);
638 462 : gel(X,i) = xi;
639 462 : if (!_0 && gequal0(xi)) _0 = i;
640 280 : else if (!_1 && gequal1(xi)) _1 = i;
641 2618 : for (j = 1; j < li; j++) Evec[vi[j]] = i;
642 : }
643 : /* add 0,1 if needed */
644 203 : if (!_0) { gel(X, i) = gen_0; _0 = i++; }
645 203 : if (!_1) { gel(X, i) = gen_1; _1 = i++; }
646 203 : l = i; setlg(X, l);
647 203 : av = avma;
648 203 : x = (double*)stack_malloc_align(l * sizeof(double), sizeof(double));
649 203 : y = (double*)stack_malloc_align(l * sizeof(double), sizeof(double));
650 840 : for (j = 1; j < l; j++)
651 : {
652 637 : GEN t = gel(X,j);
653 637 : x[j] = (j == _1)? 0: -dbllog2(gsubsg(1, t));
654 637 : y[j] = (j == _0)? 0: -dbllog2(t);
655 : }
656 840 : for (i = 1; i < l; i++)
657 1379 : for (j = i+1; j < l; j++) z = maxdd(z, x[i] + y[j]);
658 203 : b = 0;
659 840 : for (i = 1; i < l; i++)
660 : {
661 637 : GEN t = real_i(gel(X,i));
662 637 : long e = -gexpo(gsubgs(gmul2n(t,1), 1));
663 637 : b = maxss(b, e);
664 : }
665 203 : set_avma(av);
666 203 : if (z >= 2) pari_err_IMPL("polylogmult in this range");
667 203 : bitprec = prec2nbits(prec) + 64*(1 + (k >> 5));
668 203 : N = 1 + bitprec / (2 - z);
669 203 : bitprec += z * N;
670 203 : prec2 = nbits2prec(bitprec + b);
671 203 : X = gprec_wensure(X, prec2);
672 203 : get_ibin(&ibin, &ibin1, N, prec2);
673 203 : pab = get_pab(N, k);
674 203 : H = zetamult_hash(_0, _1, ibin, ibin1);
675 203 : r1 = real_1(prec2);
676 840 : for (i = 1; i < l; i++)
677 637 : if (i != _0 && i != _1)
678 231 : hash_insert(H, mkvecsmall(i), filllg1(ibin1, r1, gel(X,i), N, prec2));
679 203 : r = fillrec(H, Evec, _0, _1, X, pab, r1, N);
680 203 : if (DEBUGLEVEL) err_printf("polylogmult: k = %ld, %ld nodes\n", k, H->nb);
681 203 : 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 + EXTRAPREC64);
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 + EXTRAPREC64);
697 70 : bit += maxss(log2zeta_bound(a), 64);
698 70 : return zetamult_Akhilesh(atoe(a), bit, prec);
699 : }
700 : GEN
701 196 : zetamult(GEN s, long prec)
702 : {
703 196 : pari_sp av0 = avma, av;
704 196 : GEN z, avec, r = cgetr(prec);
705 :
706 196 : av = avma; avec = zetamultconvert_i(s,1);
707 147 : 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 252 : zetamult_interpolate(GEN s, GEN t, long prec)
714 : {
715 252 : pari_sp av = avma, av2;
716 : long i, k, l, la, bit;
717 : GEN avec, v, V;
718 :
719 252 : if (lg(s) == 1) return gen_1;
720 238 : if (!t) return zetamult(s, prec);
721 42 : avec = zetamultconvert_i(s, 1); k = zv_sum(avec);
722 42 : bit = prec2nbits(prec);
723 42 : if (bit <= 128 || k > 20 || (bit >> k) < 4)
724 42 : return zetamult_interpolate2_i(vecsmall_reverse(avec), t, prec);
725 0 : v = allstar(avec); l = lg(v); la = lg(avec);
726 0 : V = cgetg(la, t_VEC);
727 0 : for (i = 1; i < la; i++) gel(V,i) = utor(0, prec + EXTRAPREC64);
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 : GEN
739 133 : polylogmult_interpolate(GEN a, GEN z, GEN t, long prec)
740 : {
741 133 : pari_sp av = avma;
742 : GEN V, avec, A, AZ, Z;
743 : long i, la, l;
744 133 : switch(typ(a))
745 : {
746 119 : case t_VEC:
747 119 : case t_COL: a = gtovecsmall(a); break;
748 7 : case t_VECSMALL: break;
749 7 : default: pari_err_TYPE("polylogmult", a);
750 : return NULL;/*LCOV_EXCL_LINE*/
751 : }
752 126 : if (!z) return zetamult_interpolate(a, t, prec);
753 91 : switch (typ(z))
754 : {
755 0 : case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
756 0 : z = mkvec(z); break;
757 91 : case t_VEC: case t_COL: break;
758 0 : case t_VECSMALL: z = zv_to_ZV(z); break;
759 0 : default: pari_err_TYPE("polylogmult [z]", z);
760 : }
761 91 : if (lg(z) != lg(a))
762 7 : pari_err_TYPE("polylogmult [#s != #z]", mkvec2(a,z));
763 84 : if (!t) return gerepilecopy(av, zetamultevec(aztoe(a,z,prec), prec));
764 28 : avec = zetamultconvert_i(a, 1); la = lg(avec);
765 28 : AZ = allstar2(avec, z);
766 28 : A = gel(AZ, 1); l = lg(A);
767 28 : Z = gel(AZ, 2); V = zerovec(la-1);
768 182 : for (i = 1; i < l; i++)
769 : {
770 154 : pari_sp av2 = avma;
771 154 : GEN ai = gel(A,i), e = aztoe(ai, gel(Z,i), prec);
772 154 : long n = lg(ai)-1; /* > 0 */
773 154 : gel(V,n) = gerepileupto(av2, gadd(gel(V,n), zetamultevec(e, prec)));
774 : }
775 28 : return gerepileupto(av, poleval(vecreverse(V),t));
776 : }
777 :
778 : GEN
779 0 : polylogmult(GEN a, GEN z, long prec)
780 : {
781 0 : return polylogmult_interpolate(a, z, NULL, prec);
782 : }
783 :
784 : /**************************************************************/
785 : /* ALL MZV's */
786 : /**************************************************************/
787 :
788 : /* Given admissible evec w = 0e_2....e_{k-1}1, compute a,b,v such that
789 : * w=0{1}_{b-1}v{0}_{a-1}1 with v empty or admissible.
790 : * Input: binary vector evec */
791 : static void
792 945 : findabv(GEN w, long *pa, long *pb, long *pminit, long *pmmid, long *pmfin)
793 : {
794 945 : long le = lg(w) - 2;
795 945 : if (le == 0)
796 : {
797 105 : *pa = 1;
798 105 : *pb = 1;
799 105 : *pminit = 2;
800 105 : *pmfin = 2;
801 105 : *pmmid = 1;
802 : }
803 : else
804 : {
805 : long a, b, j, lv;
806 1050 : for (j = 1; j <= le; j++)
807 1050 : if (!w[j+1]) break;
808 840 : *pb = b = j;
809 1890 : for (j = le; j >= 1; j--)
810 1575 : if (w[j+1]) break;
811 840 : *pa = a = le + 1 - j;
812 840 : lv = le + 2 - a - b;
813 840 : if (lv > 0)
814 : {
815 315 : long v = fd(w, b + 1, le - a + 2), u = v + (1 << (lv-1));
816 315 : *pminit = (((1 << b) - 1) << (lv - 1)) + (v/2) + 2;
817 315 : *pmfin = (u << (a - 1)) + 2;
818 315 : *pmmid = (u >> 1) + 2;
819 : }
820 : else
821 : {
822 525 : *pminit = (1 << (b - 1)) + 1;
823 525 : *pmfin = (a == 1) ? 2 : (1 << (a - 2)) + 2;
824 525 : *pmmid = 1;
825 : }
826 : }
827 945 : }
828 :
829 : /* Returns L:
830 : * L[1] contains zeta(emptyset)_{n-1,n-1},
831 : * L[2] contains zeta({0})_{n-1,n-1}=zeta({1})_{n-1,n-1} for n >= 2,
832 : * L[m+2][n] : 1 <= m < 2^{k-2}, 1 <= n <= N + 1
833 : * contains zeta(w)_{n-1,n-1}, w corresponding to m,n
834 : * L[m+2] : 2^{k-2} <= m < 2^{k-1} contains zeta(w), w corresponding to m
835 : (code: w=0y1 iff m=1y). */
836 : static GEN
837 105 : fillL(long k, long bitprec)
838 : {
839 105 : long N = 1 + bitprec/2, prec = nbits2prec(bitprec);
840 105 : long s, j, n, m, K = 1 << (k - 1), K2 = K/2;
841 105 : GEN p1, p2, pab = get_pab(N, k), L = cgetg(K + 2, t_VEC);
842 :
843 105 : get_ibin(&gel(L,1), &gel(L,2), N, prec);
844 840 : for (m = 1; m < K2; m++)
845 : {
846 735 : gel(L, m+2) = p1 = cgetg(N+1, t_VEC);
847 83055 : for (n = 1; n < N; n++) gel(p1, n) = cgetr(prec);
848 735 : gel(p1, n) = gen_0;
849 : }
850 945 : for (m = K2; m < K; m++) gel(L, m+2) = utor(0, prec);
851 525 : for (s = 2; s <= k; s++)
852 : { /* Assume length evec < s filled */
853 : /* If evec = 0e_2...e_{s-1}1 then m = (1e_2...e_{s-1})_2 */
854 420 : GEN w = cgetg(s, t_VECSMALL);
855 420 : long M = 1 << (s - 2);
856 420 : pari_sp av = avma;
857 1995 : for (m = M; m < 2*M; m++)
858 : {
859 : GEN pinit, pfin, pmid;
860 : long comp, a, b, mbar, minit, mfin, mmid, mc;
861 1575 : p1 = gel(L, m + 2);
862 5145 : for (j = s - 1, mc = m, mbar = 1; j >= 2; j--, mc >>= 1)
863 : {
864 3570 : w[j] = mc & 1;
865 3570 : mbar = (1 - w[j]) | (mbar << 1);
866 : }
867 : /* m, mbar are dual; handle smallest, copy the other */
868 1575 : comp = mbar - m; if (comp < 0) continue; /* m > mbar */
869 945 : if (comp)
870 : {
871 630 : p2 = gel(L, mbar + 2);
872 630 : setisclone(p2); /* flag as dual */
873 : }
874 : else
875 315 : p2 = NULL; /* no copy needed if m = mbar */
876 945 : findabv(w, &a,&b,&minit,&mmid,&mfin);
877 945 : pinit= gel(L, minit);
878 945 : pfin = gel(L, mfin);
879 945 : pmid = gel(L, mmid);
880 105840 : for (n = N-1; n > 1; n--, set_avma(av))
881 : {
882 104895 : GEN t = mpmul(gel(pinit,n+1), gmael(pab, n, b));
883 104895 : GEN u = mpmul(gel(pfin, n+1), gmael(pab, n, a));
884 104895 : GEN v = gel(pmid, n+1), S = s < k ? gel(p1, n+1): p1;
885 104895 : S = mpadd(S, mpdiv(mpadd(mpadd(t, u), v), gmael(pab, n, a+b)));
886 104895 : mpaff(S, s < k ? gel(p1, n) : p1);
887 104895 : if (p2 && s < k) mpaff(S, gel(p2, n));
888 : }
889 : { /* n = 1: same formula simplifies */
890 945 : GEN t = gel(pinit,2), u = gel(pfin,2), v = gel(pmid,2);
891 945 : GEN S = s < k ? gel(p1,2): p1;
892 945 : S = mpadd(S, mpadd(mpadd(t, u), v));
893 945 : mpaff(S, s < k ? gel(p1,1) : p1);
894 945 : if (p2 && s < k) mpaff(S, gel(p2, 1));
895 945 : set_avma(av);
896 : }
897 945 : if (p2 && s == k) mpaff(p1, p2);
898 : }
899 : }
900 105 : return L;
901 : }
902 :
903 : /* bit 1 of flag unset: full, otherwise up to duality (~ half)
904 : * bit 2 of flag unset: all <= k, otherwise only k
905 : * half: 2^(k-3)+ delta_{k even} * 2^(k/2-2), sum = 2^(k-2)+2^(floor(k/2)-1)-1
906 : * full: 2^(k-2); sum = 2^(k-1)-1 */
907 : static GEN
908 105 : zetamultall_i(long k, long flag, long prec)
909 : {
910 105 : GEN res, ind, L = fillL(k, prec2nbits(prec) + 32);
911 105 : long m, K2 = 1 << (k-2), n = lg(L) - 1, m0 = (flag & 4L) ? K2 : 1;
912 :
913 105 : if (!(flag & 2L))
914 : {
915 77 : res = cgetg(n - m0, t_VEC);
916 77 : ind = cgetg(n - m0, t_VECSMALL);
917 938 : for (m = m0; m < n - 1; m++)
918 : {
919 861 : gel(res, m - m0 + 1) = m < K2 ? gmael(L, m + 2, 1) : gel(L, m + 2);
920 861 : ind[m - m0 + 1] = m;
921 : }
922 : }
923 : else
924 : { /* up to duality */
925 : long nres, c;
926 28 : if (k == 2) nres = 1;
927 28 : else if (!(flag & 2L))
928 0 : nres = (1 << (k - 2)) + (1 << ((k/2) - 1)) - 1;
929 : else
930 28 : nres = (1 << (k - 1));
931 28 : res = cgetg(nres + 1, t_VEC);
932 28 : ind = cgetg(nres + 1, t_VECSMALL);
933 350 : for (m = m0, c = 1; m < n - 1; m++)
934 : {
935 322 : GEN z = gel(L,m+2);
936 322 : if (isclone(z)) continue; /* dual */
937 182 : if (m < K2) z = gel(z,1);
938 182 : gel(res, c) = z;
939 182 : ind[c] = m; c++;
940 : }
941 28 : setlg(res, c);
942 28 : setlg(ind, c);
943 : }
944 105 : return mkvec2(res, ind);
945 : }
946 :
947 : /* fd(e, 2, lg(e)-2), e = atoe(avec) */
948 : static long
949 1876 : atom(GEN avec)
950 : {
951 1876 : long i, m, l = lg(avec);
952 1876 : if (l < 3) return 0;
953 1232 : m = 1; /* avec[1] != 0 */
954 1708 : for (i = 2; i < l-1; i++) m = (m << avec[i]) + 1;
955 1232 : return m << (avec[i]-1);
956 : }
957 : static long
958 1876 : atoind(GEN avec, long flag)
959 1876 : { return atom(avec) + (flag? 1: (1 << (zv_sum(avec) - 2))); }
960 : /* If flag is unset, L has all k1 <= k, otherwise only k */
961 : static GEN
962 644 : zetamultstar_i(GEN L, GEN avec, long flag)
963 : {
964 644 : GEN s = allstar(avec), S = gen_0;
965 644 : long i, l = lg(s);
966 2520 : for (i = 1; i < l; i++) S = gadd(S, gel(L, atoind(gel(s,i), flag)));
967 644 : return S;
968 : }
969 :
970 : /* bit 0: notstar/star
971 : * bit 1: full/half (ignored if star, always full)
972 : * bit 2: all <= k / only k
973 : * bit 3: without / with index */
974 : GEN
975 126 : zetamultall(long k, long flag, long prec)
976 : {
977 126 : pari_sp av = avma;
978 : GEN Lind, L, res;
979 : long K, k1, ct, fl;
980 :
981 126 : if (flag < 0 || flag > 15) pari_err_FLAG("zetamultall");
982 126 : if (k < 1) pari_err_DOMAIN("zetamultall", "k", "<", gen_1, stoi(k));
983 119 : if (k == 1) return cgetg(1, t_VEC);
984 112 : if (k >= 64) pari_err_OVERFLOW("zetamultall");
985 105 : if (!(flag & 1L))
986 : { /* not star */
987 49 : Lind = zetamultall_i(k, flag, prec);
988 49 : res = (flag & 8L)? Lind : gel(Lind, 1);
989 49 : return gerepilecopy(av, res);
990 : }
991 : /* star */
992 56 : fl = flag & 4L; /* 4 if k, else 0 (all <= k) */
993 56 : Lind = gerepilecopy(av, zetamultall_i(k, fl, prec)); /* full */
994 56 : L = gel(Lind, 1);
995 56 : K = 1 << (k - 2);
996 56 : res = cgetg(fl? K+1: 2*K, t_VEC);
997 196 : for (ct = 1, k1 = fl? k: 2; k1 <= k; k1++)
998 : {
999 140 : GEN w = cgetg(k1 + 1, t_VECSMALL);
1000 140 : long M = 1 << (k1 - 1), m;
1001 784 : for (m = 1; m <= M; m += 2)
1002 : {
1003 644 : pari_sp av = avma;
1004 644 : long j, mc = m;
1005 3556 : for (j = k1; j >= 1; j--) { w[j] = mc & 1; mc >>= 1; }
1006 644 : gel(res, ct++) = gerepileupto(av, zetamultstar_i(L, etoa(w), fl));
1007 : }
1008 : }
1009 56 : if (flag & 8L) res = mkvec2(res, gel(Lind, 2));
1010 56 : return gerepilecopy(av, res);
1011 : }
|
