概述
BBP算法求圆周率π,它的优点在于可以直接求出圆周率的某一位d开始的一串数字,而不依赖于第d位之前的数字,传统的算法需要一位位地求。BBP算法节省了一大部分时间和内存,但BBP所求的圆周率是用16进制表示的。
公式推导
核心公式:
第一步:将上式分解成四段求解,可以写成通项公式:
整理一下:
求第n+1位时,转换一下公式:
两边乘上:
第一个求和部分(n-k)>0,所以将分子取分母的模,目的是为了让分子分母相除后得到小数部分,
第二个求和部分(n-k)<0,无需再次取模:小数部分(如果ans<0,则ans再加上1,使他变为正数):
最后ans*16转化为16进制,取整数部分。
#include <set> #include <map> #include <cmath> #include <stack> #include <queue> #include <vector> #include <string> #include <cstdio> #include <cstring> #include <sstream> #include <iomanip> #include <iostream> #include <algorithm> #define clr(str,x) memset(str,x,sizeof(str)) #define FRER() freopen("in.txt","r",stdin); #define FREW() freopen("out.txt","w",stdout); #define MAX_INF 0x7fffffff #define INF 0x3f3f3f3f #define maxn typedef long long int ll; using namespace std; ll q_pow(ll a,ll b,ll mod) { ll res=1; while(b) { if(b&1) res=res*a%mod; a=a*a%mod; b/=2; } return res; } double bbp(int n,ll k,ll b) { double res=0; for(int i=0;i<=n;i++) res+=(q_pow(16,n-i,8*i+b)*1.0/(8*i+b)); for(int i=n+1;i<=n+1000+1;i++) res+=pow(16,n-i)*1.0/(8*i+b); return k*res; } int main() { //FRER() //FREW() int T,Case=1,n; cin>>T; while(T--) { double ans=0; cin>>n; ans=bbp(n-1,4,1) - bbp(n-1,2,4) - bbp(n-1,1,5) - bbp(n-1,1,6); ans = ans - (int)ans; if(ans < 0) ans+=1.0; ans*=16.0; int x=(int)ans; char c; if(x>=0&&x<=9) c=x+'0'; else c=x-10+'A'; printf("Case #%d: %d %c\n", Case++, n, c); } return 0; }
参考文献: