背景

Fast Multipole Method (FMM) 是求解$n$体问题中的受力所使用的技术，被称为20世纪十大算法之一。既然名头这么大，为什么我直到大三的下学期才在营销号的推送里第一次听说呢？并不是我孤陋寡闻了，而是因为这个“十大算法”其实并非正经的计算机科学研究机构评选，而是SIAM（工业与应用数学学会）选出来的，当然就是以计算数学中的算法为主了（虽然混进去了一个快速排序，还有FFT、Monte-Carlo等等应用场景更加广泛的算法）。

虽然大三才第一次接触，我第一次阅读FMM的wiki page的时候就立刻被吸引了，相见恨晚。通俗地说，FMM解决的就是“任意位置的FFT”：假设空间中有位于$x_1\dots x_n$的$n$个质点，其质量为$m_1\dots m_n$，则它们的引力势能为

如果$x_1\dots x_n$分布在一个均匀的点阵上，使用FFT就能轻松计算这一卷积。可惜，$n$体问题中的物体当然不会乖乖排列成一个矩形点阵，FFT对此就无能为力了。难道我们就只能用$O(n^2)$的暴力求和来龟速计算吗？然而先哲们注意到虽然$x_i$没有规律，可卷积核$\frac{1}{|x-y|}$有良好的性质，这就给了我们加速的机会。这个加速过程运用了许多优美的数值分析思想，这让当时沉浸在数值分析的美感中的我深感陶醉。虽然如此，本文作为一篇入门文章，仅会涉及其中的一小部分技巧，希望能通过最小的理解难度展示FMM的核心思想。

我们也不会考虑复杂的三维问题，而是使用一个一维的例子：[ZJOI2014] 力。在这个简化的例子中，$n$个电荷分别位于数轴上的$1\dots n$处，电荷为$q_1\dots q_n$，所求的场强（在静电力常量为一的单位选择下）即为

虽然这个问题与$n$体问题中的引力势能求解不同，但FMM同样也能解决；事实上，FMM能处理相当多物理问题中所需的求和。当然，这个问题中的电荷是等距分布的，所以FFT即可解决（毕竟这就是我当年学习FFT时的入门题）。然而下文将展示，FMM不仅也能完成求解，还比FFT更加优秀。

我的思考

最开始，我对于FMM到底是怎么实现的一无所知，只知道存在一种叫做“FMM”的算法，能够快速处理与$\frac{1}{(x-y)^2}$核的卷积。我决定靠自己的力量，想出一种我自己的“FMM”，这样会对理解该算法有很大的帮助。（才不会是由于wiki和一堆资料看完了，chatGPT也问了，还是啥都没懂呢）

我不知道FMM，那么我知道什么呢？我知道$\frac{1}{(x-y)^2}$作为一个解析函数，可以被泰勒展开近似。我也知道，与多项式核的卷积是平凡的：因为$(x-y)^k$可以被展开为一系列的$y^j$与$x^{k-j}$乘积之和，只要预处理出$\sum_{y} q_y y^j$，然后乘上组合数再乘上$x^{k-j}$再加起来，就可以实现与多项式核的卷积。我不熟悉的是什么呢？$\frac{1}{x^2}$在$y$处的泰勒展开收敛半径为$|y|$，必须限制在$\frac{|y|}{2}$半径内才能做到只用$\log_2{\epsilon^{-1}}$项达到所需的精度要求$\epsilon$。这就说明无论我们在哪里泰勒展开，都无法用几个展开覆盖所有的$q_y$，这是我之前没有遇到过的情况。虽然如此，只需要在$y,\frac{y}{2},\frac{y}{4},\dots$这$\log n$个地方展开，每个地方展开$O(\log \epsilon^{-1})$阶，就能覆盖到所有的$q_y$了。

虽然我已经成功得到了一个非平凡的做法，并且已经足够通过这道题了（很难想象$n\le 10^5$、时限3s的题标算会是一个log，2014年的机器有这么慢吗，现在我$O(n^2)$的暴力都跑不到三秒），但我对此并不满意。为什么呢？因为我已经知道了FMM算法的复杂度是$O(n\log \epsilon^{-1})$，而我想出来的这个做法是$O(n\log n\log \epsilon^{-1})$。因此我下面开始想怎么优化这个做法。首先，我现在查询的区间总量为$O(n \log n )$种，每个都要$O(P)$的计算量（下文记泰勒展开阶数$P=O(\log \epsilon^{-1})$），这就已经是两个log的计算量了，基本没有优化的余地。但这显然是不必要的，我们如果建一棵线段树，只查询节点所代表的那些区间，这些信息也够用了。更进一步地，我们可以仅保留大小至少为$P$的那$O(n/P)$个节点，大小$P$以内的部分暴力查询，这样维护的信息量就只有$O(n)$了，才有希望做到$O(nP)$。（这也是一个常见的优化思路了。在这个具体的例子中，小块暴力可以降低时间复杂度；并且在基本上所有的分治算法中，都可以这样做来优化常数）

