lp5395 【模板】第二类斯特林数·行

第二类斯特林数,指的是一组表示「将n个不同的元素划分为m个非空不相交集的方案数」的组合数。有时写作\(S(n,m)\)或\(\left\{\begin{matrix}n\\m\end{matrix}\right\}\)
从题意来看,我们可以容易地想到一个递推方案:每个新的元素,可以为它新开一个集合,或者放到已有的任何一个集合里面。所以我们得到一个递推式:
$$S(n,m)=S(n-1,m-1)+m*S(n-1,m)$$
然后问题是怎么求值。

首先我们有一个式子,我们称之为二项式反演的通项式。
$$ f(n)=\sum_{i=0}^nC_{n}^{i}g(n) \Leftrightarrow g(n)=\sum_{i=0}^n(-1)^{n-i}C_{i}{n}f(n)$$
举个例子:
$$\begin{matrix} f_{1}&=&g_{1}\\f_{2}&=&g_{1}&+&g_{2}\\f_{3}&=&g_{1}&+&2*g_{2}&+&g_{3}\\f_{4}&=&g_{1}&+&3*g_{2}&+&3*g_{3}&+&g_{4} \end{matrix}$$

由此可得:

$$\begin{matrix} g_{1}&=&f_{1}\\g_{2}&=&f_{2}&-&f_{1}\\g_{3}&=&f_{3}&-&2*g_{2}&-&g_{1}\\&=&f_{3}&-&2*f_{2}&+&f_{1}\\g_{4}&=&f_{4}&-&3*g_{3}&-&3*g_{2}&-&g_{1}\\&=&f_{4}&-&3*f_{3}&-&3*g_{2}&-&4*g_{1}\\&=&f_{4}&-&3*f_{3}&+&3*f_{2}&-&f_{1} \end{matrix}$$

然后我们考虑\(m^n\)的组合意义,也就是,将\(n\)个不同的球,放到\(m\)个不同的盒子(可以有空盒)的方案数。
我们不妨枚举有装东西的盒子的个数。令它为\(i\),那么选取这些盒子的方案数即是\(C_{m}^{i}\)。
选取了盒子之后,问题就转化为了将\(n\)个不同物品装到\(m\)个不同盒子里,使得每个盒子非空。这就相当于,将\(n\)个不同物品装到\(m\)个相同盒子里,使得每个盒子非空的方式——也就是\(S(n,m)\),乘上盒子的排列顺序——也就是\(m!\)。
形式化的说,就是:
$$m^n=\sum_{i=0}^{m}S(n,i)i!C_{m}^{i}$$
我们发现这个式子的形式符合二项式反演的通项式。
因而我们将\(m^n\)作为\(f\),将\(S(n,i)i!\)作为\(g\)。
那么反演得到:
$$S(n,m)=\frac{1}{m!}\sum_{i=0}^{m}(-1)^{m-i}i^nC_{m}^{i}$$
考虑\(C_{m}^{i}=\frac{m!}{i!(m-i)!}\),我们可以将上式化简得到:
$$S(n,m)=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}$$
于是我们有:
$$\left\{\begin{matrix}n\\m\end{matrix}\right\}=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}$$
考虑函数卷积的本质,对于卷积\(f=gu\),其本质是: $$f(n)=\sum_{i=0}^{n}g(i)u(n-i)$$ 因而,我们发现,上式本质上就是: $$\left\{\begin{matrix}n\\m\end{matrix}\right\}=\sum_{i=0}^{m}\frac{(-1)^{m-i}}{(m-i)!}*\frac{i^n}{i!}$$
上一个NTT卷积一下就好了。

#include<iostream>
#include<cstdio>

