快速傅里叶变换与快速数论变换从站在门外到入门

前言

litble不会FFT和NTT，是自己太蒻了。

学习数学知识需要耐心，所以也请和蒟蒻litble一样站在门外的人保持耐心来学，加油吧！
另外，litble很好奇，为什么《数学一本通》讲这些东西只讲了半页纸。

复数

参考博客，同时借了张图。

复数的基本定义

虚数单位:$i=\sqrt{-1}$

对于一个数轴上的实数a，让a×(-1)得到-a，那么a就在数轴上旋转了180°。也就是说，a×i×i相当于把a旋转180°，a×i就是把a旋转90°。于是我们就得到了一个美丽的平面：

复数：一个复数x可以表示为$x=a+bi$，当b=0的时候，这个复数为实数。我们可以把一个复数看作以上平面内的一个点(a,b)。其中，$a$叫做复数的实部，$b$叫做虚部。

共轭复数:两个实部相等，虚部为相反数的复数。

模:即复数的绝对值，用向量来表示复数时的长度，量化一下就是$\sqrt{a^2+b^2}$

复数的三角形式

如图所示，其中$\theta$叫做复数的辐角，幅角可以表示为$\theta +2k\pi$的形式，而小于180°的幅角叫做复数的幅角的主值，两个复数相等当且仅当它们的模和幅角的主值相等。

记r为复数x的模，则$x=r(\cos \theta +i\sin \theta )$叫做复数的三角形式。

复数的运算

复数的运算也满足实数的运算法则。

加减

代数式:$(a+bi)+(c+di)=(a+c)+(b+d)i$

三角式：满足矢量运算法则

乘除

代数式:$(a+bi)(c+di)=ac+bci+adi-bd=(ac-bd)+(bc+ad)i$

三角式:

$r_1(\cos {\theta }_1+i\sin {\theta }_1)r_2(\cos {\theta }_2+i\sin {\theta }_2)=$

$r_1{r}_2(\cos ({\theta }_1+{\theta }_2)+\sin ({\theta }_1+{\theta }_2))$

也就是积的模等于各复数模的积，积的幅角等于各复数幅角的和。

单位根

n次单位根即能够满足方程${\omega }^n=1$的复数

由复数的三角式乘除运算法则可以得到，有n个这样的复数，它们分布于平面的单位圆上，并平分这个单位圆。

n次单位根就是：$$e^{\frac{2\pi ki}{n}},k=0,1,2,...,n-1$$

还有一个很牛的公式，叫做欧拉公式:

$$e^{\theta i}=\cos \theta +i\sin \theta$$

就可以知道n次单位根的算数表示法了。记${\omega }_n=e^{\frac{2\pi i}{n}}$

多项式

参考博客

多项式乘法

给定多项式$A(x)={\sum }_{i=0}^n{a}_i{x}^i$和$B(x)={\sum }_{i=0}^n{b}_i{x}^i$，则它们的积是：
$$C(x)=\sum _{j+k=i,0\le j,k\le n}{a}_j{b}_k{x}^i$$

点值表示法

我们可以取若干个x取值$x_0,x_1,...$，那么多项式$A(x)$的点值表示法就是$A=(x_0,A(x_0)),(x_1,A(x_1))...$

比如说多项式$A(x)=2x^2+x-1$的点值表达式就可以是$A=(0,-1),(1,2),(2,9)$

对于一个次数从0到n-1的多项式，要取n个点才能把它确定下来。

有了点值表示法后，多项式乘法就可以化为$A=(x_0,A(x_0))...$和$B=(x_0,B(x_0))...$的乘积为$C=(x_0,A(x_0)B(x_0))...$了，现在我们只需要迅速地完成点值表示法与代数表示法之间的转化即可。

快速傅里叶变换

参考博客（没错就是数王dalao的博客）

折半引理

对于n大于0且n为奇数

{ ω n 2 k ∣ 0 ≤ k ≤ n , k ∈ Z } = { ω n / 2 k ∣ 0 ≤ k ≤ n / 2 , k ∈ Z } \{ \omega_n^{2k}|0 \leq k \leq n,k \in Z \}= \{ \omega_{n/2}^k|0 \leq k \leq n/2,k \in Z \} {ωn2k​∣0≤k≤n,k∈Z}={ωn/2k​∣0≤k≤n/2,k∈Z}
因为集合要去重嘛。

证明:

