这篇文章上次修改于 600 天前,可能其部分内容已经发生变化,如有疑问可询问作者。
卷积:$h_i=\sum _{j=0}^{i}f_j \times g_{i-j}$,其中 $h_i$ 是多项式 $H$ 的 $i$ 次项系数。其实就是多项式乘法。
FFT(快速傅里叶变换)
思想
直接朴素乘会导致太慢了,时间复杂度直逼 $\mathcal{O}(n^2)$,仍然停留在多项式的系数表达式太浅薄了,可以尝试一下点值表达式。
对于一个有 $n$ 项的多项式,只要取不同的 $n$ 个值即可确定这个多项式,如果直接把两个多项式同时取 $n$ 个值,再把相同 $x$ 值得到的结果相乘,最后转换为系数表达,时间复杂度不会有任何改进,所以可以从取值的方面来下功夫。
前置知识
1. 虚数:定义虚数单位 $i$ ,保证 $i^2=-1$,则可以把一个虚数定义为 $a+bi$,其中 $a,b$ 是实数。
2. 数轴:定义一个平面直角坐标系,横轴为实轴,纵轴为虚轴,这个坐标系定义是复数平面,任意一个点都代表一个实数或虚数,横坐标为其实部,纵坐标为其虚部。
3. 模长、幅角:接下来以 $(x,y)$ 代表从 $(0,0)$ 到 $(x,y)$ 的向量,模长记为 $|(x,y)|$ 为 $(x,y)$ 到 $(0,0)$ 的欧几里得距离,幅角定义为正实部逆时针旋转一个幅角即可与向量重合。
3. 虚数的运算:
加减法直接把实部和虚部分别相加即可。
乘法在代数意义上的运算:$(a+bi)\times(c+di)$=$ac+adi+bci-bd$=$(ac-bd)+(ad+bc)i$
几何意义:幅角相加,模长相乘。
4.单位根
定义 $\omega_n$ 为 $\omega_n^n=1$,很容易推出其中一个单位根为模长为1,幅角为$360/n$,很容易想到这个数自乘 $n$ 次,旋转 $360$度,与正实轴重合,模长为1,$\omega_n^k$也满足${(\omega_n^k)}^n=1$,证明同理。
5.单位根性质
$\omega_n^k=\omega_{n/2}^{k/2}$ ,证明:$\omega_n^k$为模长为1,幅角为$\frac{360k}{n}$度,$\omega_{n/2}^{k/2}$也是一样的。
$\omega_n^{k}=-\omega_n^{k+n/2}$, 证明:多转了180度,所以就是相反数。
$\omega_n^0=\omega_n^n=1$,由定义得。
$\omega_n^{k}=\omega_n^{n+k}$,证明:转360度转回来了。
正向FFT(系数到点值)
可以代入所有的单位根,为了统一计算,把所有的多项式补全到 $2^k$,其中k是满足 $2^k \leq n+m+1$的最小值。高次项全部补0即可。
$A(x)=(a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1})$
按照次数奇偶性分组,$A(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x+a_3x^3+…+a_{n-1}x^{n-1})$
继续定义 B(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{(n-2)/2}$
$C(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{(n-2)/2}$
则 $A(x)=B(x^2)+xC(x^2)$
将 $\omega_n^k$ 代入原式得 $A(\omega_n^k)=B(\omega_n^{2k})+\omega_n^kC(\omega_n^{2k})$,因为相比原式项数减半所以把 $\omega_n^{2k}$ 化为 $\omega_{n/2}^k$得 $A(\omega_n^k)=B(\omega_{n/2}^k)+\omega_n^kC(\omega_{n/2}^k)$
$B,C$递归处理即可。
IFFT(逆变换)
上面算出来了点值表示法,如何转成系数表达式捏。
设原多项式的点值表达式在 $\omega_n^k(0 \leq k<n)$ 得到 $y_0,y_1,…,y_{n-1}$
设多项式 $D(x)=y_0+y_1x+y_2x^2+y_3x^3+…y_{n-1}x^{n-1}$。
设在 $\omega_n^{-k}(0 \leq k<n)$的 点值表达式是 $z_0,z_1,…,z_{n-1}$。
可证得$c_k=\sum_{i=0}^{n-1}a_j(\sum^{n-1}_{i=0}(\omega_n^{j-k})^i)$
核心的就是$\sum^{n-1}_{i=0}(\omega_n^{j-k})^i$,设$S(x)=\sum^{n-1}_{i=0}x^i$
易证 $S(x)=\frac {1-1}{w_n^k-1}$
当$w_n^k!=1$时,$S(x)=0$
当$w_n^k=1$时,$S(x)=n$
代入得$c_k=na_k$
所以 $a_k=c_k/n$
优化
递归实现容易爆栈,而且极其慢,考虑优化为迭代实现。
优化为迭代实现的瓶颈在于如何划分系数,不过可以非常简单的把需要合并在一起的项放在一起。
只有一项的多项式,点值表达式和系数表达式是一样的。
因为上述按奇偶性分开,其实是按照下标二进制表达从后往前来比较的来划分的,相同的会放在一起,所以可以直接按照下表二进制反转的顺序来依次合并。
#include <iostream>
#include <cstdio>
#include <complex>
#include <cmath>
#define int long long
inline int read();
const double Pi=acos(-1.0);
std::complex <double> f[8000005],g[8000005],ans[8000005];
int r[8000005];
void fft(std::complex <double> *f,int limit,int type);
signed main(){
int n=read(),m;
std::complex <double> temp=f[0];
m=read();
int limit=1,k=0;
while(limit<(n+m+2)){
limit=(limit)<<1;
k++;
}
for(int i=0;i<=n;i++){
f[i]=read();
}
for(int i=0;i<=m;i++){
g[i]=read();
}
for(int i=0;i<=limit;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
}
fft(f,limit,1);
fft(g,limit,1);
for(int i=0;i<=limit;i++){
f[i]*=g[i];
}
fft(f,limit,-1);
for(int i=0;i<=n+m;i++){
printf("%lld ",(int)(f[i].real()/(limit)+0.5));
}
return 0;
}
inline int read(){
int x=0,f=1;char c=getchar();
while(c<'0'||c>'9'){
c=='-'?f=-1:1;
c=getchar();
}
while(c>='0'&&c<='9'){
x=(x<<3)+(x<<1)+(c^48);
c=getchar();
}
return f*x;
}
void fft(std::complex <double> *f,int limit,int type){
for(int i=0;i<=limit;i++){
if(i<r[i]){
swap(f[i],f[r[i]]);
}
}
for(int length=2;length<=limit;length=(length<<1)){
std::complex <double> Wn (std::cos(Pi*2/length),std::sin(type*Pi*2/length));
for(int i=0;i<=limit;i+=length){
std::complex <double> wn1(1,0);
for(int k=i;k<length/2+i;k++,wn1=wn1*Wn){
std::complex <double> x=f[k],y=f[k+length/2]*wn1;
f[k]=x+y;
f[k+length/2]=x-y;
}
}
}
return ;
}
NTT (快速数论变换)
思想
因为FFT涉及到复杂的复数运算,不仅很慢,不能取模,甚至有精度问题。
所有如果在取模意义下,也可以用与FFT所利用的的单位根相似的原根。
数学知识
原根的定义理所当然,定义一个数 $n$ 的原根为 $g$ ,是满足$g^n \equiv 1(mod~p)$的最小整数。
原根满足 $g^i(1 \leq i \leq n)$ 模 $p$ 两两不同。
一般取模数为 $998244353$ ,其原根是 $3$。
对于代值的数的表示仍然是$\omega_n$:
$w_n \equiv g^{\frac{p-1}{n}} (mod~p)$
$\omega_n^k=\omega_{n/2}^{k/2}$ 。
$\omega_n^{k}=-\omega_n^{k+n/2}$。
$\omega_n^0=\omega_n^n=1$。
$\omega_n^{k}=\omega_n^{n+k}$。
直接把FFT的代值改成NTT
#include <iostream>
#include <cstdio>
#include <cmath>
#define int long long
inline int read();
int a[8000005],b[8000005];
int li=1,k=0;
int r[8000005];
const int mod=998244353;
const int g=3;
const int g2=332748118;
void print(int x){
register char s[20];
register int i=0;
if(x<0){
x=-x;
putchar('-');
}
if(x==0){
putchar('0');
return;
}
while(x){
s[i++]=x%10;
x/=10;
}
while(i){
putchar(s[--i]+'0');
}
return;
}
inline int fast_power(int a, int k) //快速幂,a为底数,k为指数
{
register int res = 1;
while (k)
{
if (k & 1)
res = res * a % mod;
a = a * a % mod;
k >>= 1;
}
return res;
}
inline void ntt(int *a,int type){
for(register int i=0;i<=li;i++){
if(i<r[i]){
std::swap(a[i],a[r[i]]);
}
}
int gn;
for(register int i=1;i<li;i<<=1){
gn=fast_power(type?g:g2,(mod-1)/(i<<1));
for(register int j=0;j<li;j+=(i<<1)){
int g0=1;
for(int k=0;k<i;k++,g0=(g0*gn)%mod){
int x=a[j+k];
int y=(a[i+j+k]*g0)%mod;
a[j+k]=(x+y)%mod;
a[i+j+k]=(x-y+mod)%mod;
}
}
}
return ;
}
signed main(){
#ifdef ONLINE_JUDGE
#else
// freopen(".in","r",stdin);
// freopen(".out","w",stdout);
#endif
int n=read(),m=read();
for(register int i=0;i<=n;i++){
a[i]=read();
}
for(register int i=0;i<=m;i++){
b[i]=read();
}
while(li<(n+m+2)){
li=(li)<<1;
k++;
}
for(register int i=0;i<=li;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
}
ntt(a,1);
ntt(b,1);
for(register int i=0;i<=li;i++){
a[i]*=b[i];
a[i]%=mod;
}
ntt(a,0);
long long inv=fast_power(li,mod-2);
for(register int i=0;i<=n+m;i++){
print((a[i]*inv)%mod);
putchar(' ');
}
return 0;
}
inline int read(){
int x=0,f=1;char c=getchar();
while(c<'0'||c>'9'){
c=='-'?f=-1:1;
c=getchar();
}
while(c>='0'&&c<='9'){
x=(x<<3)+(x<<1)+(c^48);
c=getchar();
}
return f*x;
}
/*
Anything about this program:
Type:
Description:
Example:
1:
In:
Out:
More:
*/