#define Swap(A,B) (A^=B^=A^=B)
const long long P=167772161;
const long long g0=3,gi=55924054;
int L=1,R[1<<21|1];
long long invL;
inline int pw(int A,int X){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT=RT*BS%P;
		}
		BS=BS*BS%P;
		X>>=1;
	}
	return RT;
}
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	invL=P-(P-1)/L;
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}
inline void FNTT(int *A,int typ){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[R[i]],A[i]);
		}
	}
	int gn,g,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		gn=pw(~typ?g0:gi,(P-1)/i);
		for(int j=0;j<L;j+=i){
			g=1;
			for(int k=0;k<M;++k,g=1ll*g*gn%P){
				X=A[j+k],Y=1ll*g*A[M+j+k]%P;
				A[j+k]=(X+Y)%P,A[M+j+k]=(X-Y)%P;
			}
		}
	}
}
int n,a[1<<21|1],b[1<<21|1],inv[1<<21|1];
void init(){
	scanf("%d",&n);
	prpr(n+1<<1);
	inv[0]=inv[1]=1;
	for(int i=2;i<=n;++i){
		inv[i]=1ll*(P-P/i)*inv[P%i]%P;
	}
	for(int i=1;i<=n;++i){
		inv[i]=1ll*inv[i-1]*inv[i]%P;
	}
	for(int i=0;i<=n;++i){
		a[i]=1ll*pw(-1,i)*inv[i]%P;
		b[i]=1ll*pw(i,n)*inv[i]%P;
	}
	FNTT(a,1);
	FNTT(b,1);
	for(int i=0;i<L;++i){
		a[i]=1ll*a[i]*b[i]%P;
	}
	FNTT(a,-1);
	for(int i=0;i<=n;++i){
		printf("%d ",1ll*(a[i]+P)*invL%P);
	}
}

int main(){
	init();
	return 0;
}

lp4245 【模板】任意模数NTT

我们首先考虑一个很暴力的玩法,直接找一个很大很大很大的模数,然后用int128上一个朴素NTT,再将结果对题目给的模数取模,似乎就可以了。
然而仔细考虑一下发现这个做法实在太暴力了(才不是因为找不到合适的模数呢。)我们考虑另一种方案。
现在假设我们已经用好几种模数求出最终的各种值了,那么原值是多少呢?很显然可以用中国剩余定理来求出。
我们计算了一下,发现只要用三个质数,值域就足以用中国剩余定理求出最终值了。
然而我们知道,模运算没有交换律。故而我们不能用常规的中国剩余定理合并,应该需要使用一种类似于扩展中国剩余定理的递推合并方法来合并它们的解。
这就完成了任意模数的NTT。
另外要注意数组大小…

#include<iostream>
#include<cstdio>

inline void Swap(int &A,int &B){
	A^=B^=A^=B;
}

