FFT（快速傅里叶变换）

简介

$\text{FFT}$ 可以在 $O(n\log n)$ 的时间内计算两个序列的卷积/循环卷积。

其最常见的应用是计算多项式乘法。

卷积&多项式乘法

两个数列 $\bold a$ 和 $\bold b$ 的卷积 $\bold c$ 定义为：

$c_k=\sum_{i+j=k}a_ib_j$

类似地定义循环卷积：

$c_k=\sum_{(i+j)\bmod n=k}a_ib_j$

那么两个多项式 $A(x)=\sum\limits_{i=0}^{n}a_ix^i$ 和 $B(x)=\sum\limits_{i=0}^{n}b_ix^i$ 的乘积定义为：

$C(x)=A(x)\cdot B(x)=\sum_{i=0}^{n}\sum_{j=0}^{n}a_ib_jx^{i+j}=\sum_{i=0}^{2n}c_ix^i$

其中，$\bold c$ 是 $\bold a$ 和 $\bold b$ 的卷积。因此多项式乘法本质上就是计算两序列的卷积/循环卷积。

如果直接暴力计算的话，不难发现复杂度是 $O(n^2)$ 的。

怎么优化这个 $n^2$ 的复杂度啊？

用点值表示法试试！

如果用点值表示法，我们发现，如果知道了 $C(x)=A(x)B(x)$ 的足够多的点值，那么就可以计算出 $C(x)$ 的各项系数。

如果直接算出这些点值，需要在 $O(n^2)$ 的时间>_<

所以我们优化了个寂寞啊

至少要低于 $O(n^2)$ 吧qwq......

这里我们考虑一组特殊的点值。

单位根

我们定义 $\sqrt{-1}=i$，即 $i^2=-1$。

方程 $x^n-1=0$ 的 $n$ 个复根 $x$ 被称为 $n$ 次单位根，记作 $\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$。

其中，$\omega_n^k$ 定义为 $e^{\frac{2k\pi i}{n}}=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}$。而 $e^{i\theta}=\cos\theta+i\sin\theta$ 就是欧拉公式，它的证明和 $\text{FFT}$ 没啥关系，就不写啦>_<

例如，$3$ 次单位根有：$\omega_3^0=1,\omega_3^1=\dfrac{\sqrt{3}i-1}{2},\omega_3^2=\dfrac{-\sqrt{3}i-1}{2}$。

我们需要用到单位根的这两个性质：

$\omega_{n}^{k}=\omega_{nd}^{kd},d\in \mathbb N\tag{1}$

事实上，$\omega_{nd}^{kd}=\cos\dfrac{2kd\pi}{nd}+i\sin\dfrac{2kd\pi}{nd}=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}=\omega_n^k$，故 $(1)$ 成立。

$\omega_{2n}^{n+k}=-\omega_{2n}^k\tag{2}$

事实上，$\omega_{2n}^{n+k}=\omega_{2n}^n\cdot\omega_{2n}^k=\omega_2^1\cdot\omega_{2n}^k=-\omega_{2n}^k$，故 $(2)$ 成立。

离散傅里叶变换

设 $\bold a$ 是长为 $n$ 的序列，对于 $k\in [0,n)$，定义：

$b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}$

那么我们称 $\bold b$ 为 $\bold a$ 的**「离散傅里叶变换」**，即 $\text{DFT}$ ，记作 $\bold b=\mathcal{F}(\bold a)$ 。

如果我们将序列 $\bold a$ 与 $\bold b$ 的循环卷积记作 $\bold c=\bold a\ast \bold b$，那么不难发现：

$\mathcal{F}(\bold c)=\mathcal{F(\bold a\ast\bold b)=}\mathcal{F(\bold a)}\cdot\mathcal{F(\bold b)}$

其中 $\cdot$ 表示逐点乘积。

如果我们令多项式 $A(x)=\sum a_ix_i$，不难发现 $b_k$ 就是 $A(x)$ 在 $\omega_n^k$ 处的点值，即：$A(\omega_n^k)=b_k$。

也就是说，想要计算 $\bold a$ 的 $\text{DFT}$，那么我们就需要 $A(x)=\sum a_ix^i$ 的 $n$ 个点值。

同样地，只要我们知道了 $\bold a$ 的 $\text{DFT}$，我们同样可以利用这 $n$ 个点值表示出 $A(x)$，计算出 $A(x)$ 的各项系数。

但是，如果朴素地计算 $\text{DFT}$ 也就是这 $n$ 个点值，那复杂度是 $O(n^2)$ 的>_<

于是就出现一个问题：如何快速计算一个序列 $\bold a$ 的 $\text{DFT}$ 呢？

蝴蝶操作

我们考虑将多项式 $A(x)=\sum\limits_{i=0}^{n-1}$ 的 $n$ 项系数按奇偶性分类：

不妨设 $n=2m$，则

$\begin{aligned} A(x)&=\sum_{i=0}^{n-1}a_ix^i=\sum_{i=0}^{m-1}a_{2i}x^{2i}+\sum_{i=0}^{m-1}a_{2i+1}x^{2i+1}\\ &=\sum_{i=0}^{m-1}a_{2i}x^{2i}+x\sum_{i=0}^{m-1}a_{2i+1}x^{2i} \end{aligned}$

我们记

$A_0(x)=\sum_{i=0}^{m-1}a_{2i}x^{i}\\ A_1(x)=\sum_{i=0}^{m-1}a_{2i+1}x^{i}$

那么就得到了两个次数不超过 $m$ 的多项式 $A_0(x)$ 和 $A_1(x)$，并且有：

A(x)=A_0(x^2)+xA_1(x^2)\tag{#}

现在，如果我们有了 $A_0(\omega_m^k)$ 和 $A_1(\omega_m^k)$ 的值，带入 $(\#)$ 式，计算得：

$A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_{n}^kA_1(\omega_n^{2k})=A_0(\omega_m^{k})+\omega_{n}^kA_1(\omega_m^k)\\ A(\omega_n^{k+m})=A_0(\omega_n^{2k+n})+\omega_{n}^{k+m}A_1(\omega_n^{2k+n})=A_0(\omega_m^{k})-\omega_{n}^kA_1(\omega_m^{k})$

诶，只利用 $A_0(\omega_m^k)$ 和 $A_1(\omega_m^k)$ 的值，我们居然直接算出了 $A(x)$ 的两个点值！

也就是说，只要我们有了 $A_0(x)$ 和 $A_1(x)$ 在 $\omega_m^0,\omega_m^1,\cdots,\omega_m^{m-1}$ 处的点值，就可以在 $O(n)$ 的时间内计算出 $A(x)$ 在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 处的点值！

这，这不就是个递归的事嘛！

我们递归下去计算 $A_0(x)$ 和 $A_1(x)$ 的点值，每次将问题规模减半，合并子问题又是 $O(n)$ 的，那么复杂度就是：

$T(n)=2T(n/2)+O(n)=O(n\log n)$

好耶！我们成功在 $O(n\log n)$ 的时间内计算出了 $A(x)$ 在 $n$ 次单位根上的点值，也就是其系数数列 $\bold a$ 的 $\text{DFT}$ ！！

这个过程通常被称为**「蝴蝶操作」**。

现在我们就可以在 $O(n\log n)$ 的时间内计算出 $\bold a$ 和 $\bold b$ 的 $\text{DFT}$。

再根据 $\mathcal{F}(\bold a\ast \bold b)=\mathcal{F(\bold a)}\cdot\mathcal{F(\bold b)}$，我们可以在 $O(n)$ 的时间内计算出它们的卷积的 $\text{DFT}$。

但是我们要算的是 $\bold a\ast \bold b$ 而不是它的 $\text{DFT}$ 啊=_=

如果把 $\text{DFT}$ 看作点值，然后用拉格朗日插值，把 $\mathcal{F(\bold a\ast\bold b)}$ 转化为 $\bold a\ast\bold b$，那复杂度还是 $O(n^2)$ >_<

那现在的问题就是：我们知道了序列 $\bold c$ 的 $\text{DFT}$，如何在低于 $O(n^2)$ 的时间内计算出 $\bold c$？

离散傅里叶变换的逆变换

对于长为 $n$ 的序列 $\bold a$，若其 $\text{DFT}$ 为 $\bold b$，即：

$b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\tag{1}$

我们证明 $(1)$ 的逆变换，即其 $\text{IDFT}$ 为：

$a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{2}$

我们将 $(1)$ 中 $b_k$ 的表达式带入 $(2)$，得到：

$\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}=\dfrac{1}{n}\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_i\omega_n^{ij}$

交换一下求和顺序（其实就是提取公因式），得到：

$\dfrac{1}{n}\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}=\dfrac{1}{n}\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}$

