37
votes

A few years ago, it was proven that PRIMES is in P. Are there any algorithms implementing their primality test in Python? I wanted to run some benchmarks with a naive generator and see for myself how fast it is. I'd implement it myself, but I don't understand the paper enough yet to do that.

2
AKS is of great theoretical importance but its performance is terrible (Miller-Rabin is much better). There are many primality tests implemented in Python. - jfs
Sorry to disappoint you but there is a little mathematical misunderstanding here! Although we haven't formally determined if P == NP, but even if P != NP, it doesn't follow that P == Fast! - Tomasz Gandor

2 Answers

55
votes

Quick answer: no, the AKS test is not the fastest way to test primality. There are much much faster primality tests that either assume the (generalized) Riemann hypothesis and/or are randomized. (E.g. Miller-Rabin is fast and simple to implement.) The real breakthrough of the paper was theoretical, proving that a deterministic polynomial-time algorithm exists for testing primality, without assuming the GRH or other unproved conjectures.

That said, if you want to understand and implement it, Scott Aaronson's short article might help. It doesn't go into all the details, but you can start at page 10 of 12, and it gives enough. :-) There is also a list of implementations (mostly in C++) here.

Also, for optimization and improvements (by several orders of magnitude), you might want to look at this report, or (older) Crandall and Papadopoulos's report, or (older still) Daniel J Bernstein's report. All of them have fairly detailed pseudo-code that lends itself well to implementation.

-6
votes

Yes, go look at AKS test for primes page on rosettacode.org

def expand_x_1(p):
    ex = [1]
    for i in range(p):
        ex.append(ex[-1] * -(p-i) / (i+1))
    return ex[::-1]

def aks_test(p):
    if p < 2: return False
    ex = expand_x_1(p)
    ex[0] += 1
    return not any(mult % p for mult in ex[0:-1])
    print('# p: (x-1)^p for small p')
    for p in range(12):
        print('%3i: %s' % (p, ' '.join('%+i%s' % (e, ('x^%i' % n) if n else '')
                                   for n,e in enumerate(expand_x_1(p)))))

print('\n# small primes using the aks test')
print([p for p in range(101) if aks_test(p)])

and the output is:

# p: (x-1)^p for small p
  0: +1
  1: -1 +1x^1
  2: +1 -2x^1 +1x^2
  3: -1 +3x^1 -3x^2 +1x^3
  4: +1 -4x^1 +6x^2 -4x^3 +1x^4
  5: -1 +5x^1 -10x^2 +10x^3 -5x^4 +1x^5
  6: +1 -6x^1 +15x^2 -20x^3 +15x^4 -6x^5 +1x^6
  7: -1 +7x^1 -21x^2 +35x^3 -35x^4 +21x^5 -7x^6 +1x^7
  8: +1 -8x^1 +28x^2 -56x^3 +70x^4 -56x^5 +28x^6 -8x^7 +1x^8
  9: -1 +9x^1 -36x^2 +84x^3 -126x^4 +126x^5 -84x^6 +36x^7 -9x^8 +1x^9
 10: +1 -10x^1 +45x^2 -120x^3 +210x^4 -252x^5 +210x^6 -120x^7 +45x^8 -10x^9 +1x^10
 11: -1 +11x^1 -55x^2 +165x^3 -330x^4 +462x^5 -462x^6 +330x^7 -165x^8 +55x^9 -11x^10 +1x^11

# small primes using the aks test
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]