chap 1-5 数论算法-素数

1. 费马素性测试-1

上一章提到过两数互质的问题，质数（素数）一直以来都是数论中的重点内容。一个基本问题就是：如何判断一个数是否为素数？最容易想到的方法：将这个数做质因数分解，若因数只有1和它本身，则为素数，但是质因数分解效率很低；另一种方法是，依次判断是否有比它小的数整除它，本质上还是质因数分解。下面要介绍的一种高效的素性测试方法，基于一个重要的定理：

费马小定理(Fermat’s little theorem)：

若 $N$为素数，则对于任意 $1\leq a < N$，都有 $a^{N-1}\equiv 1(\mod N)$成立。

$Proof$：设集合$S$为非零整数模$N$(素数)的完全剩余系，即$S$ = {$b_1,b_2,\cdots,b_{N-1}$}，其中$b_i=i$。设$a$为一个正整数，则可以证明 {$ab_1\mod N, ab_2\mod N, \cdots, ab_{N-1}\mod N$} 仍然等价于$S$：反证法，假设存在$ab_i\equiv ab_j(\mod N),b_i\neq b_j$，由于$a$与$N$互质，可以推出$b_i\equiv b_j(\mod N)$，矛盾，故结论成立。

综上所述：$S$ = {$1,2,\cdots,N-1$} = {$a\cdot 1\mod N, a\cdot 2\mod N, \cdots, a\cdot (N-1)\mod N$}，把它们所有的元素相乘得到：$(N-1)!\equiv a^{N-1}\cdot (N-1)! (\mod N)$。又因为$(N-1)!$与$N$互质，可以推出$a^{N-1}\equiv 1(\mod N)$。▮

利用费马小定理，再加上有很快的快速幂取模算法，我们可以设计出下面这个似乎合理的素性测试方法（判断输入$N$是否为素数）：

bool primality(int N) {
	if (N == 1) return false;
	// pick a positive integer < N at random
	int a = rand() % (N - 1) + 1;
	if (modexp(a, N - 1, N) == 1) return true;
	else return false;
}

仔细分析一下，若$N$确实为素数，则对于所有的$1\leq a < N$都应该能通过测试，这没问题。但是，如果$N$为合数，是否一定意味着所有的$1\leq a < N$都不能通过测试？答案是否定的，$N$为合数时，仍然存在$1\leq a < N$能通过测试。这说明这里的算法有一定概率失败，即在$N$为合数的情况下，输出$N$为素数（伪素数）。那是否能把失败的概率降到最低？

2. 费马素性测试-2

既然当$N$为合数时，一部分$a$通过测试，一部分不通过，而且只要存在不通过的，说明$N$肯定是合数，那么我们可以多进行几次上述测试，使失败概率降低：

// k=100
bool primality2(int N) {
	if (N == 1) return false;
	for (int i = 0; i < 100; ++i) {
		int a = rand() % (N-1) + 1;
		if (modexp(a, N - 1, N) != 1) return false;
	}
	return true;
}

失败的概率具体是多少呢？

假设整数$a(1\leq a < N)$与$N$互质，且$a^{N-1}\not\equiv 1(\mod N)$，$a$将无法通过测试。假设有若干整数$b(1\leq b < N)$满足$b^{N-1}\equiv 1(\mod N)$即能通过测试，那么：

$$
(a\cdot b)^{N-1}\equiv a^{N-1}\cdot b^{N-1}\equiv a^{N-1}\not\equiv 1 (\mod N)
$$

说明，只要存在一个通过测试的$b$，就必存在一个无法通过测试的$a\cdot b$。那这些$a\cdot b \mod N$是否相等？固定$a$不变，选择不同的$b$，显然$ab_i\not\equiv ab_j (\mod N)$。由此我们得到结论：无法通过测试的元素个数$\geq$通过测试的元素个数，即最多只有一半的元素能通过测试。

从群论的角度也能证明这个结论。所有模$N$后还与$N$互质的正整数构成了乘法群$Z_N$，集合$B$ = {$b:b^{N-1}\equiv 1 (\mod N)$}是$Z_N$的子群，说明$B$中的元素个数一定整除$Z_N$的元素个数。若$B$不含有$Z_N$的所有元素，那么$B$最多含有$|Z_N|/2$个元素。

这个结论说明了以下事实：

$$
Pr(primality(N)=true,N=prime)=1
$$

$$
Pr(primality(N)=true,N\neq prime)\leq \frac{1}{2}
$$

所以改进后的$primality2(N)$的失败概率为：

$$
Pr(primality2(N)=true,N\neq prime)\leq \frac{1}{2^k}
$$

若多次测试的次数$k=100$，那么该算法失败的概率最多为$2^{-100}$，可以认为不可能失败。

3. Miller-Rabin 素性测试

