implement and verify the Miller-Rabin primality test

millerrabin.c

#include <stdio.h>
#include <stdlib.h>

int modular_pow(int base, int exponent, int modulus)
{
	if(modulus == 1) return 0;
	int result = 1;
	base %= modulus;
	while(exponent > 0){
		if(exponent & 1){
			result = (result * base) % modulus;
		}
		exponent >>= 1;
		base = (base * base) % modulus;
	}
	return result;
}

int miller_rabin_test(int n)
{
	const int bases[] = {2, 7, 61};
	int r = 0, d = n - 1, i, j;
	
	if(n == 2 || n == 3) return 1;
	if((n & 1) == 0 || n < 2) return 0;
	
	while(!(d & 1)){
		d >>= 1;
		r++;
	}
	for(j = 0; j < sizeof(bases) / sizeof(bases[0]); j++){
		int a = bases[j];
		int x = modular_pow(a, d, n);
		if(x == 1 || x == n - 1) continue;
		for(i = 1; i < r; i++){
			x = modular_pow(x, 2, n);
			if(x == n - 1) break;
		}
		if(i == r) return 0;
	}
	return 1;
}

int main()
{
	int i, j;
	for(i = 0; i < 2000000000; i++){
		if(miller_rabin_test(i)){
			for(j = 5; j * j <= i; j++){
				if(i % j == 0) break;
			}
			if(j * j <= i) printf("Fail for i=%d\n", i);
		}
	}
	return 0;
}