上面都是一些抽象的叙述，用来解释我的思路；从这里开始，我会更具体地介绍这一算法是如何实现的。首先，我们约定用$x$表示一个想要求值的位置（“汇”），而$y$表示一个产生贡献的位置（“源”）。根据基本的数学知识，

在上式中把原点平移到任何一点$c$处当然也成立（即变量代换$x\to x-c,y\to y-c$）。对于代表区间$[l,r]$的节点$p$，我们想要用区间中点$c=\frac{l+r}{2}$处的泰勒展开表达离$c$足够远的$x$处函数值；具体来说，只要$|x-c|\ge r-l\ge 2|y-c|$， 下式在精度范围内成立：

所以我们只要保存$M_{p,k}=\sum_{y\in [l,r]}q_y(y-c)^k$即可在$O(P)$时间内计算该节点对$x$的贡献。对于大小为$O(P)$的这$O(n/P)$个节点，我们用$O(nP^{-1}P^2)=O(nP)$的时间从$q_y$中暴力求和计算$M_{p,k}$即可；而对于更大的$O(n/P)$个节点，我们可以用多项式平移将它的两个子节点的$M$数组合并为自身的$M$数组，复杂度为$O(nP^{-1}P^2)=O(nP)$。什么，你说可以用FFT算多项式平移？但我们需要计算出的信息量至少为$O(nP)$，已经没有加速的可能了，所以直接暴力与组合数卷积实现多项式平移即可。

虽然这棵线段树总算是建起来了，可我失望地发现总体的复杂度其实并没有发生变化：对于一个$x$，仍然要查询$O(\log n)$个节点，每个节点上的计算量还是$O(P)$的，仍然是两个log。这说明这仍然不是FMM算法：我们还需要设计一个同时查询多个节点的方法，才能避免多次遍历线段树带来的$O(n\log n)$的代价。但这一想法也并非一无是处：由于它是一棵线段树，天生可以支持单点修改、区间修改、区间查询等功能，这应该是真正的FMM算法所无法做到的，尽管这些OI题目设计思路似乎在实际科研中没啥应用。

我的独立思考就到此为止了：我无法想出该如何同时查询$n$个节点。但这些努力并非徒劳，因为chatGPT接下来的回复总算能看懂了！虽然这些回复包含许多错误，但确实也混有有价值的信息。我在chatGPT的提示下补充查询过程，就构成了完整的FMM算法。

更快的查询

为了同时查询多个位置上的求和，我们需要注意到所求值也可以被视为一个$x$的函数，也可以被泰勒展开：这样，只需要查询一个点上的求和结果，就能得到其附近的所有$x$处的求和结果了。我们之前假设$x$离$c$足够远、$y$离$c$足够近，因为要用$c$近似$y$；而这里需要反过来，假设$x$离$c$足够近、$y$离$c$足够远，因为要用$c$近似$x$。在之前的近似公式中交换$x$和$y$就得到了

其中$c$为$x$所在节点代表的区间$[l,r)$的中点，需要$|y-c|\ge r-l$保证该近似在精度范围内成立。这里为什么改成了左闭右开区间呢？因为想要让后续的过程更简单，我们必须保证精确的二等分，即$[l,r)$的两个子节点是$[l,c)$和$[c,r)$，而非线段树里常见的划分方法。使用左闭右开区间保证了任何一点的值都不会被算两次。

由于只要$|y-c|\ge r-l$该近似公式就成立，而且已经使用了精确的二等分，只要一个节点$p’$与$p$不相邻，$p’$包含的所有$y$的贡献就都能被近似计算。此处“相邻”的定义为$[l(p’),r(p’))\cap [2l-r,2r-l)=\varnothing$。我们记这些远处节点的$y$贡献之和为$L_p$，于是计算某个叶节点$p$包含的全部$x$处的答案就只需要使用$L_p$和$[2l-r,2r-l)$中的$O(P)$个值，查询复杂度总算降到了$O(nP^{-1}\cdot P^{2})=O(nP)$。所剩的问题只有计算$L$了。

根据这个交换$x$和$y$的近似公式，远处节点的$y$贡献之和为

因此$L_p$的具体形式为$L_{p,k}=\sum_{y\in p’, p’\textrm{ is far from }p} q_y(y-c)^{-k-2}$。大致上讲，$L$是所有足够远的节点$p’$上的$M_{p’}$之和，可$M_{p’}$是$y-c_M$的多项式，而$L$是$(y-c_L)^{-1}$的多项式，该如何从$M$算出$L$呢？根据二项式定理，有

由于$|y-c_M|\le \frac12 |c_M-c_L|$，这一级数也只需截取$P=\log_2(\epsilon^{-1})$项即可达到精度要求。这就给出了从$M$计算$L$的公式，该步骤在FMM算法里被称为M2L。（之前合并两个子节点的$M$的步骤被称为M2M）

