求多项式卷积的变换

• 多项式的卷积

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

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

  $$(f⨂g{)}_x=\sum _{i=0}^x{f}_i\cdot g_{x-i}$$

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

• 快速傅里叶变换（FFT）

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

           \ \ \ \ \ \ \,       这里还有 $3\mathtt{B}\mathtt{l}\mathtt{u}\mathtt{e}1\mathtt{B}\mathtt{r}\mathtt{o}\mathtt{w}\mathtt{n}$ 的视频，可以形象理解一下：【形象展示傅里叶变换】

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

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

  $$(f⨂g{)}_x=\sum _{i=1}^x{f}_i\cdot g_{x-i}$$

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

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$ 的多项式 $g$ ,和一个长度为 $n$ 的多项式 $f$，定义$f$为：

  $$f_x=\sum _{i=1}^x{f}_{x-i}\cdot g_i$$

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

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

  $$f_x=(f⨂g{)}_x=\sum _{i=1}^x{g}_i\cdot f_{x-i}$$

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

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

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

  $$f_x=\sum _{i=1}^x{g}_i\cdot f_{x-i}$$

  $$f_x=\sum _{i=1}^x{g}_{x-i}\cdot f_i$$

  $$f_x=\sum _{i=1}^x{g}_{x-i}\cdot A_i$$

  $$f_x=g⨂A$$

           \ \ \ \ \ \ \,       注意我们这里把$g$的长度只取前$r-l$个哦，因为$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);
  }
  

           \ \ \ \ \ \ \,       然而……要是你用这个就死了，主要是死在 $\mathtt{m}\mathtt{e}\mathtt{m}\mathtt{s}\mathtt{e}\mathtt{t}$了，所以我们把 $\mathtt{m}\mathtt{e}\mathtt{r}\mathtt{g}\mathtt{e}\_\mathtt{f}\mathtt{f}\mathtt{t}$ 托下来：

  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 }^{2}n)$更快的做法吗？当然有了啊，但是又要涉及生成函数的姿势，还有多项式求逆！！！我不会啊！！！

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

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

  $$F(x)=\sum _{i=0}^n{f}_i\cdot x^i$$

  $$G(x)=\sum _{i=0}^n{g}_i\cdot x^i$$

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

  $$\begin{array}{rl}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⨂g{)}_i\end{array}$$

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

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

  $$F(x)\times G(x)\equiv F(x)-f_0{x}^0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

  $$F(x)\times G(x)\equiv F(x)-f_0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

  $$F(x)\times \left(1-G(x)\right)\equiv f_0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

  $$F(x)\equiv \frac{f_0}{1-G(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

  $$F(x)\equiv \frac{1}{1-G(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

  $$F(x)\equiv {(1-G(x))}^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

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

---

• 快速数论变换（NTT）

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

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

  • 定义：

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

  • 性质：

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

    $2$、如果$(\mathrm{m}\mathrm{o}\mathrm{d}b)$意义下有原根，那么原根恰好有$\phi (\phi (b))$个，这个数量很小。

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

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

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

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

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

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

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

           \ \ \ \ \ \ \,       下面直接给板子了,这里的模数是$998244353$，原根是$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$的多项式$f$，和一个常数$a$，我们知道可以这样运算：

  $$\frac{f}{a}=\sum _{i=1}^n\frac{f_i}{a}$$

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

  $$\frac{g}{a}=\sum _{i=1}^n\frac{g_i}{a}$$

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

  $${\left(\frac{f}{a}⨂\frac{g}{a}\right)}_x=\sum _{i=1}^x\frac{f_i}{a}\cdot \frac{g_{x-i}}{a}$$

  $$\frac{f}{a}⨂\frac{g}{a}=\frac{f⨂g}{a^2}$$

  $$f⨂g=\left(\frac{f}{a}⨂\frac{g}{a}\right)\times a^2$$

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

  $$f⨂g=\left(\frac{f}{a}⨂\frac{g}{a}\right)\times a^2+\left(f\%a⨂\frac{g}{a}\right)\times a+\left(\frac{f}{a}⨂g\%a\right)\times a+\left(f\%a⨂g\%a\right)$$

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

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=\sum _{s\subseteq O}{f}_s\cdot x^s$$

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

  • 交 $\bigcap$ （与 $\mathrm{a}\mathrm{n}\mathrm{d}$）

    $$F\underset{\bigcap }{⨂}G=\sum _{i\subseteq O}\sum _{j\subseteq O}{f}_i\times g_j[i\bigcap j=S]$$

  • 并 $\bigcup$ （或 $\mathrm{o}\mathrm{r}$）

    $$F\underset{\bigcup }{⨂}G=\sum _{i\subseteq O}\sum _{j\subseteq O}{f}_i\times g_j[i\bigcup j=S]$$

  • 对称差 $\Delta$ （异或 $\mathrm{x}\mathrm{o}\mathrm{r}$）

    $$F\underset{\Delta }{⨂}G=\sum _{i\subseteq O}\sum _{j\subseteq O}{f}_i\Delta g_j[i\Delta =S]$$

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

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

运算方式	交 $\bigcap$ （与 $\mathrm{a}\mathrm{n}\mathrm{d}$）	并 $\bigcup$ （或 $\mathrm{o}\mathrm{r}$）	对称差 $\Delta$ （异或 $\mathrm{x}\mathrm{o}\mathrm{r}$）
沃尔什变换	$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);
}