Solovay-Strassen Primality Test

[tition]

How does the ModPow function use doubles

My bad, thought it was part of std::

[Numeri]

The RaiseToPower function is saying that the value must be a modifiable lvaule when I use it like this:

RaiseToPower(a+0.0,(n-1)/2,n)

and when I do this:

RaiseToPower(a,(n-1)/2,n)

It says that the expression must have arithmetic, enum, pointer, or handle type.

What does that all mean?

[tition]

Numeri, I fixed the function a second time, now you should have no problems.

It says that the expression must have arithmetic, enum, pointer, or handle type

The original function I typed was void RaiseToPower(int& inputOutput, /*etc*/) This means that variable inputOutput is passed by reference, i.e. a pointer to the variable is passed, rather than the value of the variable. a+0.0 is stored in a temporary variable, so passing a reference to it was impossible.
[Edit:]Oops the above is wrong: the function I typed was returning void, which can not be compared to an integer. The new version of the function has the same semantics as your original ModPow. It should be faster though.


[Edit:] I think I saw an error (maybe that was what ne555 was referring to?):

if( x==0 || ModPow(a + 0.0,( n - 1 ) / 2, n) != x % n )

should be

/*Edit: ouch forgot this one:*/
if (x<0)
 x+=n;
if( x==0 || ModPow(a + 0.0,( n - 1 ) / 2, n) != x       /*nothing here*/)

[cire]

The Jacobi symbol, written (n/m) or (n/m) is defined for positive odd m as...

From: http://mathworld.wolfram.com/JacobiSymbol.html

Also:

The Jacobi symbol is a generalisation of the Legendre symbol to (a/n), where n can be any odd integer

from the wikipedia site you got the algorithm description from.

You might try something like the following, which hopefully delivers a little more information about why specific random numbers aren't returning the result you expect, although all of this information could be easily extracted by stepping through the code in a debugger.

bool SSPA(int n, int k)		//Solovay-Strassen Primality Algorithm
{
	std::cout << "Testing " << n << '\n' ;

	for(k;k>0;k--)
	{
		int a = (rand()%((n-3)/2))*2 + 3;  // pick odd numbers in the range 3 to n-1
		int x = Jacobi(a,n);

		if ( x < 0 )   // x modulo n for -x
			x += n ;

		int testval = ModPow(a, (n-1)/2, n) ;

		std::cout << "a = " << a << "     testval = " << testval << "     x= " << x << '\n' ;

		if ( x == 0 || testval != x )
			return false ;
	}
	
	return true;
}

[Numeri]

Using cire's code (minus the testing bits) the algorithm works! I guess the odd numbers / negative moduli issues were the problem.

[Numeri]

Oh... It stops working at around 46349, the 4793th prime. The numbers get too large for the functions. So, to make it accept long long unsigned int types (Is it fine to use long long unsigned int types?) I need a rand() that generates that large of numbers. Would this work?

long long unsigned int LRand()
{
   long long unsigned int L=rand();
   int shift=ceil(log(RAND_MAX+0.0)/log(2.0));
   for(int i=0;i<(64/(shift+1));i++)
   {
   L<<shift;
   L+=rand();
   }
   return L;
}

EDIT: Forgot to return a value.

[Numeri]

Yep, it works!
It even returns that the 1.0e8th prime (2,038,074,743) is prime!

Thank you so much everyone for your help!