P3803 【模板】多项式乘法（FFT）

---

$更好的阅读体验$

原题链接

P3811
AC记录：Accepted

题目大意

给定一个 $n$ 次多项式 $F(x)$，和一个 $m$ 次多项式 $G(x)$。
请求出 $F(x)$ 和 $G(x)$ 的卷积。

输入格式

第一行两个整数 $n,m$。
接下来一行 $n+1$ 个数字，从低到高表示 $F(x)$ 的系数。
接下来一行 $m+1$ 个数字，从低到高表示 $G(x)$ 的系数。

输出格式

一行 $n+m+1$ 个数字，从低到高表示 $F(x)\cdot G(x)$ 的系数。

$\mathbf{S}\mathbf{a}\mathbf{m}\mathbf{p}\mathbf{l}\mathbf{e}$ $\mathbf{I}\mathbf{n}\mathbf{p}\mathbf{u}\mathbf{t}$

1 2
1 2
1 2 1

$\mathbf{S}\mathbf{a}\mathbf{m}\mathbf{p}\mathbf{l}\mathbf{e}$ $\mathbf{O}\mathbf{u}\mathbf{t}\mathbf{p}\mathbf{u}\mathbf{t}$

1 4 5 2

$\mathbf{H}\mathbf{i}\mathbf{n}\mathbf{t}\&\mathbf{E}\mathbf{x}\mathbf{p}\mathbf{l}\mathbf{a}\mathbf{i}\mathbf{n}$
样例如图所示。
其中粉色虚线为 $F(x)$，绿色虚线为 $G(x)$，蓝色实线为 $F(x)\cdot G(x)$。

数据范围

保证输入中的系数大于等于 $0$ 且小于等于 $9$。
对于 $100\%$ 的数据，$1\le n,m\le {10}^6$。

解题思路

$\text{温馨提示}\text{此题解的字数有亿点点多，思维难度也会偏难，请读者}\text{做好准备。}$

---

题目都告诉你了，是用FFT啊。

在介绍FFT之前，你需要知道一些前置的知识。

---

前置知识

$\text{1.}$复数
复数由两部分组成：实部和虚部。设 $a,b$ 为实数，$i=\sqrt{-1}$，则形如 $a+bi$ 的数叫做复数。

$\text{2.}$复数的运算法则
设有两个复数为 $a+bi$ 和 $c+di$。

加法：
$$\begin{array}{c}(a+bi)+(c+di)\\ =a+bi+c+di\\ =(a+c)+(b+d)i\end{array}$$
其实就是把实部和虚部分别相加就可以了。

乘法：
$$\begin{array}{c}(a+bi)(c+di)\\ =ac+adi+bci+bd{i}^2\\ =ac+(ad+bc)i-bds\\ =(ac-bd)+(ad+bc)i\end{array}$$

几何定义：复数相乘，模长相乘，幅角相加。

$\text{Attention : }$以下内容默认 $n$ 为 $2$ 的正整数次幂。

$\text{3.}$单位根
在复平面（$x$ 轴代表实部，$y$ 轴代表虚部）上画一个单位圆（半径为 $1$ 的圆），以圆点为起点，圆的 $n$ 等分点为终点，做 $n$ 个向量，设幅角为正且最小的向量对应的复数为 ${\omega }_n$，称为 $n$ 次单位根。下图所示的是 $3$ 次单位根。

根据复数乘法的运算法则，不难推出剩下的 $n-1$ 个单位根为 ${\omega }_n^2,{\omega }_n^3,\cdots ,{\omega }_n^n$。
注意：${\omega }_n^0={\omega }_n^n$。

幅角：假设以逆时针为正方向，从$x$轴正半轴到已知向量的转角的有向角叫做幅角

$\text{4.}$单位根的性质

• 由欧拉公式可得，${\omega }_n^k=\cos k\frac{2\pi }{n}+i\sin k\frac{2\pi }{n}$。
• $\forall a,{\omega }_{an}^{ak}={\omega }_n^k$
  证明：
  $$\begin{array}{c}{\omega }_{an}^{ak}=\cos ak\cdot \frac{2\pi }{an}+i\sin ak\frac{2\pi }{an}\\ =\cos k\cdot \frac{2\pi }{n}+i\sin k\frac{2\pi }{n}\\ ={\omega }_n\end{array}$$
• ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
  证明：
  $$\begin{array}{c}{\omega }_n^{k+\frac{n}{2}}={\omega }_n^k\cdot {\omega }_n^{\frac{n}{2}}\\ ={\omega }_n^k\cdot \left(\cos \frac{n}{2}\cdot \frac{2\pi }{n}+i\sin \frac{n}{2}\cdot \frac{2\pi }{n}\right)\\ ={\omega }_n^k\cdot \left(\cos \pi +i\sin \pi \right)\\ ={\omega }_n^k\cdot (-1)\\ =-{\omega }_n^k\end{array}$$

接下来，就是重头戏FFT了。

---

快速傅里叶变换（FFT）

设有一个多项式 $A(x)$，他的系数为 $(a_0,a_1,a_2,\cdots ,a_{n-1})$。
那么
$$A(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1}$$
按照 $x$ 的次数的奇偶性分类，有
$$A(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1})$$
设
$$\begin{gathered} A_1(x)=a_0+a_{2}x+\cdots +a_{n-2}{x}^{\frac{n}{2}-1} \\ {A}_2(x)=a_1+a_{3}x+\cdots +a_{n-1}{x}^{\frac{n}{2}-1} \end{gathered}$$
不难得到
$$A(x)=A_1(x^2)+x{A}_2(x^2)$$
接下来就是把单位根代入了。
但是，不要着急，先把前半部分代进去。
将 ${\omega }_n^k(k<\frac{n}{2})$ 代入原式，有
$$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})$$
再把后半部分代进去。
将 ${\omega }_n^{k+\frac{n}{2}}(k<\frac{n}{2})$ 代入原式，有
$$\begin{gathered} A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2\left(k+\frac{n}{2}\right)})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2(k+\frac{n}{2})}) \\ =A_1({\omega }_n^{2k+n})-{\omega }_n^k{A}_2({\omega }_n^{2k+n}) \\ =A_1({\omega }_n^{2k}\cdot {\omega }_n^n)-{\omega }_n^k{A}_2({\omega }_n^{2k}\cdot {\omega }_n^n) \\ =A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k}) \end{gathered}$$
发现没有，这两个式子只有一个常数项不同！
由于 $k$ 在取遍 $\left[0,\frac{n}{2}-1\right]$ 的时候，$k+\frac{n}{2}$ 取遍了 $\left[\frac{n}{2},n-1\right]$。所以，我们在求第一部分的时候，也可以同时求出第二部分的值。
直接把问题缩小了一半！
直接递归求解就好了。
时间复杂度：$\Theta (n{\log }_{2}n)$。

---

快速傅里叶逆变换（IFFT）

不要以为FFT就结束了。
刚才FFT是基于 点值表示法 的。
但题目里面说的是要 系数表示法 。所以要把点值表示法转换成系数表示法，这个过程就叫 傅里叶逆变换 。
设 $(y_0,y_1,\cdots ,y_{n-1})$ 为 $(a_0,a_1,\cdots ,a_{n-1})$ 的傅里叶变换（就是点值表示法）。
又设 $(c_0,c_1,\cdots ,c_{n-1})$ 为 $(y_0,y_1,\cdots ,y_{n-1})$ 在 $({\omega }_n^0,{\omega }_n^{-1},\cdots ,{\omega }_n^{-(n-1)})$ 的点值表示，即
$$\begin{gathered} c_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i \\ =\sum _{i=0}^{n-1}{a}_j(\sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i) \end{gathered}$$
先到这里，让我们换一个地方。
设 $S(x)=\sum _{i=0}^{n-1}{x}^i$。
当 $k\ne 0$ 时，把 ${\omega }_n^k$ 代入得
$$S({\omega }_n^k)=1+{\omega }_n^k+({\omega }_n^k{)}^2+\cdots +({\omega }_n^k{)}^{n-1}(1)$$
$(1)\times {\omega }_n^k$，得
$${\omega }_n^{k}S({\omega }_n^k)={\omega }_n^k+({\omega }_n^k{)}^2+\cdots +({\omega }_n^k{)}^n(2)$$
$(2)-(1)$，得
$$\begin{gathered} {\omega }_n^{k}S({\omega }_n^k)-S({\omega }_n^k)=({\omega }_n^k{)}^n-1 \\ ({\omega }_n^k-1)S({\omega }_n^k)=({\omega }_n^n{)}^k-1 \\ S({\omega }_n^k)=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1} \\ S({\omega }_n^k)=\frac{1-1}{{\omega }_n^k-1} \\ S({\omega }_n^k)=0 \end{gathered}$$
所以，当 $k\ne 0$ 时，$S({\omega }_n^k)=0$。
那当 $k=0$ 时呢？
很显然，当 $k=0$ 时，$S({\omega }_n^k)=n$。因为此时 ${\omega }_n^k=1$，而又因为 $1^n=1,n\in \mathbb{Z}$，而这样的项共有 $n$ 项，所以整个式子的值为 $n\times 1=n$。

