【Python】计算任意位数的圆周率π（Machin Formula）

本文链接：https://blog.youkuaiyun.com/m0_65005346/article/details/125939287

本文介绍了使用马青公式计算π的原理，并提供了Python和C++两种编程语言的实现方式，详细阐述了计算过程中需要注意的精度问题及解决方案。通过代码示例，展示了如何通过循环和取模运算来逼近π的值，同时讨论了计算精度与代码中固定参数的矛盾，并提出了可能的改进方案。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

前言

今天的python课，我给学生介绍math库。其中一个孩子问，能看到PI的值吗？我马上开始了我的表演。进入命令行，导入，调用math.pi,一气呵成。谁知，这孩子并不买账！嚷嚷道：我要看小数点后100位@@@

一、马青公式

Machin公式：π=16arctan(1/5)-4arctan(1/239)

二、Python代码实现

1.理清思路

计算π 的公式已经拿到，我们有一个明确的想要精确的π的位数，那么当式子中n越大的时候，第n项就会越小，小到计算这个精确目标时可以忽略的地步。

注意1：由于程序中浮点数的有效位数有限（小数点后很长的话会被省略，很大的数如果是小数也会用科学技术法），会损失很大精度，所以用整数来计算π。

比如想计算小数点后 7 位，就计算 π * 10^7 是多少，返回 31415926。

比如想计算小数点后 3 位，程序就会给你返回整数 3141，是扩大了 10^3 的结果。

注意2: 想用整数计算 π，那除法的时候就要使用取模运算来保证结果也是整数，以达到保证精度的目的。

注意3: 用取模操作代替除法，但是取模对于除法来说本身就是有精度损失的，为了弥补这部分精度，可以先把想精确的位数增加 k，最后把结果减少 k 位。

注意4: 注意n:求 π 的数学公式中的 n，还是精确到π 的位数。

代码如下（python）：

# 计算准确圆周率的马青公式
import time
n = int(input('请输入圆周率小数点后的位数:'))
start_time = time.time()   # 计算pi计时起点
w = n + 10                 # 考虑精度，位数增加 k，结果减少 k 位
b = 10 ** w                # 用整数计算pi  如3.14 * 10^2 = 314
x1 = b * 4// 5             # 第一项前半部分
x2 = b // -239             # 第一项后半部分
sum = x1 + x2              # 第一项的值
n *= 2
for i in range(3, n, 2):
    x1 //= -25
    x2 //= -239*239
    x = (x1+x2)//i
    sum += x
mpi = sum * 4
mpi //= 10**10
end_time = time.time()
run_time = str(end_time - start_time)
mpi_str = str(mpi)
mpi_str = mpi_str[:1] + '.' + mpi_str[1:]
# f = open('mpi.txt', 'w')
# f.write(mpi_str)
# f.close()
print('运行时间: ' + run_time)
print('计算结果: ',mpi_str)
print('')

2.C++版本

代码如下（C++）：

/**********************************************************
 *使用梅钦公式π=16arctan(1/5)-4arctan(1/239)
 *通过建立数组来储存超过类型范围的小数部分
 *用long类型数组pi[ ]来储存各项（每一项有4位)
 *使用long类型term[ ]来实现对每一次泰勒展开结果的储存
 *之后通过进位等操作，实现用数组储存足够位数的π。         
 ***********************************************************/
void machin_pi(long L)
{
    //常量C用于进位退位操作，10000代表pi数组以四位进行储存
    const long C = 10000;

    //flag用于正负号的判定
    short flag;

    //梅钦公式具体应用时π=80*（1/25-1/（25*25*3）+......)+(-956)*(......)
    //xishu[]为两项展开的外系数，div[]为每次进阶指数所用的除数
    //odd为1,3，5等奇数，shang_1，r_1是进阶指数的商和余数,r_2是计算π各项的余数
    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;

    L = L / 4 + 3;

    //term[]用来完成对每个升幂项的足够长度保存
    //pi[]用来储存最后结果
    long *term = new long[L];
    long *pi = new long[L];

    while (i < L)pi[i++] = 0;

    while (k - 2)
    {
        //完成对term的初始化
        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];                //t、r_1、term[]用来储存25与57121的各幂
                term[j] = t / div[k];                 //每次进阶都用上一轮已经算好的term数组进行升阶（后一项都是前一项除以25或57121）
                r_1 = t % div[k];                     //保存余数，在下一项中×C用于退位
                
                tt = r_2 * C + term[j];               //tt、r_2、shang_1用来储存泰勒展开各项
                shang_1 = tt / odd;                   //在term[]的基础上除以奇数1、3、5等
                r_2 = tt % odd;                       //保存余数，在下一项中用于退位

                pi[j] += (flag?shang_1:-shang_1);     //用flag进行标定正负号
            }
            if (term[i])--i;                          //当term[i]==0时，后续展开不会影响这一项（这四位小数），进行++i,直接从下一项开始
            odd += 2; 
            flag=!flag;                               //泰勒展开正负号改变                               
        }
        ++k;                                          //处理完arccot5后进入下一轮处理arccot239
    }
    k = 0;

    //上述循环之后，pi[]中已经储存了各项数据
    //但pi[]中数据可能存在负数、五位数及以上等等情况
    //接下来的循环进行处理，通过进位全部转化成正四位数
    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;
        }
    }

    //输出圆周率
    printf("3.");
    for (i = 1; i < len / 4 + 1; ++i)
    {
        printf("%04ld", pi[i]);
    }
    switch (len%4)
    {
    case 0:break;
    case 1:printf("%01ld", pi[i]/ 1000); break;
    case 2:printf("%02ld", pi[i]/ 100); break;
    case 3:printf("%03ld", pi[i]/ 10); break;
    }

    //释放内存
    delete[] term;
    delete[] pi;
}