Miller Rabin Primality Test Wrong Result

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 ?

Code:

---

#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.

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...

Code:

---

               
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?

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.

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.

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?