inline long long mlt(long long A,long long X,long long P){
	long long BS=A,RT=0;
	while(X){
		if(X&1){
			RT+=BS;
			RT%=P;
		}
		BS+=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

inline long long pw(long long A,long long X,long long P){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT*=BS;
			RT%=P;
		}
		BS*=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

const long long MOD[3]={469762049,998244353,1004535809};
const long long MM=468937312667959297;
const long long g0=3,gi[3]={156587350,332748118,334845270};

int L=1,inv[3],R[1<<21|1];
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	for(int i=0;i<3;++i){
		inv[i]=MOD[i]-(MOD[i]-1)/L;
	}
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}

inline void NTT(int *A,int typ,int P){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[i],A[R[i]]);
		}
	}
	int bs,nw,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		bs=pw(~typ?g0:gi[P],(MOD[P]-1)/i,MOD[P]);
		for(int j=0;j<L;j+=i){
			nw=1;
			for(int k=0;k<M;++k,nw=1ll*nw*bs%MOD[P]){
				X=A[j+k],Y=1ll*nw*A[j+k+M]%MOD[P];
				A[j+k]=(X+Y)%MOD[P];
				A[j+k+M]=(X-Y+MOD[P])%MOD[P];
			}
		}
	}
}
int C[1<<21|1],D[1<<21|1],ans[3][1<<21|1];
inline void FNTT(int *A,int *B,int P){
	for(int i=0;i<L;++i){
		C[i]=A[i],D[i]=B[i];
	}
	NTT(C,1,P);
	NTT(D,1,P);
	for(int i=0;i<L;++i){
		C[i]=1ll*C[i]*D[i]%MOD[P];
	}
	NTT(C,-1,P);
	for(int i=0;i<L;++i){
		ans[P][i]=(int)(1ll*(C[i]+MOD[P])*inv[P]%MOD[P]);
	}
}
int a[1<<21|1],b[1<<21|1],n,m,p;
long long t[3];
void init(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=0;i<=n;++i){
		scanf("%d",a+i);
	}
	for(int i=0;i<=m;++i){
		scanf("%d",b+i);
	}
	prpr(n+m);
	FNTT(a,b,0);FNTT(a,b,1);FNTT(a,b,2);
	t[0]=pw(MOD[1]%MOD[0],MOD[0]-2,MOD[0]);t[1]=pw(MOD[0]%MOD[1],MOD[1]-2,MOD[1]);t[2]=pw(MM%MOD[2],MOD[2]-2,MOD[2]);//分别求出要用到的三个逆元。 
	long long T1,T2;
	for(int i=0;i<=n+m;++i){
		T1=(mlt(1ll*ans[0][i]*MOD[1]%MM,t[0],MM)+mlt(1ll*ans[1][i]*MOD[0]%MM,t[1],MM))%MM;//直接用CRT合并一二式。 
		T2=((ans[2][i]-T1)%MOD[2]+MOD[2])%MOD[2]*t[2]%MOD[2];//用EXCRT将第三式合并之。 
		printf("%d ",(int)(((T2%p)*(MM%p)%p+T1%p)%p));//求得最终解。 
	}
}

int main(){
	init();
	return 0;
}

lp3803 【模板】多项式乘法(FFT)(NTT实现)

程序运行效率。(评测系统为duck.ac)

我们现在已经有一种快速完成多项式卷积的方案了。然而这种方案还存在两个缺陷。
首先,是浮点型运算带来的精度误差。这在多次乘积之后会不断扩大,最终到一个无法解决的地步。
第二,就是结构体运算带来的大常数。这也使得这道题的时间限制变得比较紧张。
于是,我们考虑另辟蹊径,使用某一种大模数来完成整数多项式卷积。
快速数论变换(Fast Number Transform)就是一种这样的算法。

回到FFT的本质,FFT是一种精确选取取值以后代值优化系数表达和点值表达相互转化的方案。我们尝试在模意义下找到一些满足与单位复根性质类似的数。
一种在数论上被称为原根的数可以满足这一族性质。
首先我们有原根的定义。
形式化的说,若\(g\)模\(p\)的阶等于\(\phi(p)\),则称\(g\)为\(p\)的一个原根。
具体来说,我们认为,若\(g^x \equiv 1\ (mod\ p) \)总是有解,那么我们称\(g\)是质数\(p\)的一个原根。
对于\(\forall g\in Z,p\in P\),总是有\(g^{p-1}\equiv 1\ (mod\ p) \)
那么,对于一个\(p=t2^{k}+1\)我们有,\(g^{\frac{p-1}{2^i}}\ (0\le i\le k)\)在模$latexp$意义下始终有正整数解。这一点和单位复根类似。
通过选取满足上述条件的质数,即\(p=t2^{k}+1\),我们还有这样的性质:
消去引理:
$$g^{\frac{dk(p-1)}{dn}}\equiv g^{\frac{k(p-1)}{n}}$$
折半引理:
$$(g^{\frac{(k+\frac{n}{2})(p-1)}{n}})^2\equiv (g^{\frac{k(p-1)}{n}}g^{\frac{p-1}{2}})^2\equiv (g^{\frac{k(p-1)}{n}}(p-1))^2\equiv (g^{\frac{k(p-1)}{n}})^2$$
求和引理:
$$\sum_{j=0}^{n-1}(g^{\frac{k(p-1)}{n}})^j\equiv \frac{(g^{\frac{k(p-1)}{n}})^n-1}{g^{\frac{k(p-1)}{n}}-1}\equiv \frac{g^{k(p-1)}-1}{g^{\frac{k(p-1)}{n}}-1}\equiv \frac{1^k-1}{g^{\frac{k(p-1)}{n}}}\equiv 0$$
折半引理的一个推论:
$$\sum_{i=0}^{n-1}g^{\frac{i(p-1)}{n}}\equiv 0$$
我们惊讶地发现原根具有单位复根具有的三个定理和一个推论。这说明,我们可以将快速傅里叶变换的操作以原根的形式拓展到模意义下。
我们令:
$$\omega_{n}^{k}=g^{\frac{k(p-1)}{n}}$$
那么我们就可以简单地将快速傅里叶变换中的单位复根全部替换为原根,这样就实现了快速数论变换。
对于这里的模数\(p\),我们可以选取998244353,它是\(2^{23}*119+1\)。

