This is easy to do with PARI/GP. Here is my code
```
default(realprecision,38);
p(n) = binomial(n+2,2);
Y(k,b) = sumpos(l=0, (2*l+3)/(p(k)+p(l))^b);
Z(a,b) = sumpos(k=0, (2*k+3)/p(k)^a*Y(k,b));
print(2*Z(2,2)+4*Z(1,3))
/* 15.999999999999999999999999999893290870 */
default(realprecision,50);
print(2*Z(2,2)+4*Z(1,3))
/* 15.999999999999999999999999999999999999999987875982 */
```
It takes a few seconds for 38 decimal places and
almost a minute for 50 decimal places. You
can read the documentation for what numerical
methods are being used.