BBP算法计算圆周率（BBP Formula HDU - 6217）

 

概述

BBP算法求圆周率π，它的优点在于可以直接求出圆周率的某一位d开始的一串数字，而不依赖于第d位之前的数字，传统的算法需要一位位地求。BBP算法节省了一大部分时间和内存，但BBP所求的圆周率是用16进制表示的。

公式推导

核心公式：

 

第一步：将上式分解成四段求解，可以写成通项公式：

整理一下：

求第n+1位时，转换一下公式：

两边乘上：

第一个求和部分（n-k）>0，所以将分子取分母的模，目的是为了让分子分母相除后得到小数部分，
第二个求和部分（n-k）<0，无需再次取模：

小数部分（如果ans<0，则ans再加上1，使他变为正数）：

最后ans*16转化为16进制，取整数部分。

 

 题目链接：HDU-6217（BBP Formula）

#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;
}

 

参考文献：

•  [圆周率pi的BBP计算公式之详细证明]