算法设计与分析期末大作业

1.算法作业一

1.1.1.问题描述：

通常，股民希望通过低买高卖来获取股票的最大收益。比如一支股票，5天内价格分别为[10，11，7，10，6]。如果第1天以10元的价格买进，第2天以11元的价格卖出，则收益为1元；如果第3天以7元的价格买进，第4天以10元的价格卖出，则收益为3元。显然，3元是最大收益。
由此，提出问题：对于任意$n$天的股票价格序列，希望找出最佳的买入点和卖出点。

1.1.2.算法设计：

对于任意$n$天的股票价格，有人提出一种算法I：先找到这组股票价格中的最低值和最高值，然后从最低值出发向右股价序列中找到一个最大值；此外，从最高值出发向左股价序列中找到一个最小值。比较这两种情况下得到的收益值，选择收益大的作为买进和卖出的时间点。
对于任意$n$天的股票价格，还可以有算法II（穷举法）：对于股价序列，两两求出它们之间的差值，并记录当前的最优收益值。
对于任意$n$天的股票价格，还可以有两种分治求解策略：
（1）算法Ⅲ（顺序对分治）：先用分治算法求出该组股价序列中的（从小到大）顺序对（顺序对的概念参照PPT课件的逆序对），再来求解差值最大的顺序对，以此作为最优投资的选择。
（2）算法IV（直接分治）：假设该序列存在函数$max\_profit\_dc(price)$可求得最佳买卖点，则利用分治递归策略，不妨将序列一分为二，每个子序列可利用相同的$max\_profit\_dc(price)$来求解最佳买卖点，依此类推，直至子问题继续分解为最小规模的子问题，再进行合并回代。合并过程中，需要注意最佳买卖点不一定都在子问题内部找到，也可能在子问题的分界点上实现。

1.1.3.作业要求：

（1）上述算法I，能够保证任何情况下都能获得最佳买卖点吗？如果不能，试举一个反例。

（2）试分析穷举法和两种分治算法的算法复杂度。

（3）利用算法Ⅲ（顺序对分治）或算法Ⅳ（直接分治），编写一个程序，计算规模为$n(≥10)$的一个股价序列的最佳买卖点。数据自定义，给出结果，并附代码。

1.2.1.算法Ⅰ的分析

算法Ⅰ并不能保证任何情况下都能获得最佳买卖点的。要证明该算法的错误，只需要举出一个反例即可。

譬如题干所述实例$[10, 11, 7, 10, 6]$，有$[10, 11_{max}, 7, 10, 6_{min}]$，通过算法Ⅰ得到的最佳买卖点为$\{10,11\}$，但实际最佳买卖点为$\{7,10\}$。即算法Ⅰ的策略是错误。

1.2.2.算法Ⅱ的分析

算法Ⅱ（穷举法）通过两两对比求出差值。核心代码如下：

1
2
3
4
5
6
7
8
9
10
11

int Max=-0x3f3f3f,start=0,over=0;
for(int i=1;i<=n;i++){
    for(int j=i+1;j<=n;j++){
        if(day[j]-day[i]>Max){
            Max=day[j]-day[i];
            start=i,over=j;
        }
    }
}
cout<<"最大收益值："<<Max<<endl;
cout<<"买入日："<<start<<" "<<"卖出日："<<over;

算法Ⅱ要经历双重循环，即时间复杂度为$O(n^2)$。开辟一维数组，空间复杂度为$O(n)$。

1.2.3.算法Ⅲ、Ⅳ的分析

（1）算法Ⅳ描述

在本题中我们以算法Ⅳ为例，叙述算法思想。将题目抽象成数学模型，则为：对于给定的一个长度为$n$的序列$a$，找到一对点对$(i,j)$使得$a_j-a_i$的值最大。

这个问题是CDQ分治的一种典型例子，算法流程如下：

1、找到这个序列的中点$mid$。

2、将所有点对$(i,j)$划为三类：

• $1\le i\le mid,1\le j\le mid$的点对；
• $1\le i\le mid,mid+1\le j\le n$的点对；
• $mid+1\le i\le n,mid+1\le j\le n$的点对。

3、将$(1,n)$这个序列拆成两个序列$(1,mid)$和$(mid+1,n)$。此时第一类点对和第三类点对都在这两个序列中。

4、递归地处理第一类和第三类点对。

5、通过在每一层递归中处理第二类点对。

（2）算法Ⅲ、Ⅳ时空复杂度

