shadowice1984 的博客

shadowice1984 的博客

题解 P4457 【[BJOI2018]治疗之雨】

posted on 2018-05-09 08:31:25 | under 题解 |

非常传统的高斯消元题目……

如果熟练之后可以很快的切掉这道题?


本体题解

其实在概率题中有一个非常传统的模型

给你一堆点,点u有$P(u,v)$的概率通向点v,求从某个点S出发,经过多少步之后可以第一次到达点T,通常有一些点进入之后就不可以继续出来,比如我们一般认为T点进去之后就不会出来了

然后我们发现这个东西如果不是DAG就不可能使用dp求解,因此我们一般使用高斯消元列方程的做法解决这类问题,一般来讲我们设$E_{i}$为从i出发到达点T的期望步数

那么我们可以得到这样的式子

$E_{i}=\sum_{j=1}^{n}E_{j}P_{i,j}+1$

$E_{T}=0$

意思就是这个点的期望步数为所有后继状态的期望步数乘上对应的转移概率之和(显然这符合期望的定义),最后无论如何转移步数都会+1(或者你可以认为是所有的后继状态的概率之和为1,然后我们本来这个1在Σ里面,我们提了出来)

然后就可以由n个变量,n个方程通过高斯消元求解了


那么对于这道题我们也可以认为是这个模型

我们可以把每一个血量值i作为一个点,这样我们现在只需要求出转移的概率$P_{i,j}$就可以利用高斯消元大力列方程求解了

发现从状态i可以转移到的状态只有0~i以及i+1这些状态

那么我们我们只需要求出这个英雄受到了a点伤害的概率就可以列方程了

那么这个英雄受到a点伤害的事件可以写成一个01序列,其中受到伤害记为1,不受伤害记为0,那么由于受到了a点伤害,所以有a个1,显然这样的01序列有$C_{k}^{a}$种,每种01序列出现的概率都是

$(\frac{1}{m+1})^{a}(\frac{m}{m+1})^{k-a}$

所以最后受到k点伤害的概率就是

$C_{k}^{a}(\frac{1}{m+1})^{a}(\frac{m}{m+1})^{k-a}$

写成阶乘的形式就是

$\frac{k!}{a!(k-a)!}(\frac{1}{m+1})^{a}(\frac{m}{m+1})^{k-a}$

显然这玩意可以在我们打表了逆元之后线性的递推出来

好了现在我们有了英雄受到k点伤害的概率,让我们开始列方程吧……

以下记$P_{k}$为英雄受到恰好k点伤害的概率,$E_{i}$为血量为i的英雄血量归0的期望步数

对于i!=n的情况我们都需要分两种情况讨论,这个英雄有没有被奶到……

所以我们从i转移到j的概率应该是

$\frac{m}{m+1}P_{i-j}+\frac{1}{m+1}P_{i-j-1}$

特别的,转移到i+1的前提是没有受到伤害,所以概率应该是

$\frac{1}{m+1}P_{0}$

特别的,对于i=n的情况我们没有被奶中的可能性,因此i转移到j的概率应该是

$P_{i-j}$

另外的一点是我们的方程并没有考虑受到过量伤害的问题……,但是我们发现这应该会影响的系数是$E_{0}$的系数,但是有一个问题,我们对于$E_{i}$的定义是血量为i的英雄血量归零的期望步数……,显然$E_{0}=0$因此我们无需考虑$E_{0}$前的系数……


所以最后我们列出来的方程大概长这样

对于i=1~n-1

$E_{i}=(\sum_{j=1}^{j}(\frac{m}{m+1}P_{i-j}+\frac{1}{m+1}P_{i-j-1})E_{j})+(\frac{1}{m+1}P_{0})E_{i+1}+1$

对于 i=n

$E_{n}=(\sum_{j=1}^{n}P_{n-j}E_{j})+1$

