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

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

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;
}
  • 1
    点赞
  • 18
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
好的,我知道你想让我帮你解决关于函数模块和循环的问题。关于使用 Machin 公式计算圆周率,我可以给你一些指导。 首先,让我们了解一下 Machin 公式Machin 公式是一种用于计算圆周率公式,它的公式如下: π / 4 = 4 × arctan(1/5) - arctan(1/239) 其中,arctan 表示反正切函数。 实际上,这个公式是由丹尼尔·贝尔努利(Daniel Bernoulli)在 1655 年发现的,但由于其后被约翰·马钦(John Machin)在 1706 年重新发现,所以这个公式就被称为 Machin 公式。 接下来,我们可以将这个公式转化为计算机程序的形式。在程序中,我们需要使用反正切函数和循环来计算圆周率。 下面是一个使用 Python 语言实现 Machin 公式计算圆周率的示例代码: ```python import math def machin_formula(n): pi_sum = 0 for i in range(n): pi_sum += ((-1) ** i) / (2 * i + 1) * (4 / 5) ** (2 * i + 1) pi_sum -= ((-1) ** i) / (2 * i + 1) * (1 / 239) ** (2 * i + 1) pi = 4 * pi_sum return pi n = 1000 pi = machin_formula(n) print("圆周率的值为:", pi) print("与 math 模块中的值相差:", abs(math.pi - pi)) ``` 在这个程序中,我们定义了一个名为 `machin_formula` 的函数,它接受一个参数 `n`,表示计算的精度。函数中使用了一个循环,从 0 到 `n-1`,根据 Machin 公式计算圆周率的值。 你可以尝试修改 `n` 的值,来得到不同精度的圆周率值。当然,值越大,计算时间就越长。 最后,我们使用 `print` 函数来输出计算出的圆周率值,并与 Python 内置的 `math` 模块中的值进行比较,以验证我们的程序的正确性。 希望这个示例能够对你有所帮助,如果你有任何问题,请随时向我提问。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值