当前位置：首页 > news >正文

C++使用拉玛努金公式计算π的值

news 2025/9/14 5:14:55

拉玛努金公式是一个快速收敛计算π值的公式，公式如下：

$\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum_{k=0}^{\infty} \frac{(4k)!(1103 + 26390k)}{(k!)^4 \cdot 396^{4k}}$

根据以上公式，我们可以使用C++编写程序计算π的值，考虑到普通计算机的性能问题，我初步用该程序计算小数点后的2000位。

1. 分析步骤

包含必要的头文件，并使用GMP的高精度数值类型（mpf_class用于浮点数，mpz_class用于整数）。
设置GMP的默认浮点精度（二进制位数），我们设置为7000位二进制精度（大约相当于2000位十进制精度）。
定义计算阶乘和整数幂的辅助函数（因为GMP没有直接提供mpz_class的幂函数，我们使用循环实现）。
初始化变量：总和total_sum，当前π的近似值pi_approx，以及用于收敛检查的变量。
迭代计算级数的每一项，直到连续3次迭代小数点后2000位不再变化。
每次迭代中，计算当前项，并更新总和和π的近似值。
将π的近似值转换为字符串，并比较小数点后2000位是否连续3次相同。
输出结果和计算时间。

注意：拉马努金公式收敛非常快，通常只需要很少的迭代次数（k=0,1,2,...）即可达到很高的精度。

由于阶乘和幂运算会随着k增大而急剧增大，我们使用GMP的多精度整数来避免溢出。

由于需要使用GMP库，以ubuntu为例，先安装相关的库。

sudo apt install libgmp-dev libmpfr-dev

2. 头文件引入和命名空间

#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <gmpxx.h>
#include <cmath>using namespace std;
using namespace std::chrono;

<iostream>: 用于输入输出操作
<iomanip>: 用于格式化输出
<string>: 用于字符串操作
<chrono>: 用于计时功能
<gmpxx.h>: GMP库的C++接口，提供高精度数学运算
<cmath>: 标准数学函数库

3. 辅助函数

阶乘计算函数

mpz_class factorial(int n) {mpz_class result = 1;for (int i = 1; i <= n; ++i) {result *= i;}return result;
}

使用GMP的mpz_class类型处理大整数阶乘
通过循环乘法计算阶乘，避免使用递归可能导致的栈溢出
返回mpz_class类型，可以处理非常大的整数

整数幂计算函数

mpz_class int_pow(mpz_class base, int exponent) {mpz_class result = 1;for (int i = 0; i < exponent; ++i) {result *= base;}return result;
}

自定义整数幂函数，因为GMP库没有直接提供mpz_class的幂函数
通过循环乘法计算幂，适用于整数指数
同样返回mpz_class类型，处理大数运算

4. 主函数

精度设置

const int precision = 7000;
mpf_set_default_prec(precision);

设置GMP浮点数的默认精度为7000位二进制精度
7000位二进制精度 ≈ 7000 / log₂(10) ≈ 2100位十进制精度，足够2000位需求
mpf_set_default_prec设置所有后续创建的mpf_class对象的默认精度

初始化变量

mpf_class pi_approx(0, precision);
mpf_class total_sum(0, precision);
string last_pi_str;
int k = 0;
int convergence_count = 0;
const int convergence_threshold = 3;

pi_approx: 存储当前π的近似值
total_sum: 存储拉马努金公式中的级数和
last_pi_str: 存储上一次迭代的π值字符串表示
k: 迭代计数器
convergence_count: 收敛计数器，记录连续多少次迭代结果相同
convergence_threshold: 收敛阈值，连续3次相同则认为收敛

常量计算

mpf_class sqrt2(0, precision);
mpf_sqrt_ui(sqrt2.get_mpf_t(), 2);
mpf_class C = 2 * sqrt2 / 9801;

计算拉马努金公式中的常数部分 (2√2)/9801
mpf_sqrt_ui计算2的平方根，使用GMP库的高精度平方根函数
存储在C变量中供后续使用

主循环

while (true) {// 计算当前项mpz_class numerator = factorial(4 * k) * (1103 + 26390 * k);mpz_class denominator = int_pow(factorial(k), 4) * int_pow(396, 4 * k);mpf_class term = mpf_class(numerator) / mpf_class(denominator);// 添加到总和total_sum += term;// 计算当前π的近似值pi_approx = 1 / (C * total_sum);// 转换为字符串并检查收敛// ...
}

