FFT & NTT 学习笔记

标签：原根 int FFT NTT 笔记 Complex limit omega

FFT

FFT 是一种高效实现 DFT 和 IDFT 的方式，可以在 \(O(n \log n)\) 的时间内求多项式的乘法。

多项式的点值表示

不同于用每项的系数来表示一个多项式，我们知道对于给定的 \(n+1\) 个点值，可以确定唯一的 \(n\) 次多项式。这种用点值表示多项式的方法叫点值表示法。
如果知道 \(F(x)\) 和 \(G(x)\) 的点值表示，求出 \(F(x)\times G(x)\) 的点值表示是 \(O(n+m)\) 的。复杂度瓶颈变成了如何快速转换多项式的两种表示形式。
把系数转成点值的算法叫 离散傅里叶变换(DFT)，相应地，把点值转成系数的算法叫 逆离散傅里叶变换(IDFT)。

DFT

前置知识：单位根
好像高中数学的复数里没讲这玩意，当然也有可能是我学习偷工减料所以不知道！
OIWiki-复数-单位根

已经说过，要选取 \(n\) 个点并求出它们对应 \(F(x)\) 的值。先观察 \(F(x)\) 的性质：

\[F(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 \]

把奇偶项分开，得

\[F(x)=(a_0+a_2x^2+a_4x^4+a_6x^6)+(a_1x+a_3x^3+a_5x^5+a_7x^7)\\ 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)\\ F(x)=G(x^2)+x\times H(x^2) \]

似乎看起来对每个 \(x\) 求 \(G\) 和 \(H\)，要求的还是 \(n\) 个值。但这时单位根的性质就能用上了：

令 \(x\) 分别为 \(\omega_n^1,\omega_n^2,...,\omega_n^n\)，那么有：

\[\begin{aligned} F(\omega_n^k)&=G(\omega_n^{2k})+x\times H(\omega_n^{2k})\\ &=G(\omega_{n/2}^k)+x\times H(\omega_{n/2}^k) \end{aligned} \]

根据 \(\omega_n^k=-\omega_n^{k+n/2}\)，且 \(G,H\) 均为偶函数，可以得到：

\[F(\omega_n^k)=G(\omega_{n/2}^{k})-x\times H(\omega_{n/2}^{k}) \]

发现这事好的，因为只用求前一半的 \(G\) 和 \(H\) 的值就可以了，把它们两个分别递归处理即可。
但是 FFT 只能对次数为 \(2^m\) 的多项式这么做，所以如果多项式不是恰好 \(2^m\) 项，在后面补一些系数为 \(0\) 的项即可。

IDFT

想完成多项式乘法，还要把点值表示重新转换回系数表示，这个过程叫 IDFT。
构造方法实在过于神秘，以至于菜猫猫并不能理解这个思路的由来，所以这里直接给出做法：
设对 \(F(x)\) 进行 DFT 后得到点值表示 \(y_0,y_1,...y_{n-1}\)，构造 \(G(x)=y_0+y_1x+...+y_{n-1}x^{n-1}\)。
将 \(x=\omega_n^0,\omega_n^{-1}...\omega_n^{-(n-1)}\) 代入，得：

\[\begin{aligned} G(\omega_n^{-k})&=\sum\limits_{i=0}^{n-1} y_i \omega_n^{-ki}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} a_j\times (\omega_n^i)^j\times \omega_n^{-ki}\\ &=\sum_{j=0}^{n-1}a_j \sum_{i=0}^{n-1}(\omega_n^{j-k})^i \end{aligned} \]

那么当且仅当 \(j=k\)，\(G=a_k\times n\)。因此对 \(G(x)\) 代入 \(x=\omega_n^0,\omega_n^{-1}...\omega_n^{-(n-1)}\) 求 DFT，结果除以 \(n\) 即为多项式的各项系数。
实现上可以设定一个参数为 \(1\) 或 \(-1\) 传入 DFT 函数，减小码量。

实现-递归版

众所周知 c++ 自带的 complex 类常数极大，建议手写。
至此，我们可以写出朴素的递归版本 FFT 了：

const int N=4e6+5;
typedef double db;
int n,m;
db pai=acos(-1.0);
struct Complex{
	db x,y;
}a[N],b[N];
Complex operator +(Complex a,Complex b) {return {a.x+b.x,a.y+b.y};}
Complex operator -(Complex a,Complex b) {return {a.x-b.x,a.y-b.y};}
Complex operator *(Complex a,Complex b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(int limit,Complex a[],int op)
{
	if(limit==1) return;
	Complex a1[(limit>>1)+5],a2[(limit>>1)+5];
	for(int i=0;i<=limit;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];
	FFT(limit>>1,a1,op),FFT(limit>>1,a2,op);
	Complex Wn={cos(2.0*pai/limit),sin(2.0*pai/limit)*op},w={1,0};
	for(int i=0;i<(limit>>1);i++,w=w*Wn)
	{
		a[i]=a1[i]+w*a2[i];
		a[i+(limit>>1)]=a1[i]-w*a2[i];
	}
}
int main()
{
	n=read(),m=read();
	for(int i=0;i<=n;i++) a[i].x=read();
	for(int i=0;i<=m;i++) b[i].x=read();
	int w=1;while(w<=m+n) w<<=1;
	FFT(w,a,1),FFT(w,b,1);
	for(int i=0;i<=w;i++) a[i]=a[i]*b[i];
	FFT(w,a,-1);
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/w+0.5));
	return 0;
}

