计算圆周率π可以通过随机数、泰勒展开、割圆术、投针术等方法,本节主要介绍前两项的c++实现。
方法一:随机数算法
即,蒙特卡罗法,是先在平面区域均匀的(随机的)撒点,并按照点分布在不同区域的数量,来代表不同区域面积(这样就可以模拟计算圆形区域的面积,从而可求π),用整数的比例关系来代替面积的比例关系。
如下图所示,
其中,半径r=1的圆的面积为: S圆
S圆 =π * r * r = π * 1 * 1 = π
圆的外接正方形(边长为2)的面积为:S正
S正 = 2 * 2 = 4
那么,我们向正方形内部随机撒点,假设这些点落入正方形内的数量为N(总的点数),落入圆内的点数为n。
我们有下面的比例关系:
随着撒点数量的增大,点的数量会越来越能够逼近其所在区域的面积。
上面的等式就变成了,n是已知的,N也是已知的,S正也是已知的,求S圆的值。
此时,
实现思路:
上图正方形内部的点的坐标,假设为(x, y)。
那么,x和y的值均为0~1之间的值。
我们不断的创造出0~1之间的小数,来模拟一个点在正方形中的坐标x的值
我们不断的创造出0~1之间的小数,来模拟一个点在正方形中的坐标y的值
假设我们一共撒点N个,通过x,y值,我们就可以计算出x,y这个点在圆的内部还是在圆的外部。
点(x,y)在圆的内部,就是:
所以,我们就可以计算出π的近似值了。
代码如下:
#include<iostream>
#include <ctime>
using namespace std;
int main()
{
double x, y, r, pi;
int n = 0, N = 1e7;
srand(time(NULL));
for (int i = 0; i < N; i++)
{
x = rand() / double(RAND_MAX);//取值为0-1之间的随机数
y = rand() / double(RAND_MAX);
r = sqrt(x * x + y * y);
if (r < 1)
{
//(1) your code
}
}
pi = n * 4.0 / N;
cout <<pi << endl;
}
随着N取值增大,通过程序模拟,得出如下数据
因为是随机数,所以同一个N值也可能输出不同的结果,精确度并不是随着N增大一直正相关。
方法二:泰勒展开
S0 = 1
S1 = S0×(1 /(2×1+1))= 1/3
S2 = S1×(2 /(2×2+1))= 1/3 × 2/5 = 2/15
S3 = S2×(3 / (2×3+1))= 2/15 × 3/7 = 6/105
代码如下:
#include <iostream>
#include <iomanip>
using namespace std;
int main()
{
cout.precision(20);//小数部分20位有效数字
int n = 100;//迭代次数
double S = 1;// S0=1
double Sum = S;
for (int m = 1; m < n; m++)
{
//下面的代码中int和double在一起运算,运算结果类型为double类型
S = m / (2 * m + 1.0) * S;
//(2) your code
}
double pi = 2 * Sum;
cout << "pi=" << pi << endl;
return 0;
}
输出如下:
圆周率的实际值为(20位小数):3.14159265358979323846
可以看出,迭代100次之后,小数点后面的14位都是正确的。
隐式类型转换
上面的泰勒展开求π的代码中:
S = m / (2 * m + 1.0) * S;
- 1.0 是一个 double 类型;
- 整数 m 是 int 类型;
- 2 * m 的结果是 int类型;
- 2 * m + 1.0 的结果是 double 类型;这里发生了隐式类型转换
- m / (2 * m + 1.0) 的结果是 double 类型; 这里再次发生了隐式类型转换
可以看出,隐式类型转换是很方便的。
但是,如果上面的1.0 写成了int类型的 1 ,
m / (2 * m + 1)
就会出问题。这会导致难以察觉的逻辑错误。
这会导致两个整数相除,分母始终大于分子,从而结果为 0。 参考整数除法。
如需答案和答疑:请私信。