#!/usr/bin/env python3
"""
Rigorous chance-coincidence enumeration v2 for m_p/m_e.

Honest accounting: the DCT formula 12·153 + 4μ² + 1/144 = 1836.152842
matches CODATA 1836.152673 at relative error 9.2 × 10⁻⁸ ≈ 92 ppb,
not 9 ppb (the Tier 4 audit reported "9 ppb" but the actual gap is
~92 ppb). We use this honest tolerance.

We expand the search:
  formula = a·b ± c·X ± 1/d ± e/f
  where X ∈ {μ², 1/φⁿ, mixed √5 terms, 1/240, 1/120, ...}
  and a,b,c,d,e,f range over the small-integer building blocks.

We also test what fraction of the search space hits the target at:
  • 92 ppb (the actual DCT match)
  • 1 ppm
  • 10 ppm
  • 100 ppm
  • 0.1%
This calibrates the chance probability honestly across precision tiers.
"""
import math

phi = (1 + math.sqrt(5))/2
target = 1836.152673
DCT = 12*153 + (7-3*math.sqrt(5))/2 + 1/144

print(f"Target (CODATA m_p/m_e): {target}")
print(f"DCT prediction:          {DCT}")
print(f"DCT − target:            {DCT - target:+.6e}")
print(f"Relative agreement:      {(DCT - target)/target:+.3e} = {(DCT-target)/target*1e9:+.1f} ppb")
print()

# Building blocks (an honest ENUMERATIVE basis)
SMALL_INTS = sorted(set([
    *range(1, 51),
    5, 8, 16, 24, 120, 600, 720, 1200, 153, 144, 240, 248,
    12, 30, 20, 31, 154, 72, 9, 25, 36, 49, 81, 100,
]))

CONT_BLOCKS = {
    "mu2":     (7 - 3*math.sqrt(5))/8,    # μ²
    "4mu2":    (7 - 3*math.sqrt(5))/2,    # = 1/φ⁴
    "phi-1":   1/phi,
    "phi-2":   1/phi**2,
    "phi-4":   1/phi**4,
    "phi-6":   1/phi**6,
    "phi-8":   1/phi**8,
    "0":       0.0,
    "sqrt5":   math.sqrt(5),
}

# Target window tiers — count hits at each
tiers = {
    "1 ppb":      1e-9 * target,
    "10 ppb":     1e-8 * target,
    "92 ppb":     9.2e-8 * target,
    "1 ppm":      1e-6 * target,
    "10 ppm":     1e-5 * target,
    "100 ppm":    1e-4 * target,
}

# =====================================================================
# Sweep
# =====================================================================
hits = {tname: 0 for tname in tiers}
hits_examples = {tname: [] for tname in tiers}
checks = 0

# Find a*b near target ± 5 (gives leading term within reach)
ab_pairs = []
for a in SMALL_INTS:
    if a == 0: continue
    for b in SMALL_INTS:
        if b == 0: continue
        if abs(a*b - target) < 5:
            ab_pairs.append((a, b))
print(f"(a,b) leading-term pairs with a*b within ±5 of target: {len(ab_pairs)}")
for (a,b) in ab_pairs[:10]:
    print(f"  {a}·{b} = {a*b}")
print()

# For each leading pair, sweep two correction terms:
#    leading + s1*c*X + s2*e/f
for (a, b) in ab_pairs:
    leading = a * b
    for Xname, Xval in CONT_BLOCKS.items():
        for c in SMALL_INTS:
            term2_pos = c * Xval
            for s1 in (+1, -1):
                term2 = s1 * term2_pos
                for e in SMALL_INTS:
                    for f in SMALL_INTS:
                        if f == 0: continue
                        for s2 in (+1, -1):
                            term3 = s2 * e / f
                            val = leading + term2 + term3
                            checks += 1
                            err = abs(val - target)
                            for tname, twin in tiers.items():
                                if err < twin:
                                    hits[tname] += 1
                                    if len(hits_examples[tname]) < 5:
                                        sgn1 = "+" if s1>0 else "-"
                                        sgn2 = "+" if s2>0 else "-"
                                        formula = f"{a}·{b} {sgn1} {c}·{Xname} {sgn2} {e}/{f}"
                                        hits_examples[tname].append((formula, val, (val-target)/target*1e9))

print(f"Total formulas evaluated: {checks:,}")
print()

# =====================================================================
# Report
# =====================================================================
print("="*78)
print(f" CHANCE-COINCIDENCE TABLE (search basis: {len(ab_pairs)} a·b pairs × cX × e/f)")
print("="*78)
print(f"{'Window':>10s} {'#Hits':>8s} {'P(chance)':>12s} {'lnB':>7s} {'σ-eq':>6s}")
for tname, twin in tiers.items():
    h = hits[tname]
    p = h / checks if h > 0 else 1/checks
    if h > 0:
        lnB = math.log(checks/h)
        sigma_eq = math.sqrt(2*lnB) if lnB > 0 else 0
    else:
        lnB = math.log(checks)  # upper bound
        sigma_eq = math.sqrt(2*lnB)
    print(f"{tname:>10s} {h:>8d} {p:>12.3e} {lnB:>7.2f} {sigma_eq:>6.2f}")

print()
print("Examples of formulas hitting the target at each tier:")
for tname in tiers:
    print(f"\n  {tname} window:")
    if hits_examples[tname]:
        for (form, val, ppb) in hits_examples[tname][:3]:
            print(f"    {form:50s} = {val:.6f}   ({ppb:+.2f} ppb)")
    else:
        print(f"    [no hits in {checks:,} formulas]")

# =====================================================================
# DCT formula leverage
# =====================================================================
print()
print("="*78)
print(" DCT FORMULA POSITION IN THE TIER STRUCTURE")
print("="*78)
DCT_diff_ppb = (DCT - target)/target * 1e9
print(f"  DCT formula: 12·153 + 4μ² + 1/144 = {DCT:.6f}")
print(f"  Δ from CODATA: {DCT_diff_ppb:+.2f} ppb")
print(f"  Tier hit: {'92 ppb' if abs(DCT_diff_ppb) < 92 else '100 ppm or worse'}")