求多项式卷积的变换

• 多项式的卷积

  $\ \ \ \ \ \ \,$      在泛函分析中，卷积、旋积或摺积(英语：Convolution)是通过两个函数f 和g 生成第三个函数的一种数学算子，表征函数f 与g经过翻转和平移的重叠部分的面积。

  $\ \ \ \ \ \ \,$      对于一个长度为 $n$n 的多项式 $f$f ,和一个长度为 $m$m 的多项式 $g$g，那么他们的卷积就可以表达为：

  $(f \bigotimes g)_x = \sum_{i=0}^{x}f_i\cdot g_{x-i}$(f⨂g)x​=i=0∑x​fi​⋅gx−i​

  $\ \ \ \ \ \ \,$      其中暴力求的复杂度会是$O(nm)$O(nm)的，在很多情况下是不允许有这么多复杂度的，下面我们来引入一些快速变换：

• 快速傅里叶变换（FFT）

  $\ \ \ \ \ \ \,$      这个算法算是非常常见的对卷积的处理了，讲起来也非常繁琐，留个链接【傅里叶变换(FFT)学习笔记】 【VictoryCzt_FFT学习笔记】。

  $\ \ \ \ \ \ \,$      这里还有 $\tt 3Blue1Brown$3Blue1Brown 的视频，可以形象理解一下：【形象展示傅里叶变换】

  $\ \ \ \ \ \ \,$      然后就贴板子吧：

  $\ \ \ \ \ \ \,$      依然是对于一个长度为 $n$n 的多项式 $f$f ,和一个长度为 $m$m 的多项式 $g$g，那么他们的卷积就可以表达为：

  $(f \bigotimes g)_x = \sum_{i=1}^{x}f_i\cdot g_{x-i}$(f⨂g)x​=i=1∑x​fi​⋅gx−i​

  $\ \ \ \ \ \ \,$      使用快速傅里叶变换（FFT）可以达到$O(n\log n)$O(nlogn)的优秀复杂度。

struct cpx{
	double r,i;
	inline cpx operator *(const cpx&x)const{return (cpx){r*x.r-i*x.i,r*x.i+i*x.r};}
	inline cpx operator +(const cpx&x)const{return (cpx){r+x.r,i+x.i};}
	inline cpx operator -(const cpx&x)const{return (cpx){r-x.r,i-x.i};}
}a[N],b[N];
int m,n,R[N];
void FFT(cpx*a,int f,int la){
	int n=la;
	for(register int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);
	for(register int i=1;i<n;i<<=1){
		cpx wn=(cpx){cos(pi/i),f*sin(pi/i)};
		for(register int j=0;j<n;j+=(i<<1)){
			cpx w=(cpx){1,0};
			for(register int k=0;k<i;++k,w=w*wn){
				cpx x=a[j+k],y=w*a[j+k+i];
				a[j+k]=x+y;a[j+k+i]=x-y;
			}
		}
	}
	if(f==-1)
	for(register int i=0;i<n;i++)a[i].r/=n;
}
int merge_fft(cpx *a,cpx *b,int la,int lb){
	int n=la,m=lb;
	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
	for(register int i=0;i<n;i++)
	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	FFT(a,1,n);FFT(b,1,n);
	for(register int i=0;i<=n;i++)a[i]=a[i]*b[i];
	FFT(a,-1,n);
	return m;
}

