Editing Openai/68ec50da-cf00-8005-b5f6-b683506e5853 (section)

==== The idea is to multiply the factors up to a large cutoff BBB and bound the remaining tail. ====
For 0<x≤1/20<x\le 1/20<x≤1/2 we have −log⁡(1−x)≤x+2x2-\log(1-x)\le x+2x^2−log(1−x)≤x+2x2.
Since 2/p2≤2/192<0.0062/p^2 \le 2/19^2<0.0062/p2≤2/192<0.006, for primes p>17p>17p>17 this gives

∣log⁡ ⁣ ⁣∏p>B ⁣(1−2p2)∣≤2 ⁣ ⁣∑p>B ⁣1p2+8 ⁣ ⁣∑p>B ⁣1p4≤2B+83(B−1)3,\bigl|\log\!\!\prod_{p>B}\!(1-\tfrac{2}{p^2})\bigr|
\le 2\!\!\sum_{p>B}\!\frac{1}{p^2}+8\!\!\sum_{p>B}\!\frac{1}{p^4}
\le \frac{2}{B}+\frac{8}{3(B-1)^3},logp>B∏(1−p22)≤2p>B∑p21+8p>B∑p41≤B2+3(B−1)38,
(using ∑n>Bn−2≤1/B\sum_{n>B}n^{-2}\le 1/B∑n>Bn−2≤1/B and ∑n>Bn−4≤1/(3(B−1)3)\sum_{n>B}n^{-4}\le 1/\bigl(3(B-1)^3\bigr)∑n>Bn−4≤1/(3(B−1)3)).
Thus the full product lies in

[P(B) e−(2/B+8/(3(B−1)3)) , P(B)],\Bigl[P(B)\,e^{-\bigl(2/B+8/(3(B-1)^3)\bigr)}\ ,\ P(B)\Bigr],[P(B)e−(2/B+8/(3(B−1)3)) , P(B)],
where P(B)=∏17<p≤B, p≡1(4)(1−2/p2)P(B)=\prod_{17<p\le B,\ p\equiv1(4)}(1-2/p^2)P(B)=∏17<p≤B, p≡1(4)(1−2/p2) is the partial product.

<syntaxhighlight lang="python">import math

def primes_up_to(n: int) -> list[int]:
    """Simple sieve of Eratosthenes, returns all primes ≤ n."""
    if n < 2:
        return []
    sieve = bytearray(b"\x01") * (n + 1)
    sieve[:2] = b"\x00\x00"
    m = int(n**0.5)
    for i in range(2, m + 1):
        if sieve[i]:
            step = i
            start = i * i
            sieve[start:n + 1:step] = b"\x00" * ((n - start) // step + 1)
    return [i for i in range(2, n + 1) if sieve[i]]

def product_mod4eq1_excluding_17(B: int) -> tuple[float, float, float, int]:
    """
    Returns:
      (partial_product P(B), rigorous_lower_bound, rigorous_upper_bound, count_of_primes_used)
    for ∏_{p>17, p≡1 (mod 4), p≤B} (1 - 2/p^2), together with a rigorous
    lower/upper bound for the infinite product (tail bounded analytically).
    """
    primes = primes_up_to(B)
    selected = [p for p in primes if p % 4 == 1 and p > 17]
    # Sum logs for numerical stability
    s = 0.0
    for p in selected:
        s += math.log1p(-2.0 / (p * p))
    partial = math.exp(s)
    # Tail bound: |log Tail| ≤ 2/B + 8/(3(B-1)^3)
    tail_log_bound = 2.0 / B + 8.0 / (3.0 * (B - 1) ** 3) if B > 1 else float("inf")
    lower = partial * math.exp(-tail_log_bound)  # full product ≥ this
    upper = partial                               # full product ≤ this
    return partial, lower, upper, len(selected)

def product_with_target_error(tol=2e-7, start_B=1_000_000, max_B=50_000_000):
    """
    Increase B until the rigorous interval [lower,upper] has width ≤ tol.
    Returns (value_midpoint, lower, upper, B, primes_used).
    """
    B = start_B
    while True:
        partial, lower, upper, k = product_mod4eq1_excluding_17(B)
        if upper - lower <= tol or B >= max_B:
            return (0.5 * (lower + upper), lower, upper, B, k)
        B = min(max_B, B * 10)

if __name__ == "__main__":
    val, low, up, B, k = product_with_target_error(tol=2e-7)
    print(f"cutoff B = {B:,}, primes used = {k}")
    print(f"product ≈ {val:.12f}  (rigorous interval [{low:.12f}, {up:.12f}])")
    print(f"absolute error ≤ {up - low:.3e}")

</syntaxhighlight>

===== Using B=10,000,000B=10{,}000{,}000B=10,000,000 (ten million), I obtain: =====
* partial product P(B)=0.991 161 303 471 1718P(B)=0.991\,161\,303\,471\,1718P(B)=0.9911613034711718,
* rigorous interval for the infinite product [ 0.991 161 105 238 931, 0.991 161 303 471 1718 ][\,0.991\,161\,105\,238\,931,\ 0.991\,161\,303\,471\,1718\,][0.991161105238931, 0.9911613034711718],

so the true value is 0.9911612…0.9911612\ldots0.9911612… with absolute error <2×10−7<2\times10^{-7}<2×10−7.

(You can tighten the bound by increasing <code>start_B</code> or by lowering <code>tol</code>.)