多项式相关&&生成函数相关&&一些题目（updating...）

2022-05-17 14:59:15

文章目录

多项式的运算

多项式的加减法，数乘
多项式乘法
多项式求逆
多项式求导
多项式积分
多项式取对
多项式取exp
多项式开方
多项式的除法/取模
分治FFT

生成函数
相关题目

多项式的运算

多项式的加减法，数乘

这个大家应该都会吧

不推了。

直接上代码：

	friend inline poly operator+(const poly&a,const poly&b){
		poly ret(max(a.deg(),b.deg()));
		for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
		for(ri i=0;i<=b.deg();++i)ret[i]=add(ret[i],b[i]);
		return ret;
	}
	friend inline poly operator-(const poly&a,const poly&b){
		poly ret(max(a.deg(),b.deg()));
		for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
		for(ri i=0;i<=b.deg();++i)ret[i]=dec(ret[i],b[i]);
		return ret;
	}
	friend inline poly operator*(const int&a,const poly&b){
		poly ret(b.deg());
		for(ri i=0;i<=b.deg();++i)ret[i]=mul(a,b[i]);
		return ret;
	}

多项式乘法

这个大家应该都会吧

还是推一推吧。

已知的：

A(x)=∑i=0naixiA(x)=\sum_{i=0}^na_ix^iA(x)=∑i=0naixi

B(x)=∑i=0mbixiB(x)=\sum_{i=0}^mb_ix^iB(x)=∑i=0mbixi

要求的：

C(x)=A(x)B(x)=∑i=0n+m(∑j=0min{i,n}aj∗bi−j)xiC(x)=A(x)B(x)=\sum_{i=0}^{n+m}(\sum_{j=0}^{min\{i,n\}}a_j*b_{i-j})x^iC(x)=A(x)B(x)=∑i=0n+m(∑j=0min{i,n}aj∗bi−j)xi

显然直接暴力做是O(n2)O(n^2)O(n2)的，考虑如何优化。

那么我们使用fftfftfft或者nttnttntt来实现点值表示法和系数表示法之间的快速转化。

为了方便起见，我们将A,BA,BA,B的最高次数统一成一个222的幂（对于超过n/mn/mn/m的项的系数看成0即可）

所谓的系数表示法就是我们平常用的那种。

而点值表示法，就是把这个多项式理解成一个函数，用这个函数上的若干个点的坐标来描述这个多项式:

f(x)=(x0,y0),(x1,y1),...,(xn,yn)=p0,p1,p2,...,pnf(x)=(x_0,y_0),(x_1,y_1),...,(x_n,y_n)=p_0,p_1,p_2,...,p_nf(x)=(x0,y0),(x1,y1),...,(xn,yn)=p0,p1,p2,...,pn

假设我们已经将A,BA,BA,B两个函数转化成了点值表示，于是就可以马上求出CCC的点值表示：

A(x)=(xa,0,ya,0),(xa,1,ya,1),...,(xa,n−1,ya,n−1)=pa,0,pa,1,...,pa,n−1A(x)=(x_{a,0},y_{a,0}),(x_{a,1},y_{a,1}),...,(x_{a,n-1},y_{a,n-1})=p_{a,0},p_{a,1},...,p_{a,n-1}A(x)=(xa,0,ya,0),(xa,1,ya,1),...,(xa,n−1,ya,n−1)=pa,0,pa,1,...,pa,n−1

B(x)=(xb,0,yb,0),(xb,1,yb,1),...,(xb,n−1,yb,n−1)=pb,0,pb,1,...,pb,n−1B(x)=(x_{b,0},y_{b,0}),(x_{b,1},y_{b,1}),...,(x_{b,n-1},y_{b,n-1})=p_{b,0},p_{b,1},...,p_{b,n-1}B(x)=(xb,0,yb,0),(xb,1,yb,1),...,(xb,n−1,yb,n−1)=pb,0,pb,1,...,pb,n−1

那么C(x)=pa,0∗pb,0,pa,1∗pb,1,...,pa,n−1∗pb,n−1=pc,0,pc,1,...,pc,n−1C(x)=p_{a,0}*p_{b,0}, p_{a,1}*p_{b,1},...,p_{a,n-1}*p_{b,n-1}=p_{c,0},p_{c,1},...,p_{c,n-1}C(x)=pa,0∗pb,0,pa,1∗pb,1,...,pa,n−1∗pb,n−1=pc,0,pc,1,...,pc,n−1