• 分治FFT

  $\ \ \ \ \ \ \,$      对于一个长度为 $n$n 的多项式 $g$g ,和一个长度为 $n$n 的多项式 $f$f，定义$f$f为：

  $f_x=\sum_{i=1}^xf_{x-i}\cdot g_i$fx​=i=1∑x​fx−i​⋅gi​

  $\ \ \ \ \ \ \,$      其中，$f_0=1$f0​=1。

  $\ \ \ \ \ \ \,$      很容易看出这个式子是卷积形式:

  $f_x=(f \bigotimes g)_x = \sum_{i=1}^{x}g_i\cdot f_{x-i}$fx​=(f⨂g)x​=i=1∑x​gi​⋅fx−i​

  $\ \ \ \ \ \ \,$      我们想到上面提到的FFT算法，但是这里是$f$f函数有卷之前算到的自己，所以我们第一次$f$f的有效范围是$1$1，卷一次，可以得到$f$f有效范围是$2$2，，再一次是$4$4，就像倍增一样，所以我们的复杂度就应该是$O(n\log ^2 n)$O(nlog2n)。

  $\ \ \ \ \ \ \,$      具体怎么做？

  $\ \ \ \ \ \ \,$      假设对于$x$x的取值区间$l$l到$r$r，我们知道了$l$l到 $\frac{l+r}{2}$2l+r​ 的$f_x$fx​的值，那么我们令一个长度为 $\frac{r-l}{2}$2r−l​ 的多项式 $A$A 等于这知道的 $\frac{r-l}{2}$2r−l​ 项，我们现在想知道前面这一段对于后面 $\frac{l+r}{2}$2l+r​ 到 $r$r 这一段的贡献是多少：

  $f_x= \sum_{i=1}^{x}g_i\cdot f_{x-i}$fx​=i=1∑x​gi​⋅fx−i​

  $f_x= \sum_{i=1}^{x}g_{x-i}\cdot f_i$fx​=i=1∑x​gx−i​⋅fi​

  $f_x= \sum_{i=1}^{x}g_{x-i}\cdot A_i$fx​=i=1∑x​gx−i​⋅Ai​

  $f_x= g \bigotimes A$fx​=g⨂A

  $\ \ \ \ \ \ \,$      注意我们这里把$g$g的长度只取前$r-l$r−l个哦，因为$x$x的范围限制，我们需要保证卷出来的长度是我们需要的。然后贴板子：

  void cdq_FFT(int l,int r){
    	if(l==r)return;
    	int mid=(l+r)>>1;
    	cdq_FFT(l,mid);
    memset(A,0,sizeof(A));memset(G,0,sizeof(G));
    	for(int i=l;i<=mid;i++)A[i-l].r=(double)f[i];
    	for(int i=1;i<=r-l;i++)G[i-1].r=(double)g[i];
    	merge_fft(A,G,mid-l,r-l-1);
    	for(int i=mid+1;i<=r;i++)f[i]+=(int)A[i-l-1].r;
    	cdq_FFT(mid+1,r);
  }
  

  $\ \ \ \ \ \ \,$      然而……要是你用这个就死了，主要是死在 $\tt memset$memset了，所以我们把 $\tt merge\_fft$merge_fft 托下来：

  void cdq_FFT(int l,int r){
    	if(l==r)return;
    	int mid=(l+r)>>1;
    	cdq_FFT(l,mid);
    	int L=0,n=1;
    	for(;n<=r-l+1;n<<=1)L++;
    	for(int i=0;i<n;i++)
    	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    	for(int i=0;i<=n;i++)
    	  A[i].r=G[i].r=A[i].i=G[i].i=0
    	for(int i=l;i<=mid;i++)A[i-l].r=(double)f[i];
    	for(int i=1;i<=r-l;i++)G[i-1].r=(double)g[i];
    	FFT(A,1,n);FFT(G,1,n);
    	for(register int i=0;i<=n;i++)A[i]=A[i]*G[i];
    	FFT(A,-1,n);
    	for(int i=mid+1;i<=r;i++)f[i]+=(int)A[i-l-1].r;
    	cdq_FFT(mid+1,r);
  }
  

  $\ \ \ \ \ \ \,$      至于洛谷那道模板题P4721还是算了吧，要取模，建议把FFT换成NTT，（NTT模数明示）

  $\ \ \ \ \ \ \,$      那么还有比$O(n\log ^2n)$O(nlog2n)更快的做法吗？当然有了啊，但是又要涉及生成函数的姿势，还有多项式求逆！！！我不会啊！！！

  $\ \ \ \ \ \ \,$      先占坑:

  $\ \ \ \ \ \ \,$      我们令$F$F为$f$f的生成函数，$G$G为$g$g的生成函数：

  $F(x)=\sum_{i=0}^{n}f_i\cdot x^i$F(x)=i=0∑n​fi​⋅xi

  $G(x)=\sum_{i=0}^{n}g_i\cdot x^i$G(x)=i=0∑n​gi​⋅xi

  $\ \ \ \ \ \ \,$      那么就有：

  $\begin{aligned}F(x)\times G(x) &= \left(\sum_{i=0}^{n}f_i\cdot x^i\right)\cdot \left(\sum_{i=0}^{n}g_i\cdot x^i\right)\\&=\sum_{i=0}^{n}x^i\sum_{j=0}^{i}g_j\cdot f_{i-j}\\&=\sum_{i=0}^{n}x^i(f \bigotimes g)_i\end{aligned}$F(x)×G(x)​=(i=0∑n​fi​⋅xi)⋅(i=0∑n​gi​⋅xi)=i=0∑n​xij=0∑i​gj​⋅fi−j​=i=0∑n​xi(f⨂g)i​​

  $\ \ \ \ \ \ \,$      也就是说两个函数的生成函数的积，就是他们卷积的生成函数~~@贝尔级数~~。

  $\ \ \ \ \ \ \,$      我们知道的$f_{mid+1,r}= f_{l,mid} \bigotimes g$fmid+1,r​=fl,mid​⨂g，当然了，$f_0$f0​是不满足这个卷积的，所以有：

  $F(x)\times G(x)\equiv F(x)-f_0x^0\ \ ({\rm mod}\ x^n)$F(x)×G(x)≡F(x)−f0​x0  (mod xn)

  $F(x)\times G(x)\equiv F(x)-f_0\ \ ({\rm mod}\ x^n)$F(x)×G(x)≡F(x)−f0​  (mod xn)

  $F(x)\times \left(1-G(x)\right)\equiv f_0\ \ ({\rm mod}\ x^n)$F(x)×(1−G(x))≡f0​  (mod xn)

  $F(x)\equiv \frac{f_0}{1-G(x)}\ \ ({\rm mod}\ x^n)$F(x)≡1−G(x)f0​​  (mod xn)

  $F(x)\equiv \frac{1}{1-G(x)}\ \ ({\rm mod}\ x^n)$F(x)≡1−G(x)1​  (mod xn)

  $F(x)\equiv {(1-G(x))}^{-1}\ \ ({\rm mod}\ x^n)$F(x)≡(1−G(x))−1  (mod xn)

  $\ \ \ \ \ \ \,$      然后多项式求逆啊_，并不会~，复杂度$O(n\log n)$O(nlogn)，就不贴代码了啊。

