Deterministic Miller-Rabin implementation

Question

I'm trying to implement a primality check function with a deterministic Miller-Rabin algorithm but the results are not always correct: when checking first 1,000,000 numbers it only founds 78,495 instead of 78,498.

This is obtained using [2, 7, 61] as a base which, according to wikipedia, should always be correct for values up to 4,759,123,141.
The interesting thing is that the 3 missing primes are exactly the ones componing the base (2, 7 and 61).

Why is this happening? The code I'm using is the following:

T modular_power(T base, T exponent, T modulo) {
    base %= modulo;
    T result = 1;

    while (exponent > 0) {
        if (exponent % 2 == 1)
            result = (result * base) % modulo;
        base = (base * base) % modulo;
        exponent /= 2;
    }

    return result;
}

bool miller_rabin(const T& n, const vector<T>& witnesses) {
    unsigned int s = 0;
    T d = n - 1;
    while (d % 2 == 0) {
        s++;
        d /= 2;
    }

    for (const auto& a : witnesses) {
        if (modular_power<T>(a, d, n) == 1)
            continue;

        bool composite = true;
        for (unsigned int r = 0; r < s; r++) {
            if (modular_power<T>(a, (T) pow(2, r) * d, n) == n - 1) {
                composite = false;
                break;
            }
        }

        if (composite)
            return false;
    }

    return true;
}

bool is_prime(const T& n) {
    if (n < 4759123141)
        return miller_rabin(n, {2, 7, 61});
    return false; // will use different base
}

Further - M-R asserts: `{n in odd | n >= 3}`. So one simple rejection in `is_prime`: `if ((n & 0x1) == 0) return (n == 2);` ... discards all evens as composite, except `(2)`. The deterministic test explicitly *passes* the case: `if ((a %= n) == 0)` (since this test tells us nothing), and further ensures that `(0 < a < n)` from this point. — Brett Hale, Jan 14 '18 at 17:16
My [deterministic test](https://stackoverflow.com/a/35825375/906839) for *all* unsigned 64-bit values. You might drop back to `uint64_t` if you're only testing 32-bit values. — Brett Hale, Jan 14 '18 at 17:21
@BrettHale I'm using this for [Project Euler](https://projecteuler.net/) problems and they require pretty big numbers so `uint64_t` might not be enough. Your tips are very good, I had already implemented the even check and will probably use a couple more of your optimisations, thanks for the help! — Becks, Jan 14 '18 at 18:17

score 3 · Accepted Answer · answered Jan 13 '18 at 17:50

Miller-Rabin indeed does not work when the base and the input are the same. What happens in that case is that a^d mod n is zero (because a mod n is zero, so this is really raising zero to some irrelevant power), and the rest of the algorithm is unable to "escape" from zero and concludes that you're dealing with a composite.

As a special case of that, Miller-Rabin never works with an input of 2, because there is no base that can be selected. 2 itself is useless, so is 1, that leaves nothing.

Deterministic Miller-Rabin implementation

1 Answers1