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))))
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)