分子计算: factorial(4 * k) * (1103 + 26390 * k)
- 计算拉马努金公式中当前项的分子部分
分母计算: int_pow(factorial(k), 4) * int_pow(396, 4 * k)
- 计算拉马努金公式中当前项的分母部分
- 使用自定义的int_pow函数计算幂
项计算: 将分子除以分母得到当前项的值
总和更新: 将当前项添加到总和中
π值计算: 根据拉马努金公式计算π的近似值
- π = 1 / (C * total_sum)

字符串转换和收敛检查

// 将π值转换为字符串
mp_exp_t exp;
string pi_str = pi_approx.get_str(exp, 10, 0);// 格式化字符串，确保有小数点
if (exp <= 0) {pi_str = "0." + string(-exp, '0') + pi_str;
} else {pi_str = pi_str.substr(0, exp) + "." + pi_str.substr(exp);
}// 提取小数点后2000位
size_t dot_pos = pi_str.find('.');
string decimal_part;
if (pi_str.length() > dot_pos + 2000) {decimal_part = pi_str.substr(dot_pos + 1, 2000);
} else {decimal_part = pi_str.substr(dot_pos + 1);decimal_part.append(2000 - decimal_part.length(), '0');
}// 检查是否收敛
if (!last_pi_str.empty()) {if (last_pi_str == decimal_part) {convergence_count++;if (convergence_count >= convergence_threshold) {break;}} else {convergence_count = 0;}
}last_pi_str = decimal_part;

字符串转换: 使用get_str方法将高精度浮点数转换为字符串
- exp参数返回指数值，用于确定小数点的位置
- 10表示使用十进制
- 0表示不限制输出长度
格式化字符串: 根据指数值添加小数点
- 如果指数为负，需要在前面补零
- 如果指数为正，在相应位置插入小数点
提取小数部分: 提取小数点后的2000位数字
- 如果字符串长度不足，用零填充
收敛检查: 比较当前迭代与上一次迭代的小数部分
- 如果相同，增加收敛计数器
- 如果连续3次相同，退出循环
- 如果不同，重置收敛计数器

进度输出和安全措施

// 每10次迭代输出一次进度
if (k % 10 == 0) {auto current_time = high_resolution_clock::now();auto elapsed = duration_cast<seconds>(current_time - start_time);cout << "迭代次数: " << k << ", 已用时: " << elapsed.count() << "秒" << endl;
}// 防止无限循环（安全措施）
if (k > 300) {cout << "达到最大迭代次数，退出循环" << endl;break;
}

每10次迭代输出当前进度和已用时间
设置最大迭代次数为300，防止无限循环
拉马努金公式收敛非常快，通常只需要几次迭代

结果输出

// 输出π值
mp_exp_t exp;
string pi_str = pi_approx.get_str(exp, 10, 0);// 格式化字符串，确保有小数点
if (exp <= 0) {pi_str = "0." + string(-exp, '0') + pi_str;
} else {pi_str = pi_str.substr(0, exp) + "." + pi_str.substr(exp);
}size_t dot_pos = pi_str.find('.');
string integer_part = pi_str.substr(0, dot_pos);
string decimal_part = pi_str.substr(dot_pos + 1, 2000);cout << "\nπ值（小数点后2000位）:" << endl;
cout << integer_part << "." << decimal_part << endl;

最终计算完成后，再次格式化π值字符串
分离整数部分和小数部分
输出整数部分和小数点后的2000位数字

5. 完整代码

