浅谈快速傅里叶变换（FFT）

多项式乘法

我们定义两个次多项式$F(x),G(x)$，有：

$$F(x)=\sum f_i\times x^i\\ G(x)=\sum g_i\times x^i$$

其中$\{f_i\},\{g_i\}$表示两个多项式的系数

显然当一个多项式的系数被确定后这个多项式就被确定了，我们把这称之为多项式的系数表达

两个多项式可以相乘得到一个结果，比如：

$\quad(x^2+2x+3)(2x^2+3x+4)\\ =2x^4+3x^3+4x^2+4x^3+6x^2+8x+6x^2+9x+12\\ =2x^4+7x^3+16x^2+17x+12$

我们令$H(x)=F(x)\times G(x)$，因为$H(x)$中的$x^k$只能由$F(x)$中的$x^i$和$G(x)$中的$x^{k-i}$相乘得到，所以我们有：

$$h_k=\sum\limits_{i=0}^{k}f_i\times g_{k-i}$$

同样，这里的$\{h_i\}$表示$H(x)$的系数

可以观察到按照这个定义式暴力计算的复杂度是$O(n^2)$级别的，很多时候难以接受，所以我们需要一种高效的算法

FFT

快速傅里叶变换(Fast Fourier Transfromation)简称FFT，可以实现在$O(n\log n)$时间里计算多项式乘法

对于一个多项式$F(x)$，我们把它看作一个关于$x$的函数，显然可以得到它的图像。带入一个值$x_0$便可以得到一个对应的$y_0$，$(x_0,y_0)$便构成了一个有序数对表示点的坐标

结论：$n$个点的坐标可以唯一确定一个$n-1$次函数

感性证明：

对于一个$n-1$次函数，我们带入$n$个值得到有$n$个方程的方程组，此时有$n$个未知数（系数），解方程即可求出系数，系数确定函数（多项式）也随之确定

所以当我们有了$n$对有序实数时，原多项式$F(x)$的表达式就唯一确定，这也就意味着我们可以用$n$个有序实数对等效表示一个多项式，我们将这一组有序实数对称为这个多项式的点值表达

我们现在给出两个多项式$F(x),G(x)$，令$H(x)=F(x)\times G(x)$

