BJOI2018 治疗之雨

在洛谷上发了份题解，链接：https://www.luogu.org/blog/user22151/solution-p4457

我不是很会矩阵做这道题(不会卡常)。

但是我们可以用本质差不多的方法，但是不用矩阵，可能常数会小一点？反正目前洛谷开O2不开O2、bz都是最快的……

用$dp[i]$表示还剩$i$点血的时候，期望还要有多少轮操作。

用$f[x]$表示$K$次中有$i$次打到英雄的概率。

$f[x]=C(K,x)(frac{1}{m+1})^x(frac{m}{m+1})^{(k-x)}$

首先对于$m=0$的情况特判。

对于$m!=0$情况，我们发现有

$dp[0]=0$

$dp[1]=a_{1,2}dp[2]+a_{1,1}dp[1]+1$

...

$dp[i]=a_{i,i+1}dp[i+1]+a_{i,i}dp[i]+...+a_{i,1}dp[1]+1$

...

$dp[n]=a_{n,n}dp[n]+a_{n,n-1}dp[n-1]+...+a_{n,1}dp[1]+1$

其中$a_{i,j}$表示英雄$i$点血在一轮操作后变成$j$点血的概率，这个可以根据$f$算出来，不再赘述。

所以说我们可以得到$n$个等式，第$i$个式子是$dp[i]$等于多少多少

$n$个未知数，$n$个等式(方程)，很优秀的样子。

但是我们发现，第$i$个等式右边有$dp[i+1]$，所以不能直接递推，感觉很恼火。

我们把第$i$个式子右边的$dp[i+1]$放到等式左边来，然后把$dp[i]$放到等式右边，那么这个式子就变成了$dp[i+1]$等于多少多少。

注意第$n$个式子例外，因为和$dp[n+1]$无关，仍是$dp[n]$等于多少多少。

所以这样我们可以用第$i$个式子算$dp[i+1]$？

但是我们算不出$dp[1]$呀，也没法直接往后递推。

那我们先设$dp[1]=X$，那么我们可以根据第$1$~$n-1$个式子，把$dp[i]$表示成$AX+B$的形式，递推算出$dp[n]=A_1X+B_1$，然后根据第$n$个式子算出$dp[n]=A_2X+B_2$。

那么$A_1X+B_1=A_2X+B_2$，解一元一次方程直接算出$X$。

注意特判解不出$X$的情况输出-1。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(register int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(register int i=(a);i>=(b);--i)
const int maxn=3000+7;
const ll mod=1e9+7;
ll Td,n,P,m,K,X;

char cc; ll ff;
template<typename T>void read(T& aa) {
    aa=0;cc=getchar();ff=1;
    while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
    if(cc=='-') ff=-1,cc=getchar();
    while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
    aa*=ff;
}

ll qp(ll x,ll k) {
    if(x<=1) return x;
    ll rs=1;
    while(k) {
        if(k&1) rs=rs*x%mod;
        k>>=1; x=x*x%mod;
    }
    return rs;
}

ll mo(ll x) {if(x>=mod) x-=mod; return x;}

struct T{
    ll x,y;
    T(ll x=0.,ll y=0.):x(x),y(y){}
    inline T operator + (const T& b) const{return T(mo(x+b.x),mo(y+b.y));}
    inline T operator - (const T& b) const{return T(mo(x-b.x+mod),mo(y-b.y+mod));}
    inline T operator * (const T&b) const{return T((x*b.y+y*b.x)%mod,y*b.y%mod);}
}dp[maxn],f[maxn],o,r,t,inv,finv;

ll work() {
    if(K==0) return -1;
    if(K==1&&n!=1) return -1;
    dp[0]=T(0,0);
    ll x; T p;
    For(i,1,n-1) {
        dp[i]=T(0,1)+dp[max(i+1-K,(ll)0)];
    }
    dp[n]=T(0,1)+dp[max(n-K,(ll)0)];
    return dp[P].y;
}	

ll solve() {
    if(P==0) return 0;
    if(n==1&&K==0) return -1;
    ll x; T p;
    o=T(0,m+1); r=T(0,m); inv=T(0,qp(qp(m+1,K),mod-2));
    f[0]=T(0,1); x=min(K,n);
    For(i,1,x) f[i]=f[i-1]*T(0,(K-i+1)*qp(i,mod-2)%mod);
    For(i,0,x) f[i]=f[i]*inv*T(0,qp(m,K-i)); 
    if(m==0) return work();
    finv=T(0,qp(f[0].y,mod-2));
    if(f[0].y==1||f[0].y==0) return -1;
//	cerr<<clock()<<"
";
    //f[i]: probability of get i of K , finv : inv of f[0]
    dp[0]=T(0,0);
    dp[1]=T(1,0);
    //dp[i+1]=((m+1)(dp[i]-1)-m*sum(f[j]*dp[i-j])-sum(f[j]*dp[i-j+1]))/f[0]
    For(i,1,n-1) {
        dp[i+1]=o*dp[i]-o;
        x=min((ll)i,K);
        p=T(0,0);
        For(j,0,x) p=p+f[j]*dp[i-j];
        dp[i+1]=dp[i+1]-p*r;
        x=min((ll)i+1,K);
        For(j,1,x) dp[i+1]=dp[i+1]-f[j]*dp[i-j+1];
        dp[i+1]=dp[i+1]*finv;
    }
//	cerr<<clock()<<"
";
    //(1-f[0])*dp[n]=sum(f[j]*dp[n-j])+1;
    t=T(0,1); x=min(n,K);
    For(j,1,x) t=t+f[j]*dp[n-j];
    t=t*T(0,qp(mo(mod+1-f[0].y),mod-2));
    t=t-dp[n];
    if(t.x==0&&dp[P].x!=0) return -1;
    if(t.x) X=(mod-1)*t.y%mod*qp(t.x,mod-2)%mod;
    return (X*dp[P].x%mod+dp[P].y)%mod;
}

int main() {
    read(Td);
    while(Td--) {
        read(n); read(P); read(m); read(K);
        printf("%lld
",solve());
    }
//	cerr<<clock()<<"
";
    return 0;
}