---

• 快速数论变换（NTT）

  $\ \ \ \ \ \ \,$      由于FFT是在复平面上操作变换，会用到很多的三角函数，丢精度非常严重，当有模数的时候且模数为NTT模数，形如$a\cdot 2^b+1$a⋅2b+1，且为素数的时候，常见的是$998244353$998244353，我们推荐使用NTT来操作。

  - 【原根】%%%大巨佬orz

  • 定义：

    $\ \ \ \ \ \ \,$      考虑方程$a^{x}\equiv 1({\rm mod}\ b)$ax≡1(mod b)，根据欧拉定理，当$a$a与$b$b互质的时候，一定有解$x=\varphi(b)$x=φ(b)，如果$\varphi(b)$φ(b)既为该方程的最小解，那么$a$a就称为$({\rm mod}\ b)$(mod b)的原根。

  • 性质：

    $\ \ \ \ \ \ \,1$      1、若$a$a就为$({\rm mod}\ b)$(mod b)的原根，那么集合$\{a^1,a^2,...,a^{b-1}\}=\{0,1,2,...,b-1\}${a1,a2,...,ab−1}={0,1,2,...,b−1}。也就是说，$a$a的任意次方遍历了$({\rm mod}\ b)$(mod b)的剩余系。

    $\ \ \ \ \ \ \,2$      2、如果$({\rm mod}\ b)$(mod b)意义下有原根，那么原根恰好有$\varphi(\varphi(b))$φ(φ(b))个，这个数量很小。

    $\ \ \ \ \ \ \,3$      3、原根一般都不大，我们可以暴力从$2$2到$b-1$b−1枚举$a$a，判断$a^{\varphi(b)}\equiv 1({\rm mod}\ b)$aφ(b)≡1(mod b)是否当且仅当$x=\varphi(b)$x=φ(b)时成立，注意$x$x的解仅有可能是$\varphi(b)$φ(b)的约数，还有快的求法，可以用欧拉定理和裴蜀定理证明，可是我不会。

  $\ \ \ \ \ \ \,$      当然了，对于原根的理解我们只需要理解到它是$a^{x}\equiv 1({\rm mod}\ b)$ax≡1(mod b)的最小解，而且可以暴力求就行了2333就行了。

  $\ \ \ \ \ \ \,$      那么如何用到快速变换中去呢？

  $\ \ \ \ \ \ \,$      考虑FFT中的单位根$x^n=1$xn=1，现在在取模的意义下，就变成了$x^n\equiv 1\ \ ({\rm mod}\ P)$xn≡1  (mod P)，所以我们用原根$gn$gn去替换掉FFT中的单位根$wn$wn就好了呀，现在的问题是，原根具体怎么求：

  $\ \ \ \ \ \ \,$      对于一个NTT模数的原根，是有特殊求法的，这就是为什么一般NTT会有NTT模数这个限制：

  $gn=g^{\frac{P-1}{2^{i+1}}}$gn=g2i+1P−1​

  $\ \ \ \ \ \ \,$      原根表在这里：【FFT用到的各种素数】

  $\ \ \ \ \ \ \,$      下面直接给板子了,这里的模数是$998244353$998244353，原根是$3$3：

