Implementing Miller-Rabin in C

Question

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?

Possibly because your `modpow()` overflows computing `a^b` *before* reducing modulo-`m`. You'll need to reduce `a` and `c` as you go. If you write some tests for `modpow()`, you'll have greater confidence in it. — Toby Speight, Mar 21 '17 at 16:59

score 6 · Accepted Answer · answered Mar 21 '17 at 18:18

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;
}

score 1 · Answer 2 · answered Mar 22 '17 at 10:44

There are problems in your modpow() function. You probably want to use unsigned types for the arguments and the result (what would a negative m mean, anyway?) Secondly, it overflows because it tries to compute a^b before reducing it modulo-m. You need to reduce a and c as you go.

The best way to deal with this is to write some tests:

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

#include <stdio.h>
unsigned test(unsigned a, unsigned b, unsigned m, unsigned expected)
{
    unsigned actual = modpow(a, b, m);
    if (actual == expected)
        return 0;
    printf("modpow(%u, %u, %u) returned %u; expected %u\n",
           a, b, m, actual, expected);
    return 1;
}


int main()
{
    return test(0, 0, 2, 1)
        +  test(1005, 16, 100, 25)
        ;
}

The first test passes (but raises a question - what result do you want when m < 2?); the second one fails:

modpow(1005, 16, 100) returned 49; expected 25

Let's modify modpow() to reduce the result at each step:

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

It now passes! We can make another failing test:

int main()
{
    return test(0, 0, 2, 1)
        +  test(1005, 16, 100, 25)
        +  test(100000005, 16, 1000000000, 587890625)
        ;
}

modpow(100000005, 16, 1000000000) returned 919214529; expected 587890625

Now we need to calculate the multiplication using a larger type:

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

Once we have sufficient confidence in the modpow() function, we can debug the rest of the algorithm.

Note that if your integer sizes are different to mine, you'll need to use different values in the tests to replicate the results. I chose large numbers ending in 005 because we know the last two digits are invariant 25 regardless of the power. You might find that dc -e '???|p' helps in generating test cases (supply the three arguments on its stdin, and it will print the expected value).

Implementing Miller-Rabin in C

2 Answers2