二、概率算法
概率算法分类：
- 数字算法：求近似解
- Monte Carlo算法：解判定问题或只有一个正确解的问题。答案未必正确。
- Las Vegas算法：决不返回错误答案。
- Sherwood算法：总是给出正确答案。
2.1 数字算法
找到数字问题的近似解
\(\Pi\)值计算：投飞镖模拟，单位圆和外切正方形。
\[\frac{落在圆内}{总数} = \frac{k}{n} =\frac{\pi r^2}{4r^2} = \frac{\pi}{4} \\ \Rightarrow \pi = 4 \frac{k}{n} \]Darts(n){ k=0; for i=1 to n do { x = uniform(0,1); y = uniform(0,1); if(x*x+y*y<=1) then k++; } return 4k/n; }
数字积分：Monte Carlo积分，但不是本课中的MC算法。
计算积分\(S=\int_0^!f(x)dx\)
HitorMiss算法：单位面积正方形投飞镖，落入积分函数下方的为k，\(\frac{k}{n}=\frac{S}{1}\)
HitorMiss(f, n){ k=0; for i=1 to n do { x = uniform(0,1); y = uniform(0,1); if(y<=f(x)) then k++; } return k/n; }Crude算法：积分区间上随机取点，取算术平均值乘区间宽度。\(\int_a^bf(x)dx=(b-1)\frac{1}{n}\sum_{i=1}^{n}f(x_i)\)
Crude(f,n,a,b){ sum = 0; for i=1 to n do{ x = uniform(a,b); sum = sum + f(x); } return (b-a)*sum/n; }相对而言，Crude算法方差更小，但是HitorMiss算法在相同时间迭代次数更多。
解题思路：\(\pi=4\int_0^1\sqrt{1-x^2}dx\)
Trapezoid确定性算法：梯形算法，将区间划分为n-1个子区间
\(S=\delta*(f(a+\delta)+f(a+2\delta)+\cdots+\frac{f(a)+f(b)}{2})\)
Trapezoid(f,n,a,b){ // 假设n>=2 delta = (b-a)/(n-1); sum = (f(a)+f(b))/2; for(x=a+delta; x<=b-delta; x+=delta){ sum = sum + f(x); } return sum * delta; }确定性算法有时求不出解，概率算法适用于高维定积分。
概率计数：估计整数值，如集合的大小或者集合中不同对象数目的估计。
求集合的大小：有放回的随机抽取元素，k是第1次出现重复元素之前所选出的数目。
\(n足够大时,k=\beta\sqrt{n}=\sqrt{\frac{\pi}{2}n},即n=\frac{2k^2}{\pi}\)
SetCount(X){ k = 0; S = {}; a = uniform(X); do{ k++; S = S U {a}; a = uniform(X); }while(a不属于S) return 2k**2/pi }多重集合不同对象数目估计：将单词映射到长度为m的位串，\(\pi(y,b)=\begin{cases}i,y[i]=b的最小i \\ k+1,y[i]!=b \end{cases}\)
WordCount(){ y[1...m+1] = 0; for 磁带上的每个单词 x do { i = pi(h(x), 1); // x的散列值中等于1的最小位置，表示x是以 00...01....开头的 y[i] = 1; } return pi(y, 0); // 返回y中等于0的最小位置 } 2^(k-2)<=n<=2^k，更准确的n=2^k/1.54703
2.2 Sherwood算法
总是给出正确的答案，平滑不同输入实例的执行时间。
用途：平均性能较优，但是最坏性能较差时才使用。比如选择和排序，对元素的相对顺序敏感。
基本过程：
- 将被解实例变换到一个随机实例。（预处理）
- 用确定算法解随机实例，得到一个解。
- 将此解变换为原实例的解。（后处理）
RH(x){ // 用Sherwood算法计算f(x) n = length(x); // x的size为n r = uniform(A[n]); // 随机选取一个元素 y = u(x,r); // 将原实例x转化为随机实例y s = f(y); // 用确定算法求y的解s return v(r,s); // 将s的解还原为x的解 }
选择和排序的Sherwood算法：将输入实例打乱，无需后处理。在调用确定性算法前，先洗牌
Shuffle(T){
n = length(T);
for i=1 to n-1 do {
// 在T[1...n]中随机选取一个元素放在T[i]上
j = uniform(1...n);
sqap(T[i],T[j]);
}
}
离散对数计算：\(a=g^x mod \ p, 计算x=log_{g,p}a\)
简单算法：
- 计算\(g^x mod \ p,0 \le x\le p-1\)，x有限个取值，因为离散对数\(<g>\)是个循环群
- 验证\(a=g^x mod \ p\)是否成立
dlog(g,a,p){ // 当这样的算法不存在时，算法返回p x = 0; y = 1; do{ x++; y = y*g; // 计算y=g**x }while(a!=y%p && x!=p); return x; }定理：
- \(log_{g,p}(st \% p)=(log_{g,p}s+log_{g,p}t)\%(p-1)\)
- \(log_{g,p}(g^r \% p)=r,0\le r \le p-2\)
Sherwood改造算法：
dlogRH(g, a, p){ r = uniform(0..p-2); b = g^r mod p; // 求幂模b=g**r%p c = ba mod p; // 定理1 y = dlog(g,c,p); // 使用确定性算法求log_{g,p}c, y=r+x return (y-r) mod (p-1); // 求x }\[\begin{align} y &= log_{g,p}(ba \% p) = (log_{g,p}b+log_{g,p}a)\%(p-1) \\ &= (log_{g,p}(g^r \% p)+log_{g,p}a)\%(p-1) \\ &= (r+x)\%(p-1) \end{align} \]
搜索有序表：数值数组和指针数组构成有序的静态链表。
时间为\(O(n)\)的确定算法：直接顺序搜索，平均时间\(O(\frac{n}{2})\)
Search(x, i){ while x > val[i] do i = ptr[i]; return i; } A(x){ return Search(x, head); }时间为\(O(n)\)的概率算法：从随机的一个点左右开始顺序搜索，平均时间\(O(\frac{n}{3})\)
D(x){ i = uniform(1..n); y = val[i]; case{ x<y: return Search(x, head); // 从头开始搜索 x>y: return Search(x, ptr[i]); // 从下一个位置开始搜索 otherwise: return i; // case3 x==y } }平均时间为\(O(\sqrt{n})\)的确定算法：在前\(\left \lfloor \sqrt{n} \right \rfloor\)个数中找不大于x的最大整数的下标，从该数开始向后查找。
B(x){ // 设x在val[1..n]中 i = head; max = val[i]; // max初值是表val中的最小值 for j=1 to [sqrt(n)] do { // 在前根号n个数中找不大于x的最大整数y的下标 y = val[j]; if max < y <= x then { i = y; max = y; } } return Search(x, i); }
2.3 Las Vegas算法
特点：
- 更有效率的算法
- 时间上界可能不存在
- 可能找不到解，但如果找到一定是正确的
- 成功的概率随执行时间的增加而增加
算法的一般形式：\(LV(x,y,success),x为输入,y为输出\)
约定记号：
- \(p(x),算法成功的概率\)
- \(s(x),算法成功时的期望时间\)
- \(e(x),算法失败时的期望时间\)
一般会多次运行LV算法，知道成功
Obstinate(x){ repeat LV(x,y,success); until success; return y; }找到正确解的期望时间：
\[\begin{align} t(x) &= p(x)s(x)+(1-p(x))(e(x)+t(x)) \\ &=s(x)+\frac{1-p(x)}{p(x)}e(x) \end{align} \]
2.3.1 八皇后问题
要求：
- 行不冲突：行号相等
- 列不冲突：列号相等
- 主对角线不冲突：行列号之差相等
- 辅对角线不冲突：行列号之和相等
传统回溯法：分别决定每一行放在哪个位置
i = j = 1;
while i ≤ 8 do { // 当前行号i≤ 8
检查当前行i：从当前列j起向后逐列试探，寻找安全列号；
if 找到安全列号 then {
放置皇后，将列号记入栈，并将下一行置为当前行(i++);
j=1; //当前列置为1
} else {
退栈回溯到上一行，即i--;
移去该行已放置的皇后，以该皇后所在列的下一列作为当前列;
}
}
Las Vegas算法：
约定：
- \(try[1...8]表示第i个皇后放在(i,try[i])上\)
- \(try[1...k]称为k-promising\)，当且仅当满足k皇后要求
- \(\forall i!=j \in [1,k],try[i]-try[j]\notin \{i-j,0,j-i\}\)
基本步骤：
- 遍历改行所有的列位置是否可行，\(\frac{1}{nb}\)概率下选用该列
- 如果有解则放，没解则退出
QueensLV(success){
// 贪心的LV算法，所有皇后都是随机放置
col, diag45, diag135 = {}; // 列、两个对角线初始化为空集
k = 0; // 行号
repeat
nb = 0; // 计数器，nb值为(k+1)th皇后的open位置总数
// 遍历列号，试探(k+1,i)是否安全
for i=1 to 8 do {
if (i!∈col) and (i-k-1!∈diag45) and (i+k+1!∈diag135) then {
// 列i对(i+1)th皇后可用，但不马上放置
nb = nb + 1;
if uniform(1..nb)=1 then // 或许放在第i列
j = i;
}
}
if(nb>0) then {
// nb=0时无安全位置
// 在所有nb个安全位置上，(k+1)th皇后选择位置j列的概率为1/nb
k = k+1;
try[k] = j;
col = col ∪ {j};
diag45 = diag45 ∪ {j-k};
diag135 = diag135 ∪ {j+k};
}
until nb=0 or k=8; // 当找不到合适的位置，或者try是8-promising时，退出
success = (nb>0);
}
算法12行，在所有可能的列位置上，选中某一列的概率都是\(\frac{1}{nb}\)，从最后一项往前倒推。
改进算法：联合回溯法和LV算法，指定使用LV算法只确定前几行的位置。
...上述算法...
until nb=0 or k=stepVegas; // 已随机放好stepVegas个皇后
if nb>0 then
backtrace(k, col, diag45, diag135, success);
else
success = false;
2.3.2 模p平方根
\(x \equiv y^2 (mod \ p),x\in[1,p-1],p为奇素数\)，则
- x为模p的二次剩余
- y为x模p的平方根
定理：
- 任何一个模p的二次剩余有两个不同的平方根 \(b,(p-b)\)
- 1到p-1之间的整数恰有一半时模p的二次剩余
- \(\forall x \in [1,p-1],p为奇素数,有x^{(p-1)/2} \equiv \pm1(mod \ p)\)
- \(x是模p的二次剩余，当且仅当x^{(p-1)/2} \equiv 1(mod \ p)\)
如何判定x是否为模p的二次剩余？计算\(x^{\frac{p-1}{2}}(mod \ p)\)是否为1。
如何计算x模p的两个平方根？
-
\(if \ p\equiv 3 (mod \ 4)\)，则平方根为 \(\pm x^{(p+1)/4}\)
-
\(if \ p\equiv 1 (mod \ 4)\)，只能使用LV算法
例子：\(p=53\equiv1(mod \ 4),x=7\)，求7模53的平方根。
基本步骤：
- 随机取\(a\in[1,p-1]\)
- 计算c和d，\((a+\sqrt{x})^{(p-1)/2} \equiv c+d\sqrt{x} \ (mod \ p)\)
- 若\(d=0\)，则返回假
- 若\(c=0\)，则找到y，\(dy \equiv 1 \ (mod \ p)\)
rootLV(x, p, y, success){ // 计算y
a = uniform(1..p-1); // 并不知道a取多少
if a**2 % p == x % p then { // 可能性很小
success = true; y =a;
}else{
计算c和d，使得0<=c,d<=p-1，(a+sqrt(x))**((p-1)/2) % p == c+d*sqrt(x) % p;
if d==0 then success = false;
else{ // c=0
success = true;
计算y使得，d*y % p == 1, 1<=y<=p-1 // 修改Euclid算法可求
}
}
}
算法成功的概率大于0.5，因此平均调用两次即可求得x的平方根。
2.3.3 整数的因数分解
整数的素因数分解，\(n=p_1^{m_1}p_2^{m_2}\cdots p_k^{m_k}\)
朴素split算法：找到n的一个非平凡因数
split(n){
// n是素数返回1，否则返回找到的n的最小非平凡因数
for i=2 to [sqrt(n)] do
if n % i == 0 then return i;
return 1;
}
合数的二次剩余：上节的二次剩余要求p为奇素数，则恰有两个不同的平方根。若\(n=pq,p和q为不同素数\)，则二次剩余恰有4个平方根。
Dixon因数分解算法：
基本步骤：
- 找到与n互素的整数a和b，使得 \(a^2 \equiv b^2 \ (mod \ n) 且 a !\equiv \pm b \ (mod \ n)\)
- 则n和a+b的最大公因子是一个非平凡因子，\(x=gcd(n, a+b)\)。
Dixon(n,x,success){ // 找合数n的某一非平凡因子
if n是偶数 then {
x = 2; success = true;
}else{
for i=2 to [log_{2}n] do {
if n**(1/i) 是整数 then {
x = n**(1/i); success = true; return;
}
}
找到a, b，使得 a**2 % n == b**2 % n;
if a % n == +-b % n then success = false;
else {
x = gcd(a+b, n); // 最大公因子
success = true;
}
}
}
k-平滑：一个整数x的所有素因子均在前k个素数中。
确定a和b的取值：提前选定n和k。
随机选择\(x\in[1,n-1]\)。
- 若x不与n互素，则找到了非平凡因子
- 否则，\(y=x^2 \ mod \ n\)，若y是k平滑的，则将x和y的因数分解保存在表里。
- 一直重复，直到选择了k+1个互不相同的整数x。
在k+1个因式分解式中找到一个非空子集\(\{y_1y_2\cdots y_k\}\)，使得相应的因式分解的积中前k个素数的指数均为偶数（包含0），如\(y_1y_2y_4y_8=2^63^45^47^011^213^217^2\)
令\(S=\{y_{s_1},y_{s_2},\cdots,y_{s_k}\},a=\prod_{y_i\in S} x_i,b=\sqrt{\prod_{y_i\in S} y_i}\)，\(if \ a !\equiv \pm b \ (mod \ n)\)，则只需要求a+b和n的最大公因子即可。
时间分析：k的取值将会影响 \(x^2 \ mod \ n\)是k平滑的可能性，以及因数分解y时的时间。通常取如下时间
2.4 Monte Carlo算法
存在某些问题无法找到有效的算法获得正确解。MC以高概率找到正确解。
约定：
- p-正确：一个MC算法以不小于p的概率返回一个正确的解
- 相容的、一致的：一个MC算法对同一实例返回相同的正确解。
多次调用同一个算法，选择出现次数最多的解，以增加一致的，p-正确算法成功的概率。
有偏算法：
- 偏真算法：MC算法返回true时的解一定正确，只有返回false时才有可能产生错误的解。
- 偏\(y_0\)算法：算法返回\(y_0\)时总是正确的，返回非\(y_0\)时p概率正确。
重复调用一个一致的，p-正确的，偏\(y_0\)的MC算法k次，得到一个 \(1-(1-p)^k\) -正确的算法。
2.4.1 主元素问题
一个具有n个元素的数组，若某一元素的个数大于\(\frac{n}{2}\)，则其为主元素。
判断T是否有主元素：随机设置一个x，然后看看有多少个
maj(T){
i = uniform(1..n);
x = T[i];
k = 0;
for j=1 to n do
if T[j] == x then
k = k+1;
return k>n/2;
}
偏真，0.5-正确的算法
重复调用降低错误率：
maj2(T){
if maj(T) then return true; // 1次成功
else return maj(T); // 调用2次
}
偏真，0.75-正确的算法
算法改进：控制算法出错概率小于\(\varepsilon > 0\)，应调用次数为\(k=\left \lceil lg\frac{1}{\varepsilon} \right \rceil\)
// e 为错误概率
majMC(T, e){
k = [lg(1/e)]; // 上取整
for i=1 to k do
if maj(T) then return true;
return false;
}
// O(nlg(1/e))
2.4.2 素数测定
判定一个整数是否为素数，尚未找到有效的确定性算法或LV算法
简单的概率算法：
prime(n){
d = uniform(2..[sqrt(n)]); // 下取整
return (n mod d)!=0;
}
Fermat小定理：\(if \ n 是素数,则\forall a \in [1,n-1],有a^{n-1} mod \ n = 1\)
变换命题：\(if \ \exist a \in [1,n-1],使a^{n-1}mod \ n != 1, 则n不是素数\)
素性测定算法：使用Fermat小定理，偏假的。即算法返回假，一定为假。
Fermat(n){
a = uniform(1..n-1);
if a**(n-1) mod n == 1 then return true; // 未必正确
else return false; // 一定正确
}
伪素数和素性伪证据：符合Fermat小定理的合数，称为以a为底的伪素数，a称为n的素性伪证据。
将Fermat测试重复多次，也不能降低误差。
强伪素数：n是大于4的奇数，s和t是整数，使得\(n-1=2^st,t为奇数\)，集合\(B(n)\)满足如下条件：
- n为素数，则 \(\forall a \in [2,n-2], a \in B(n)\)
- n为合数，\(if \ a \in B(n)\)，则称n为一个以a为底的强伪素数，a为n素性的强伪证据。
Btest(a, n){
// n为奇数，2<=a<=n-2, 返回true则代表a属于B(n)。则n为强伪素数或素数
s=0; t=n-1; // t开始为偶数
repeat
s++;t/=2;
until t mod 2 = 1; //n-1=(2**s)t, t为奇数
x = a**t mod n;
if x=1 or x=n-1 then return true; // 满足条件1或2
for i=1 to s-1 do{
x = x**2 mod n;
if x=n-1 then return true; // 满足条件2
}
return false;
}
0.75-正确
强伪证据的数目比伪证据数目少很多。
\(强伪证据 \Rightarrow 伪证据\)
Miller-Rabin测试：0.75-正确，偏假的
MillRab(n){ // 奇n>4，返回真时表示素数，假时表示合数
a = uniform(2..n-2);
return Btest(a,n); // 测试n是否为强伪素数
}
RepeatMillRob：重复实验k次，\(1-4^{-k}\)-正确的MC算法
RepeatMillRob(n,k){
for i=1 to k do
if MillRob(n) == false then return false; // 一定是合数
return true;
}
若给出出错概率不超过\(\varepsilon\)，则重复次数为 \(k=\left \lceil \frac{1}{2}lg\frac{1}{\varepsilon} \right \rceil\)
确定n素性时间为：\(O(lg^3n \ lg\frac{1}{\varepsilon})\)
2.4.3 矩阵乘法验证
验证\(矩阵AB=C\)
设 \(X_{1 \times n}=(a_1,a_2,\cdots,a_n),a_i\in\{0,1\}\)，则问题转化为 \(XAB=XC\)
基本步骤：
- 计算\(XA\)
- 计算\(XA\)与\(B\)的乘积
- 计算\(XC\)
GoodProduct(A,B,C,n){
for i=1 to n do
x[i] = uniform(0..1);
if XAB=XC then return true;
else return false;
}
偏假的，0.5-正确
返回假一定为假
改进算法：重复k次
RepeatGoodProduct(A,B,C,n,k){
for i=1 to k do // 重复k次
if GoodProduct(A,B,C,n) == false then return false;
return true;
}
偏假的，\(1-2^{-k}\)-正确
给出错误概率：即$2^{-k}=\varepsilon $
GP(A,B,C,n,e){
k = [lg(1/e)]; // 上取整
return RepeatGoodProduct(A,B,C,n,k);
}
// O(n**2log(1/e))