UVALive8201-BBP公式计算圆周率

在1995年，Simon Plouffe 发现了一种特殊的求和方法来表示某些常数。两年后，David H. Bailey 和 Peter Borwein 在他们的一篇论文中详细介绍了这一发现，并将其命名为 Bailey-Borwein-Plouffe (BBP) 公式。这一公式的出现令人震惊，因为它提供了一种计算 π 的第 n 个十六进制数字的方法，而无需计算前 n-1 个数字。

BBP 公式如下：

π = Σ (k=0 到 ∞) 16^(-k) * (4/(8k + 1) - 2/(8k + 4) - 1/(8k + 5) - 1/(8k + 6))

这个问题要求你计算 π 的第 n 个十六进制数字，紧随小数点之后。例如，π 的十六进制表示为 3.243F6A8885A308D313198A2E...，其中第 1 个数字是 2，第 11 个数字是 A，第 15 个数字是 D。

输入

输入的第一行包含一个整数 T (1 ≤ T ≤ 32)，表示测试用例的数量。接下来的每一行包含一个整数 n (1 ≤ n ≤ 100000)，表示需要计算的十六进制数字的位置。

输出

对于每个测试用例，输出一行，格式为“Case #X: N D”，其中 X 是测试用例的编号，N 是给定的 n 值，D 是计算得到的十六进制数字（0-9 或 A-F）。

样例输入

5
1
11
111
1111
11111

样例输出

Case #1: 1 2
Case #2: 11 A
Case #3: 111 D
Case #4: 1111 A
Case #5: 11111 E

解决方案

为了计算一个十进制小数 x 的第 n 位十六进制小数，可以先将 x 转换为十六进制，然后将小数点右移 n 位。这样，小数点左侧的第一个整数位就是所需的答案。然而，由于 n 可能非常大，直接计算会导致精度不足。因此，可以将 x 乘以 16^n，即将小数点右移 n 位。

对于 BBP 公式，可以将其分解为四个部分：

π = 4Σ1 - 2Σ2 - Σ3 - Σ4

其中：

Σ1 = Σ (k=0 到 ∞) 16^(-k) / (8k + 1)

Σ2 = Σ (k=0 到 ∞) 16^(-k) / (8k + 4)

Σ3 = Σ (k=0 到 ∞) 16^(-k) / (8k + 5)

Σ4 = Σ (k=0 到 ∞) 16^(-k) / (8k + 6)

以 Σ1 为例，为了计算第 n 位小数，可以将小数点右移 n-1 位：

Σ1 = Σ (k=0 到 n-1) 16^(n-1-k) / (8k + 1) + Σ (k=n 到 ∞) 16^(n-1-k) / (8k + 1)

由于整数部分对最终结果没有影响，可以通过取模去除整数部分，保留小数部分：

Σ1 = Σ (k=0 到 n-1) [16^(n-1-k) mod (8k + 1)] / (8k + 1) + Σ (k=n 到 ∞) 16^(n-1-k) / (8k + 1)

前半部分可以通过快速幂算法轻松实现，而后半部分是一个收敛级数，随着 k 的增加，对应的项逐渐趋近于零。因此，可以将无穷大取为一个较大的数，如 500。

以下是实现代码：

#include 
#include 
#include 
#include 
#include 
#define re register
#define il inline
#define ll long long
defined ld long double
using namespace std;
const ll MAXN = 1e2 + 5;
const ll TABLE = 26;
const ld INF = 1e9;
const ld EPS = 1e-9;

// 快速幂模
ll powmod(ll a, ll n, ll md) {
    ll ans = 1;
    while (n) {
        if (n & 1) {
            ans = (ans * a) % md;
        }
        a = (a * a) % md;
        n >>= 1;
    }
    return ans;
}

// BBP Formula
ll BBP(ll n) {
    ld ans = 0;
    for (re ll i = 0; i         ll k = n - 1 - i;
        ll a = 8 * i + 1;
        ll b = a + 3;
        ll c = a + 4;
        ll d = a + 5;
        ans += 4 * powmod(16, k, a) / (ld)a - 2 * powmod(16, k, b) / (ld)b - powmod(16, k, c) / (ld)c - powmod(16, k, d) / (ld)d;
    }
    for (re ll i = n; i <= n + 10; ++i) {
        ll k = n - 1 - i;
        ll a = 8 * i + 1;
        ll b = a + 3;
        ll c = a + 4;
        ll d = a + 5;
        ans += 4.0 * powl(16.0, (ld)k) / a - 2.0 * powl(16.0, (ld)k) / b - 1.0 * powl(16.0, (ld)k) / c - 1.0 * powl(16.0, (ld)k) / d;
    }
    ans -= (ll)ans;
    ans = ans <0 ? ans + 1 : ans;
    return ((ll)(ans * 16)) % 16;
}

int main() {
    ios::sync_with_stdio(false);
    int T;
    cin >> T;
    for (re int i = 1; i <= T; ++i) {
        int n;
        cin >> n;
        ll ans = BBP(n);
        cout <<"Case #" <        cout <    }
    return 0;
}