算法Ⅲ（顺序对分治）和算法Ⅳ（直接分治）满足递归方程

$T(n)= \begin{cases} O(1),\qquad \qquad \,n<2\\ 2T(\frac{n}{2})+O(n),n\ge2 \end{cases}$

当$n\ge2$时，有$T(n)=n+2T(\frac{n}{2})=n+2(\frac{n}{2})+4T(\frac{n}{4})=...=nlogn$。即算法Ⅲ、Ⅳ的时间复杂度为$O(nlogn)$。

空间复杂度为$O(n)$。

（3）实例应用

例如，对于实例$[1,2,5,7]$，在第一层递归中分为$[1,2]$和$[5,7]$，左半部分第一类点对得到的最大收益是1，右半部分第三类点对得到的最大收益是2。再通过两边循环得到的第二类点对最大收益为7-1=6。则合并所有情况取最大值，最大收益为6。

1.3.算法Ⅳ运行结果

（1）

Input:

5

10 11 7 10 6

Output:

最大收益值：3

买入日：第3天 卖出日：第4天

（2）

Input:

10

4 23 7 23 3 4 34 23 10 9

Output:

最大收益值：31

买入日：第5天 卖出日：第7天

1.4.算法Ⅳ代码附录

代码附录链接：https://paste.org.cn/BfrkIv8CW7

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

#include<bits/stdc++.h>
using namespace std;

const int N=1e5+7;
const int inf=0x3f3f3f3f;
int a[N];

int max_profit_dc(int l, int r,int& Min,int& Max){
	
    if(l==r){
    	Min=l,Max=r;
        return 0;
    }
    
    int mid=(l+r)>>1;
	int x1_Min=0,x1_Max=0,x3_Min=0,x3_Max=0,x2_Min=inf,x2_Max=-inf;
	
    int x1=max_profit_dc(l,mid,x1_Min,x1_Max); //递归处理第一类点对的情况 
    int x3=max_profit_dc(mid+1,r,x3_Min,x3_Max); //递归处理第三类点对的情况 
    
    Min=mid,Max=mid+1;
    
    for(int i=mid+1;i<=r;i++){ //找到mid~r最大的值 
    	if(a[i]>x2_Max){
    		Max=i;
    		x2_Max=a[i];
		}
	}
    for(int i=l;i<=mid;i++){ //找到l~mid最小的值 
    	if(a[i]<x2_Min){
    		Min=i;
    		x2_Min=a[i];
		}
	}
    int x2=x2_Max-x2_Min; //在每一层递归中单独处理第二类点对的情况 
    if(x2<x3){ //如果第二类的最大收益值小于第三类 
    	Min=x3_Min;
    	Max=x3_Max;
	}
	if(a[Max]-a[Min]<x1){ //如果第二、三类的最大收益值小于第一类 
		Max=x1_Max;
		Min=x1_Min;
	}
    return max(max(x1,x3),x2);
}
int main(){
	int n;
	cout<<"请输入这只股票价格序列的长度"<<endl;
	cin>>n;
	cout<<"请分别输入这只股票不同日期的价格"<<endl;
	for(int i=1;i<=n;i++)
		cin>>a[i];
	int Max=0,Min=0;
	cout<<"最大收益值："<<max_profit_dc(1,n,Min,Max)<<endl;
	cout<<"买入日：第"<<Min<<"天 "<<"卖出日：第"<<Max<<"天"<<endl;
}

2.算法作业二

2.1.1.问题描述：

给定两个字符串$x$和$y$，如果要将字符串$x$替换为$y$，转换过程有三个基本操作：插入（Ins）字符C、删除（Del）字符C、将字符C替换（Rep）为字符C’。如何组合这些操作从而使得转换的代价最小？

2.1.2.算法设计：

（1）本问题是一个优化问题，可采用动态规划法来求解。

（2）**该问题在文本编辑、基因检测等领域有重要应用。**对于单纯的字符串转换代价问题，可假设单次插入、删除、替换操作具有相同的代价（=1）。

（3）定义子问题：

考虑将问题的后缀作为子问题，也即分别取输入字符串$x$和$y$的后面部分作为子问题。不妨设函数DP是用于求解该问题的策略，则子问题可表示为DP$(x[i:],y[j:])$，其中，$0≤i≤len(x)$, $0≤j≤len(y)$，$len()$表示求长度的函数。易知子问题个数为$len(x)×len(y)$。

（4）建立子问题的递归关系。

