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. >>

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]


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


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


Is there any way to calculate products like $C$ to high precision in Mathematica? Thanks.

Zhenhua Liu

Posted 2015-08-15T11:25:27.107

Reputation: 153


Page 11 of http://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]: https://oeis.org/A065463

– Patrick Stevens – 2015-08-15T12:22:32.677

1Did not have any luck with the inverse symbolic calculator though. – bobbym – 2015-08-15T15:34:40.027

Thanks, @PatrickStevens! @bbgodfrey, – Zhenhua Liu – 2015-08-16T00:32:14.303

@bbgodfrey, I'll try my best to help, though I'm still a novice in using Mathematica. – Zhenhua Liu – 2015-08-16T00:41:15.150

@bobbym, Thanks. It seems no analytic expression is available. – Zhenhua Liu – 2015-08-16T00:41:21.707



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]]

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]]

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]]

J. M.'s ennui

Posted 2015-08-15T11:25:27.107

Reputation: 115 520

Thanks a lot! That's great help to me. – Zhenhua Liu – 2015-08-16T00:42:27.557

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. – yarchik – 2020-07-06T17:54:33.763

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. – Vaclav Kotesovec – 2020-07-08T20:20:40.030


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}]


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)));
 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}]


Vaclav Kotesovec

Posted 2015-08-15T11:25:27.107

Reputation: 2 114

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? – yarchik – 2020-07-06T17:57:36.677

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. – Vaclav Kotesovec – 2020-07-06T18:56:24.910

For application of Riemann zeta prime function see (for example) http://numbers.computation.free.fr/Constants/Miscellaneous/constantsNumTheory.html

– Vaclav Kotesovec – 2020-07-06T19:05:54.923

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. – yarchik – 2020-07-06T20:11:00.777

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. – J. M.'s ennui – 2020-07-09T04:06:58.957