「BJOI2018」治疗之雨（概率+高斯消元转递推）

https://loj.ac/problem/2513

这一类问题做多了现在看到都是秒。

(O(n^2))预处理(g[i])表示k轮后，第一个数恰好少了(i)的概率。

设(f[i])表示(i)的期望经过次数，那么(f[i]=[i=p]+sum_{j=i-1}^n f[j]*系数(j->i))。

这个系数可以用(g)来(O(1))算出，高斯消元后，(Ans=sum_{i=1}^n f[i])

根据树上随机游走那一套递推，我们从右往左做，可以一直得到(f[i]=k*f[i-1]+b)的形式

到(i=1)时，因为(f[0]=0)，所以(f[1]=b)，然后又从可以左往右推出(f[i])的值。

时间复杂度：(O(T*n^2))

Code

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int mo = 1e9 + 7;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

const int N = 1505;

int T;
ll n, p, m, k;
ll m1;

ll g[N];

ll C(ll n, ll m) {
	if(n < m) return 0;
	ll s = 1;
	fo(i, n - m + 1, n) s = s * i % mo;
	ll v = 1;
	fo(i, 1, m) v = v * i % mo;
	return s * ksm(v, mo - 2) % mo;
}

ll calc(int i, int j) {
	if(j <= 0) return 0;
	ll s = 0;
	ll v = j == n ? 1 : (m * m1 % mo);
	if(i <= j) {
		v = v * g[j - i] % mo;
		s = v;
	}
	v = j == n ? 0 : m1;
	if(v) {
		v = v * g[j + 1 - i] % mo;
		s = (s + v) % mo;
	}
	return s;
}

struct P {
	ll x, y;
	P(ll _x = 0, ll _y = 0) {
		x = _x, y = _y;
	}
};

P operator * (P a, P b) {
	return P(a.x * b.x % mo, (a.y + b.y * a.x) % mo);
}
P operator * (P a, ll b) {
	return P(a.x * b % mo, a.y * b % mo);
}

P f[N];
ll s[N];

void work() {
	scanf("%lld %lld %lld %lld", &n, &p, &m, &k);
	if(k == 0) {
		pp("-1
"); return;
	}
	if(m == 0 && k == 1 && n > 1) {
		pp("-1
"); return;
	}
	m1 = ksm(m + 1, mo - 2);
	fo(i, 0, n) {
		g[i] = C(k, i) * ksm(m1, i) % mo * ksm(m1 * m % mo, k - i) % mo;
	}
	fd(i, n, 1) {
		ll a0 = 1, b0 = 0;
		a0 = (a0 - calc(i, i)) % mo;
		P c = P(1, 0);
		fo(j, i + 1, n) {
			c = f[j] * c;
			P d = c * calc(i, j);
			a0 = (a0 - d.x) % mo;
			b0 = (b0 + d.y) % mo;
		}
		if(i == p) b0 ++;
		f[i] = P(calc(i, i - 1), b0) * ksm(a0, mo - 2);
	}
	ll ans = 0;
	fo(i, 1, n) {
		s[i] = (s[i - 1] * f[i].x + f[i].y) % mo;
		ans = (ans + s[i]) % mo;
	}
	ans = (ans % mo + mo) % mo;
	pp("%lld
", ans);
}

int main() {
	scanf("%d", &T);
	fo(ii, 1, T) {
		work();
	}
}