Comprobación de números primos con el FIPS186

Que no se diga, aquí está el código del algoritmo y una comprobación para ver el número de falsos positivos y negativos:

---

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
import random

""" FIPS186
Step 1. Set i = 1 and n > or = to 50.
Step 2. Set w = the integer to be tested, w = 1 + 2 a m, where m is odd and 2 a is m, where m is odd and 2 is the largest power of 2 dividing w - 1

Step 3. Generate a random integer b in the range 1 < b < w
Step 4. Set j = 0 and z = b**m mod w mod w
Step 5. If j = 0 and z = 1, or if z = w - 1, go to step 9
Step 6. If j > 0 and z = 1, go to step 8
Step 7. j = j + 1. If j < a, set z = z 2 mod w and go to step 5 mod w and go to step 5

Step 8. w is not prime. Stop.
Step 9. If i < n, set i = i + 1 and go to step 3. Otherwise, w is probably prime
"""

from math import sqrt
# Comprobación real
def rchk(n):
    for i in xrange(2,int(sqrt(n))+1):
        if not n % i:
            return False
    return True
     
# Especificada en el FIPS186
def chkprime(w, n = 50):
    m = w-1
    a = 0
    while not m & 1:
        m >>= 1
        a += 1
    i = 0
    while True:
        b = random.randint(2, w-1)
        j = 0
        while True:
            z = (b**m) % w
            if ( j == 0 and z == 1) or \
              z == w-1:
                if (i<n):
                    i += 1
                    break

                else:
                    return True

            if j > 0 and z == 1:
                return False

            j += 1
            if j < a:
                z = ( z * z ) % w

            else:
                return False

fn = 0
fp = 0
pr = 0
inicio = 10
fin = 10000
import sys
for i in xrange(inicio,fin):
    c = chkprime(i)
    r = rchk(i)
    if not i%100:
        sys.stdout.write("\r"+str(i)+" "+str(pr))
        sys.stdout.flush()
    if c:
        pr += 1
    if c != r:
        if c:
            fp += 1
            pr -= 1
        else:
            fn += 1
            pr += 1

print ""
print "Primos:", pr
print "Falsos negativos:", fn
print "Falsos positivos:", fp

---

Por ejemplo:

---

Primos: 1225
Falsos negativos: 608
Falsos positivos: 0

---

ps: se me pasó poner la licencia, como siempre está bajo WTFPL

Saludos