【c语言随笔4】数论之概率性素数性测试(Miller-Robin Algorithms)

算法思想： 利用素数的性质
若n为大于2的素数,则有 n-1 = (2^k)*q, k>0, q为奇数, 且 1<a<(n-1)
则满足以下二者之一的条件:
1. a^q mod n = 1;
2. 存在j在[1, k]区间内的正整数, a^{(2^(j-1))*q} mod n = n-1
若满足二者条件之一未必为素数
若二者都不满足，必为合数

注意：
本文不作证明，可以参考其他博主博客；
当输入值过大，代码中可能存在数值上溢情况，可将int型换为long 或 long long型，请自行酌情分析。

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

// 快速模幂指数算法：奇数减一，偶数除二； --> 原理：模运算的乘法性质
int quickModPow(int base, int power, int n){
	int r; // result

	r = 1;
	while(power > 0){
		if(power & 1){ // 等价于power%2 == 1
			r = r * base % n;
		}
		power >>= 1; // power/=2
		base = (base * base) % n;
	}

	// return result
	return r;
}

// Miller-Robin Algorithms
int TEST(int n){ // n为要测试素性的数
	// 算法思想： 利用素数的性质
	// 若n为大于2的素数,则有 n-1 = (2^k)*q, k>0, q为奇数, 且 1<a<(n-1)
	// 则满足以下二者之一的条件
	// 1. a^q mod n = 1;
	// 2. 存在j在[1, k]区间内的正整数, a^((2^(j-1))*q) mod n = n-1
	// 若满足二者条件之一未必为素数
	// 若二者都不满足，必为合数
	int k, q, a; 
	int i, t;

	if((!(n%2) && n != 2) || n <= 0) return 0;
	else if(n == 2 || n == 1) return 1;

	/* n-1 = q*(2^k), q is odd number */
	q = (n-1)/2;
	k = 1;
	while(!(q%2)){
		q = q/2;
		k++;
	}
	// printf("k: %d, q: %d\n", k, q);

	// a random in {2,..., n-2}, 1<a<n-1
	a = 0;
	while(!a || a == 1){
		a = rand()%(n-1);
	}
	// printf("a: %d\n", a);

	// possibility 1
	t = quickModPow(a, q, n);
	if(t == 1) return 1;

	// possibility 2
	for (i = 0; i < k; ++i)
	{	
		t = pow(2, i);
		t = t * q;
		t = quickModPow(a, t, n);
		if (t == n-1)
		{
			return 1;
		}
	}

	// n is composite number
	return 0;
}


int main(int argc, char const *argv[])
{
	int i;

	// 测试10次，失败概率为10E-6
	for (i = 0; i < 10; ++i)
	{
		// 测试素数1999
		if(!TEST(1999)){
			printf("not prime!\n");
			break;
		}
	}
	if (i == 10)
	{
		printf("prime!\n");
	}

	return 0;
}

额外知识补充（素数分布）： 任意整数n附近，平均每 ln(n) 个整数中会有一个素数，由于素数必不是偶数(2除外)，故要寻找n附近的一个素数大概需要 0.5ln(n) 次。