Random Prime Generator/Modular Arithmetic

[Bryan Olson]

Tuvas wrote:
> Okay, now I get the correct number of 561 pseudoprimes, 5, so I can
> assume that it is indeed working right.

Your code now looks right, so I'm guessing "5" was a typo,
perhaps because "5" is just below "8" on the numeric keypad.


You can simplify your loop like:

def is_strong_pseudo_prime(a, n):
     # check for n == 2 and/or other input validation omitted
     if pow(a, n - 1, n) != 1:
         return False
     x = n - 1
     while x % 2 == 0:
         x = x / 2
         y = pow(a, x, n)
         if y == n - 1:
             return True
         elif y != 1:
             return False
     return True


For time efficiency, I prefer a slightly longer version that does
all the square-root-checking in the course of the one modular
exponentiation.

def is_mr_pseudoprime(n, w):
     """ Pass positive integers n and w, with w < n.
         Return true iff n is prime or a strong base-w pseudo-prime.
     """
     assert n == long(n) and w == long(w) and 1 < w < n
     power = n - 1
     pow2 = 0
     while power % 2 == 0:
         power /= 2
         pow2 += 1
     r = pow(w, power, n)
     for _ in xrange(pow2):
         prev = r
         r = r * r % n
         if r == 1:
             return prev in (1, n - 1)
     return r == 1


-- 
--Bryan