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

概述

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

公式推导

核心公式：

$\pi ={\sum }_{k=0}^{\infty }\left\{\frac{1}{16^k}(\frac{4}{8k-1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6})\right\}$

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

$S_i={\sum }_{k=0}^{\infty }\left\{\frac{1}{16^k}\ast \frac{d}{8k+j}\right\}$

整理一下：
$S_i=d\ast {\sum }_{k=0}^{\infty }\left\{\frac{1}{16^k\ast (8k+j)}\right\}$
求第n+1位时，转换一下公式：

$S_i=d\ast {\sum }_{k=0}^n\left\{\frac{1}{16^k\ast (8k+j)}\right\}+d\ast {\sum }_{k=n+1}^{\infty }\left\{\frac{1}{16^k\ast (8k+j)}\right\}$

两边乘上$16^n$：

$16^n\ast S_i=d\ast {\sum }_{k=0}^n\left\{\frac{16^{n-k}}{(8k+j)}\right\}+d\ast {\sum }_{k=n+1}^{\infty }\left\{\frac{16^{n-k}}{(8k+j)}\right\}$

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

$t_i=d\ast {\sum }_{k=0}^n\left\{\frac{16^{n-k}mod(8k+j)}{(8k+j)}\right\}+d\ast {\sum }_{k=n+1}^{\infty }\left\{\frac{16^{n-k}}{(8k+j)}\right\}$

小数部分（如果ans<0，则ans再加上1，使他变为正数）：
$ans=\sum t_i-(int)\sum t_i$

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

试题链接 ：http://acm.hdu.edu.cn/showproblem.php?pid=6217

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

6
5
Case #1: 5 6
11
Case #2: 11 A
11
Case #3: 11 A
111
Case #4: 111 D
25
Case #5: 25 0
6
Case #6: 6 A

Process returned 0 (0x0)   execution time : 3031.679 s
Press any key to continue.