pythonnumpyprimesprimality-test

Numpy AKS primality function


I'm implementing the AKS-primality test using Numpy. Specifically, I'm using the polynomial capabilities of Numpy to find the coefficients of the equation (x - 1)^n - (x^n - 1) and then returning True if all those coefficients are divisible by the prime candidate.

def isPrime(n):
    if n < 2:
        return False
    poly1 = np.polynomial.polynomial.polypow([1, -1], n)
    poly2 = np.zeros(n + 1, dtype=int)
    poly2[0] = 1
    poly2[-1] = -1
    coefficients = np.polysub(poly1, poly2)
    divisibility_test = lambda x : x % n != 0
    non_divisibles = filter(divisibility_test, coefficients)
    try:
        _ = next(non_divisibles)
        return False
    except StopIteration:
        return True

I recognize I don't have the most performant solution here. I'm confused why this implementation produces correct answers only for inputs below 59. It fails to recognize any primes greater than 53.

Edit: An easy way to demonstrate the results looks like this:

print(list(filter(isPrime, range(100))))

Solution

  • To use python's arbitrary precision integers you can specify dtype=object to the appropriate methods when constructing coefficient lists, i.e.

    poly1 = np.polynomial.polynomial.polypow(np.array([1, -1], dtype=object), n)
    poly2 = np.zeros(n + 1, dtype=object)