背景
圆周率被广泛应用在众多的实际场景中,然而圆周率的精确值一直是数学界悬而未决的问题。随着计算机的广泛普及,使用计算机去求解圆周率吸引了越来越多的关注。
代码运用了蒙特卡洛法、割圆法、无穷级数法(格雷戈里-莱布尼茨级数)、梅钦法。
精度表示
误差的绝对值取对数再取负数
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;
}