I call shenanigans on this.

    import timeit
    def test(M,n):
        values = [i * M for i in range(1, n + 1)]
        s = set(values)
        sum(v in s for v in values)
    M = (1 << 61) - 1
    for n in [1000, 2000, 4000, 8000, 16000]:
        print(f"M=2^61-1, {n=:5d} ->", timeit.timeit(lambda: test(M,n), number=3))
    for n in [1000, 2000, 4000, 8000, 16000]:
        print(f"M=1,      {n=:5d} ->", timeit.timeit(lambda: test(1,n), number=3))

Magically, when you stop using BIGINTS as the set members and just use regular ints, there is no such quadratic explosion.

    M=2^61-1, n= 1000 -> 0.13144792200182565
    M=2^61-1, n= 2000 -> 0.48016051898594014
    M=2^61-1, n= 4000 -> 2.058760045998497
    M=2^61-1, n= 8000 -> 7.843778470996767
    M=2^61-1, n=16000 -> 40.01485426299041
    M=1,      n= 1000 -> 0.000748768012272194
    M=1,      n= 2000 -> 0.0015750699967611581
    M=1,      n= 4000 -> 0.003037029004190117
    M=1,      n= 8000 -> 0.006664915999863297
    M=1,      n=16000 -> 0.012693285010755062

The runtime is being spent hashing bigints, comparing candidate bigint(s) against reference bigints, and summing bigints. And there's also some set lookups.