高精度黑科技——二进制分裂（Binary Splitting）

gitxiaozheng · 2025-08-10 17:37:43 · 算法·理论

什么是二进制分裂？

二进制分裂（Binary Splitting）可以十分快速地解决求一类无理数小数点后 n 位问题，这一类无理数的计算公式通常可以化成如下形式：

\sum^{\infin}_{i=1} \frac{P(i)}{R(i)} \prod^{i}_{j=1}\frac{R(j)}{Q(j)}

与其他在这方面的算法相比，该算法有以下优势： - **适用范围广**。我们可以用它轻松解决求 $\pi$，$e$ 和一系列著名的无理数小数点后 $n$ 位的问题。 - **更好理解**。理解这个做法，你只需要了解 $\sum$ 和 $\prod$ 的基础知识，拥有分式通分的能力以及拆式子的毅力。 - **相对好写**。（注意只是相对于其他这方面的算法，肯定没有暴力好写） - 为什么不用好写的暴力呢？自然是因为它**效率很高**。二进制分裂效率同样优于很多非暴力算法，是解决这一类问题最快的算法之一。二进制分裂是 2024 年 5 月提出来的，在 [numberworld](https://www.numberworld.org/y-cruncher/internals/binary-splitting.html) 上面有原文，原文有二进制分裂的三种形式，这里主要介绍一般形式和其中一个特殊形式，同时原文作者也附有[更多用二进制分裂计算著名无理数的实例](https://www.numberworld.org/y-cruncher/internals/binary-splitting-library.html)。 ## 具体流程我们直接讨论一般性问题的解法：计算无理数 $a$ 后 $m$ 位的值，公式如下所示： $$ a=\sum^{\infin}_{i=1} \frac{P(i)}{R(i)} \prod^{i}_{j=1}\frac{R(j)}{Q(j)} $$ ### 前置定义我们先定义部分和 $S(l,r)$ 为： $$ S(l,r)=\sum^{r}_{i=l} \frac{P(i)}{R(i)} \prod^{i}_{j=l}\frac{R(j)}{Q(j)} $$ 考虑将 $S(l,r)$ 写成类似于 $\frac{A}{B}$ 的形式。因此我们可以定义 $$ P(l,r)=\sum^{r}_{i=l} P(i) \cdot (\prod^{i-1}_{j=l}R(j)) \cdot (\prod^{r}_{j=i+1}Q(j)) $$ $$ Q(l,r)=\prod^{r}_{i=l}Q(i) $$ $$ R(l,r)=\prod^{r}_{i=l}R(i) $$ 因此我们便将 $P(l,r)$ 视作 $S(l,r)$ 的分子，$Q(l,r)$ 视作部分和的分母，$R(l,r)$ 为辅助计算的函数，即现在 $S(l,r)$ 可以表示为： $$ S(l,r)=\frac{P(l,r)}{Q(l,r)} $$ ::::info[如果你不知道为什么 $A(x)$ 为 $S(l,r)$ 的分子] 我们把 $S(l,r)$ 直接化成类似于 $\frac{A}{B}$： $$ \begin{aligned} S(l,r)&= \sum^{r}_{i=l} \frac{P(i)}{R(i)} \cdot \frac{\prod^{i}_{j=l}R(j)}{\prod^{i}_{j=l}Q(j)}\\ &=\sum^{r}_{i=l} \frac{P(i) \cdot \prod^{i-1}_{j=l}R(j)}{\prod^{i}_{j=l}Q(j)} \\ &=\frac{\sum^{r}_{i=l} P(i) \cdot (\prod^{i-1}_{j=l}R(j)) \cdot (\prod^{r}_{j=i+1}Q(j))}{\prod^{r}_{i=l}Q(i)}\\ \end{aligned} $$ 这样子一一对应上去，便能得出 $P(x)$ 的形式和含义。 :::: 假设我们要快速求 $S(1,r)$，那么我们现在就要考虑如何快速算出 $P(1 ,r),Q(1,r),R(1,r)$。 ### 递归过程这里开始就是二进制分裂的核心部分，其思想是把整个求和区间掰成两半，分成 $[l,m]$ 和 $[m+1,r]$ 两部分递归处理，并将各部分乘积结果合并。以下取 $m=\lfloor \frac{l+r}{2}\rfloor$。那么每次递归，我们采取如下的合并式子计算即可： $$ \begin{cases} P(l,r)=P(l,m) \cdot Q(m+1,r)+R(l,m) \cdot P(m+1,r)\\ Q(l,r)=Q(l,m) \cdot Q(m+1,r)\\ R(l,r)=R(l,m) \cdot R(m+1,r) \end{cases} $$ ::::info[具体推导过程]{open} - 对于 $Q(l,r)$ 和 $R(l,r)$ 的合并，将掰开的两部分结果相乘即可。这没什么好说的。 $$ \begin{aligned} Q(l,r)&=\prod^{m}_{i=l}Q(i) \cdot \prod^{r}_{i=m+1}Q(i) \\ &=Q(l,m) \cdot Q(m+1,r) \end{aligned} $$ $$ \begin{aligned} R(l,r)&=\prod^{m}_{i=l}R(i) \cdot \prod^{r}_{i=m+1}R(i) \\ &=R(l,m) \cdot R(m+1,r) \end{aligned} $$ - 对于 $P(l,r)$，这一堆东西看起来不好直接得出合并的式子，考虑将 $S(l,r)$ 作如下拆分（注意细微差别）： $$ \begin{aligned} S(l,r)&= \sum^{\color{red} m}_{i=l} \frac{P(i)}{R(i)} \prod^{i}_{j=l}\frac{R(j)}{Q(j)}+(\prod^{m}_{j=l}\frac{R(j)}{Q(j)}) \cdot (\sum^{\color{red} r}_{i=m+1} \frac{P(i)}{R(i)} \prod^{i}_{\color{red} j=m+1}\frac{R(j)}{Q(j)})\\ &=S(l,m)+(\prod^{m}_{j=l}\frac{R(j)}{Q(j)}) \cdot S(m+1,r)\\ &=S(l,m)+\frac{R(l,m)}{Q(l,m)} \cdot S(m+1,r) \\ &=\frac{P(l,m)}{Q(l,m)}+\frac{R(l,m)}{Q(l,m)} \cdot \frac{P(m+1,r)}{Q(m+1,r)} \\ &=\frac{P(l,m) \cdot Q(m+1,r)+R(l,m) \cdot P(m+1,r)}{Q(l,m)\cdot Q(m+1,r)}\\ &=\frac{P(l,m) \cdot Q(m+1,r)+R(l,m) \cdot P(m+1,r)}{Q(l,r)} \\ \end{aligned} $$ 由于我们知道 $S(l,r)=\frac{P(l,r)}{Q(l,r)}$，于是自然而然的得到了 $P(l,r)$ 的合并式子，即 $P(l,r)=P(l,m) \cdot Q(m+1,r)+R(l,m) \cdot P(m+1,r)$。 :::: ### 递归边界很明显当 $l=r$ 时结束递归。此时： $$ \begin{cases} P(l,r)=P(l) \\ Q(l,r)=Q(l) \\ R(l,r)=R(l) \end{cases} $$ 至此，基本流程就这么结束了，我们快速算出了 $P(1,r),Q(1,r),R(1,r)$ 的值，直接

S(l,r)=\frac{P(l,r)}{Q(l,r)}

## 代码框架 ~~人生苦短，我用 Python！~~ 这是一个用 Python 实现的代码框架。 ```python def P(x): # return val def Q(x): # return val def R(x): # return val # P(x),Q(x),R(x)为关于x的易于计算的函数。 def bs(l,r): if l==r: return (P(l),Q(l),R(l)) m=(a+b)//2 P1,Q1,R1=bs(l,m) P2,Q2,R2=bs(m+1,r) P=P1*Q2+P2*R1 Q=Q1*Q2 R=R1*R2 return (P,Q,R) # P1=P(l,m) # Q1=Q(l,m) # R1=R(l,m) # P2=P(m+1,r) # Q2=Q(m+1,r) # R2=R(m+1,r) # P,Q,R=(P(l,r),Q(l,r),R(l,r)) ``` ## 时间复杂度假设我们要计算 $S(1,n)$。那么，这个算法如果使用 FFT 实现高精度乘法，一般复杂度为 $O(n \log^3 n)$。 ::::info[分析与证明]{open} 递归次数显然是为 $O(\log n)$，递归到第 $k$ 层要合并 $2^k$ 个区间，计算三次乘法。合并一次区间时间复杂度为 $O(M(d))$，所以第 $k$ 层合并所有区间时间复杂度为 $O(2^k \cdot M(d))$，$d$ 为做高精度乘法时两个大数的位数。具体分析 $O(M(d))$： - 使用 FFT 加速乘法，则单次高精度乘法时间开销 $M(d)=O(d \log d)$。 - 由于 $P(x),Q(x),R(x)$ 为有理函数，$x\leq n$ 时其值位数为 $O(\log n)$，区间长度 $L=\frac{n}{2^k}$，因此得大数位数 $d=O(L \log n)$。最后得到合并单区间时间复杂度为： $$O(M(L \log n))=O(L \log n \cdot \log(L \log n))=O(L \log^2 n)=O(\frac{n}{2^k} \log^2 n)$$ （你要知道 $L \leq n$ 的时候，$\log(L \log n) \le \log(n \log n)$。）于是显然第 $k$ 层合并所有区间时间复杂度为 $O(n \log^2 n)

最后得到这个算法的最终复杂度为 O(n \log^3 n)。 ::::

例题

1.计算 e

对应洛谷题目：https://www.luogu.com.cn/problem/P1729。

是不是感觉很容易了呢......

设我们计算 m 项级数的和后可以使 e 取到小数点后第 n 位，e_m=\sum^{m}_{i=0} \frac{1}{i!}，来点简单的数学推导：

\begin{aligned} e_m &=1+\sum^{m}_{i=1} \frac{1}{i!} \\ &=1+\sum^{m}_{i=1}\frac{1}{1}\prod^{i}_{j=1}\frac{1}{j} \end{aligned}

不妨直接定义：

\begin{cases} S(l,r)=\sum^{r}_{i=l} \frac{1}{i!}\\ P(i)=1\\ Q(i)=i \\ R(i)=1 \end{cases}

直接按照上面的方式计算就完事了。

::::success[参考代码]

由于数据范围小，因此没有用 python 中自带 FFT 的 Decimal 库来加速运算，因为 FFT 常数不小会显得更慢。

代码没必要如此格式化，只是方便理解而已。

import sys,math
sys.set_int_max_str_digits(0)  # 注意到要计算小数点后10000位

def P(x):
    return 1
def Q(x):
    return x
def R(x):
    return 1
def bs(l,r):
    # 核心代码
    if l==r:
        return (P(l),Q(l),R(l))
    m=(l+r)//2
    P1,Q1,R1=bs(l,m)
    P2,Q2,R2=bs(m+1,r)
    QT=Q1*Q2
    RT=R1*R2
    PT=P1*Q2+R1*P2
    return (PT,QT,RT)
def main():
    n=int(input())
    # 多计算一些
    m=n+50

    PT,QT,RT=bs(1,m)

    fenzi=PT-QT
    fenmu=QT

    ansnum=math.floor((fenzi*(10**n))//fenmu)

    ansstr=str(ansnum)
    print("2.")

    for i in range(0,len(ansstr),50):
        l=ansstr[i:i+50]
        output=[l[k:k+10]for k in range(0,len(l),10)]
        print(' '.join(output))
if __name__=="__main__":
    main()

这里是提交记录，可见虽然用 python 实现，却一点也不虚。 ::::

::::info[思考：是否存在方法二？] 辅助函数 R(x) 都等于1了，那么 R(l,r) 不用说也等于 1，那么我们代入就能得到新的合并式子：

\begin{cases} P(l,r)=P(l,m) \cdot Q(m+1,r)+ P(m+1,r)\\ Q(l,r)=Q(l,m) \cdot Q(m+1,r)\\ \end{cases}

新的递推边界自然为

\begin{cases} P(l,r)=P(l) \\ Q(l,r)=Q(l) \\ \end{cases}

于是我们就得到了二进制分裂的一个特殊形式。

代码特别好写，在原通用形式上删掉有关 R(x) 的定义就可以了，这里只给出核心代码和提交记录，在数据范围小的情况下，这也并没有快多少。

def bs(l,r):
    if l==r:
        return (P(l),Q(l))
    m=(l+r)//2
    P1,Q1=bs(l,m)
    P2,Q2=bs(m+1,r)
    QT=Q1*Q2
    PT=P1*Q2+P2
    return (PT,QT)

这种特殊形式也是 numberworld 原文作者所采用的方法，不过只有计算 e 和 \frac{1}{e} 才有用到。 ::::

2.计算 \pi

对应洛谷题目：https://www.luogu.com.cn/problem/P1727。

我们直接考虑使用 Chundnovsky 公式：

\begin{aligned} \frac{1}{\pi}&=12 \sum_{k=0}^{\infty}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k+\frac{3}{2}}}\\ &=\frac{1}{426880\sqrt{10005}}\sum_{k=0}^{\infty}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k}} \end{aligned}