我们考虑后面的那个小可爱 $\sum\limits_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}$

当 $j=k$ 时，有：

$\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\sum_{i=0}^{n-1}\omega_n^{ik-ik}=\sum_{i=0}^{n-1}1=n$

当 $j\neq k$ 时，则有：

$\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\sum_{i=0}^{n-1}\omega_n^{i(j-k)}$

我们发现这东西是个类似于等比数列求和的形式，于是用一下等比数列求和公式：

$\sum_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum_{i=0}^{n-1}\left(\omega^{j-k}\right)^i=\dfrac{1-\left(\omega^{j-k}_{n}\right)^n}{1-\omega_n^{j-k}}=\dfrac{1-\left(\omega_n^n\right)^{j-k}}{1-\omega_n^{j-k}}=\dfrac{1-1}{1-\omega_n^{j-k}}=0$

于是只要 $j\neq k$，$\dfrac{1}{n}a_j\sum\limits_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_{n}^{-ki}$ 整个就是个 $0$，不会对答案产生任何影响。

因此这个式子就被化简成了：

$\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}=\dfrac{1}{n}\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{ij}\cdot\omega_n^{-ki}=\dfrac{1}{n}a_k\cdot n=a_k$

即

$a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{2}$

证毕。

类似地，我们也可以证明 $(1)$ 和 $(2)$ 是等价的。

我们称 $(2)$ 为**「离散傅里叶变换的逆变换」**，即 $\text{IDFT}$。如果朴素地按照表达式计算 $\text{IDFT}$，那么复杂度还是 $O(n^2)$ 的。

但实际上，观察 $\text{IDFT}$ 的表达式，对比一下 $\text{DFT}$ 的表达式：

$b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\tag{DFT}\\$

$a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\tag{IDFT}$

实际上只需要将计算 $\text{DFT}$ 时用到的「蝴蝶操作」中的 $\omega_n$ 全部替换成 $\omega_n^{-1}$ ，类似地递归下去，同样可以在 $O(n\log n)$ 的时间内计算 $\text{IDFT}$。唯一不同的是算完之后需要乘上 $\dfrac{1}{n}$ 。

其实再多看两眼一下 $\text{IDFT}$ 的表达式，可以发现只需要翻转一下 $\bold b$ 的后 $n-1$ 项，再做一次 $\text{DFT}$ 然后除一下，也相当于做完了 $\text{IDFT}$。

所以 $\text{IDFT}$ 同样可以在 $O(n\log n)$ 的时间内被计算。

总结

总结一下我们的 $\text{FFT}$ 流程：

