High precision calculation of infinite product involving prime numbers

Question

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.

Answer 1

Using the formula given in the arXiv preprint Patrick linked to for the "carefree constant" gives:

Exp[NSum[(-1)^k PrimeZetaP[k] (1 - LucasL[k])/k, {k, 2, ∞}, Compiled -> False,
         Method -> "AlternatingSigns", NSumTerms -> 20, WorkingPrecision -> 30]]
   0.704442200999165592738713909247

Note that this agrees with the result in the OEIS up to twenty digits. The only speed-limiting part of this is the calculation of the prime zeta function.

---

In fact, Vaclav's answer can also be modified, so that one can exploit the convergence acceleration capabilities of NSum[].

The key identity to use for this (whose proof I leave as an exercise) is that

$$\sum_{n=2}^\infty \frac{(-1)^n \left(1-L_n\right)}{n p^n}=\log\left(1-\frac1{p(1+p)}\right)$$

where $L_n$ is a Lucas number. (This should also give a hint on how the Lucas number showed up in my previous solution.)

With that,

5/6 Exp[NSum[((-1)^n (1 - LucasL[n]) (PrimeZetaP[n] - 2^-n))/n, {n, 2, ∞},
             Method -> "AlternatingSigns", NSumTerms -> 55, WorkingPrecision -> 125]]
   0.704442200999165592736603350326637210188586431417098049414226842591097056682006778536808244145693135370271359151436811784885404

Of course, we can take out more terms the same way as in Vaclav's answer:

Product[1 - 1/(p (p + 1)), {p, Prime[Range[3]]}]
Exp[NSum[((-1)^n (1 - LucasL[n]))/n (PrimeZetaP[n] -
                                     Evaluate[Sum[p^-n, {p, Prime[Range[3]]}]]),
         {n, 2, ∞}, Method -> "AlternatingSigns",
         NSumTerms -> 45, WorkingPrecision -> 135]]
   0.7044422009991655927366033503266372101885864314170980494142268425910970566820067785368082441456931337676420607204592721529533500243226539

Answer 2

This is my program that is also based on PrimeZetaP, but is much faster:

 $MaxExtraPrecision = 1000; Clear[f]; f[p_] := (1 - 1/(p*(p + 1))); 
 Do[c = Rest[CoefficientList[Series[Log[f[1/x]], {x, 0, m}], x]]; 
 Print[f[2] * Exp[N[Sum[Indexed[c, n]*(PrimeZetaP[n] - 1/2^n), {n, 2, m}], 100]]], {m, 100, 1000, 100}]

 (*
 0.7044422009991655927366033503228685017899862452605546001261786165281699328975676993338639277737095544
 0.7044422009991655927366033503266372101885864314170980494113214008058012515341318822970728422740691116
 0.7044422009991655927366033503266372101885864314170980494142268425910970566820067785338281514179409974
 0.7044422009991655927366033503266372101885864314170980494142268425910970566820067785368082441456931338
 0.7044422009991655927366033503266372101885864314170980494142268425910970566820067785368082441456931338
 0.7044422009991655927366033503266372101885864314170980494142268425910970566820067785368082441456931338
 *)

A more detailed description of the method

$$\prod_{p}f(p) = \prod_{p}\left(1-\frac{1}{p(p+1)}\right) = \exp\left( \sum_{p}\log\left(1-\frac{1}{p(p+1)}\right)\right)$$

We expand the function in a power series:

 Normal[Series[Log[f[1/x]], {x, 0, 10}]] /. x -> 1/p // TraditionalForm

$$\log\left(1-\frac{1}{p(p+1)}\right)=-\frac{1}{p^2}+\frac{1}{p^3}-\frac{3}{2 p^4}+\frac{2}{p^5}-\frac{17}{6 p^6}+\frac{4}{p^7}-\frac{23}{4 p^8}+...$$

$$\sum_{p}\log\left(1-\frac{1}{p(p+1)}\right)=\sum _{n=2}^{\infty } c_n P(n)$$

where $P(n)$ is the PrimeZetaP function.

The rate of convergence can be even faster if we explicitly compute the first few terms of the product (in the program above for $p = 2$). The following is an example where we separately calculate the terms corresponding to $p = 2,3,5$:

 $MaxExtraPrecision=1000; Clear[f]; f[p_]:=(1-1/(p*(p+1)));
 Do[c=Rest[CoefficientList[Series[Log[f[1/x]],{x,0,m}],x]];
 Print[f[2]*f[3]*f[5]*Exp[N[Sum[Indexed[c,n]*(PrimeZetaP[n] - 1/2^n - 1/3^n - 1/5^n), {n,2,m}],100]]], {m,100,500,100}]

 (*
 0.70444220099916559273660335032663721018858643141709804941422684259077579047803922512750258917649590736
 0.70444220099916559273660335032663721018858643141709804941422684259109705668200677853680824414569313377
 0.70444220099916559273660335032663721018858643141709804941422684259109705668200677853680824414569313377
 0.70444220099916559273660335032663721018858643141709804941422684259109705668200677853680824414569313377
 0.70444220099916559273660335032663721018858643141709804941422684259109705668200677853680824414569313377
 *)