我们将数列$\{x_0,x_1\dots x_{2n}\}$分别带入$F(x),G(x)$，得到两组有序数对$(x_0,y_0),(x_1,y_1)\dots (x_{2n},y_{2n})$和$(x_0,y_0^{'}),(x_1,y_1^{'})\dots (x_{2n},y_{2n}^{'})$

显然$F(x_0)\times G(x_0)=y_0\times y_0^{'}$，而这个值就等于$H(x_0)$

将两组有序实数对都这样相乘，我们会得到一组新的有序实数对$(x_0,y_0^{''}),(x_1,y_1^{''})\dots (x_{2n},y_{2n}^{''})$，而这组有序实数对就和$H(x)$对应，如此便可以求出$H(x)$的系数，从而确定$H(x)$的表达式。这就是FFT的核心思想：

将多项式转变为点值表达，相乘求出乘积多项式的点值表达，最后再转变为系数表达

将系数表达转化为点值表达的过程称为DFT，将点值表达转化为系数表达的过程称为IDFT

复数

我们定义形如$a+bi$的数为复数，其中$i=\sqrt{-1}$。我们定义$a$为该复数的实部，$bi$为该复数的虚部

所有实数都可以在一条坐标轴上找到唯一一个对应的点，我们再添加一条竖着的数轴表示虚数，即可使所有复数同这个新形成的平面直角坐标系里的每个点一一对应：

比如上图中这个点$A$即可表示复数$3+2i$

复数可以进行四则运算：

复数相加：实部和虚部分别相加

如：$(3+2i)+(5+4i)=8+6i$

复数相减：取相反数再相加

复数相乘：像多项式一样相乘

如：

$$\begin{equation} \begin{aligned} (3+2i)(5+4i)&=15+12i+10i+8i^2\\ & =7+22i \end{aligned} \end{equation}$$

注意这里$i^2=-1$

复数相除在FFT里暂时用不到，故暂不说明

这里我们需要注意复数相乘在刚才的平面直角坐标系里面的体现：

（该图来自command_block大佬/bx/bx）

我们定义该坐标系中一个点到原点的距离为这个复数的模长，这条连线与$x$轴正半轴的夹角为这个复数的辐角

观察我们可以发现：复数相乘，模长相乘，辐角相加，证明很简单，代入计算即可

单位根

我们做DFT选择带入的值时，一般都会考虑一些常见好算的整数带入。理论上DFT的时候可以带入任意值，但是带入整数后的IDFT过程不好处理（详细证明见后），而傅里叶想到了带入单位根

我们定义方程$x^n=1$的复数解为$n$次单位根

根据“复数相乘，模长相乘，辐角相加”可以得到所有单位根的模长都是$1$，换句话说就是所有单位根都在单位圆上。因此模长不会对结果造成影响，我们只考虑辐角

我们画出一个单位圆：

容易发现辐角大小为$0,\dfrac{1}{n}\times 2\pi,\dfrac{2}{n}\times 2\pi\dots \dfrac{n-1}{n}\times 2\pi$的复数都是$n$次单位根

由代数基本定理：$n$次方程在复数域内有且只有$n$个根可知以上$n$个复数即为所有单位根

显然所有单位根平分单位圆

我们用$\omega_n^i$来表示逆时针看第$i$个单位根，比如$\omega_n^0$就表示$(1,0)$这个点，所有$n$次单位根为$\omega_n^0,\omega_n^1\dots \omega_n^{n-1}$

当然我们也可以定义$\omega_n^{-k}$为顺时针旋转$\dfrac{k}{n}$个圆周得到的单位根

令$\omega_n^k$的辐角为$\theta$，则该坐标显然为$(\cos\theta,\sin\theta)$

单位根有几个性质：

1.$\omega_n^0=1$对于任意$n$都成立

显然

2.$(\omega_n^1)^k=\omega_n^k$

证明：

根据“复数相乘，模长相乘，辐角相加”可知，每次乘上一个$\omega_n^1$就相当于逆时针转动了$\dfrac{1}{n}$个圆周，转动$k$次就到达了第$k$个单位根

3.$\omega_n^i \times \omega_n^j=\omega_n^{i+j}$

证明：

$\omega_n^i \times \omega_n^j=(\omega_n^1)^i\times (\omega_n^1)^j=(\omega_n^1)^{i+j}=\omega_n^{i+j}$

4.$(\omega_n^i)^j=\omega_n^{ij}$

证明：

$(\omega_n^i)^j=\omega_n^i\times \omega_n^i\dots \times \omega_n^i=\omega_n^{i+i+\dots +i}=\omega_n^{ij}$

5.$\omega_{pn}^{pi}=\omega_n^i$，其中$p$为常数

证明：

由于所有单位根平分单位圆，所以$\omega_{pn}^{pi}$的辐角就是$\dfrac{pi}{pn}\times 2\pi=\dfrac{i}{n}\times 2\pi$，与$\omega_n^i$辐角相同

6.若$n$为偶数，则$\omega_n^{i+n/2}=-\omega_n^i$

证明：

$\omega_n^{i+n/2}=\omega_n^i\times \omega_n^{n/2}$，显然$\omega_n^{n/2}$辐角为$\pi$，乘上它相当于旋转$\pi$，到达关于原点的对称点，此时显然结果为$-\omega_n^i$，不懂可以画图看看

DFT

我们将一个$n-1$次多项式$F(x)$（这里我们让$n=2^k,k\in Z$，若不成立则在高位补$0$）按照奇偶分成两块：

$$F(x)=(f_0\times x^0+f_2\times x^2+\dots +f_{n-2}\times x^{n-2})+(f_1\times x^1+f_3\times x^3+\dots +f_{n-1}\times x^{n-1})$$

我们再定义两个多项式$F_L(x),F_R(x):$

$$F_L(x)=f_0\times x^0+f_2\times x^1+f_4\times x^2+\dots +f_{n-2}\times x^{n/2-1}\\ F_R(x)=f_1\times x^0+f_3\times x^1+f_5\times x^2+\dots +f_{n-1}\times x^{n/2-1}$$

显然有：

$$F(x)=F_L(x^2)+xF_R(x^2)$$

我们将$\omega_n^k(k<\dfrac{n}{2})$带入上式：

$$\begin{equation} \begin{aligned} F(\omega_n^k)&=F_L((\omega_n^k)^2)+\omega_n^kF_R((\omega_n^k)^2)\\ &=F_L(\omega_n^{2k})+\omega_n^kF_R(\omega_n^{2k})\\ &=F_L(\omega_{n/2}^k)+\omega_n^kF_R(\omega_{n/2}^k) \end{aligned} \end{equation}$$

我们再把$\omega_n^{k+n/2}(k<\dfrac{n}{2})$代入上式：

$$\begin{equation} \begin{aligned} F(\omega_n^{k+n/2})& =F_L((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}F_R((\omega_n^{k+n/2})^2)\\ &=F_L(\omega_n^{2k+n})-\omega_n^{k}F_R(\omega_n^{2k+n})\\ &=F_L(\omega_n^{2k})-\omega_n^{k}F_R(\omega_n^{2k})\\ &=F_L(\omega_{n/2}^k)-\omega_n^{k}F_R(\omega_{n/2}^k) \end{aligned} \end{equation}$$

发现什么？这两个式子只有一个符号的差距

所以说如果我们知道了$F_L(x)$和$F_R(x)$在$\omega_{n}^0,\omega_{n}^1\dots \omega_{n}^{n-1}$的点值表示，代入上式即可在$O(n)$时间里求出$F(x)$的点值表示

我们再次观察上述过程，看到第一步：

我们将一个$n-1$次多项式$F(x)$（这里$n=2^k,k\in Z$）按照奇偶分成两块

这个过程是什么？分治

所以我们同样可以对$F_L(x)$和$F_R(x)$进行分治，递归处理，直到只剩下一项

$Code:$

#include <bits/stdc++.h>
using namespace std;
#define MAXN 100005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], tp[MAXN];

void dft(Complex *f, int len) {
    if (len == 1) return;
    Complex *fl = f, *fr = f + (len >> 1);
    int i;

    for (i = 0; i < len; ++i) tp[i] = f[i];
    for (i = 0; i < (len >> 1); ++i) {
        fl[i] = tp[i << 1];
        fr[i] = tp[i << 1 | 1];
    }

    dft(fl, len >> 1);
    dft(fr, len >> 1);

    Complex w(cos(2 * pi / len), sin(2 * pi / len)), buf(1, 0);
    for (i = 0; i < (len >> 1); ++i) {
        tp[i] = fl[i] + buf * fr[i];
        tp[i + (len >> 1)] = fl[i] - buf * fr[i];
        buf = buf * w;
    }

    for (i = 0; i < len; ++i) f[i] = tp[i];
}

IDFT

IDFT就是把点值表示再转化为系数表示，也就是DFT的逆过程

我们令$DFT(F)$得到的点值数列为$G$，那么有：

$$f_k=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}g_i(\omega_n^{-k})^i$$

证明：

显然根据DFT定义有：

$$g_k=\sum\limits_{i=0}^{n-1}f_i(\omega_n^k)^i$$

所以有：

$$\begin{equation} \begin{aligned} \sum\limits_{i=0}^{n-1}g_i(\omega_n^{-k})^i&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}f_j(\omega_n^i)^j(\omega_n^{-k})^i\\ &=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}f_j\omega_n^{ij-ik}\\ &=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}f_j\omega_n^{i(j-k)} \end{aligned} \end{equation}$$

