FFT&NTT(以及扩展)

预备知识：用于NTT

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

前言

FT: 傅里叶变换:

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

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

FFT: 快速傅里叶变换

NTT: 快速数论变换

$\$

谈及核心思想

1.单位根:

构造$\omega_n$为$n$阶单位根(不知道$\omega_n$的值域)，满足性质$\omega_n^n=\omega_n^0=1$

对于$2|n$,$\omega _n^{\frac{n} {2} }=-1$

显然$\omega_n$满足一个非常简单的性质：折半引理 $\begin{aligned} \forall 2|i\and 2|n , \omega_n^i=\omega_{\frac{n} {2} }^{\frac{i} {2} }\end{aligned}$

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

2.多项式与点值式的转化

一个$n$阶多项式最普通的表示就是$F(x)=\sum_{i=0}^{n-1} a_ix^i$

然而，多项式也可以用$n$个互不相关的点表示，即$(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})$

两者可以互相转化

对于同$x_i$的点值，两个多项式卷积时，其$y_i$可以直接对应相乘

FFT/NTT的核心过程是

多项式$\longrightarrow$ 点值式$\longrightarrow$点值式对应相乘$\longrightarrow$多项式

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

3.分治思想

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

$\$

FFT的单位根

$(x,y)$指复数$i=\sqrt{-1},(x,y)=x+yi$

基本运算$(x,y)+(a,b)=(x+a,y+b),(x,y)\cdot (a,b)=(ax-by,ay+bx)$

FFT的单位根是:$\omega_n$=$(cos(\frac{2\pi} {n}),sin(\frac{2\pi} {n}))$

而$\omega _n^i=(cos(\frac{2\pi} {n}\cdot i),sin(\frac{2\pi} {n}\cdot i))$ (展开发现就是三角函数求和公式)

显然满足单位根的性质

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

$\$

NTT

相信您已经了解了原根的一些性质，$\text{NTT}$的单位根常用原根构造

$\text{NTT}$的单位根实际有较大的局限性，对于质数$P$只能构造出$n|P-1,\omega_n=g^{\frac{P-1} {n} }$

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

通常我们$P$取$998244353$，$2^{23}|(P-1)$，它的一个原根是3

实际上，为了满足下面分治需要，构造的模数通常满足$P-1=s\cdot 2^t$的$t$较大，这类模数我们常称作$\text{NTT}$模数

$\$ $\$

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

$\$

多项式转点值式

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

对于点值式，我们构造的点横坐标为$x_i=\omega_n^i$

具体目标是对于函数$F(x)$，求出在$x_0,x_1,\cdots ,x_{n-1}$上的函数值

即求出$F(x_i)=a_0\omega_n^0+a_1\omega_n^{i}+a_2\omega_n^{2i}+\cdots$

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

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

令$m=\frac{n} {2}$

$\begin{aligned} 2|i,G(x_i)=a_0\omega_{m}^0+a_2\omega_{m}^{\frac{i} {2} }+a_4\omega_{m}^{\frac{i} {2}\cdot 2}+\cdots\end{aligned}$ $\begin{aligned} 2|i,H(x_i)=a_1\omega_{m}^0+a_3\omega_{m}^{\frac{i} {2} }+a_5\omega_{m}^{\frac{i} {2}\cdot 2}+\cdots\end{aligned}$

更简洁的描述为

$i<m,G(x_i’)=a_0\omega_{m}^0+a_2\omega_{m}^{i}+a_4\omega_{m}^{2i}+\cdots$

$i<m,H(x_i’)=a_1\omega_{m}^0+a_3\omega_{m}^{i}+a_5\omega_{m}^{2i}+\cdots$

由于$G(x’_i),H(x’_i)$计算的是$[0,m-1]$项，而求$F(x_i)$时用到的是$0,2,4,\cdots$项，实际需要访问$G(x^2_i),H(x^2_i)$

和$F(x_i)$的式子比较，我们得到合并的式子为

$F(x_i)=G(x^2_i)+x_i H(x^2_i)$

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

$F(x_i)=G(x’_i)+x_i H(x’_i)$

注意$x_i=x’_{i\mod m}$

为了保证复杂度，尽量使得每次分治的子问题都分为两部分，这样的复杂度为$O(n\log n)$

