FFT 和 NTT

FFT:https://www.luogu.com.cn/problem/P3803
NTT:https://www.luogu.com.cn/problem/P4245

日常抄 oi-wiki

FFT

多项式乘法对多项式的系数向量进行卷积：\(\sum_{i=1}^{n+m}\sum_{j=1}^i f(j)g(i-j)\)
DFT 可以认为是将多项式的系数表示法变成点值表示法，然后将点值分别相乘即可得到 \(f\times g\) 的点值表示，再用 IDFT 变换回去即可
而 FFT 便是快速实现 DFT 和 IDFT 的过程

定义 \(x^n=1\) 复数域的 \(n\) 个解为 \(n\) 次单位复根，表示为 \(w_n\)。有 \(w_n=e^{\frac{2\pi i}{n}}=\cos(\frac{2\pi}{n})+i\sin(\frac{2\pi}{n})\)，也就是把复平面下的单位圆均分成了 \(n\) 份
那么 \(x^n=1\) 的解就可以表示为 \(w_n^k,k=0,1,2,\cdots,n-1\)
有性质：

• \(w_n^n=1\)
• \(w_n^k=w_{2n}^{2k}\)
• \(w_{2n}^k=-w_{2n}^{k+n}\)

DFT 在处理时将多项式按照次数的奇偶分治来求在每个 \(w_n^k\) 处的取值，比如有多项式：\(f(x)=\sum_{i=0}^7 a_ix^i\)
那么按照奇偶分类：\(f(x)=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6)\)
然后分别定义 \(g(x)=a_0+a_2x+a_4x^2+a_6x^3),h(x)=a_1+a_3x+a_5x^2+a_7x^3\)
所以原式变成了 \(g(x^2)+x\cdot h(x^2)\)
如果分治的求出 \(g(x),h(x)\) 的点值，那么：

\[\begin{align} \operatorname{DFT}(f(w_n^k))&=\operatorname{DFT}(g((w_n^k)^2))+w_n^k\cdot \operatorname{DFT}(h((w_n^k)^2))\\ &=\operatorname{DFT}(g(w_n^{2k}))+w_n^k\cdot \operatorname{DFT}(h(w_n^{2k}))\\ &=\operatorname{DFT}(g(w_{n/2}^{k}))+w_n^k\cdot \operatorname{DFT}(h(w_{n/2}^{k}))\\ \end{align} \]

\[\begin{align} \operatorname{DFT}(f(w_n^{k+n/2}))&=\operatorname{DFT}(g(w_n^{2k+n}))+w_n^{k+n/2}\cdot \operatorname{DFT}(h(w_n^{2k+n}))\\ &=\operatorname{DFT}(g(w_{n/2}^{k}))-w_n^k \cdot \operatorname{DFT}(h(w_{n/2}^k)) \end{align} \]

然后现在已经有了一个递归的做法，需要多项式的长度为 \(n=2^k\)，所以只要对输入的数补齐到第一个大于等于它的 2 的幂次就行
然后带入 \(w_n^0,w_n,w_n^2,\cdots,w_n^{n-1}\) 求值即可

更常用的是迭代版的
考虑分治的过程，当次数为 8 时，分治过程如下：

step 1: 0,1,2,3,4,5,6,7
step 2: 0,2,4,6; 1,3,5,7
step 3: 0,4; 2,6; 1,5; 3,7;
step 4: 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\) 为第 \(i\) 为要交换到的地方，那么讨论 \(i\) 的最后一个二进制位是 \(0\) 或 \(1\) 即可得知 \(rev_i\) 的最高二进制位，然后再根据 \(rev_{i/2}\) 的值进行递推：

for(reg int i=0;i<n;i++) rev[i]=rev[i>>1]>>1,rev[i]|=(i&1)?(n>>1):0;

然后提前交换就有了迭代版的 FFT

---

然后就是它的逆变换，将点值表示转化为系数表示
\(f(x)\) 组成的向量可以由系数组成的向量乘一个关于带入的值（也就是 \(w_n^k\)）的矩阵得出，而 FFT 的过程实际上就是快速求了这样的一个矩阵和向量的乘积，具体来说就是：

