This is easy to do with [PARI/GP](https://pari.math.u-bordeaux.fr/).
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](https://pari.math.u-bordeaux.fr/doc.html)
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.

More details on Monien summation are in the paper
[Gaussian Summation: An Exponentially Converging Summation Scheme](https://arxiv.org/abs/math/0611057)
by Hartmut Monien.