专栏文章

多项式总结

GGreat_Influence
算法·理论参与者 42已保存评论 50

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
50 条
当前快照
1 份
快照标识符
@mhz5sk8e
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文
此处多项式算法，主要针对基于快速傅里叶变换 (FFT)(FFT) 或者快速数论变换 (NTT)(NTT) 的多项式演算，主要用来解决各种母函数、排列组合、字符串题目。此处仅介绍相关底层算法。文章可能不定期更新。
[TOC]

1.多项式乘法

点我
一切算法的基础。直接做为 O(n2)O(n^2) ，利用 FFTFFT 或者 NTTNTT 来实现可以做到 O(nlog⁡n)O(n\log n)。
配合预处理单位根可以提升精度和速度。

代码:

CPP

inline void FFT(comp*F,int tl)
{
	rep(i,0,Len)if(i<rev[i])swap(F[rev[i]],F[i]);
	for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
		for(register int j=0;j<Len;j+=i)rep(k,0,ii)
		{
			comp x=F[j+k+ii]*g[t][k];
			F[j+k+ii]=F[j+k]-x;
			F[j+k]=F[j+k]+x;
		}
	if(tl==-1)
	{
		reverse(F+1,F+Len);
		rep(i,0,Len)F[i]=F[i]/Len;
	}
}

inline void NTT(int*F,int typ)
{
	Rep(i,1,Len-1)if(i<rev[i])swap(F[i],F[rev[i]]);
	for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
		for(register int j=0;j<Len;j+=i)rep(k,0,ii)
		{
			register int tt=(uint64)F[j+k+ii]*g[t][k]%mod;
			F[j+k+ii]=ad(F[j+k],mod-tt);
			F[j+k]=ad(F[j+k],tt);
		}
	if(typ==-1)
	{
		reverse(F+1,F+Len);
		register uint64 invn=power(Len,mod-2);
		rep(i,0,Len)F[i]=invn*F[i]%mod;
	}
}

int X[MAXN],Y[MAXN];

inline void mul(int *a,int *b,int *c,int lenl,int lenr)
{
	if(lenl<=300/lenr)
	{
		Rep(i,0,lenl+lenr)c[i]=0;
		Rep(i,0,lenl)Rep(j,0,lenr)c[i+j]=(c[i+j]+a[i]*b[j])%mod;
		return;
	}
	for(Len=2;Len<=lenl+lenr;Len<<=1);
	Rep(i,0,lenl)X[i]=a[i];
	Rep(i,0,lenr)Y[i]=b[i];
	rep(i,lenl+1,Len)X[i]=0;
	rep(i,lenr+1,Len)Y[i]=0;
	NTT(X,1),NTT(Y,1);
	rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
	NTT(X,-1);
	Rep(i,0,lenl+lenr)c[i]=X[i];
}
注意 FFTFFT 的复数类最好自己手写，比系统自带的快多了。
以下几个操作大多基于取模来实现，存在一个特殊的模数，满足其减 11 后存在较多 22 的因数。无说明一般为 998244353998244353 。

2.多项式求逆

点我
也算是基础吧。

推导过程:

设要求 FF 关于 xnx^n 的逆，则
设
fi(x)≡F−1(x)(modx2i)f_i(x)\equiv F^{-1}(x)\pmod{x^{2^i}}
假设已知 fif_i ，要求 fi+1f_{i+1} ，那么
F(x)fi(x)≡F(x)fi+1(x)(modx2i)F(x)f_i(x)\equiv F(x)f_{i+1}(x)\pmod{x^{2^i}}
[fi(x)−fi+1(x)]≡0(modx2i)[f_i(x)-f_{i+1}(x)]\equiv 0\pmod{x^{2^i}}
两边平方，则
[fi(x)−fi+1(x)]2≡0(modx2i+1)[f_i(x)-f_{i+1}(x)]^2\equiv 0\pmod{x^{2^{i+1}}}
fi2(x)−2fi(x)fi+1(x)+fi+12(x)≡0(modx2i+1)f_i^2(x)-2f_i(x)f_{i+1}(x)+f_{i+1}^2(x)\equiv 0\pmod{x^{2^{i+1}}}
同乘FF,得
fi+1(x)−2fi(x)+F(x)fi2(x)≡0(modx2i+1)f_{i+1}(x)-2f_i(x)+F(x)f_i^2(x)\equiv 0\pmod{x^{2^{i+1}}}
fi+1(x)≡2fi(x)−F(x)fi2(x)(modx2i+1)f_{i+1}(x)\equiv 2f_i(x)-F(x)f_i^2(x)\pmod{x^{2^{i+1}}}
利用这个公式直接倍增就好了。时间复杂度利用主定理可以证明为 O(nlog⁡n)O(n\log n) 。

代码:

CPP

inline void Inv(int*F,int*G,int ln)
{
	Iv[0]=power(F[0],mod-2);
	rep(i,0,ln)X[i]=Y[i]=0;
	for(register int Ln=2;Ln<=ln;Ln<<=1)
	{
		rep(i,0,Ln)X[i]=F[i];
		rep(i,0,(Ln>>1))Y[i]=Iv[i];
		Len=Ln,calrev();
		NTT(X,1),NTT(Y,1);
		rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
		NTT(X,-1);
		rep(i,0,(Ln>>1))X[i]=0;
		NTT(X,1);
		rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
		NTT(X,-1);
		rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
	}
	rep(i,0,ln)G[i]=Iv[i];
}

3.多项式除法

点我
已知函数FF和GG，求QQ和RR，使得
F=QG+RF=QG+R
保证 Deg(R)<Deg(G)<Deg(F)Deg(R)<Deg(G)<Deg(F) 。

推导过程:

