求多项式卷积的变换

$\ \ \ \ \ \ \ \,$关于导数和牛顿迭代的复习笔记：

多项式的卷积

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

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

$(f \bigotimes g)_x = \sum_{i=0}^{x}f_i\cdot g_{x-i}$

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

快速傅里叶变换（FFT）

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

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

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

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

$(f \bigotimes 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}^xf_{x-i}\cdot g_i$

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

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

$f_x=(f \bigotimes 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 \bigotimes 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);
}

$\ \ \ \ \ \ \,$然而……要是你用这个就死了，主要是死在 $\tt memset$了，所以我们把 $\tt 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)$更快的做法吗？当然有了啊，但是又要涉及生成函数的姿势，还有多项式求逆！！！我不会啊！！！

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

$\ \ \ \ \ \ \,$我们令$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{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_{mid+1,r}= f_{l,mid} \bigotimes g$，当然了，$f_0$是不满足这个卷积的，所以有：

$F(x)\times G(x)\equiv F(x)-f_0x^0\ \ (mod \ x^n)$

$F(x)\times G(x)\equiv F(x)-f_0\ \ (mod \ x^n)$

$F(x)\times \left(1-G(x)\right)\equiv f_0\ \ (mod \ x^n)$

$F(x)\equiv \frac{f_0}{1-G(x)}\ \ (mod \ x^n)$

$F(x)\equiv \frac{1}{1-G(x)}\ \ (mod \ x^n)$

$F(x)\equiv {(1-G(x))}^{-1}\ \ (mod \ x^n)$

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

---

快速数论变换（NTT）

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

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

• 定义：

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

• 性质：

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

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

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

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

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

  $\ \ \ \ \ \ \,$考虑FFT中的单位根$x^n=1$，现在在取模的意义下，就变成了$x^n\equiv 1\ \ (mod \ 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，你可以用这个方法过洛谷的模板，但是写起来非常复杂。下面还有一种方法:

MTT

Mettaton

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

$\ \ \ \ \ \ \,$（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} \bigotimes \frac{g}{a}\right)_x = \sum_{i=1}^{x}\frac{f_i}{a}\cdot \frac{g_{x-i}}{a}$

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

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

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

${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)$

$\ \ \ \ \ \ \,$就差不多这样了，取值$a$的时候需要保证用FFT时不会爆炸，不过也容易看出，这个算法精度特别特别卡，卡死了，需要开 $long\ double$ 还有一些特别的奇技淫巧，下面给出模板，$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$ （与 $and$）

  $F\bigotimes_{\bigcap} G=\sum_{i\subseteq O}\sum_{j\subseteq O}f_i\times g_j[i\bigcap j=S]$

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

  $F\bigotimes_{\bigcup} G=\sum_{i\subseteq O}\sum_{j\subseteq O}f_i\times g_j[i\bigcup j=S]$

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

  $F\bigotimes_{\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$ （与 $and$）	并 $\bigcup$ （或 $or$）	对称差 $\Delta$ （异或 $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);
}