# High precision calculation of infinite product involving prime numbers

I'm recently studying some topics in analytic number theory and I have encountered results involving the infinite product $$C=\prod_{p}\left(1-\frac{1}{p(p+1)}\right)$$ where $p$ denotes calculating the product over all prime numbers.

Or more intuitively,$$C=\left(1-\frac{1}{2\cdot3}\right)\cdot\left(1-\frac{1}{3\cdot4}\right)\cdot\left(1-\frac{1}{5\cdot6}\right)\cdot\left(1-\frac{1}{7\cdot8}\right)\cdot\left(1-\frac{1}{11\cdot12}\right)\cdots$$ Mathematica can show that $C\approx0.704442.$

However, I need a high precision result, roughly accurate to at least 20 digits, so that I could use the Inverse Symbolic Calculator to make an educated guess of any possible analytic expression for $C$ in forms like $\dfrac{\zeta(2)\zeta(3)}{\zeta(6)}$.

When I use NProduct, I get

NProduct[1 - 1/(Prime[k] (Prime[k] + 1)), {k, 1, Infinity},
AccuracyGoal -> 20, PrecisionGoal -> 20]

During evaluation of In[14]:= Prime::intpp: Positive integer argument expected in Prime[15.]. >>
During evaluation of In[14]:= Prime::intpp: Positive integer argument expected in Prime[14.]. >>
During evaluation of In[14]:= Prime::intpp: Positive integer argument expected in Prime[13.]. >>
During evaluation of In[14]:= General::stop: Further output of Prime::intpp will be suppressed during this calculation. >>
0.705024


However, when I use N[Product[]], the calculation was so slow that I had to abort the calculation after minutes of calculating.

N[Product[1 - 1/(Prime[k] (Prime[k] + 1)), {k, 1, 100000}], 50]

0.70444223707595873775750824971600704569326374380066

N[Product[1 - 1/(Prime[k] (Prime[k] + 1)), {k, 1, 1000000}], 50]

0.70444220359808164073982578152558851839613830115938

N[Product[1 - 1/(Prime[k] (Prime[k] + 1)), {k, 1, 10000000}], 50]

$Aborted  Is there any way to calculate products like$C$to high precision in Mathematica? Thanks. • Page 11 of arxiv.org/pdf/0903.2514v2.pdf gives the expansion, as$Q^{(1)}_1\$. It's the "carefree constant", according to Wikipedia. (Thanks to [OEIS] for this.) [OEIS]: oeis.org/A065463 – Patrick Stevens Aug 15 '15 at 12:22
• Did not have any luck with the inverse symbolic calculator though. – bobbym Aug 15 '15 at 15:34
• Thanks, @PatrickStevens! @bbgodfrey, – Zhenhua Liu Aug 16 '15 at 0:32
• @bbgodfrey, I'll try my best to help, though I'm still a novice in using Mathematica. – Zhenhua Liu Aug 16 '15 at 0:41
• @bobbym, Thanks. It seems no analytic expression is available. – Zhenhua Liu Aug 16 '15 at 0:41

Exp[NSum[(-1)^k PrimeZetaP[k] (1 - LucasL[k])/k, {k, 2, ∞}, Compiled -> False,