我们分类讨论：

$\bullet$ 若$j=k$

则有：

$$\begin{equation} \begin{aligned} \sum\limits_{i=0}^{n-1}g_i(\omega_n^{-k})^i&=\sum\limits_{i=0}^{n-1}f_k\omega_n^0\\&=nf_k \end{aligned} \end{equation}$$

$\bullet$ 若$j\not=k$

令$p=j-k$

则有：

$$\begin{equation} \begin{aligned} \sum\limits_{i=0}^{n-1}g_i(\omega_n^{-k})^i&=\sum\limits_{i=0}^{n-1}f_{k+p}\omega_n^{ip}\\ &=nf_{k+p}\omega_n^p\sum\limits_{i=0}^{n-1}\omega_n^i \end{aligned} \end{equation}$$

令$S=\sum\limits_{i=0}^{n-1}\omega_n^i$

则有$\omega_n^1S=\sum\limits_{i=1}^{n}\omega_n^i$

所以有：

$$(\omega_n^1-1)S=\omega_n^n-\omega_n^0\\ S=\dfrac{\omega_n^n-\omega_n^0}{\omega_n^1-1}=0$$

所以$(6)$式的值为$0$，对结果不产生贡献

综上可知，$(4)$式的值为$nf_k$，故原等式成立

由此我们再代入$\omega_n^{0}, \omega_n^{-1} \dots \omega_n^{-n+1}$即可完成IDFT

