多项式

FFT

首先，多项式有两种表示方法。

系数表示法

就是普通的表示多项式的方法

用\(A(n)\)表示\(n\)次多项式，则

\(A(n)=\sum\limits_{i=0}^{n}a_ix^i\)

形象一点\(\Downarrow\)

即：

但时间复杂度是 \(d^2\) ,不太优。

点值表示法

我们可以发现，
\(d\) 次多项式都可以用 \(d+1\) 个点表示出来

则：

所以说，\(FFT\)的思路是

即：

1. 计算点值\((DFT)\)
   系数 \(\Rightarrow\) 点值

2. 将点值乘起来

3. 用点值还原成多项式的系数\((IDFT)\)
   点值 \(\Rightarrow\) 系数

---

DFT

为了方便计算我们不妨假设 \(n = 2^k\) . 记

\(ω_n = e^{\frac{2\pi i}{n}}\)

单位圆上的点：\(z=\cos\theta+i\sin\theta(0\le\theta<2\pi)\)

单位根\(\omega_n^k=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n}\)

\(\omega\) 具有周期性 \(\omega_n^{k+n}=\omega_n^k\) , 对称性 \(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\) , 折半性 \(\omega_n^{2k}=\omega_{n/2}^k\)

点值位置的选取为 \(n\) 个 \(n\) 次单位根

\(ω^0_n ,ω^1_n ,\dots,ω^{n-1}_n\)

记

\(A(x) =∑\limits^{n-1}_{i=0}x^i = A_1 (x^2 ) + xA_2 (x^2 )\)

其中

\(A_1(x) = a_0 + a_2x + a_4x^2 + ··· + a_{n-2}x^{\frac{n}{2} -1}\)

\(A_2(x) = a_1 + a_3x + a_5 x^2 + ··· + a _{n-1}x^{\frac{n}{2} -1}\)

代入 \(x = ω^k_n\):

\(A(ω^k_n ) = A_1 ((ω^k_n )^2 ) + ω_k^n A_2 ((ω^k_n )^2 )= A_1 (ω^k_{n/2} ) + ω^k_n A_2 (ω^k_{n/2})\)

另一方面, 代入 \(x = ω^{k+\frac{n}{2}}\)有:

\(A(ω^{k+\frac{n}{2}}_n) = A_1 (ω^k_{n/2}) -ω^k_n A_2 (ω^k_{n/2})\)

也就是说，如果我们已经求出了\(A_1 ,A_2\) 在 \(n/2\) 次单位根的取值，那么我们可以通过 \(O(n)\) 计算出 \(A\) 在 \(n\) 次单位根的取值。

根据主定理

\(T(n) = 2T(n/2) + O(n) = O(n\log n)\)

通过分治我们解决了 \(DFT\) 的问题

IDFT

从另一个角度看 \(DFT\), 我们实际上是把一个单位根矩阵\(W = [ω^{ij}]\) 左乘在了系数向量上, 那么实现 \(IDFT\) 实际上我们需要找到这个单位根矩阵的逆。