#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <gmpxx.h>
#include <cmath>using namespace std;
using namespace std::chrono;// 计算阶乘的辅助函数
mpz_class factorial(int n) {mpz_class result = 1;for (int i = 1; i <= n; ++i) {result *= i;}return result;
}// 计算整数幂的辅助函数
mpz_class int_pow(mpz_class base, int exponent) {mpz_class result = 1;for (int i = 0; i < exponent; ++i) {result *= base;}return result;
}int main() {// 设置精度（二进制位数）// 我们需要约2000位十进制精度，所以设置约2000*log2(10) ≈ 6640位二进制精度const int precision = 7000;mpf_set_default_prec(precision);cout << "开始计算π值（使用拉马努金公式）..." << endl;cout << "目标精度：小数点后2000位" << endl;cout << string(50, '=') << endl;auto start_time = high_resolution_clock::now();// 初始化变量mpf_class pi_approx(0, precision);mpf_class total_sum(0, precision);string last_pi_str;int k = 0;int convergence_count = 0;const int convergence_threshold = 3; // 连续3次不变则认为收敛// 常量mpf_class sqrt2(0, precision);mpf_sqrt_ui(sqrt2.get_mpf_t(), 2);mpf_class C = 2 * sqrt2 / 9801;// 迭代计算直到收敛while (true) {// 计算当前项mpz_class numerator = factorial(4 * k) * (1103 + 26390 * k);mpz_class denominator = int_pow(factorial(k), 4) * int_pow(396, 4 * k);mpf_class term = mpf_class(numerator) / mpf_class(denominator);// 添加到总和total_sum += term;// 计算当前π的近似值pi_approx = 1 / (C * total_sum);// 将π值转换为字符串mp_exp_t exp;string pi_str = pi_approx.get_str(exp, 10, 0);// 格式化字符串，确保有小数点if (exp <= 0) {pi_str = "0." + string(-exp, '0') + pi_str;} else {pi_str = pi_str.substr(0, exp) + "." + pi_str.substr(exp);}// 提取小数点后2000位size_t dot_pos = pi_str.find('.');string decimal_part;if (pi_str.length() > dot_pos + 2000) {decimal_part = pi_str.substr(dot_pos + 1, 2000);} else {decimal_part = pi_str.substr(dot_pos + 1);decimal_part.append(2000 - decimal_part.length(), '0');}// 检查是否收敛if (!last_pi_str.empty()) {if (last_pi_str == decimal_part) {convergence_count++;if (convergence_count >= convergence_threshold) {break;}} else {convergence_count = 0;}}last_pi_str = decimal_part;k++;// 每10次迭代输出一次进度if (k % 10 == 0) {auto current_time = high_resolution_clock::now();auto elapsed = duration_cast<seconds>(current_time - start_time);cout << "迭代次数: " << k << ", 已用时: " << elapsed.count() << "秒" << endl;}// 防止无限循环（安全措施）if (k > 300) {cout << "达到最大迭代次数，退出循环" << endl;break;}}auto end_time = high_resolution_clock::now();auto total_time = duration_cast<seconds>(end_time - start_time);cout << string(50, '=') << endl;cout << "计算完成！总迭代次数: " << k << endl;cout << "总用时: " << total_time.count() << "秒" << endl;// 输出π值mp_exp_t exp;string pi_str = pi_approx.get_str(exp, 10, 0);// 格式化字符串，确保有小数点if (exp <= 0) {pi_str = "0." + string(-exp, '0') + pi_str;} else {pi_str = pi_str.substr(0, exp) + "." + pi_str.substr(exp);}size_t dot_pos = pi_str.find('.');string integer_part = pi_str.substr(0, dot_pos);string decimal_part = pi_str.substr(dot_pos + 1, 2000);cout << "\nπ值（小数点后2000位）:" << endl;cout << integer_part << "." << decimal_part << endl;return 0;
}

6. 编译运行

g++ -o pi_calculator pi_calculator.cpp -lgmp -lgmpxx -std=c++11

7. 程序运行结果

开始计算π值（使用拉马努金公式）...
目标精度：小数点后2000位
==================================================
迭代次数: 10, 已用时: 0秒
迭代次数: 20, 已用时: 0秒
迭代次数: 30, 已用时: 0秒
迭代次数: 40, 已用时: 0秒
迭代次数: 50, 已用时: 0秒
迭代次数: 60, 已用时: 0秒
迭代次数: 70, 已用时: 0秒
迭代次数: 80, 已用时: 0秒
迭代次数: 90, 已用时: 0秒
迭代次数: 100, 已用时: 0秒
迭代次数: 110, 已用时: 0秒
迭代次数: 120, 已用时: 0秒
迭代次数: 130, 已用时: 0秒
迭代次数: 140, 已用时: 0秒
迭代次数: 150, 已用时: 0秒
迭代次数: 160, 已用时: 0秒
迭代次数: 170, 已用时: 0秒
迭代次数: 180, 已用时: 0秒
迭代次数: 190, 已用时: 0秒
迭代次数: 200, 已用时: 0秒
迭代次数: 210, 已用时: 0秒
迭代次数: 220, 已用时: 0秒
迭代次数: 230, 已用时: 0秒
迭代次数: 240, 已用时: 0秒
迭代次数: 250, 已用时: 0秒
==================================================
计算完成！总迭代次数: 253
总用时: 0秒总

π值（小数点后2000位）:
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086
51328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344
61284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209
20962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611
79310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217
98609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409
01224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998
37297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959
09216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131
51557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160
35637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521
04752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219
92458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988
18347977535663698074265425278625518184175746728909777727938000816470600161452491921732172147723501414419735685
48161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671
79049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848
86269456042419652850222106611863067442786220391949450471237137869609563643719172874677646575739624138908658326
4599581339047802759009