题解：P1727 计算π

· · 题解

以下部分内容来自 https://www.numberworld.org，有兴趣的读者可以去了解更多内容。

截止至 2025 年 7 月 15 日，本做法为非打表做法最优解，总用时仅 98ms，记录，甚至 10^6 位的计算也可以在 2s 内完成，如果有人有更优的做法也请告诉我。

我们都知道， Chudnovsky 公式是一个非常快计算 \pi 的公式，每计算一项可以增加约 14.18 位十进制精度。公式如下：

\dfrac 1\pi=\dfrac{\sqrt{10005}}{4270934400}\Large \sum\limits^\infty_{k=0}\normalsize(-1)^k\dfrac{(6k)!}{(k!)^3(3k)!}\dfrac{13591409+545140134k}{640320^{3k}}

这个公式很复杂，我们显然不能一项一项进行计算，这会使时间复杂度达到 O(n^2\log n)，再加上常数的影响，无法体现这个公式的优越性。

假设我们要计算以下形式的式子：

\large\frac{P(0,n)}{Q(0,n)}=\sum\limits^n_{k=1}\frac{P(k)}{R(k)}\prod\limits^k_{i=1}\frac{R(i)}{Q(i)} $P(l,r)$ 为计算 $l$ 到 $r$ 项（$k$ 从 $l$ 到 $r$）时 $P$ 的值，$Q,R$ 同理。 取 $m=\lfloor\frac{l+r}2\rfloor P(l,r)=P(l,m)Q(m,r)+P(m,r)R(l,m) Q(l,r)=Q(l,m)Q(m,r) R(l,r)=R(l,m)R(m,r)

注意上面的是 m 而不是 m+1，也就是区间中点处会在两侧同时出现。

边界条件（由于端点处的特殊性，若在初始时有 l<r，则不会出现 l=r 的情况）：

P(r-1,r)=P(r) Q(r-1,r)=Q(r) R(r-1,r)=R(r)

这种算法叫做 binary splitting。

对于所有 k\ge1，需要满足 R(k)<Q(k)，级数才可以较快的收敛。

我们把 Chudnovsky 公式改写成上面式子的形式（把阶乘和指数运算写到乘积里面），得：

P(r)=(-1)^r(13591409+545140134r)(2r-1)(6r-5)(6r-1) Q(r)=10939058860032000r^3 R(r)=(2r-1)(6r-5)(6r-1) \pi=\dfrac{1}{\sqrt{10005}}\dfrac{4270934400Q(0,n)}{P(0,n)+13591409Q(0,n)}

我们采用上面的方法分治计算，由于 Q 的三次项系数比 R 的三次项系数大很多，所以会以很快的速度收敛，这可以得出这个级数的误差约为 (\frac{10939058860032000}{2\times6\times6})^{-n}=151931373056000^{-n}，\lg 151931373056000\approx14.18，这就是文章开头所说的每计算一项增加的十进制精度的由来。

接下来我们分析这个算法的时间复杂度，我们要精确到 N 位数字，就要计算 O(N) 项级数。

如果我们把 P 和 Q 用高精度整数计算出来，那么时间复杂度 T(n)=2T(n/2)+O(n\log n)，即 T(n)=O(n\log^2n)，这里的位数 n=O(N\log N)（n! 有 O(n\log n) 位），所以复杂度为 O(N\log^3N)。

但是如果我们在计算中使用高精度浮点数，只保留 N 位精度，多余的精度舍去，这会使计算快上不少，再加上别的地方常数的影响，这个算法甚至比理论复杂度为 O(n\log^2n) 的 AGM 算法要快上不少。

参考实现：

void _pi_impl(int l,int r,BigFloat &P,BigFloat &Q,BigFloat &R){
    if (l+1==r){
        P=(BigFloat(13591409)+BigFloat(545140134)*r)*(R=BigFloat(2*r-1)*(6*r-5)*(6*r-1));
        if (r&1) P.s.sign^=1;
        Q=BigFloat(10939058860032000ll)*r*r*r;
        return;
    }
    int mid=(l+r)>>1;
    _pi_impl(l,mid,P,Q,R);
    BigFloat P2,Q2,R2;
    _pi_impl(mid,r,P2,Q2,R2);
    P*=Q2;
    P+=P2*R;
    Q*=Q2,R*=R2;
}
// Chudnovsky
BigFloat pi(){ // eps=O(151931373056000^-n)
    BigFloat p,q,r;
    // auto st=clock();
    _pi_impl(0,ceil(BigFloat_prec*digit/std::log10(151931373056000)),p,q,r);
    // std::cout<<"main:"<<double(clock()-st)/CLOCKS_PER_SEC*1000<<"ms\n";
    // st=clock();
    auto tmp=sqrt(BigFloat(10005));
    // std::cout<<"sqrt:"<<double(clock()-st)/CLOCKS_PER_SEC*1000<<"ms\n";
    return BigFloat(4270934400ll)*q/(tmp*(p+q*13591409));
}

高精度模板可以去这里自取，同时也欢迎大家来查错。