如图所示子问题的结构，为了将字符串$x$转换为$y$，对于字符$x[i]$有以下三种可能的操作：

• 删除字符$x[i]$。如图b：DP$(x[i:],y[j:])=$DP$(x[i+1:],y[j:])+cost(Delx[i])$；

• 插入字符$y[j]$。如图c：DP$(x[i:],y[j:])=$DP$(x[i:],y[j+1:])+cost(Insy[j])$；

• 替换字符$x[i]$为$y[j]$。如图d：DP$(x[i:],y[j:])=$DP$(x[i+1:],y[j+1:])+cost(Repy[j])$。

其中，$cost()$函数为执行一次操作的代价，Del、Ins、Rep表示删除、插入、替换操作。于是，DP$(x[i:],y[j:])$就等于以上三者中的最小值。由此可建立完整的递归关系（从略）。递归初始条件为DP$(len(x),len(y)) = 0$。

2.1.3.作业要求

（1）理解以上动态规划算法思想，写出完整的递归关系式。

（2）对输入x=“book”和y=“broke”，利用上述动态规划思想，人工求解将$x$转换为$y$的最小转换代价（可列表求解）。

（3）将以上算法付诸代码实现，给出三个输入实例的分析结果，并附代码。

2.2.算法分析

（1）算法描述

不妨记执行一次操作的代价为一个单位，对于当前的$x_i$和$y_j$，子问题的递归关系统共有四种情况：

（1）当$x_i$等于$y_j$时，有$dp_{i,j}=dp_{i+1,j+1}$

（2）当$x_i$不等于$y_j$时，有三种可能

• 在字符串$x$中删除$x_i$，则有$dp_{i,j}=dp_{i+1,j}+1$
• 在字符串$x$中插入$y_j$，相当于在字符串$y$中删除$y_j$，则有$dp_{i,j}=dp_{i,j+1}+1$
• 在字符串$x$中将$x_i$替换为$y_j$，相当于在字符串$x$中删除$x_i$后又在字符串$x$中插入$y_j$，则有$dp_{i,j}=dp_{i+1,j+1}+1$

同时，也要注意边界情况：当字符串$x$或$y$为空的时候，后缀子问题的值应当为子问题的长度。即$dp_{i,len(y)}=len(x)-i$，或$dp_{len(x),j}=len(y)-j$。

综合所有情况，得到后缀子问题的状态转移方程为：

$\begin{cases} dp_{i,len(y)}=len(x)-i,\qquad \qquad \qquad \qquad \qquad j=len(y)\\ dp_{len(x),j}=len(y)-j,\qquad \qquad \qquad \qquad \qquad i=len(x)\\ dp_{i,j}=dp_{i+1,j+1},\qquad \qquad \qquad \qquad \qquad \qquad \;x_i=y_j\\ dp_{i,j}=min\{dp_{i+1,j},\;dp_{i,j+1},\;dp_{i+1,j+1}\}+1,\;\;\;else \end{cases}$

（2）时空复杂度

算法涉及两层循环嵌套，时间复杂度为$O(nm)$。空间复杂度为$O(n+m)。$其中$n=len(x)$，$m=len(y)$。

（3）实例应用

人工求解将$x=book$转换为$y=broke$的最小转换代价，列表如下。其中箭头表示当前数值从哪个状态转移（转移方向不唯一）。

2.3.运行结果

以下是三个实例的测试结果。

（1）

Input:

$x=aaaa,\; y=aa$

Output:

2

状态转移列表如下，答案出口为$dp_{0,0}$，即【】中所示。

（2）

Input:

$x=Jingyi,\;y=Jophie$

Output:

5

状态转移列表如下，答案出口为$dp_{0,0}$，即【】中所示。

（3）

Input:

$x=thank,\;y=thinking$

Output:

4

状态转移列表如下，答案出口为$dp_{0,0}$，即【】中所示。

2.4.代码附录

代码附录链接：https://paste.org.cn/6bwe4lKfu4

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#include<bits/stdc++.h>
using namespace std;
int main(){
	string x,y;
	cin>>x>>y;
	vector<vector<int>>dp(x.size()+1,vector<int>(y.size()+1,0));
	for(int i=0;i<=x.size();i++) dp[i][y.size()]=x.size()-i;
	for(int j=0;j<=y.size();j++) dp[x.size()][j]=y.size()-j;
	for(int i=x.size()-1;i>=0;i--){
		for(int j=y.size()-1;j>=0;j--){
			if(x[i]==y[j]){
				dp[i][j]=dp[i+1][j+1];
			} else{
				dp[i][j]=min(dp[i+1][j+1],min(dp[i+1][j],dp[i][j+1]))+1;
			}
		}
	}
	cout<<dp[0][0]<<endl;
}

