Miller-Rabin算法&乘法逆元

这周是另一个素性测试的算法，以及密码学中常用的求乘法逆元的算法、

素数的两个性质

任意素数n可表示为n = 2^k+1，k>=0，q为奇数
- 特殊情况：n = 2时， 2 = 2⁰+1
n是素数，a是小于n的正整数，则a² mod n = 1 当且仅当a mod n = 1或 a mod n = n-1
- 充分性：a² mod n = (a mod n) * (a mod n) = 1
- 必要性：a² mod n =1 , 则a² -1 mod n = 0;
- 即(a+1)(a-1) mod n =0
- ∵ n为素数，∴ (a+1) mod n =0 或 (a-1) mod n =0;
- 即得出a mod n = 1或 a mod n = n-1

Miller-Rabin算法

n = 2^k + 1 (k > 0 , q为奇数) 是大于2的素数，a是大于1且小于n-1的整数，则如下两个条件之一成立

a^q mod n = 1
存在j (j>=1 && j<=k)，满足a^{2^j-1q} mod n = n-1

原理

费马小定理

a^n-1 mod n = a ^{2^kq} mod n = 1

序列：a^q mod n，a^2q mod n，a^4q mod n，......，a ^{2^k-1q} mod n，a ^{2^kq} mod n

该序列有几个特点：

后一项恰为前一项的平方
最后一项为1

所以，可以得出Miller-Rabin算法的结论；

该序列所有项均为1
- 则a^q mod n = 1
序列存在一项为n-1，使之后所有项为1
- 存在j (j>=1 && j<=k)，满足a^{2^j-1q} mod n = n-1

Miller-Rabin算法实现

确定整数k和q ，满足n = 2^kq + 1
随机选择整数a ，满足a>1且 a<n-1
若a^q mod n = 1，返回不确定
若存在j (j>=1 && j<=k)，满足a^{2^j-1q} mod n = n-1，返回不确定
返回 "合数"

通过Miller-Rabin素性测试的数不一定是素数；无法通过MillerRabin素性测试的数一定不是素数（必要条件）

代码实现

std::string miller_rabin_prime_test(unsigned int n, unsigned int a) {
	// TODO: 对n进行非确定性素性测试

	//先特判下
	if (a <= 1 || a >= n - 1)	return "error";
	if (n == 2)	return "uncertain";			//2返回uncertain,虽然我们知道2是素数
	if (n & 1 == 0)	return "not_prime";		//除2外的偶数，不是素数
	if (n == 1)	return "not_prime";			//1不是素数

	//将n分解为2^k * q + 1
	int k = 1;
	int q = 1;
	size_t t = n - 1;
	while (((t / 2) & 1) == 0)	//除2后仍为偶数
	{
		k++;
		t /= 2;
	}
	q = t;
	
	//判定1
	if (mod_exp(a, q, n) == 1) return "uncertain";

	//判定2
	size_t temp = mod_exp(a, q, n);
	for (int j = 1; j <= k; j++)
	{
		if (temp % n == n - 1)	return "uncertain";
		temp = (temp * temp) % n;
	}

	return "not_prime";
}

乘法逆元

存在条件：a和m的最大公约数为1（a与m互素），记为gcd(a,m) = 1

若a与m互素，则存在整数k1，k2，满足k1*a + k2*m = 1
两边同取mod m得：k1*a = 1 mod m
即k1是a关于模m的乘法逆元

欧几里得算法

辗转相除，求最大公约数

比如gcd(a,b)：若b！=0

则gcd(b,a%b)

欧几里得扩展算法

在欧几里得算法的基础上，计算辗转相除的系数

实现代码：

int ex_gcd(int a, int m, int &x, int &y) {
	if (m == 0) {
		x = 1; y = 0;
		return a; // 到达递归边界返回上一层
	}
	int r = ex_gcd(m, a % m, x, y);
	int t = x; x = y; y = t - a / m * y;
	return r; // 得到最大公约数
}

解释：

r = x*m + y*(a%m)
∴ r = x*m + y*(a - a/m *m)
r = y*a + (x - a/m *y)*m
所以需要有t......

Gridea