本文已参与「新人创造礼」活动， 一起开启创造之路。

快速傅里叶改换 (Fast Fourier Transform) ，即利用核算机核算离散傅里叶改换（DFT）的快速核算办法，简称 FFT，于1965年由 J.W.库利和 T.W.图基提出。咱们在这儿只讨论其在算法竞赛中的运用。

问题定义

关于两个多项式 f(x)=∑i=0n−1aixif(x)=\sum_{i=0}^{n-1}a^i\times x^if(x)=∑i=0n−1​aixi 和 g(x)=∑i=0m−1bixig(x)=\sum_{i=0}^{m-1}b^i\times x^ig(x)=∑i=0m−1​bixi，求他们相乘后形成的多项式 $f(x)g(x)$。

问题分析

明显咱们能够以O(n2)O(n^2)O(n2)的复杂度暴力核算上式的每一项的系数。

for (int i=0;i<=n;++i)
    	for (int j=0;j<=m;++j) c[i+j]=a[i]*b[j];

考虑优化。

方才咱们定义了 f(x)=∑i=0n−1aixif(x)=\sum_{i=0}^{n-1}a^i\times x^if(x)=∑i=0n−1​aixi，咱们将其的 $n$ 个系数次序罗列下来，就能够仅有确认这个多项式了。这种表明办法称为多项式的系数表明方式，如 $f(x)$ 能够表明为

通过线性代数的相关常识咱们知道，假如咱们有该多项式所描述的线上的 $n+1$ 个点，也能够仅有确认这个多项式，这种表明办法称为多项式的点值表明方式。所以 $f(x)$ 也能够记为

假如咱们分别在 $f(x)$ 和 $g(x)$ 上确认 $n+m-1$ 个点

则 $f(x)g(x)$ 的点值表明方式如下

假如咱们能够做到快速转化多项式的系数表明和点值表明，咱们就能够用 $O(n)$ 的时刻复杂度核算 $f(x)g(x)$ 的成果了！可是假如咱们给多项式带入 $n+m-1$ 个点去求点值表明，变量的每个取值都需求 $O(n)$ 的核算对应的函数值，总的时刻复杂度是 O(n2)O(n^2)O(n2) 等级，乃至比暴力做常数大得多（

而 FFT 能够完成在 $O(nlogn)$ 等级完成多项式点值表明与系数表明的相互转化。

DFT

DFT 的任务是把一个多项式从系数表明转化为点值表明。

咱们考虑将多项式从系数表明转化为点值表明的时刻瓶颈是什么。明显咱们必须要最终得到 $n$ 个点，而每个自变量去求对应函数值时又不得不进行 $O(n)$ 等级的运算，所以咱们期望对不同自变量求解中产生的进程量能够互相利用。加之即使自变量的值为 2，在指数级运算中数值增大的速度也难以承受。简单发现 $1$ 是个很不错的挑选，为了凑够 $n$ 个点，咱们把视界转向了复平面。

单位根

