Efficiently summing reciprocal of polynomial over primes. For fun, I have been trying to calculate the sum of the reciprocals of the cube-full numbers. I have managed to show that the limit is equal to $$\frac{\zeta(3)\zeta(4)\zeta(5)}{\zeta(8)\zeta(10)}\prod_{p\ \mbox{prime}}\left(1-\frac1{(p^4+1)(p^5+1)}\right)\approx1.3397841535.$$ This product converges pretty fast, because of the order $9$ polynomial $f(x)=(x^4+1)(x^5+1)$ in the numerator. By simply taking the primes up to $10^8$, I already got $64$ digits of precision. $$1.3397841535743472465991525865148860527752422497881828066630150676$$ However, this method requires exponential time to calculate more digits. I was wondering if there is a faster, or even a polynomial time algorithm to calculate more digits.
One thing I tried is to take the logarithm of the product. $$\log\left(\prod_{p\ \mbox{prime}}\left(1-\frac1{f(p)}\right)\right)=\sum_{p\ \mbox{prime}}\log\left(1-\frac1{f(p)}\right)$$ By taking the Taylor series of the natural logarithm, we get $$\log\left(1-\frac1{f(p)}\right)=\sum_{n=1}^\infty\frac{-1}{n(f(p))^n}.$$ By absolute convergence, we can interchange the sums to obtain $$\sum_{n=1}^\infty\frac{-1}n\sum_{p\ \mbox{prime}}\frac1{(f(p))^n}.$$ For all $n$, of course $(f(p))^n$ is a polynomial, so the question becomes how we can efficiently sum the reciprocal of a polynomial over the primes. Is there some sort of analog for the Euler-Maclaurin formula for primes?
 A: $$\log(1-\frac{x^{9}}{(1+x^4)(1+x^5)}) = \sum_{m\ge 0} a_m x^{9+m}$$
$$\log\prod_p (1-\frac{p^{-9}}{(1+p^{-4})(1+p^{-5})}) =\sum_p \sum_m a_m p^{-9-m}=\sum_m a_m \sum_p p^{-9-m}$$ $$=\sum_m a_m \sum_l \frac{\mu(l)}{l} \log\zeta(l(9+m)) = \sum_{n\le N}(\sum_{l|n}\frac{ a_{n/l-9} \mu(l)}{l} )\log\zeta(n)+O(2^{-N})$$
A: You can use the Taylor expansion of $ 1 - \frac{x^9}{(1+ x^4)(1+x^5)}$ to pull out as many zeta factors as you need
\begin{eqnarray}
&& \prod_{p} \left(1 - \frac{1}{(p^4 + 1)(p^5 + 1)} \right) \\
 &=& \prod_{p} \left( 1 - \frac{1}{p^9} + \frac{1}{p^{13}} + \frac{1}{p^{14}} - \frac{1}{p^{17}} - \frac{1}{p^{18}} + O(\frac{1}{p^{{19}}}) \right) \\
&=& \prod_{p} \left(1 - \frac{1}{p^9}\right) \left(1 + \frac{1}{p^{13}} + \frac{1}{p^{14}} - \frac{1}{p^{17}} - \frac{1}{p^{18}} + O(\frac{1}{p^{{19}}})\right) \\
&& \dots \\
&=& \prod_{p} \left(1 - \frac{1}{p^9}\right)\left(1 + \frac{1}{p^{13}}\right)\left(1 + \frac{1}{p^{14}}\right)\left(1 - \frac{1}{p^{17}}\right)\left(1 - \frac{1}{p^{18}}\right) \left( 1 + O(\frac{1}{p^{{19}}})\right) \\
&& \frac{\zeta(13)\zeta(14)}{\zeta(9)\zeta(26)\zeta(28)\zeta(17)\zeta(18)}\prod_{p}  \left(1 + O(\frac{1}{p^{{19}}})\right) \\
\end{eqnarray}
In general you can do this for euler products of the form $\prod_p f(p)$ where $f$ has a convergent expansion in $p^{-1}$. The convergence of the above formula can be slow, but a solution is to compute the product explicitly for primes less than some $N$, and then use the formula for the remainder, being sure to adjust the zeta functions for the missing factors.
