马青公式计算圆周率程序

马青公式(梅钦公式、Machin)计算圆周率是这样的：

用这公式每增加计算一项，就可以增加约1.39位圆周率准确值，计算百万以下级别精度的圆周率，速度还算满意。

马青公式级数展开后，可以看做两大项，编程可以开两个线程分别同时对这两大项进行计算。为了精简，在这里采用单线程编程。

级数中的分数，分母增长很快，一下子就超出int、long甚至__int64或long long的数值范围，看似要用数组来表示分母才行。其实不必，因为规律性很强，我们对一个分式，不断除以同一个低精度数(25、57121)，就可以得到所有分数的值。

我们采用在程序运行中读入计算位数的形式编程，这样的好处是编译好的程序，可以很方便计算各种不同精度圆周率。

用马青公式编程计算圆周率，相对于网上流传的“四行代码计算圆周率”会复杂些，不过也没复杂多少。看到一些网友用马青公式编程，洋洋洒洒上百行，看都看晕了。

下面是笔者用马青公式计算圆周率最精简的C++代码(万进制)：

#include<stdio.h>
#include<iostream>
using namespace std;
int main(void)
{       //本程序每四位数输出结果，如果请求计算的位数不是4的整数倍，最后输出可能会少1~3位
	long a[2]={956,80},b[2]={57121,25},i=0,j,k,p,q,r,s=2,t,u,v,N,M=10000;
	printf("%9cMachin%6cpi=16arctan(1/5)-4arctan(1/239)\nPlease input a number.\n",32,32);
	cin>>N,N=N/4+3;
	long *pi=new long[N],*e=new long[N];
	while(i<N)pi[i++]=0;
	while(--s+1)
	{
		for(*e=a[k=s],i=N;--i;)e[i]=0;
		for(q=1;j=i-1,i<N;e[i]?0:++i,q+=2,k=!k)
			for(r=v=0;++j<N;pi[j]+=k?u:-u)u=(t=v*M+(e[j]=(p=r*M+e[j])/b[s]))/q,r=p%b[s],v=t%q;
	}
	while(--i)(pi[i]=(t=pi[i]+s)%M)<0?pi[i]+=M,s=t/M-1:s=t/M;
	for(cout<<"3.";++i<N-2;)printf("%04ld",pi[i]);
	delete []pi,delete []e,cin.ignore(),cin.ignore();
	return 0;
}

上面代码计算一万位耗时0.1秒(CPU：Core(i5-2500)3.30 GHz  内存：DDR3 1333MHz)，速度还算过得去。想要更快就得借助运算库、多线程编程或汇编了。