由此得：

\pi=\frac{426880\sqrt{10005}}{\sum_{k=0}^{\infty}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k}}}

~~有一种前面 1+1 后面造飞机的感觉。~~

看上去这么丑的一堆式子，直接计算肯定不行的，二进制分裂就派上用场了！

::::info[注]{open}

本例可能存在多种构造方式，这里主要介绍 numberworld 原文作者的构造方式，仅供参考。

::::

仍然设我们计算 m 项级数的和后可以使 \pi 取到小数点后第 n 位。

设 \pi_m=\frac{426880\sqrt{10005}}{\sum_{k=0}^{m}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k}}}。

那么

\begin{aligned} \pi_m&=\frac{426880\sqrt{10005}}{13591409+\sum_{k=1}^{m}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k}}}\\ &=\frac{426880\sqrt{10005}}{13591409+\sum_{k=1}^{m}(13591409+545140134k)\cdot(-1)^k\cdot\frac{ (6k)!}{640320^{3k}(3k)!(k!)^3} } \end{aligned}

令 F(a) 等于右边带有阶乘的分式，即：

F(a)=\frac{(6a)!}{640320^{3a}(3a)!(a!)^3}

显然 F(0)=1。

由于我们知道：

\begin{aligned} F(a)&=\prod^{a}_{j=1}\frac{F(j)}{F(j-1)}\\ \end{aligned}

