Miller-Rabin deterministic primality test (C)

Question

I have found a program in C that implemented a deterministic variant of the Miller-Rabin primality test here. However, the modified code (which can be seen below) doesn't work when dealing with numbers bigger than 2^32, even though I use unsigned long long data type to store my numbers. Which should be able to hold integers up to 2^64. Where does it go wrong?

In short: my problem is that my code correctly determines if a number is prime or not but only if it is smaller than 2^32, which shouldn't be the case since I can store numbers up to 2^64

unsigned long long power(unsigned long long a, unsigned long long n, unsigned long long mod)
{
    unsigned long long power = a;
    unsigned long long result = 1;

    while (n)
    {
        if (n & 1)
            result = (result * power) % mod;
        power = (power * power) % mod;
        n >>= 1;
    }
    return result;
}

// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1
unsigned long long witness(unsigned long long n, unsigned long long s, unsigned long long d, unsigned long long a)
{
    unsigned long long x = power(a, d, n);
    unsigned long long y;

    while (s) {
        y = (x * x) % n;
        if (y == 1 && x != 1 && x != n-1)
            return 0;
        x = y;
        --s;
    }
    if (y != 1)
        return 0;
    return 1;
}

/*
 * if n < 1,373,653, it is enough to test a = 2 and 3;
 * if n < 9,080,191, it is enough to test a = 31 and 73;
 * if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
 * if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
 * if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
 * if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
 * if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
 */

int is_prime_mr(unsigned long long n)
{
    if (((!(n & 1)) && n != 2 ) || (n < 2) || (n % 3 == 0 && n != 3))
        return 0;
    if (n <= 3)
        return 1;

    unsigned long long d = n / 2;
    unsigned long long s = 1;
    while (!(d & 1)) {
        d /= 2;
        ++s;
    }
    unsigned long long a1 = 4759123141;
    unsigned long long a2 = 1122004669633;
    unsigned long long a3 = 2152302898747;
    unsigned long long a4 = 3474749660383;
    if (n < 1373653)
        return witness(n, s, d, 2) && witness(n, s, d, 3);
    if (n < 9080191)
        return witness(n, s, d, 31) && witness(n, s, d, 73);
    if (n < a1)
        return witness(n, s, d, 2) && witness(n, s, d, 7) && witness(n, s, d, 61);
    if (n < a2)
        return witness(n, s, d, 2) && witness(n, s, d, 13) && witness(n, s, d, 23) && witness(n, s, d, 1662803);
    if (n < a3)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11);
    if (n < a4)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13);
    return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);
}

int main()
{
  unsigned long long in1 = 4294967291;
  unsigned long long in2 = 4294967311;
  int a = is_prime_mr(in1); //4294967291 is the last prime below 2^32, should return 1 and does so correctly
  printf("%d\n",a);
  int b = is_prime_mr(in2); //4294967311 is the first prime after 2^32, should also return 1 but returns 0
  printf("%d",b);
  return 0;
}

Answer 1

As @Michael Dorgan suggested, the problem lay in the multiplication of two 64-bit integers which was resolved using a different implementation of powermod.

//computes (a*b) (mod n)
unsigned long long mul_mod(unsigned long long a, unsigned long long b, unsigned long long m)
{
   unsigned long long d = 0, mp2 = m >> 1;
   if (a >= m) a %= m;
   if (b >= m) b %= m;
   for (int i = 0; i < 64; ++i)
   {
       d = (d > mp2) ? (d << 1) - m : d << 1;
       if (a & 0x8000000000000000ULL)
           d += b;
       if (d >= m) d -= m;
       a <<= 1;
   }
   return d;
}
//computes (a^b) (mod m)
unsigned long long powermod(unsigned long long a, unsigned long long b, unsigned long long m)
{
    unsigned long long r = m==1?0:1;
    while (b > 0) {
        if (b & 1)
            r = mul_mod(r, a, m);
        b = b >> 1;
        a = mul_mod(a, a, m);
    }
    return r;
}

My code now works for integers up to 2^64.

Answer 2

Here's a simple, deterministic implementation of the Miller-Rabin primality test in C, designed for 64-bit unsigned integers. This implementation avoids overflows, has no external dependencies, and is compliant with C99 standards.

Implementation

typedef unsigned long long int ulong;

ulong mul_mod(ulong a, ulong b, const ulong mod) {
    ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
    for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
    return res % mod;
}

ulong pow_mod(ulong n, ulong exp, const ulong mod) {
    ulong res = 1; // return (n ^ exp) % mod.
    for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
    return res;
}

int is_prime(ulong n) {
    // Perform a Miller-Rabin test, it should be a deterministic version.
    static const ulong bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    static const int n_bases = (int) sizeof(*bases);
    for (int i = 0; i < n_bases; ++i)
        if (n % bases[i] == 0)
            return n == bases[i];
    if (n < bases[n_bases - 1] * bases[n_bases - 1])
        return 1 < n;
    // Depending on the size of the number, we don't need to test all the bases.
    int lim = n < 2152302898747 ? n < 25326001 ? n < 2047 ? 1 : n < 1373653 ? 2 : 3 : n < 3215031751 ? 4 :
        5 : n < 341550071728321 ? n < 3474749660383 ? 6 : 7 : n < 3825123056546413051 ? 9 : 12, res = 1, a = 0;
    ulong b, c;
    for (b = n - 1; ~b & 1; b >>= 1, ++a);
    for (int i = 0; i < lim && res; ++i)
        if (c = pow_mod(bases[i], b, n), c != 1) {
            for (int d = a; d-- && (res = c + 1 != n);)
                c = mul_mod(c, c, n);
            res = !res;
        }
    return res;
}

It's straightforward, efficient, and free from external dependencies.

Source: github.com/michel-leonard/C-MathSnip/blob/main/miller_rabin.c