FFT及NTT

FFT——快速傅里叶变换(附NTT——快速数论变换)

FFT是一种能在O(nlogn)时空复杂度内求解两个函数卷积的优秀算法。

算法思想（DFT）：

对于函数 \(f(x)=\Sigma_{i=0}^{n-1}c_ix^{ik}\) 构造一个函数，即k次离散傅里叶级数 \(h(w^k)\) ，其中 \(w\) 是单位元

\[h(w^k)=\Sigma_{j=0}^{n-1}c_jw^{jk}(w=e^{\frac{2i\pi}{n}}) \]

由于单位元 \(w\) 具有性质

\[\Sigma_{j=0}^{n-1}w^{j}=0 \]

所以我们求卷积时，可以直接将 \(f(x1)\) 的 \(h(x1)\) 与 \(f(x2)\) 的 \(h(x2)\) 相乘得到 \(h(x3)\)，即为 \(f(x1)\) 与 \(f(x2)\) 的卷积 \(f(x3)\) 的k次离散傅里叶级数。

单位元的求法也很简单,可借用实复坐标系帮助求出：

\[w=cos(\frac{\pi}{n})+sin(\frac{\pi}{n})*i; \]

现在我们需要将 \(h(x3)\) 转换回 \(f(x3)\) ，相当于对 \(h(x3)\) 进行一次逆傅里叶变换，类似的有结论：

\[c_k=\frac{1}{n}g(w^{-k})=\frac{1}{n}\Sigma_{j=0}^{n-1}h(w^{j})*w^{-jk} \]

证明：由 \(h(x)\) 的定义得，

\[c_k=\frac{1}{n}\Sigma_{j=0}^{n-1}h(w^{j})*w^{-jk}=\frac{1}{n}\Sigma_{j=0}^{n-1}w^{-jk}\Sigma_{p=0}^{n-1}c_pw^{jp}=\frac{1}{n}\Sigma_{j=0}^{n-1}\Sigma_{p=0}^{n-1}c_pw^{j(p-k)} \]

1、\(p=k\) 时，\(c_kw^{j(k-k)}=c_k\to\Sigma_{j=0}^{n-1}c_kw^{j(k-k)}=nc_k\)

2、\(p\ne k\) 时，\(\Sigma_{j=0}^{n-1}c_pw^{j(p-k)}=0\to\Sigma_{j=0}^{n-1}\Sigma_{p=0}^{n-1}c_pw^{j(p-k)}=0 (p\ne k)\)

综上，\(\Sigma_{j=0}^{n-1}\Sigma_{p=0}^{n-1}c_pw^{j(p-k)}=nc_k\)，证毕。

算法优化（FFT）：

显然，如果我们直接枚举 \(j\) 和 \(p\) ，时间复杂度为 \(O(n^2)\) ，与暴力复杂度相同，且常数更大（那我们为什么要学这个？），所以需要优化。

我们再次仔细观察公式

\[h(w^k)=\Sigma_{j=0}^{n-1}c_jw^{jk}(w=e^{\frac{i\pi}{n}}) \]

我们可以发现

\[h(w^k)=\Sigma_{j=0}^{(n-1)/2}c_{j*2}w^{2*j*k}+\Sigma_{j=0}^{n/2-1}c_{j*2+1}w^{(2*j+1)*k} \]

所以可以构造两个函数

\[he(w^k)=\Sigma_{j=0}^{(n-1)/2}c_{j*2}w^{j*k},ho(w^k)=\Sigma_{j=0}^{n/2-1}c_{2*j+1}w^{j*k}; \]

原式改写为

\[h(w^k)=he(w^{2*k})+w*ho(w^{2*k}) \]

由于 \(he(x^{2*k})\) 与 \(ho(x^{2*k})\) 与 \(h(x^k)\) 同构，所以可以一起求解。

但是我们发现，如果每次都对两边求解，那么总时间复杂度 \(O(n^2)\) ，依然达不到要求。

继续观察

\[h(w^k)=he(w^{2*k})+w*ho(w^{2*k}),h(w^{k+n/2})=he(w^{2*k})+w^{n/2+1}*ho(w^{2*k}); \]

由欧拉公式

\[e^{i\pi}=1 \]

再由单位元定义

\[w^{n/2}=e^{\frac{i\pi}{2}}=-1 \]

所以

\[h(w^k)=he(w^{2*k})+w*ho(w^{2*k}),h(w^{k+n/2})=he(w^{2*k})-w*ho(w^{2*k}); \]

由于 \(h(w^k)\) 与 \(h(w^{k+n/2})\) 可以一起求出，所以可以只对一边求解，时间复杂度降为O(nlogn)

所以有如下递归写法

