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

本文介绍了一个复杂的概率问题“治疗之雨”,通过动态规划模型和高斯消元法求解第一个数变为最小值0的期望操作次数。详细解析了转移概率的计算、DP状态转移方程的推导,以及如何利用高斯消元法高效解决该问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

<更新提示>


<正文>

治疗之雨

Description

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

在不为最大值的数中等概率随机选择一个(如果没有则不操作),把它加一;

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

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

Input Format

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

Output Format

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

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

否则,可以证明答案一定为有理数,那么请输出答案模 1000000007 1000000007 1000000007的余数,即设答案为 a b \frac{a}{b} ba a , b a,b a,b为互质的正整数),你输出的整数为 x x x,那么你需要保证 0 ≤ x < 1000000007 0 \leq x < 1000000007 0x<1000000007 a ≡ b x   m o d   1000000007 a \equiv bx\ mod\ 1000000007 abx mod 1000000007

Sample Input
2
2 1 1 1
2 2 1 1
Sample Output
6
8
Solution

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

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

p i = ( k i ) m k − i ( m + 1 ) k p_i=\frac{\binom{k}{i}m^{k-i}}{(m+1)^k} pi=(m+1)k(ik)mki

组合意义: k k k次减小作用在 m + 1 m+1 m+1个数上共有 ( m + 1 ) k (m+1)^k (m+1)k种方案,强制现在 i i i次作用在第一个数上,方案数就是 ( k i ) m k − i \binom{k}{i}m^{k-i} (ik)mki

q i q_i qi为一回合数值变小超过 i i i的概率,然后我们考虑转移:

f i = m m + 1 ( q i + ∑ j = 0 i − 1 p j ( f i − j + 1 ) ) + 1 m + 1 ( q i + 1 + ∑ j = 0 i p j ( f i + 1 − j + 1 ) ) 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}^ip_j(f_{i+1-j}+1)\right) fi=m+1m(qi+j=0i1pj(fij+1))+m+11(qi+1+j=0ipj(fi+1j+1))

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

然后考虑化简:

f i = m m + 1 ( 1 + ∑ j = 0 i − 1 p j f i − j ) + 1 m + 1 ( 1 + ∑ j = 0 i p j f i + 1 − j ) 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}^ip_jf_{i+1-j}\right) fi=m+1m(1+j=0i1pjfij)+m+11(1+j=0ipjfi+1j)

我们发现 ∑ j = 0 x p j + q x + 1 = 1 \sum_{j=0}^xp_j+q_{x+1}=1 j=0xpj+qx+1=1,于是就可以取出求和号中的项消掉 q i q_i qi q i + 1 q_{i+1} qi+1,得到上式。

边界情况: f n = 1 + ∑ j = 0 n − 1 p j f i − j f_n=1+\sum_{j=0}^{n-1}p_jf_{i-j} fn=1+j=0n1pjfij

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

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

值得注意的是,当消元时出现形容 0 x = a 0x=a 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;
}


<后记>

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值