#include<iostream>
#include<cstdio>
#define Swap(A,B) (A^=B^=A^=B)
const long long P=998244353;
const long long g0=3,gi=332748118;
//*****非常重要:gi是332748118!!!!!!! 
int L=1,R[1<<21|1];
long long invL;
inline int pw(int A,int X){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT=RT*BS%P;
		}
		BS=BS*BS%P;
		X>>=1;
	}
	return RT;
}
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	invL=P-(P-1)/L;
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}
inline void FNTT(int *A,int typ){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[R[i]],A[i]);
		}
	}
	int gn,g,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		gn=pw(~typ?g0:gi,(P-1)/i);
		for(int j=0;j<L;j+=i){
			g=1;
			for(int k=0;k<M;++k,g=1ll*g*gn%P){
				X=A[j+k],Y=1ll*g*A[M+j+k]%P;
				A[j+k]=(X+Y)%P,A[M+j+k]=(X-Y)%P;
			}
		}
	}
}
int n,m,a[1<<21|1],b[1<<21|1];
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%d",a+i);
	}
	for(int i=0;i<=m;++i){
		scanf("%d",b+i);
	}
	prpr(n+m);
	FNTT(a,1);
	FNTT(b,1);
	for(int i=0;i<L;++i){
		a[i]=1ll*a[i]*b[i]%P;
	}
	FNTT(a,-1);
	for(int i=0;i<=n+m;++i){
		printf("%d ",1ll*(a[i]+P)*invL%P);
	}
}

int main(){
	init();
	return 0;
}

lp3803 【模板】多项式乘法(FFT)(非递归版本)

快速傅里叶变换(FFT)的一种较高效的实现

*由于快速傅里叶变换(FFT)常常被用在通信领域,故而FFT的实现效率成为了一个需要被考察的对象。

事实上,观察FFT的实现过程,我们发现它没有非递归不可的理由。

具体来说,我们发现,如果我们想要把FFT改造为递推模式,最大的障碍就是每一次合并操作。
在快速离散傅里叶变换中,每一次合并操作都会导致两个距离为\(\frac{n}{2}\)的点的值合并在一起。因此,如果将原序列传入,那么得到的就是一个乱序的队列。
我们考虑将序列以什么顺序传入。
我们发现,每一次这样的合并操作本质上都是将两个部分「翻转」了。
以0~7为例:
0 1 2 3 4 5 6 7
0 4 2 6 1 5 3 7
将它们写成二进制:
000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111
发现变换之后每个数的二进制位都倒转了。
故而我们可以通过倒转二进制位的方式来预处理原序列,使得FFT之后的序列顺序是正常的。
但是不知道是什么原因,这个代码跑出来的速度竟然也是4500ms…这就很让人挠头了。
于是向PinkRabbit巨神请教了几个经验:
1.不用中间数组,直接在原数组上操作。
2.调整赋值操作,使得修改尽量少。
3.将蝴蝶变换预处理为一个数组。
综合几点,就成功地卡到了2500ms
当然还是有待卡常,毕竟这还是一个比较危险的时间。

