#!/usr/bin/env python3
"""
ERDŐS PRIME DIVISIBILITY CONJECTURE — Complete Verification Suite
=================================================================
Verifies EVERY computational claim in the paper. Run this script;
if it prints "ALL CLAIMS VERIFIED" at the end, every computation
in the paper is independently confirmed.

Claims verified:
  1. Zero counterexamples for n ≤ 4400, i ≤ 15  (§5.7, Appendix A)
  2. Exactly 1 fully obstructed triple: (10,3,5)  (§5.5, Appendix A)
  3. B(i) table: B(3)=9, B(5)=81, B(7)=4375      (§5.6)
  4. Appendix B: 1 FO, 8 Band Escape, 1 Loneliness Escape  (Appendix B)
  5. j_0(i) thresholds for i=4..9                  (§5.6 table)
  6. Brun–Titchmarsh: 2i/ln(i) < i-1 for i ≥ 10   (§5.4 Corollary)
  7. i=3: only (10,3,5) is FO, for ALL j ≤ 10000   (§5.5)
  8. i=4..9: zero FO triples, j ≤ 1000              (§5.6)

Runtime: ~10-20 minutes depending on hardware.
"""
import sys
import time
import math
from collections import defaultdict

# ═══════════════════════════════════════════════════════════════════
# PRIMITIVES
# ═══════════════════════════════════════════════════════════════════

def sieve_primes(limit):
    is_prime = [True] * (limit + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(limit**0.5) + 1):
        if is_prime[i]:
            for j in range(i*i, limit + 1, i):
                is_prime[j] = False
    return [i for i in range(2, limit + 1) if is_prime[i]]

def vp(n, p):
    """p-adic valuation of n."""
    if n == 0: return float('inf')
    v = 0
    while n % p == 0: n //= p; v += 1
    return v

def vp_binom(n, k, p):
    """p-adic valuation of C(n,k) via Legendre."""
    val = 0; pp = p
    while pp <= n:
        val += n // pp - k // pp - (n - k) // pp
        pp *= p
    return val

def is_smooth(n, bound, primes):
    """Is n smooth with respect to primes <= bound?"""
    rem = n
    for p in primes:
        if p > bound: break
        while rem % p == 0: rem //= p
    return rem == 1

def primepi(x, primes):
    """Count primes <= x."""
    count = 0
    for p in primes:
        if p > x: break
        count += 1
    return count

PRIMES = sieve_primes(100000)

# ═══════════════════════════════════════════════════════════════════
# CORE FUNCTIONS
# ═══════════════════════════════════════════════════════════════════

def find_witness(n, i, j):
    """Find prime p >= i with p | gcd(C(n,i), C(n,j)), or None."""
    for p in PRIMES:
        if p < i: continue
        if p > n: break
        if vp_binom(n, i, p) >= 1 and vp_binom(n, j, p) >= 1:
            return p
    return None

def is_fully_obstructed(n, i, j):
    """
    Check if (n,i,j) is fully obstructed:
    - i-product is j-smooth
    - every prime > i in i-product is tame (in (j-i,j]), lonely, v=1
    - no Band Escape (no prime in (i, j-i])
    """
    for k in range(i):
        if not is_smooth(n - k, j, PRIMES):
            return False

    for p in PRIMES:
        if p <= i: continue
        if p > j: break
        for k in range(i):
            if (n - k) % p == 0:
                if p <= j - i:
                    return False          # Band escape
                for k2 in range(i, j):
                    if (n - k2) % p == 0:
                        return False      # Not lonely
                if vp(n - k, p) >= 2:
                    return False          # Bridge applies
                break
    return True

# ═══════════════════════════════════════════════════════════════════
# CLAIM 1 & 2: Main verification + fully obstructed classification
# ═══════════════════════════════════════════════════════════════════

