# 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 Aug 15, 2015 at 12:22 • Did not have any luck with the inverse symbolic calculator though. Aug 15, 2015 at 15:34 • Thanks, @PatrickStevens! @bbgodfrey, – user32407 Aug 16, 2015 at 0:32 • @bbgodfrey, I'll try my best to help, though I'm still a novice in using Mathematica. – user32407 Aug 16, 2015 at 0:41 • @bobbym, Thanks. It seems no analytic expression is available. – user32407 Aug 16, 2015 at 0:41 ## 2 Answers 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  • Thanks a lot! That's great help to me. – user32407 Aug 16, 2015 at 0:42 • Can you please explain the connection between the original and your expression. Which formulas from the preprint are you using? I think it is better to type the equations here because external links might disappear, even from arXiv. Jul 6, 2020 at 17:54 • Both methods are almost identical in the use of PrimeZetaP, the difference is that the formula with Lucas numbers is applicable only to this one product, while my decomposition in the power series is applicable to any convergent product over all prime numbers. Jul 8, 2020 at 20:20 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
*)

• I have the same comment to you, as to the accepted answer. Can you please explain the connection between your implementation and the original product. Also, number like 100 and 1000 look a bit arbitrary. What do they illustrate, the rate of convergence? Can you explicitly estimate the rate? Jul 6, 2020 at 17:57
• This program is general and can be used for any convergent product over all primes. Experimentally: if we want an accuracy of 100 decimal places, about 400 terms of the series are needed, for 200 digits at least 800 terms and for 300 digits at least 1100 terms are needed. Jul 6, 2020 at 18:56
• For application of Riemann zeta prime function see (for example) numbers.computation.free.fr/Constants/Miscellaneous/… Jul 6, 2020 at 19:05
• Thank you, I have see your link, but there are some formatting issues for equations. Therefore, it would be nice to have some relevant formulas in the post. Jul 6, 2020 at 20:11
• Using the Lucas number result in my answer, you could just do c = LinearRecurrence[{-2, 0, 1}, {0, -2, 3}, m]/Range[m] in both of your snippets. Jul 9, 2020 at 4:06