然后再把C(x)C(x)C(x)还原成系数表达式即可。

注意：我们需要保证xai,xbix_{ai},x_{bi}xai,xbi互不相同）

现在就只用考虑如何实现点值表示和系数表示的互换了。

也就是如何用更少的计算次数来求出nnn个不同的xxx值对应的yyy值。

考虑有两个具有特殊性质的东西：

单位根保证了wn0,wn1,...,wnn−1w_n^0,w_n^1,...,w_n^{n-1}wn0,wn1,...,wnn−1是互不相同的并且有wnij=(wni)jw_n^{ij}=(w_n^i)^jwnij=(wni)j

而原根在模数为质数ppp的时候也有g0,g1,g2,...,gn−1g^0,g^1,g^2,...,g^{n-1}g0,g1,g2,...,gn−1是互不相同的并且gij≡(gi)jmod  pg^{ij}\equiv (g^i)^j \mod pgij≡(gi)jmodp

这满足了我们上面的性质，因此我们考虑将wn0,wn1,...wnn−1w_n^0,w_n^1,...w_n^{n-1}wn0,wn1,...wnn−1作为x0,x1,...xn−1x_0,x_1,...x_{n-1}x0,x1,...xn−1带入求点值。

然后要用到两个引理：

折半引理：wnk∗2=wn2kw_n^{k*2}=w_{\frac n2}^kwnk∗2=w2nk(n为偶数 )
消去引理：wnk=−wnk+n2w_n^{k}=-w_n^{k+\frac n2}wnk=−wnk+2n

这两个引理可以画个单位圆简单证明（建议各位神犇自己简要证明一下）

然后利用按照下标的奇偶性来进行分治处理：

f(x)=∑i=0naixif(x)=\sum_{i=0}^na_ix^if(x)=∑i=0naixi

=>f(wnk)=∑i=0naiwnikf(w_n^k)=\sum_{i=0}^na_i w_n^{ik}f(wnk)=∑i=0naiwnik

=>f(wnk)=∑i=0n2−1a2iwn2ik+wnk∑i=0n2−1a2i+1wn2ikf(w_n^k)=\sum_{i=0}^{\frac n2-1}a_{2i}w_n^{2ik}+w_n^k\sum_{i=0}^{\frac n2-1}a_{2i+1}w_n^{2ik}f(wnk)=∑i=02n−1a2iwn2ik+wnk∑i=02n−1a2i+1wn2ik

同时又有：

f(wnk+n2)=∑i=0n2−1a2iwn2ik−wnk∑i=0n2−1a2i+1wn2ikf(w_n^{k+\frac n2})=\sum_{i=0}^{\frac n2-1}a_{2i}w_n^{2ik}-w_n^k\sum_{i=0}^{\frac n2-1}a_{2i+1}w_n^{2ik}f(wnk+2n)=∑i=02n−1a2iwn2ik−wnk∑i=02n−1a2i+1wn2ik

这一步需要用到引理

所以我们只要算出两个重新分配了系数的多项式的值就可以了。

显然一直分下去只有logloglog层。

于是总时间复杂度O(nlogn)O(nlogn)O(nlogn)

注意到递归的时间复杂度很慢，我们可以预处理最后一层的系数的位置然后用迭代的方式还原回去来优化常数

代码：

vector<int>pos;
int lim,tim;
inline void init(int up){
	lim=1,tim=0;
	while(lim<=up)lim<<=1,++tim;
	pos.resize(lim),pos[0]=0;
	for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}

原理：

我们发现分治完之后相当于将iii和iii对应的二进制数在limlimlim下反转之后对应的新二进制数i′i'i′这两个位置的系数ai,ai′a_i,a_i'ai,ai′交换了位置。

这就成功的实现了系数转点值。

下面来看点值转系数：

注：下面的图片有几点没写清楚：

最后的时候ei,i=0,ei,j∣i̸=je_{i,i}=0,e_{i,j|i\not=j}ei,i=0,ei,j∣i̸=j
那个InI_nIn指的就是对角矩阵。

以上是fftfftfft的证明，nttnttntt同理。

代码：

fft:

