Miller-Rabin Primality test FIPS 186-3 implementation

Question

Im trying to implement the Miller-Rabin primality test according to the description in FIPS 186-3 C.3.1. No matter what I do, I cannot get it to work. The instructions are pretty specific, and I dont think I missed anything, and yet Im getting true for non-prime values.

What did I do wrong?

template <typename R, typename S, typename T>
T POW(R base, S exponent, const T mod){
    T result = 1;
    while (exponent){
        if (exponent & 1)
            result = (result * base) % mod;
        exponent >>= 1;
        base = (base * base) % mod;
    }
    return result;
}



// used uint64_t to prevent overflow, but only testing with small numbers for now
bool MillerRabin_FIPS186(uint64_t w, unsigned int iterations = 50){
    srand(time(0));
    unsigned int a = 0;
    uint64_t W = w - 1; // dont want to keep calculating w - 1
    uint64_t m = W;
    while (!(m & 1)){
        m >>= 1;
        a++;
    }

    // skipped getting wlen
    // when i had this function using my custom arbitrary precision integer class,
    // and could get len(w), getting it and using it in an actual RBG 
    // made no difference 

    for(unsigned int i = 0; i < iterations; i++){
        uint64_t b = (rand() % (W - 3)) + 2; // 2 <= b <= w - 2
        uint64_t z = POW(b, m, w);
        if ((z == 1) || (z == W))
            continue;
        else
            for(unsigned int j = 1; j < a; j++){
                z = POW(z, 2, w);
                if (z == W)
                    continue;
                if (z == 1)
                    return 0;// Composite
            }
    }
    return 1;// Probably Prime
}

this:

std::cout << MillerRabin_FIPS186(33) << std::endl;
std::cout << MillerRabin_FIPS186(35) << std::endl;
std::cout << MillerRabin_FIPS186(37) << std::endl;
std::cout << MillerRabin_FIPS186(39) << std::endl;
std::cout << MillerRabin_FIPS186(45) << std::endl;
std::cout << MillerRabin_FIPS186(49) << std::endl;

is giving me:

0
1
1
1
0
1

Principled

Can we see POW? I see a mistake that could declare some primes as composite, but nothing jumps out for the other way round. For what values are you getting wrong results? – Trainor

Where is your definition of POW? – Monodic

pow is fine, but ill put it up in any case – Principled

BTW, your random numbers are not uniformly distributed -- that modulo-operation skews the results. – Anyone

What happens if result * base or base * base overflows? – Anyone

they shouldnt, since the values are so small. also, normally, i would be using a custom arbitrary precision integer, so they wouldnt overflow. i just changed it to uint64_t to double check that the math in the custom class wasnt wrong – Principled

Answer 1

The only difference between your implementation and Wikipedia's is that you forgot the second return composite statement. You should have a return 0 at the end of the loop.

Edit: As pointed out by Daniel, there is a second difference. The continue is continuing the inner loop, rather than the outer loop like it's supposed to.

for(unsigned int i = 0; i < iterations; i++){
    uint64_t b = (rand() % (W - 3)) + 2; // 2 <= b <= w - 2
    uint64_t z = POW(b, m, w);
    if ((z == 1) || (z == W))
        continue;
    else{
        int continueOuter = 0;
        for(unsigned int j = 1; j < a; j++){
            z = POW(z, 2, w);
            if (z == W)
                continueOuter = 1;
                break;
            if (z == 1)
                return 0;// Composite
        }
        if (continueOuter) {continue;}
    }
    return 0; //This is the line you're missing.
}
return 1;// Probably Prime

Also, if the input is even, it will always return probably prime since a is 0. You should add an extra check at the start for that.

Monodic

One hopes a primality test isn't being used on even numbers. :) – Anyone

Oh c'mon, this certainly doesn't deserve a downvote... it's a good point. :) – Anyone

Why the downvote? It's a legitimate issue and I had no idea what numbers were being tested at the time I wrote that. – Monodic

now im getting 0 0 0 0 0 0 and i added the {} to the else, so its not returning 0 after 1 iteration – Principled

@calccrypto: you can't just add the break and the return 0, because now whenever you break the return 0 is triggered. Set a flag (say contloop) to false before the loop. Instead of just break, set "contloop = 1; break;}, and then only return 0 if contloop is false. [Or something.] – Durfee

Answer 2

In the inner loop,

        for(unsigned int j = 1; j < a; j++){
            z = POW(z, 2, w);
            if (z == W)
                continue;
            if (z == 1)
                return 0;// Composite
        }

you should break; instead of continue; when z == W. By continueing, in the next iteration of that loop, if there is one, z will become 1 and the candidate is possibly wrongly declared composite. Here, that happens for 17, 41, 73, 89 and 97 among the primes less than 100.

Trainor

Aargh, just as I was about to hit send, too. I think both this and the return 0 if this loop makes it all the way through are necessary. – Durfee

Wow, I can't believe I missed that. The continue is only continuing the inner loop, not the outer loop like it's supposed to. – Monodic