Product over the primes 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
 A: 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$
A: 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.