此时我们再反观一开始的问题：为什么要代入单位根？因为单位根的性质使得DFT和IDFT变得十分简单，代入别的数难以达到这个效果

这个IDFT的过程可以被理解为范德蒙德矩阵求逆，这也衍生出单位根反演，在此不多做赘述

$Code:$（和DFT几乎一样）

#include <bits/stdc++.h>
using namespace std;
#define MAXN 100005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], tp[MAXN];

void idft(Complex *f, int len) {
    if (len == 1) return;
    Complex *fl = f, *fr = f + (len >> 1);
    int i;

    for (i = 0; i < len; ++i) tp[i] = f[i];
    for (i = 0; i < (len >> 1); ++i) {
        fl[i] = tp[i << 1];
        fr[i] = tp[i << 1 | 1];
    }

    idft(fl, len >> 1);
    idft(fr, len >> 1);

    Complex w(cos(2 * pi / len), -sin(2 * pi / len)), buf(1, 0);
    for (i = 0; i < (len >> 1); ++i) {
        tp[i] = fl[i] + buf * fr[i];
        tp[i + (len >> 1)] = fl[i] - buf * fr[i];
        buf = buf * w;
    }

    for (i = 0; i < len; ++i) f[i] = tp[i];//除以n的事放到主函数里面做
}

综合在一起，我们就得到了最初版本的FFT：

#include <bits/stdc++.h>
using namespace std;
#define MAXN 100005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], a[MAXN], b[MAXN], tp[MAXN];

void dft(Complex *f, int len) {
    if (len == 1) return;
    Complex *fl = f, *fr = f + (len >> 1);
    int i;

    for (i = 0; i < len; ++i) tp[i] = f[i];
    for (i = 0; i < (len >> 1); ++i) {
        fl[i] = tp[i << 1];
        fr[i] = tp[i << 1 | 1];
    }

    dft(fl, len >> 1);
    dft(fr, len >> 1);

    Complex w(cos(2 * pi / len), sin(2 * pi / len)), buf(1, 0);
    for (i = 0; i < (len >> 1); ++i) {
        tp[i] = fl[i] + buf * fr[i];
        tp[i + (len >> 1)] = fl[i] - buf * fr[i];
        buf = buf * w;
    }

    for (i = 0; i < len; ++i) f[i] = tp[i];
}

void idft(Complex *f, int len) {
    if (len == 1) return;
    Complex *fl = f, *fr = f + (len >> 1);
    int i;

    for (i = 0; i < len; ++i) tp[i] = f[i];
    for (i = 0; i < (len >> 1); ++i) {
        fl[i] = tp[i << 1];
        fr[i] = tp[i << 1 | 1];
    }

    idft(fl, len >> 1);
    idft(fr, len >> 1);

    Complex w(cos(2 * pi / len), -sin(2 * pi / len)), buf(1, 0);
    for (i = 0; i < (len >> 1); ++i) {
        tp[i] = fl[i] + buf * fr[i];
        tp[i + (len >> 1)] = fl[i] - buf * fr[i];
        buf = buf * w;
    }

    for (i = 0; i < len; ++i) f[i] = tp[i];
}