先考虑一种操作:
FR(x)=xDeg(F)F(1x)F_R(x)=x^{Deg(F)}F(\frac{1}{x})
拆开发现该操作的作用就是将FF的系数反位。
然后开始化式子。
为了简便，设Deg(F)=n,Deg(G)=mDeg(F)=n,Deg(G)=m，那么可知Deg(Q)=n−m,Deg(R)=m−1Deg(Q)=n-m,Deg(R)=m-1。
∵F(x)=G(x)Q(x)+R(x)\because F(x)=G(x)Q(x)+R(x)
∴F(1x)=G(1x)Q(1x)+R(1x)\therefore F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})
∴xnF(1x)=xnG(1x)Q(1x)+xnR(1x)\therefore x^nF(\frac{1}{x})=x^nG(\frac{1}{x})Q(\frac{1}{x})+x^nR(\frac{1}{x})
∴xnF(1x)=xmG(1x)xn−mQ(1x)+xn−m+1∗xm−1R(1x)\therefore x^nF(\frac{1}{x})=x^mG(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}*x^{m-1}R(\frac{1}{x})
∴FR(x)=GR(x)QR(x)+xn−m+1∗RR(x)\therefore F_R(x)=G_R(x)Q_R(x)+x^{n-m+1}*R_R(x)
将这个式子放在模n−m+1n-m+1次方下，可以发现RRR_R的所有系数正好全部被模掉，失去影响。
∴FR(x)=GR(x)QR(x)(modxn−m+1)\therefore F_R(x)=G_R(x)Q_R(x)\pmod{x^{n-m+1}}
∴QR(x)=FR(x)GR−1(x)(modxn−m+1)\therefore Q_R(x)=F_R(x)G_R^{-1}(x)\pmod{x^{n-m+1}}
多项式求逆即可。至于RR，直接R=F−QGR=F-QG就可以求出。时间复杂度 O(nlog⁡n)O(n\log n) 。

代码:

CPP
static int X[MAXN],Y[MAXN];

inline void Div(int *a,int n,int *b,int m)
{
	memcpy(X,a,sizeof X);memcpy(Y,b,sizeof Y);
	reverse(b,b+m+1);reverse(X,X+n+1);
	Rep(i,n-m+1,n)X[i]=0;
	Inv(b,n-m+1);
	while(Len<=(n<<2))Len<<=1;
	calrev(Len);
	mul(X,b);reverse(X,X+n-m+1);
	Rep(i,n-m+1,n)X[i]=0;
	mul(Y,X);
	Rep(i,0,m-1)b[i]=ad(a[i],mod-Y[i]);
	memcpy(a,X,sizeof X);
}

4.多项式开方

点我
开始有点麻烦了。依旧采用倍增来求。
已知FF，求GG使得G2≡F(modxn)G^2\equiv F\pmod{x^n}
假设求出 T02≡F(modxt)T_0^2\equiv F\pmod{x^{t}} ，要求 T2≡F(modx2t)T^2\equiv F\pmod{x^{2t}} 。
有
T02−T2≡0(modxt)T_0^2-T^2\equiv 0\pmod{x^t}
(T0+T)(T0−T)≡0(modxt)(T_0+T)(T_0-T)\equiv 0\pmod{x^t}
因为存在平方，所以存在2个解。设 T−T0≡0(modxt)T-T_0\equiv 0\pmod{x^{t}} ，则再平方有
T′2−2TT+T2≡0(modx2t)T'^2-2T_T+T^2\equiv 0\pmod{x^{2t}}
又因为 T2≡F(modx2t)T^2\equiv F\pmod{x^{2t}} ，则
T′2−2TT0+F≡0(modx2t)T'^2-2TT_0+F\equiv 0\pmod{x^{2t}}
T≡T02+F2T0(modx2t)T\equiv \frac{T_0^2+F}{2T_0}\pmod{x^{2t}}
直接倍增计算即可。时间复杂度依旧 O(nlog⁡n)O(n\log n) 。

代码:

