FFT&NTT(以及扩展)

预备知识:用于NTT

NTT/FFT其实本质相同,用途是快速求解 多项式乘积

前言

FT: 傅里叶变换:

这是一个工程上的概念,可以简述为:一个周期性的信号波段可以用 若干个正弦曲线 的带权和表示

DFT: 离散傅里叶变换,这是傅里叶变换在离散情况下的变种

FFT: 快速傅里叶变换

NTT: 快速数论变换

 

谈及核心思想

1.单位根:

构造ωnn阶单位根(不知道ωn的值域),满足性质ωnn=ωn0=1

对于2|n,ωnn2=1

显然ωn满足一个非常简单的性质:折半引理 2|i\and2|n,ωni=ωn2i2

ωn实际上是一个在幂次上呈现n元循环的数值

2.多项式点值式的转化

一个n阶多项式最普通的表示就是F(x)=i=0n1aixi

然而,多项式也可以用n互不相关的点表示,即(x0,y0),(x1,y1),,(xn1,yn1)

两者可以互相转化

对于同xi的点值,两个多项式卷积时,其yi可以直接对应相乘

FFT/NTT的核心过程是

多项式 点值式点值式对应相乘多项式

而用单位根来构造快速的多项式与点值式的转化

3.分治思想

用于降低多项式与点值式转换的复杂度

 

FFT的单位根

(x,y)指复数i=1,(x,y)=x+yi

基本运算(x,y)+(a,b)=(x+a,y+b),(x,y)(a,b)=(axby,ay+bx)

FFT的单位根是:ωn=(cos(2πn),sin(2πn))

ωni=(cos(2πni),sin(2πni)) (展开发现就是三角函数求和公式)

显然满足单位根的性质​

(实际上可以发现,这个说是点值其实就是信号序列的三角函数表示)

 

NTT

相信您已经了解了原根的一些性质,NTT的单位根常用原根构造

NTT的单位根实际有较大的局限性,对于质数P只能构造出n|P1,ωn=gP1n

计算在模意义下就能满足单位根的性质

通常我们P998244353223|(P1),它的一个原根是3

实际上,为了满足下面分治需要,构造的模数通常满足P1=s2tt较大,这类模数我们常称作NTT模数

  

以上部分均为基础知识,相对来说应该不会太难,下面是主要难点

 

多项式转点值式

接下来我们考虑如何将多项式转化为点值式

对于点值式,我们构造的点横坐标为xi=ωni

具体目标是对于函数F(x),求出在x0,x1,,xn1上的函数值

即求出F(xi)=a0ωn0+a1ωni+a2ωn2i+

接下来就是核心的分治思想,注意,这里的分治是子问题严格等大

对于当前问题,分成两部分子问题求解(实际是可以分成多部分的,但是这个是特殊情况暂时不予讨论),即求解

m=n2

2|i,G(xi)=a0ωm0+a2ωmi2+a4ωmi22+ 2|i,H(xi)=a1ωm0+a3ωmi2+a5ωmi22+

更简洁的描述为

i<m,G(xi)=a0ωm0+a2ωmi+a4ωm2i+

i<m,H(xi)=a1ωm0+a3ωmi+a5ωm2i+

由于G(xi),H(xi)计算的是[0,m1]项,而求F(xi)时用到的是0,2,4,项,实际需要访问G(xi2),H(xi2)

F(xi)的式子比较,我们得到合并的式子为

F(xi)=G(xi2)+xiH(xi2)

带入折半引理,实际等价于

F(xi)=G(xi)+xiH(xi)

注意xi=ximodm

为了保证复杂度,尽量使得每次分治的子问题都分为两部分,这样的复杂度为O(nlogn)

附:实际上,分为d个子问题时,每次合并的复杂度为O(nd),因此复杂度为

保证每次分治为两个严格等大的子问题,可以从一开始就把n扩充为2的幂次

1
2
int N=1;
while(N<=n+m) N<<=1;