于是我们直接大力硬算一波 \frac{F(a)}{F(a-1)}，这个读者可以自己去算。结果是

\frac{F(a)}{F(a-1)}=\frac{(2a-1)(6a-1)(6a-5)}{10939058860032000a^3}

于是：

F(a)=\prod^{a}_{j=1}\frac{(2j-1)(6j-1)(6j-5)}{10939058860032000j^3}

反代回去（这里只拿分母出来看）：

\begin{aligned} &13591409+\sum_{k=1}^{m}(13591409+545140134k)\cdot(-1)^k\cdot\frac{ (6k)!}{640320^{3k}(3k)!(k!)^3}\\ &=13591409+\sum_{k=1}^{m}(13591409+545140134k)\cdot(-1)^k\cdot F(k)\\ &=13591409+\sum_{k=1}^{m}(13591409+545140134k)\cdot(-1)^k\cdot \prod^{k}_{j=1}\frac{(2j-1)(6j-1)(6j-5)}{10939058860032000j^3}\\ &=13591409+\sum_{k=1}^{m}\frac{(13591409+545140134k)\cdot(-1)^k\cdot(2k-1)(6k-1)(6k-5)}{(2k-1)(6k-1)(6k-5)}\cdot \prod^{k}_{j=1}\frac{(2j-1)(6j-1)(6j-5)}{10939058860032000j^3}\\ \end{aligned}

