Other suggestion: square-free integers

[stla]

from gmpy2 import mpz, next_prime, f_mod


def squarefreenums(n, Nmax=None):
    """Returns the first square-free integers. A square-free integer is an 
    integer which is divisible by no other perfect square than 1.
    
    Parameters
    ----------
    n : int or None
        A positive integer. It is ignored if `Nmax` is not `None`, otherwise 
        the function returns the first `n` square-free integers.
    Nmax : int or None
        A positive integer. If it is not `None`, the function returns the first 
        square-free integers lower than `Nmax`

    Returns
    -------
    list
        A list of `mpz` integers, the first `n` square-free integers or the 
        first square-free integers lower than `Nmax`.
    
    Examples
    --------
    >>> squarefreenums(n=15)
    [mpz(1), mpz(2), mpz(3), mpz(5), mpz(6), mpz(7), mpz(10), mpz(11), mpz(13), mpz(14), mpz(15), mpz(17), mpz(19), mpz(21)]
    >>> squarefreenums(Nmax=16)
    [mpz(1), mpz(2), mpz(3), mpz(5), mpz(6), mpz(7), mpz(10), mpz(11), mpz(13), mpz(14), mpz(15)]
    >>> # this should go to 6/pi² when Nmax -> inf:
    >>> len(squarefreenums(n=None, Nmax=1000)) / 1000
    0.608
    >>> from math import pi
    >>> 6 / pi**2
    0.6079271018540267
    
    """
    if n is None and Nmax is None:
        raise ValueError("Specify an integer for either `n` or `Nmax`.")
    n = n if Nmax is None else Nmax
    out = [None]*n
    j = 0
    out[j] = mpz(1)
    primes = next_prime(0)
    prodprimes = primes
    i = mpz(2)
    while j < n-1:
        if f_mod(prodprimes, i) == 0:
            j += 1
            out[j] = i
            if(i > Nmax):
                return out[:j]
        nprime = next_prime(primes)
        prodprimes = prodprimes * nprime
        primes = nprime
        i += 1
    return out[:j]