CPP
inline void Sqrt(int*F,int*G,int ln)
{
	G[0]=1;
	for(register int Ln=2;Ln<=ln;Ln<<=1)
	{
		mul(G,G,ExX,Ln>>1,Ln>>1);
		rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
		rep(i,Ln>>1,Ln)ExX[i]=0;
		Inv(G,ExY,Ln>>1);
		Len=Ln,calrev();
		NTT(ExX,1),NTT(ExY,1);
		rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
		NTT(ExX,-1);
		rep(i,0,Ln>>1)G[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
	}
}

5.多项式积分和求导

这就不要例题了吧...
设F(x)=∑i=0naixiF(x)=\sum\limits_{i=0}^n a_ix^i，则
定义:
∫F(x)=∑i=0naii+1xi+1\int F(x)=\sum\limits_{i=0}^n \frac{a_i}{i+1}x^{i+1}
dFdx=∑i=1niaixi−1\frac{dF}{dx}=\sum\limits_{i=1}^n ia_ix^{i-1}
直接O(n)O(n)计算就可以了。

代码:

CPP
inline void Direv(int *a,int n)
{Rep(i,1,n)a[i-1]=(ll)a[i]*i%mod;a[n]=0;}
	
inline void Inter(int *a,int n)
{Repe(i,n,1)a[i]=(ll)a[i-1]*power(i,mod-2)%mod;a[0]=0;}

6.多项式取 ln⁡\ln

模板题点我
做完模板可以点我
ln⁡\ln 的定义为 exp⁡\exp 的逆变换。具体定义请参考下面的 exp⁡\exp 。
定义式为:
ln⁡(1−x)=−∑i=1∞xii\ln(1-x)=-\sum_{i=1}^\infty \frac{x^i}{i}
存在推导式:
ln⁡F=∫F′F≡∫F′F−1\ln F=\int \frac{F'}{F}\equiv \int F'F^{-1}
直接计算即可。

代码:

CPP
inline void Ln(int *a,int n)
{
	memcpy(X,a,sizeof X);
	Inv(X,n);Direv(a,n);
	while(Len<=(n<<2))Len<<=1;
	calrev(Len);
	mul(a,X);
	Inter(a,n);
	Rep(i,n+1,Len)a[i]=0;
}

7.多项式牛顿迭代

题目请见 WC2019WC2019 汪乐平的课件。
采用倍增实现。
已知GG，求多项式FF使得G(F(x))=0G(F(x))= 0。
设已知F0F_0满足G(F0(x))≡0(modxt)G(F_0(x))\equiv 0\pmod{x^t}，要求G(F(x))≡0(modx2t)G(F(x))\equiv 0\pmod{x^{2t}}。
存在递推式
F(x)≡F0(x)−G(F0(x))G′(F0(x))(modx2t)F(x)≡F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^{2t}}
其中，G′(F(x))=dGdFG'(F(x))=\frac{\mathbb{d}G}{\mathbb{d}F}，即 GG 对 FF 的导数。
直接计算即可。

7.1推导过程:

函数 GG 对 FF 进行泰勒展开，得:
G(F)=G(F0)+G′(F0)(F−F0)+G′′(F0)(F−F0)2+⋯G(F)=G(F_0)+G'(F_0)(F-F_0)+G''(F_0)(F-F_0)^2+\cdots
因为F≡F0(modxt)F\equiv F_0\pmod{x^t}，因此F−F0≡0(modxt)F-F_0\equiv 0\pmod{x^t}，因此(F−F0)k≡0(modx2t),k≥2(F-F_0)^k\equiv 0\pmod{x^{2t}},k\ge 2。
因此后面几项在 mod x2t\bmod x^{2t}意义下全部为00，直接忽略。
因此得到关系式
G(F)=G(F0)+G′(F0)(F−F0)G(F)=G(F_0)+G'(F_0)(F-F_0)
整理即得上式。
(代码依 GG 不同而不同)

7.2.其他式子利用牛顿迭代推导的方式

7.2.1 求逆

考虑对方程F(x)G(x)−1=0F(x)G(x)-1=0 作牛顿迭代，得:
δ(F)=GF−1\delta(F)=GF-1
δ′(F)=G\delta'(F)=G
可以得到递推式: F=F0−δ(F0)δ′(F0)F=F_0-\frac{\delta(F_0)}{\delta'(F_0)}
=F0−G0F0−1G0=F_0-\frac{G_0F_0-1}{G_0}
=F0−F0(G0F0−1)=F_0-F_0(G_0F_0-1)
=2F0−G0F02=2F_0-G_0F_0^2

7.2.2 开方

考虑对方程F2−G=0F^2-G=0作牛顿迭代，得:
δ(F)=F2−G\delta(F)=F^2-G
δ′(F)=2F\delta'(F)=2F
可以得到递推式:
F=F0−δ(F0)δ′(F0)F=F_0-\frac{\delta(F_0)}{\delta'(F_0)} =F0−F02−G02F0=F_0-\frac{F_0^2-G_0}{2F_0} =F02+G02F0=\frac{F_0^2+G_0}{2F_0}

7.2.3 exp⁡\exp

见下文。

7.3 牛顿法的通用优化技巧

使用牛顿法，我们通常需要多次求出F0−G(F0)G′(F0)F_0-\frac{G(F_0)}{G'(F_0)}。有时为了方便，我们会将前面的F0F_0和后面部分共同计算。
可以发现，因为F≡F0(modxn2)F\equiv F_0\pmod{x^{\frac{n}{2}}}，所以后面部分一定有长度为n2\frac{n}{2}的部分全部都是00。此时，我们可以将后半部分位移至前面，然后再进行操作。通常可以将多项式乘法的长度减少一半，时间复杂度会大大优化。
例如，多项式开方的公式为F=F0+12F0−1(G−F02)F=F_0+\frac{1}{2}F_0^{-1}(G-F_0^2)。应用以上方法可以知道，(G−F02)(G-F_0^2)有一半都是00。那么，如果当前层长度为nn，从n2\frac{n}{2}转移过来，此时我们做多项式乘法只需要做长度为nn的而不是2n2n的。同理，求逆也只需要处理到nn。这样可以省下近一半常数。

8.多项式 exp⁡\exp

点我
由泰勒展开可得 ex=∑i=0∞xii!e^x=\sum_{i=0}^\infty \frac{x^i}{i!}。因此，我们定义 ∑i=0∞F(x)ii!=exp⁡F\sum_{i=0}^\infty \frac{F(x)^i}{i!}=\exp{F}。介于 ee 在膜域下无定义，一般默认F(x)F(x) 常数项为 00 。
已知函数 FF ，求 G≡exp⁡FG\equiv \exp{F} 。保证 FF 常数项为 00 。
可以知道 exp⁡(ln⁡P)=P\exp(\ln P)=P 。
所以 ln⁡G−F=0\ln G-F=0 。
直接套用上方多项式牛顿迭代的公式，可得
G=G0(1−ln⁡G0+F)G=G_0(1-\ln G_0+F)
倍增即可。时间复杂度 O(nlog⁡n)O(n\log n) 。
该算法常数巨大，经验算大约在 70∼8070\sim 80 左右。因此可以粗略认为是 O(nlog⁡2n)O(n\log^2 n) 的。
最近发现常数过大还没有分治 FFTFFT 跑得快。

代码:

CPP
static int ExX[MAXN],ExY[MAXN];

inline void Exp(int *F,int *G,int ln)
{
	G[0]=1;
	for(register int Lx=2;Lx<=ln;Lx<<=1)
	{
		rep(i,Lx>>1,Lx)G[i]=0;
		Ln(G,ExX,Lx);
		rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
		rep(i,0,Lx>>1)ExY[i]=G[i];
		rep(i,Lx>>1,Lx)ExX[i]=0;
		Len=Lx;
		calrev();
		NTT(ExY,1),NTT(ExX,1);
		rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
		NTT(ExX,-1);
		rep(i,0,Len>>1)G[i+(Len>>1)]=ExX[i];
	}
}

9.多项式幂次

例题见10.拉格朗日反演。
直接上式子吧...
当 F[0]=1F[0]=1 时:
Fk=exp⁡(kln⁡F)F^k=\exp(k\ln F)
否则，设最低项为 atxta_tx^t ,则
Fk=atkxtkexp⁡(kln⁡Fatxt)F^k=a_t^kx^{tk}\exp(k\ln \frac{F}{a_tx^t})
直接按照式子计算即可。时间复杂度O(nlog⁡n)O(n\log n) 。

注意:

直接求 ln⁡\ln 再 exp⁡\exp 的方法仅适用于常数项为1的情况(因为多项式 exp⁡\exp 计算时默认常数项为 11 )。因此当常数项不为 11 时需要将多项式的常数项强制变成 11 。多判一下即可(不过大部分时间都不用管)。

代码:

CPP
inline void Pow(int*F,int*G,int ln,ll k)
{
	int lst=ln+1;
	Rep(i,0,ln)if(F[i]){lst=i;break;}
	if(ln&&(__int128)lst*k>ln)
	{memset(G,0,sizeof(int)*(ln+1));return;}
	int md=ln-k*lst,bs=F[lst],iv=power(bs,mod-2);
	Rep(i,lst,ln)G[i]=(ll)G[i]*iv%mod;
	Ln(G+lst,G,md);
	k%=mod;
	Rep(i,0,md)G[i]=(ll)G[i]*k%mod;
	Exp(G,G,md);
	bs=power(bs,k);
	Repe(i,md,0)G[i+lst*k]=(ll)G[i]*bs%mod;
	memset(G,0,sizeof(int)*(lst*k));
}

9.5.多项式 kk 次方根

在膜意义下，存在一下转化:
Fk=F1k=Fk−1\sqrt[k]F=F^{\frac{1}{k}}=F^{k^{-1}}
直接套用上面的幂次即可。复杂度为 O(nlog⁡n)O(n\log n) 。
代码就咕了。

10.拉格朗日反演

点我
已知 GG ，求函数 FF 使得 F[G(x)]=xF[G(x)]=x 。
经过复杂的推理，存在关系式:
[xn]F(x)=1n[xn−1](xG(x))n[x^n] F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n
其中[xn][x^n]表示多项式xnx^n项系数。
直接按照式子模拟即可。时间复杂度 O(nlog⁡n)O(n\log n) 。

代码:

CPP
inline int Lagrange(int *a,int n)
{
	Rep(i,0,n-1)a[i]=a[i+1];
	Inv(a,n);
	pow(a,n,n);
	return (ll)a[n-1]*power(n,mod-2)%mod;
}

10.5 扩展拉格朗日反演

其实就是另一个式子。
[xn]H[F(x)]=1n[xn−1]H′(x)(xG(x))n[x^n]H[F(x)]=\frac{1}{n}[x^{n-1}]H'(x)(\frac{x}{G(x)})^n

11.分治FFT

点我
已知 F[0]=1F[0]=1 ,求 FF 满足 F[i]=∑j=1iF[i−j]∗G[j]F[i]=\sum_{j=1}^i F[i-j]*G[j]
其实哪怕是F[i]=aF[i−b]+c+∑j=1idF[i−j]∗G[j]F[i]=aF[i-b]+c+\sum_{j=1}^i dF[i-j]*G[j]也是可以解决的。
考虑利用分治。
假设我们求出了l→midl\to mid的答案，要求这些点对mid+1→rmid+1\to r的影响，那么对右半边点xx的贡献为:
wx=∑i=lmidf[i]g[x−i]w_x=\displaystyle\sum_{i=l}^{mid} f[i]g[x-i]
这部分可以利用卷积来快速计算。计算完以后，答案直接加到答案数组就可以了。
需要注意的是，如果要求左边点对右边点的影响，首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边，在求出中间，然后在递归右边。
时间复杂度 O(nlog⁡2n)O(n\log^2n) 。
有一个简单卡常方法。求一层的时候，我们会做长度为 md−l+1+r−l+1≈1.5(r−l+1)md-l+1+r-l+1\approx 1.5(r-l+1) 的卷积，但是这没有必要。一个显然的观察是我们只会使用在 md−l+1∼r−l+1md-l+1\sim r-l+1 内的值，而前半截内容显然是没有用的。因此我们可以使用长度为 r−l+1r-l+1 的卷积。一般情况下可以快 50%50\% 。
代码(求 exp⁡\exp):
CPP
void cdq_FFT(int l,int r)
{
	if(l==r)
	{
		if(!l)F[l]=1;
		else F[l]=(ll)F[l]*iv[l]%mod;
		return;
	}
	int md=(l+r)>>1;cdq_FFT(l,md);
	for(Len=2;Len<=r-l;Len<<=1);
	calrev();
	memset(X,0,sizeof(int)*Len);
	memset(Y,0,sizeof(int)*Len);
	memcpy(X,F+l,sizeof(int)*(md-l+1));
	memcpy(Y,a,sizeof(int)*(r-l));
	NTT(X,1),NTT(Y,1);
	rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
	NTT(X,-1);
	Rep(i,md+1,r)F[i]=ad(F[i],X[i-l-1]);
	cdq_FFT(md+1,r);
}

11.5.特殊的分治FFT

点我
先假设你会前面的推式子部分。
你现在得到了一个长得很丑的式子，大概长这样: fn=(n−1)fn−1+∑i=1n−2(i−1)fifn−if_n=(n-1)f_{n-1}+\sum_{i=1}^{n-2}(i-1)f_if_{n-i} 换句话说，你得到了一个类似于F=F+F∗FF=F+F*F 的式子。这个式子利用上面的分治 FFTFFT 方法当然是求不出的了。
这时候，你就需要对你的分治 FFTFFT 进行魔改了。魔改结果大概长这样:
CPP
void cdqFFT(int *f,int l,int r)
{
    if(l==r)
    {
        f[l]=ad(f[l],(uint64)(l-1)*f[l-1]%mod);
        return;
    }
    int mid=(l+r)>>1;
    cdqFFT(f,l,mid);
    for(Len=r-l+1;Len>(Len&-Len);Len+=(Len&-Len));
    Rep(i,0,Len)A[i]=B[i]=0;
    Rep(i,l,mid)A[i-l]=(uint64)f[i]*(i-1)%mod,B[i-l]=f[i];
    calrev(Len);mul(A,B);
    Rep(i,mid+1,r)if(i>=l*2)f[i]=ad(f[i],A[i-2*l]);
    if(l^2)
    {
        Rep(i,0,Len)A[i]=B[i]=0;
        Rep(i,2,min(l-1,r-l))A[i-2]=f[i];
        Rep(i,l,mid)B[i-l]=f[i];
        mul(A,B);
        Rep(i,mid+1,r)if(i>=l+2)
            f[i]=ad(f[i],(uint64)(i-2)*A[i-l-2]%mod);
    }
    cdqFFT(f,mid+1,r);
}
稍微解释一下。
首先，前面的 (n−1)fn−1(n-1)f_{n-1} 明显是可以在递归最底层在加上的。这个很好处理。
问题就在于后面的 ∑i=1n−2(i−1)fifn−i\sum_{i=1}^{n-2}(i-1)f_if_{n-i}。这部分很难求。但是你可以分两种情况来进行讨论:
  1. 所计算区间左端点是全局左端点。原来怎么做就怎么做。但是注意不要在卷积中用到还没有计算完毕的项。
  2. 所计算区间不是全局左端点。可以知道会存在一些项之前并没有计算过。我们将这些项合并一起计算。具体来说就是把 (i−1)fifn−i(i-1)f_if_{n-i} 和 (n−i−1)fn−ifi(n-i-1)f_{n-i}f_i 一起计算。那么这两项加在一起就变成了(n−2)fifn−i(n-2)f_if_{n-i} 。直接加上即可。
其余部分和一般分治FFTFFT相同。注意处理几个边界情况。
注意，我程序内部有一个 if(l⊗2)if(l\otimes 2) ，这是判断左端点是否为全局左端点的。因为这道题的全局左端点是 22 。

12.常系数齐次线性递推

点我
主要就是推式子。
推导过程和代码见线性递推。

13.多项式求三角函数

经过审稿管理员提示，我发现原来这东西还是有用的。。。
已知F(x)F(x)，求p[F(x)]p[F(x)]。其中， p(x)p(x) 为某一个三角函数 (sin⁡,cos⁡,tan⁡,sec⁡,csc⁡,cot⁡)(\sin,\cos,\tan,\sec,\csc,\cot) 。
首先，可以知道，如果求出了sin⁡\sin 和cos⁡\cos ，那么其它的三角函数也可以轻易推出。那么问题就是如何求出这两个三角函数。
根据欧拉定理可知:
ew4x=cos⁡(x)+w4sin⁡(x)e^{w_4x}=\cos(x)+w_4\sin(x)
其中， w4w_4 为四次单位根，也就是常用的 ii 。
那么同理可得
ew4F(x)=cos⁡[F(x)]+w4sin⁡[F(x)]e^{w_4F(x)}=\cos[F(x)]+w_4\sin[F(x)]
将F(x)F(x) 的系数都乘以 w4w_4 然后 exp⁡\exp 就可以了。时间复杂度 O(nlog⁡n)O(n\log n) 。
注意以上的内容是 FFTFFT 的求法。如果要用 NTTNTT 的话也不是不行。
由
cos⁡(x)=∑i=0∞[i mod 2=0](−1)i2xii!\cos(x)=\sum_{i=0}^\infty \frac{[i\bmod 2=0](-1)^\frac{i}{2}x^i}{i!}
cosh⁡(x)=ex+e−x2=∑i=0∞[i mod 2=0]xii!\cosh(x)=\frac{e^x+e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod 2=0]x^i}{i!}
可得
cos⁡(x)=ew4x+e−w4x2\cos(x)=\frac{e^{w_4x}+e^{-w_4x}}{2}
cos⁡(F)=12(ew4F−e−w4F)\cos(F)=\frac{1}{2}(e^{w_4F}-e^{-w_4F})
因为 w4w_4 在常用 NTTNTT 模数下一般都有定义，且实系数函数的 cos⁡\cos 肯定也是实系数函数，因此这个公式计算出来肯定不含复系数。
至于 sin⁡\sin ，可以利用类似的方法来计算。即:
sin⁡(x)=∑i=0∞[i mod 2=1](−1)i−12xii!\sin(x)=\sum_{i=0}^\infty \frac{[i\bmod2=1](-1)^{\frac{i-1}{2}}x^i}{i!}
sinh⁡(x)=ex−e−x2=∑i=0∞[i mod 2=1]xii!\sinh(x)=\frac{e^x-e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod2=1]x^i}{i!}
推得:
sin⁡(x)=w432(ew4x−e−w4x)\sin(x)=\frac{w_4^3}{2}(e^{w_4x}-e^{-w_4x})
同理得:
sin⁡F=w432(ew4F−e−w4F)\sin F=\frac{w_4^3}{2}(e^{w_4F}-e^{-w_4F})

代码:

CPP
static int Sn[MAXN],Cs[MAXN];

inline void Cos(int*F,int*G,int ln)
{
	Rep(i,0,ln)Cs[i]=(ll)F[i]*w4%mod;
	Exp(Cs,Cs,ln),Inv(Cs,G,ln);
	Rep(i,0,ln)G[i]=((ll)Cs[i]+G[i])*inv2%mod;
}

const int w43div2=(ll)w4*w4%mod*w4%mod*inv2%mod;

inline void Sin(int*F,int*G,int ln)
{
	Rep(i,0,ln)Sn[i]=(ll)F[i]*w4%mod;
    Exp(Sn,Sn,ln),Inv(Sn,G,ln);
	Rep(i,0,ln)G[i]=((ll)Sn[i]+mod-G[i])*w43div2%mod;
}

14.多项式反三角函数

这个东西和上面的没有半毛钱关系。。。
根据积分表可以知道:
arcsin⁡x=∫11−x2\arcsin x = \int \frac{1}{\sqrt{1-x^2}}
arctan⁡x=∫11+x2\arctan x = \int \frac{1}{1+x^2}
那么同理可得:
arcsin⁡F=∫F′1−F2\arcsin F = \int \frac{F'}{\sqrt{1-F^2}}
arctan⁡F=∫F′1+F2\arctan F = \int \frac{F'}{1+F^2}
直接计算即可。
代码:
CPP
inline void Arcsin(int*F,int*G,int ln)
{
	memcpy(Sn,F,sizeof(int)*(ln+1));
	for(Len=2;Len<=ln*2;Len<<=1);calrev();
	NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
	NTT(Sn,-1);
	Sn[0]=1;rep(i,ln+1,Len)Sn[i]=0;
	Rep(i,1,ln)Sn[i]=mod-Sn[i];
	Sqrt(Sn,Sn,ln);
	Inv(Sn,Sn,ln);
	Deriv(F,G,ln);
	mul(G,Sn,G,ln,ln);
	Rep(i,ln+1,ln<<1)G[i]=0;
	Inter(G,G,ln);
}

inline void Arctan(int*F,int*G,int ln)
{
	memcpy(Sn,F,sizeof(int)*(ln+1));
	for(Len=2;Len<=ln*2;Len<<=1);calrev();
	NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
	NTT(Sn,-1),Sn[0]=1;
	rep(i,ln+1,Len)Sn[i]=0;Inv(Sn,Sn,ln);
	Deriv(F,G,ln),mul(G,Sn,G,ln,ln);
	Rep(i,ln+1,ln<<1)G[i]=0;
	Inter(G,G,ln);
}

15.多项式多点求值

点我
已知函数 FF 和数组 {ai}\{a_i\} ，求 {F(ai)}\{F(a_i)\} 。
直接计算是 O(nm)O(nm) 的。此时，我们可以考虑降幂。
构造一个函数 Dl,r=∏i=lr(x−ai)D_{l,r}=\prod_{i=l}^r (x-a_i) 。
可以发现，∀k∈[l,r],Dl,r(ak)=0\forall k\in [l,r],D_{l,r}(a_k)=0 。
因此，∀k∈[l,r],F(ak)=F(ak) mod Dl,r(ak)\forall k\in [l,r],F(a_k)=F(a_k)\bmod D_{l,r}(a_k) 。
于是，设 Fl,r(x)F_{l,r}(x) 表示 ala_l 到 ara_r 在原函数上的点构成的多项式，则:
Fl,mid(x)=Fl,r(x) mod Dl,mid(x)F_{l,mid}(x)=F_{l,r}(x)\bmod D_{l,mid}(x)
直接分治即可。时间复杂度为 O(nlog⁡2n)O(n\log^2 n) 。注意 DD 的计算可以直接分治乘法每次乘起来，不影响复杂度。

代码:

CPP
static int solv[18][2][MAXN];

inline void calc(int*a,int l,int r,int lev,int dir)
{
	if(l==r)
	{
		solv[lev][dir][0]=mod-a[l];
		solv[lev][dir][1]=1;return;
	}
	int md=(l+r)>>1;
	calc(a,l,md,lev+1,0),calc(a,md+1,r,lev+1,1);
	mul(solv[lev+1][0],solv[lev+1][1],solv[lev][dir],md-l+1,r-md);
}

static int pol[18][MAXN],Q[MAXN];

void getnum(int*a,int*ans,int l,int r,int lev,int n)
{
	if(n>r-l)
	{
		calc(a,l,r,0,0);
		Div(pol[lev-(lev>0)],solv[0][0],Q,pol[lev],n,r-l+1);
		n=r-l;
	}
	else if(lev){Rep(i,0,n)pol[lev][i]=pol[lev-1][i];}
	if(l==r)
	{
		ans[l]=pol[lev][0];
		return;
	}
	int md=(l+r)>>1;
	getnum(a,ans,l,md,lev+1,n);
	getnum(a,ans,md+1,r,lev+1,n);
}

16.多项式快速插值

点我
给定二维平面上点集{(xi,yi)}\{(x_i,y_i)\}(保证xix_i互不相同)，求函数F(x)F(x)，满足所有点均在该函数上。
核心内容是优化拉格朗日插值。
具体过程见多项式快速插值。

Ex1.MTT

点我还有我
上面的所有函数基于模数为一个 NTTNTT 模数(即 a∗2k+1a*2^k+1 ，其中kk比较大)。如果不是的话，可以将多项式乘法换成以下的函数:
CPP
inline void mul(int*F,int*G)
{
	rep(i,0,Len)A[i]=comp(F[i]>>15,F[i]&32767),
		C[i]=comp(G[i]>>15,G[i]&32767);
	DBLFFT(A,B,1),DBLFFT(C,D,1);
	rep(i,0,Len)
	{
		register comp t=A[i];
		A[i]=B[i]*C[i]+A[i]*D[i];
		B[i]=B[i]*D[i];
		C[i]=C[i]*t;
	}
	DBLFFT(A,B,-1),FFT(C,-1);
	rep(i,0,Len)
	{
		register int lz=(((ll)floor(C[i].re+0.5))%mod*73741817
			 +((ll)floor(A[i].re+0.5))%mod*32768+(ll)floor(B[i].re+0.5))%mod;
		F[i]=i<=k?lz:0;
	}
}
其中 DBLFFTDBLFFT 是共轭优化过的 FFTFFT 。下面一章有代码。
简单解释其实就是将系数拆成 215a+b2^{15}a+b ，然后利用等式 (215A+B)(215C+D)=230AC+215(BC+AD)+BD(2^{15}A+B)(2^{15}C+D)=2^{30}AC+2^{15}(BC+AD)+BD 直接计算。原理为 doubledouble 的精度可以承受住 101410^{14} ，因此能够支撑 215∗215∗105=1073741824000002^{15}*2^{15}*10^5=107374182400000 以内的数字运算。
具体还是看 20162016 毛爷爷写的国家队论文吧。
顺便一提，配合共轭优化后可以做到 44 次 DFTDFT ，基本和一般的乘法没什么区别了。

Ex2.BlueStein算法

这个算例题吧
已知多项式 F(x)F(x) ，求 ∀i∈[0,k),Ai=F(wki)\forall i\in[0,k),A_i=F(w_k^i) 。其中 wkw_k 为 kk 次单位根。
被这个东西送到了退役边缘。。。只怪我自己没有认真看论文。。。
因为 wkt=wkt mod kw_k^t=w_k^{t\bmod k} ，因此可以认为 FF 的最高次幂为 k−1k-1 。
设 FF 的 xtx^t 次方系数为 ftf_t 。
那么，将式子写出，得:
F(wkt)=∑i=0k−1fiwktiF(w_k^t)=\sum_{i=0}^{k-1}f_iw_k^{ti}
=∑i=0k−1fiwk−(k−t)i=\sum_{i=0}^{k-1}f_iw_k^{-(k-t)i}
设这个东西为 ansk−tans_{k-t} ，则 F(wkt)=ansk−tF(w_k^t)=ans_{k-t} 。
运用套路，得:
ti=(t+i2)−(t2)−(i2)ti=\binom{t+i}{2}-\binom{t}{2}-\binom{i}{2}
因此有:
anst=∑i=0k−1fiwk(i−t2)−(−t2)−(i2)ans_t=\sum_{i=0}^{k-1} f_iw_k^{\binom{i-t}{2}-\binom{-t}{2}-\binom{i}{2}}
=wk−(−t2)∑i=0k−1fiwk−(i2)∗wk(−(t−i)2)=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{-(t-i)}{2}}
因为我们每一项都要卷到 k−1k-1 ，因此考虑将系数直接后移 k−1k-1 位，可以得到新的式子:
anst=wk−(−t2)∑i=0k−1+tfiwk−(i2)∗wk(k−1−(k−1+t−i)2)ans_t=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1+t} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{k-1-(k-1+t-i)}{2}}
我们设 Ai=fiwk−(i2)A_i=f_iw_k^{-\binom{i}{2}} ， Bi=wk(k−1−(k−1+t−i)2)B_i=w_k^{\binom{k-1-(k-1+t-i)}{2}} ，则
anst=wk−(−t2)[k−1+t](A∗B)ans_t=w_k^{-\binom{-t}{2}}[k-1+t](A*B)
时间复杂度 O(nlog⁡n)O(n\log n) 。
这个过程被称为 CZTCZT 。 ICZTICZT 只需要将 F1F_1 到 Fk−1F_{k-1} 的部分翻转，然后整体除掉 kk 即可。
算出这个东西后，我们就可以将他当做以 wkw_k 为公比的插值。因此如果将多个式子相卷，得到的结果将会是这些式子在膜 xkx^k 意义下的循环卷积。有的时候会有用。
代码:
CPP
inline void Blue_Stein(int*F,int*G,int ln)
{
	rep(i,0,ln)A[i]=(ll)F[i]*power(ome,ln-(ll)i*(i-1)/2%ln,mod)%mod;
	Rep(i,0,2*ln-2)B[i]=power(ome,ln+(ll)(ln-1-i)*(ln-2-i)/2%ln,mod);
	mul(A,B,A,ln-1,2*ln-2);
	reverse(A+1,A+k);
	rep(i,0,ln)G[i]=(ll)A[i+k-1]*power(ome,ln-(ll)(-i)*(-i-1)/2%ln,mod)%mod;
}
注意使用的时候大部分情况下模数都不会很友好，因此结合上面的 MTTMTT 或者直接上 FFTFFT 是一般情况。

Ex3.常数优化

根据经验可以发现，多项式方面的常数差异非常大，有的时候别人跑过了但是你没跑过就凉了。因此，常数优化是非常重要的。限于篇幅就只推荐这篇和这篇了。后面那篇的详细内容和不详细证明可以参考毛爷爷的国家队论文。
此处放出共轭优化的代码。该代码可以同时对 FF 和 GG 进行 DFTDFT 或者 IDFTIDFT 操作。
CPP

static comp Z[MAXN];

inline void DBLFFT(comp*F,comp*G,int tl)
{
	if(tl==1)
	{
		memcpy(Z,F,sizeof(comp)*Len);
		FFT(Z,1);
		memcpy(F,Z,sizeof(comp)*Len);
		reverse(F+1,F+Len);
		rep(i,0,Len)F[i].im=-F[i].im;
		rep(i,0,Len)G[i]=(Z[i]-F[i])/2,F[i]=(Z[i]+F[i])/2
			,swap(G[i].re,G[i].im),G[i].im=-G[i].im;
	}
	else
	{
		rep(i,0,Len)F[i]=comp(F[i].re-G[i].im,F[i].im+G[i].re);
		FFT(F,-1);
		rep(i,0,Len)G[i]=comp(F[i].im,0),F[i]=comp(F[i].re,0);
	}
}

Ex4.快速子集变换

点我
和 FFTFFT 的形式差不多。代码和内容在这里。

Ex4.5快速子集运算

限于篇幅在另一篇文章作说明。

总集:帕秋莉的超级多项式

已知 n−1n-1 次多项式 FF 和 kk ，求
[(1+ln⁡(1+1exp⁡(∫1F)))k]′\left[\left(1+\ln(1+\frac{1}{\exp(\int\frac{1}{\sqrt{F}})})\right)^k\right]'
这道题汇集了几乎所有的多项式基础操作，有兴趣的可以尝试做一下，练练手。

代码:

CPP
#include<bits/stdc++.h>
#define Rep(i,a,b) for(register int i=(a);i<=(b);++i)
#define Repe(i,a,b) for(register int i=(a);i>=(b);--i)
#define rep(i,a,b) for(register int i=(a);i<(b);++i)
#define pb push_back
#define mx(a,b) (a>b?a:b)
#define mn(a,b) (a<b?a:b)
typedef unsigned long long uint64;
typedef unsigned int uint32;
typedef long long ll;
using namespace std;

namespace IO
{
	const uint32 Buffsize=1<<15,Output=1<<24;
	static char Ch[Buffsize],*S=Ch,*T=Ch;
	inline char getc()
	{
		return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++);
	}
	static char Out[Output],*nowps=Out;
	
	inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;}

	template<typename T>inline void read(T&x)
	{
		x=0;static char ch;T f=1;
		for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1;
		for(;isdigit(ch);ch=getc())x=x*10+(ch^48);
		x*=f;
	}

	template<typename T>inline void write(T x,char ch='\n')
	{
		if(!x)*nowps++='0';
		if(x<0)*nowps++='-',x=-x;
		static uint32 sta[111],tp;
		for(tp=0;x;x/=10)sta[++tp]=x%10;
		for(;tp;*nowps++=sta[tp--]^48);
		*nowps++=ch;
	}
}
using namespace IO;

