bzoj5292:[Bjoi2018]治疗之雨

传送门

首先设\(m\)为随从个数，\(k\)为暗影打击装甲的个数，\(p\)为剩余生命值，\(n\)为生命上限

然后考虑每个回合受到伤害，设\(A_i\)为每个回合被攻击\(i\)次的概率

\[A_i=C^{i}_{k}*（\frac{1}{m+1}）^i*(\frac{m}{m+1})^{k-i}\\ \]

然后考虑答案，我们设\(F_i\)为当前有\(i\)滴血被干掉的期望回合数，再设\(g_i\)为当前回合打出超过\(i\)滴血的概率，容易得出

\[F_i=\frac{1}{m+1}(g_{i+1}+\sum_{j=0}^{i}A_{j}(F_{i+1-j}-1))+\frac{m}{m+1}(g_i+\sum_{j=0}^{i-1}A_j(F_{i-j}+1))\\ F_i=\frac{1}{m+1}(1+\sum_{j=0}^{i}A_{j}F_{i+1-j})+\frac{m}{m+1}(1+\sum_{j=0}^{i-1}A_jF_{i-j})\\ \]

好像没有什么问题，但是还得考虑满血的情况，这个时候是回不上血的

所以最后得出来

\[i!=n:\\ F_i=\frac{1}{m+1}(1+\sum_{j=0}^{i}A_{j}F_{i+1-j})+\frac{m}{m+1}(1+\sum_{j=0}^{i-1}A_jF_{i-j})\\ i==n:\\ F_n=1+\sum_{j=0}^{n-1}A_jF_{n-j} \]

然后考虑概率期望的通常做法：高斯消元

由于\(n=1500\)，高斯消元显然会TLE

然后我们观察一下这个方程组，发现第\(i\)个方程只有前\(i+1\)个系数不为\(0\)，所以考虑每次消元都只留下\(i,i+1\)，然后最后一个方程只有第\(n\)个位置上系数不为\(0\)，然后解出\(F_n\)，递推出答案就好了，那么复杂度就可以优化到\(O(n^2)\)了

提示：bzoj太慢，mod太多会TLE，请优化

代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
#include<map>
using namespace std;
void read(int &x) {
    char ch; bool ok;
    for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
    for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
#define rg register
const int maxn=1.5e3+10,mod=1e9+7;
int n,T,m,p,k,a[maxn][maxn],fac[maxn],inv[maxn],A[maxn],facc[maxn];
int mul(int x,int y){return 1ll*x*y-1ll*x*y/mod*mod;}
int mi(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        b>>=1,a=1ll*a*a%mod;
    }
    return ans;
}
void prepare()
{
    fac[0]=inv[0]=1;
    for(rg int i=1;i<=1500;i++)fac[i]=1ll*fac[i-1]*i%mod;
    inv[1500]=mi(fac[1500],mod-2);;
    for(rg int i=1499;i;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod;
}
int C(int x,int y)
{
    int ans=1;
    if(y<=1500)ans=1ll*ans*fac[y]*inv[y-x]%mod;
    else ans=1ll*ans*facc[x]%mod;
    return 1ll*ans*inv[x]%mod;
}
void gauss()
{
    int y[maxn];
    for(rg int i=1;i<=n;i++)
    {
        for(rg int j=1;j<i;j++)
        {
            int t=1ll*a[i][j]*y[j]%mod;
            a[i][j+1]=(a[i][j+1]-mul(t,a[j][j+1])+mod)%mod;
            a[i][j]=0,a[i][n+1]=(a[i][n+1]-mul(t,a[j][n+1])+mod)%mod;
        }
        y[i]=mi(a[i][i],mod-2);
    }
    int ans=mul(a[n][n+1],y[n]);
    for(rg int i=n-1;i>=p;i--)ans=1ll*(a[i][n+1]-mul(a[i][i+1],ans)+mod)%mod*y[i]%mod;
    printf("%d\n",ans?(ans+mod)%mod:-1);
}
int main()
{
    read(T),prepare();
    while(T--)
    {
        memset(a,0,sizeof a);
        read(n),read(p),read(m),read(k);int res=0;
        if(!k||(!m&&k==1)){puts("-1");continue;}
        if(!m)
        {
            while(p>0){if(p<n)++p;p-=k,++res;}
            printf("%d\n",res),res=0;continue;
        }
        for(rg int i=1;i<=n;i++)a[i][i]=a[i][n+1]=mod-1;
        int ny=mi(m+1,mod-2),u=mul(ny,m);memset(A,0,sizeof A);facc[0]=1;
        if(k>1500)for(rg int i=1;i<=n;i++)facc[i]=mul(facc[i-1],(k-i+1));
        for(rg int i=0;i<=min(k,n);i++)A[i]=mul(C(i,k),mul(mi(ny,i),mi(u,k-i)));
        for(rg int i=1;i<n;i++)
        {
            for(rg int j=0;j<=i;j++)a[i][i+1-j]=(a[i][i+1-j]+mul(ny,A[j]))%mod;
            for(rg int j=0;j<i;j++)a[i][i-j]=(a[i][i-j]+mul(u,A[j]))%mod;
        }
        for(rg int i=0;i<n;i++)a[n][n-i]=(a[n][n-i]+A[i])%mod;
        gauss();
    }
}