虽然$L$已经可以计算了，我们已经可以批量查询$O(P)$个$x$上的求和了，可不幸的是算法的复杂度并没有发生任何改变：每个叶节点要查询$O(\log n)$个节点的$M$值，需要$O(P^2\log n )$的计算量来计算单个$L$，总计算量还是$O(nP^{-1}\cdot P^2\log n )=O(nP\log n)$。真正精妙的是，$L$总是可以继承父节点的值：只要某区间离父节点足够远，它离子节点就一定足够远！不过要注意一下，父节点的$L$是$(y-c_f)^{-1}$的多项式，而子节点的$L$是$(y-c)^{-1}$的多项式，我们还需要一个“有理式平移”的过程：二项式定理给出

由于$|c_f-c|\le \frac12 |y-c_f|$，我们还是只需要截取$P$项。因此从父节点的$L$可以在$O(P^2)$时间内算出子节点的$L$，这一步骤被称为L2L。我们有了M2L和L2L，可到底哪些区间里的$M$是要通过M2L贡献进$L$的呢？毕竟大部分的区间已经通过父节点的$L$纳入贡献了。为了弄清这一点，我模拟了$[0,8)$上的算法执行流程：

[0,8)
[0,4)                   [4,8)
[0,2)       [2,4)       [4,6)       [6,8)
[0,1) [1,2) [2,3) [3,4) [4,5) [5,6) [6,7) [7,8)

L[0,2): M[4,8)
L[2,4): M[6,8)
L[4,6): M[0,2)
L[6,8): M[0,4)

L[0,1): M[2,4) + M[4,8) = L[0,2) + M[2,4)
L[1,2): M[4,8) + M[3,4) = L[0,2) + M[3,4)
L[2,3): M[0,1) + M[6,8) + M[4,6) = M[0,1) + L[2,4) + M[4,6)
L[3,4): M[0,2) + M[6,8) + M[5,6) = M[0,2) + L[2,4) + M[5,6)
...

模拟之后这一问题的答案就变得明了了：使用二叉堆的存储格式时（即$p$的左儿子位于$2p$，右儿子位于$2p+1$），

• 如果$p$是左儿子，添加$p/2+1$和$p-2$处的$M$值（整数除法；仅未越界时适用，下同）；
• 如果$p$是右儿子，添加$p/2-1$和$p+2$处的$M$值。

这就是完整的一维FMM算法了，总结如下：

1. 使用精确的二等分建立线段树，仅建立那些大小至少为$T=O(P)$的节点。预处理出$2P$以内的组合数。
2. 使用M2M过程自下而上计算出所有节点上的$M$值。
3. 使用L2L和M2L过程自上而下计算出所有节点上的$L$值；使用叶节点上的$L$值和近邻处的暴力求值计算答案。

结果

[ZJOI2014] 力一题的数据范围为$0\le q_i\le 10^9$，绝对误差要求是$10^{-2}$。这说明理论上应该要取$P=\lceil \log_2(10^{11})\rceil =37$才能通过；事实上我自己想出来的两个log的实现确实是取的$P=37$。但真正的$O(nP)$的FMM算法仅取$P=32$就够了，我猜这是由于相邻的节点间的贡献在完整的FMM算法里是通过暴力计算的，只有间隔较远的点计算精度和$P$有关，而由于平方反比，间隔$\ge P$的点贡献较小。

而关于小块暴力的阈值$T$，我实测下来$T\in [P,4P]$的用时都几乎差不多，于是选择了$T=2P$。

完整的代码见fmm_1d.cpp。该代码使用未优化的IO，但已是此题的最快代码之一了；而使用了和此前的本题最优解相同的快速IO之后，FMM就成为了遥遥领先的最优解。为什么$O(n\log \epsilon^{-1})$会吊打$O(n\log n)$呢？因为此题需要两次卷积，即至少4次的FFT，而且至少是$2n$点的FFT，如果$n$略大于$2^k$就会是接近$4n$点的FFT。虽然FMM也不能说常数很小，但比起8到16倍的$n\log n$还是好多了。

虽然此题已经完满地解决了，但这只是FMM的冰山一角。仍有许多的问题尚未解决：

• 这只是1维的情形，科研中用到的2维和3维又该如何处理呢？
• 我使用了一点的泰勒展开来近似核函数$(x-y)^{-2}$，可wiki上面明明说是用Chebyshev节点上的插值啊？别的地方也有说用球谐展开的？
• 如果把$(x-y)^{-2}$替换为更复杂的核，本文中的核心技巧$(x-y)^{-2}=x^{-2}(1-y/x)^{-2}$将核函数转换为一元函数进行分析不就行不通了？
• FFT有其逆运算IFFT，FMM作为和FFT很接近的线性运算，其逆运算是什么？换句话说，如果我们知道$n$体问题中势能场的$n$个点值，能快速求出这$n$个物体的质量吗？

为了进行科研，我需要解决所有的这些问题，这将是后续的博文的内容，见FMM进阶应用：计算二维点对距离和。

• Previous
  杀戮尖塔无限合集
• Next
  FMM进阶应用：计算二维点对距离和