这样n-1个变量n-1个方程我们就可以解出$E_{m}$了,注意这里我们直接代入了$E_{0}=0$从而消掉了矩阵的一行一列

做完了?


等等,$n=1500$你跟我说高斯消元?

所以我们不能进行$O(n^3)$的高斯消元,我们观察一下这个矩阵,会发现它大概长这样

1 1 0 0 0 0 ……

1 1 1 0 0 0 ……

1 1 1 1 0 0 ……

1 1 1 1 1 0 ……

会发现它的主对角线之上只有一个变量的系数不为1,因此我们可以两行两行的向下消,最后削成这个样子

1 1 0 0 0 0 ……

0 1 1 0 0 0 ……

0 0 1 1 0 0 ……

0 0 0 1 1 0 ……

0 0 0 0 1 1 ……

然后我们发现最后一行没有那个多余的1,所以我们从最后一行开始倒着消元就可以把这个矩阵削成单位矩阵了~

复杂度$O(n^2)$

另外注意判断诸多无解的情况,比如什么k=0,m=0,k=1之类的

还有什么m=0的时候也需要注意判掉

上代码~

// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
using namespace std;const int N=1510;typedef long long ll;const ll mod=1e9+7;
inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod){if(p&1){r=r*a%mod;}}return r;}
int n;int p;ll m;ll k;ll mp[N][N];ll f[N];ll inv1;ll inv;ll res;int T;
inline void solve()
{
    scanf("%d%d%lld%lld",&n,&p,&m,&k);
    if(k==0){printf("-1\n");return;}//特判一些特殊情况
    if(m==0)
    {
        if(k==1){printf("-1\n");return;}
        for(;p>0;){if(p<n){p++;}p-=k;res++;}printf("%lld\n",res);res=0;return;
    }inv=po(m+1,mod-2);inv1=po(m,mod-2);f[0]=po(m*inv%mod,k);
    for(int i=1;i<=min((ll)(n+1),k);i++)f[i]=f[i-1]*inv1%mod*po(i,mod-2)%mod*(k-i+1)%mod;//递推f数组
    for(int i=1;i<n;i++)for(int j=1;j<=i;j++)mp[i][j]=(f[i-j]*m+f[i-j+1])%mod*inv%mod;//列方程
    for(int i=1;i<n;i++)mp[i][i+1]=f[0]*inv%mod;//i+1
    for(int i=1;i<n;i++)(mp[i][i]+=mod-1)%=mod;//1~i
    for(int i=1;i<=n;i++){mp[i][n+1]=mod-1;}//常数项
    for(int i=1;i<=n;i++){mp[n][i]=f[n-i];}(mp[n][n]+=mod-1)%=mod;
    for(int i=1;i<=n;i++)//高斯消元
    {
        ll inv=po(mp[i][i],mod-2);mp[i][i]=1;(mp[i][n+1]*=inv)%=mod;
        if(i!=n){(mp[i][i+1]*=inv)%=mod;}//两行两行的向下消
        for(int j=i+1;j<=n;j++)
        {
            ll mult=mp[j][i];mp[j][i]=0;(mp[j][i+1]+=mod-mult*mp[i][i+1]%mod)%=mod;
            (mp[j][n+1]+=mod-mult*mp[i][n+1]%mod)%=mod;
        }
    }
    for(int i=n;i>1;i--){(mp[i-1][n+1]+=mod-mp[i-1][i]*mp[i][n+1]%mod)%=mod;mp[i-1][i]=0;}//倒着消上去
    printf("%lld\n",mp[p][n+1]);
}
inline void clear()//清空
{
    for(int i=1;i<=min((ll)(n+1),k);i++){f[i]=0;}
    for(int i=1;i<=n;i++)for(int j=1;j<=n+1;j++)mp[i][j]=0;
}
int main(){scanf("%d",&T);for(int i=1;i<=T;i++){solve();clear();}return 0;}//拜拜程序~