程序计算精确圆周率π的方法

wabil

已于 2024-07-16 22:37:19 修改

阅读量1.4w

点赞数 3

分类专栏： VC 文章标签：算法

于 2016-10-27 16:24:21 首次发布

本文链接：https://blog.csdn.net/wabil/article/details/52946977

版权

VC 专栏收录该内容

10 篇文章 0 订阅

订阅专栏

圆周率π=3.14159265358……无穷多位[事实上是超越数], 自古以来，人们就对π的精确计算充满兴趣。从最早的几何方法到现代的计算机算法，计算π的精度逐步提高。本文将介绍快速计算圆周率的算法以及它们在计算机上的应用。

首先给一种王炸级别的获取方法：python+mpmath库【pip install mpmath】

from mpmath import mp
mp.dps = 1000 
print( mp.pi)   ## 输出1000位圆周率

Python 运行效果
这里提前为大家轻松列举1000位圆周率，其实mpmath内部的算法也是用c语言实现的。

下面我们选一种关于π比较简单, 并且收敛非常快的级数。
π /2= 1 + 1/3 + (1/3)(2/5) + (1/3)(2/5)(3/7) +(1/3)(2/5)(3/7)(4/9) + …
调整下：
π = 2 + 2/3 + (2/3)(2/5)+(2/3)(2/5)(3/7)+(2/3)(2/5)(3/7)(4/9) + …

void main()
{
	double x = 2, z = 2;
	int a = 1, b = 3;
	while (z > 0)
	{   //when b->max, z->0
		z = z * a / b;
		x += z;
		a++;
		b += 2;
	}
	printf("π=%.14f\n", x);
}

执行效果： π=3.1415926535898 ；这个程序如此简单它可以计算到小数点后14位，因受限于double的精度。【感兴趣的可以使用GMP大数库，MPFR高精度库再做尝试】

把上面的程序改造下, 将double类型的x，z变成数组保存每位小数，让它精确到小数点后面 1000 位再测试试试！！


void main()
{
	const int ARRSIZE = 1010;
	const int DISPCNT = 1000;	 // 定义数组大小，显示位数
	char x[ARRSIZE], z[ARRSIZE]; // x[0] x[1] . x[2] x[3] x[4] .... x[ARRSIZE-1]
	int a = 1, b = 3, c, d, i;
	int Run = 1, Cnt = 0;
	memset(x, 0, ARRSIZE);
	memset(z, 0, ARRSIZE);
	x[1] = 2;
	z[1] = 2;

	while (Run && (++Cnt < 100000))
	{
		d = 0;
		i = ARRSIZE - 1;
		do
		{
			c = z[i] * a + d;
			z[i] = c % 10;
			d = c / 10;
		} while (i--);

		d = 0;
		i = 0;
		do
		{
			c = z[i] + d * 10;
			z[i] = c / b;
			d = c % b;
		} while (i++ < ARRSIZE);

		Run = 0;
		for (i = ARRSIZE - 1; i > 0; i--)
		{
			c = x[i] + z[i];
			x[i] = c % 10;
			x[i - 1] += c / 10;
			Run |= z[i];
		}
		a++;
		b += 2;
	}
	printf("calc %d times\r\n", Cnt);
	printf("Pai=%d.\r\n", x[1]);
	for (i = 0; i < DISPCNT; i++)
	{
		if (i && ((i % 100) == 0))
			printf("\r\n");
		printf("%d", (int)x[i + 2]);
	}
}

以下是运行的效果，可以看到经过3343轮计算后前面的999位数字都是对的！
在这里插入图片描述
这下心里就有底了, 是不是改变数组大小就可以计算更多位数呢？答案是肯定的。
如果把定义数组大小和显示位数改为:
const int ARRSIZE = 10100;
const int DISPCNT = 10000;
执行的结果表明计算33538轮后圆周率精度就可以达到惊人的10000位！
初步结论，这个算法可以通过计算大概3.3轮就能输出一位正确的数值。

趣味扩展
基于以上算法，其实我们就可以对π的一些属性展开研究，如统计π每个位数数字出现的概率，等等

算法提高精度的原理:
以上程序的原理是利用数组把计算结果保存起来, 其中数组每一项保存10进制数的一位,
小数点定位在数组第1个数和第2个数之间, 即小数点前面2位整数, 其余都是小数位。
利用电脑模拟四则运算的笔算方法来实现高精度的数据计算,没想到最原始的方法竟然是精度最高的。

利用这个原理我们就可以扩展到其他常量的计算如自然底数e=2.71828…任意位数

#include <stdio.h>
#include <string.h>

// e = 1+1/1!+1/2!+1/3!+…+1/n! = 1+1/2*(1+1/3*(1+…*(1+1/n))…)));
int main()
{
#define TIMES 100 /// 输出e有效位数
	int n, e, d, i, f[TIMES + 1];
	i = TIMES;
	while (i--)
		f[i] = 10;
	printf("e=2.");
	i = 0;
	e = 0;
	while (i++ < TIMES)
	{
		d = 0;
		for (n = TIMES; n > 0; n--)
		{
			d = d + f[n] * 10;
			f[n] = d % n; // 余数
			d = d / (n);  // 整数
		}
		printf("%d", (e + d / 10) % 10);
		e = d % 10;
	}
	sleep(1);
	printf("\n\n");
	return 0;
}

用python+mpmath打印下答案对比

import mpmath
mpmath.mp.dps = 1000
# 计算自然对数的底数e的值，1000位
print(mpmath.e)

作为文章结束，给出一个某年Obfuscated C Contest佳作选录: 800位精度圆周率

#include <stdio.h>
long a=10000, b, c=2800, d, e, f[2801], g; 
main(){
for(;b-c;)f[b++]=a/5; 
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a) 
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}