#include<iostream>
#include<cstdio>
#include<cmath>
#define PI 3.141592653589793238462

struct Complex{
	double r;
	double i;
	Complex(double ri=0.0,double ii=0.0):r(ri),i(ii){}
	inline Complex operator+(const Complex &B)const{
		return (Complex){r+B.r,i+B.i}; 
	}
	inline Complex operator-(const Complex &B)const{
		return (Complex){r-B.r,i-B.i};
	}
	inline Complex operator*(const Complex &B)const{
		return (Complex){r*B.r-i*B.i,r*B.i+i*B.r};
	}
}a[1<<21|1],b[1<<21|1],s[1<<21|1];
int r[1<<21|1];
inline void FFT(Complex *A,int N,int typ){
	//ButterflyTransform
	for(int i=0;i<N;++i){
		(r[i]>i)?std::swap(A[i],A[r[i]]),0:0;
	}
	Complex bs,nw,x,y;
	for(int i=2,mid;i<=N;i<<=1){
		bs=Complex(cos(2.0*PI/i),typ*sin(2.0*PI/i));
		mid=i>>1;
		for(int j=0;j<N;j+=i){
			nw=Complex(1.0,0.0);
			for(int k=0;k<mid;++k,nw=nw*bs){
				x=A[j+k],y=(A[mid+j+k]*nw);
				A[j+k]=x+y;
				A[mid+j+k]=x-y;
			}
		}
	}
}
int n,m,cnt=1,byt=0;;
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%lf",&a[i].r);
	}
	for(int i=0;i<=m;++i){
		scanf("%lf",&b[i].r);
	}
	while(cnt<=n+m){
		cnt<<=1;
		++byt;
	}
	for(int i=0;i<cnt;++i){
		r[i]=(r[i>>1]>>1)|((i&1)<<(byt-1));
	}
	FFT(a,cnt,1);
	FFT(b,cnt,1);
	for(int i=0;i<cnt;++i){
		a[i]=a[i]*b[i];
	}
	FFT(a,cnt,-1);
	//注意这里typ应该是-1 
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(a[i].r/cnt+0.5));
		//注意输出应该强转格式。 
	}
	puts("");
}

int main(){
	init();
	return 0; 
}

lp3803 【模板】多项式乘法(FFT)

*快速傅里叶变换是一种被广泛运用于各个领域算法,尤其是通信领域。它的主要用途是计算两个函数的卷积。在这里我们讨论的是用于多项式乘法的快速傅里叶变换。

一个\(n\)次多项式函数\(f\)有两种表达方式:点值表达法和系数表达法。
系数表达法:
系数表达法是一个多项式函数最直接的表达法。
对于一个多项式函数,它总是可以表达为这样一个形式:
$$f(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+…a_{n}x^{n}$$
点值表达法:
对于一个多项式函数,我们可以得到它的一个取值集合:
$$\{(x_{0},f(x_{0})),(x_{1},f(x_{1})),…(x_{n-1},f(x_{n-1}))\}$$
这个取值集合本身也代表了这个多项式函数。因为,当我们选取好一个值以后,总是可以通过一系列计算利用这个取值集合得到这个值关于这个多项式函数的映射。

暴力计算两个多项式的\(f,g\)的卷积的复杂度根据其表达法的不同有着不同的复杂度。
暴力计算两个系数表达的多项式的卷积的复杂度是\(O(n^2)\)的,具体来说就是将它们像竖式乘法一样相乘。形如:
$$fg=\sum_{i=0}^{n}\sum_{i=0}^{m}a_{i}b_{j}x^{i+j}$$
暴力计算两个取值相同的点值表达的多项式的卷积的复杂度是\(O(n)\)的,具体来说就是将它们的每个点值表达式的每一位的结果相乘。
$$fg={(x_{0},f(x_{0})g(x_{0})),(x_{1},f(x_{1})g(x_{1})),…(x_{n+m},f(x_{n+m})g(x_{n+m})))}$$

