快速傅里叶变换FFT/快速数论变换NTT

前言

其实这是3个月前的内容。。。
因为各种原因拖到了现在就是懒而已

FFT

FFT全称快速傅里叶变换（fast Fourier transform），用来解决多项式乘法之类的问题
当然主要还是用来搞卷积
搞♂积
设有两个多项式A(x)和B(x)，要求C(x)=A(x)B(x)
暴力算时间复杂度是$O(n^2)$，用FFT可以优化到$O(n\log n)$

点值与插值

点值

顾名思义，点值就是一个函数上某个点的值
点值对(x,y)就表示在函数上x对应的值y

一个函数的最高项为n，则需要n+1个x不同的点值对就可以把原函数还原出来
那么{(x0,y0),(x2,y2)…(xn,yn)}就是函数的点值表达，把函数转换成点值表达的运算就是点值运算

插值

插值就是点值的逆运算，相当于给出若干个点，要求一条穿过所有点的函数
可以解方程，但一般是用拉格朗日插值法
原来我之前一直用的都是O(2^n)的插值法
拉格朗日♂插值法
暴力的做法就是先把A和B做点值运算，然后把AB的函数值相乘，之后插值运算还原出C
为了保证C的唯一性，A和B都要代入至少C的项数个值

---

点值和插值的复杂度都是$O(n^2)$，过程上并不能优化，但可以选择带入特定的值x

复数

形如$a+bi$的数就是复数，其中i是复数单位，$i=\sqrt{-1}$
a为复数的实部，b为复数的虚部

复平面

复数一般用平面直角坐标系来表示

x轴表示实部，y轴表示虚部
那么一个复数可以用向量来表示，设复数Z的幅角为Z到x正半轴的夹角，Z的模长|z|表示从原点到Z的距离

那么Z可以用三角函数来表示

之后考虑两个复数相乘，就是这样

总之，就是模长相乘，幅角相加

n次单位复数根

就是方程${\omega }^n=1$的解，可以得出一共有n个解
因为结果等于1，所以模长一定为1
所以等于从一个点开始转n次，每次转相同的角度，最后转回初始位置

根据感性理解，n个解一定是从(1,0)开始，把复平面平均分成n份且模长为1的点

这是n=8的情况，其它情况类似

性质

n次单位复数根有一些有用的性质

首先如上图可得，n次单位复数根是n个一循环
也就是${\omega }^i={\omega }^j$ $i\equiv j(modn)$

然后还有一个群的性质：
${\omega }^i̸={\omega }^j$ $imodn̸=jmodn$

消去引理：
${\omega }_{dn}^{di}={\omega }_n^i$ （2|n）

折半引理：
$({\omega }_n^i{)}^2=({\omega }_n^{i+\frac{n}{2}}{)}^2={\omega }_n^{2i}$
因为
$({\omega }_n^i{)}^2={\omega }_n^{2i}$
$({\omega }_n^{i+\frac{n}{2}}{)}^2={\omega }_n^{2i+n}$
因为n次单位复数根的性质，${\omega }_n^{2i}={\omega }_n^{2i+n}$
得证

