# Random Prime Generator/Modular Arithmetic

Bryan Olson fakeaddress at nowhere.org
Mon Mar 6 16:24:35 CET 2006

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

```