\[\begin{bmatrix} f(w_n^0)\\ f(w_n)\\ f(w_n^2)\\ \vdots\\ f(w_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1&1&1&\cdots&1\\ 1&w_n&w_n^2&\cdots&w_n^{n-1}\\ 1&w_n^2&w_n^4&\cdots&w_n^{2(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&w_n^{n-1}&w_n^{2(n-1)}&\cdots&w_n^{(n-1)^2} \end{bmatrix} \times \begin{bmatrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n-1} \end{bmatrix} \]

现在知道等号左边的向量和乘的矩阵，要反推右边的系数向量，那么只要给中间的矩阵求逆即可
那么直接套用 \(O(n^3)\) 求逆的方法
其实这个矩阵的逆矩阵可以直接构造出来，每一位取倒数再除以 \(n\) 即可

然后只要给左边的 \(f(x)\) 向量用 \(w_n^{-k}=\cos(\frac{2\pi k}{n})-i\sin(\frac{2\pi k}{n})\) 构成的矩阵再做一次 FFT 就是逆变换了
最终在结果上再除以 \(n\) 就行

const double PI=acos(-1);
struct Complex{
	double a,b;
}a[N];
inline Complex operator + (Complex x,Complex y){return (Complex){x.a+y.a,x.b+y.b};}
inline Complex operator - (Complex x,Complex y){return (Complex){x.a-y.a,x.b-y.b};}
inline Complex operator * (Complex x,Complex y){return (Complex){x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a};}
int rev[N];
inline void init(reg int n){
	for(reg int i=0;i<n;i++) rev[i]=rev[i>>1]>>1,rev[i]|=(i&1)?(n>>1):0;
}
inline void fft(int n,Complex *a,int type){
	for(reg int i=0;i<n;i++)if(rev[i]<i) std::swap(a[i],a[rev[i]]);
	Complex wn,w,o;
	for(reg int h=1;h<n;h<<=1){
		wn.a=cos(2*PI/(h<<1));wn.b=sin((type?2:-2)*PI/(h<<1));
		for(reg int i=0;i<n;i+=h<<1){
			w.a=1;w.b=0;
			for(reg int j=i;j<i+h;j++,w=w*wn){
				o=w*a[j+h];
				a[j+h]=a[j]-o;a[j]=a[j]+o;
			}
		}
	}
	if(!type) for(reg int i=0;i<n;i++) a[i].a/=n,a[i].b/=n;
}

NTT

如果求得的系数需要取模就不能直接 FFT 了
所以考虑到对于质数 \(p=qn+1,n=2^m\)，可以取其原根 \(g\)，满足 \(g^{qn}\equiv 1 \pmod p\)，然后就可以将 \(g_n=g^q\) 等价的看作 \(w_n\)
更进一步，\(g_k=g^{(p-1)/k}\) 可以等价的看作 \(w_k\)
发现它也满足之前说的 \(w_n\) 的性质：\(g_n^n\equiv 1\pmod p,g_{2n}^{2k}\equiv g_{n}^k\pmod p,g_{2n}^k\equiv g_{2n}^{n+k}\pmod p\)
所以就如果选取的质数有足够多的 \(2\) 的因数，就可以直接把 FFT 中的 \(w_n\) 换成 \(g^{(p-1)/n}\) 来计算
满足要求的质数有：998244353,1004535809,469762049

对于任意模数 NTT，可以直接选取以上三个质数分别做一次普通 NTT，然后再 CRT 合并，值域一般够用

放一个任意模数 NTT 的代码

#define G 3
#define mod1 998244353
#define mod2 1004535809
#define mod3 469762049
#define N 270006
int n,m,p;
inline int power(reg int a,reg int b,reg int mod){
	reg int ans=1;
	while(b){
		if(b&1) ans=(long long)ans*a%mod;
		b>>=1;a=(long long)a*a%mod;
	}
	return ans;
}
const int inv1=power(mod1,mod2-2,mod2),inv2=power((long long)mod1*mod2%mod3,mod3-2,mod3);
inline int Mod(reg int a,reg int mod){return a>=mod?a-mod:a;}
int rev[N];
int f1[N],f2[N],f3[N],g1[N],g2[N],g3[N];
inline void init(reg int n){
	for(reg int i=0;i<n;i++) rev[i]=rev[i>>1]>>1,rev[i]|=(i&1)?(n>>1):0;
}
inline void ntt(reg int n,int *a,int type,int mod){
	for(reg int i=0;i<n;i++)if(rev[i]<i) std::swap(a[i],a[rev[i]]);
	for(reg int h=1;h<n;h<<=1){
		int gn=power(G,(mod-1)/(h<<1),mod),g,o;
		if(!type) gn=power(gn,mod-2,mod);
		for(reg int i=0;i<n;i+=h<<1){
			g=1;for(reg int j=i;j<i+h;j++,g=(long long)g*gn%mod){
				o=(long long)g*a[j+h]%mod;
				a[j+h]=Mod(a[j]-o+mod,mod);a[j]=Mod(a[j]+o,mod);
			}
		}
	}
	if(!type){
		int inv=power(n,mod-2,mod);
		for(reg int i=0;i<n;i++) a[i]=(long long)a[i]*inv%mod;
	}
}
inline int get(long long a,long long b,long long c){
	long long x=(b-a+mod2)%mod2*inv1%mod2*mod1+a;
	return ((c-x%mod3+mod3)%mod3*inv2%mod3*mod1%p*mod2%p+x)%p;
}
int main(){
	n=read();m=read();p=read();
	for(reg int i=0;i<=n;i++) f1[i]=f2[i]=f3[i]=read()%p;
	for(reg int i=0;i<=m;i++) g1[i]=g2[i]=g3[i]=read()%p;
	int max=1;while(max<=n+m) max<<=1;
	init(max);
	ntt(max,f1,1,mod1);ntt(max,f2,1,mod2);ntt(max,f3,1,mod3);
	ntt(max,g1,1,mod1);ntt(max,g2,1,mod2);ntt(max,g3,1,mod3);
	for(reg int i=0;i<max;i++)
		f1[i]=(long long)f1[i]*g1[i]%mod1,f2[i]=(long long)f2[i]*g2[i]%mod2,f3[i]=(long long)f3[i]*g3[i]%mod3;
	ntt(max,f1,0,mod1);ntt(max,f2,0,mod2);ntt(max,f3,0,mod3);
	for(reg int i=0;i<=n+m;i++) printf("%d ",get(f1[i],f2[i],f3[i]));
	return 0;
}