Numerical Evaluation of MZVs

DISCLAIMER: Check the evaluations using other software, before trusting them. No guarantee that the code is correct or free of silly mistakes. Gp/pari implements both zeta and interpolated zetas, as in order to verify answers.

Mathematica implementation from the talk by Steven Charlton in Computing Multiple Zeta Seminar. Method originally given by Don Zagier in gp/pari. $MinPrecision = 100; zetaM[M_, trorder_, s_List] := Module[{asym, i, d, vec}, (     asym = Table[0, {Length[s] + 1}, {trorder + 1}];      asym1, 1 = 1.0`1;      For[d = 1, d <= Length[s], d++,       For[i = 1, i <= trorder, i ++,        asymd + 1, i + 1 =          1/i (If[i + 1   < sd, 0, asym d, i - s[[d + 2]]] - Sum[asymd + 1, k + 1 Binomial[i, i + 1 - k], {k, 0, i - 1}])        ]];      vec = Sum[asymAll, i + 1 / M^(i), {i, 0, trorder}];      For[i = M, i >= 1, i--,       For[r = Length[s] + 1, r >= 2, r--,         vecr = vecr + vecr - 1*1/i^sr - 1;         ]; ];      vec-1    )];

Code to evaluate interpolated \( \zeta^t(s_1,\ldots,s_r) \). Main idea is that for \[ \zeta^t_{\gg M}(s_1,\ldots,s_r) = \sum_{i_1 \geq \cdots \geq i_r > M} \frac{t^{\delta_{i_1 = i_2} + \cdots + \delta_{i_{r-1} = i_r}}}{i_1^{s_1} \cdots i_r^{s_r}} \,, \] we have the recursion \[ \zeta^t_{\gg M-1}(s_1,\ldots,s_r) = \zeta^t_{\gg M}(s_1,\ldots,s_r) + \frac{1}{M^{s_r}}   \zeta^t_{\gg M}(s_1,\ldots,s_{r-1}) + \frac{t}{M^{s_r + s_{r-1}}}  \zeta^t_{\gg M}(s_1,\ldots,s_{r-2}) + \cdots + \frac{t^{r-1}}{M^{s_r + \cdots + s_1}}  \zeta^t_{\gg M}(\emptyset) \] Then use the same idea as in the talk to find the asymptotic series, and use the recursion to explicitly sum the head of the series. $MinPrecision = 100; zetaInterpM[M_, trorder_, s_List] := Module[{ k, r, i}, (   asym = Table[0, {Length[s] + 1}, {trorder + 1}];    asym1, 1 = 1.0`1;    For[d = 1, d <= Length[s], d++,     For[i = 1, i <= trorder, i ++,      asymd + 1, i + 1 = 1/i ( Sum[If[i + 1 - Total[sd - k ;; d] + 1 < 1, 0, t^(k) asym (d - 1) - k + 1,                i + 1 - Total[s[[d - k ;; d] + 1]]], {k, 0, d - 1}] -            Sum[asymd + 1, k + 1 Binomial[i, i + 1 - k], {k, 0, i - 1}]) // Expand;     ]     ];    vec = Sum[asymAll, i + 1 /M^(i), {i, 0, trorder}] // Expand;    For[i = M, i >= 1, i--,     For[r = Length[s] + 1, r >= 2, r--,       vecr = vecr + (Sum[ t^(r - 1 - k )/i^(Total[sk ;; r - 1])* veck - 1 + 1, {k, 1, r - 1} ] // Expand)      ];      ];    vec-1    )];