附:实际上，分为$d$个子问题时，每次合并的复杂度为$O(n\cdot d)$，因此复杂度为

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

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

附:$d$个子问题时，设子问题答案为$G_j(x_i)$，则合并的式子为

$\begin{aligned} F(x_i)=\sum_{j=0}^{d-1}x_i^jG_j(x_i^d)=\sum_{i=0}^{d-1}x_i^jG_j(x'_{i\mod \frac{n} {d} })\end{aligned}$

点值式转多项式

一个简单的性质：单位根反演 $\sum_{j=0}^{n-1}\omega_n^{ij}= \left\{\begin{aligned} \frac{\omega_n^{in}-1} {\omega_n^i-1}=0 && i\ne 0\\ n && i=0\end{aligned} \right.$

设点值式对应$y_i$的序列为$b_i$

则$n\cdot a_i=\sum_{j=0}^{n-1}\omega_n^{-ij} b_j$

证明

$\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j=\sum_{j=0}^{n-1} \omega_n^{-ij}(\sum_{k=0}^{n-1}a_k\omega_n^{jk})\end{aligned}$ $\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j= \sum_{k=0}^{n-1}a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)} \end{aligned}$

由上面的式子，发现只有$k-i=0$时右边的求和式有值，所以上式成立

因此点值式转多项式直接把系数改为$\omega_n^{-i}$即可

$\$ $\$

Tips:

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

因此，如果乘法过后的多项式产生了超过$>n$的项$x^i$，会溢出到$x^{i\mod n}$

2.点值式并不是不满足除法，只是除法得到的多项式并不一定是一个$n$元以内的多项式，除了恰好整除的情况，得到的通常是一个无穷级数的式子，如$\begin{aligned} \frac{1} {1-x}=\frac{1-x^{\infty} } {1-x}=\sum_{i=0}^{\infty}x^i\end{aligned}$

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

$\$

代码实现与优化

模板题传送门

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

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

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);
}

由于$(\omega_n)^{\frac{n} {2} }=-1$，所以还可以简化为

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，最后输出要取整

蝴蝶优化

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

可以看到，分治时我们按照$i \mod 2$分成两组，然后继续分

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

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

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$:分治子问题大小为$2i$，$l$:合并区间的左端点为$l$，右端点为$l+2i$，$j$枚举合并位置

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;
            }
        }
    }
}

事实上我们还有更快的写法，就是将$\omega_n^i$预处理出来

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

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

$\$

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

由于每次求得点值是$\omega_n^{-i}=\omega_n^{n-i}$

所以可以直接用多项式转点值式的函数，最后把$[1,n-1]$这一段翻转，每个数除掉$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为例)

带入上面已经提到的优化，无预处理系数的$\text{NTT}$大概是这样的

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;
    }
}

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

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$，我们只需要访问$\omega_n^{0},\omega_n^{1},\cdots,\omega_n^{\frac{n} {2}-1}$

那么对于分治长度$n$，我们在$w$数组的第$\frac{n} {2}$ ~ $n-1$项依次存储这些值

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

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且精度较高)，不需要预处理

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背包也是可以的)

我们可以用NTT$n\log n$合并两个大小为$n$的背包

分治时，每次合并两个分治子问题，总共的时间就是$\sum size\log n$

每个背包的$size$会被计算$\log n$次，所以总共复杂度是$n \log ^2 n$

$\$

2.CDQ+NTT

模板题传送门

对于形如$dp_i=\sum_{j=0}^{i-1}dp_jg_{i-j}$的$dp$转移(就是dp转移与差值有关)

由于求$dp_i$时，需要保证$dp_0,dp_1,\cdots,dp_{i-1}$才能卷积，这个限制，我们可以用CDQ分治解决

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

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

算法流程

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$分治+$FFT$：HDU-5730 题解

7.$CDQ$+NTT/降次前缀和优化$dp$：HDU-5332 题解

8.容斥+$MTT$：HDU-6088 题解

9.图上$dp$：

联通图个数：BZOJ-3456 题解

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

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

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

$\$ $\$ $\$

orangejuice's blog

按ctrl会跳搜索，右上角选择day/night mode 似乎修好了搜索