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.
