calgorithmprimesprimality-test

Implementing Miller-Rabin in C


I'm trying to implement the Miller-Rabin primality test in C99, but I'm coming across some problems getting it to work. I crafted a small test-set to verify whether or not the implementation works, here's how I'm checking for primes

int main() {
    int foo[11] = {0, 1, 2, 3, 4, 7, 28, 73, 125, 991, 1000};
    for (int i = 0; i < 11; i++) {
        printf("%s; ", isprime(foo[i], 5000) ? "Yes" : "No");
    }
    return 0;
}

From the numbers listed, the expected output would be

No; No; Yes; Yes; No; Yes; No; Yes; No; Yes; No;

However, as implemented , the output I get is the following:

No; No; Yes; Yes; No; Yes; No; No; No; No; No;

Here's how I wrote the algorithm

int randint (int low, int up){
    return rand() % (++up - low)+low;
}

int modpow(int a, int b, int m) {
    int c = 1;
    while (b) {
        if (b & 1) {
            c *= a;
        }
        b >>= 1;
        a *= a;
    }
    return c % m;
}

bool witness(int a, int s, int d, int n) {
    int x = modpow(a,d,n);
    if(x == 1) return true;
    for(int i = 0; i< s-1; i++){
        if(x == n-1) return true;
        x = modpow(x,2,n);
    }
    return (x == n- 1);
}

bool isprime(int x, int j) {
    if (x == 2) {
        return true;
    }
    if (!(x & 1) || x <= 1) {
        return false;
    }
    int a = 0;
    int s = 0;
    int d = x - 1;

    while (!d&1){
        d >>=1;
        s+=1;
    }
    for(int i = 0; i < j; i++){
        a = randint(2, x-1);
        if(!witness(a,s,d,x)){
            return false;
        }
    }

    return true;
}

What am I doing wrong? Why is the test failing for "large" primes, but working for very small ones? How may I fix this?


Solution

  • With Visual Studio 2015 Community edition I found two problems. First the line:

    while (!d&1){
    

    needs to be:

    while (!(d&1)){
    

    Secondly, as mentioned in the comments your modpow function is overflowing. Try:

    int modpow(int a, int d, int m) {
        int c = a;
        for (int i = 1; i < d; i++)
            c = (c*a) % m;
        return c % m;
    }