I'm trying to estimate the product $$\prod_{p\lt q\lt r\lt s}1-\frac{24}{(pqrs)^2}$$ where $p,q,r,s$ are primes.

This is for the purpose of calculating the density of Sloane's A070284 [1]. The idea is that there are 4! congruence classes mod $(pqrs)^2$ such that $n,n+1,n+2,n+3$ are each congruent to 0 mod the square of one of the primes. A number is not the first of four such numbers exactly when it is not in any of these congruence classes for any quadruple of distinct primes.

I feel that there should be some way to accelerate the calculation, possibly using the prime zeta function. At the least, this is useful for the sum approximation which is far easier to calculate with dynamic programming (making it essentially linear rather than quartic, if you can spare $O(n)$ memory). Unfortunately there are enough congruence classes that the cancellation is an important part of the problem, so the sum approximation is poor.

Any suggestions for the calculation, or different approaches to the problem, would be appreciated.

[1] http://oeis.org/classic/A070284

  • $\begingroup$ For what it's worth, A070284(100000) = 46541275, which suggests that the product is around 0.9978. $$\prod_{p\lt q\lt r\lt s\lt1000}1-\frac{24}{(pqrs)^2}\approx0.9975$$ which is perhaps close? $\endgroup$
    – Charles
    Jul 23, 2010 at 5:34
  • $\begingroup$ As far as I can see from, for example, [W. Butske, L.M. Jaje and D.R. Mayernik, Math. Comp. 69 (2000) 407--420], it's not a big deal to compute sums/products over primes up to $10^7$---$10^9$. You don't say how close you wish to be the constant... Do you wish to prove it is 1?! :-) $\endgroup$ Jul 23, 2010 at 8:12
  • 2
    $\begingroup$ There are 664579 primes less than $10^7$, so calculating the formula to $10^7$ would take ${664579\choose4}\approx8\times10^{21}$ calculations of the inner loop. Implemented efficiently the inner loop would take four multiprecision operations (addition, subtraction, squaring, and division) for a total of perhaps a thousand cycles. $8\times10^{24}$ cycles at 3.0 GHz would take ~ 85,000 core-years. $\endgroup$
    – Charles
    Jul 23, 2010 at 14:32
  • $\begingroup$ I calculate the product up to 7000 as 0.9974814. $\endgroup$
    – Charles
    Jul 23, 2010 at 21:34
  • $\begingroup$ btw, do you know lower bound? $\endgroup$
    – joro
    Apr 14, 2013 at 8:27

2 Answers 2


Modulo errors got unconditionally slightly smaller upper bound than your conjecture. prodeulerrat from Cohen's pari script can compute $\prod_{\text{p prime}}{1-1/O(x^2)}$.

For fixed $p < q < r < L$ define $g(p,q,r)=\prod_{\text{s prime}}{1-24/(pqrs)^2}$

$g$ is efficiently computable by the script and contains the contribution of $s > L$ times extra factors $ < L$. Compute the extra factors the slow way for each triple to get the exact product for $p < q < r < L < s$ for fixed $p,q,r$. Upper bound (since all terms are $<1$) for your product is the product over the exact products times the product $< L$. .

For $L=200$ got $0.9974805703680839724432166987$


Rewrite the product as $$\prod_{p\lt q\lt r}f\left(r,\frac{24}{(pqr)^2}\right)$$ with $$f(r,\alpha)=\prod_{s\gt r}1-\frac{\alpha}{s^2}$$ where $p,q,r,s$ are always taken over the primes.

Now $f$ can be approximated by $$f(r,\alpha)\approx\left(\prod_{s=N(r)}^{L}1-\frac{\alpha}{s^2}\right)-\frac{1}{\alpha}P_{s\ge L}(2)$$ where $N(x)$ is the least prime greater than $x$ and $$P_{\ge L}(z)=\sum_{s\ge L}\frac{1}{s^2}=P(z)-\sum_{s\lt L}\frac{1}{s^2}=\sum_{k=1}^\infty\frac{\mu(k)}{k}\log\zeta(kz)-\sum_{s\lt L}\frac{1}{s^2}.$$

Let $F=\prod_{s=N(r)}^{L}1-\alpha/s^2$ denote the finite part and $I=P_{s\ge L}(2)/\alpha$ denote the infinite part (so the above is $f\approx F-I$). Then a better approximation is $$f(r,\alpha)\approx Fe^{-I}$$ which takes into account the cancellation amongst the terms.

Now a table of values of $P_{\ge L}(2)$ can be constructed readily in quasilinear time, allowing a good approximation for $$\prod_{p\lt q\lt r\lt L,r\lt s}1-\frac{24}{(pqrs)^2}$$ to be calculated for 'reasonable' values of $L$ (say, L = 10000). That is, the number of terms calculated does not improve, but the final prime is taken from $N(r)$ to $\infty$ rather than $N(r)$ to $L$, at the cost of a cubic number of table lookups, divisions, and subtractions -- essentially, 'for free'.

This, together with the smallness of the terms, may give a reasonable precision. Standard tricks (calculating the product $P$ of $1-x_i$ as a sum $x_i(1-S_i)$ with $P=1-S_\infty$, storing partial products rather than recalculating) improve the constant factor several-fold.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.