3.算法作业三

3.1.1.问题描述：

假设零钱系统的币值为$\{1, p, p^2, … , p^n\},\;p>1$，且每个钱币的重量都等于1，设计一个最坏情况下时间复杂度最低的算法，使得对任何钱数$y$，该算法得到的零钱个数最少。

3.1.2.作业要求：

（1）说明该问题的贪心算法策略，并举出一个实例。

（2）参考教材定理4.6，试简单证明贪心策略的最优性。

3.2.算法分析：

（1）贪心策略

将本题题意抽象为数学模型，则对于序列$\{1, p, p^2, … , p^n\}, p>1$，有约束函数：

$\sum_{i=0}^mx_ip^i=y,\qquad x_i\in N$

且其目标函数为

$min\{\sum_{i=1}^mx_i\}$

假如对进制转换足够敏感的话，很容易发现这个问题的本质上是将$y$表示成$p$进制后，将每个数位加起来的和就是最少的零钱个数。即按照币值从大到小排列零钱，从币值大的$p^k$开始贪心，每种钱尽量多用。如果剩余零钱小于该币值，再考虑用下一种钱币。

（2）实例说明

对于$y=58,p=7$，则用$p$进制来表示$y$则有$y=112_{(7)}$，则零钱个数$w=1+1+2=4$。

如果用贪心的思想，即对于$y=58,p=7$，首先取一个$p^2=49$，剩下$58-49=9$元；再取一个$p^1=7$，剩下$9-7=2$元；再取两个$p^0=1$元，正好取完。则有$p^2+p^1+2p^0=y$，最少的零钱个数为$w=1+1+2=4$。

（3）贪心策略的证明

记$F_k(y)$表示在动态规划算法（最优解）中用前$k$种硬币，总钱数为$y$的最小重量。记$G_k(y)$表示在贪心算法种用前$k$种硬币，总钱数为$y$的最小重量。

下面通过数学归纳法来证明对于任意$k$，有$F_k(y)\equiv G_k(y)$。

（1）当$k=1$时显然成立。

（2）假设$k\gt1$时成立，即$G_k(y)=F_k(y)$且存在$\theta=p,\;\delta=0$​，则有

$\theta p^k-\delta=p^{k+1}$

且满足$0\le\delta\lt p,\theta\in N$，则由教材定理4.6可得

$G_{k+1}(y)=F_{k+1}(y)$

即$k+1$时也成立。则对于任意的$k$，均有$F_k(y)\equiv G_k(y)$，证毕。

（4）时空复杂度

时间复杂度为$O(n)$，其中$n$为$y$位数的长度。空间复杂度为$O(1)$，仅为常数级别。

4.算法作业四

4.1.1.问题描述：

许多计算机科学问题可利用图来建立模型，且最后问题的求解过程往往就是在图上进行遍历和搜索。

4.1.2.作业要求：

（1）请查阅有关资料，举1~2个可转换为图搜索问题的应用场景，并简要说明其中图搜索的机制（不少于200字）。

（2）回溯法和分支限界方法也常在图的搜索算法中。请简述回溯法和分支限界方法的算法思想和算法设计步骤。

4.2.图搜索的应用场景与搜索机制

图的搜索其实可以应用在大学选课的场景中。比如大学课程中有：《微积分》，《线性代数》，《离散数学》，《概率论与统计学》，《算法导论》，《机器学习》等等。当我们想要学习《机器学习》的时候，就必须先学会《微积分》和《概率论与统计学》。这些课程就相当于几个顶点$u$, 顶点之间的有向边$(u,v)$就相当于学习课程的顺序。这样的选课场景我们通常称其解决方案为”拓扑排序“。下面将介绍如何将这个过程抽象出来，用算法来实现。

但是如果某一天排课的老师打瞌睡了，说想要学习《算法导论》，还得先学《机器学习》，而《机器学习》的前置课程又是《算法导论》。在这里，《算法导论》和《机器学习》间就出现了一个闭环，显然我们现在没办法弄清楚到底需要学什么了，于是我们也没办法进行拓扑排序了。因而如果有向图中存在环路，那么我们就没办法进行拓扑排序了。即拓扑排序还有一个应用场景是，在有向图中判环。

