9章【同余】

蒙特卡罗算法【和一个模型联系起来】

舍伍德算法【】

拉斯维加斯算法【最不靠谱，思路：把该算法和固有算法相结合】

随机数

学习目标：

利用数据序列的随机性和概率 分布等点，设计解决问题的算法或提高 已有算法的效率。

随机性 (randomness) ：

某一事件集合中的各个事件所表现出来的不确定性。可能遵循某个概率分布。

 随机序列 (random sequence)/ 随机变量序列：

如果用 X 1， X 2 ……X n代表随机变量，这些随机变量如果按照顺序出现，就形成了随机序列。

(1) 序列中的每个变量都是随机的；

(2) 序列本身就是随机的。

对产生随机数的数学方法要求：

(1) 产生的数列具有随机样本的统计学性质

(2) 产生的数列要有足够的周期

(3) 产生数列要速度快，占用计算机内存少

线性同余法

(1) 同余：

设 a, b, m 为整数， m ＞ 0 ，若 a-b 为 m的整数倍，则称a 与b关于模m 同余；记为 a≡b(mod m) 。

(2) 同余式满足 :

反射性（ reflection ） , 即a ≡ a ( mod n )

对称性（ symmetry ） ,即由 a ≡ b ( mod n ), 可得 b ≡ a ( mod n )

传递性（ transitivity）,即由 a ≡ b ( mod n ), b ≡ c ( mod n ), 可得a ≡ c ( mod n )

(3) 剩余类 ：对于正整数 n ，把全体整数按 n 的余数分成 n 类，

每一类数的全体称为模 n 的一个剩余类（或称同余类）。

(4) 完全剩余系：

在模 n的每一个剩余类中各取一个元素，则这 n 个数就组成模 n的一个完全剩余系。

在完全剩余系中，每个成员对于模 n 的余数不同。

(5) 缩剩余系：

在任意一个模 n 的完全剩余系，仅保留与 n互素的那些数（共φ (n)个数，其中，φ(n)为欧拉函数），则这φ (n) 个数组成模 n的一个缩剩余系（或称既约剩余系、简化剩余系，简称缩系）。

(6) 欧拉函数： 对正整数 n ，欧拉函数是小于或等于 n的正整数中与 n 互质的数的数目 , 因此 φ(1)=1 。

(7)线性同余发生器(Linear Congruence Generator, LCG)

$\begin{Bmatrix} a_0=d & \\ an=(ba_{n-1}+c)mod m & n=1,2... \end{}$

d是初始值【种子】

统一术语：an-1是乘数，c是增量，m是模数

当d不变时，产生的数列永远不变；结果可预测

其中， b≥0, c ≥0, m ＞ 0, d≤m 。 m, c 互质， m为机器所能取得的大数

LCG（线性同余发生器）改善

【例9-1】用线性同余法生成一个随机序列。

计算模型：

选择系统时间作为种子，令b=1194211693L，c=12345L，m=2147483647(2^31-1)，则有

系统函数：

srand(time(NULL));

rand()

算法设计与描述
输入:种子值
输出：伪随机序列
//伪随机序列 const unsigned long b = 1194211693L; const unsigned long c = 12345L; static unsigned long rand_seed; void mysrand(unsigned long seed) { rand_seed = seed; } long myrand() { rand_seed = (rand_seed * b + c) % ((1 << 31) - 1);//左移相当于×2 return rand_seed; } double myrandf() { return (double)rand_seed / ((1 << 31) - 1); } int main() { int i; mysrand(time(0));//time(0)是种子数 srand(time(0));//系统函数 for (i = 0; i < 50; i++) { myrand(); cout<<"rand_seed:" << rand_seed;//或者是cout<<myrand(); cout <<"\tmyrandf():"<< myrandf(); cout<<"\trand():"<<(double)rand()/RAND_MAX<<endl; } return 0; }

代码：

 #include<iostream>
#include<math.h>
#include<ctime>
#include<cstdlib>
using namespace std;
//伪随机序列

const unsigned long b = 1194211693L;
const unsigned long c = 12345L;
static unsigned long rand_seed;

void mysrand(unsigned long seed)
{
	rand_seed = seed;
}
long myrand()
{
	rand_seed = (rand_seed * b + c) % ((1 << 31) - 1);//左移相当于×2
	return rand_seed;
}

double myrandf()
{
	return (double)rand_seed / ((1 << 31) - 1);
}