附:d个子问题时,设子问题答案为Gj(xi),则合并的式子为

F(xi)=j=0d1xijGj(xid)=i=0d1xijGj(ximodnd)

点值式转多项式

一个简单的性质:单位根反演 j=0n1ωnij={ωnin1ωni1=0i0ni=0

设点值式对应yi的序列为bi

nai=j=0n1ωnijbj

证明

j=0n1ωnijbj=j=0n1ωnij(k=0n1akωnjk) j=0n1ωnijbj=k=0n1akj=0n1ωnj(ki)

由上面的式子,发现只有ki=0时右边的求和式有值,所以上式成立

因此点值式转多项式直接把系数改为ωni即可

  

Tips:

1.由于单位根的循环特性,溢出会直接溢出到本来的式子里

因此,如果乘法过后的多项式产生了超过>n的项xi,会溢出到ximodn

2.点值式并不是不满足除法,只是除法得到的多项式并不一定是一个n元以内的多项式,除了恰好整除的情况,得到的通常是一个无穷级数的式子,如11x=1x1x=i=0xi

真正要求除法,通常是求前n项的结果,即需要用到多项式乘法逆

 

代码实现与优化

模板题传送门

然后我们得到一份优美的代码(FFT)

(Complex是C++库自带的复数,M_PI是C++自带π常量)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void FFT(int n,Complex *a,int f) {
if(n==1) return;
Complex tmp[N];
int m=n/2;
rep(i,0,m-1) tmp[i]=a[i<<1],tmp[i+m]=a[i<<1|1]; // 按照奇偶分类
memcpy(a,tmp,sizeof(Complex) * n);
FFT(m,a,f),FFT(m,a+m,f); // 分两半,算g(x),h(x)
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0); // w=x^1,e=x^i
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m]; // f(x_i)=g(x_i)+e*h(x_i)
e=e*w;
}
rep(i,m,n-1) {
tmp[i]=a[i-m]+e*a[i];
e=e*w;
}
memcpy(a,tmp,sizeof(T)*n);
}

由于(ωn)n2=1,所以还可以简化为

1
2
3
4
5
6
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0);
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m];
tmp[i+m]=a[i]-e*a[i+m];
e=e*w;
}

由于用了double,最后输出要取整

蝴蝶优化

我们加一点优化,取代递归的分治过程

可以看到,分治时我们按照imod2分成两组,然后继续分

这个过程中,实际上我们就是将i的二进制位前后翻转

所以我们可以暴力处理出i分治底层的位置

1
2
3
4
5
6
7
rep(i,0,n-1) {
int x=i,s=0;
for(int j=1;(j<<c)<=n;++j) {
s=(s<<1)|(x&1);
x>>=1;
} // s就是最终位置
}

当然也是有O(n)处理方法的

1
2
3
int N=1,c=-1;
while(N<=n+m) N<<=1,c++;
rep(i,1,N-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);

(建议自己模拟一下)

有了这个翻转数组,我们可以直接从分治底层开始解决整个问题,每次合并操作完全相同

每次分治问题的大小,依次合并每一个子问题区间即可

为了在一个数组上完成操作,还需要注意合并顺序

代码解释i:分治子问题大小为2il:合并区间的左端点为l,右端点为l+2ij枚举合并位置

1
2
3
4
5
6
7
8
9
10
11
12
13
14
void FFT(int n,Complex *a){
for(int i=1;i<n;i<<=1) {
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n));
for(int l=0;l<n;l+=i*2) {
Complex e(1,0);
for(int j=l;j<l+i;++j,e=e*w) {
Complex t=a[j+m]*e; // a'[j]=a[j]+e*a[j+m]
// a'[j+i]=a[j]-e*a[j+m]
a[j+m]=a[j]-t;
a[j]=a[j]+t;
}
}
}
}

事实上我们还有更快的写法,就是将ωni预处理出来

(注意这个FFT的预处理很考验double精度,不能每次都直接累乘上去,隔几个就要重新调用依次三角函数)