因而将场景抽象出数学模型，即为我们在一个有向无环图中，将图中的顶点以线性的方式进行排序，使得对任何的顶点$u$到$v$的有向边$(u,v)$ ，都有$u$在$v$的前面。

我们可以通过广度优先搜索或者深度优先搜索来实现拓扑排序。广度优先的思路其实就是对每个入度为0的且未被访问过的节点进行广度优先搜索。

在搜索过程中，只要搜索了$u$与$v$之间的边，那就将$v$的入度减1，相当于删除边的操作。入度为零就代表着它的前置条件已经完全满足。一个节点只有当其入度为0且未被访问过，才对其进行广度优先搜索。具体算法流程请看代码。

代码提取链接：https://paste.org.cn/wDu2ox1MNW

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

#include<bits/stdc++.h>
using namespace std;

const int N=1e5+7;

vector<int>g[N];
vector<int>ans;
int v,e,indeg[N];
bool vis[N];

void bfs(int u){
    queue<int>q;
    vis[u]=true;
    q.push(u);
    while(!q.empty()){
        int t=q.front();
        q.pop();
        ans.push_back(t);
        for(auto x:g[t]){
            indeg[x]--;//将入度减一
            if(!indeg[x]&&!vis[x]){//从入度为0且未被访问过的节点压入队列中
                q.push(x);
                vis[x]=true;
            }
        }
    }
}

int main(){
    cin>>v>>e;
    for(int i=0;i<e;++i){
    	int u,v;
        cin>>u>>v;
        g[u].push_back(v);//邻接表建有向图
        indeg[v]++;//记录入度
    }
    for(int i=0;i<v;++i){
        if(!indeg[i]&&!vis[i])
            bfs(i);//从入度为0且未被访问过的节点开始搜索
    }
    for(int x:ans)
        cout<<x<<endl;
}

4.3.1.回溯 & 分支限界的思想

在屈婉玲教授的《算法设计与分析》中，回溯与分支限界法都是基于深度优先搜索的基础上进行优化改进的。回溯算法实际上一个类似枚举的搜索尝试过程，主要是在搜索尝试过程中寻找问题的解，当发现已不满足求解条件时，就“回溯”返回，尝试别的路径。回溯法是一种优选搜索法，按选优条件向前搜索，以达到目标。但当探索到某一步时，发现原先选择并不优或达不到目标，就退回一步重新选择，这种走不通就退回再走的技术为回溯法，而满足回溯条件的某个状态的点称为“回溯点”。分支限界法则是在回溯算法的基础上，设定一个代价函数对搜索分支进行剪枝优化，从而降低算法的复杂度。

4.3.1.回溯 & 分支限界的算法设计步骤

（1）回溯 & 分支限界的一般步骤：

1、针对所给问题，定义问题的解空间；

2、确定易于搜索的解空间结构；

3、以深度优先方式搜索解空间，并在搜索过程中确定界函数避免无效搜索。

（2）回溯 & 分支限界的实例

以TSP旅行商问题为例，下图为一个实例。

其回溯法的解空间为：

分支限界法的解空间为：

不妨设巡回路线从1开始，解向量为$<i_0=1,\ i_1,\ldots,i_{n-1}>$，其中$ i_1,\ldots,i_{n-1}$为${2,\ldots,n}$的一个排列。搜索树为上图所示。其中代价函数为$F$，界函数（当前得到的最短巡回路线长度)为$B$，且有

$F=\sum_{j=1}^{k}c_j+l_{i_k}+\sum_{i_j\notin B} l_{i_j}$

由此，我们便确定好了解空间和剪枝函数（界函数）。以下为分支限界法的实例代码，分支限界与回溯的不同之处主要在于多了一个界函数的剪枝优化。

代码提取链接：https://paste.org.cn/hqffCgqfri

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

