【洛谷6633】[ZJOI2020] 抽卡(多项式毒瘤题)

点此看题面

  • 一个卡池中有\(n\)张卡,标号为\(a_{1\sim n}\)
  • 每轮随机从中抽取一张卡,求期望多少轮后能抽出\(k\)张编号连续的卡。
  • \(n\le2\times10^5\)

大毒瘤题,感觉自己真是啥也不会。

翻遍了全网题解依然有一些地方不是非常明白,对应地方的解释就先坑着或做个标记吧。

主要问题有两个,一是不太明白初始得出的生成函数为什么长这样,二是我自始至终就没有搞懂式子里的\(u\)到底是一个什么玩意。

初始的生成函数构造

假设我们已经抽出了\(w\)张不同的卡,抽出下一张不同的卡的概率就是\(\frac{n-w}n\),那么期望轮数就是\(\frac n{n-w}\),与具体抽出了哪些卡并没有关系。

因此我们只需要考虑抽出\(w\)张不同的卡之后仍然没有\(k\)张编号连续的卡的状态数即可。

我们先分开讨论每一个极长连续段,容易发现段与段之间是相互独立、互不干扰的。

对于一个长度为\(s\)的极长连续段,其中选出\(w\)张牌的方案数是:

\[[x^{n+1-w}](\sum_{i=0}^{k-1}x^i)^w \]

这里考虑的是我们在末尾添加一张一定不会被选中的辅助卡,那么每张没被选择的卡之前就可以放连续\(0\sim k-1\)张被选中的卡。

(但按照我的理解,感觉没被选中的卡数\(n+1-w\)才应该是生成函数的幂次,因为我们是考虑在每张没被选择的卡之前放选择的卡,而被选中的卡数\(w\)应该是最后需要的项的次数。不知道为什么题解里和我的想法是反着的,这里不太明白先做个标记。)

然后我们将没被选中的卡也加到生成函数当中,并添加一个辅助元\(u\)表示被选择的\(0\)的数量。

于是设\(G(x)=\sum_{i=1}^kx^i\),也就是此时每一块对应的生成函数。

由于\((n+1-w)+w=n+1\),那么我们所需要的就是:

\[[x^{n+1}]\sum(uG(x))^i=[x^{n+1}]\frac1{1-uG(x)} \]

(我也不知道\(u\)是什么意思,“表示被选择的\(0\)的数量”已经是我能找到的最详尽的解释了。)

扩展拉格朗日反演

首先给出拉格朗日反演和扩展拉格朗日反演的两个公式,设\(G(F(x))=x\),则满足:

