圆周率的多种数值计算方法以及绘制数据曲线

背景
圆周率被广泛应用在众多的实际场景中，然而圆周率的精确值一直是数学界悬而未决的问题。随着计算机的广泛普及，使用计算机去求解圆周率吸引了越来越多的关注。

代码运用了蒙特卡洛法、割圆法、无穷级数法（格雷戈里-莱布尼茨级数）、梅钦法。

精度表示
误差的绝对值取对数再取负数

level_acc=−log⁡(|pi−π|)

level_acc用于输出精度、作为绘图坐标

最后排序按照level_acc从高到低保存在本地txt文件，以及保存从0到循环次数的所有数据到本地csv文件，从而可以画出精度-速度，和精度-循环次数的曲线，从而通过可视化的手段清晰地看出各种方法的优劣。

#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
#include <random>
#include <algorithm> // for std::sort
#include <vector>
#include <chrono>    // for timing
using namespace std;
using namespace std::chrono;

const double REAL_PI = 3.14159265358979323846;// 圆周率的真实值

struct PiMethodResult {
    //定义了一个结构体PiMethodResult，用于存储每种方法的名称、计算出的圆周率值、误差和执行时间。
    string methodName; //方法名字
    double calculatedPi; //计算结果
    double error; //误差
    long long executionTime; // 添加运行时间成员变量
};

// 割圆法
double Geyuan(int n, long long& duration) {
    //n表示迭代次数
    auto start = high_resolution_clock::now(); // 记录开始时间
    double tmp = 0.5, s = 1; 
    //初始化两个变量tmp和s。tmp用于存储每次迭代计算的中间值，初始值为0.5；s用于累积每次迭代的结果，初始值为1。
    for (int i = 1; i <= n; ++i) {
        //开始一个循环，从1迭代到n，每次迭代都会更新tmp和s的值。
        tmp = sqrt((tmp + 1) / 2);
        //在每次迭代中，计算tmp的新值。这是割圆法的一个关键步骤，通过迭代逼近圆的边长。
        s *= tmp;
    }
    double pi = 1.5 * sqrt(3) / s;
    //计算圆周率π的近似值。根据割圆法的公式，这里使用了1.5乘以根号3再除以s。
    auto stop = high_resolution_clock::now(); // 记录结束时间
    duration = duration_cast<microseconds>(stop - start).count(); // 计算执行时间
    return pi;
}

// 蒙特卡洛法
double montecarlo_pi(int epochs, long long& duration) { //epochs表示随机采样的次数
    auto start = high_resolution_clock::now();
    std::default_random_engine random_engine;//创建一个默认的随机数引擎random_engine，用于生成随机数。
    std::uniform_real_distribution<double> distribution(0, 1);
    //定义一个均匀分布的随机数生成器distribution，它会生成[0, 1]区间内的随机实数。
    int count = 0;
    for (int i = 0; i < epochs; i++) {
        //开始一个循环，从0迭代到epochs - 1，每次迭代都会生成两个随机数并计算它们是否落在单位圆内。
        double x = distribution(random_engine);
        double y = distribution(random_engine);
        double distance = sqrt(x * x + y * y);
        if (distance <= 1.0)
            count++;
    }
    //根据蒙特卡洛法的原理，圆周率π的近似值可以通过4* (落在圆内的点数 / 总点数)来估计
    double PI = 4.0 * static_cast<double>(count) / epochs;
    auto stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start).count();
    return PI;
}
// 莱布尼茨级数
double Leibniz(int n, long long& duration) {//n表示级数的项数
    auto start = high_resolution_clock::now();
    double pi = 0;//初始化一个变量pi，用来累积级数的各项，初始值为0。
    for (int i = 0; i < n; i++) {
        pi += 4.0 * pow(-1, i) / (2 * i + 1);
        //在每次迭代中，计算莱布尼茨级数的第i+1项，公式为4 * (-1)^i / (2 * i + 1)，并累加到pi中。这里使用了pow函数来计算(-1)^i。
    }
    auto stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start).count();
    return pi;
}