int n, m;
signed main() {
    scanf("%d%d", &n, &m);
    int i;

    for (i = 0; i <= n; ++i) scanf("%lf", &a[i].x);
    for (i = 0; i <= m; ++i) scanf("%lf", &b[i].x);
    for (m += n, n = 1; n <= m; n <<= 1)
        ;

    dft(a, n), dft(b, n);
    for (i = 0; i < n; ++i) f[i] = a[i] * b[i];
    idft(f, n);
    for (i = 0; i <= m; ++i) printf("%d ", (int)(f[i].x / n + 0.5));//在这里除以n
    return 0;
}

由于DFT和IDFT的代码差不多，所以我们可以把他们写到一起：

#include <bits/stdc++.h>
using namespace std;
#define MAXN 5000005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], a[MAXN], b[MAXN], tp[MAXN];

void fft(Complex *f, int len, int flag) {
    if (len == 1) return;
    Complex *fl = f, *fr = f + (len >> 1);
    int i;

    for (i = 0; i < len; ++i) tp[i] = f[i];
    for (i = 0; i < (len >> 1); ++i) {
        fl[i] = tp[i << 1];
        fr[i] = tp[i << 1 | 1];
    }

    fft(fl, len >> 1, flag);
    fft(fr, len >> 1, flag);

    Complex w(cos(2 * pi / len), flag * sin(2 * pi / len)), buf(1, 0);
    for (i = 0; i < (len >> 1); ++i) {
        tp[i] = fl[i] + buf * fr[i];
        tp[i + (len >> 1)] = fl[i] - buf * fr[i];
        buf = buf * w;
    }

    for (i = 0; i < len; ++i) f[i] = tp[i];
}

int n, m;
signed main() {
    scanf("%d%d", &n, &m);
    int i;

    for (i = 0; i <= n; ++i) scanf("%lf", &a[i].x);
    for (i = 0; i <= m; ++i) scanf("%lf", &b[i].x);
    for (m += n, n = 1; n <= m; n <<= 1)
        ;

    fft(a, n, 1), fft(b, n, 1);
    for (i = 0; i < n; ++i) f[i] = a[i] * b[i];
    fft(f, n, -1);
    for (i = 0; i <= m; ++i) printf("%d ", (int)(f[i].x / n + 0.5));
    return 0;
}

以上代码可以通过洛谷P3803 【模板】多项式乘法（FFT）

优化

我们可以观察到上述代码使用递归实现FFT的速度并不尽如人意，所以我们要优化

首先观察分治的过程：

我们发现分治到最后某个数的下标就是原来的下标在二进制意义下翻转的结果，所以我们可以用一个类似数位dp的东西来预处理出每个数在分治序列中的下标，从而避免了递归处理

具体地说，我们令$r_i$为递归后原数列中$f_i$所在的位置。我们把$i$拆成最后一位和前面所有位两部分，显然$r_i$就是最后一位加上前面所有位的翻转

所以有：

r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1))

同时我们观察这部分代码：

for (i = 0; i < (len >> 1); ++i) {
    tp[i] = fl[i] + buf * fr[i];
    tp[i + (len >> 1)] = fl[i] - buf * fr[i];
    buf = buf * w;
}

for (i = 0; i < len; ++i) f[i] = tp[i];

拷贝数组的常数会比较大，所以我们可以一边做一边拷贝：

for (i = 0; i < (len >> 1); ++i) {
    Complex t = buf * fr[i];
    f[i + (len >> 1)] = fl[i] - t;
    f[i] = fl[i] + t;//注意赋值的顺序
    buf = buf * w;
}

同时还可以迭代实现：