#include<bits/stdc++.h>
using namespace std;
const int kMax=107;
vector<pair<int,int>>g[kMax]; 
vector<vector<int>>ans;
vector<int>temp;
vector<int>vis,l;
int n,m,B=0x3f3f3f3f;
void dfs(int fa,int sum,int count){
    if(count==n){
        if(fa==1&&sum<=B){
            ans.push_back(temp);//存储当前最优解 
            B=sum;//更新界函数 
        }
        return ;
    }
    for(auto i:g[fa]){
       if(!vis[i.first]){
            int F=sum+l[fa];//求解代价函数F
            for(int i=2;i<=n;i++){
                if(!vis[i]){
                    F+=l[i];
                }
            }
            if(F>B){ //如果代价函数F大于界函数B，则回溯到父节点。 
                return ;
            }
            vis[i.first]=true;
            temp.push_back(i.first);
            dfs(i.first,sum+i.second,count+1);
            temp.pop_back();//回溯
            vis[i.first]=false;
        }
    }
}
int main(){
    cin>>n>>m;
    vis.resize(n+1,false);
    l.resize(n+1,0x3f3f3f3f);
    while(m--){
        int u,v,w;cin>>u>>v>>w;
        g[u].push_back({v,w});//建图 
        g[v].push_back({u,w});
    }
    for(int i=1;i<=n;i++){
        for(auto it:g[i]){
            l[i]=min(l[i],it.second);//预处理每个节点所连边的最小权值
        }
    }
    temp.push_back(1);
    dfs(1,0,0);//搜索入口 
    cout<<"最短巡回路径有"<<ans.size()<<"条，分别为"<<endl; 
    for(auto i:ans){
        for(auto it:i){
            cout<<it<<" ";
        }cout<<endl;
    }
    cout<<"最短巡回路径长度为"<<B<<endl;
}
/*

4 6
1 2 5
1 3 9
1 4 4
3 4 7
2 4 2
2 3 13

*/

5.算法作业五

5.1.1.问题描述

根据高等数学知识，积分$\int{e^{x^2}dx}$没有初等函数形式的原函数。但对于定积分$\int_0^1{e^{x^2}dx}$，我们可以通过设计随机数值方法计算其积分的近似值。

5.1.2.作业要求

（1）采用线性同余法，编写一个随机数生成的程序 。

（2）应用自编随机数生成程序，分别在平面区域$[0,1]×[0,e]$内取$100、1000、10000、100000$个随机点，设计随机数值算法计算定积分$\int_0^1{e^{x^2}dx}$的近似值，并列表分析近似值的收敛性。

5.2.算法分析

由线性同余法在$O(n)$的时间复杂度产生$n$个随机数序列，其核心代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

#include<bits/stdc++.h>
using namespace std;

unsigned long long seed=chrono::system_clock::now().time_since_epoch().count();
mt19937 rnd(seed);

const int mod=1e9+7,B=rnd(),C=rnd();
int main(){
	int n;
	cin>>n;
	unsigned long long x=(B*rnd()+C)%mod;
	for(int i=1;i<=n;i++){
		x=(B*x+C)%mod;
		cout<<x<<" ";
	}
}

本算法计算定积分$\int_0^1{e^{x^2}dx}$的核心思想是蒙特卡罗法。要得到平面区域$I=[0,1]×[0,e]$内的随机点，则对于随机数$x$，可做变换$x^*=\frac{x}{mod}$，$y^*=\frac{xe}{mod}$，从而得到随机点$(x^*,y^*)$。我们记落在积分区域内点的数量为$cnt$，总的点个数为$n$，则得到估算的积分近似值为$\frac{cnt}{n}$。

由自适应辛普森积分，我们首先得到精确值$\int_0^1{e^{x^2}dx}≈1.462651745907$。在本题，我们通过列表法比较近似值来的绝对误差和方差来评价近似值的精确性。

5.3.运行结果 & 结果分析

我们发现，随着$n$的增大，$\epsilon$更加小，即精确性越高，随机算法更加稳定。

同时，根据大数定律，我们也了解到，由于每次测试的简单样本都是独立同分布的，随着测试次数的增多，算术平均结果更加趋近于精确值。

5.4.代码附录

代码提取链接：https://paste.org.cn/v23DKWNouR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

#include<bits/stdc++.h>
#define double long double
using namespace std;

unsigned long long seed=chrono::system_clock::now().time_since_epoch().count();
mt19937 rnd(seed);

const unsigned long long mod=1e9+7;
const unsigned long long B=rnd(),C=rnd();
const double eps=1e-6;
const double E=2.718281828459;

int main(){
	int n,cnt=0;
	double x,y;
	cin>>n;
	unsigned long long x1=(B*rnd()+C)%mod,y1=(B*rnd()+C)%mod;
	for(int i=1;i<=n;i++){
		x=x1*1.0/mod;
		y=y1*E/mod;
		if(pow(E,pow(x,2))-y>eps)
			cnt++;
		x1=(B*x1+C)%mod;
		y1=(B*y1+C)%mod;
	}
	cout<<fixed<<setprecision(18)<<cnt*1.0/n*E<<endl;
}