Deterministic MillerRabin implementation
I'm trying to implement a primality check function with a deterministic MillerRabin 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
}
1 answer

MillerRabin 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, MillerRabin 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.