回到原式，继续考虑
$$c_k=\sum _{i=0}^{n-1}{a}_j(\sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i)$$
由上面的式子可以得出，当 $j\ne k$ 时，$\sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i=0$。当 $j=k$ 时，$\sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i=n$。因此，
$$\begin{gathered} c_k=n{a}_k \\ {a}_k=\frac{1}{n}{c}_k \end{gathered}$$
由于 $(c_0,c_1,\cdots ,c_{n-1})$ 为 $(y_0,y_1,\cdots ,y_{n-1})$ 在 $({\omega }_n^0,{\omega }_n^{-1},\cdots ,{\omega }_n^{-(n-1)})$ 的点值表示，所以相当于又做了一次FFT。

---

迭代实现

尽管你用了上面这么多字之精华写出来了一串代码，但是，你交到洛谷上面去，你还是会$\text{WA }\text{77pts}$。这时候就需要迭代来实现了。
盗用一下某位大佬的图

发现了吗？我们实际上依此操作的元素，实际上就是原来序列的下标的二进制反转得到的！
所以现在的难题就是：如何得到操作后的序列呢？
我们像，对于某个数 $x$ 来说，他的原序列，是由 $\frac{x}{2}$ 左移一位再加上最后一位的特判得到的。那么在操作后的序列中也一样，他是由 $\frac{x}{2}$ 右移一位再加上最后一位（也就是第一位）的特判得到的。因此，我们就可以在 $\Theta (n)$ 的时间内求出 $[1,n]$ 操作后的序列所对应的数了。这里给出转移公式：
设数 $x$ 在 $k$ 位二进制下操作后的序列所对应的数为 $r_x$，有
$$r_x=\lfloor \frac{r_{\frac{x}{2}}}{2}\rfloor +[x\mathrm{m}\mathrm{o}\mathrm{d}2=1]\cdot 2^{k-1}$$
其中 $[x]$ 代表当 $x$ 成立时 $[x]=1$，不成立时 $[x]=0$。

上代码

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

struct Complex{
    Complex(double a=0.00,double b=0.00):real(a),imag(b){}
    double real,imag;
};

const double    pi=acos(-1.00);
Complex         a[4000010],b[4000010];
int             resort[4000010];
int             n,m,lim,dig;

Complex operator + (Complex a,Complex b)
{
    return Complex(a.real+b.real,a.imag+b.imag);
}

Complex operator - (Complex a,Complex b) {
    return Complex(a.real-b.real,a.imag-b.imag);
}

Complex operator * (Complex a,Complex b)
{
    double real,imag;
    real=a.real*b.real-a.imag*b.imag;
    imag=a.real*b.imag+a.imag*b.real;
    return Complex(real,imag);
}

void FFT(Complex *c,int state)
{
    for(int i=0; i<lim; i++)
        if(i<resort[i])
            swap(c[i],c[resort[i]]);
    for(int i=1; i<lim; i<<=1)
    {
        Complex W1n(cos(pi/i),state*sin(pi/i));
        for(int Size=i<<1,j=0; j<lim; j+=Size)
        {
            Complex W(1.00,0.00);
            for(int k=0; k<i; k++,W=W*W1n)
            {
                Complex x=c[j+k],y=W*c[j+i+k];
                c[j+k]=x+y;
                c[j+i+k]=x-y;
            }
        }
    }
    return;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    /* Code */
    cin>>n>>m;
    for(int i=0; i<=n; i++)
        cin>>a[i].real;
    for(int i=0; i<=m; i++)
        cin>>b[i].real;
    lim=1,dig=0;
    while(lim<=n+m)
        lim<<=1,dig++;
    for(int i=0; i<lim; i++)
        resort[i]=(resort[i>>1]>>1)|((i&1)<<(dig-1));
    FFT(a,1);
    FFT(b,1);
    for(int i=0; i<lim; i++)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0; i<=n+m; i++)
        cout<<(int)(a[i].real/lim+0.5)<<" ";
    return 0;
}

完美切题$\sim$