上述关于概率的推导存在一个前提：假设整数$a(1\leq a < N)$与$N$互质，且$a^{N-1}\not\equiv 1(\mod N)$。但是有一类数，很少但是存在，叫做Carmichael数，它的特点是：对于所有与$N$互质的整数$a(1\leq a < N)$，$a^{N-1}\equiv 1(\mod N)$，比如561。对这些数进行费马素性测试，将无法保证上述$\frac{1}{2^k}$的错误率。

有一种叫Miller-Rabin素性测试的算法可以解决这个问题。假设需要被测试的整数为$N$，则$N-1$可以被表示为$2^tu$的形式。计算下面的序列：

$$
a^u\mod N,a^{2u}\mod N,\cdots,a^{2^tu}=a^{N-1}\mod N
$$

若$a^{N-1}\not\equiv 1(\mod N)$，根据费马小定理，$N$一定是合数；但是如果$a^{N-1}\equiv 1(\mod N)$，$N$有可能是伪素数，需要进行下面的检验：在上述序列中，如果存在某个元素（非第一个）的值等于1但是它前一个元素的值不等于$(-1)\mod N$（即$N-1$），那么$N$就是合数。

$Proof$：假设$x^2=1+kN$，则$(x+1)(x-1)=kN$，说明$x\pm1\equiv 0 (\mod N)$，即$x\equiv\pm 1(\mod N)$。这里的$x^2$就是上述等于1的元素，$x$为它前一个元素。▮

实际上，Miller-Rabin素性测试就是在费马素性测试中加入了额外的检测。代码如下：

// Miller-Rabin
bool primality3(int N) {
	if (N == 1) return false;
	for (int i = 0; i < 50; ++i) {
		int a = rand() % (N - 1) + 1;

		int t = 0;
		int u = N - 1;
		while (u % 2 == 0) {
			u /= 2;
			t++;
		}

		int prev = -1;
		int curr = modexp(a, u, N);
		int val = 0;
		for (int i = 0; i < t; ++i) {
			prev = curr;
			curr = modexp(curr, 2, N);
			if (curr == 1 && prev == -1) val = prev;
		}

		if (curr != 1) return false;
		else {
			if (val == N - 1) return false;
		}

	}
	return true;
}

查阅相关证明，可知对于所有的$N$（包含Carmichael number）：

$$
Pr(primality3(N)=true,N\neq prime)\leq \frac{1}{4^k}
$$

错误率比以前更低，时间复杂度为$O(kn^3)$。这也是实际应用中，比较常用的素性检测方法。

4. 生成随机素数

由于素数的诸多性质，很多算法会提前建立好素数表，那么该如何高效地生成$n-bit$素数呢？一个简单的想法就是：随机取一个数，进行素性测试，若通过测试则输出；若不通过，则重复以上过程。

int gen_random_prime(int lb, int ub) {
	int N = rand() % (ub - lb) + lb;
	while(!primality3(N)) N = rand() % (ub - lb) + lb;
	return N;
}

测试结果：

gen_random_prime_test
Input the lower bound: 1
Input the upper bound: 100
Input the num of primes: 10
Primes generated:
17
61
13
11
17
7
37
53
47
83
time = 6ms
======================================

比较棘手的问题是，我们不确定会重复多少次素性测试，那该如何分析这个算法的复杂度？借助下面的定理：

拉格朗日素数定理：

设 $\pi (x)$为$\leq x$的素数的个数，则$\pi (x)\approx x/\ln x $，或 $\lim\limits_{x\to+\infty}\frac{\pi(x)}{x\ln x}=1 $。

可以推出：$n-bit$的所有数中随机取一个，大约有$p=\frac{1}{n}$的概率是素数。令随机变量$X$表示该算法需要执行循环的次数，则：

$$
E[X]=\sum_{i=1}^{\infty}i\cdot P[X=i]=\sum_{i=1}^{\infty}i\cdot (1-p)^{i-1}p=p\cdot \frac{d}{dp}(\sum_{i=1}^{\infty}-(1-p)^i)=p\cdot \frac{1}{p^2}=\frac{1}{p}
$$

注意到$p=\frac{1}{n}$，则该算法平均需要$O(n)$次循环之后停止，故平均复杂度是$O(kn^4)$。这种情况下，算法执行的computer steps不确定，故只能求出平均的step。

5. 随机算法

本章所讲的算法都具有一定的随机性。令人惊讶的是，对于一些问题的求解，最快的算法都是随机化的。我们不需要算法一定输出正确的结果，我们只需要算法以高概率输出正确结果。

随机算法一般分为两类：

1. Las Vegas算法：永远输出正确的结果。如生成随机素数。
2. Monte Carlo算法：输出有出错的概率。如素性测试。

随机算法在很多方面都有应用（后面会讲到）：Hashing、排序、中位数、min-cut问题等等。

Quick Link: 算法笔记整理