void fft(Complex *f, int flag) {
    int i, j, k;
    for (i = 0; i < n; ++i)
        if (i < r[i]) swap(f[i], f[r[i]]);//预处理位置

    for (int len = 2; len <= n; len <<= 1) {//由小到大合并而不是由大到小递归
        Complex w(cos(2 * pi / len), flag * sin(2 * pi / len));
        for (j = 0; j < n; j += len) {
            Complex buf(1, 0);
            for (k = j; k < j + (len >> 1); ++k) {
                Complex t = buf * f[k + (len >> 1)];
                f[k + (len >> 1)] = f[k] - t;
                f[k] = f[k] + t;
                buf = buf * w;
            }
        }
    }
}

完整版（也是最常见的）FFT代码：

#include <bits/stdc++.h>
using namespace std;
#define MAXN 5000005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], a[MAXN], b[MAXN], tp[MAXN];

int n, m;
int r[MAXN];
void fft(Complex *f, int flag) {
    int i, j, k;
    for (i = 0; i < n; ++i)
        if (i < r[i]) swap(f[i], f[r[i]]);

    for (int len = 2; len <= n; len <<= 1) {
        Complex w(cos(2 * pi / len), flag * sin(2 * pi / len));
        for (j = 0; j < n; j += len) {
            Complex buf(1, 0);
            for (k = j; k < j + (len >> 1); ++k) {
                Complex t = buf * f[k + (len >> 1)];
                f[k + (len >> 1)] = f[k] - t;
                f[k] = f[k] + t;
                buf = buf * w;
            }
        }
    }
}

signed main() {
    scanf("%d%d", &n, &m);
    int i, cnt = 0;

    for (i = 0; i <= n; ++i) scanf("%lf", &a[i].x);
    for (i = 0; i <= m; ++i) scanf("%lf", &b[i].x);
    for (m += n, n = 1; n <= m; n <<= 1) ++cnt;

    for (i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (cnt - 1));

    fft(a, 1), fft(b, 1);
    for (i = 0; i < n; ++i) f[i] = a[i] * b[i];
    fft(f, -1);
    for (i = 0; i <= m; ++i) printf("%d ", (int)(f[i].x / n + 0.5));
    return 0;
}

可以看到优化后效果还是很明显的

例题

洛谷P1919 【模板】A*B Problem升级版

给出两个正整数$A,B$，求$A\times B$的值
$1\le A,B\le 10^{1000000}$

我们可以把一个$n$位正整数看成一个$n-1$次的$x=10$的多项式，然后直接乘起来就可以了

#include <bits/stdc++.h>
using namespace std;
#define MAXN 3000005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f[MAXN], a[MAXN], b[MAXN], tp[MAXN];

int n, m;
int r[MAXN];
void fft(Complex *f, int flag) {
    int i, j, k;
    for (i = 0; i < n; ++i)
        if (i < r[i]) swap(f[i], f[r[i]]);

    for (int len = 2; len <= n; len <<= 1) {
        Complex w(cos(2 * pi / len), flag * sin(2 * pi / len));
        for (j = 0; j < n; j += len) {
            Complex buf(1, 0);
            for (k = j; k < j + (len >> 1); ++k) {
                Complex t = buf * f[k + (len >> 1)];
                f[k + (len >> 1)] = f[k] - t;
                f[k] = f[k] + t;
                buf = buf * w;
            }
        }
    }
}

int idx = 0;
int ans[MAXN];
char s1[MAXN], s2[MAXN];
signed main() {
    scanf("%s%s", s1, s2);
    int i, cnt = 0;
    n = strlen(s1), m = strlen(s2);
    --n, --m;

    for (i = n; i >= 0; --i) a[n - i].x = (s1[i] - 48) * 1.0;
    for (i = m; i >= 0; --i) b[m - i].x = (s2[i] - 48) * 1.0;

    for (m += n, n = 1; n <= m; n <<= 1) ++cnt;

    for (i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (cnt - 1));

    fft(a, 1), fft(b, 1);
    for (i = 0; i < n; ++i) f[i] = a[i] * b[i];
    fft(f, -1);

    for (i = 0; i <= n; ++i) {
        ans[i] += (int)(f[i].x / n + 0.5);
        if (ans[i] >= 10) {
            if (i == n) ++n;
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10;
        }
    }

    while (!ans[n] && n) --n;
    ++n;
    while (n--) printf("%d", ans[n]);
    return 0;
}