终于化成了二进制分裂可以应用的形式！

现在就非常清晰了！

\begin{cases} S(1,m)=\sum_{k=1}^{m}\frac{(-1)^k (6k)!}{(3k)! (k!)^3} \cdot \frac{13591409+545140134k}{640320^{3k}}\\ P(a)=(13591409+545140134a)\cdot(-1)^a\cdot(2a-1)(6a-1)(6a-5) \\ Q(a)=10939058860032000a^3 \\ R(a)=(2a-1)(6a-1)(6a-5) \end{cases}

由此得到

\pi_m=\frac{426880\sqrt{10005}}{13591409+S(1,m)}

还没完！由于最后求出 S(1,m) 需要一次高精度除法，带回去又需要一次高精度除法，很慢，所以不妨直接代入 S(1,m)=\frac{P(1,m)}{Q(1,m)}，得到：

\begin{aligned} \pi_m&=\frac{426880\sqrt{10005}}{13591409+\frac{P(1,m)}{Q(1,m)}}\\ &=\frac{426880\sqrt{10005}\cdot Q(1,m)}{13591409 \cdot Q(1,m)+P(1,m)} \end{aligned}

完成！

::::info[注]{open} numberworld 原文作者最后把分子的 \sqrt{10005} 移到分母去了，变成：

