【洛谷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\)张牌的方案数是：
这里考虑的是我们在末尾添加一张一定不会被选中的辅助卡，那么每张没被选择的卡之前就可以放连续\(0\sim k-1\)张被选中的卡。
（但按照我的理解，感觉没被选中的卡数\(n+1-w\)才应该是生成函数的幂次，因为我们是考虑在每张没被选择的卡之前放选择的卡，而被选中的卡数\(w\)应该是最后需要的项的次数。不知道为什么题解里和我的想法是反着的，这里不太明白先做个标记。）
然后我们将没被选中的卡也加到生成函数当中，并添加一个辅助元\(u\)表示被选择的\(0\)的数量。
于是设\(G(x)=\sum_{i=1}^kx^i\)，也就是此时每一块对应的生成函数。
由于\((n+1-w)+w=n+1\)，那么我们所需要的就是：
（我也不知道\(u\)是什么意思，“表示被选择的\(0\)的数量”已经是我能找到的最详尽的解释了。）
扩展拉格朗日反演
首先给出拉格朗日反演和扩展拉格朗日反演的两个公式，设\(G(F(x))=x\)，则满足：
这里我们需要用的是第二个式子，即扩展拉格朗日反演的式子，此时的\(H(x)=\frac1{1-ux}\)。
直接套用公式得到：
（注意，在这个式子中依然有\(u\)的存在，但我仍然不知道\(u\)的作用，代码中相关部分的实现也让我很迷惑。）
牛顿迭代求\(F(x)\)
感觉我对牛顿迭代掌握还不错，至少这一部分都能自己推出来。
现在我们的核心问题就只有求\(F(x)\)了，由于最终需要的是\(\frac x{F(x)}\)，因此设\(f(x)=\frac{F(x)}x\)。
由于\(G(F(x))=x\)，所以：
变形去分母，然后把所有项移到等号一边，得到：
代入\(f(x)\)得到：
发现这种式子可以使用牛顿迭代，它的公式是：
这里的\(g(x)=x^{k+1}f_0(x)^{k+1}-(1+x)xf_0(x)+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;//输出答案
}
