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).