1. 利用 $\text{DFT}$，在 $O(n\log n)$ 的时间内分别计算出 $\bold a$ 和 $\bold b$ 的 $\text{DFT}$。
2. 把算出来的出来的这两个 $\text{DFT}$ 相乘，在 $O(n)$ 的时间内得到 $\bold c$ 的 $\text{DFT}$。
3. 利用 $\text{IDFT}$，在 $O(n\log n)$ 的时间内，把 $\bold c$ 的 $\text{DFT}$ 转化为 $\bold c$ ，就计算出了 $\bold a$ 和 $\bold b$ 的卷积。

综上，我们在 $O(n\log n)$ 的时间内计算出了两个序列的卷积！

FFT 代码实现&注意事项&技巧（递归版）

typedef complex<double> comp;
#define pi acos(-1)//不要手打 pi QAQ，最好用 acos(-1)
void dft(comp *a,int p){
    if(p==1)return;//递归边界
    
    int q=p>>1;
    comp a1[q],a2[q];
    for(int i=0;i<q;i++){
        a1[i]=a[i<<1];
        a2[i]=a[(i<<1)|1];
    }//按奇偶分类
    
    dft(a1,q);dft(a2,q);//分治
    
    for(int i=0;i<q;i++){
        comp omega=comp(cos(2*pi*i/p),sin(2*pi*i/p));
        a[i]=a1[i]+omega*a2[i];
        a[i+q]=a1[i]-omega*a2[i];//照搬我们推出来的式子qwq
    }
}

void idft(comp *a,int N){
    dft(a,N);
    reverse(a+1,a+N);
    for(int i=0;i<=N;i++)a[i]=(a[i]+comp(0.5,0))/((double)N);//输出需要是整数
}

注意到如果我们按照上述的操作进行 $\text{FFT}$，只能处理长度为 $2^k$ 的数列的卷积，其中 $k\in \mathbb N^*$。

因此，我们一般会选取一个最小的 $2^k\ge n+m$ 作为序列长度，多出来的地方就看做 $0$，然后再进行 $\text{FFT}$。

所以，实际上，在计算卷积的时候，我们通常会：

int k=1;for(;k<=n+m;k<<=1);

然后再：

dft(x,k);
dft(y,k);
for(int i=0;i<=k;i++)z[i]=x[i]*y[i];
idft(z,k);

最后才是：

for(int i=0;i<=n+m;i++)cout<<(int)z[i].real()<<" ";

此时，由于我们寻找的是最小的 $2^k\ge n+m$ ，因此数组并不能直接开成 $10^6$ 这么大，而要开到最小的 $2^k\ge 10^6+10^6$ ，大概是 $2^{21}=2097152$ qwq。

方便起见，开成 $2100000$ 也是没有问题的QAQ。

整个代码十分简洁易懂（所以窝并没有加多少注释），就是把我们之前推导的式子照搬过来了qwq。

迭代实现 FFT

然而递归版常数大到爆炸，尤其是我们递归的时候居然还要递归一下开一遍内存，很容易空间爆炸>_<

（虽然递归版开 O2 就能直接艹过板子题）

而且，$\text{FFT}$ 作为一个经常要使用的算法，如果它的常数过大导致整个程序时间爆炸是很不好的qwq。

能不能把 $\text{FFT}$ 改写成迭代形式啊？

画出递归树：（图片咕了）

在这个递归树中，如果我们能够快速地处理出叶子节点所构成的序列 $A$（本例中为 $A=\{a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7\}$，那么就可以自底向上地实现整个 $\text{FFT}$ 操作。

具体地，我们现在有了一个序列 $A$，就是这个递归树的叶子节点所构成的序列。

首先我们从前往后，两个两个地取出 $A$ 中元素，然后对于每两个元素，算出它们的 $\text{DFT}$，这样就有了叶子节点上一层的 $n/2$ 个 $\text{DFT}$。

然后用这 $n/2$ 个 $\text{DFT}$ ，可以再算出来再往上一层的 $n/4$ 个 $\text{DFT}$，直到最后算出来递归树的「根」处的 **整个序列的 $\text{DFT}$ **。

