Miller Rabin Primality Test Wrong Result

[Peter_APIIT]

Hello to all, i try to code miller rabin algorithm with C++ but the program doesn't did what it suppose to be.

What wrong with my program ?

#include <bitset>
#include <vector>
#include <limits>
#include <algorithm>
#include <functional>
#include <utility>

#include <ctime>

using namespace std;

const unsigned short NITE = 6;

typedef unsigned long long ulong;

void Userinput(ulong&);
bool MillerRabinPrimeTest(const ulong&);
ulong ComputeModularExponentiation(const ulong&, const ulong&, const ulong&);
ulong ComputeExponentiation(const ulong&, const ulong&);

// =========================================
int main()
{
	ulong number(0);

	while (1)
	{
		Userinput(number);
		cout << MillerRabinPrimeTest(number);
	}
	
	return 0;
}
// =========================================
void Userinput(ulong& number)
{
	cout << "Enter Number : ";
	cin >> number;

	while (cin.fail())
	{
		cin.clear();
		cin.ignore(numeric_limits<int>::max(), '\n');

		cout << "\nEnter Number : ";
		cin >> number;
	}
}
// =========================================
bool MillerRabinPrimeTest(const ulong& number)
{
	ulong a(0), x(0), y(0), tempNumber(0), s(1);
	vector<int> primeFactor(10);
	bool IsPrime(false);

	tempNumber = number - 1;

	primeFactor[0] = 2;
	primeFactor[1] = 3;
	primeFactor[2] = 5;
	primeFactor[3] = 7;
	primeFactor[4] = 11;
	primeFactor[5] = 13;
	primeFactor[6] = 19;
	primeFactor[7] = 23;
	primeFactor[8] = 29;
	primeFactor[9] = 31;
	
	if (number > 2)
	{
			// Write them as 2 ^ s * d
		while (tempNumber % 2 == 0)
		{
			// tempNumber is d
			tempNumber /= 2;
		}

		// compute s 
		while (ComputeExponentiation(2, s) * tempNumber != number - 1)
		{
			++s;
		}

		for (int loop = 0;loop < NITE;++loop)
		{
			srand((unsigned int)time(0));
			// rand() % upperBound + startnumber
			 a = rand() % (number - 2) + 2;
		//	a = primeFactor[loop];
			x = ComputeModularExponentiation(a, tempNumber, number);

			for (int r = 0; r < s-1;++r)
			{
				y = (x * x) % number;
				if (x == 1 || y == -1)
				{
					return true;
				}
			
			}

		}
	}

	return IsPrime;
}
// =========================================
ulong ComputeModularExponentiation(
			const ulong& a, const ulong& d, 
			const ulong& n)
{
	enum {NBITS = numeric_limits<ulong>::digits };
	string bitExponent = bitset<NBITS>((unsigned long)d).to_string(); 
	typedef string::size_type strType;
	strType MSSB = bitExponent.find_first_of('1'); 
	ulong result(a % n);

	MSSB += 1;

	while (MSSB < NBITS)
	{
		result *= result;

		if (bitExponent[MSSB] == '1')
		{
			result = result * a % n; 
		}
		++MSSB;
	}

	return result;
}
// =========================================
ulong ComputeExponentiation(const ulong& base, const ulong& exponent)
{
	enum {NBITS = numeric_limits<ulong>::digits };
	string bitExponent = bitset<NBITS>((unsigned long)exponent).to_string(); 
	typedef string::size_type strType;
	strType MSSB = bitExponent.find_first_of('1'); 
	ulong result(base);

	MSSB += 1;

	while (MSSB < NBITS)
	{
		result *= result;

		if (bitExponent[MSSB] == '1')
		{
			result *= base; 
		}
		++MSSB;
	}

	return result;
}

Please explain where my program goes wrong.

Thanks.

[NMTop40]

I don't know the test but there are a few things funny about your code sample. I guess some of the one-letter identifiers may be "standard" within the mathematical realm of that formula, however...

		
for (int r = 0; r < s-1;++r)
{
    y = (x * x) % number;
    if (x == 1 || y == -1)
    {
        return true;
    }
}

x and number never change during this loop so y will be the same each time it loops so it will either be -1 or not and I don't get the point of looping.

Just use a simple if.

In addition this looks like it could be never terminating:

What exactly is ComputeExponentation and the modular version doing?

[kempofighter]

It is very helpful when the OP tells us what the program is supposed to do vs. what it is actually doing. I have never heard of the miller rabin algorithm and don't have any idea what the expected output should be. Give us a sample of the expected output vs. the actual output so that we don't have to spend a bunch of time debugging the program to learn. Many people simply won't do that and rightly so.

[Peter_APIIT]

What exactly is ComputeExponentation and the modular version doing?

The ComputeExponentation raise a base number to power of exponent and the modulus version does % after the compute exponentiation.


Sample Input
2
3
5
6
7

Sample Output
Prime
Prime
Prime
Composite
Prime

Thanks.

[NMTop40]

One guess for a start is that MSSB is reading from the wrong side, but do not convert to a string and use find_first_of, why not just use std::bitset's in-built functions?