【算法设计技巧】随机化算法

随机化算法(randomized algorithm)：在算法期间，随机数至少有一次用于决策。该算法的运行时间不只依赖于特定的输入，而且还依赖于所出现的随机数。在前面的数据结构中，随机化算法隐式地用在了完美散列和通用散列中。

随机数发生器

需要一种方法来生成随机数。实际上，真正的随机性在计算机上是不可能生成的，因为这些数将依赖于算法，从而不可能是随机的。因此，产生伪随机数就足够了。

产生随机数的最简单方法是线性同余数发生器(linear congruential generator)，它于1951年由 Lehmer 首先描述：数 x₁, x₂, … 的生成满足

    x_i+1 = Ax_i mod M

为了开始这个序列，必须给出 x₀ 的某个值，这个值叫作种子(seed)。如果 x₀ = 0，那么这个序列远不是随机的，但如果 A 和 M 选择得正确，则任何其他的 1 ≤ M 都是同等有效的。如果 M 是素数，那么 x_i 就绝不会是 0。例如，如果 M=11，A=7，而 x₀ = 1，则生成的数为 7, 5, 2, 3, 10, 4, 6, 9, 8, 1, 7, 5, … 在 M -1 = 10个数之后序列将重复。因此，这个序列的周期为 M-1。如果 M 是素数，那么总存在对 A 的一些选择能够给出全周期(full period) M - 1。如果 M 选择得很大，比如 31 比特的素数，那么对于大部分的应用来说周期应该是非常大的。Lehmer 建议使用 31 个比特的素数 M = 2³¹ - 1 = 2 147 483 647。

在许多库中的发生器基于函数

    x_i+1 = (Ax_i + C) mod 2^B

其中，B 的选取要匹配机器整数的比特位数，而 A 和 C 均为奇数。但是，这些发生器总是产生在奇偶之间交替的 x_i 的值——很难具有理想的性质。事实上，低 k 位是以周期 2^k循环的。

C++ 所提供的发生器包括线性同余发生器，带有类模板linear_congruential_engine，允许对参数 A、C和M进行类型说明，

template<typename UnsignedType, UnsignedType A, UnsignedType C, UnsignedType M>
class linear_congruential_engine;

同时还包括 typedef，

typedef linear_congruential_engine<unsigned int, 48271, 0, 2147483647> minstd_rand0;

此外，库中还提供了一个基于更新算法的发生器，叫作 Mersenne Twister，同时提供了 typedef mt19937。

使用C++11 随机数工具的类

#include<chrono>
#include<random>
#include<functional>
using namespace std;

class UniformRandom
{
public:
	UniformRandom(int seed=currentTimeSeconds()):generator{seed}{}

	//返回一个伪随机的 int
	int nextInt()
	{
		static uniform_int_distribution<unsigned int>distribution;
		return distribution(generator);
	}

	//返回一个在范围[0, high)之间的伪随机的 int
	int nextInt(int high)
	{
		return nextInt(0, high - 1);
	}

	//返回一个在范围[low...high]之间的伪随机 int
	int nextInt(int low, int high)
	{
		uniform_int_distribution<int>distribution(low, high);
		return distribution(generator);
	}

	//返回一个在范围[0...1) 之间的伪随机数 double
	double nextDouble()
	{
		static uniform_real_distribution<double> distribution(0, 1);
		return distribution(generator);
	}
	
private:
	mt19937 generator;
};

int currentTimeSeconds()
{
	auto now = chrono::high_resolution_clock::now().time_since_epoch();
	return (chrono::duration_cast<chrono::seconds>(now)).count();
}

跳跃表

素性测试

这一节将考查一个大数是否是素数的问题。

下面将给出一个可以测试素性的多项式时间算法。如果这个算法宣称一个数不是素数，则可以肯定这个数不是素数。若该算法宣称一个数是素数，则这个数将以高的概率，但不是 100%地肯定是素数。算法的关键是费马的著名定理。

费马小定理(Fermat’s Lesser Theorem)：如果 P 是素数，且 0＜A＜P，那么 A^P-1 ≡ 1 (mod P)。

例如，67 是素数，因此 2⁶⁶ ≡ 1 (mod 67)。这提出了测试一个数 N 是否是素数的算法：只要检验一下是否 2^N-1 ≡ 1 (mod N)。如果不成立则可以肯定 N 不是素数，反之 N 很可能是素数。例如，满足该等式但不是素数的最小的 N = 341。

如果 P 是素数且 0＜X＜P，那么 X² ≡ 1 (mod P) 仅有的两个解为 X = 1，P - 1 。

因此，如果在计算 A^N-1 (mod N) 的任一时刻发现违背了该定理，那么可以断言 N 肯定不是素数。

/**
* 实现基本的素性测试的函数
* 如果witness不返回1，那么n肯定是一个合数
* 通过计算 a^i (mod n)并随时查看 1 的非平凡的平方根来做到的
*/
HugeInt witness(const HugeInt& a, const HugeInt& i, const HugeInt& n)
{
	if (i == 0)
		return 1;

	HugeInt x = witness(a, i / 2, n);
	if (x == 0) //如果n递归地为合数，则停止
		return 0;

	//若我们找到1 的一个非平凡的平方根，则n不是素数
	HugeInt y = (x * x) % n;
	if (y == 1 && x != 1 && x != n - 1)
		return 0;

	if (i % 2 != 0)
		y = (a * y) % n;

	return y;
}

/**
* 在随机化素性测试中所查询witnes的次数
*/
const int TRIALS = 5;

/**
* 随机化素性测试
* 调整TRIALS以增加可信度
* n是要测试的数
* 如果返回false，则n肯定不是素数
* 如果返回true，则n很可能是素数
*/
bool isPrime(const HugeInt& n)
{
	Random r;

	for (int counter = 0; counter < TRIALS; ++counter)
		if (witness(r.randomHugeInt(2, n - 1), n - 1, n) != 1)
			return false;
	return true;
}

素性测试的随机化算法很重要，因为这些算法一直比那些最好的非随机化算法要明显地快。虽然其可能偶尔会产生错误的结果，但其发生的机会可以限制到足够小，可以忽略不计。