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);
Z(a,b) = sumnum(k=0, (2*k+3)/p(k)^a*Y(k,b));
default(realprecision,57);
print(2*Z(2,2)+4*Z(1,3))
/* 16.0000000000000000000000000000000000000000000000000000000 */
It takes four 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 better.