$({\omega }_n^{k+\frac{n}{2}}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}{\omega }_n^n=({\omega }_n^k{)}^2={\omega }_{n/2}^k$

如果上式有哪里不懂的话，我猜你八成是忘了${\omega }_n=e^{\frac{2\pi i}{n}}$这件事了。

快速傅里叶变换

即将多项式快速转化为点值表示法。

首先，对原多项式做奇偶性划分：

$$A_0(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{x}{2}-1}$$

$$A_1(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{x}{2}-1}$$

所以$A(x)=A_0(x^2)+x{A}_1(x^2)$

我们把这个变成点集表达式的时候，可以取一些特殊的x来加速，比如，复数${\omega }_n^k$
所以$A({\omega }_n^k)=A_0(({\omega }_n^k{)}^2)+{\omega }_n^k{A}_1(({\omega }_n^k{)}^2)$

由于$(e^{\frac{2\pi i}{n}}{)}^{2k}=(e^{\frac{4\pi i}{n}}{)}^k$

所以$A({\omega }_n^k)=A_0({\omega }_{n/2}^k)+{\omega }_n^k{A}_1({\omega }_{n/2}^k)$
又由折半引理，同时${\omega }_n^{k+\frac{n}{2}}=(e^{\frac{2\pi i}{n}}{)}^{k+\frac{n}{2}}=(e^{\frac{2\pi i}{n}}{)}^k{e}^{\pi i}={\omega }_n^k{\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$

所以得到$A({\omega }_n^{k+n/2})=A_0({\omega }_{n/2}^k)-{\omega }_n^k{A}_1({\omega }_{n/2}^k)$
当k取遍$[0,\frac{n}{2}-1]$时，$k+\frac{n}{2}$就会取遍$[\frac{n}{2},n-1]$

这样我们可以利用分治，迅速将多项式转化成点值表示法了，复杂度为$O(nlo{g}_{2}n)$

逆离散快速傅里叶变换

即将点值表示法快速转化为多项式

这相当于是一个解方程的过程，方程可以写作矩阵(1)，下图摘自树王dalao的博客。

所以呢$e_{ij}={\sum }_{k=0}^{n-1}{d}_{ik}{v}_{kj}={\sum }_{k=1}^n{\omega }_n^{-ik}{\omega }_n^{kj}={\sum }_{k=1}^n{\omega }_n^{k(j-i)}$

然后当i=j时，显然$e_{ij}=n$

当i不等于j时，显然$e_{ij}={\sum }_{k=0}^{n-1}{\omega }_n^{(j-i)k}$，然后根据等比数列随便搞一搞得到：

$$e_{ij}=\frac{({\omega }_n^n{)}^{j-i}-1}{{\omega }_n^{j-i}-1}=0$$

所以$\frac{1}{n}{E}_n=I_n$($I_n$表示n×n的单位矩阵)，也就是$\frac{1}{n}D=V^{-1}$
那么：

而这个D矩阵真是完美，可以和类似于快速傅里叶变换的方法搞定，就大大减小了代码复杂度

蝶形算法

参考博客

例题：洛谷P3803

又名蝴蝶操作，Cooley-Tukey算法（其实就是个迭代）。

这个名字非常优美，这个算法也很优美。

我们手动模拟一个递归过程，每次把偶数位放在前面，奇数位放在后面处理。第一行和第二行表示当前位置的数的二进制形式。那么（下表依旧摘抄自数王dalao）：

$$\begin{array}{ccccccccc}& 000& 001& 010& 011& 100& 101& 110& 111\\ & 0& 1& 2& 3& 4& 5& 6& 7\\ & 0& 2& 4& 6& 1& 3& 5& 7\\ & 0& 4& 2& 6& 1& 5& 3& 7\\ & 0& 4& 2& 6& 1& 5& 3& 7\\ & 000& 100& 010& 110& 001& 101& 011& 111\end{array}$$

发现位置最后的数，是其二进制反过来后得到的数。其实证明也很容易，我们因为按照奇偶划分，会把末尾是1的放后面，末尾是0的放前面。然后把倒数第二位是0的放前面，是1的放后面…算法步骤：

1. 将二进制反序后的序列处理出来。

2. 枚举当前处理的区间长度$2i$（至于为什么要乘2，因为代码里乘了2，此处i不是虚数单位）

3. 根据欧拉公式计算${\omega }_{2i}$

4. 枚举每个区间开头j

5. 枚举该区间的每个偶数位那一半的$a_{j+k}$和奇数位那一半的$a_{i+j+k}$

6. 令$t1=a_{j+k}$，$t2={\omega }_{2i}^k{a}_{i+j+k}$，那么更改$a_{j+k}=t1+t2$，$a_{i+j+k}=t1-t2$。由于这个操作的流程像蝴蝶的翅膀，所以该算法称为蝶形算法。

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define db double
const int N=2650000;//注意空间
const db pi=3.1415926535897384626;
struct com{db r,i;}a[N],b[N];//复数
int n,m,len,rev[N];
com operator + (com A,com B) {return (com){A.r+B.r,A.i+B.i};}
com operator - (com A,com B) {return (com){A.r-B.r,A.i-B.i};}
com operator * (com A,com B) {return (com){A.r*B.r-A.i*B.i,A.i*B.r+A.r*B.i};}
com operator / (com A,db  B) {return (com){A.r/B,A.i/B};}
void FFT(com *a,int x) {
	for(int i=0;i<n;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);//防止两次交换
	for(int i=1;i<n;i<<=1) {//步骤2
		com wn=(com){cos(pi/i),x*sin(pi/i)};//步骤3
		for(int j=0;j<n;j+=(i<<1)) {//枚举区间开头，步骤4
			com w=(com){1,0},t1,t2;
			for(int k=0;k<i;++k,w=w*wn)//步骤5
				t1=a[j+k],t2=w*a[j+k+i],a[j+k]=t1+t2,a[j+k+i]=t1-t2;
		}
	}
	if(x==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;
}
int main() {
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i) scanf("%lf",&a[i].r);
	for(int i=0;i<=m;++i) scanf("%lf",&b[i].r);
	m=n+m; for(n=1;n<=m;n<<=1) ++len;
	for(int i=0;i<n;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));//二进制下的翻转，步骤1
	FFT(a,1),FFT(b,1);
	for(int i=0;i<=n;++i) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<=m;++i) printf("%d ",(int) (a[i].r+0.5));
	return 0;
}