//复数定义
struct Complex{
	double x,y;
	friend inline Complex operator+(const Complex&a,const Complex&b){return (Complex){a.x+b.x,a.y+b.y};}
	friend inline Complex operator-(const Complex&a,const Complex&b){return (Complex){a.x-b.x,a.y-b.y};}
	friend inline Complex operator*(const Complex&a,const Complex&b){return (Complex){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
}a[N],b[N];
//fft
inline void fft(Complex a[],int type){
	for(int i=0;i<lim;++i)if(i<pos[i])swap(a[i],a[pos[i]]);
	for(int mid=1;mid<lim;mid<<=1){
		Complex w_n=(Complex){cos(pi/mid),type*sin(pi/mid)};
		for(int j=0,len=mid<<1;j<lim;j+=len){
			Complex w=(Complex){1,0};
			for(int k=0;k<mid;++k,w=w*w_n){
				Complex a0=a[j+k],a1=w*a[j+k+mid];
				a[j+k]=a0+a1,a[j+k+mid]=a0-a1;
			}
		}
	}
}

ntt：

inline void ntt(vector<int>&a,int type){
	for(ri i=0;i<lim;++i)if(i<pos[i])swap(a[i],a[pos[i]]);
	for(ri mid=1,wn,mult=(mod-1)/2,typ=type==1?3:(mod+1)/3;mid<lim;mid<<=1,mult>>=1){
		wn=ksm(typ,mult);
		for(ri j=0,len=mid<<1;j<lim;j+=len)for(ri w=1,a0,a1,k=0;k<mid;++k,w=mul(w,wn)){
			a0=a[j+k],a1=mul(w,a[j+k+mid]);
			a[j+k]=add(a0,a1),a[j+k+mid]=dec(a0,a1);
		}
	}
	if(type==-1)for(ri i=0,inv=ksm(lim,mod-2);i<lim;++i)a[i]=mul(a[i],inv);
}

多项式相乘：

	friend inline poly operator*(const poly&a,const poly&b){
		int n=a.deg(),m=b.deg();
		init(n+m),A.resize(lim),B.resize(lim);
		poly ret(lim-1);
		for(ri i=0;i<=n;++i)A[i]=a[i];
		for(ri i=0;i<=m;++i)B[i]=b[i];
		for(ri i=n+1;i<lim;++i)A[i]=0;
		for(ri i=m+1;i<lim;++i)B[i]=0;
		ntt(A,1),ntt(B,1);
		for(ri i=0;i<lim;++i)A[i]=mul(A[i],B[i]);
		return ntt(A,-1),ret.a=A,ret;
	}

一道板题：uoj#34多项式乘法

 fftfftfft写法

 nttnttntt写法

多项式求逆

定义：

对于一个多项式A(x)A(x)A(x)，如果存在一个多项式B(x)B(x)B(x)，满足B(x)B(x)B(x)的次数小于等于A(x)A(x)A(x)且A(x)B(x)≡1mod  xnA(x)B(x)≡1 \mod x^nA(x)B(x)≡1modxn，那么我们称B(x)为A(x)A(x)A(x)在模xnx^nxn意义下的逆元，简单记作A−1(x)A^{−1}(x)A−1(x)
求法：

n=1?n=1?n=1?那不就是ccc的逆元么。

n>1?n>1?n>1?我们令B(x)=A−1(x)B(x)=A^{-1}(x)B(x)=A−1(x)

那么有A(x)B(x)≡1mod  xnA(x)B(x)\equiv 1 \mod x^nA(x)B(x)≡1modxn

然后可以用类似倍增的方法求。

假设我们已经知道C(x)C(x)C(x)满足A(x)C(x)≡1mod  xn2A(x)C(x)\equiv 1\mod x^{\frac n2}A(x)C(x)≡1modx2n（这里的n2\frac n22n都是向上取整）

显然A(x)B(x)≡1mod  xn2A(x)B(x)\equiv 1\mod x^{\frac n2}A(x)B(x)≡1modx2n是成立的。

我们将两式相减：

A(x)(B(x)−C(x))≡0mod  xn2A(x)(B(x)-C(x))\equiv 0\mod x^{\frac n2}A(x)(B(x)−C(x))≡0modx2n

所以B(x)−C(x)≡0mod  xn2B(x)-C(x)\equiv 0\mod x^{\frac n2}B(x)−C(x)≡0modx2n

然后将两边平方：

B2(x)−2B(x)C(x)+C2(x)≡0mod  xn2B^2(x)-2B(x)C(x)+C^2(x)\equiv 0\mod x^{\frac n2}B2(x)−2B(x)C(x)+C2(x)≡0modx2n

=>B2(x)−2B(x)C(x)+C2(x)≡0mod  xnB^2(x)-2B(x)C(x)+C^2(x)\equiv 0\mod x^nB2(x)−2B(x)C(x)+C2(x)≡0modxn

这一步很关键，请神犇们仔细思考原因

然后两边同时乘上A(x)A(x)A(x)

=>B(x)−2C(x)+A(x)C2(x)≡0mod  xn)B(x)-2C(x)+A(x)C^2(x)\equiv 0\mod x^n)B(x)−2C(x)+A(x)C2(x)≡0modxn)

