#!/usr/bin/env python3
"""tp_bound.py: a rigorous enclosure of the twin-prime constant
    TP = prod_{p>2} (1 - 1/(p-1)^2)
from the partial product over 2 < p <= P (mpmath.iv, outward rounding) and the tail bound
    1 >= prod_{p>P} (1 - 1/(p-1)^2) >= 1 - sum_{p>P} 1/(p-1)^2 >= 1 - sum_{n>=P} 1/n^2 > 1 - 1/(P-1),
(p > P prime => p - 1 >= P, distinct p give distinct p - 1; prod(1-x_i) >= 1 - sum x_i for x_i in [0,1]).
Output: tp_bound.out (lower and upper ends). Positive control: the published value 0.66016181584686957392... must lie inside."""
import sys, numpy as np
from mpmath import iv, mp
iv.dps = 30; mp.dps = 30
P = int(sys.argv[1]) if len(sys.argv) > 1 else 10**7
s = np.ones(P + 1, dtype=bool); s[:2] = False
for i in range(2, int(P ** 0.5) + 1):
    if s[i]: s[i * i::i] = False
primes = np.nonzero(s)[0]
prod = iv.mpf(1)
for p in primes:
    p = int(p)
    if p == 2: continue
    prod *= 1 - iv.mpf(1) / (p - 1) ** 2
lower = prod * (1 - iv.mpf(1) / (P - 1))
TP_lo = mp.make_mpf(lower._mpi_[0]); TP_hi = mp.make_mpf(prod._mpi_[1])
known = mp.mpf('0.66016181584686957392781211001455577843262336028473341331944')
with open('tp_bound.out', 'w') as f:
    f.write(f'P = {P}, odd primes used = {len(primes) - 1}\n')
    f.write(f'TP in [{mp.nstr(TP_lo, 20)}, {mp.nstr(TP_hi, 20)}]\n')
    f.write(f'TP_lo_str {mp.nstr(TP_lo, 25)}\n')
    f.write(f'positive control: published 0.66016181584686957392... inside: {TP_lo <= known <= TP_hi}\n')
print(open('tp_bound.out').read())