\[[x^n]G(x)=\frac1n[x^{n-1}](\frac{x}{F(x)})^n\\ [x^n]H(G(x))=\frac1n[x^{n-1}]H'(x)(\frac{x}{F(x)})^n \]

这里我们需要用的是第二个式子,即扩展拉格朗日反演的式子,此时的\(H(x)=\frac1{1-ux}\)

直接套用公式得到:

\[[x^{n+1}]\frac1{1-uG(x)}=\frac1{n+1}[x^n]\frac{u}{(1-ux)^2}(\frac x{F(x)})^{n+1} \]

(注意,在这个式子中依然有\(u\)的存在,但我仍然不知道\(u\)的作用,代码中相关部分的实现也让我很迷惑。)

牛顿迭代求\(F(x)\)

感觉我对牛顿迭代掌握还不错,至少这一部分都能自己推出来。

现在我们的核心问题就只有求\(F(x)\)了,由于最终需要的是\(\frac x{F(x)}\),因此设\(f(x)=\frac{F(x)}x\)

由于\(G(F(x))=x\),所以:

\[\sum_{i=1}^{k}F(x)^k=\frac{F(x)^{k+1}-F(x)}{F(x)-1}=x \]

变形去分母,然后把所有项移到等号一边,得到:

\[F(x)^{k+1}-(1+x)F(x)+x=0 \]

代入\(f(x)\)得到:

\[x^{k+1}f(x)^{k+1}-(1+x)xf(x)+x=0 \]

发现这种式子可以使用牛顿迭代,它的公式是:

\[f(x)=f_0(x)-\frac{g(x)}{g'(x)} \]

这里的\(g(x)=x^{k+1}f_0(x)^{k+1}-(1+x)xf_0(x)+x\),代入得到:

\[f(x)=f_0(x)-\frac{x^{k+1}f_0(x)^{k+1}-(1+x)xf_0(x)+x}{x^{k+1}(k+1)f_0(x)^k-(1+x)x}=f_0(x)-\frac{x^{k}f_0(x)^{k+1}-(1+x)f_0(x)+1}{x^{k}(k+1)f_0(x)^k-(1+x)} \]

这样一来\(f(x)\)就可以求出了,然后只要再做一遍多项式求逆就能求出\(\frac x{F(x)}\)了。

分治\(NTT\)

发现我们现在求出的只是每一个极长连续段的生成函数。

因此我们最后还要再做一次分治\(NTT\)把这些生成函数全部卷起来。

然后根据最终得到的多项式计算答案就行了。

(我的代码不知道为什么跑得特别慢,是卡着时限过的。。。)

代码:\(O(nlog^2n)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 200000
#define X 998244353
#define IC(x,y) (1LL*IFac[x]*Fac[y]%X*Fac[(x)-(y)]%X)
using namespace std;
int n,k,a[N+5],s[N+5],f[N+5],ff[N+5];vector<int> F[N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
int Fac[N+5],IFac[N+5];I void Init_F()//预处理阶乘和阶乘逆元
{
	RI i;for(Fac[0]=i=1;i<=n;++i) Fac[i]=1LL*Fac[i-1]*i%X;
	for(IFac[n]=QP(Fac[n],X-2),i=n;i;--i) IFac[i-1]=1LL*Fac[i]*i%X;
}
namespace Poly//多项式板子
{
	#define Init(n) P=1,L=0;W(P<=2*(n)) P<<=1,++L;\
		for(i=0;i^P;++i) A[i]=B[i]=0,R[i]=((R[i>>1]>>1)|((i&1)<<L-1));
	#define Init_(n) P=1,L=0;W(P<=(n)) P<<=1,++L;\
		for(i=0;i^P;++i) A[i]=B[i]=0,R[i]=((R[i>>1]>>1)|((i&1)<<L-1));
	int PR=3,P,L,R[N<<2],A[N<<2],B[N<<2];
	I void NTT(int* s,CI op)//NTT
	{
		RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(x=s[i],s[i]=s[R[i]],s[R[i]]=x);
		for(i=1;i^P;i<<=1) for(U=QP(QP(PR,op),(X-1)/(i<<1)),j=0;j^P;j+=i<<1) for(S=1,k=0;
			k^i;++k,S=1LL*S*U%X) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%X))%X,s[i+j+k]=(x-y+X)%X;
	}
	I void Mul(CI n,int* a,int* b)//多项式乘法(针对数组)
	{
		RI i;Init(n);for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),i=0;i<=n;++i) a[i]=1LL*A[i]*t%X;
	}
	I void Mul(vector<int>& a,vector<int>& b)//多项式乘法(针对vector)
	{
		RI i,n=a.size()-1,m=b.size()-1;Init(n+m);
		for(i=0;i<=n;++i) A[i]=a[i];for(i=0;i<=m;++i) B[i]=b[i];
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),a.clear(),i=0;i<=n+m;++i) a.push_back(1LL*A[i]*t%X);
	}
	I void Inv(CI n,int* a,int* b)//多项式求逆
	{
		if(!n) return (void)(b[0]=QP(a[0],X-2));RI i;Inv(n>>1,a,b);
		Init(n);for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=(2LL*B[i]-1LL*A[i]*B[i]%X*B[i]%X+X)%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),i=0;i<=n;++i) b[i]=1LL*A[i]*t%X;
	}
	int p[N+5];I void Ln(CI n,int* a,int* b)//多项式Ln
	{
		RI i;for(i=0;i<=n;++i) b[i]=0;Inv(n,a,b);
		Init(n-1);for(i=0;i<=n-1;++i) A[i]=1LL*a[i+1]*(i+1)%X,B[i]=b[i];
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),b[0]=i=0;i^n;++i) b[i+1]=1LL*A[i]*t%X*QP(i+1,X-2)%X;
	}
	int q[N+5];I void Exp(CI n,int* a,int* b)//多项式Exp
	{
		if(!n) return (void)(b[0]=1);RI i;Exp(n>>1,a,b);
		Ln(n,b,p),Init(n);for(i=0;i<=n;++i) A[i]=b[i],B[i]=(!i-p[i]+a[i]+X)%X;
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),i=0;i<=n;++i) b[i]=1LL*A[i]*t%X;
	}
	int r[N+5];I void Pow(CI n,int* a,CI k)//多项式快速幂
	{
		Poly::Ln(n,a,r);for(RI i=0;i<=n;++i) a[i]=0,r[i]=1LL*r[i]*k%X;Poly::Exp(n,r,a);
	}
	int a[N+5],b[N+5],w[N+5];I void Newton(CI n,int* f)//牛顿迭代
	{
		if(!n) return (void)(f[0]=1);RI i;Newton(n>>1,f);
		for(i=0;i<=n;++i) a[i]=b[i]=w[i]=0;a[0]=1,b[0]=X-1,b[1]=X-1;//注意每次调用前清空
		for(i=0;i<=(n>>1);++i) a[i]=(a[i]-f[i]+X)%X,a[i+1]=(a[i+1]-f[i]+X)%X;
		for(i=0;i<=(n>>1);++i) w[i]=f[i];
		for(Pow(n,w,k),i=0;i<=n-k;++i) b[k+i]=(1LL*(k+1)*w[i]+b[k+i])%X;
		for(Mul(n,w,f),i=0;i<=n-k;++i) a[k+i]=(w[i]+a[k+i])%X;
		for(i=0;i<=n;++i) w[i]=0;Inv(n,b,w),Mul(n,a,w);//a,b分别表示分子和分母的多项式
		for(i=0;i<=n;++i) f[i]=(f[i]-a[i]+X)%X;//更新f数组
	}
	I void Solve(CI l,CI r)//分治NTT
	{
		if(l==r) return;RI i,mid=l+r>>1;Solve(l,mid),Solve(mid+1,r);Mul(F[l],F[mid+1]);
	}
}
int main()
{
	RI i;for(scanf("%d%d",&n,&k),Init_F(),i=1;i<=n;++i) scanf("%d",a+i);
	sort(a+1,a+n+1),Poly::Newton(n,ff),Poly::Inv(n,ff,f);//预处理
	RI c=0,t=1;for(a[n+1]=1e9,i=2;i<=n+1;++i) a[i]^(a[i-1]+1)&&(s[++c]=t,t=0),++t;//求出每个极长连续段长度
	RI j,v;for(i=1;i<=c;++i)//分别考虑每个极长连续段
	{
		for(j=0;j<=s[i];++j) ff[j]=f[j];Poly::Pow(s[i],ff,s[i]+1);
		for(v=QP(s[i]+1,X-2),F[i].push_back(1),//求出该极长连续段对应生成函数
			j=1;j<=s[i];++j) F[i].push_back(1LL*(s[i]-j+1)*ff[j]%X*v%X);//(这里的系数我也不太明白)
	}
	RI ans=0;for(Poly::Solve(1,c),i=0;i^n;++i)//枚举选出i张不同的卡片
		ans=(1LL*F[1][i]*n%X*QP(n-i,X-2)%X*IC(n,i)+ans)%X;//统计答案
		//n/(n-i)表示对期望轮数的贡献,1/C(n,i)即这i张不同的卡片是对应非法状态的概率
	return printf("%d\n",ans),0;//输出答案
}
posted @ 2021-01-15 21:17  TheLostWeak  阅读(34)  评论(0编辑  收藏