于是B(x)≡2C(x)−A(x)C2(x)mod  xnB(x)\equiv2C(x)-A(x)C^2(x)\mod x^nB(x)≡2C(x)−A(x)C2(x)modxn

乘法可以用fft/nttfft/nttfft/ntt加速，因为每次递归的时候多项式最高次项都减少一半，所以总复杂度仍然是O(nlogn)O(nlogn)O(nlogn)

代码：

	vector<int>A,B;
	inline poly poly_inv(poly a,const int&k){
		a=a.extend(k);
		if(k==1)return poly(0,ksm(a[0],mod-2));
		poly f0=poly_inv(a,(k+1)>>1);
		return (2*f0-((f0*f0.extend(k))*a).extend(k)).extend(k);
	}

一道板题：洛谷4238

多项式求逆板

多项式求导

默认大家都会函数求导

这个多项式求导属于最简单的那一种

对于一个多项式f(x)=∑i=0naixif(x)=\sum_{i=0}^na_ix^if(x)=∑i=0naixi

它求导的结果f′(x)=∑i=0n−1(i+1)ai+1xif'(x)=\sum_{i=0}^{n-1}(i+1)a_{i+1}x^if′(x)=∑i=0n−1(i+1)ai+1xi

于是直接模拟即可。

代码：

	inline poly poly_direv(poly  a){
		poly ret(a.deg()-1);
		for(ri i=0;i<=ret.deg();++i)ret[i]=mul(a[i+1],i+1);
		return ret;
	}

多项式积分

相当于是多项式求导的逆运算。

对于一个多项式f(x)=∑i=0naixif(x)=\sum_{i=0}^na_ix^if(x)=∑i=0naixi

它求导的结果∫f(x)dx=∑i=1n+1ai−1ixi\int\mathrm f(x)dx=\sum_{i=1}^{n+1}\frac{a_{i-1}}ix^i∫f(x)dx=∑i=1n+1iai−1xi

可以看出来一个函数的导函数积分起来等价于自己。

代码：

	inline poly poly_inter(poly a){
		poly ret(a.deg()+1);
		for(ri i=1;i<=ret.deg();++i)ret[i]=mul(Inv[i],a[i-1]);
		return ret;
	}

多项式取对

我们令g(x)=lnf(x)g(x)=lnf(x)g(x)=lnf(x)

那么根据链式法则求导知：

g′(x)=f′(x)f(x)g'(x)=\frac{f'(x)}{f(x)}g′(x)=f(x)f′(x)

我们已经会多项式求逆和多项式积分，多项式求导了，于是成功解决了多项式取对。

代码：

	inline poly poly_ln(poly a,int len){
		poly ret=a.poly_direv(a);
		return ret=ret*a.poly_inv(a,len),ret.poly_inter(ret);
	}

多项式取exp

这个时候我们要提到一个重要的方法：牛顿迭代法

假设我们要求h(B(x))≡A(x)mod  xnh(B(x))\equiv A(x) \mod x^nh(B(x))≡A(x)modxn中的B(x)B(x)B(x)

现在考虑构造有一个以多项式为变量的函数g(f)=h(f)−Ag(f)=h(f)-Ag(f)=h(f)−A

那么要求的就是g(f)g(f)g(f)模xnx^nxn意义下的零点。

假设已经求出了$g(f)g(f)g(f)模x⌊x2⌋x^{\left\lfloor\frac x 2\right\rfloor}x⌊2x⌋的零点f0f_0f0

那么现在g(fa)=g(f0)+g′(f0)(fa−f0)+g′′(f0)2(fa−f0)2+...g(f_a)=g(f_0)+g'(f_0)(f_a-f_0)+\frac{g''(f_0)}2(f_a-f_0)^2+...g(fa)=g(f0)+g′(f0)(fa−f0)+2g′′(f0)(fa−f0)2+... 泰勒展开

由于f0f_0f0是模x⌊x2⌋x^{\left\lfloor\frac x 2\right\rfloor}x⌊2x⌋的零点，所以有：

