1

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))))
mttpgn
  • 327
  • 6
  • 17
  • 1
    I think this has to do with integer overflow. Raising something to the 59'th power will result in very large values. – Kevin Apr 02 '21 at 04:11
  • @Kevin Same result for dtype=float – mttpgn Apr 02 '21 at 09:21
  • Try `dtype=object`. – President James K. Polk Apr 02 '21 at 13:43
  • @PresidentJamesK.Polk That similarly limits the prime output to no more than the number 53. While dtype is an optional argument to zeros, it isn't for polypow. – mttpgn Apr 02 '21 at 15:15
  • Like this: `poly1 = np.polynomial.polynomial.polypow(np.array([1, -1], dtype=object), n)` – President James K. Polk Apr 02 '21 at 21:36
  • @PresidentJamesK.Polk If you put that as an answer I'll accept – mttpgn Apr 02 '21 at 22:34
  • I understand the question is more about getting what you have working, but a pedantic note that you're implementing the old exponential-time method, rather than AKS. See the algorithm part of your linked Wikipedia page. The cool thing about AKS is both limiting the size and number of coefficients to vastly less than the naive method. – DanaJ Apr 03 '21 at 09:04
  • 1
    @DanaJ The question was about why Numpy was exhibiting certain behavior. "Working" here helps me see why I was running into the constraint I was. Learning the answer shouldn't imply this implementation is totally algorithmially correct yet. – mttpgn Apr 03 '21 at 14:45
  • @mttpgn That's why I phrased my comment the way I did. You did start your comment saying you were implementing AKS. Imagine if you said you were implementing Quicksort and proceeded to show a Bubblesort algorithm. Wouldn't you expect someone to note that your code isn't the right algorithm? – DanaJ Apr 04 '21 at 02:43

1 Answers1

1

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)
President James K. Polk
  • 40,516
  • 21
  • 95
  • 125