尝试用程序计算Π(3.141592653......)

1. $\pi$

$\pi$的重要性或者地位不用多说，有时候还是很好奇，精确地$\pi$值是怎么计算出来的。研究$\pi$的精确计算应该是很多数学家计算机科学家努力的方向，很多算式看起来和$\pi$没有一点关系，但是最后的表达式中却包含$\pi$。
作为一个国家定义的从事信息技术的“新时代农民工”，要验证某些式子相对来说就容易得多了。

2. 用微积分来计算$\pi$

2.1 原理

小学时候就知道了圆的面积公式$S=\pi r^2$,到了高中学习了微积分，方程与曲线之后可以知道：
$$x^2+y^2=r^2$$
表示一个圆，如果$y>=0$则表示上半圆，这个时候上半圆可以表达为一个函数:
$$y=\sqrt{(r^2-x^2)}$$
当它是一个单位圆时（r=1）图形如下所示：

根据从微积分的角度来看曲线$y(x)与x=0，x\in [-r,r]$围成的区域面积，就是圆面积的一半，即上半圆的面积。$令圆面积为S,则有：$
$$S=\pi r^2=2{\int }_{-r}^r\sqrt{(r^2-x^2)}dx$$
现在的工作就变成计算上面右侧的积分式了。积分式用化曲为直（马师傅称为：接 化 发），把圆分成许多直立的小矩形，再求和，小矩形的高为$y_i=\sqrt{(r^2-x_i^2)},宽为\Delta x,\Delta x=\frac{1}{n}$,n表示把圆的面积分为n个小矩形的面积来近似，n越大则近似程度越好，而且当n趋近于无穷时就和积分式完全一致了。
为了求$\pi$那我们只需要用程序帮我们把圆分成一个一个的小矩形再求和就行了，设定圆的半径，每个矩形的宽度，就可以计算积分式的近似值了。

下面的程序就来计算下式
$$\pi =\frac{2}{r^2}{\int }_{-r}^r\sqrt{(r^2-x^2)}dx$$

2.2 代码

#include <cmath>
#include <iostream>


/**
 * @brief calculate sqrt(r^2-x^2)
 * 
 * @param x  : 
 * @return constexpr long double 
 */
constexpr inline long double IntegralFormula(const long double radius,const long double x)
{
    return std::sqrt(radius*radius-x*x);
} 

/**
 * @brief 计算曲线y=sqrt(r^2-x^2),在 x in [-r,r]的积分
 * 
 * @param radius  : 
 * @param delta_x  : 
 * @return constexpr long double 
 */
constexpr long double Integral(const long double radius,const long double delta_x)
{
    long double sum = 0.0;
    for (double x = -radius; x <= radius; x += delta_x)
     {
        sum += delta_x * IntegralFormula(radius,x);
     }
     return sum;
}

/**
 * @brief 
 * 
 * @param radius  : 半径
 * @param bin_num  : 分成小矩形的数目
 * @return constexpr long double 
 */
constexpr long double CalculatePi(const long double radius,double bin_num)
{
    long double delta_x = 2 * radius / bin_num;
    long double integral = Integral(radius,delta_x);
    long double pi = integral * 2.0 / (radius * radius);

    return pi;
}

/**
 * @brief 计算pi值，用单位圆的上半圆做积分，y=sqrt(r^2-x^2),x in [-r,r]
 * 
 */
void Test()
{
    constexpr long double true_pi = 3.141592653589793;
    printf("True PI: %.15Lf\n",true_pi);
    for(double radius = 10; radius <= 100; radius*=10)
        for(double bin_num = 100000; bin_num <= 10000000; bin_num*=10)
        {
            printf("radius:%f, \t bin_num: %f,\t pi: %.20Lf \n",radius,bin_num,CalculatePi(radius,bin_num));
        }
    
}

void Test2()
{
    double radius = 10000;
    double bin_num = 100000;
    printf("radius:%f, \t bin_num: %f,\t pi: %.20Lf \n",radius,bin_num,CalculatePi(radius,bin_num));
}

int main()
{
    Test();
    Test2();
    return 0;
}

2.3 结果

输出

True PI: 3.141592653589793
radius:10.000000,        bin_num: 100000.000000,         pi: 3.14159254840509465393
radius:10.000000,        bin_num: 1000000.000000,        pi: 3.14159265028029712060 
radius:10.000000,        bin_num: 10000000.000000,       pi: 3.14159265331958147799 
radius:100.000000,       bin_num: 100000.000000,         pi: 3.14159254846500206543
radius:100.000000,       bin_num: 1000000.000000,        pi: 3.14159265024260935645 
radius:100.000000,       bin_num: 10000000.000000,       pi: 3.14159265331696711159 
radius:10000.000000,     bin_num: 100000.000000,         pi: 3.14159254840753097960

2.4 分析

这重方式用long double的计算精度，但是最后计算的精度最高只到"3.141592653",小数点后9位。分再多的小矩形也没用了，因为计算时候的截断误差已经足够影响计算精度了。不知道有没有什么好的策略保持计算精度的。