Solovay-Strassen Primality Test

[Numeri]

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!

#include "stdafx.h"
#include <time.h>
#include <iostream>
#include <conio.h>
#include <stdlib.h>

using namespace std;

int ModPow(int base, int exp, int mod)
{
	int Ans=1;
	for(exp;exp>=1;exp--)
	{
		Ans*=base;
		Ans%=mod;
	}
	return Ans;
}

int GCD(int a, int b)
{
	if( a == 1 && b == 0 )
	{
		return 1;
	}
	else
	{
		int Intermediate;
		while(b!=0)
		{
			Intermediate=a;
			a=b;
			b=Intermediate % a;
		}
		return a;
	}
	return a;
}



int Jacobi(int a, int n)
{
  int Mul=1, Temp;
  if(n%2==0)
	  return 2;
  if(GCD(a,n)!=1)
	  return 0;
  while(a!=0)
  {
	  while(a%2==0)
	  {
		  a/=2;
		  if(n%8==3 || n%8==5)
		  {
			  Mul=-Mul;
		  }  
	  }
	  Temp=a;
	  a=n;
	  n=Temp;
	  a%=n;
	  if(((a%4)==3)&&((n%4)==3))
	  {
		  Mul=-Mul;
	  }
	}
  if(n==1)
	  return Mul;
  else if(n==0)
	  return 0;
}

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 )
		{
			return false;
		}
	}
	
	return true;
}

int _tmain(int argc, _TCHAR* argv[])
{
	srand(time(NULL));
	int N,K;
Start:
	cin>>N>>K;
	cout<<"\n\n";
	if(SSPA(N,K))
	{
		cout<<"Probably Prime\n\n\n";
	}
	else
        {
		cout<<"Not Prime\n";
        }
	
		cout<<"\n\n\n";
		goto Start;
	

		
	
	_getche();
	return 0;
	
}

[cire]

just keep running in to problems.

And these problems are?

[Numeri]

Sorry, I meant to type that in.

No matter what number you input, prime or not, it returns not prime.

Thanks!

[Numeri]

Well, at least for most prime numbers it returns not prime. For a few seemingly random primes (like 19), it actually works.

[tition]

Hi Numeri, trying to get a hang of your code.

First remark:
you shouldn't be calling gcd with argument zero.

int GCD(int a, int b)
{
	if( a == 1 && b == 0 )
	{
		return 1;
	} else
	int 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.

[Numeri]

Thanks, I'll change that.

[ne555]

If the user calls gcd with a=1 and b=0, he/she *should* get division by zero.

He should not. And you are not doing that either.

@OP:
_ std::swap()

_ Jacobi can be 0, -1 or 1.
When computing the modulus of a negative number shit happens. Make it a positive equivalent

_ Your Jacobi is incorrect

/*	  Temp=a;
	  a=n;
	  n=Temp;*/
	  swap(a,n);
	  a%=n; //too early
	  if(((a%4)==3)&&((n%4)==3))
	  {
		  Mul=-Mul;
	  }

[Numeri]

I tested the Jacobi and it was returning the correct answer. Where should I put the a%=n;? The std::swap() is very handy.

Thanks!

[ne555]

It is not

$ Jacobi(3,7)
your program >> 1
correct >> -1

As it is, you are checking with a modified value. The check should be done with the originals.

	  if(((a%4)==3)&&((n%4)==3))
	  {
		  Mul=-Mul;
	  }
 	  swap(a,n);
	  a%=n; //here 

[Numeri]

Oops. Maybe I didn't check enough values.

Thanks for that!

[tition]

ne555 is right: your Jacobi symbol is buggy.

Jacobi(2,3)=-1,
your program returns 1.

@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.

[Numeri]

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:

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:

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)
		{
			return false;
		}
	}
	
	return true;
}

[tition]

It still doesn't work though

What are the symptoms?

[Numeri]

Same as before, for all inputs of N and K, it returns not prime. (once I saw it return prime for some primes but not all, but not anymore)

Thanks!

[ne555]

me wrote:

_ Jacobi can be 0, -1 or 1.
When computing the modulus of a negative number shit happens. Make it a positive equivalent

[Numeri]

What do you mean make it a positive equivalent? Sorry, I don't understand.

Thanks though!

[tition]

Why not:

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;
			return false;
		}
	}
	
	return true;
}

and then let us test the bad values that fail your test?

[Numeri]

N = 5:
3,5
2,5.

N=7:
5,7
3,7

N=11:
7,11
8,11
6,11
2,11

N=13:
5,13
7,13
11,13
2,13
6,13
8,13

N=17:
5,17
11,17
7,17
10,17
6,17
3,17

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.

It seems pretty random.

[tition]

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:

int RaiseToPower(int input, unsigned int 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;
}

[Numeri]

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

I will try your code, Thanks!