位逆序置换

递归的常数太大了，接下来讲讲怎么不递归。
考虑每次把系数按奇偶分开的过程，如果 \(i\) 的二进制翻转后为 \(j\)，那么下标为 \(i\) 的数操作到最后就会跑到位置 \(j\) 上。
这启发我们使用迭代实现，直接从底层倒推即可。

const int N=4e6+5;
const double pi=acos(-1);
int n,m,limit=1;
struct Complex {double x,y;}a[N],b[N];
il Complex operator +(Complex a,Complex b) {return {a.x+b.x,a.y+b.y};}
il Complex operator -(Complex a,Complex b) {return {a.x-b.x,a.y-b.y};}
il Complex operator *(Complex a,Complex b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
int to[N];
void FFT(Complex *a,int tp)
{
	for(int i=0;i<limit;i++) if(i<to[i]) swap(a[i],a[to[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		Complex Wn={cos(pi/mid),sin(pi/mid)*tp};
		for(int j=0;j<limit;j+=(mid<<1))
		{
			Complex w={1,0};
			for(int k=0;k<mid;k++,w=w*Wn)
			{
				Complex x=a[j+k],y=w*a[j+mid+k];
				a[j+k]=x+y,a[j+mid+k]=x-y;
			}
		}
	}
}
int main()
{
	n=read(),m=read();
	for(int i=0;i<=n;i++) a[i].x=read();
	for(int i=0;i<=m;i++) b[i].x=read();
	int k=0; while(limit<=n+m) limit<<=1,k++;
	for(int i=0;i<limit;i++) to[i]=(to[i>>1]>>1)|((i&1)<<(k-1));
	FFT(a,1),FFT(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/limit+0.5));
	return 0;
}

NTT

前置知识-原根

阶：使得 \(a^n\equiv1\pmod{m}\) 的最小正整数 \(n\) 称为 \(a\) 模 \(m\) 的阶。
原根：若 \(\gcd(g,m)=1\)，且 \(\delta_m(g)=\varphi(m)\)，则 \(g\) 是 \(m\) 的原根。
存在定理 一个数 \(m\) 存在原根当且仅当 \(m=2,4,p^{\alpha},2p^{\alpha}\)，其中 \(p\) 为奇素数，\(\alpha \in \mathbb{N}^{*}\)。
判定定理 设 \(m\ge 3,\gcd(a,m)=1\)，则 \(a\) 是模 \(m\) 的原根的充要条件是，对于 \(\varphi (m)\) 的每个质因数 \(p\)，都有 \(a^{\frac{\varphi(m)}{p}}\not\equiv 1\pmod m\)。
原根个数 若一个数有原根，则它原根的个数为 \(\varphi(\varphi(m))\)。
以上结论的证明：link

NTT 用来解决模意义下的多项式乘法问题，其运算均为整数，在常数和精度方面均高于 FFT。
由原根的性质，我们把 FFT 中所有 \(\omega_n\) 替换成 \(g^{\frac{p-1}{n}}\) 即可。
常见模数是 \(998244353\)，它的原根是 \(3\)。

为什么我的 ntt 跑得比 fft 慢好多啊 /dk

#define int long long
const int N=4e6+5,mod=998244353;
il int qpow(int n,int k=mod-2)
{
	int res=1;
	for(;k;n=n*n%mod,k>>=1) if(k&1) res=res*n%mod;
	return res;
}
int n,m,a[N],b[N],inv=qpow(3),limit=1,to[N];
il void NTT(int *a,int tp)
{
	for(int i=0;i<limit;i++) if(i<to[i]) swap(a[i],a[to[i]]);
	for(int len=1;len<limit;len<<=1)
	{
		int Wn=qpow(tp>0?3:inv,(mod-1)/(len<<1));
		for(int i=0;i<limit;i+=(len<<1))
			for(int j=0,w=1;j<len;j++,w=w*Wn%mod)
			{
				int x=a[i+j],y=w*a[i+len+j]%mod;
				a[i+j]=(x+y)%mod,a[i+len+j]=(x-y+mod)%mod;
			}
	}
}
signed main()
{
	n=read(),m=read();
	for(int i=0;i<=n;i++) a[i]=read();
	for(int i=0;i<=m;i++) b[i]=read();
	int k=0; while(limit<=m+n) limit<<=1,k++;
	for(int i=0;i<limit;i++) to[i]=(to[i>>1]>>1)|((i&1)<<(k-1));
	NTT(a,1),NTT(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i]%mod;
	NTT(a,-1);
	for(int i=0;i<=n+m;i++) printf("%lld ",a[i]*qpow(limit)%mod);
	return 0;
}

完结撒花！

标签：原根,int,FFT,NTT,笔记,Complex,limit,omega
From： https://www.cnblogs.com/ying-xue/p/17676005.html