void NTT(int *a,int f,int la){
	int n=la;
	for(register int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);
	for(register int i=1;i<n;i<<=1){
		int gn=power(3,(mod-1)/(i<<1));
		for(register int j=0;j<n;j+=(i<<1)){
			int g=1;
			for(register int k=0;k<i;++k,g=1ll*g*gn%mod){
				int x=a[j+k],y=1ll*g*a[j+k+i]%mod;
				a[j+k]=(x+y)%mod;a[j+k+i]=(x-y+mod)%mod;
			}
		}
	}
	if(f==-1){
		reverse(a+1,a+n);
		int inv=power(n,mod-2);
		for(register int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
	}
}
int merge_ntt(int *a,int *b,int la,int lb){
	int n=la,m=lb;
	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
	for(register int i=0;i<n;i++)
	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	NTT(a,1,n);NTT(b,1,n);
	for(register int i=0;i<=n;i++)a[i]=1ll*a[i]*b[i]%mod;
	NTT(a,-1,n);
	return m;
}

$\ \ \ \ \ \ \,$      那么不是NTT模数的时候怎么办啊？我们不妨用CRT合并，具体在这里可以看到的【由CRT合并引入的各种扩展求余和同余问题】。这个方法叫三模NTT，你可以用这个方法过洛谷的模板，但是写起来非常复杂。下面还有一种方法:

• Mettaton

  

  $\ \ \ \ \ \ \,$      OH，YEAH~

  $\ \ \ \ \ \ \,$      （MTT就是镁塔顿啊，雾）

  $\ \ \ \ \ \ \,$      好吧让我们重新打开一下：

• MTT

  $\ \ \ \ \ \ \,$      对于一个长度为$n$n的多项式$f$f，和一个常数$a$a，我们知道可以这样运算：

  $\frac{f}{a}=\sum_{i=1}^{n}\frac{f_i}{a}$af​=i=1∑n​afi​​

  $\ \ \ \ \ \ \,$      同理，对于一个长度为$n$n的多项式$g$g:

  $\frac{g}{a}=\sum_{i=1}^{n}\frac{g_i}{a}$ag​=i=1∑n​agi​​

  $\ \ \ \ \ \ \,$      那么有：

  $\left(\frac{f}{a} \bigotimes \frac{g}{a}\right)_x = \sum_{i=1}^{x}\frac{f_i}{a}\cdot \frac{g_{x-i}}{a}$(af​⨂ag​)x​=i=1∑x​afi​​⋅agx−i​​

  $\frac{f}{a} \bigotimes \frac{g}{a}= \frac{{f}\bigotimes {g}}{a^2}$af​⨂ag​=a2f⨂g​

  ${{f}\bigotimes {g}}= \left(\frac{f}{a} \bigotimes \frac{g}{a}\right)\times a^2$f⨂g=(af​⨂ag​)×a2

  $\ \ \ \ \ \ \,$      所以我们想，要是我们把每一项的系数拆了，最后合并不就是我们想要的答案了吗？于是就有了上面那个这个式子，但是很多时候并不能整除啊，就有了下面这个式子:

  ${{f}\bigotimes {g}}=\left(\frac{f}{a} \bigotimes \frac{g}{a}\right)\times a^2+\left({f}\%{a} \bigotimes \frac{g}{a}\right)\times a+\left(\frac{f}{a} \bigotimes {g}\%{a}\right)\times a+\left({f}\%{a} \bigotimes {g}\%{a}\right)$f⨂g=(af​⨂ag​)×a2+(f%a⨂ag​)×a+(af​⨂g%a)×a+(f%a⨂g%a)

  $\ \ \ \ \ \ \,$      就差不多这样了，取值$a$a的时候需要保证用FFT时不会爆炸，不过也容易看出，这个算法精度特别特别卡，卡死了，需要开 $\rm long\ double$long double 还有一些特别的奇技淫巧，下面给出模板，$a=2^{15}$a=215，这个板子因为精度问题是过不了洛谷的模板题的，看看就好：

void merge(cpx *a,cpx *b,int L,int len,int val,int *Ret){
	for(int i=0;i<=L;i++)f[i]=a[i]*b[i];
	FFT(f,-1,L);
	for(int i=0;i<len;i++)
	Ret[i]=(Ret[i]+1ll*(long long)(f[i].r+0.5)%mod*val%mod)%mod;
}
void merge_mtt(int *A,int *B,int la,int lb,int *Ret){
	int M=32768,L,l=0,len=la;
	for(int i=0;i<(len<<1);++i)a[i]=b[i]=c[i]=d[i]=(cpx){0,0};
	for(int i=0;i<=len;i++)a[i].r=A[i]/M,b[i].r=A[i]%M;
	for(int i=0;i<=len;i++)c[i].r=B[i]/M,d[i].r=B[i]%M;
	for(L=1;L<=len;L<<=1)++l;
	for(int i=0;i<L;++i)R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
	FFT(a,1,L);FFT(b,1,L);FFT(c,1,L);FFT(d,1,L);
	for(int i=0;i<len;++i)Ret[i]=0;
	merge(a,c,L,len,M*M%mod,Ret);
	merge(b,c,L,len,M%mod,Ret);
	merge(a,d,L,len,M%mod,Ret);
	merge(b,d,L,len,1,Ret);
	for(int i=0;i<len;++i)Ret[i]=(Ret[i]+mod)%mod;
}

---

• 快速沃尔什变换（FWT）

  $\ \ \ \ \ \ \,$      首先引入一种生成函数集合幂级数$F$F，其中$f$f 是集合中的元素：

  $F=\sum_{s\subseteq O}f_s\cdot x^s$F=s⊆O∑​fs​⋅xs

  $\ \ \ \ \ \ \,$      集合幂级数的集合卷积：

  • 交 $\bigcap$⋂ （与 $\rm and$and）

    $F\bigotimes_{\bigcap} G=\sum_{{i\subseteq O}}\sum_{{j\subseteq O}}f_i\times g_j[i\bigcap j=S]$F⋂⨂​G=i⊆O∑​j⊆O∑​fi​×gj​[i⋂j=S]

  • 并 $\bigcup$⋃ （或 $\rm or$or）

    $F\bigotimes_{\bigcup} G=\sum_{{i\subseteq O}}\sum_{{j\subseteq O}}f_i\times g_j[i\bigcup j=S]$F⋃⨂​G=i⊆O∑​j⊆O∑​fi​×gj​[i⋃j=S]

  • 对称差 $\Delta$Δ （异或 $\rm xor$xor）

    $F\bigotimes_{\Delta} G=\sum_{{i\subseteq O}}\sum_{{j\subseteq O}}f_i\Delta g_j[i\Delta=S]$FΔ⨂​G=i⊆O∑​j⊆O∑​fi​Δgj​[iΔ=S]

  $\ \ \ \ \ \ \,$      而快速沃尔什变换（FWT）就是快速求出集合幂级数的集合卷积的算法，证明在这里看【FWT(Fast Walsh Transformation)快速沃尔什变换学习笔记】，下面给出结论和模板代码：

  $\ \ \ \ \ \ \,$      其中$x=a_{i+j}$x=ai+j​，$y=a_{i+j+bit}$y=ai+j+bit​:

运算方式	交 $\bigcap$⋂ （与 $\rm and$and）	并 $\bigcup$⋃ （或 $\rm or$or）	对称差 $\Delta$Δ （异或 $\rm xor$xor）
沃尔什变换	$\ \ a_{i+j}=x+y\ \ $	$\ \ a_{i+j+bit}=x+y\ \ $	$\ \ a_{i+j}=x+y\ \ $,$,\ \ a_{i+j+bit}=x-y\ \ $
沃尔什逆变换	$\ \ a_{i+j}=x-y\ \ $	$\ \ a_{i+j+bit}=y-x\ \ $	$\ \ a_{i+j}=\frac{x+y}{2}\ \ $,$,\ \ a_{i+j+bit}=\frac{x-y}{2}\ \ $

void FWT(int *a,int n,int f,int type){
  for(int i=1;i<n;i<<=1)
  for(int j=0;j<n;j+=(i<<1))
  for(int k=0;k<i;k++){
    int p=a[j+k],q=a[j+k+i];
    if(type==1){//and &
      a[j+k]=(1ll*f*q+p+mod)%mod;
    }
    if(type==2){//or |
      a[j+k+i]=(1ll*f*p+q+mod)%mod;
    }
    if(type==3){//xor ^
      a[j+k]=(p+q)%mod;
      a[j+k+i]=(p-q+mod)%mod;
      if(f==-1){
        a[j+k]=1ll*a[j+k]*inv2%mod;
        a[j+k+i]=1ll*a[j+k+i]*inv2%mod;
      }
    }
  }
}
void merge_fwt(int *a,int *b,int lim,int type){
  FWT(a,lim,1,type);FWT(b,lim,1,type);
  for(int i=0;i<lim;i++)a[i]=(long long)a[i]*b[i]%mod;
  FWT(a,lim,-1,type);
}