def verify_claim_1_and_2(n_max=4400, i_max=15):
    """Verify zero counterexamples and find all FO triples."""
    print("=" * 70)
    print(f"CLAIM 1: Zero counterexamples for n <= {n_max}, i <= {i_max}")
    print(f"CLAIM 2: Identify all fully obstructed triples")
    print("=" * 70)

    total = 0; fails = 0; fo_triples = []
    t0 = time.time()

    for n in range(3, n_max + 1):
        for j in range(2, n // 2 + 1):
            if n < 2 * j: continue
            for i in range(1, min(j, i_max + 1)):
                total += 1
                w = find_witness(n, i, j)
                if w is None:
                    fails += 1
                    print(f"  *** COUNTEREXAMPLE: ({n},{i},{j}) ***")
                elif is_fully_obstructed(n, i, j):
                    fo_triples.append((n, i, j, w))
        if n % 500 == 0:
            print(f"  n={n}: {total:,} triples, {fails} fails, "
                  f"{time.time()-t0:.0f}s", file=sys.stderr)

    elapsed = time.time() - t0
    print(f"  Triples checked: {total:,}")
    print(f"  Counterexamples: {fails}")
    print(f"  Fully obstructed: {len(fo_triples)}")
    for (n, i, j, w) in fo_triples:
        print(f"    ({n}, {i}, {j}) -> witness p={w}")
    print(f"  Time: {elapsed:.0f}s")

    ok1 = (fails == 0)
    ok2 = (len(fo_triples) == 1 and fo_triples[0][:3] == (10, 3, 5))
    print(f"  CLAIM 1: {'PASS' if ok1 else 'FAIL'}")
    print(f"  CLAIM 2: {'PASS' if ok2 else 'FAIL'}")
    return ok1, ok2

# ═══════════════════════════════════════════════════════════════════
# CLAIM 3: B(i) table
# ═══════════════════════════════════════════════════════════════════

def verify_claim_3():
    print("\n" + "=" * 70)
    print("CLAIM 3: B(i) table")
    print("=" * 70)

    expected = {3: 9, 5: 81, 7: 4375}
    ok = True

    for bound in [2, 3, 5, 7, 11]:
        max_pair = 0
        for x in range(2, 200001):
            if is_smooth(x, bound, PRIMES) and is_smooth(x - 1, bound, PRIMES):
                max_pair = x
        S = [p for p in PRIMES if p <= bound]
        status = ""
        if bound in expected:
            if max_pair == expected[bound]:
                status = "  matches paper"
            else:
                status = f"  MISMATCH (paper says {expected[bound]})"
                ok = False
        print(f"  S = primes <= {bound:>2} = {S}: "
              f"B = {max_pair:>6}  (pair: {max_pair-1}, {max_pair}){status}")

    print(f"  CLAIM 3: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# CLAIM 4: Appendix B — the 10 "hard" triples
# ═══════════════════════════════════════════════════════════════════

def verify_claim_4():
    print("\n" + "=" * 70)
    print("CLAIM 4: Appendix B — 10 hard triples classification")
    print("  Expected: 1 FO, 8 Band Escape, 1 Loneliness Escape")
    print("=" * 70)

    triples = [(10,3,5), (16,3,7), (16,3,8), (22,3,11), (26,3,13),
               (27,3,13), (28,3,13), (27,4,13), (28,4,13), (28,5,13)]

    fo_count = 0; band_count = 0; lonely_count = 0; other_count = 0

    for (n, i, j) in triples:
        w = find_witness(n, i, j)
        fo = is_fully_obstructed(n, i, j)

        has_band = False
        for p in PRIMES:
            if p <= i or p > j - i: continue
            for k in range(i):
                if (n - k) % p == 0:
                    has_band = True; break
            if has_band: break

        # Check loneliness escape: tame prime with 2nd multiple in (j-i)-block
        has_lonely_esc = False
        if not fo and not has_band:
            for p in PRIMES:
                if p <= i or p > j: continue
                if p <= j - i: continue  # would be band
                hit_i = False
                for k in range(i):
                    if (n - k) % p == 0: hit_i = True; break
                if not hit_i: continue
                for k2 in range(i, j):
                    if (n - k2) % p == 0:
                        has_lonely_esc = True; break
                if has_lonely_esc: break

        if fo:
            fo_count += 1; label = "FULLY OBSTRUCTED"
        elif has_band:
            band_count += 1; label = "Band Escape"
        elif has_lonely_esc:
            lonely_count += 1; label = "Loneliness Escape"
        else:
            other_count += 1; label = "other escape"

        print(f"  ({n:>2},{i},{j:>2}): witness={w}, {label}")

    ok = (fo_count == 1 and band_count == 8 and lonely_count == 1)
    print(f"  Fully obstructed: {fo_count} (expected 1)")
    print(f"  Band Escape: {band_count} (expected 8)")
    print(f"  Loneliness Escape: {lonely_count} (expected 1)")
    print(f"  CLAIM 4: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# CLAIM 5: j_0(i) thresholds for i=4..9
# ═══════════════════════════════════════════════════════════════════

def verify_claim_5():
    print("\n" + "=" * 70)
    print("CLAIM 5: j_0(i) thresholds (section 5.6 table)")
    print("=" * 70)

    expected_j0 = {4: 6, 5: 6, 6: 7, 7: 8, 8: 9, 9: 10}
    ok = True

    for i in range(4, 10):
        target = i - 2
        exc = []
        for j in range(i + 1, 10001):
            count = primepi(j, PRIMES) - primepi(j - i, PRIMES)
            if count > target:
                exc.append(j)
        j0 = (max(exc) + 1) if exc else (i + 1)
        matches = (j0 == expected_j0[i])
        if not matches: ok = False
        exc_str = str(exc) if len(exc) <= 5 else f"{exc[:3]}... ({len(exc)} values)"
        print(f"  i={i}: j_0 = {j0:>6} (expected {expected_j0[i]}), "
              f"exc j: {exc_str}  {'matches' if matches else 'MISMATCH'}")

    print(f"  CLAIM 5: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# CLAIM 6: Brun–Titchmarsh: 2i/ln(i) < i-1 for i >= 10
# ═══════════════════════════════════════════════════════════════════

def verify_claim_6():
    print("\n" + "=" * 70)
    print("CLAIM 6: 2i/ln(i) < i-1 for all i >= 10")
    print("=" * 70)

    ok = True
    for i in range(3, 21):
        lhs = 2 * i / math.log(i)
        rhs = i - 1
        holds = lhs < rhs
        if i >= 10 and not holds:
            ok = False
        indicator = "  <-- threshold" if i == 10 else ""
        print(f"  i={i:>2}: 2i/ln(i) = {lhs:>7.3f}, i-1 = {rhs:>2}, "
              f"{'holds' if holds else 'fails'}{indicator}")

    # Bulk check i=10..10000
    for i in range(10, 10001):
        if 2 * i / math.log(i) >= i - 1:
            print(f"  FAILS at i={i}!")
            ok = False

    if ok:
        print(f"  Verified for i = 10..10000: always holds.")
    print(f"  CLAIM 6: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# CLAIM 7: i=3 — only (10,3,5) is FO, for ALL j <= 10000
# ═══════════════════════════════════════════════════════════════════

def verify_claim_7():
    print("\n" + "=" * 70)
    print("CLAIM 7: i=3, only (10,3,5) is fully obstructed (j <= 10000)")
    print("=" * 70)

    i = 3; fo_list = []; t0 = time.time()

    for j in range(4, 10001):
        for n in range(2 * j, min(2 * j + 200, 5 * j + 1)):
            smooth = True
            for k in range(i):
                if not is_smooth(n - k, j, PRIMES):
                    smooth = False; break
            if not smooth: continue
            if is_fully_obstructed(n, i, j):
                w = find_witness(n, i, j)
                fo_list.append((n, i, j, w))
        if j % 2000 == 0:
            print(f"  j={j}: {len(fo_list)} FO, "
                  f"{time.time()-t0:.0f}s", file=sys.stderr)

    print(f"  FO triples with i=3, j <= 10000: {len(fo_list)}")
    for (n, i, j, w) in fo_list:
        print(f"    ({n}, {i}, {j}) -> witness p={w}")

    ok = (len(fo_list) == 1 and fo_list[0][:3] == (10, 3, 5))
    print(f"  Time: {time.time()-t0:.0f}s")
    print(f"  CLAIM 7: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# CLAIM 8: i=4..9 — zero FO triples, j <= 1000
# ═══════════════════════════════════════════════════════════════════

def verify_claim_8():
    print("\n" + "=" * 70)
    print("CLAIM 8: i=4..9, zero fully obstructed triples (j <= 1000)")
    print("=" * 70)

    ok = True; t0 = time.time()

    for i in range(4, 10):
        fo_count = 0
        for j in range(i + 1, 1001):
            for n in range(2 * j, min(2 * j + 200, 5 * j + 1)):
                smooth = True
                for k in range(i):
                    if not is_smooth(n - k, j, PRIMES):
                        smooth = False; break
                if not smooth: continue
                if is_fully_obstructed(n, i, j):
                    fo_count += 1
                    w = find_witness(n, i, j)
                    print(f"    ({n},{i},{j}) is FO! witness={w}")

        status = "zero FO" if fo_count == 0 else f"{fo_count} FO -- UNEXPECTED"
        print(f"  i={i}: {status}")
        if fo_count > 0: ok = False

    print(f"  Time: {time.time()-t0:.0f}s")
    print(f"  CLAIM 8: {'PASS' if ok else 'FAIL'}")
    return ok

# ═══════════════════════════════════════════════════════════════════
# MAIN
# ═══════════════════════════════════════════════════════════════════

def main():
    print("+" + "=" * 68 + "+")
    print("|  ERDOS PRIME DIVISIBILITY CONJECTURE — FULL VERIFICATION          |")
    print("|  Checks every computational claim in the paper.                   |")
    print("+" + "=" * 68 + "+")
    print()

    results = {}
    t_start = time.time()

    # Fast claims first
    results[3] = verify_claim_3()
    results[4] = verify_claim_4()
    results[5] = verify_claim_5()
    results[6] = verify_claim_6()
    # Medium claims
    results[8] = verify_claim_8()
    results[7] = verify_claim_7()
    # Slow claim (main verification)
    r1, r2 = verify_claim_1_and_2()
    results[1] = r1
    results[2] = r2

    # Final report
    elapsed = time.time() - t_start
    print("\n" + "=" * 70)
    print("FINAL REPORT")
    print("=" * 70)
    all_pass = True
    for claim in sorted(results.keys()):
        status = "PASS" if results[claim] else "FAIL"
        if not results[claim]: all_pass = False
        print(f"  Claim {claim}: {status}")

    print(f"\n  Total time: {elapsed:.0f}s ({elapsed/60:.1f} min)")
    print()

    if all_pass:
        print("  +========================================+")
        print("  |   ALL CLAIMS VERIFIED SUCCESSFULLY     |")
        print("  +========================================+")
    else:
        print("  +========================================+")
        print("  |   SOME CLAIMS FAILED — SEE ABOVE       |")
        print("  +========================================+")

    return 0 if all_pass else 1

if __name__ == "__main__":
    sys.exit(main())