（13）幂模m与逐次平方法

逐次平方法

如何计算

$5^{100000000000000}(\mathrm{m}\mathrm{o}\mathrm{d}12830603)$

呢？如果12830603是素数，你会设法使用费马小定理，即使不是素数，也可以利用欧拉公式。事实上，$12830603=3571\cdot 3593$且

$\varphi (12830603)=\varphi (3571)\cdot \varphi (3593)=3570\cdot 3592=12823440$

欧拉公式告诉我们，对任何a与m，若$gcd(a,m)=1$，则

$a^{\varphi (m)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$

所以可利用事实

$100000000000000=7798219\cdot 12823440+6546640$

来“简化”我们的问题：

$5^{100000000000000}=(5^{12823440}{)}^{7798219}\cdot 5^{6546640}$
$\equiv 5^{6546640}(\mathrm{m}\mathrm{o}\mathrm{d}12830603)$

现在“只”需计算5的$6546640$次幂，然后用模$12930603$进行简化。不幸的是，数$5^{6546640}$有400多万位数，即使用计算机计算也是很困难的。

用来计算$a^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$的一个巧妙想法叫做逐次平分法。

逐次平分计算$a^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$

1.将k表成2的幂次和：

$k=u_0+u_1\cdot 2+u_3\cdot 2^2+\cdots +u_r\cdot 2^r$

其中每个$u_i$是0或1。（这种表示式叫做k的二进制展开。）

2.使用逐次平分法制作模m的a的幂次表。

$a^1\equiv A_0(\mathrm{m}\mathrm{o}\mathrm{d}m)$
$a^2\equiv (a^1{)}^2\equiv A_0^2\equiv A_1(\mathrm{m}\mathrm{o}\mathrm{d}m)$
$a^4\equiv (a^2{)}^2\equiv A_1^2\equiv A_2(\mathrm{m}\mathrm{o}\mathrm{d}m)$
$a^8\equiv (a^4{)}^2\equiv A_2^2\equiv A_3(\mathrm{m}\mathrm{o}\mathrm{d}m)$
$\cdots$
$a^{2r}\equiv (a^{2r-1}{)}^2\equiv A_{r-1}^2\equiv A_r(\mathrm{m}\mathrm{o}\mathrm{d}m)$

注意要计算表的每一行，仅需要取前一行最末的数，平方它然后用模m简化。也注意到表有$r+1$行，其中r是第1步中k的二进制展开式中2的最高指数。

3.乘积

$A_0^{u_0}\cdot A_1^{u_1}\cdot A_2^{u_2}\cdots A_r^{u_r}(\mathrm{m}\mathrm{o}\mathrm{d}m)$

同余于$a^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$。注意到所有$u_i$是0或1，因此这个数实际上是$u_i$等于1的那些$A_i$的乘积。
证明

$a^k=a^{u_0+u_1\ast 2+u_2\ast 2^2+\cdots +u_r\ast 2^r}$
$=a^{u_0}\cdot (a^2{)}^{u_1}\cdot (a^{2^2}{)}^{u_2}\cdots (a^{2^r}{)}^{u_r}$
$=A_0^{u_0}\cdot A_1^{u_1}\cdot A_2^{u_2}\cdots A_r^{u_r}(\mathrm{m}\mathrm{o}\mathrm{d}m)$

代码实现

int power(int e,int n) {//无模运算
	if (e == 1 || n == 0)return 1;
	return n % 2 == 1 ? e * power(e*e, n / 2) : power(e*e, n / 2);
}
int power(int e, int n, int m) {//加入了取模运算
	if (e == 1 || n == 0)return 1;
	return n % 2 == 1 ? (e * power((e*e)%m, n / 2, m)) % m : power((e*e)%m, n / 2, m);
}
int power(int e, int n, int m) {//位运算优化
	if (e == 1 || n == 0)return 1;
	return n & 1 ? (e * power((e*e) % m, n >> 1, m)) % m : power((e*e) % m, n >> 1, m);
}

快速证明合数

之前我们要证明一个数是合数，需要试着用不超过$\sqrt{n}$的每个数去除，查看能否找到因数。这样显然效率很低。
有了费马小定理与逐次平方法，我们就有办法快速的证明一个数是合数，而不用求出任何它的因数。
假设我们想证明m为合数。
（1）我们取一个小于m的数a。
（2）计算$gcd(m,a)$，如果$gcd(m,a)\ne 1$那么说明我们已经找到了m的一个因数，得证。
（3）如果$gcd(m,a)=1$，那么计算

$a^{m-1}(\mathrm{m}\mathrm{o}\mathrm{d}m)$

如果m是素数那么一定有$a^{m-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$，所以如果$a^{m-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$，那么m一定不是素数即是个合数，这样我们就不用求任何因数证明了一个数是合数。
其中$a^{m-1}(\mathrm{m}\mathrm{o}\mathrm{d}m)$可以通过逐次平方法快速求解。

我们可以通过改变a的值，如$a=5,7,11\cdots$发现m使$a^{m-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$是不是可以断言m是一个素数呢？很有可能，可惜不一定。

注意通过费马小定理不能结论性的证明一个素数，因为的确存在合数使得任意a，$a^{m-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$这样的m被称为卡米歇尔数。最小的卡米歇尔数是561，之后我们会进一步研究卡米歇尔数。

代码实现

#include<iostream>
using namespace std;

int n;

int Eucild(int a, int b)
{
	int r = a % b;
	while (r)
	{
		a = b;
		b = r;
		r = a % b;
	}
	return b;
}

int power(int e, int n, int m) {
	cout << e << ' ' << n << ' ' << m << endl;
	if (e == 1 || n == 0)return 1;
	return n & 1 ? (e * power((e*e) % m, n >> 1, m)) % m : power((e*e) % m, n >> 1, m);
}

bool isCombinde(int n) {
	int a = 2;
	if (n < a)return 0;
	int m = Eucild(n, a);
	if (m != 1)return 1;
	else if (power(a, n - 1, n) != 1)return 1;
	else return 0;
}

int main() {
	
	cin >> n;

	if (isCombinde(n))cout << n << "一定是合数" << endl;
	else cout << n << "可能是素数" << endl;

	system("pause");
	return 0;
}