在一些应用情况,我们需要用比\(O(n^2)\)更优的复杂度计算两个多项式函数的系数表达法卷积。
很容易可以想到的是将一个多项式函数的系数表达式转换为点值表达式,计算卷积后再转换回去,这样可能就能够以一个可以接受的复杂度计算两个多项式的系数表达式。
复杂度的瓶颈在于点值表达式和系数表达式的相互转换。
快速傅里叶变换(Fast Fourier Transform)泛指是一种通过选取巧妙的取值,使得点值表达式可以快速地和系数表达式相互转换的算法。其中,将系数表达式转换为点值表达式的转换被称为快速离散傅里叶变换,将点值表达式转换为系数表达式的转换被称为快速逆离散傅里叶变换。
为了完成快速傅里叶变换,我们需要了解我们即将选取的一组取值——单位复根,以及它的性质。

方程\(x^n=1\)的所有解构成了\(n\)个单位复根,记作\(\omega_{n}\),其中的第\(i\)项记作\(\omega_{n}^{i}\)
多项式的单位复根具有一些奇妙的性质:
消去引理:
$$\omega_{dn}^{dk}=e^{\frac{2\pi dki}{dn}}=e^{\frac{2\pi ki}{n}}=\omega_{n}^{k}$$
折半引理:
$$(\omega_{n}^{k+\frac{n}{2}})^2=(e^{\frac{2\pi (k+\frac{n}{2})i}{n}})^2=(e^{\frac{2\pi ki}{n}}*(e^{\pi i}))^2=(-e^{\frac{2\pi ki}{n}})^2=(-\omega_{n}^{k})^2=(\omega_{n}^{k})^2$$
求和引理:
$$\sum_{j=0}^{n-1}(\omega_{n}^{k})^j=\frac{(\omega_{n}^{k})^n-1}{\omega_{n}^{k}-1}=\frac{(e^{\frac{2\pi ki}{n}})^n-1}{\omega_{n}^{k}-1}=\frac{0}{\omega_{n}^{k}-1}=0\ (k>0)$$
折半引理的一个推论:
$$\sum_{k=0}^{n-1}\omega_{n}^{k}=0$$

现在我们考虑如何通过选取单位复根来完成快速地将函数表达在系数表达法和取值表达法之间转换。
我们不由自主地选择分治来解决这个问题。为了实现奇偶分治,我们需要将一个函数的系数表达式扩充到\(2^k-1\)次。
然后我们分别考虑FFT和逆FFT
我们作出了如下定义:
$$f(x)=\sum_{i=0}^{(2^n)-1}a_{i}x^{i}$$
$$f_0(x)=\sum_{i=0}^{2^{(n-1)}}a_{2i}x^{i}$$ $$f_1(x)=\sum_{i=0}^{2^{(n-1)}}a_{2i+1}x^{i}$$
则:
$$f(x)=f_0(x^2)+xf_1(x^2)$$
首先我们考虑如何取值——下面我们证明,当我们已经拥有了一个函数的奇偶分治后的函数的点值表达式,我们可以在线性时间内得到这个函数的点值表达式。
对于某一个取值\(f(\omega_{n}^{k})\ (k\le\frac{n}{2})\),我们有:
$$f(\omega_{n}^{k})=f_0(\omega_{n}^{2k})+\omega_{n}^{k}f_1(\omega_{n}^{2k})$$
$$f(\omega_{n}^{k})=f_0(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}f_1(\omega_{\frac{n}{2}}^{k})$$
$$f(\omega_{n}^{k+\frac{n}{2}})=f_0(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}f_1(\omega_{n}^{2k+n})$$
$$f(\omega_{n}^{k+\frac{n}{2}})=f_0(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}f_1(\omega_{\frac{n}{2}}^{k})$$
通过这样地递归取值,我们可以在\(O(nlog_2n)\)的时间内将一个多项式从它的系数表达法转化为点值表达法。