void file(void)
{
	FILE*DSD=freopen("polynomial.in","r",stdin);
	FILE*CSC=freopen("polynomial.out","w",stdout);
}

const int MAXN=1<<21;

namespace poly
{
	const int mod=998244353,gen=3;
	
	static int g[23],iv[MAXN];

	inline int power(int u,int v)
	{
		register int sm=1;
		for(;v;v>>=1,u=(uint64)u*u%mod)if(v&1)
			sm=(uint64)sm*u%mod;
		return sm;
	}

	inline void predone()
	{
		Rep(i,1,21)g[i]=power(gen,(mod-1)>>i);
		iv[1]=1;
		Rep(i,2,1e5)iv[i]=mod-(uint64)mod/i*iv[mod%i]%mod;
	}

	static int Len,rev[MAXN];

	inline void calrev()
	{
		int II=log(Len)/log(2)-1;
		Rep(i,1,Len-1)rev[i]=rev[i>>1]>>1|(i&1)<<II;
	}

	inline int ad(int u,int v){return(u+=v)>=mod?u-mod:u;}

	inline void NTT(int*F,int typ)
	{
		Rep(i,1,Len-1)if(i<rev[i])swap(*(F+i),*(F+*(rev+i)));
		for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
		{
			register int wn=*(g+t);
			for(register int j=0;j<Len;j+=i)
			{
				register int w=1;
				rep(k,0,ii)
				{
					register int tt=(uint64)*(F+j+k+ii)*w%mod;
					*(F+j+k+ii)=ad(*(F+j+k),mod-tt);
					*(F+j+k)=ad(*(F+j+k),tt);
					w=(uint64)w*wn%mod;
				}
			}
		}
		if(typ==-1)
		{
			reverse(F+1,F+Len);
			register uint64 invn=power(Len,mod-2);
			rep(i,0,Len)*(F+i)=invn**(F+i)%mod;
		}
	}