求和引理：
对于任意正整数n和非负整数k，且n不是k的倍数
则有
${\sum }_{j=0}^{n-1}({\omega }_n^k{)}^j=0$
证明可以用等比数列求和
${\sum }_{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{1-({\omega }_n^k{)}^n}{1-{\omega }_n^k}=\frac{1-({\omega }_n^n{)}^k}{1-{\omega }_n^k}=\frac{1-1^k}{1-{\omega }_n^k}=0$

DFT

DFT全称离散傅里叶变换（Discrete Fourier transform）
定义多项式A(x)的离散傅里叶变换DFT(A)，n为A的项数
则$DFT(A{)}_k=A({\omega }_n^k)$，就是依次代入${\omega }_n^{0k-1}$，求A的点值表达

分治策略

暴力求DFT还是$O(n^2)$，所以要用到n次单位复数根的性质

（为了方便处理，首先把AB都扩大到2的幂大小）
设当前分治到多项式A(x)，$A(x)=a_0+a_{1}x+...+a_{n-1}{x}^{n-1}$
定义两个多项式A0(x),A1(x)，其中
$A0(x)=a_0+a_{2}x+...+a_{n-2}{x}^{\frac{n}{2}-1}$
$A1(x)=a_1+a_{3}x+...+a_{n-1}{x}^{\frac{n}{2}-1}$
然后A(x)就可以表示成
$A(x)=A0(x^2)+xA1(x^2)$
（注意这里指的是每一个不同的x）
所以就变成分治求A0(x)和A1(x)了

因为折半引理，当x取了平方之后，$({\omega }_n^i{)}^2=({\omega }_n^{i+\frac{n}{2}}{)}^2$，所以实际上代入的x就会减少一半
然后就只用代入一半的x就行了

这样每次分成另两个子问题，但每个子问题的规模会减小一半，所以每向下一层，求解的总数一定为n，一共向下$\log n$层
所以时间复杂度为$O(n\log n)$

IDFT

IDFT（inverse discrete Fourier transform）是DFT的逆运算
DFT操作转换成矩阵就是这样：

其中y=DFT(A)
可以看出，$y=V\ast a$，其中V是中间的矩阵
那么$a=y\ast V^{-1}$，就是乘上V的逆矩阵

${V^{-1}}_{i,j}={\omega }_n^{-ji}/n$
证明：
要证明该逆矩阵成立，只需要证出$V\ast V^{-1}=I$，$I$是单位矩阵
$V\ast V^{-1}$的第(i,j)项为
${\sum }_{k=0}^{n-1}{\omega }_n^{ik}{\omega }_n^{-jk}/n$
$=\frac{1}{n}{\sum }_{k=0}^{n-1}{\omega }_n^{ik}{\omega }_n^{-jk}$
$=\frac{1}{n}{\sum }_{k=0}^{n-1}{\omega }_n^{(i-j)k}$
如果i=j，那么结果为1；否则根据求和引理，结果为0
证毕。

所以只需要把${\omega }_n^{ij}$变成${\omega }_n^{-ij}$，最后的答案再/n就行了

非递归FFT

其实我前面都是在扯淡
FFT用递归实现的话，空间/常数巨大且不好写

所以就要用到非递归FFT

数字代表当前位置的x，圆圈代表值
左边就是把0~n-1二进制位反过来
等于是说每次把最后一位为0的放在上，把最后一位为1的放在下，最后二进制就会被翻转

大概看一下这个图也能找出规律吧。。。
每次枚举一个块，之后把这个块从中间分开，做若干次变换

就像这样，做完一层后做下一层
然后就没了
感觉FFT挺简单

FFT最后可能有精度误差，所以要判一下
（C++可以用operator来定义复数）

例题

可以看一下这里

code

裸的多项式乘法
输入n和m表示多项式的最高项数，然后输入两个多项式，求两个多项式的乘积

//Started Date: 13/06/2018 18:23
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
using namespace std;

struct C
{
	double a,b;
	C (double _a=0,double _b=0) {a=_a;b=_b;};
}; //复数
C operator+(C x,C y) {return C(x.a+y.a , x.b+y.b);}
C operator-(C x,C y) {return C(x.a-y.a , x.b-y.b);}
C operator*(C x,C y) {return C(x.a*y.a-x.b*y.b , x.a*y.b+x.b*y.a);}
C a[1000001],b[1000001],A[1000001],B[1000001];

int c[1000001];
int i,j,k,l,len,n,m,N;

void init()
{
	scanf("%d%d",&n,&m);
	fo(i,0,n) scanf("%lf",&a[i].a);
	fo(i,0,m) scanf("%lf",&b[i].a);
	
	len=ceil(log(max(n,m)+1)/log(2))+1;
	N=pow(2,len);
	fo(i,0,N-1)
	{
		k=0;
		j=i;
		fo(l,1,len)
		{
			k=(k<<1)+(j&1);
			j>>=1;
		}
		c[i]=k;
	}
}

void fft(C *a,int type)
{
	int i,j,k,l,s,S;
	C a1,a2;
	
	s=N;S=1;
	fo(i,1,len)
	{
		int S2=S;
		s>>=1;S<<=1;
		
		fo(j,0,S2-1)
		{
			C w(cos(2*M_PI*j/S*type) , sin(2*M_PI*j/S*type));
			
			fo(k,0,s-1)
			{
				l=j+k*S;
				a1=a[l];a2=w*a[l+S2];
				
				a[l]=a1+a2;
				a[l+S2]=a1-a2;
			}
		}
	}
}

int main()
{
	init();
	
	fo(i,0,N-1) A[i]=a[c[i]],B[i]=b[c[i]];
	fo(i,0,N-1) a[i]=A[i],b[i]=B[i];
	
	fft(a,1);
	fft(b,1);
	
	fo(i,0,N-1)
	a[i]=a[i]*b[i];
	fo(i,0,N-1) A[i]=a[c[i]];
	fo(i,0,N-1) a[i]=A[i];
	
	fft(a,-1);
	
	j=N-1;
	while (abs(a[j].a)<0.0000001) j--;
	fo(i,0,j) printf("%0.0lf ",a[i].a/N);printf("\n");
}

NTT

用来搞有模数的情况

实际上不一定要用n次单位复数根来定义${\omega }_n$，只需要满足上面所说的性质
主要是群的性质
${\omega }^i̸={\omega }^j$ $imodn̸=jmodn$

那么在模p的意义下，有一个很好的代替品——原根
设v是p的原根，则满足当0< i< j< p时，$v^i̸=v^j(modp)$
所以可以用$v^{(p-1)\ast i/n}$来代替${\omega }_n^i$（因为一共有p-1个不同的值）

但是要用原根还要满足$(p-1)/n$为整数，所以要用到特殊的质数
即形如a*2^k+1这样的质数（因为n是2的幂），最常用的就是998244353和1004535809，原根都是3
然后就用原根的幂替换n次单位复数根，其它和FFT一样

code

高精度乘法，输入a和b输出a*b

#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
using namespace std;

const int P=998244353;
const int G=3;

long long a[262144],b[262144],A[262144];
long long ans[262144];

char x[262144];
char y[262144];
int c[262144];
int i,j,k,l,len,n,m,N;

long long qpower(long long a,int b)
{
	long long ans=1;
	
	while (b)
	{
		if (b&1) ans=ans*a%P;
		a=a*a%P;
		b>>=1;
	}
	
	return ans;
}

void init()
{
	scanf("%s",x);
	scanf("%s",y);
	n=strlen(x)-1;
	m=strlen(y)-1;
	
	fo(i,0,n) a[i]=x[n-i]-48;
	fo(i,0,m) b[i]=y[m-i]-48;
	
	len=ceil(log(max(n,m)+1)/log(2))+1;
	N=pow(2,len);
	fo(i,0,N-1)
	{
		k=0;
		j=i;
		fo(l,1,len)
		{
			k=(k<<1)+(j&1);
			j>>=1;
		}
		c[i]=k;
	}
}

void ntt(long long a[],int type)
{
	int i,j,k,l,s,S;
	long long a1,a2;
	
	fo(i,0,N-1) A[i]=a[c[i]];
	fo(i,0,N-1) a[i]=A[i];
	
	s=N;S=1;
	fo(i,1,len)
	{
		int S2=S;
		s>>=1;S<<=1;
		
		long long W=qpower(G,((P-1)/S*type+P-1));
		
		fo(j,0,s-1)
		{
			long long w=1;
			
			fo(k,0,S2-1)
			{
				l=j*S+k;
				a1=a[l];a2=w*a[l+S2]%P;
				
				a[l]=(a1+a2)%P;
				a[l+S2]=(a1-a2)%P;
				
				w=w*W%P;
			}
		}
	}
}

int main()
{
	init();
	ntt(a,1);
	ntt(b,1);
	
	fo(i,0,N-1) a[i]=a[i]*b[i];
	
	ntt(a,-1);
	
	j=qpower(N,P-2);
	fo(i,0,N-1) ans[i]=(a[i]*j)%P,ans[i]+=(ans[i]<0)?P:0;
	fo(i,0,N-1)
	if (ans[i]>9)
	{
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
	
	j=N-1;
	while (j && !ans[j]) j--;
	fd(i,j,0) printf("%lld",ans[i]);printf("\n");
}

后记

其实FFT和NTT把板子背下来就行了
任意模数NTT就用中国剩余定理搞搞高级的做法还不会