primes.py

from math import floor
import numpy as np

def primes_upto(n):
    maxdiv = floor(n ** 0.5) + 1
    divs = np.arange(2, maxdiv)
    nums = np.arange(2, n + 1)
    rem = np.mod.outer(nums, divs)
    mask = np.less.outer(nums, divs ** 2)
    isprime = np.all(mask | rem, axis=1)
    return np.argwhere(isprime).ravel() + 2

if __name__ == '__main__':
    print(primes_upto(100))