// 梅钦公式
double machin_pi(long L, long long& duration) { //L表示级数的项数
    auto start = high_resolution_clock::now();
    const long C = 10000;//定义一个常量C，用于后续计算中的乘除操作，这里取值为10000，用于将数值转换为四位数字。
    short flag;//定义一个标志变量flag，用于控制级数中正负项的交替。
    long xishu[2] = { 80, -956 }, div[2] = { 25, 57121 }, odd, shang_1, r_1, r_2, i = 0, j, t, tt, k = 0, len = L;
    //定义两个数组xishu和div，存储梅钦公式中的特定系数。
    L = L / 4 + 3;
    long* term = new long[L];
    long* pi = new long[L];
    while (i < L) pi[i++] = 0;
    while (k - 2) {
        flag = 1;
        term[0] = xishu[k];
        i = 1;
        while (i < L) {
            term[i] = 0;
            ++i;
        }
        for (i = 0, odd = 1; i < L; ++i) {
            j = i;
            for (r_1 = r_2 = 0;j < L; ++j) {
                t = r_1 * C + term[j];
                term[j] = t / div[k];
                r_1 = t % div[k];
                tt = r_2 * C + term[j];
                shang_1 = tt / odd;
                r_2 = tt % odd;
                pi[j] += (flag ? shang_1 : -shang_1);
            }
            if (term[i])--i;
            odd += 2;
            flag = !flag;
        }
        ++k;
    }
    //这部分代码实现了梅钦公式的计算逻辑，通过迭代计算term和pi数组的值。
    k = 0;
    while (--i) {
        pi[i] += k;
        if (pi[i] < 0) {
            k = pi[i] / C - 1;
            pi[i] %= C;
            pi[i] += C;
        }
        else {
            k = pi[i] / C;
            pi[i] %= C;
        }
    }
    double PI = 3.0;
    for (i = 1; i < len / 4 + 1; ++i) {
        PI += pi[i] / pow(10, 4 * i);
    }
    auto stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start).count();
    delete[] term;
    delete[] pi;
    //释放之前动态分配的内存。
    return PI;
}

int main() {
    // 获取用户输入的循环次数
    int numIterations;
    cout << "请输入计算圆周率的最大循环次数：";
    cin >> numIterations;

    // 打开CSV文件进行写入
    ofstream outFileCsv("C:\\Users\\opiumpunk\\Desktop\\output.csv");
    if (!outFileCsv.is_open()) {
        cerr << "无法打开文件 output.csv" << endl;
        return 1;
    }

    // 写入CSV文件表头
    outFileCsv << "循环次数,割圆法_圆周率计算结果,割圆法_误差,割圆法_运行时间 (微秒),"
        << "蒙特卡洛法_圆周率计算结果,蒙特卡洛法_误差,蒙特卡洛法_运行时间 (微秒),"
        << "莱布尼茨级数_圆周率计算结果,莱布尼茨级数_误差,莱布尼茨级数_运行时间 (微秒),"
        << "梅钦公式_圆周率计算结果,梅钦公式_误差,梅钦公式_运行时间 (微秒)" << endl;

    // 打开TXT文件进行写入
    ofstream outFileTxt("C:\\Users\\opiumpunk\\Desktop\\output.txt");
    if (!outFileTxt.is_open()) {
        cerr << "无法打开文件 output.txt" << endl;
        return 1;
    }

    // 写入TXT文件表头
    outFileTxt << "按误差从高到低排序的结果：" << endl;
    outFileTxt << "--------------------------" << endl;

    for (int i = 0; i <= numIterations; ++i) {
        long long duration;

        // 割圆法
        duration = 0;
        double pi1 = Geyuan(i, duration);
        double error1 = -log(fabs(pi1 - REAL_PI));
        outFileCsv << i << "," << pi1 << "," << error1 << "," << duration << ",";

        // 蒙特卡洛法
        duration = 0;
        double pi2 = montecarlo_pi(i, duration);
        double error2 = -log(fabs(pi2 - REAL_PI));
        outFileCsv << pi2 << "," << error2 << "," << duration << ",";

        // 莱布尼茨级数
        duration = 0;
        double pi3 = Leibniz(i, duration);
        double error3 = -log(fabs(pi3 - REAL_PI));
        outFileCsv << pi3 << "," << error3 << "," << duration << ",";

        // 梅钦公式
        duration = 0;
        double pi4 = machin_pi(i, duration);
        double error4 = -log(fabs(pi4 - REAL_PI));
        outFileCsv << pi4 << "," << error4 << "," << duration << endl;

        // 如果是输入的循环次数，保存到TXT文件中
        if (i == numIterations) {
            outFileTxt << "计算次数: " << i << endl;
            outFileTxt << "割圆法:\t圆周率计算结果 = " << fixed << setprecision(10) << pi1
                << ", 误差 = " << fixed << setprecision(10) << error1 << ", 运行时间 = " << duration << " microseconds" << endl;
            outFileTxt << "蒙特卡洛法:\t圆周率计算结果 = " << fixed << setprecision(10) << pi2
                << ", 误差 = " << fixed << setprecision(10) << error2 << ", 运行时间 = " << duration << " microseconds" << endl;
            outFileTxt << "莱布尼茨级数:\t圆周率计算结果 = " << fixed << setprecision(10) << pi3
                << ", 误差 = " << fixed << setprecision(10) << error3 << ", 运行时间 = " << duration << " microseconds" << endl;
            outFileTxt << "梅钦公式:\t圆周率计算结果 = " << fixed << setprecision(10) << pi4
                << ", 误差 = " << fixed << setprecision(10) << error4 << ", 运行时间 = " << duration << " microseconds" << endl;
        }
    }

    // 关闭文件
    outFileCsv.close();
    outFileTxt.close();

    cout << "结果已保存到 output.csv 和 output.txt 文件中。" << endl;

    return 0;
}