\(W'= [ω^{-ij}]\).乘起来是一个单位阵的常数倍

于是

\(W^{-1} =\frac{W'}n\)

\(\omega_n^k=\cos \theta + i\sin\theta\)

\(\omega_n^k=\frac{1}{\cos \theta+i\sin \theta}=\cos \theta -i\sin \theta=\omega_n^{-k}\)

逆变换带入单位根的倒数等于虚部变为负号。

然后我们就发现, \(IDFT\) 就是把 \(DFT\) 里面的 \(ω\) 换成 \(ω^{-1}\) ，最后除以 $n $就行了。

推导过程：

\(A(\omega_n^{-k})=\sum\limits_{i=0}^{n-1}f(\omega_n^i)\omega_n^{-ik}\)

\(=\sum\limits_{i=0}^{n-1}\omega_n^{-ik}\sum\limits_{j=0}^{n-1}a_j(\omega_n^i)^j\)

\(=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j\omega_n^{i(j-k)}\)

记\(S(\omega_n^a)=\sum\limits_{i=0}^{n-1}(\omega_n^a)^i\)

当\(a\bmod n=0\) 时，\(S(\omega_n^a)=n\)

当\(a\bmod n\not=0\) 时，\(S(\omega_n^a)=n\)

\(S(\omega_n^a)=\sum\limits_{i=0}^{n-1}(\omega_n^a)^i\)

\(\omega_n^aS(\omega_n^a)=\sum\limits_{i=1}^{n}(\omega_n^a)^i\)

\(S(\omega_n^a)=\frac{(\omega_n^a)^n-(\omega_n^a)^0}{\omega_n^a-1}=0\)

即：

\(S(\omega_n^a)=\left\{\begin{matrix} n,\qquad a=0\\ 0,\qquad a\not=0 \end{matrix}\right.\)

代回可得：

\(A(\omega_n^{-k})=\sum\limits_{j=0}^{n-1}a_jS(\omega_n^{j-k})=a_k\times n\)

\(Code:\)

递归实现

#include<iostream>
#include<cmath>
#include<cstdio>
#define db double
using namespace std;
const int N=4e6+10;
const db pi=acos(-1.0);
int n,m;
struct complex{
	db x,y;
}a[N],b[N];
complex operator + (complex a, complex b) {return complex{a.x+b.x,a.y+b.y};}
complex operator - (complex a, complex b) {return complex{a.x-b.x,a.y-b.y};}
complex operator * (complex a, complex b) {return complex{a.x*b.x-a.y*b.y,b.x*a.y+b.y*a.x};}
void FFT(int lmt,complex *a,int type){
	if(lmt==1) return ;
	complex a1[lmt],a2[lmt];
	for(int i=0;i<=lmt;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];
	FFT(lmt>>1,a1,type);
	FFT(lmt>>1,a2,type);
	complex Wn={cos(2.0*pi/lmt),type*sin(2.0*pi/lmt)},w=complex{1,0};
	for(int i=0;i<(lmt>>1);++i){
		a[i]=a1[i]+w*a2[i],a[i+(lmt>>1)]=a1[i]-w*a2[i];
		w=w*Wn;
	}
}
int main(){
	cin>>n>>m;
	for(int i=0;i<=n;++i) cin>>a[i].x;
	for(int i=0;i<=m;++i) cin>>b[i].x;
	int lmt=1;while(lmt<=n+m) lmt <<= 1;
	FFT(lmt,a,1);
	FFT(lmt,b,1);
	for(int i=0;i<=lmt;++i) a[i]=a[i]*b[i];
	FFT(lmt,a,-1);
	for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x / lmt + 0.5));
	return 0;
}

迭代实现

蝴蝶变换：每个位置的元素下标做了二进制翻转，这个变换称为位逆序变换（即蝴蝶变换

可以预处理出翻转后的数。

设\(n=2^k,k\) 表示二进制的长度

设 \(R[x]\) 表示长度为 \(k\) 的二进制数 \(x\) 翻转后的数，要求出 \(R[0],R[1]\dots R[n-1]\)

求\(R[x]\)时 ，发现\(R[\frac{x}{2}]\)已求出，所以求\(R[]\)的时候，先把 \(x\) 右移一位（去掉个位上的数），然后翻转，再右移（去掉高位补的0），如果 \(x\) 的个位是 0，翻转后最高位就是0，如果是1，则为1，因此还要加上\(\frac{n}{2}=2^{k-1}\)

即：\(R[x]=\frac{R[\frac{x}{2}]}{2}+[x\)&\(1]\times \frac{n}{2}\)

#include<iostream>
#include<cstdio>
#include<cmath>
#define db double
#define il inline
using namespace std;
const int N=4e6+10;
const int pi=acos(-1.0);
struct complex{
	db x,y;
}a[N],b[N];
int n,m,limit=1,l,r[N<<1];
int read(){
	int f=0,x=0;char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f|=(ch=='-');
	for(; isdigit(ch);ch=getchar()) x=(x<<1)+(x<<3)+(ch^48);
	return f?-x:x;
}
complex operator + (complex a, complex b) {return complex{a.x+b.x,a.y+b.y};}
complex operator - (complex a, complex b) {return complex{a.x-b.x,a.y-b.y};}
complex operator * (complex a, complex b) {return complex{a.x*b.x-a.y*b.y,b.x*a.y+b.y*a.x};}
void FFT(complex *a,int type){
	for(int i=0;i<limit;++i) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int m=2;m<=limit;m<<=1){
		complex Wn={cos(2.0*pi/m),type*sin(2.0*pi/m)};
		for(int i=0;i<limit;i+=m){
			complex w={1,0};
			for(int j=0;j<m/2;++j){
				complex x =a[i+j],y=w*a[i+j+(m>>1)];
				a[i+j]=x+y,a[i+j+(m>>1)]=x-y;
				w=w*Wn;
			}
		}
	}
	return ;
}
int main(){
	n=read();m=read();
	for(int i=0;i<=n;++i) cin>>a[i].x;
	for(int i=0;i<=m;++i) cin>>b[i].x;
	while(limit<=n+m) limit<<=1,l++;
	for(int i=0;i<limit;++i) r[i]=(r[i>>1]>>1|(i&1)<<(l-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

前置知识：欧拉定理、阶、原根。

\(NTT\) 是 \(FFT\) 在数论基础上的实现。

数论变换时一种计算卷积的快速算法。

类似于 \(FFT\)，只是 \(FFT\) 利用的单位圆和单位根，\(NTT\) 用的是原根

我的理解：

\(a^n\equiv1(\bmod p)\) 中 最小的正整数\(n\)即为阶，记\(\delta_p(a)\)：\(a\) 模 \(p\) 的阶。

而 \(\delta_p(a)=\varphi(p)\) , 则称 \(a\) 为模 \(p\) 的一个原根。

为了多次二分，模数 \(p\) 应选择形如\(q\times 2^k+1\)的质数，其中，\(q\)为奇素数，\(k\) 为整数。

因 \(p\) 是质数，故 \(\varphi(p)=p-1\) , 则\(g^0,g^1,\dots g^n-1\) 在模 \(p\) 意义下两两不同，从中选取对称的\(n(2^b)\)个值：

\(g_n^0=1,g_n^1=g^{\frac{p-1}{n}},\dots,g_n^k=g^{\frac{p-1}{n}k},\dots,g_n^{n-1}\)

以下证明均省略\(\bmod p\)

\(g\) 是模 \(p\) 的原根，\(g_n^k=g^{\frac{p-1}{n}k}\)

原根的性质：

1. 指数性

   \(g_n^kg_n^m=g^{\frac{p-1}{n}k}g^{\frac{p-1}{n}m}=g_n^{k+m}\)

2. 周期性

   \(g_n^{k+n}=g_n^kg_n^n=g_n^k(g^{p-1})=g_n^kg^{\varphi(p)}=1\)

3. 对称性

   \(g_n^{k+\frac{n}{2}}=-g_n^k\)

4. 折半性

   \(g_n^{2k}=g_{n/2}^k\)

原根的倒数 \(g_n^{-k}\) 就是原根 \(g_n^k\) 模 \(p\) 的逆元。

剩下的感觉和 \(FFT\) 差不多，但是要注意精度和模数。

\(Code:\)

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=3e6;
const int mod=998244353;
const int g=3;
ll a[N],b[N];
ll qpow(ll a,ll x){
	ll ans=1;
	while(x){
		if(x&1) ans=ans*a%mod;
		a=(a*a)%mod;
		x>>=1;
	}
	return ans%mod;
}
int n,m,limit=1,r[N],l;
void NTT(ll *a,int type){
	for(int i=0;i<limit;++i) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int m=1;m<limit;m<<=1){
		ll Wn=qpow(g,(mod-1)/(m<<1));
		if(type==-1) Wn=qpow(Wn,mod-2);
		for(int i=0;i<limit;i+=(m<<1)){
			ll w=1;
			for(int j=0;j<m;++j,w=(w*Wn)%mod){
				ll x=a[i+j],y=w*a[i+j+m]%mod;
				a[i+j]=(x+y)%mod;a[i+j+m]=(x-y+mod)%mod;
			}
		}
	}
	return ;
}
int main(){
	cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i];
	for(int i=0;i<=m;i++) cin>>b[i];
	while(limit<=m+n) limit<<=1,++l;
	for(int i=0;i<limit;++i) r[i]= (r[i>>1]>>1 | (i&1)<<(l-1));
	NTT(a,1);
	NTT(b,1);
	for(int i=0;i<=limit;++i) a[i]=(a[i]*b[i])%mod;
	NTT(a,-1);
	ll inv=qpow(limit,mod-2);
	for(int i=0;i<=n+m;++i) cout<<a[i]*inv%mod<<" ";
	return 0;
}