	static int X[MAXN],Y[MAXN],Iv[MAXN];

	inline void mul(int *a,int *b,int *c,int lenl,int lenr)
	{
		if((ll)lenl*lenr<=300)
		{
			Rep(i,0,lenl+lenr)X[i]=0;
			Rep(i,0,lenl)Rep(j,0,lenr)
				X[i+j]=(X[i+j]+(ll)a[i]*b[j])%mod;
			Rep(i,0,lenl+lenr)c[i]=X[i];
			return;
		}
		for(Len=2;Len<=lenl+lenr;Len<<=1);
		calrev();
		Rep(i,0,lenl)X[i]=a[i];
		Rep(i,0,lenr)Y[i]=b[i];
		rep(i,lenl+1,Len)X[i]=0;
		rep(i,lenr+1,Len)Y[i]=0;
		NTT(X,1),NTT(Y,1);
		rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
		NTT(X,-1);
		Rep(i,0,lenl+lenr)c[i]=X[i];
		rep(i,lenl+lenr+1,Len)c[i]=0;
	}

	inline void Inv(int*F,int*G,int ln)
	{
		Iv[0]=power(F[0],mod-2);
		for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
		{
			rep(i,ln+1,Ln)F[i]=0;
			rep(i,0,Ln)X[i]=F[i],Y[i]=0;
			rep(i,0,(Ln>>1))Y[i]=Iv[i];
			Len=Ln,calrev();
			NTT(X,1),NTT(Y,1);
			rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
			NTT(X,-1);
			rep(i,0,(Ln>>1))X[i]=0;
			NTT(X,1);
			rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
			NTT(X,-1);
			rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
		}
		Rep(i,0,ln)G[i]=Iv[i];
	}