信任看到这儿的大家都知道什么是复平面吧_(:3

在复平面中的单位圆上的点模长均为 $1$，如下图所示

则单位圆上一切的点能够表明为 $(cos(),isin())$。

$n$ 次单位根是指 n 次幂为 1 的复数。

由欧拉公式

当 $x=$，且 $n,k\in N+$

所以 $n$ 次单位根恰有 $n$ 个且第 $k$ 个是

由上式咱们能够发现，单位根实际上便是把 $(1,0)$ 作为其间一个等分点，把复平面上的单位圆进行 $n$ 等分，然后从原点指向 $n$ 个等分点形成的向量。其间幅角为正且最小的是 e2in=cos(2n)+isin(2n)e^\frac{2i\pi}{n}=cos(\frac{2\pi}{n})+i\times sin(\frac{2\pi}{n})en2i​=cos(n2​)+isin(n2​)，咱们把它记作 n\omega_nn​。

由三角函数的性质和根本数学运算规律能够得到单位根具有以下性质：

1. 共有 $n$ 个不同的 $n$ 次单位根，其间第 $k$ 个便是 nk\omega_n^knk​。
2. 对一切的整数 $p$，都满意 nk=pnpk\omega_n^k=\omega_{pn}^{pk}nk​=pnpk​。
3. nk+n2=−nk\omega_n^{k+\frac{n}{2}}=-\omega_n^{k}nk+2n​​=−nk​
4. nk+n=nk\omega_n^{k+n}=\omega_n^{k}nk+n​=nk​
5. ∑i=0n−1ni=0\sum_{i=0}^{n-1}\omega_n^i=0∑i=0n−1​ni​=0

前四条都比较明显，咱们证明一下第五条性质：

首先令

S=1+x+x2+x3+…+xn−1S=1+x+x^2+x^3+…+x^{n-1}S=1+x+x2+x3+...+xn−1

则

xS=x+x2+x3+…+xn−1+xnxS=x+x^2+x^3+…+x^{n-1}+x^nxS=x+x2+x3+...+xn−1+xn

两式作差能够得到

xS−S=(x+x2+x3+…+xn−1+xn)−(1+x+x2+x3+…+xn−1)=xn−1=(x−1)(1+x+x2+x3+…+xn−1)\begin{aligned} xS-S&=(x+x^2+x^3+…+x^{n-1}+x^n)-(1+x+x^2+x^3+…+x^{n-1} )\\ &=x^n-1\\ &=(x-1)\times(1+x+x^2+x^3+…+x^{n-1}) \end{aligned}xS−S​=(x+x2+x3+...+xn−1+xn)−(1+x+x2+x3+...+xn−1)=xn−1=(x−1)(1+x+x2+x3+...+xn−1)​

由 xn−1=(x−1)(1+x+x2+x3+…+xn−1)x^n-1=(x-1)\times(1+x+x^2+x^3+…+x^{n-1})xn−1=(x−1)(1+x+x2+x3+...+xn−1) 能够得到

• 当 $x=1$ 时，S=n。
• 否则，若 xn=1x^n=1xn=1，则 S=0。

得证 ∑i=0n−1ni=0\sum_{i=0}^{n-1}\omega_n^i=0∑i=0n−1​ni​=0。

那么咱们把单位根作为自变量有什么用途呢？

FFT

关于一个多项式 f(x)=∑i=0n−1aixif(x)=\sum_{i=0}^{n-1}a^i\times x^if(x)=∑i=0n−1​aixi，不失一般性的，咱们能够把 $n$ 视为 2k2^k2k，即整个多项式一共有 2k2^k2k 项，高位补 0。咱们先对其进行一定的转化

令

则

将 $n$ 次单位根带入后能够得到

当 k≥n2k\ge\frac{n}{2}k≥2n​ 时，令 k′=k−n2k’=k-\frac{n}{2}k′=k−2n​，能够得到

假如咱们递归地做下去，咱们只要知道两个多项式 f1(x)f_1(x)f1​(x) 和 f2(x)f_2(x)f2​(x) 在 (n20,n21,…,n2n2−1)(\omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,…,\omega_\frac{n}{2}^{\frac{n}{2}-1})(2n​0​,2n​1​,...,2n​2n​−1​) 处的函数值，就能够求出 $f(x)$ 在 (n0,n1,…,nn−1)(\omega_n^0,\omega_n^1,…,\omega_n^{n-1})(n0​,n1​,...,nn−1​) 处的函数值了。当 $n=1$ 时，1=1\omega_1=11​=1，则函数值便是多项式系数表明的系数。

总体复杂度是 $O(nlogn)$。

IDFT

接下来咱们考虑怎样把多项式从点值表明转化回系数表明。

假如咱们把 DFT 表明为矩阵方式，就能够得到

其间 aia_iai​ 是原多项式 $i$ 次方项的系数，yjy_jyj​ 是自变量 nj\omega_n^jnj​ 对应的函数值。 令

则咱们只需求找到 \Omega 的一个逆矩阵 −1\Omega^{-1}−1，就能够将多项式由点值表明转化为系数表明了。怎样找 \Omega 的逆矩阵呢？

咱们方才证明了 $n$ 次单位根的和为 $0$，遭到这一特性的启示，咱们发现

所以咱们能够得到 \Omega 的一个逆矩阵

所以将一个多项式在分治的进程中乘上的单位根变为其共轭复数，分治完的每一项除以 $n$ 即可得到原多项式的每一项系数。

递归完成

#include <iostream>
#include <cstdio>
#include <cmath>
#include <vector>
using namespace std;
const double Pi=acos(-1.0);
struct asdf {
    long double x,y;
    asdf operator + (const asdf b) const { return (asdf){x+b.x,y+b.y}; }
	asdf operator - (const asdf b) const { return (asdf){x-b.x,y-b.y}; }
	asdf operator * (const asdf b) const { return (asdf){x*b.x-y*b.y,x*b.y+y*b.x}; }
	asdf operator * (const int b) const {return (asdf){x,y*b}; }
};
int N,M,l,tot=1;
void FFT(int len,vector<asdf> &a,int qwq)
{
	if (len==1) return;
	vector<asdf> a1(len>>1),a2(len>>1);
	for (int i=0;i<len;++i)
		if (i%2) a2[i>>1]=a[i];
		else a1[i>>1]=a[i];
	FFT(len>>1,a1,qwq);
	FFT(len>>1,a2,qwq);
	asdf Wn={cos(2.0*Pi/len),qwq*sin(2.0*Pi/len)};
	asdf w={1,0};
	for (int i=0;i<(len>>1);i++)
	{
		a[i]=a1[i]+w*a2[i];
		a[i+(len>>1)]=a1[i]-w*a2[i];
		w=w*Wn;
	}
}
int main()
{
    int N,M;
    scanf("%d%d",&N,&M);
    for (;tot<=N+M;l++) tot<<=1;
    vector<asdf> a(tot+1),b(tot+1);
    for (int i=0;i<=N;i++) scanf("%Lf",&a[i].x);
    for (int i=0;i<=M;i++) scanf("%Lf",&b[i].x);
    FFT(tot,a,1),FFT(tot,b,1);
    for (int i=0;i<=tot;i++) a[i]=a[i]*b[i];
    FFT(tot,a,-1);
    for (int i=0;i<=N+M;i++)
    {
    	printf("%lld ",(long long)(a[i].x/tot+0.5));
    }
    return 0;
}

非递归完成

可是这样核算常数很大，咱们试图用非递归办法完成。

考虑一下奇偶分隔时发生了什么。比方一个长度为 $8$ 的多项式的分治进程为

假如咱们能够直接将数组变成最终一步的样子，就能够自下而上递推核算了。

观察上述进程能够发现，第一次奇偶分列时，下标最低位是 $0$ 的系数被咱们分去了前半段。第2次整个序列长度减小了一半，能够视为将一切元素的下标整除 $2$ 后再将下标最低位是 $0$ 的系数分给前半段。以此类推。

所以系数 aia_iai​ 最终到达的方位是 $i$ 的二进制位逆序后的成果。这个东西能够递推预处理。令 $r[i]$ 表明改换后 $i$ 应该在的方位，多项式的长度为 $tot$，即多项式系数的下标为 $[0,tot-1]$，不失一般性将 $tot$ 视为 $2$ 的幂且令 n=log2(tot)n=log_2(tot)n=log2​(tot)，分析如下：

假定咱们当前枚举到了 i=∑j=0n−1sj2ji=\sum_{j=0}^{n-1}s_j\times 2^ji=∑j=0n−1​sj​2j，其二进制表明为

i=(sn−1,sn−2,…,s1,s0)2i=(s_{n-1},s_{n-2},…,s_1,s_0)_2i=(sn−1​,sn−2​,...,s1​,s0​)2​

将其右移一位得到

i>>1=(0,sn−1,sn−2,…,s1)2i>>1=(0,s_{n-1},s_{n-2},…,s_1)_2i>>1=(0,sn−1​,sn−2​,...,s1​)2​

那么咱们已经求得了将 $i>>1$ 翻转后的成果

r[i>>1]=(s1,s2,…,sn−2,sn−1,0)2r[i>>1]=(s_1,s_2,…,s_{n-2},s_{n-1},0)_2r[i>>1]=(s1​,s2​,...,sn−2​,sn−1​,0)2​

所以咱们只需求将 $r[i>>1]$ 右移一位再或上 s0<<(n−1)s_{0}<<(n-1)s0​<<(n−1)，就能够得到

r[i]=(s0,s1,s2,…,sn−2,sn−1)2r[i]=(s_0,s_1,s_2,…,s_{n-2},s_{n-1})_2r[i]=(s0​,s1​,s2​,...,sn−2​,sn−1​)2​

综上，搬运方程为：

下面给出非递归完成的代码。

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double Pi=acos(-1.0);
struct asdf
{
    double x,y;
}a[10000010],b[10000010];
asdf operator + (asdf a,asdf b)
{
	return (asdf){a.x+b.x,a.y+b.y};
}
asdf operator - (asdf a,asdf b)
{ 
	return (asdf){a.x-b.x,a.y-b.y};
}
asdf operator * (asdf a,asdf b)
{ 
	return (asdf){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
int N,M,l,r[10000010],tot=1;
void FFT(asdf *A,int qwq)
{
    for (int i=0;i<tot;i++) 
        if (i<r[i]) swap(A[i],A[r[i]]);
    for (int mid=1;mid<tot;mid<<=1)
	{
        asdf Wn=(asdf){cos(Pi/mid),qwq*sin(Pi/mid)};
        for (int R=mid<<1,j=0;j<tot;j+=R)
        {
            asdf w=(asdf){1,0};
            for (int k=0;k<mid;k++,w=w*Wn) 
            {
                asdf x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y,A[j+mid+k]=x-y;
            }
        }
    }
}
int main()
{
    int N,M;
    scanf("%d%d",&N,&M);
    for (int i=0;i<=N;i++) scanf("%lf",&a[i].x);
    for (int i=0;i<=M;i++) scanf("%lf",&b[i].x);
    for (;tot<=N+M;l++) tot<<=1;
    for (int i=0;i<tot;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1),FFT(b,1);
    for (int i=0;i<=tot;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for (int i=0;i<=N+M;i++) printf("%d ",(int)(a[i].x/tot+0.5));
    return 0;
}

预处理优化精度

队友用我之前写的上面那个板子在代码源提交大 WA 特 WA，发现精度丢失首要出现在咱们求第 $k$ 个 $n$ 次单位根时使用了乘法。

假如咱们将其改成每次使用单位根时用用公式直接核算，需求屡次核算三角函数值，很简单 TLE。

咱们能够直接预估程序需求处理的项数最多的多项式有多少项，并将其补成 $2$ 的若干次幂，记为 $FFTN$。直接预处理好 $FFTN$ 次单位根，每次使用时就能够通过性质

对一切的整数 $p$，都满意 nk=pnpk\omega_n^k=\omega_{pn}^{pk}nk​=pnpk​

快速且精度较高地求出咱们需求的单位根。

代码如下：

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double Pi=acos(-1.0);
const int FFTN=262144; 
struct asdf {
    long double x,y;
    asdf operator + (const asdf b) const { return (asdf){x+b.x,y+b.y}; }
	asdf operator - (const asdf b) const { return (asdf){x-b.x,y-b.y}; }
	asdf operator * (const asdf b) const { return (asdf){x*b.x-y*b.y,x*b.y+y*b.x}; }
	asdf operator * (const int b) const {return (asdf){x,y*b}; }
}a[10000010],b[10000010],nw[10000010];
int N,M,l,r[10000010],tot=1;
void FFT(asdf *A,int qwq)
{
    for (int i=0;i<tot;i++) 
        if (i<r[i]) swap(A[i],A[r[i]]);
    for (int mid=1;mid<tot;mid<<=1)
	{
        for (int R=mid<<1,j=0;j<tot;j+=R)
        {
        	int t=FFTN/2/mid;
            for (int k=0;k<mid;k++) 
            {
                asdf x=A[j+k],y=nw[t*k]*qwq*A[j+mid+k];
                A[j+k]=x+y,A[j+mid+k]=x-y;
            }
        }
    }
}
int main()
{
	for (int i=0;i<FFTN;++i)
		nw[i]=(asdf){cosl(2*Pi*i/FFTN),sinl(2*Pi*i/FFTN)};
    int N,M;
    scanf("%d%d",&N,&M);
    for (int i=0;i<=N;i++) scanf("%Lf",&a[i].x);
    for (int i=0;i<=M;i++) scanf("%Lf",&b[i].x);
    for (;tot<=N+M;l++) tot<<=1;
    for (int i=0;i<tot;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1),FFT(b,1);
    for (int i=0;i<=tot;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for (int i=0;i<=N+M;i++) printf("%lld ",(long long)(a[i].x/tot+0.5));
    return 0;
}

其他

虽然上面的代码已经比较优秀了，可是咱们还能够对常数进行进一步的优化，日后遇到卡常的题目我再来进行弥补。