假设我们已经有了这个序列 $A$ 的初始值，那么就能写出这样的代码：

typedef complex<double> comp;
comp A[MAXN];
void DFT(comp *a,int N){    //N 是 2 的正整数次幂。
    Pre(a,N);
    //假装我们已经通过一些神奇操作在低于 O(n log n) 的时间内预处理出了 A[] qwq。
    int n=log2(N);
    for(int i=1;i<=n;i++){
        int m=1<<i;         //m=2^i，就是递归时这一层中每组的长度
        for(int j=0;j<N;j+=m){
            for(int k=0;k<m/2;k++){
                comp omega=comp(cos(2*k*pi/m),sin(2*k*pi/m));
                comp x=omega*A[j+k+m/2],y=A[j+k];
                A[j+k]=y+x;A[j+k+m/2]=y-x;
            }
        }
    }
    //算完之后的 A 就是 a 的 DFT 了qwq。
    //看上去是三层的循环qwq
    //但其实第二层和第三层合起来才只是一个 O(n)，因此实际上总复杂度是 O(n log n)。
}

所以，到底怎么在一个较好的时间复杂度内确定 $A$ 啊？

直接模拟递归过程显然不行qwq......虽然是 $O(n\log n)$ 的，但是递归下去常数不还是爆炸>_<

位逆序置换

我们仔细观察一下递归的过程，还以 $7$ 次多项式的八个系数为例：

• 一开始，这八个系数好好地放在那里：$\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}$。
• 我们把这个系数按照奇偶分了个类：$\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\}$。
• 又把每一个按照奇偶分了个类：$\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\}$。
• 最后，变成了：$\{a_0\},\{a_4\},\{a_2\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\}$。

观察 $(0,4,2,6,1,5,3,7)$ 的二进制表示：$(000,100,010,110,001,101,011,111)$

诶，把这个二进制表示翻转过来：

000 100 010 110 001 101 011 111
翻转后: 000 001 010 011 100 101 110 111

这，这不就是 $(0,1,2,3,4,5,6,7)$ 的二进制表示嘛！

实际上这其中的道理是很好解释的：考虑我们分组时的分组依据。

第一次我们将奇偶分了个组，把偶数放到前面，把奇数放到后面，也就相当于把二进制表示中末尾为 $0$ 的放到前面，末位为 $1$ 的放到后面。

再想想把 $(0,1,2,3,4,5,6,7)$ 的二进制表示：显然前面四个的首位是 $0$，后面四个的首位是 $1$，毕竟这个序列是递增的嘛QAQ。

也就是说，我们在末位上把 $0$ 全都提到了前头，把 $1$ 全都扔到了后头，而 $(0,1,2,3,4,5,6,7)$ 则是在首位上把 $0$ 提到前头，把 $1$ 扔到后头，一模一样啊QAQ！

那同理，我们后面的操作也都是把每组中的 $0$ 提到前头，把 $1$ 扔到后头；这恰巧相当于 $(0,1,2,3,4,5,6,7)$ 这一顺序排列在每组中的操作。

如果我们令 $\text{rev}(n)$ 为 「将 $n$ 的二进制表示翻转过来之后得到的值」，根据上面的结论，我们就有：

$A_{\text{rev}(i)}=a_i$

那我们就可以这么处理序列 $A$

for(int i=0;i<n;i++)A[rev(i)]=a[i];

其中 rev(n) 是一个计算「将 $n$ 的二进制翻转后得到的值」的函数。

计算 rev(n) 显然可以做到 $O(\log n)$，毕竟 $n$ 在二进制下的位数也一共只有 $O(\log n)$ 个，随便你怎么搞都行嘛qwq。

这样的话就做到了 $O(n\log n)$ 计算序列 $A$。

实际上，如果 $n$ 在二进制表示下的后 $L$ 位（本例中 $L=3$），那么将 $n$ 的这 $L$ 位翻转，其实就相当于：把前面的 $L-1$ 位拽出来，翻转一下，再把最后一位搞到最前头去qwq！

