Miller-Rabin primality test issue in C++

Question

I've been trying to implement the algorithm from wikipedia and while it's never outputting composite numbers as primes, it's outputting like 75% of primes as composites.

Up to 1000 it gives me this output for primes:

3, 5, 7, 11, 13, 17, 41, 97, 193, 257, 641, 769

As far as I know, my implementation is EXACTLY the same as the pseudo-code algorithm. I've debugged it line by line and it produced all of the expected variable values (I was following along with my calculator). Here's my function:

bool primeTest(int n)
{
    int s = 0;
    int d = n - 1;

    while (d % 2 == 0)
    {
        d /= 2;
        s++;
    }

    // this is the LOOP from the pseudo-algorithm
    for (int i = 0; i < 10; i++)
    {
        int range = n - 4;
        int a = rand() % range + 2;
        //int a = rand() % (n/2 - 2) + 2;
        bool skip = false;
        long x = long(pow(a, d)) % n;

        if (x == 1 || x == n - 1)
            continue;

        for (int r = 1; r < s; r++)
        {
            x = long(pow(x, 2)) % n;

            if (x == 1)
            {
                // is not prime
                return false;
            }
            else if (x == n - 1)
            {
                skip = true;
                break;
            }
        }

        if (!skip)
        {
            // is not prime
            return false;
        }
    }

    // is prime
    return true;
}

Any help would be appreciated D:

EDIT: Here's the entire program, edited as you guys suggested - and now the output is even more broken:

bool primeTest(int n);

int main()
{
    int count = 1;     // number of found primes, 2 being the first of course
    int maxCount = 10001;
    long n = 3;
    long maxN = 1000;
    long prime = 0;

    while (count < maxCount && n <= maxN)
    {
        if (primeTest(n))
        {
            prime = n;
            cout << prime << endl;
            count++;
        }

        n += 2;
    }

    //cout << prime;
    return 0;
}

bool primeTest(int n)
{
    int s = 0;
    int d = n - 1;

    while (d % 2 == 0)
    {
        d /= 2;
        s++;
    }

    for (int i = 0; i < 10; i++)
    {
        int range = n - 4;
        int a = rand() % range + 2;
        //int a = rand() % (n/2 - 2) + 2;
        bool skip = false;
        //long x = long(pow(a, d)) % n;
        long x = a;
        for (int z = 1; z < d; z++)
        {
            x *= x;
        }

        x = x % n;

        if (x == 1 || x == n - 1)
            continue;

        for (int r = 1; r < s; r++)
        {
            //x = long(pow(x, 2)) % n;
            x = (x * x) % n;

            if (x == 1)
            {
                return false;
            }
            else if (x == n - 1)
            {
                skip = true;
                break;
            }
        }

        if (!skip)
        {
            return false;
        }
    }

    return true;
}

Now the output of primes, from 3 to 1000 (as before), is:

3, 5, 17, 257

I see now that x gets too big and it just turns into a garbage value, but I wasn't seeing that until I removed the "% n" part.

Answer 1

The likely source of error is the two calls to the pow function. The intermediate results will be huge (especially for the first call) and will probably overflow, causing the error. You should look at the modular exponentiation topic at Wikipedia.

Answer 2

Source of problem is probably here:

x = long(pow(x, 2)) % n;

pow from C standard library works on floating point numbers, so using it is a very bad idea if you just want to compute powers modulo n. Solution is really simple, just square the number by hand:

x = (x * x) % n