原根

小L学会FFT后，自以为非常牛逼，于是要dalao给她道题水水。dalao随手给了她一道…
卡精度的题。

在快速傅里叶变换里，我们利用${\omega }_n$实现了对于消去引理和折半引理的满足。那么，有没有什么整数也满足呢？

原根就满足。

what is 原根 you may ask

定义P的原根为满足$g^{\varphi (P)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$的整数g，请记住这个定义。

对于一个质数 P,如果它存在原根，那么$g^i̸=g^j$。由于这条性质，所以我们可以利用原根来生成点集表示法。

快速数论变换

参考博客

现在我们用$g^{\frac{\varphi (P)}{n}}$来代替${\omega }_n$进行计算。

对于P的要求：要求$\frac{\varphi (P)}{n}$为整数，又因为当P为质数时，$\varphi (P)=P-1$，所以$P=k{2}^q+1$，其中$2^q\ge n$

原根是否可以代替单位根，取决于FFT的引理，原根是否也满足，以折半引理为例,在模P的意义下：

$$(g^{\frac{\varphi (P)}{n}}{)}^2=g^{\frac{2\varphi (P)}{n}}=g^{\frac{\varphi (P)}{n/2}}$$

$$(g^{\frac{\varphi (P)}{n}}{)}^{k+\frac{n}{2}}=g^{\frac{\varphi (P)}{2}}(g^{\frac{\varphi (P)}{n}}{)}^k=-(g^{\frac{\varphi (P)}{n}}{)}^k$$

所以是满足的，其他的东西推导也类似啦。

那么怎么求原根呢？用bsgs打表。除非题目给你模数，否则可以一劳永逸地选择几个记住：一位dalao的归纳（998244353是个好东西）

如果P不是个质数，就用中国剩余定理玩合并。

但是以上两个过程都比较痛苦就是了…给出NTT模板的代码，例题还是上面那道：

#include<bits/stdc++.h>
using namespace std;
const int G=3,mod=(119<<23)+1,N=2650000;
int n,m,len,a[N],b[N],rev[N];
int ksm(int x,int y) {
	int re=1;
	while(y) {
		if(y&1) re=1LL*re*x%mod;
		x=1LL*x*x%mod,y>>=1;
	}
	return re;
}
void NTT(int *a,int x) {
	for(int i=0;i<n;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1;i<n;i<<=1) {
		int gn=ksm(G,(mod-1)/(i<<1));
		for(int j=0;j<n;j+=(i<<1)) {
			int t1,t2,g=1;
			for(int k=0;k<i;++k,g=1LL*g*gn%mod) {
				t1=a[j+k],t2=1LL*g*a[j+k+i]%mod;
				a[j+k]=(t1+t2)%mod,a[j+k+i]=(t1-t2+mod)%mod;
			}
		}
	}
	if(x==1) return;
	int ny=ksm(n,mod-2); reverse(a+1,a+n);//在FFT里面，我们已经利用x在步骤3中完成了翻转。
	for(int i=0;i<n;++i) a[i]=1LL*a[i]*ny%mod;
}
int main() {
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i) scanf("%d",&a[i]);
	for(int i=0;i<=m;++i) scanf("%d",&b[i]);
	m=n+m;for(n=1;n<=m;n<<=1)++len;
	for(int i=0;i<n;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	NTT(a,1),NTT(b,1);
	for(int i=0;i<n;++i) a[i]=1LL*a[i]*b[i]%mod;
	NTT(a,-1);for(int i=0;i<=m;++i) printf("%d ",a[i]);
	return 0;
}

应用

bzoj4503

将某一元素翻转来构成卷积，是一种常用手段。那么我们将B串翻转，字母全部转化为数字，"?"转化为0之后，我们得到A从位置j开始能和B匹配的要求是：

$$C_{j+m-1}=\sum _{i=0}^{m-1}(a_{j+i}-b_{m-i-1}{)}^2\ast b_{m-1-i}=0$$

展开发现是三个卷积的和，分别用FFT或者NTT搞一搞，查找为0的C值即可。

bzoj4827/洛谷P3723

假设增加的值为C,差异值公式：

$$\sum _{i=0}^{n-1}(a_i+C-b_i{)}^2=\sum _{i=0}^{n-1}{a}_i+\sum _{i=0}^{n-1}{b}_i+n{C}^2-2\sum _{i=0}^{n-1}{a}_i{b}_i+2C\sum _{i=0}^{n-1}(a_i-b_i)$$

假定C为一个常数，那么这个式子有一堆的常数，主要矛盾是求最大的$2{\sum }_{i=0}^{n-1}{a}_i{b}_i$，那么将B手链翻转后用FFT或NTT搞一搞即可。

最后，可以枚举C或者利用二次函数搞定和C有关的最优值。

bzoj3527/洛谷P3338

可得$$E_i=\sum _{j=1}^{i-1}\frac{q_j}{(i-j{)}^2}-\sum _{j=i+1}^n\frac{q_j}{(i-j{)}^2}$$
令$a_i=q_i$，$b_i=\frac{1}{i^2}$，可以发现$\frac{q_j}{(i-j{)}^2}=a_j{b}_{i-j}$
原式就是两个卷积相减，后一个卷积需要反向来构造。

bzoj3771

首先是一个基本的生成函数，比如说有一把斧头的权值为$w_i$，那么就令多项式里面的$x^{w_i}$的权值加1，这样弄出来的一个多项式称为A。

比如选两把斧头，那么就是$A\ast A$，k次项的系数就是方案数，不过还要去除都选同一把斧头的情况。

那么设把斧头权值都乘以2得到的多项式为B，乘以3得到的为C，根据容斥原理，答案就是求：
$$A+\frac{A^2-B}{2}+\frac{A^3-3AB+2C}{6}$$

bzoj4332

首先令t(0)=0，$t(i)=O{i}^2+Si+U$，然后设$g(i,j)$为i个小朋友，每人都有糖的开心程度积（小朋友个数还是设为n吧，显然m小于n的时候n取m即可）。

$$g(i,j)=\sum _{k=0}^{j}g(i-1,k)t(j-k)$$

发现这是一个卷积的形式，可得：$g(i)=t^i$

我们要求的答案是${\sum }_{i=1}^{n}g(i,m)$

于是设$f(n)={\sum }_{i=1}^{n}g(i)$，那么$f(n,m)$就是答案。

考虑分治求解，$f(n)=f(\frac{n}{2})+{\sum }_{i=\frac{n}{2}+1}^{n}g(i)=f(\frac{n}{2})+t^{\frac{n}{2}}f(\frac{n}{2})$

对于奇数的n，可以首先求出$f(n-1)$然后再算。

取模就每次算完多项式乘法后模一下。

复杂度$O(nlo{g}^{2}n)$