g(fa)≡g(f0)+g′(f0)(fa−f0)≡0mod  xng(f_a)\equiv g(f_0)+g'(f_0)(f_a-f_0)\equiv0\mod x^{n}g(fa)≡g(f0)+g′(f0)(fa−f0)≡0modxn

所以移项后发现fa=f0−g(f0)g′(f0)f_a=f_0-\frac{g(f_0)}{g'(f_0)}fa=f0−g′(f0)g(f0)

这个东西有什么用呢？

我们简单举个例子：

比如说多项式求逆，可以构造g(f)=1f−Ag(f)=\frac1f-Ag(f)=f1−A，算出来fa=2f0−Af02f_a=2f_0-Af_0^2fa=2f0−Af02省去了之间的麻烦。

再比如说现在要求的多项式取expexpexp:

构造g(f)=ef−Ag(f)=e^f-Ag(f)=ef−A，算出来fa=f0(1−lnf0+A)f_a=f_0(1-lnf_0+A)fa=f0(1−lnf0+A)

然后就做完了。

代码：

	inline poly poly_exp(poly a,const int&k){
		a=a.extend(k);
		if(k==1)return poly(0,1);
		poly f0=poly_exp(a,(k+1)>>1).extend(k);
		poly delt=a-poly_ln(f0,k);
		++delt[0];
		return (f0*delt).extend(k);
	}

多项式开方

直接使用上面所说的牛顿迭代的结论，令g(x)=f2−Ag(x)=f^2-Ag(x)=f2−A，带入得到:fa=f02+A2f0f_a=\frac{f_0^2+A}{2f_0}fa=2f0f02+A

然后用多项式求逆搞一搞即可。

代码：

	inline poly poly_sqrt(poly a,const int&k){
		a=a.extend(k);
		if(k==1)return poly(0,1);
		poly f0=poly_sqrt(a,(k+1)>>1).extend(k);
		return (((f0*f0).extend(k)+a)*poly_inv((2*f0),k)).extend(k);
	}

多项式的除法/取模

现在有两个多项式:

A(x)=∑i=0naixi,B(x)=∑i=0mbixi,n>mA(x)=\sum_{i=0}^na_ix^i,B(x)=\sum_{i=0}^mb_ix^i,n>mA(x)=∑i=0naixi,B(x)=∑i=0mbixi,n>m

要求出C(x)=∑i=0n−mcixi,D(x)=∑i=0tdixi,d<mC(x)=\sum_{i=0}^{n-m}c_ix^i,D(x)=\sum_{i=0}^td_ix^i,d<mC(x)=∑i=0n−mcixi,D(x)=∑i=0tdixi,d<m，满足A(x)=B(x)C(x)+D(x)A(x)=B(x)C(x)+D(x)A(x)=B(x)C(x)+D(x)，其中C(x)C(x)C(x)类比商,D(x)D(x)D(x)类比余数。

感觉想法比较神奇。

对于一个多项式f(x)f(x)f(x)，我们定义一个fR(x)f^R(x)fR(x)表示将这个多项式的系数翻转之后得到的新多项式，如f(x)=2x3+3x2+x+5f(x)=2x^3+3x^2+x+5f(x)=2x3+3x2+x+5时，fR(x)=5x3+x2+3x+2f^R(x)=5x^3+x^2+3x+2fR(x)=5x3+x2+3x+2

然后可以惊奇的发现：fR(x)=xnf(1x)f^R(x)=x^nf(\frac1x)fR(x)=xnf(x1)

然后就有A(1x)=B(1x)C(1x)+D(1x)A(\frac 1x)=B(\frac 1x)C(\frac 1x)+D(\frac 1x)A(x1)=B(x1)C(x1)+D(x1)

所以xnA(1x)=xnB(1x)C(1x)+xmD(1x)x^nA(\frac 1x)=x^nB(\frac 1x)C(\frac 1x)+x^mD(\frac 1x)xnA(x1)=xnB(x1)C(x1)+xmD(x1)

所以AR(x)=BR(x)CR(x)+DR(x)∗xn−m+t,t∈Z且t>1A^R(x)=B^R(x)C^R(x)+D^R(x)*x^{n-m+t},t\in Z且t>1AR(x)=BR(x)CR(x)+DR(x)∗xn−m+t,t∈Z且t>1

因此AR(x)≡BR(x)CR(x)mod  xn−m+1A^R(x)\equiv B^R(x)C^R(x)\mod x^{n-m+1}AR(x)≡BR(x)CR(x)modxn−m+1

所以CR(x)≡AR(x)(BR(x))−1mod  xn−m+1CR(x)\equiv A^R(x) (B^R(x))^{-1} \mod x^{n-m+1}CR(x)≡AR(x)(BR(x))−1modxn−m+1

D(x)=A(x)−B(x)C(x)D(x)=A(x)-B(x)C(x)D(x)=A(x)−B(x)C(x)

于是我们靠nttnttntt和多项式求逆解决了这一问题。

代码：

	friend inline poly operator/(const poly&a,const poly&b){
		poly ta=a,tb=b;
		int len=1,up=a.deg()-b.deg();
		reverse(ta.a.begin(),ta.a.end()),reverse(tb.a.begin(),tb.a.end());
		ta.extend(up),tb.extend(up);
		while(len<=up)len<<=1;
		tb=tb.poly_inv(tb,len).extend(up);
		return ta=(ta*tb).extend(up),reverse(ta.a.begin(),ta.a.end()),ta;
	}
	friend inline poly operator%(const poly&a,const poly&b){return a-a/b;}

总的代码（缺三角函数以后更新）：

#include<bits/stdc++.h>
#define ri register int
using namespace std;
inline int read(){
    int ans=0;
    char ch=getchar();
    while(!isdigit(ch))ch=getchar();
    while(isdigit(ch))ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar();
    return ans;
}
typedef long long ll;
const int mod=998244353;
int n,lim,tim,m;
vector<int>A,B,pos,Inv;
#define add(a,b) ((a)+(b)>=mod?(a)+(b)-mod:(a)+(b))
#define dec(a,b) ((a)>=(b)?(a)-(b):(a)-(b)+mod)
#define mul(a,b) ((ll)(a)*(b)%mod)
inline int ksm(int a,int p){int ret=1;for(;p;p>>=1,a=mul(a,a))if(p&1)ret=mul(ret,a);return ret;}
inline void ntt(vector<int>&a,const int&type){
    for(ri i=0;i<lim;++i)if(i<pos[i])swap(a[i],a[pos[i]]);
    for(ri mid=1,wn,mult=(mod-1)/2,typ=type==1?3:(mod+1)/3;mid<lim;mid<<=1,mult>>=1){
        wn=ksm(typ,mult);
        for(ri j=0,len=mid<<1;j<lim;j+=len)for(ri w=1,a0,a1,k=0;k<mid;++k,w=mul(w,wn)){
            a0=a[j+k],a1=mul(w,a[j+k+mid]);
            a[j+k]=add(a0,a1),a[j+k+mid]=dec(a0,a1);
        }
    }
    if(type==-1)for(ri i=0,inv=ksm(lim,mod-2);i<lim;++i)a[i]=mul(a[i],inv);
}
inline void init(const int&up){
    lim=1,tim=0;
    while(lim<=up)lim<<=1,++tim;
    pos.resize(lim),pos[0]=0;
    for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}
struct poly{
    vector<int>a;
    inline int deg()const{return a.size()-1;}
    poly(int k,int x=0){a.resize(k+1),a[k]=x;}
    inline int&operator[](const int&k){return a[k];}
    inline const int&operator[](const int&k)const{return a[k];}
    inline poly extend(const int&k){poly ret=*this;return ret.a.resize(k),ret;}
    friend inline poly operator+(const poly&a,const poly&b){
        poly ret(max(a.deg(),b.deg()));
        for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
        for(ri i=0;i<=b.deg();++i)ret[i]=add(ret[i],b[i]);
        return ret;
    }
    friend inline poly operator-(const poly&a,const poly&b){
        poly ret(max(a.deg(),b.deg()));
        for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
        for(ri i=0;i<=b.deg();++i)ret[i]=dec(ret[i],b[i]);
        return ret;
    }
    friend inline poly operator*(const int&a,const poly&b){
        poly ret(b.deg());
        for(ri i=0;i<=b.deg();++i)ret[i]=mul(a,b[i]);
        return ret;
    }
    friend inline poly operator*(const poly&a,const poly&b){
        int n=a.deg(),m=b.deg();
        init(n+m),A.resize(lim),B.resize(lim);
        poly ret(lim-1);
        for(ri i=0;i<=n;++i)A[i]=a[i];
        for(ri i=0;i<=m;++i)B[i]=b[i];
        for(ri i=n+1;i<lim;++i)A[i]=0;
        for(ri i=m+1;i<lim;++i)B[i]=0;
        ntt(A,1),ntt(B,1);
        for(ri i=0;i<lim;++i)A[i]=mul(A[i],B[i]);
        return ntt(A,-1),ret.a=A,ret;
    }
    inline poly poly_inv(poly a,const int&k){
        a=a.extend(k);
        if(k==1)return poly(0,ksm(a[0],mod-2));
        poly f0=poly_inv(a,(k+1)>>1);
        return (2*f0-((f0*f0.extend(k))*a).extend(k)).extend(k);
    }
    inline poly poly_direv(poly a){
        poly ret(a.deg()-1);
        for(ri i=0;i<=ret.deg();++i)ret[i]=mul(a[i+1],i+1);
        return ret;
    }
    inline poly poly_inter(poly a){
        poly ret(a.deg()+1);
        for(ri i=1;i<=ret.deg();++i)ret[i]=mul(Inv[i],a[i-1]);
        return ret;
    }
    inline poly poly_ln(poly a,int len){
        poly ret=a.poly_direv(a);
        return ret=ret*a.poly_inv(a,len),ret.poly_inter(ret);
    }
    inline poly poly_exp(poly a,const int&k){
        a=a.extend(k);
        if(k==1)return poly(0,1);
        poly f0=poly_exp(a,(k+1)>>1).extend(k),delt=a-poly_ln(f0,k);
        ++delt[0];
        return (f0*delt).extend(k);
    }
    inline poly poly_sqrt(poly a,const int&k){
        a=a.extend(k);
        if(k==1)return poly(0,1);
        poly f0=poly_sqrt(a,(k+1)>>1).extend(k);
        return (((f0*f0).extend(k)+a)*poly_inv((2*f0),k)).extend(k);
    }
    friend inline poly operator/(const poly&a,const poly&b){
		poly ta=a,tb=b;
		int len=1,up=a.deg()-b.deg();
		reverse(ta.a.begin(),ta.a.end()),reverse(tb.a.begin(),tb.a.end());
		ta.extend(up),tb.extend(up);
		while(len<=up)len<<=1;
		tb=tb.poly_inv(tb,len).extend(up);
		return ta=(ta*tb).extend(up),reverse(ta.a.begin(),ta.a.end()),ta;
	}
	friend inline poly operator%(const poly&a,const poly&b){return a-a/b;}
};
int main(){
	...
    return 0;
}

分治FFT

一个听起来挺高大上的东西，然而很简单。

前置知识：cdqcdqcdq分治，fftfftfft,nttnttntt

考虑这样一个转移式子f0=0,fi=∑j=1ifi−jgjf_0=0,f_i=\sum_{j=1}^if_{i-j}g_jf0=0,fi=∑j=1ifi−jgj，其中ggg数组已知，让你求fff数组。

容易观察到，这个转移式是一个卷积的形式。

然而并不能直接fftfftfft因为fff的转移跟自身有关，当然你可以使用我们马上下面讲的生成函数秒掉。

生成函数做法：

我们对f,gf,gf,g构造生成函数F(x),G(x)F(x),G(x)F(x),G(x)

那么F(x)−f0=G(x)F(x)F(x)-f_0=G(x)F(x)F(x)−f0=G(x)F(x)

所以F(x)=f0G(x)−1F(x)=\frac{f_0}{G(x)-1}F(x)=G(x)−1f0，直接上多项式求逆即可，时间复杂度为O(nlogn)O(nlog_n)O(nlogn)吊打分治FFT。

新的做法：分治FFT

假设现在要求的fff值的下标为[l,r][l,r][l,r]，我们可以利用类似cdqcdqcdq分治的思想，先递归求出[l,mid][l,mid][l,mid]这一段的fff值，然后考虑[l,mid][l,mid][l,mid]对[mid+1,r][mid+1,r][mid+1,r]的贡献，最后递归[mid+1,r][mid+1,r][mid+1,r]即可。

那么现在要考虑的就只有[l,mid][l,mid][l,mid]对[mid+1,r][mid+1,r][mid+1,r]的贡献啦！

对于fi,i∈[mid+1,r],fj,j∈[l,mid],fi+=fj∗gi−jf_i,i\in[mid+1,r],f_j,j\in[l,mid],f_i+=f_j*g_{i-j}fi,i∈[mid+1,r],fj,j∈[l,mid],fi+=fj∗gi−j，然后把所有的放在一起考虑就成了一个卷积的形式。

于是我们将fl,l+1,...,midf_{l,l+1,...,mid}fl,l+1,...,mid构成的多项式和g0,1,...,r−lg_{0,1,...,r-l}g0,1,...,r−l构成的多项式乘起来更新fmid+1,mid+2,...,rf_{mid+1,mid+2,...,r}fmid+1,mid+2,...,r即可。

然而时间复杂度为O(nlogn2)O(nlog_n^2)O(nlogn2)

代码：

inline void cdqFFT(poly&a,poly&b,int l,int r){
    if(l==r)return;
    int mid=l+r>>1;
    cdqFFT(a,b,l,mid),init(2*(r-l+1));
    for(ri i=0;i<lim;++i)A[i]=B[i]=0;
    for(ri i=l;i<=mid;++i)A[i-l]=a[i];
    for(ri i=0;i<=r-l;++i)B[i]=b[i];
    ntt(A,1),ntt(B,1);
    for(ri i=0;i<lim;++i)A[i]=mul(A[i],B[i]);
    ntt(A,-1);
    for(ri i=mid+1;i<=r;++i)a[i]=add(a[i],A[i-l]);
    cdqFFT(a,b,mid+1,r);
}

生成函数

生成函数简介（百度百科）
前置知识：泰勒展开（fromfromfrom 百度百科）

我对于生成函数的理解：生成函数就相当于对一个集合的表示：

一般生成函数(Ordinary Generating Function) 也就是大家常说的OGFOGFOGF:

F(x)=∑i=0∞aixiF(x)=\sum_{i=0}^{\infty}a_ix^iF(x)=∑i=0∞aixi

指数生成函数(Exponential Generating Function) 也就是大家常说的EGFEGFEGF:

F(x)=∑i=0∞aixii!F(x)=\sum_{i=0}^{\infty}a_i\frac{x^i}{i!}F(x)=∑i=0∞aii!xi

它们可以帮助我们处理一些组合问题。

两个经常用到的公式：

1+x+x2+...=11−x1+x+x^2+...=\frac{1}{1-x}1+x+x2+...=1−x1

1+x+x3+...+xn=1−xn+11−x1+x+x^3+...+x^n=\frac{1-x^{n+1}}{1-x}1+x+x3+...+xn=1−x1−xn+1小学的等比数列求和公式

以及我们的泰勒展开公式。

解决组合问题的时候我们通常将xix^ixi的系数看成值为iii的数被凑出的方案数

事实上，我认为加法原理和乘法原理在生成函数上同样对应了具体的运算。

比如说现在我们要从集合A1,A2,A3,..,AnA_1,A_2,A_3,..,A_nA1,A2,A3,..,An中选一个值为iii的数出来问有多少种选法（加法原理），那么考虑这nnn个生成函数的和CCC，CCC中xix^ixi的系数就是答案。

再比如说我们要从A1,A2,A3,...,AnA_1,A_2,A_3,...,A_nA1,A2,A3,...,An中各选一个数加起来，问加和为iii的方案数（乘法原理），那么考虑这nnn个生成函数的积CCC，CCC中xix^ixi的系数就是答案。

是不是有点感觉了？

我们举一个例子：正整数集N={1,2,3,4,...}N=\{1,2,3,4,...\}N={1,2,3,4,...}，元素的大小定义为它的数值，定义SEQ(A)SEQ(A)SEQ(A) 是由AAA的元素排成的序列组成的集合，一个序列的大小定义为其元素大小总和，现在让我们求一求SEQ(N)SEQ(N)SEQ(N)。

SEQ(N)={正整数有序拆分}={0,1,1+1,2,1+1+1,1+2,2+1,3}SEQ(N) = \{正整数有序拆分\}=\{0,1,1+1,2,1+1+1,1+2, 2+1,3\}SEQ(N)={正整数有序拆分}={0,1,1+1,2,1+1+1,1+2,2+1,3}

考虑构造一个函数N(x)=x+x2+...=x1−xN(x)=x+x^2+...=\frac x{1-x}N(x)=x+x2+...=1−xx

于是SEQ(N)=1+N(x)+N2(x)+...=11−N(x)=1+x1−2x=1+x+2x+4x2+...SEQ(N)=1+N(x)+N^2(x)+...=\frac 1{1-N(x)}=1+\frac x{1-2x}=1+x+2x+4x^2+...SEQ(N)=1+N(x)+N2(x)+...=1−N(x)1=1+1−2xx=1+x+2x+4x2+...

是不是感觉挺好用的

那么来一道板题：bzoj3028

我的做法

再来一道给各位熟悉一下泰勒展开：poj3734

生成函数做法

码农公寓