Solovay-Strassen Primality Test

Pages: 12
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!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#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;
	
}

Last edited on
just keep running in to problems.


And these problems are?
Sorry, I meant to type that in.

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

Thanks!
Well, at least for most prime numbers it returns not prime. For a few seemingly random primes (like 19), it actually works.
Hi Numeri, trying to get a hang of your code.

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;
	} 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.
Last edited on
Thanks, I'll change that.
Last edited on
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
1
2
3
4
5
6
7
8
9
/*	  Temp=a;
	  a=n;
	  n=Temp;*/
	  swap(a,n);
	  a%=n; //too early
	  if(((a%4)==3)&&((n%4)==3))
	  {
		  Mul=-Mul;
	  }
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!
Last edited on
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.
1
2
3
4
5
6
	  if(((a%4)==3)&&((n%4)==3))
	  {
		  Mul=-Mul;
	  }
 	  swap(a,n);
	  a%=n; //here 
Last edited on
Oops. Maybe I didn't check enough values.

Thanks for that!
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.
Last edited on
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)
		{
			return false;
		}
	}
	
	return true;
}
It still doesn't work though


What are the symptoms?
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!
me wrote:
_ Jacobi can be 0, -1 or 1.
When computing the modulus of a negative number shit happens. Make it a positive equivalent
What do you mean make it a positive equivalent? Sorry, I don't understand.

Thanks though!
Why not:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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?
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.





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, 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;
}
Last edited on
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!
Pages: 12