const double pi=acos(-1.0);//π
struct node//复数
{
	double x,y;//x是实部，y是虚部
	node (double xx=0,double yy=0)
	{
		x=xx;
		y=yy;
	}
};
node operator + (node a,node b) {return node(a.x+b.x,a.y+b.y);}
node operator - (node a,node b) {return node(a.x-b.x,a.y-b.y);}
node operator * (node a,node b) {return node(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}//复数加，减，乘法
void dfs(node *a,int len)
{
    if(len==0)
        return;
   	node a1[(len>>1)+1],a2[(len>>1)+1];
    for(int i=0;i<=len;i+=2)
    {
        a1[i>>1]=a[i];
        a2[i>>1]=a[i+1];
    }
    dfs(a1,len>>1);
    dfs(a2,len>>1);
    node w(cos(pi/len),sin(pi/len)),wn(1,0);
    for(int i=0;i<=(len>>1);i++)
    {
        a[i]=a1[i]+wn*a2[i];
        a[i+(len>>1)]=a1[i]-wn*a2[i];
        wn=wn*w;
    }
}

此时，该算法的时空复杂度已经达到 \(O(nlogn)\) ,已经可以通过大部分题目，但是观察到递归写法多开临时数组且常数较大，遇到某些（毒瘤）数据很强的题，有爆栈或超时的风险，所以继续优化。

考虑观察最后序列序号与原序号。

0	1	2	3	4	5	6	7
0	4	2	6	1	5	3	7

我们再把他们写成二进制

000	001	010	011	100	101	110	111
000	100	010	110	001	101	011	111

不难发现，两组序号互为二进制反转之后的数。所以我们可以直接跳过划分下标奇偶的步骤，直接先将序列排成最后的样式，再用非递归方式计算。

计算二进制反转数的式子：

\[rev[i]=(rev[i>>1]>>1)|((i\&1)\&\&(1<<len)) \]

解释：1、如果 \(x\) 是偶数，那么 \(x\) 可以写成 \((x>>1)<<1\) 的形式，对应到反转数上就是 \(rev[i>>1]>>1\)

​ 2、如果 \(x\) 是奇数，那么 \(x\) 可以写成 \(((x>>1)<<1)|1\) 的形式，对应到反转数上就是 \((rev[i>>1]>>1)|(1<<len)\)

\((len是所有数二进制最高位的下一位的位数)\)

直接贴出代码：

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);//π
struct node
{
	double x,y;
	node (double xx=0,double yy=0)
	{
		x=xx;
		y=yy;
	}
};
node a[4000005],b[4000005];
node operator + (node a,node b) {return node(a.x+b.x,a.y+b.y);}
node operator - (node a,node b) {return node(a.x-b.x,a.y-b.y);}
node operator * (node a,node b) {return node(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
int n,m,limit=1,rev[4000005];
void FFT(node *a,int type)//快速傅里叶变换
{
	for(int i=0;i<limit;i++)
		if(i<rev[i])//防止变回去
			swap(a[i],a[rev[i]]);//将序列变成最终样式
	for(int k=2;k<=limit;k<<=1)
	{
		node wn(cos(pi/(k>>1)),type*sin(pi/(k>>1)));
		for(int i=0;i<limit;i+=k)
		{
			node w(1.0,0.0);
			for(int j=i;j<i+(k>>1);j++)
			{
				node x=a[j],y=a[j+(k>>1)];
				a[j]=x+w*y;
				a[j+(k>>1)]=x-w*y;
				w=w*wn;
			}
		}
	}
}
int main()
{
	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);
	while(limit<=n+m)
		limit<<=1;//最高位的下一位
	for(int i=0;i<limit;i++)
	{
		if((i&1)==0)
			rev[i]=rev[i>>1]>>1;
		else
			rev[i]=(rev[i>>1]>>1)|(limit>>1);
	}
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<limit;i++)
		a[i]=a[i]*b[i];
	FFT(a,-1);//快速傅里叶逆变换
	for(int i=0;i<=n+m;i++)
		printf("%lf ",a[i].x/limit);//千万不要忘了1/n这个系数
} 

最后，我们来证明为什么到最后所有数虚部都变成0

证明：

​ 不需要证明。

​ 考虑两个原函数的系数皆为实数，它们卷积的系数也必为实数。

​ 证毕。

题目模板：

P3803 [模板] 多项式乘法 (FFT)

NTT——FFT的进阶版

考虑在FFT中，我们选用 \(e^{\frac{i\pi}{n}}\) 这个复数来作为单位元，这种方法固然巧妙，但缺点就是容易丢失精度，对于取余类问题也很不好处理，复数运算的常数也较大，所以我们想到在整数中寻找这个 \(w\)

考虑在整数域中，有什么具有与 \(w\) 相同的性质呢？

（快速翻找中。。。。。。）（掏出一本积灰的笔记本。。。。。。）（打开。。。）

显然我们可以想到原根（对于原根的资料请自行查找，本文重点不在于此）

考虑原根的几条性质，都和 \(w\) 要求的完全相同。

所以我们直接用原根替换掉上述代码中的 \(w\) 就OK了。

如果是取余类问题，直接取模数的原根即可；对于非取余类且数字上界较小的问题（如1e9），可以直接用一个比较大的质数当做模数（实际上根本没有取模）（推荐998244353，原根是3）

但是由于单位元的指数是 \((mod-1)/((mod-1)/2^k)\) 有可能有一些(毒瘤) 质数，由于 \(mod-1\) 的 2 的质因子个数不够多，出现指数是小数的情况，就会出错（依然推荐998244353,998244352的2的质因子个数足足23个，支持2e6左右的数据）（差评1e9+7，1e9+6的2的质因子个数只有1，分分钟WA掉）。如果有些 核癌氪氢 的题目，规定这样的数作模数，就会很难办。（不过因为本人是一名蒟蒻，还没遇见过这种题目，评论大佬知道的戳我一下）