Accelerating convergence for some double sums I am interested in the following general double sums, for integers $a\geq 1$ and $b\geq 2$,
$$Z(a,b) = \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^a} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^b},$$
which are converging very slowly. For these sums, there is also an alternative expression as an iterated integral in dimension $a+b$, similar to multiple zeta values.
I would like more particularly to find $Z(2,2)$ and $Z(1,3)$ with precision as large as possible. Using the naive summation, I could only obtain 4 decimal digits, namely $Z(2,2) \simeq 4.7058$ and $Z(1,3) \simeq 1.6470$.
It is known that $2 Z(2,2) + 4 Z(1,3) = 16$.

What would be a smart way to accelerate the convergence, in general and in the special case using maybe the previous formula ?

 A: With Mathematica I can first sum the series over $\ell$ to get a closed-form expression in terms of polygamma functions,

$$Z(2,2)= \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^2} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^2}$$
$$\qquad=\sum_{k\geq 0}\frac{16 (2 k+3) }{(k+1)^2 (k+2)^2 \sqrt{-4 k (k+3)-7}}\left(\psi ^{(1)}\left[\frac{1}{2} \left(3-\sqrt{-4 k (k+3)-7}\right)\right)-\psi ^{(1)}\left(\frac{1}{2} \left(\sqrt{-4 k (k+3)-7}+3\right)\right)\right],$$

and then obtain a high-precision result using the Wynn epsilon method.. To fifteen decimal places I find
$$Z(2,2)=4.705905644127174\cdots$$
NSum[Z22,{k,0,Infinity},NSumTerms->2000,WorkingPrecision->250,Method->"WynnEpsilon"]

As a test, I also computed $Z(1,3)$ to fifteen decimal places, which gave
$$Z(1,3)=1.647047177936412\cdots.$$
And thus I find
$$2Z(2,2)+4Z(1,3)=15.999999999999996\cdots$$
correct to fifteen decimal places.
A: This is easy to do with PARI/GP.
Here is my code
p(n) = binomial(n+2,2);
Y(k,b) = sumnum(l=0, (2*l+3)/(p(k)+p(l))^b,sumtable);
Z(a,b) = sumnum(k=0, (2*k+3)/p(k)^a*Y(k,b));
default(realprecision,57);
sumtable = sumnuminit();
print(2*Z(2,2)+4*Z(1,3))
/* 16.0000000000000000000000000000000000000000000000000000000 */

It takes 1.361 CPU seconds for 57 decimal places. The
documentation
has some information about the methods being used and there
are several summation functions other than sumnum. For
example, I originally used sumpos but sumnum is faster.
Thanks to Henri Cohen for his comment to use sumnuminit
to speed up the calculation of Y(k,b). Thanks to Jorge
Zuniga for his comment that replacing sumnum with summonien
is much faster. And especially thanks to the developers of
PARI/GP who provided these fast numerical summation functions.
Some details are in the arXiv paper
Gaussian Summation: An Exponentially Converging Summation Scheme
by Hartmut Monien. (published in 2010 in
Mathematics of Computation).