\pi_m=\frac{1}{\sqrt{10005}} \cdot \frac{4270934400 \cdot Q(1,m)}{13591409 \cdot Q(1,m)+P(1,m)}

计算效率上没什么差别。 :::: ::::warning[特别注意：对于 m 的取值]{open}

如果直接 m=n+k （k 你可以自取正整数，为了保证精度用的。）会慢亿点，然而 Chundnovsky 公式每计算一项可以增加约 14.18 位精度，因此我们应取 m=\lfloor\frac{n}{14}\rfloor+k。

~~这是写给那些像我一样写完后才发现这个事的小傻瓜的。~~ ::::

::::success[参考代码]

同样的，由于数据范围小，不用 Decimal。

import sys,math
sys.set_int_max_str_digits(0)  # 注意到要计算小数点后10000位

def P(x):
    p=-1
    if x%2==0:
        p=1
    return (13591409+545140134*x)*(6*x-5)*(6*x-1)*(2*x-1)*p
def Q(x):
    return 10939058860032000*x*x*x
def R(x):
    return (6*x-5)*(6*x-1)*(2*x-1)
def bs(l,r):
    # 核心代码
    if l==r:
        return (P(l),Q(l),R(l))
    m=(l+r)//2
    P1,Q1,R1=bs(l,m)
    P2,Q2,R2=bs(m+1,r)
    QT=Q1*Q2
    RT=R1*R2
    PT=P1*Q2+R1*P2

    return (PT,QT,RT)
def main():
    n=int(input())
    # 多计算一些
    m=n//14+2

    PT,QT,RT=bs(1,m)

    fenzi=QT

    fenmu=(PT+13591409*QT)

    ansnum=math.floor((426880*fenzi*math.isqrt(10005*10**(n+n)))//(fenmu))

    ansstr=str(ansnum)
    ansstr=ansstr[1:n+1]
    print("3.")

    for i in range(0,len(ansstr),50):
        l=ansstr[i:i+50]
        output=[l[k:k+10]for k in range(0,len(l),10)]
        print(' '.join(output))
if __name__=="__main__":
    main()

提交记录 ::::

这道重要的例题告诉了我们：

对于含有阶乘的分式 F(a)，我们一般先 F(a)=\prod^{a}_{i=1}\frac{F(a)}{F(a-1)}，然后手动爆算出 \frac{F(a)}{F(a-1)}，最后反代回去
如果我们面对 \sum^{\infin}_{\color{red}i=0} T（T 为含字母的代数式），即下标从 0 开始的情况的时候，先将 i=0 代入 T，算出结果 T_0，然后将原式子写成 T_0+\sum^{\infin}_{i=1}T 的形式，总之一定要让下标从 1 开始！
高精度除法很慢，所以当目标能化为关于 S 的分式的时候（比如这道例题），很多情况下都要将 S(1,m)=\frac{P(1,m)}{Q(1,m)} 代入这个分式，转换为关于 P(1,m) 和 Q(1,m) 的分式，避免算出 S(1,m) 再计算。

这也是解决一般二进制分裂问题所必须要注意的。

写在最后

恭喜你看完了这一篇总字数 11k 的文章。
这个算法最先在洛谷上出现是这两位大佬的题解：题解1，题解2。本文章参考了其中的一小部分。其中@bcdmwSjy 大神用 FFT 和 C++ 完成了二进制分裂在计算 \pi 这道题的代码，砍下了这道题非打表最优解！Orz！
初三编程菜鸡，如果文章有什么你不能理解的地方或者表述有纰漏或不严谨的地方，欢迎在评论区中指出！