圆周率Pi的求解

字节跳动面试时，算法部分给我丢了这个开胃小菜；当时使用的是蒙特卡洛模拟，效果还不错。

1. 蒙特卡洛模拟

蒙特卡洛模拟用于模拟一个随机过程中不同结果可能出现的概率。那怎么使用蒙特卡洛模拟来求解圆周率Pi呢？

1.1 思想

蒙特卡洛模拟的思想是很简单的，比如：在以[0,0]，[1,1]为顶点的正方形内进行随机投点，如果该随机点到圆心[0,0]的距离小于1，则说明该点位于以[0,0]为圆心，1为半径的1/4圆内；当试验次数达到无穷时，位于圆内点的数目和总投点数目之比无限接近圆和正方形面积之比，即pi/4。

蒙特卡洛模拟
蒙特卡洛模拟

1.2 求解代码

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 <iostream>
#include <cstdlib>
#include <ctime>
#include <vector>

using namespace std;

class Solution{
public:
    float calcPi(const int numOfPoints){
        // 计算随机点到原点的距离平方
        auto getDist = [](const vector<float>& point){
            return point[0] * point[0] + point[1] * point[1];
        };
        int legal_num = 0;
        for(int i = 1; i<= numOfPoints; ++i){
            vector<float> point = getNextRandom();
            if(getDist(point) <= 1.0f){
				++legal_num;
            }
        }
        // 面积之比： Pi/4 ≈ legal_num/numOfPoints
        return 4.0*(float)legal_num/(float)numOfPoints;
    }
    vector<float> getNextRandom(){
    	int x = rand() % 100 + 1;
        int y = rand() % 100 + 1;
        return {x/100.0f, y/100.0f};
    }
};

int main(){
    srand(time(0));
    Solution s;
    const int numOfPoints = 10* 10000;
    std::cout << s.calcPi(numOfPoints) << std::endl;
    return 0;
}

/*
	out:
	3.11112
*/

从这里的输出结果可以看到，使用蒙特卡洛求解的结果精度明显不够，不过思想简单算法实现容易。

我当时面试时，输出的结果差不多就是这个样子；面试官看后说，有没有办法使得结果更加精确一些呢？

我说：可以增大样本点的个数；或者实验多组数据求取平均值；但是感觉这个模拟方法不太好。

面试官又问：还有优化的方法么？这个随机数的生成可以优化一下么？

我想：随机数的生成？什么意思？难道是在说随机数的底层生成算法上解决？虽然有接触过线性同余法，但是这水平更不不够呀。难道是说，想要保证随机数不重复？使用哈希？没想明白，就说不会了。

有想法的朋友可以告知我一下哈~

2. 级数求解

在高中或者在大学数学中，我们可能会看到这样一个公式：（通过arctanx级数求和得到）

Pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + 1/13 ...

其实这是一个交错级数（以0为极限，必收敛），是数学家莱布尼茨发现的计算圆周率Pi的公式，它的通项表达为如下：

莱布尼茨公式
莱布尼茨公式

2.1 代码

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
#include <iostream>
using namespace std;

class Solution{
public:
    double calcPi(const unsigned int n){
        double sum = 0.0, symbol = 1.0;
        for(int i = 0; i<=n; ++i){
            sum += symbol/(2 * i + 1);
            symbol = -symbol;
        }
        return 4 * sum;
    }
};

int main(){
    Solution s;
    const unsigned int n = 10 * 10000;
    std::cout << s.calcPi(n) << std::endl;
    return 0;
}

/*
	out:
	3.1416
*/

从输出结果可以看出，使用莱布尼茨级数来计算Pi值精度更好，值得一提的是该算法的求解效率大大优于蒙特卡洛模拟。

好了，关于圆周率Pi的求解方法暂时介绍到这里，大家有什么更好的方法，可以在下面畅所欲言~