int main()
{
	int i;
	mysrand(time(0));//time(0)是种子数 
	srand(time(0));//系统函数 
	for (i = 0; i < 50; i++)
	{
		myrand();
		cout<<"rand_seed:" << rand_seed;//或者是cout<<myrand(); 
		cout <<"\tmyrandf():"<< myrandf();
		cout<<"\trand():"<<(double)rand()/RAND_MAX<<endl; 
	}
	return 0;
}

蒙特卡罗算法(Monte Carlo，MC)

它将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的解。

时间： 20世纪40年代美国在第二次世界大战

项目：研制原子弹的“曼哈顿计划”

作者：S.M.乌拉姆和J.冯·诺伊曼

命名：数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo

起源：1777年，法国数学家布丰（Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法求圆周率π

统计模拟—蒙特卡罗算（MC）

1.数值计算

【例9-2】利用随机函数计算圆周率π：

可能得不到正确结果 1.随机数符合概率（正态分布）2.点数足够（样本量足够多）

精度不高

(1) 数学模型：

将n根飞镖随机投向一正方形的靶子，计算落入此正方形的内切圆中的飞镖数目k。

假定飞镖击中方形靶子任一点的概率相等。

设圆的半径为r，圆的面积s_1=πr^2 , 正方形面积s_2=4r^2

由随机序列均匀分布的假设可知落入圆中的飞镖和正方形内的飞镖平均比为：k:n=πr^2:4r^2

由此可知：π=4k/n

取单位圆，取第一象限，则有长方形：0<=x<=1，0<=y<=1；圆形：x2+y2<=1

(2) 算法描述：

Step 1: 设定总点数

Step 2: 取随机数x, y

Step 3: 判断x2+y2<=1,是则圆点数加1

Step 4: 判断所取点数达到总点数，是，计算π值，否，重复执行step 2.

(3) 算法分析：O(n)

代码：

#include<iostream>
#include<ctime>
#include<cstdlib>
using namespace std;


int main()
{
	int n;//随机数个数=所有点数 
	int inside=0;//位于圆中的点数 
	double x,y;
	
	cout<<"输入随机数的个数：";
	cin>>n;
	
	srand((unsigned int)time(NULL)) ; //初始化种子数 
	for(int i=0;i!=n;i++)
	{
		x=rand()*1.0/RAND_MAX;
		y=rand()*1.0/RAND_MAX;
		if(x*x+y*y<=1)
			inside++;
	}
	printf("%6.5lf",4.0*inside/n);
	system("PAUSE>NUL");
	
	return 0;
}

【例9-3】

设 f(x) 是 [0,1] 上连续函数，且 0 ≤ f(x) ≤ 1。需要计算的积分为

，积分值等于图中的面积G

(1) 数学模型

单位正方形内均匀地作投点试验，则随机点落在曲线下面的概率为：

假设向单位正方形内随机地投入 n 个点 (xi, yi) 。

如果有 m 个点落入 G 内，则随机点落入 G 内的概率 m/n≈Pr

(2) 算法分析：

投 n 个点，判断其在 G 中的个数，时间复杂 度 O(n)

