Hello to all, i have developed a miller rabin primality test program but return me wrong result all the time.

I don't know what wrong with it after few days of debug.

Code:

#include <iostream>
#include <sstream>
#include <string>
#include <bitset>
#include <vector>
#include <limits>
#include <algorithm>
#include <functional>
#include <utility>

#include <ctime>

using namespace std;

// =========================================
/*
Miller Rabin PT
= Further enhancement by testing n using first
  7 primes(2, 3, 5, 7, 11, 13, 17, 19)
= depends on correctness of 
		Extended Riemann Hypothesis
= If N is an odd composite number then 
			the number of witnesses to the 
			compositeness of N is at least
			(N - 1) / 2.

			Example 

			1. n = 15
						(15 - 1) / 2 = 7 composite number
						4, 6, 8, 9, 10, 12, 14
			2. n  = 9
						(9 - 1) / 2 = 4 composite number
						4, 6, 8, 9
= Test 32 bit n using 2, 7, and 61 only
		Add 3 and 5 for correctness


Algorithm Background
Let p be number to test
x is a nontrivial square root of 1 mod p. Then:
X ^ 2 congruence 1 (mod p)
(x-1)(x+1) congruence 0 (mod p)

Since x is nontrivial, x is neither 
1 nor −1 mod p, and 
therefore both x−1 and x+1 are coprime to p

If p does not divide (x-1) or (x+1) 
	but product of(x-1)(x+1) then
	IsNotPrime
	X=1;
	X=-1;


Algorithm Steps
Let n be odd prime then n-1 is even rewrite
them as 2 ^ s * d(Odd)

Therefore, 

1)a ^ d = 1 (mod n) or 

2)a ^ 2r * d = -1 (mod n) for some 0 <= r <= s-1

If a is chosen uniformly at random and n is prime, 
these formulas hold with probability 1/4.

Fermat Little Theorem used 
to proove these two formula.

a ^ n-1 = 1 (mod n)

if Sqrt(a ^ n-1) must be either 1 or -1 == -1 
	2) true
else
	a ^ 2^0 * d 
	a^d != -1 mod(n) 
	1) true

If find a ^ d != 1 (mod n) or 
		a ^ 2r * d != -1 (mod n)
		a is a witness for the compositeness of n 



Example
n = 221

n-1 
= 221 - 1   
= 220
= 2 ^ s * d
= 2 ^ 2 * 55

Randomly select a < n

a = 174


Test the two formula 

a ^ d = 1 (mod 221)
a ^ 2r *d = -1 (mod 221)

174 ^ 55  
= 47 not congruence to 1 (mod 221)

174 ^ 2(1) * 55
= 174 ^ 110 (mod 221)
= 220 (n-1)


Choose a = 137
137 ^ 55 (mod 221)
= 188 not congruence to 1 (mod 221)

137 ^ 2(1) * 55 (mod 221)
= 205 not congruence to -1 (mod 221)

Therefore, 221 is not prime



Solovay–Strassen primality test




*/
// =========================================

const unsigned short NITE = 50;

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)
	{
		cout << "\n";
		Userinput(number);
		cout << boolalpha << 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);
	bool IsPrime(false);
	
	tempNumber = number - 1;
	
	if (number > 2)
	{
			// Write them as 2 ^ s * d
		while (tempNumber % 2 == 0)
		{
			// tempNumber is d
			tempNumber /= 2;
		}

		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);

			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;
}
// =========================================



// =========================================


// =========================================



// =========================================




// =========================================


// =========================================

Thanks for any comments.

Member Avatar for iamthwee

Hello to all, i have developed a miller rabin primality test program but return me wrong result all the time.

I don't know what wrong with it after few days of debug.

<code snipped>

Thanks for any comments.

This didn't compile without problems for me.

This didn't compile without problems for me.

Remove all those comments at the top, then it will compile -- at least it did with vc++ 2008 express. I don't have time this morning to debug the code. My suggestion to the OP is learn to use your compiler's debugger and find out why the program fails to work correctly. That's what you will have to do if you ever get a professional job coding.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.