I am working on the Solovay-Strassen primality test and just keep running in to problems. If you don't know what it is, look at http://en.wikipedia.org/wiki/Solovay-Strassen_primality_test.
The algorithm is on the wikipedia page in pseudocode, and I am fairly certain I got it to C++ accurately. I also found a working algorithm for the Jacobi symbol. I have checked through the code and made sure that the Jacobi symbol function, GCD function, and the Modular Power function are correct. I can't find anything wrong with the Solovay-Strassen part either. I am using Visual Studio 2010 on Windows 7. Here is the code, thanks!
First remark: you shouldn't be calling gcd with argument zero.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
int GCD(int a, int b)
{
if( a == 1 && b == 0 )
{
return 1;
} elseint Intermediate;
while(b!=0)
{
Intermediate=a;
a=b;
b=Intermediate % a;
}
return a;
}
[Edit:] ne555 has a valid point: Do your gcd like that. If the user calls gcd with a=1 and b=0, he/she *should* get division by zero. You could also throw an exception or do some other error handling, but make sure it displays a fat error message.
@ne555: you are right I don't check whether b is zero, but gcd is not defined for
b==0. I guess I don't raise alarm just to make the code simpler. But I should.
static int gcd(int a, int b)
{ int temp;
while(b>0)
{ temp= a % b;
a=b;
b=temp;
}
return a;
}
I thought about it awhile, and I am not sure what I crossed has merit or not. They say, think first then speak, but nevermind.
Alright, now Jacobi(3,7) returns -1, as does Jacobi(2,3).
It still doesn't work though, so is it my SSPA() that is wrong now? It is the actual algorithm. Wikipedia says that it is this in pseudocode:
1 2 3 4 5 6 7
Inputs: n, a value to test for primality; k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
repeat k times:
choose a randomly in the range [2,n − 1]
x ← Jacobi(a,n)
if (x = 0) or (a^( (n-1)/2 ) is not congruent to x mod n) then return composite
return probably prime
In C++ I did:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
bool SSPA(int n, int k) //Solovay-Strassen Primality Algorithm
{
int a;
int x;
for(k;k>0;k--)
{
a=rand()%(n-3)+2;
x=Jacobi2(a,n);
if(x==0||ModPow(a+0.0,(n-1)/2,n)!=x%n)
{
returnfalse;
}
}
returntrue;
}
bool SSPA(int n, int k) //Solovay-Strassen Primality Algorithm
{
int a;
int x;
for(k;k>0;k--)
{
a=rand()%(n-3)+2;
x=Jacobi(a,n);
if( x==0 || ModPow(a + 0.0,( n - 1 ) / 2, n) != x % n )
{
std::cout << "The evil values that fail the test are: " << a << " and " << n;
returnfalse;
}
}
returntrue;
}
and then let us test the bad values that fail your test?
At first it seemed like it was primes that were bad, but then there are composites later, but the primes still occurred more frequently than the composites, and the numbers 5 and 7 occurred most frequently.
You are using a ModPow function which appears to deal with doubles. Can you print the output of ModPow as well? It is possible that ModPow gives you some nasty approximation number, and x%n is an integer converted to double, so the comparison of an approximation to a non-approximation kills you.
Also, try printing out all relevant variables to the test-fail in human-readable form.
Also, if you want to raise to power using integers only, you might wanna try something like this:
[Edit:]sorry gotta fix the code, one moment please...
[Edit:]One more edit for clarity:
...Think the code is fine now, didn't really test:
1 2 3 4 5 6 7 8 9 10 11 12 13
int RaiseToPower(int input, unsignedint thePower, int theModulus)
{ if (thePower==0)
return 1;
int squares=input%theModulus;
input=1;
while (thePower>0)
{ if (thePower%2==1)
input=(squares*input)%theModulus;
squares=(squares*squares)%theModulus;
thePower/=2;
}
return input;
}
How does the ModPow function use doubles? There isn't a single double declared in it.
For the test-fail, I meant N is the number that is entered as a prime, and the numbers below it are the bad numbers it returned. I formatted it that way because it was more readable to me... :D