Machin(梅钦/马青)公式计算圆周率π

编程基础 • 2024-12-20 23:30 • 阅读 137

作为一名计算机的初学者，因为老师作业的要求去完成圆周率的计算，因此突然产生了兴趣，想尝试自己用梅钦公式来完成这个任务。上网找了一些资料，也看见了短短几行就完成任务的代码，实在佩服，不过那样的代码实在对于我这个初学者来说太艰涩了，于是我自己大致了解了之后，自己完成了一遍，作为自己的第一次博客，来记录自己的思绪。

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

总体想法来说，是通过建立数组来储存超过类型范围的小数部分，这里我用long类型数组pi[ ]来储存各项（每一项有4位），使用long类型term[ ]来实现对每一次升幂结果的储存，之后通过进位等操作，实现用数组储存足够位数的π。

/ *使用梅钦公式π=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; }

今天的文章 Machin(梅钦/马青)公式计算圆周率π分享到此就结束了，感谢您的阅读。

Machin(梅钦/马青)公式计算圆周率π

相关推荐