算法设计与描述
输入：投射的点数n
输出：面积G
void getG() { int n;//随机数个数=所有点数 double x,y,m=0; cout<<"输入随机数的个数："; cin>>n; srand((unsigned int)time(NULL)) ; //初始化种子数 for(int i=0;i!=n;i++) { x=rand()1.0/RAND_MAX; y=rand()1.0/RAND_MAX; if(y<=f(x)) m=m+1; } printf("f(x)面积=%6.5lf",m/n); system("PAUSE>NUL"); }

PS：

MC算法在一般情况下可以保证对问题的所有实例都以高概率给出正确解，但是通常无法判定一个具体的解是否正确。

【p正确(p-correct)】

设如果一个MC算法对于问题的任一实例得到正确解释的概率不小于p，p是一个实数，且 1/2≤p＜1，称MC算法是p正确的，且称p - 1/2是该算法的优势。

【一致的(consistent)】

如果对于同一实例，蒙特卡罗算法不会给出2个不同的正确解答。

【偏真】

设MC(x)是解某个判定问题D的蒙特卡罗算法，当MC(x)返回true时解总是正确的，仅当它返回false时有可能产生错误的解，称为偏真，反之为偏假。

偏真-说真就是真，说假可能是假。

偏假-说假就是假，说真可能是真。

【例9-4】素数测试：判断给定整数是否为素数。

梅森Mersenne 2^p－1(p为素数)

算法1：穷举

模型：

(1) n mod i ==0;

(2) 穷举范围为 2 ~ n^( 1/2) 。//根号n

时间复杂度：

O(n^( 1/2) );

若 n 为 m 位 10 进制数， 则为 O(10^( m/2) )

实现：

int prime( int n)

{

for(i=2; i<= n^(1/2) ; i++)

if(n mod i == 0)

return 0;

return 1;

}

算法2：取2~n^(1/2)随机数，判断其是否为n的因子

模型：

(1)d=random(2, n^( 1/2) );

(2) n mod d ==0 。

时间复杂度： O(1)

实现：

int prime( int n)

{

d=random(2, sqrt(n));

if(n mod d== 0) return 0;

return 1;

}

分析： 找到 n 的非平凡因子， n 为合数，算法完全正确，偏假

例： n=61*43, 根号n ≈51, prime 在 2~51内随机选一整数 d 。

成功： d=43 ，算法返回 false(概率为2%) ，结果正确。

失败： d≠43 ，算法返回 true( 概率98%),结果错误

算法3：

 问题分析

(1)Wilson定理：充要条件

对于给定的正整数 n, 判定 n 是一个素数的充要条件是 (n-1)!≡ -1(mod n) 。

//（n-1）!和n互质

//n的完全剩余系的相乘和n互质

证明：取集合 A={1, 2, 3, …… ， n-1}: 则 A 构成模 n 完全剩余系。

∀ i ∈ A，存在 j ∈ A ，使得： (i*j)≡1(mod n) ，则对于A中不能两两配对的成员x ，

考虑 x 2 ≡1(mod n) -> x^ 2 -1≡0(mod n) -> (x-1)(x+1) ≡0(mod n)

解得 x≡1(mod n) 或 x≡n-1(mod n)

其余两两配对：故而 (n-1)! ≡1*(n-1) ≡-1(mod n)

(2) Fermat小定理：——必要条件，非充分

341 符合但是合数 2^340 mod 341 =1 ,341=11*31

Fermat小定理是必要条件

满足 Fermat 小定理而不是素数的合数叫 Carmichael 数

10亿个自然数中素数50847534个， Carmichael数5597，

运用Fermat小定理出错的几率为0.00011

【引理 1 】若 a,b,c 为任意 3 个整数 ,m 为正整数，且 (m,c)=1 ，则当 a·c≡b·c(mod m) 时，

有 a≡b(mod m) 。

证明： a·c≡b·c(mod m) 可得 ac–bc≡0(mod m) 可得 (a-b)·c≡0(mod m) 。因为 (m,c)=1

即 m,c 互质， c 可以约去， a– b≡0(mod m) 可得 a≡b(mod m)

【引理 2 】设 m 是一个整数且 m>1 ， b 是一个整数且 (m,b)=1 。如果

a1,a2,a3,a4,…am 是模 m 的一个完全剩余系，则

b·a1,b·a2,b·a3,b·a4,…b·am 也构成模 m 的一个完全剩余系。

证明 : 若存在 2 个整数 b·ai 和 b·aj 同余即 b·ai≡b·aj(mod m)..(i>=1 && j>=1) ，

根据引理 1 则有 a[i]≡a[j](mod m) 。根据完全剩余系的定义可知这是不可能

的，因此不存在 2 个整数 b·a[i] 和 b·a[j] 同余。

问题分析

(1)Wilson 定理： 对于给定的正整数 n, 判定 n 是一个素数的充要条件是 (n-1)!≡ -1(mod n) 。【n-1的阶乘和n互质】

(2) Fermat 小定理： 如果 n 是一个素数 , 且 (0<a<n), 则 a^（ n-1） ≡1(mod n) //这个数可能很大

(3) 利用 Fermat 小定理构造快速判断素数方法

• 二次探测定理： 如果 n 是一个素数 , 且 0<a<n, 则方程 a*a ≡1(mod n) 的解为 a=1, n-1 。

a*a % n = 1

(a-1)*(a+1) % n = 0

因此 a=1或n-1

它是F小定理的一个特例

• 二次探测定理的应用：若 a≠1 且 a≠n-1 ，则必为合数。

偏假算法【说假就是假，说真不一定是真】

(3) 利用 Fermat 小定理构造快速判断素数方法

计算模型：

计算模型

a=Random(2, n-1) //从2~n-1中随机取一个数

令m=n-1

d=a^2 mod n if m mod 2==0

d=(a*(a^2 mod n)) mod n if m mod 2==1

若d==1且 a<>1 and a<>n-1，则n是合数

（4）算法设计

算法实现：

int power(unsigned int a,unsigned int b,unsigned int n)
{
	unsigned int y=1,m=b,z=a;
	while(m>0)
	{
		while(m%2==0)
		{
			//此时是偶数 
			int x=z;
			z=z*z%n;
			m=m/2;
			if( (z==1)&&(x!=1)&&(x!=(n-1)))//判断是否为合数 
				return 1;
		}
		m--;
		//此时是奇数 
		y=y*z%n;
	}
	if(y==1)
		return 0;
	return 1;
}
int prime(unsigned int n)
{
	int i,a;
	srand(time(NULL));
	for(i=1;i<log(n)/log(10) ;i++)
	{
		a=rand()%n;
		if(power(a,n-1,n))
			return 0;
	}
}

算法分析：

代码：

#include<iostream>
#include<ctime>
#include<cstdlib>
#include<math.h> 
//#include<random>
using namespace std;
//算法1：穷举 
int primeQiongJu(int n )
{
	//从2到根号n 
	for(int i=2;i<=sqrt(n);i++)
	{
		if(n%i==0)return 0;//合数 
	}
	return 1;//素数 
}

//算法2：随机数 ——多试几次，会提高正确率 
int random01(int x,int y)
{
	srand(time(NULL));
	int z=rand()%(y-x)+x;
	
	
	return z;
}
int primeSuiJi(int n)
{
//	int d=random(2,sqrt(n));//编译器不支持

	int d=random01(2,sqrt(n));//从2-根n中随机取一个数 
	if(n%d==0)
		return 0;//非素数 
	return 1;
} 

//03F费尔玛小定理
int power(unsigned int a,unsigned int b,unsigned int n)
{
	//a是随机数 n-1 n 
	//初始化
	unsigned int y=1,m=b,z=a,x ;
	//m是指数 
	while(m>0)
	{
		while(m%2==0)//偶数 
		{
			x=z;
			z=z*z%n;
			m=m/2;
			//判断 x是否为合数
			if( (z==1)&&(x!=1)&&(x!=(n-1)) ) 
				return 1;//合数 
		}
		m--;
		y=y*z%n;// 
		 
	}
	//退出时 m<=0
	if(y==1)
		return 0;//素数 
	else return 1;//合数 
}

int primeF(int n)
{
	int a; 
	//种子数
	srand(time(NULL));
	//为了提高正确率，根据n的位数
	//多运行log10(n)次算法 
	for(int i=1;i<log10(n);i++)
	{
		a=rand()%n;// 0-(n-1)
		if(power(a,n-1,n))
			return 0;//合数 
	} 
	return 1;//素数 
} 


int main()
{
	int n=53;
	if(primeQiongJu(n))
	{
		cout<<"01.穷举："<<n<<"is prime"<<endl;
	}
	if(primeSuiJi(n))
	{
		cout<<"02.随机数："<<n<<"is prime"<<endl;
	}
	if(primeF(n))
	{
		cout<<"03.随机数："<<n<<"is prime"<<endl;
	}
	
	return 0;
}

 舍伍德算法

Sherwood 基本思想

设 A 是一个确定性算法，当它的输入实例为 x时所需的计算时间记为 T A (x) 。

设 X n 是算法 A 所有输入规模为 n的实例的全体，

则当问题的输入规模为 n 时，算法A所需的平均时间为

【例9-5】设A是含有n个元素的集合，从A中选出第k小的元素，其中1≤k≤n。

快速排序分析

最坏情况：

T(1)=θ(1), T(n)=T(n-1)+T(1)+θ(n) -> O(n^ 2 ) ；

有 序且每次都取第一个元素，每次都分为 1 元素和 n-1 个元素

最好情况：

T(1)=θ(1) , T(n)=2T(n/2)+θ(n) -> O(nlogn)，每次都恰好二分

平均情况： O(nlogn) ，

查找第 k 小元素的时间复杂度至少： O(n)

计算模型:

(1) 数列存储于数组 a 中，下标从 0 开始；

(2)pivot=a[left + rand() % (right - left)], j=right, 一趟快排后，

则有以下三个结论。

 j-left=k-1 ，则分界数据就是选择问题的答案。

 j-left>k-1 ，则选择问题的答案继续在左子集中找，问

题规模变小了。

 j-left<k-1 ，则选择问题的答案继续在右子集中找，问

题变为选择第 k-left-1 小的数，问题的规模也变小。

算法设计与描述
输入：数列a[n] 和k 值
输出：第k小元素
#include<iostream> using namespace std; int select(int a[],int left,int right,int k) { int i,j,pivot,t; if(left>=right)//结束递归 return a[left]; i=left; j=left+rand()%(right-left); swap(a[i],a[j]); pivot=a[left]; j=right+1; while(true) { while(a[++i]<pivot); while(a[--j]>pivot); if(i>=j)break; t=a[i];a[i]=a[j];a[j]=t; } if(j-left+1==k)return pivot; a[left]=a[j]; a[j]=pivot; if(k<(j-left+1))select(a,left,j-1,k); else select(a,j+1,right,k-j-1+left); } int main() { }

代码：

#include<iostream>
#include<ctime> 
#include<cstdlib>
using namespace std;

int select(int a[],int left,int right,int k)
{
	int i,j,pivot,t;
	
	if(left>=right)//结束递归 
		return a[left];
	
	i=left;
	j=left+rand()%(right-left);
	swap(a[i],a[j]);
	pivot=a[left];
	j=right+1;
	while(true)
	{
		while(a[++i]<pivot);
		while(a[--j]>pivot);
		if(i>=j)break;
		t=a[i];a[i]=a[j];a[j]=t;
	}
	if(j-left+1==k)return pivot;
	a[left]=a[j];
	a[j]=pivot;
	if(k<(j-left+1))select(a,left,j-1,k);
	else select(a,j+1,right,k-j-1+left);
	
}

int main()
{
	int a[]={1,2,3,4,5,6,7,8,9};
	int k=3;
	cout<<select(a,0,8,3);
	 
}

拉斯维加斯算法

随机生成答案并检测 — 拉斯维加斯算法

算法思想[应用范畴]： 用随机序列代替有规律的枚举，判断枚举随机结果是否正确。

特征： 得到的解都是正确的，但是有时找不到解。求得正确解的概率也依赖于算法所用的时间。

原理： 通常采用 bool 型方法来表示拉斯维加斯算法。当算法

找到一个解时返回 true ，否则 false 。当返回 false 时，说明未得

到解，那么可再次独立调用该算法，在时间允许的情况一直

运算到出解为止。

【例9-6】n皇后问题

问题分析：

对于 n 皇后问题的任何一个解而言，每一个皇后在棋盘上的位置无任可规律，不具有系统性，像是随机放置的。适合 Las Vegas 。

模型：

(1) 数据结构： try[1]~try[8] 存放摆放列坐标， (i, try[i]) 表示存放位置；//i是行，try[i]是列

(2) 为第 i 行选取位置： try[i]=Random(1,8);

(3) 检测：当 try[i]≠try[k] and abs(try[i]- try[k])≠abs(i-k) 时，位置正确

算法设计：

(1) 选取位置 try[i]=random(1,8)

(2) 检测位置是否正确 check(try, i) ，正确转步骤 (3) ，不正确且未足 1000 次，转 (1);

(3) i=8 输出正确结果， i ≠8 且大于 1000 次，无解。

时间复杂度： O(n^ 2 )

算法实现：

代码：正确率太低了

#include<iostream>
#include<cstdlib>
#include<ctime>
using namespace std;
//n皇后问题
int check(int a[],int n)
{
	int i;
	for(i=1;i<n;i++)
	{
		//位于对角线 或者 位于同一列 
		if(abs(a[i]-a[n])==abs(i-n) || (a[i]==a[n]))
			return 0;//该种选择不成立 
	}
	return 1; 
} 
//
int lv(int ty[])
{
	int i;
	srand(time(NULL));
	for(i=1;i<=8;i++)
	{
		ty[i]=rand()%8+1;//从1到8
		if(!check(ty,i)) 
			return 0;//不成立 
	}
	return 1;
}

void queens()
{
	int ty[9]={0};//存储列坐标
	int success=0;
	int n=0;
	int i;
	while(!success && n<1000000 )
	{
		success=lv(ty);
		if(success)
		{
			cout<<"success:"<<n<<endl;
			for(i=1;i<=8;i++)
				cout<<ty[i]<<"\t";
			cout<<endl;
			success=0;
		}
		n++;
		for(i=1;i<=8;i++)
		{
			ty[i]=0;
		}
//		cout<<endl<<"n="<<n;
	} 
	
}
int main()
{
	queens();
	return 0;
}

算法改进：先在棋盘中随机放置若干行，后继算法用回溯，随机放置皇后越多，回溯时间越少。