#!/usr/bin/env python3
"""
Erdos Problem 835 -- exact reduction + divisibility sieve.

835(k): can the k-subsets of [2k] be (k+1)-coloured so that every (k+1)-subset
        sees all k+1 colours on its k facets?

Reduction (proved in the note):
  835(k) is YES
    <=>  chi(J(2k,k)) = k+1               (Johnson graph chromatic number)
    <=>  A(2k,4,k) = C_k  AND a perfect tiling exists
    <=>  the middle layer C([2k],k) partitions into (k+1) Steiner systems S(k-1,k,2k)
         (a "large set" LS[k+1](k-1,k,2k)).

Here C_k = Catalan(k) = binom(2k,k)/(k+1) is the exact packing bound on any single
colour class (a distance-4 constant-weight code): each block covers k of the
binom(2k,k-1) many (k-1)-sets, so a class has <= binom(2k,k-1)/k = C_k codewords,
with equality IFF it is a Steiner system S(k-1,k,2k).

RIGOROUS CHAIN used to certify NO:
  no S(k-1,k,2k)  =>  A(2k,4,k) <= C_k - 1
                  =>  chi >= ceil( binom(2k,k) / (C_k - 1) ) > k+1  =>  835(k) is NO.

This script:
  (1) computes the Steiner divisibility conditions for S(k-1,k,2k), k=2..K,
  (2) tests the equivalence  [S(k-1,k,2k) divisibility holds]  <=>  [k+1 is prime],
  (3) confirms the strict bound gap ceil(binom(2k,k)/(C_k-1)) = k+2 when no Steiner
      system exists, i.e. one missing codeword already forces chi >= k+2.
All arithmetic is exact Python big-int; no floating point.
"""
from math import comb, isqrt
import json

def is_prime(n:int)->bool:
    if n < 2: return False
    if n % 2 == 0: return n == 2
    for d in range(3, isqrt(n)+1, 2):
        if n % d == 0: return False
    return True

def catalan(k:int)->int:
    # exact integer; binom(2k,k) is divisible by k+1
    b = comb(2*k, k)
    assert b % (k+1) == 0
    return b // (k+1)

def steiner_lambdas(k:int):
    """
    Necessary divisibility for a Steiner system S(t,K,v) with t=k-1, K=k, v=2k:
    for every 0<=i<=t, lambda_i = C(v-i, t-i) / C(K-i, t-i) must be a positive integer.
    Here C(K-i,t-i) = C(k-i, (k-1)-i) = C(k-i,1) = k-i.
    Returns (all_integer:bool, list of (i, numerator, denom, exact_or_None)).
    """
    t, K, v = k-1, k, 2*k
    rows = []
    ok = True
    for i in range(0, t+1):
        num = comb(v-i, t-i)
        den = comb(K-i, t-i)     # == k-i
        assert den == k-i
        if num % den == 0:
            rows.append((i, num, den, num//den))
        else:
            rows.append((i, num, den, None))
            ok = False
    return ok, rows

def chrom_lower_from_missing_one(k:int)->int:
    """If A(2k,4,k) <= C_k - 1, chi >= ceil( binom(2k,k)/(C_k-1) )."""
    N = comb(2*k, k)
    Ck = catalan(k)
    A = Ck - 1
    # ceil division
    return -(-N // A)

def main():
    K = 300
    survivors = []
    mismatches = []
    table = []
    for k in range(2, K+1):
        ok, rows = steiner_lambdas(k)
        prime = is_prime(k+1)
        # locate first failing i (for the certificate)
        first_fail = next((r for r in rows if r[3] is None), None)
        table.append({
            "k": k, "k_plus_1": k+1, "kplus1_prime": prime,
            "steiner_S(k-1,k,2k)_divisible": ok,
            "first_failing_i": (None if ok else first_fail[0]),
            "first_failing_lambda": (None if ok else f"{first_fail[1]}/{first_fail[2]}"),
        })
        if ok:
            survivors.append(k)
        if ok != prime:
            mismatches.append((k, ok, prime))

    print("="*74)
    print("STEINER DIVISIBILITY SIEVE for S(k-1,k,2k), k = 2..%d" % K)
    print("="*74)
    print("Claim under test:  S(k-1,k,2k) satisfies design divisibility  <=>  k+1 is prime")
    print("Mismatches found:", mismatches if mismatches else "NONE  (equivalence holds on 2..%d)"%K)
    print()
    print("Survivors (divisibility passes) up to %d:" % K)
    print("  k        :", survivors)
    print("  k+1      :", [k+1 for k in survivors])
    print("  all k+1 prime? ->", all(is_prime(k+1) for k in survivors))
    print()

    # Show the clean odd-k obstruction: i=k-2 gives lambda=(k+2)/2
    print("Worked obstruction for a few composite k+1 (first failing design equation):")
    for k in [3,5,7,8,9,11,13,14,15]:
        ok, rows = steiner_lambdas(k)
        ff = next((r for r in rows if r[3] is None), None)
        reason = f"lambda_{ff[0]} = C({2*k-ff[0]},{(k-1)-ff[0]})/{ff[2]} = {ff[1]}/{ff[2]} (not integer)"
        print(f"  k={k:3d} (k+1={k+1:3d}, {'prime' if is_prime(k+1) else 'composite'}): NO  -- {reason}")
    print()

    # The "one missing codeword forces chi >= k+2" fact
    print("Strict-gap check: if A(2k,4,k) <= C_k - 1 then chi >= ceil(binom(2k,k)/(C_k-1)):")
    for k in [3,5,7,8,9,10,16,112]:
        lb = chrom_lower_from_missing_one(k)
        print(f"  k={k:4d}:  ceil( binom(2k,k) / (C_k - 1) ) = {lb}   (= k+2 ? {lb==k+2})")
    print()

    # Frontier: surviving k up to 500 (candidates that the divisibility sieve cannot kill)
    surv500 = [k for k in range(3, 501) if steiner_lambdas(k)[0]]
    print("Open YES-candidates (sieve cannot kill) with 3<=k<=500  ==  {k : k+1 prime}:")
    print("  first several k:", surv500[:12], "...")
    print("  smallest:", surv500[0], " (k+1 =", surv500[0]+1, ", prime:", is_prime(surv500[0]+1), ")")

    with open("sieve_835_table.json","w") as f:
        json.dump({"K":K,"survivors":survivors,"mismatches":mismatches,"table":table}, f, indent=1)
    print("\nWrote sieve_835_table.json")

if __name__ == "__main__":
    main()