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.