然后考虑如何快速地将点值表达法转化为系数表达法。
我们将点值表达法获得的取值作为系数,设出了一个函数\(h\),它的每一项的系数都是原函数的点值表达式中的取值:
$$f'(x)=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{ij}x^{j})$$
据超图灵机所说,我们将\(\omega_{n}^{-k}\)代入这个函数,可以得到如下式子:
$$f'(\omega_{n}^{-k})=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{ij}\omega_{n}^{-kj})$$
$$f'(\omega_{n}^{-k})=\sum_{i=0}^{2^n-1}(a_{i}\sum_{j=0}^{2^n-1}\omega_{n}^{j(i-k)})$$
由求和引理可得,当且仅当\(i-k=0\)时,\(\omega_{n}^{j(i-k)}=1\)。其他时候均有\(\omega_{n}^{j(i-k)}=0\)
故有:
$$f'(\omega_{n}^{-k})=a_{k}\sum_{j=0}^{2^n-1}\omega_{n}^{0}=na_{k}$$
故而,我们可以通过将函数\(f’\)通过取每个单位复根的倒数来求得原函数的系数。
而这一过程也可以同样可以通过上述的分治思想来优化。

通过快速离散傅里叶变换将函数从系数表达式转化成点值表达式,然后将点值表达式卷积起来,再将卷积以后的点值表达式通过快速离散傅里叶逆变换转化为系数表达式,这样就在\(O(nlog_2n)\)的时间复杂度内完成了计算多项式函数的系数表达式的卷积。

#include<iostream>
#include<cstdio>
#include<cmath>
#define MID (LEN>>1)
#define PI 3.141592653589793238462

struct Complex{
	double r;
	double i;
	Complex(double rIN=0.0,double iIN=0.0):r(rIN),i(iIN){}
	inline Complex operator+(const Complex &B)const{
		return (Complex){r+B.r,i+B.i};
	}
	inline Complex operator-(const Complex &B)const{
		return (Complex){r-B.r,i-B.i};
	}
	inline Complex operator*(const Complex &B)const{
		return (Complex){r*B.r-i*B.i,r*B.i+i*B.r};
	}
}a[4000005],b[4000005],c[4000005];

inline void FFT(int LEN,Complex *A,Complex *B,int typ){
	if(LEN==1){
		B[0]=A[0];
		return;
	}
//	奇偶项分治 
	for(int i=0;i<MID;++i){
		B[i]=A[i<<1];
		B[i+MID]=A[i<<1|1];
	}
	for(int i=0;i<LEN;++i){
		A[i]=B[i];
	}
	FFT(MID,A,B,typ),FFT(MID,A+MID,B+MID,typ);
	Complex bs(cos((PI*2.0)/LEN),typ*sin((PI*2.0)/LEN)),nw(1.0,0.0);
	for(int i=0;i<MID;++i){
		A[i]=B[i]+nw*B[i+MID];
		A[i+MID]=B[i]-nw*B[i+MID];
		nw=nw*bs;
	}
	for(int i=0;i<LEN;++i){
		B[i]=A[i];
	}
}


int n,m,cnt;
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%lf",&a[i].r);
	}
	for(int i=0;i<=m;++i){
		scanf("%lf",&b[i].r);
	}
	cnt=1;
	while(cnt<=n+m){
		cnt<<=1;
	}
	FFT(cnt,a,c,1);
	FFT(cnt,b,c,1);
	for(int i=0;i<cnt;++i){
		a[i]=a[i]*b[i];
	}
	FFT(cnt,a,c,-1);
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(a[i].r/cnt+0.5));
	}
}

int main(){
	init();
	return 0;
}