Miller-Rabin and Primes

For those of you who do not know, Miller-Rabin is a probabilistic algorithm for testing the primality of numbers. This is different from a deterministic algorithm since with Miller-Rabin there is a change that even if the test says a number is prime that it is in fact composite. However, we can make the probability of this happening so small that we are not worried. On the other hand, if the algorithm says that a number is composite then there is no doubt that it is so.

Here is how the algorithm works. You are given a number and wish to find out if it is prime. So we choose a random integer . First we compute , if this is not 1 then we know immediately that is not prime. If this is in fact 1, we then proceed to calculate where is an odd integer. Next we calculate , if this is we choose another random integer. If it is not then we proceed with the following "cycle". Let and continue to iterate in the following manner, , until . If none of these , then is composite. But if one does satisfy this criteria then is a probable prime.

Now for the reasoning behind the algorithm. From Fermat's Little Theorem we know that where is a prime and . Using this and a few other techniques we can see that we must have that . Now if is a prime number then is a field and thus there are only 2 square roots of 1, namely 1 and negative 1. Thus if none of the is equivalent to negative 1 modulo then we have found a third square root and can immediately conclude that is not prime.

It can be shown (not it this post) that after running the algorithm with "testers" (the random integers) which all point to being prime then the probability that is composite is . Also the algorithm has runtime of where is the number of testers. It should be noted that this can be improved to if we use the fast Fourier transform to do the multiplications.

Well below is my C++ code for the Miller-Rabin algorithm.

//a plays the role of N and b is the tester
bool MR(long long a, long long b)
{
if(gcd(a,b)!=1)
return false;
unsigned long long Y;
int B=0; //a-1=M(2^B)
long long M=a-1;
while(M%2==0)
{
B++;
M=M/2;
}
Y=RS(b,M,a); //I use repeated squaring, instead of FFT, to calculate b^M
if(Y==1 || Y==a-1)
return true;
int j;
for(j between 0 and M) //this is not real code it just did this so that it would show up in blogger
{
Y=RS(Y,2,a);
if(Y==a-1)
return true;
}
return false;
}

Yes I know this looks ugly (thanks a lot blogger) but it is still readable. Below is the code for the repeated squaring function.

unsigned long long RS(long long a, long long b, long long c) // a^b mod c
{
if(b==0)
return 1;
if(b%2==0)
return (RS(a,b/2,c)*RS(a,b/2,c))%c;
if(b%2==1)
return (a*RS(a,(b-1)/2,c)*RS(a,(b-1)/2,c))%c;
}

This algorithm calculates in time.

Well if any of you are compelled to try using this here are a few cool "easy" problems you can try.
Prime Generator
Prime or Not
(Note: Using Miller-Rabin might not be the best method to solve both of these problems)