当然如果自己写复数会更快

 

关于点值式转多项式的优化

由于每次求得点值是ωni=ωnni

所以可以直接用 多项式转点值式的函数, 最后把[1,n1]这一段翻转,每个数除掉n即可

 

对于加减运算取模的优化

三目运算

1
2
a+=b,a=a>=P?a-P:a;
a-=b,a=a<0?a+P:a;

逻辑运算优化(原理是逻辑预算会在第一个确定表达式值的位置停下)

1
2
a+=b,((a>=P)&&(a-=P));
a-=b,((a<0)&&(a+=P));
 

关于系数预处理优化(以NTT为例)

带入上面已经提到的优化,无预处理系数的NTT大概是这样的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int l=0;l<n;l+=i*2) {
int e=1;
for(int j=l;j<l+i;++j,e=1ll*e*w%P) {
int t=1ll*a[j+i]*e%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}

一种简单的预处理是,每次对于每个分治大小,预处理依次系数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=1;j<i;++j) e[j]=1ll*e[j-1]*w%P;
//for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
//这个版本是沿用上一次预处理的结果,实际(只有)用这种预处理方法可以极大程度上加强FFT的精度
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}

另一种是在一开始就把所有的系数用一个数组存下来,具体过程可以描述为

对于每个分治长度n,我们只需要访问ωn0,ωn1,,ωnn21

那么对于分治长度n,我们在w数组的第n2 ~ n1项依次存储这些值

优化:我们只需要对于最大的分治长度处理,剩下的部分发现可以直接用折半引理访问得到

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
const int N=1<<21;
int a[N],w[N];
void Init(){
w[N>>1]=1;
int t=qpow(3,(P-1)/N);
rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,(N>>1)-1,1) w[i]=w[i<<1];
}
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}

三份代码在duck.ac上的评测结果表明,不预处理系数将近慢一倍

单组数据来看,预处理系数会慢一点

多组来看,预处理系数会快

实际差距不大,都可以使用

但是在某些层面来说,下面这份板子才是最好的(适用NTT,FFT且精度较高),不需要预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
 

拓展

1.分治+NTT

常用于处理多个计数背包的快速合并 (实际无权值01背包也是可以的)

我们可以用NTTnlogn合并两个大小为n的背包

分治时,每次合并两个分治子问题,总共的时间就是sizelogn

每个背包的size会被计算logn次,所以总共复杂度是nlog2n

 

2.CDQ+NTT

模板题传送门

对于形如dpi=j=0i1dpjgijdp转移(就是dp转移与差值有关)

由于求dpi时,需要保证dp0,dp1,,dpi1才能卷积,这个限制,我们可以用CDQ分治解决

对于当前分治区间[L,R]

依次考虑[L,mid]内部转移,[L,mid][mid+1,R]的转移(用FFT/NTT解决),[mid+1,R]内部转移

算法流程

1
2
3
4
5
6
7
void Solve(l,r){
if(l==r) return;
mid=(l+r)>>1;
Solve(l,mid);
(l,mid)->(mid+1,r);
Solve(mid+1,r);
}
 

3.MTT(任意模数NTT)

 

4.$n$元点值式

 

练习建议:

1.高精度乘法

2.简单应用:HDU-4609 题解

3.卷积构造模板: BZOJ-3527 题解

4.拓展卷积构造:HDU-5885 题解

5.构造卷积的应用:HDU-6061 题解

6.CDQ分治+FFTHDU-5730 题解

7.CDQ+NTT/降次前缀和优化dpHDU-5332 题解

8.容斥+MTTHDU-6088 题解

9.图上dp

联通图个数:BZOJ-3456 题解

带环联通图个数:HDU-5552 题解

森林数量和带限制森林数量:HDU - 5279 题解

10.点分治+FFT:CodeChef-PRIMEDST 题解

   

更多应用和优化参见毛啸2016论文

(如:两次FFT做卷积,4次FFT做MTT。。。)