BJOI2018 治疗之雨 高斯消元 期望dp

<更新提示>

---

<正文>

治疗之雨

Description

你现在有$m+1$个数：第一个为$p$，最小值为$0$，最大值为$n$；剩下$m$个都是无穷，没有最小值或最大值。你可以进行任意多轮操作，每轮操作如下：

在不为最大值的数中等概率随机选择一个（如果没有则不操作），把它加一；

进行$k$次这个步骤：在不为最小值的数中等概率随机选择一个（如果没有则不操作），把它减一。

现在问期望进行多少轮操作以后第一个数会变为最小值$0$。

Input Format

输入包含多组数据。 输入第一行包含一个正整数$T$，表示数据组数。 接下来$T$行 ，每行$4$个非负整数 $n,p,m,k$（含义见题目描述），表示一次询问。

Output Format

输出$T$行，每行一个整数，表示一次询问的答案。

如果无论进行多少轮操作，第一个数都不会变为最小值$0$，那么输出$-1$；

否则，可以证明答案一定为有理数，那么请输出答案模$1000000007$的余数，即设答案为$\frac{a}{b}$（$a,b$为互质的正整数），你输出的整数为$x$，那么你需要保证 $0\le x<1000000007$且$a\equiv bxmod1000000007$。

Sample Input

2
2 1 1 1
2 2 1 1

Sample Output

6
8

Solution

首先，我们肯定考虑$dp$模型，设$f_i$代表第一个数值为$i$时变为$0$的期望操作次数，然后我们要考虑转移。

转移肯定要考虑当前轮数值的变化，但是我们并不知道这一轮数值会变小多少，因为概率是未知的。我们考虑预处理概率$p_i$，$p_i$代表一回合数值变小$i$的概率：

$$p_i=\frac{\left(\genfrac{}{}{0ex}{}{k}{i}\right)m^{k-i}}{(m+1{)}^k}$$

组合意义：$k$次减小作用在$m+1$个数上共有$(m+1{)}^k$种方案，强制现在$i$次作用在第一个数上，方案数就是$\left(\genfrac{}{}{0ex}{}{k}{i}\right)m^{k-i}$。

令$q_i$为一回合数值变小超过$i$的概率，然后我们考虑转移：

$$f_i=\frac{m}{m+1}\left(q_i+\sum _{j=0}^{i-1}{p}_j(f_{i-j}+1)\right)+\frac{1}{m+1}\left(q_{i+1}+\sum _{j=0}^i{p}_j(f_{i+1-j}+1)\right)$$

组合意义：有$\frac{m}{m+1}$的概率第一个数没有加一，然后我们枚举这一次第一个数变小了多少，根据全期望公式可知就是对应的概率乘上期望。反之，有$\frac{1}{m+1}$的概率第一个数加一了，同理根据全期望公式枚举这一次减少的多少即可。

然后考虑化简：

$$f_i=\frac{m}{m+1}\left(1+\sum _{j=0}^{i-1}{p}_j{f}_{i-j}\right)+\frac{1}{m+1}\left(1+\sum _{j=0}^i{p}_j{f}_{i+1-j}\right)$$

我们发现${\sum }_{j=0}^x{p}_j+q_{x+1}=1$，于是就可以取出求和号中的项消掉$q_i$和$q_{i+1}$，得到上式。

边界情况：$$f_n=1+\sum _{j=0}^{n-1}{p}_j{f}_{i-j}$$

发现转移成环，于是考虑高斯消元。$O(n^3)$的模意义高消可以得到$70$分。

我们注意这$n$个方程的增广矩阵，就是下三角矩阵多了一条斜线，每行向下消元，就只用减两个系数，可以实现$O(n^2)$消元。

值得注意的是，当消元时出现形容$0x=a$的情况，直接输出无解即可，不必要一开始就讨论什么情况下无解。

参考代码如下：

#include <bits/stdc++.h>
using namespace std;
const int N = 1520 , Mod = 1e9+7;
int n,m,p,k,inv[N],P[N],a[N][N],b[N],f[N];
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Add(int &a,int b) { a = add( a , b ); }
inline void Sub(int &a,int b) { a = sub( a , b ); }
inline void Mul(int &a,int b) { a = mul( a , b ); }
inline int quickpow(int a,int b) { int res = 1; for (b%=(Mod-1);b;Mul(a,a),b>>=1) if ( 1 & b ) Mul(res,a); return res; }
inline int inverse(int x) { return quickpow( x , Mod-2 ); }
inline void init(void)
{
    inv[0] = inv[1] = 1;
    for (int i=2;i<=1500;i++)
        inv[i] = mul( Mod - Mod/i , inv[Mod%i] );
}
inline void input(void) { scanf("%d%d%d%d",&n,&p,&m,&k); }
inline void build(void)
{
    memset( P , 0x00 , sizeof P );
    int val = inverse( quickpow(m+1,k) ) , C = 1;
    for (int i=0;i<=min(n,k);i++)
    {
        P[i] = mul( val , quickpow(m,k-i) );
        Mul( P[i] , C ) , Mul( C , mul( inv[i+1] , k-i ) );
    }
    int invm = inverse(m+1);
    memset( a , 0x00 , sizeof a );
    memset( b , 0x00 , sizeof b );
    for (int i=1;i<=n;i++) b[i] = 1;
    for (int i=1;i<n;i++)
    {
        for (int j=0;j<i;j++)
            Sub( a[i][i-j] , mul( mul( m , invm ) , P[j] ) ),
            Sub( a[i][i-j+1] , mul( invm , P[j] ) );
        Sub( a[i][1] , mul( invm , P[i] ) );
        Add( a[i][i] , 1 );
    }
    for (int i=0;i<n;i++) Sub( a[n][n-i] , P[i] );
    Add( a[n][n] , 1 );
}
inline bool Gauss(void)
{
    for (int i=1;i<=n;i++)
    {
        int inva = inverse(a[i][i]);
        if ( a[i][i] == 0 ) return false;
        for (int j=i+1;j<=n;j++)
        {
            int rate = mul( a[j][i] , inva );
            Sub( a[j][i+1] , mul( a[i][i+1] , rate ) );
            Sub( b[j] , mul( b[i] , rate ) );
        }
    }
    for (int i=n-1;i>=1;i--)
    {
        int rate = mul( a[i][i+1] , inverse( a[i+1][i+1] ) );
        Sub( b[i] , mul( b[i+1] , rate ) );
    }
    for (int i=1;i<=n;i++) f[i] = mul( b[i] , inverse( a[i][i] ) );
    return true;
}
int main(void)
{
    init();
    int T; scanf("%d",&T);
    while ( T --> 0 )
    {
        input();
        build();
        if ( !Gauss() ) puts("-1");
        else printf("%d\n",f[p]);
    }
    return 0;
}

---

<后记>