	static int ExX[MAXN],ExY[MAXN],Op[MAXN];

	inline void Div(int*F,int*G,int*Q,int*R,int lenf,int leng)
	{
		Rep(i,0,lenf)ExX[i]=F[lenf-i];
		Rep(i,0,leng)ExY[i]=G[leng-i];
		Rep(i,leng+1,lenf-leng)ExY[i]=0;
		Inv(ExY,ExY,lenf-leng);
		Rep(i,lenf-leng+1,lenf)ExX[i]=0;
		Rep(i,lenf-leng+1,leng)ExY[i]=0;
		mul(ExX,ExY,ExY,lenf-leng,lenf-leng);
		Rep(i,lenf-leng+1,(lenf-leng)<<1)ExY[i]=0;
		Rep(i,0,(lenf-leng)>>1)swap(ExY[i],ExY[lenf-leng-i]);
		mul(ExY,G,ExX,lenf-leng,leng);
		assert(F[leng]==ExX[leng]);
		Rep(i,0,leng-1)R[i]=ad(F[i],mod-ExX[i]);
		Rep(i,0,lenf-leng)Q[i]=ExY[i];
	}

	const int inv2=(mod+1)>>1;

	inline void Sqrt(int*F,int*G,int ln)
	{
		Op[0]=sqrt(F[0]);
		for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
		{
			mul(Op,Op,ExX,Ln>>1,Ln>>1);
			rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
			rep(i,Ln>>1,Ln)Op[i]=ExX[i]=0;
			Inv(Op,ExY,Ln>>1);
			Len=Ln,calrev();
			NTT(ExX,1),NTT(ExY,1);
			rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
			NTT(ExX,-1);
			rep(i,0,Ln>>1)Op[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
		}
		Rep(i,0,ln)G[i]=Op[i];
	}

	inline void Deriv(int*F,int*G,int ln)
	{Rep(i,1,ln)G[i-1]=(uint64)F[i]*i%mod;G[ln]=0;}

	inline void Inter(int*F,int*G,int ln)
	{Repe(i,ln,1)G[i]=(uint64)F[i-1]*iv[i]%mod;G[0]=0;}

	static int LnX[MAXN];

	inline void Ln(int*F,int*G,int ln)
	{
		Deriv(F,LnX,ln),Inv(F,G,ln);
		for(Len=2;Len<=ln<<1;Len<<=1);
		rep(i,ln+1,Len)LnX[i]=G[i]=0;
		calrev();
		NTT(LnX,1),NTT(G,1);
		rep(i,0,Len)G[i]=(uint64)G[i]*LnX[i]%mod;
		NTT(G,-1);
		Inter(G,G,ln);
	}

	inline void Exp(int *F,int *G,int ln)
	{
		Op[0]=1;
		for(register int Lx=2;Lx>>1<=ln;Lx<<=1)
		{
			rep(i,Lx>>1,Lx)Op[i]=0;
			Ln(Op,ExX,Lx);
			rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
			rep(i,0,Lx>>1)ExY[i]=Op[i];
			rep(i,Lx>>1,Lx)ExX[i]=ExY[i]=0;
			Len=Lx;
			calrev();
			NTT(ExY,1),NTT(ExX,1);
			rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
			NTT(ExX,-1);
			rep(i,0,Lx>>1)Op[i+(Len>>1)]=ExX[i];
		}
		Rep(i,0,ln)G[i]=Op[i];
	}
}
using poly::power;
using poly::Len;
using poly::calrev;
using poly::NTT;
using poly::mod;
using poly::predone;
using poly::Inv;
using poly::Inter;
using poly::Deriv;
using poly::Ln;
using poly::Exp;
using poly::Sqrt;
using poly::Div;
using poly::ad;

static int n,k,F[MAXN];

int main()
{
	file();
	predone();
	read(n),read(k);
	k%=mod;
	Rep(i,0,n-1)read(F[i]);
	Sqrt(F,F,n-1);
	Inv(F,F,n-1);
	Inter(F,F,n-1);
	Exp(F,F,n-1);
	Inv(F,F,n-1);
	++F[0];
	Ln(F,F,n-1);
	++F[0];
	Ln(F,F,n-1);
	Rep(i,0,n-1)F[i]=(ll)F[i]*k%mod;
	Exp(F,F,n-1);
	Deriv(F,F,n-1);
	Rep(i,0,n-1)write(F[i],' ');
	flush();
	return 0;
}
推荐这种东西写成一个 namespacenamespace ，方便调试和复板子。

评论

共 50 条评论，欢迎与作者交流。

正在加载评论...