Looking for Primes

A post in the TopCoder forum led me to this problem. The crux of the problem is to check whether a number is prime, the number being in the range [3,2^32]. It was obvious the naive prime check would time out. At first I thought simply generating all primes under sqrt(2^32) using sieve and using those to test for primality would be enough. This method usually works for problems with numbers in this range. Code (in Python) for both methods follows.

"Naive" Algorithm

from math import *
def is_prime(n):
limit=int(sqrt(n))
if n==2:
 # 2 is the only even prime
 return True
 elif (n&1)==0:
 # n is not 2 and it is even, so it can't be a prime
 return False
 else:
 for i in xrange(3,limit,2): # check only odd divisors
 if n%i==0:
 # n is divisible i which is neither n nor 1
 # it can't be prime
 return False

  return True

Faster Method using Sieve of Eratosthenes

def generate_primes(max_value):
 #run sieve
prime=[True]*max_value # assume every number is prime
prime[0]=False # 0 and 1 are not primes
prime[1]=False
limit=int(sqrt(max_value)) # we only need to look unitl sqrt(n)

 #all even numbers over 2 are non-prime
 for i in xrange(4,max_value,2);
prime[i]=False

 prime_list=[2] # start off list with one prime - 2 
 for(i in xrange(3,limit,2)): # check all odd numbers
 if prime[i]:
prime_list.append(i) # add i to our list
 # cross out all multiples of i
 for j in xrange(i*i,max_value)
prime[j]=False

 # add all primes greater than sqrt to our list
 if limit&1==0:
limit+=1
prime_list.extend([i for i in xrange(limit,max_value,2) if prime[i]])

 return prime_list

def is_prime(prime_list, n):
limit=int(sqrt(n))
 # check all primes below sqrt(n)
 for p in prime_list[:n]:
 if (n%p)==0:
 return False
 #no prime divides it, so it must be a prime
 return True

The simpler method checks all odd numbers in [3,floor(sqrt(n))] for divisibility, while the faster method only checks the primes in the same range. This results in a significant speedup, and is usually fast enough for 32 bit primes in most problems.

It turns out that there are in excess of 6000 primes less than sqrt(2^32), and the maximum distance between a number and it nearest prime in [3,2^32] can be larger than 200. So we can have 200 prime checks, each having to use the %-operator 6000 times giving us a runtime in the order of 10^5 for each case. Notice that there are 10^5 cases, thus our program will run for (~10^5*10^5=)~10^10 iterations. In other words, in excess of 20 minutes on the SPOJ servers while the time limit is 5 seconds.

So the usual prime checks are not fast enough. Time to bring in the heavy artillery. The only practical fast methods of checking big primes use some form of probabilistic check. That means that there is a non-zero probability that your answer is wrong. But usually, if you know what you are doing, that chance is very very low. Since Java already has a ready-made probablyPrime method in its BigInteger class, I thought it would be a good idea to just use that than code something from scratch. I guess the problemsetters had thought of that too, and my simple Java solution timed out. The huge input/output involved in 10^5 cases was probably the dominant factor, but I didn’t think it would be worth it to try the weird optimized Java I/O code. There was nothing more to do than to code a probabilistic prime tester myself.

Probabilistic prime testers, despite their grandiose sounding names, are actually quite simple. The Miller-Rabin test is generally thought to be faster in practice. It uses that the fact (Fermat’s Little Theorem) that

for all primes p and all bases a
a^(p-1)=1 (mod p)

So, if a number is prime a^(p-1)%p will always be 1. If a number n gives a^(n-1)%n!=1, it most certainly isn’t prime. Thus simply checking a large enough number of random values of a gives us an answer which has a high probability of correctness. After all, if n us composite and we check enough values of a, sooner or later we will get a result that is not 1.

Unfortunately, there is a catch. There are certain rare composite numbers that pass this test for all bases. To guard against that, we use another property of primes.

For any prime p, if
x^2=1 (mod p)
x has to be either 1 or -1 (mod p)
i.e. There are no non-trivial roots of unity (mod p) if p is a prime

Thus, if n passes Fermat’s test for a reasonable number of bases, and there is no x other than 1 and p-1 for which (x^2)%p=1, n is probably prime. With a large enough sampling for the values of the base, this compound test gives an answer that is correct with high probability.

Computing the value of a^(p-1)%p is easy using exponentiation by squaring, effectively taking log(p) time and doing it while keeping all intermediate values less than p. Care still has to be taken to avoid overflows while squaring. Checking for non-trivial roots is slightly trickier. Actually, the fact that we are using exponentiation by squaring comes to our rescue. As we exponentiate, whenever we get an even exponent, we simply square the result of exponent/2. Now, suppose the result of a^(exponent/2)%p is x. We use that to compute a^(exponent)%p=square(a^(exponent/2)%m)=x*x. At this stage, if x^2%p is 1, x is a root of unity (mod p). So, x has to be either 1 or p-1. If x is anything else, we can be sure p is not a prime. Thus simply putting this check inside our modular exponentiation code completes or test. Although I used C for this particular problem for speed, I’m pasting Python Code for it. I don’t like posting code for SPOJ problems and all the other code in this post is in Python anyway.

Miller-Rabin Primality Test

from random import *

def prime(n,k): # check if n is prime
limit=k # k samples
q=n-1
 if n==2: return True
 elif ((n&1)==0): return False
 elif (n<2): return False

  for i in xrange(limit):
a=randrange(1,n)
m=q;y=1
z=a
 #modular exponentiation
 while m>0:
 while (m&1)==0:
 #even exponent
 # a^m%n=square(a^(m/2)%n)
x=z;z=(z*z)
 if (z>=n):
z%=n
 #check for roots of unity
 if (z==1 and x!=1 and x!=q):
 return False
m>>=1
 #odd exponent
 # a^m%n=(a*a^(m1))%n
m-=1;y=(y*z)
 if y>=n:
y%=n

  if y!=1:
 # certainly not prime
 return False
 #survived k iterations, n is probably prime
 return True

Usually around 20 iterations (i.e. the value of k) is enough. At SPOJ I experimented , successfully, with much lower values. The probability of an incorrect answer is 4^(-k).

(moving this post from LiveJournal messed up the indentaion in the code, which I am too lazy to fix)

Further Reading:
Fermat’s Little Theorem (Wikipedia)
Primality Test (Wikipedia)
Prime Numbers, Factorizations and Euler Function (Topcoder)
Primality Testing : Non-deterministic Algorithms (TopCoder)

Advertisement