题意
你现在有 (m+1) 个数：第一个为 (p) ，最小值为 (0) ，最大值为 (n) ；剩下 (m) 个都是无穷，没有最小值或最大值。你可以进行任意多轮操作，每轮操作如下：
在不为最大值的数中等概率随机选择一个（如果没有则不操作），把它加一；
进行 (k) 次这个步骤：在不为最小值的数中等概率随机选择一个（如果没有则不操作），把它减一。
现在问期望进行多少轮操作以后第一个数会变为最小值 (0)。
(1 leq p leq n leq 1500) ，(0 leq m, k leq 10^9) 。
Solution
显然我们只用考虑第一个数的变化。
设 (P_x) 表示一次操作 (-x) 概率，即 (k) 次中选出 (x) 次，剩余 (k-x) 次分配到其他 (m) 个数。
[P_x=frac{inom{k}{x}m^{k-x}}{(m+1)^k}
]
设 (Q_{x,y}) 表示一次操作从 (x) 变到 (y) 的概率，不难推出
[Q_{x,y}=egin{equation}
left{
egin{array}{lr}
0 & x=n 且 y=n+1 &\
P_{x-y} & x=n\
frac{1}{m+1}P_0 & y=x+1\
frac{m}{m+1}P_{x-y}+frac{1}{m+1}P_{x-y+1} & ext{otherwise}
end{array}
ight.
end{equation}
]
设 (f(i)) 表示从 (i) 变到 (0) 的期望操作数。
[egin{align}
& f(i)=1+sum_{j=1}^{i+1} Q_{i,j}f(j) & (1le i<n) \
& f(n)=1+sum_{j=1}^nQ_{n,j}f(j)
end{align}
]
我们可以 (n^2) 消元上述式子。将 ((1)) 变形可得
[egin{align}
f(i+1)=frac{f(i)-sum_{j=1}^iQ_{i,j}f(j)-1}{Q_{i,i+1}}
end{align}
]
可以利用 ((3),(4)) 得到 (f(n)) 关于 (f(1)) 的两个方程，解出 (f(1)) 后带入即可。
Code
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int Mod=1e9+7,N=1505;
inline int mul(int x, int y) { return 1ll*x*y%Mod; }
inline int po(int x, int y)
{
int r=1;
while(y)
{
if(y&1) r=mul(r,x);
x=mul(x,x), y>>=1;
}
return r;
}
int f[N],x[N],y[N],fac[N],inv[N],n,p,m,k,iv;
inline int calcp(int x, int y)
{
if(x==n&&y==n+1) return 0;
if(x==n) return f[x-y];
if(y==x+1) return mul(iv,f[0]);
return (mul(mul(m,iv),f[x-y])+mul(f[x-y+1],iv))%Mod;
}
int main()
{
int T; scanf("%d",&T);
fac[0]=inv[0]=1;
for(int i=1;i<=1500;++i) fac[i]=mul(i,fac[i-1]);
inv[1500]=po(fac[1500],Mod-2);
for(int i=1499;i;--i) inv[i]=mul(inv[i+1],i+1);
while(T--)
{
scanf("%d%d%d%d",&n,&p,&m,&k);
if(!k||(!m&&k==1))
{
puts("-1");
continue;
}
if(!m)
{
int ans=0;
while(p>0)
{
if(p<n) ++p;
p-=k,++ans;
}
printf("%d
",ans);
continue;
}
int inv0=po(m,Mod-2),inv1=po(po(m+1,k),Mod-2);
iv=po(m+1,Mod-2);
for(int i=0,j=1,g=po(m,k),l=min(n,k);i<=l;++i)
{
f[i]=mul(mul(mul(j,inv[i]),g),inv1);
j=mul(j,k-i),g=mul(g,inv0);
}
x[1]=1,y[1]=0;
for(int i=2;i<=n;++i)
{
x[i]=x[i-1],y[i]=(Mod+y[i-1]-1)%Mod;
for(int j=1;j<i;++j)
{
int tmp=calcp(i-1,j);
x[i]=(x[i]+Mod-mul(x[j],tmp))%Mod;
y[i]=(y[i]+Mod-mul(y[j],tmp))%Mod;
}
int tmp=po(calcp(i-1,i),Mod-2);
x[i]=mul(x[i],tmp),y[i]=mul(y[i],tmp);
}
int nx=0,ny=1;
for(int i=1;i<=n;++i)
{
int tmp=calcp(n,i);
nx=(nx+mul(x[i],tmp))%Mod;
ny=(ny+mul(y[i],tmp))%Mod;
}
int tmp=mul((ny-y[n]+Mod)%Mod,po((x[n]-nx+Mod)%Mod,Mod-2));
printf("%d
",(mul(tmp,x[p])+y[p])%Mod);
memset(f,0,sizeof(int)*(min(n,k)+1));
}
}