也就是说，如果我们知道了 $\text{rev}(\texttt{n>>1})$ 即 $\text{rev}(n/2)$ 的值，那么就可以直接 $O(1)$ 算出来 $\text{rev}(n)$ 的值。

来构造一下这个递推式吧qwq。

首先要把前 $L-1$ 位单拉出来并翻转：rev[i>>1]。

注意，这里我们相当于把前面的 $L-1$ 全部右移一位，把最后一位「顶掉」，这会导致最前面空出来一个 $0$。

那翻转之后这个最前面的 $0$ 就被翻到了最后一位，也就是说，我们需要再右移一位，把这个 $0$ 顶掉，即：rev[i>>1]>>1。

现在，这个 rev[i>>1]>>1 就是 rev[i] 的后面 $L-1$ 位翻转得到的结果了，而且位置是没有问题的qwq。

现在只需要把 i 的最后一位搞到前头去就行啦！

i 的最后一位，就是 i&1；搞到前头去就相当于左移 $L-1$ 位（这样它变成了第一位），再和刚才的 rev[i>>1]>>1 「或」一下即可qwq。

那么递推式就是：rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1))

因此，我们就利用递推做到了 $O(n)$ 预处理 $\text{rev}(i)$，然后 $O(n)$ 计算序列 $A$ qwq！

FFT 代码实现&注意事项&技巧（非递归版）

此时我们需要算出这些数在二进制下的位数，也就是把这一步改一下：

int k=1,L=0;
for(;k<=(n+m);k<<=1,L++)//L 就是我们需要考虑的位数。

其他的都基本没有什么区别qwq。

typedef complex<double> comp;

comp A[MAXN];
int rev[MAXN],k,L;

void init(){
	rev[0]=0;
	for(int i=0;i<k;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
}//预处理 rev[]

void dft(comp *a,int N){    //N 是 2 的正整数次幂。
	init();
	for(int i=0;i<k;i++)A[rev[i]]=a[i];
    int n=log2(N);
    for(int i=1;i<=n;i++){
        int m=1<<i;         //m=2^i，就是递归时这一层中每组的长度
        for(int j=0;j<N;j+=m){
            for(int l=0;l<m/2;l++){
                comp omega=comp(cos(2*l*pi/m),sin(2*l*pi/m));
                comp x=omega*A[j+l+m/2],y=A[j+l];
                A[j+l]=y+x;A[j+l+m/2]=y-x;
            }
        }
    }
    for(int i=0;i<k;i++)a[i]=A[i];
}

void idft(comp *a,int N){
    dft(a,N);
    reverse(a+1,a+N);
    for(int i=0;i<=N;i++)a[i]=(a[i]+comp(0.5,0))/((double)N);//输出需要是整数
}

最后对比一下递归版和非递归版：

（上面是递归版开 O2卡过，下面是非递归版不开 O2）

我知道我对照实验变量不唯一，不管了不管了，就当是突出一下非递归版吊打递归版>_<

其中非递归版最慢的一个点跑了 $1.61s$，递归版是 $1.78s$ qwq。

用 FFT 切点题？

只过掉板子题怎么行！一个蓝题根本不够嘛！说好的爆切蓝紫黑呢？

1. P1919 【模板】A*B Problem升级版（FFT快速傅里叶）

高精度乘法？

把一个数的 $10^i$ 位置上的数看成 $a_i$，就有了一个多项式 $f(x)=\sum a_ix^i$，取 $x=10$ ，这个多项式的值就是这个数。

因此本质上还是计算多项式乘法的过程。

2. P3338 [ZJOI2014]力

给一个长为 $n$ 的序列 $\{q_n\}$，定义

$F_i=\sum_{i>j}\dfrac{q_iq_j}{(i-j)^2}-\sum_{i<j}\dfrac{}{(i-j)^2}$