- 一个卡池中有(n)张卡，标号为(a_{1sim n})。
- 每轮随机从中抽取一张卡，求期望多少轮后能抽出(k)张编号连续的卡。
- (nle2 imes10^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
]
这里考虑的是我们在末尾添加一张一定不会被选中的辅助卡，那么每张没被选择的卡之前就可以放连续(0sim 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
",ans),0;//输出答案
}