[ZJOI2014]力

给出$n$个数$\{q_i\}，$令$F_j=\sum\limits_{i<j}\dfrac{q_i\times q_j}{(i-j)^2}-\sum\limits_{i>j}\dfrac{q_i\times q_j}{(i-j)^2},E_i=\dfrac{F_i}{q_i}$，对于每个$1\le i\le n$求出$E_i$
$1\le n\le 10^5$

$$\begin{equation} \begin{aligned} E_j&=\dfrac{F_j}{q_j}\\ &=\dfrac{\sum\limits_{i<j}\dfrac{q_i\times q_j}{(i-j)^2}-\sum\limits_{i>j}\dfrac{q_i\times q_j}{(i-j)^2}}{q_j} \\&=\sum\limits_{i<j}\dfrac{q_i}{(i-j)^2}-\sum\limits_{i>j}\dfrac{q_i}{(i-j)^2}\\ &=\sum\limits_{i<j}\dfrac{q_i}{(j-i)^2}-\sum\limits_{i>j}\dfrac{q_i}{(i-j)^2} \end{aligned} \end{equation}$$

令$f_i=q_i,g_i=\dfrac{1}{i^2}$

$$\begin{equation} \begin{aligned} E_j&=\sum\limits_{i=0}^{j-1}f_ig_{j-i}-\sum\limits_{i=j+1}^{n}f_ig_{i-j} \end{aligned} \end{equation}$$

我们把序列翻转一下，令新序列的$f_i=f^{'}_i$

$$\begin{equation} \begin{aligned} E_j&=\sum\limits_{i=0}^{j-1}f_ig_{j-i}-\sum\limits_{i=0}^{j-1}f^{'}_ig_{j-i} \end{aligned} \end{equation}$$

发现前后两个和式都是卷积的形式，直接FFT即可

#include <bits/stdc++.h>
using namespace std;
#define MAXN 1000005

const double pi = acos(-1);
class Complex {
   public:
    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }
    double x, y;
    Complex operator+(Complex const &B) const {
        return Complex(x + B.x, y + B.y);
    }
    Complex operator-(Complex const &B) const {
        return Complex(x - B.x, y - B.y);
    }
    Complex operator*(Complex const &B) const {
        return Complex(x * B.x - y * B.y, x * B.y + y * B.x);
    }
} f1[MAXN], f2[MAXN], a[MAXN], b[MAXN], tp[MAXN];

int n, m;
int r[MAXN];
void fft(Complex *f, int flag) {
    int i, j, k;
    for (i = 0; i < n; ++i)
        if (i < r[i]) swap(f[i], f[r[i]]);

    for (int len = 2; len <= n; len <<= 1) {
        Complex w(cos(2 * pi / len), flag * sin(2 * pi / len));
        for (j = 0; j < n; j += len) {
            Complex buf(1, 0);
            for (k = j; k < j + (len >> 1); ++k) {
                Complex t = buf * f[k + (len >> 1)];
                f[k + (len >> 1)] = f[k] - t;
                f[k] = f[k] + t;
                buf = buf * w;
            }
        }
    }
}

signed main() {
    scanf("%d", &n);
    int i, cnt = 0;

    for (i = 1; i <= n; ++i) {
        scanf("%lf", &tp[i].x);
        a[i].x = (double)(1.0 / i / i);
        b[i].x = tp[i].x;
    }
    reverse(tp + 1, tp + n + 1);

    int t = n;
    for (n = 1; n < (t << 1); n <<= 1) ++cnt;

    for (i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (cnt - 1));

    fft(a, 1), fft(b, 1), fft(tp, 1);
    for (i = 0; i <= n; ++i) {
        f1[i] = a[i] * b[i];
        f2[i] = a[i] * tp[i];
    }

    fft(f1, -1), fft(f2, -1);
    for (i = 1; i <= t; ++i) {
        printf("%lf\n", (f1[i].x / n) - (f2[t - i + 1].x / n));
    }

    return 0;
}