[BJOI2018]治疗之雨 [概率DP]

传送门

首先设k次中打中i次的概率为 p[i]

设 f[i] 为 还剩 i 滴血挂掉的期望

发现跟 f[i+1] 和 f[i] 有关, 已经可以高斯消元了, 但是 n <= 1500, 显然需要优化

好像有个高斯消元利用这道题的性质优化成 n^2 的, 但是我不会...

我们发现, 如果把第一个式子的 f[i+1] 提出来, f[i+1] 就只和f[i], f[i-1]... 有关, 就可以转移了

但是我们不知道 f[1], 类似套路, 每个 f[i] 都可以表示成 A[i] * f[1] + B[i], 然后转移

对于最后一个式子, 可以将 f[n] 表示为另一个有关 f[1]的式子, 两个方程解一下就可以得出f[1]

比较讨厌的是特殊情况, 需要特判一下情况

对于 , 需要但对做, 并在这个小类中还有特殊情况, 对于, 显然打不死

对于 

对于 

对于方程解出来系数为0 也要特判

另外, 由于数太大, 无法预处理, 所以组合数直接乘

---

#include<bits/stdc++.h>
#define N 50050
using namespace std;
const int Mod = 1000000007, up = 1500;
typedef long long ll;
ll mul(ll a, ll b){ return (a*b) % Mod;}
ll add(ll a, ll b){ return (a+b) % Mod;}
ll power(ll a, ll b){ 
	ll ans = 1; for(;b;b>>=1){if(b&1) ans = mul(ans, a); a = mul(a, a);} 
	return ans;
}
ll fac[N], inv[N];
struct Node{
	ll x, y; 
	Node(ll _x = 0, ll _y = 0){ x = _x, y = _y;}
	Node operator + (const Node &a){ return Node(add(x, a.x), add(y, a.y));}
	Node operator - (const Node &a){ return Node(add(x, Mod - a.x), add(y, Mod - a.y));}
	Node operator * (const ll &a){ return Node(mul(x, a), mul(y, a));}
}f[N];
int T, n, p, m, k; ll P[N];
ll C(int n, int m){ 
	ll up = 1, down = 1;
	for(int i=1; i<=m; i++) down = mul(down, i);
	for(int i=n; i>=n-m+1; i--) up = mul(up, i);
	return mul(up, power(down, Mod - 2));
}
ll Solve(){
	if(k == 0) return -1;
	if(k == 1 && n != 1) return -1;
	f[0] = Node(0, 0);
	ll x; Node tmp;
	for(int i=1; i<n; i++){
		f[i] = Node(0, 1) + f[max(i + 1 - k, 0)];
	} f[n] = Node(0, 1) + f[max(n - k, 0)];
	return f[p].y;
}
int main(){
	scanf("%d", &T);
	while(T--){
		scanf("%d%d%d%d", &n, &p, &m, &k);
		
		if(p == 0){ printf("0\n"); continue;}
		if(n == 1 && k == 0){ printf("-1\n"); continue;}
		
		ll invm = power(m + 1, Mod - 2);
		for(int i=0; i<=min(n, k); i++){
			P[i] = mul(C(k, i), mul(power(invm, i), power(mul(invm, m), k-i)));
		}
		if(m == 0){ printf("%lld\n", Solve()); continue;}
		if(P[0] == 0 || P[0] == 1){ printf("-1\n"); continue;}
		f[0] = Node(0, 0);
		f[1] = Node(1, 0); 
		ll invp = power(P[0], Mod - 2);
		for(int i=1; i<n; i++){
			f[i+1] = f[i] * (m + 1) - Node(0, m+1);
			Node tmp(0, 0);
			for(int j=0; j<=min(k, i); j++){
				tmp = tmp + f[i - j] * P[j];
			} f[i+1] = f[i+1] - tmp * m;
			for(int j=1; j<=min(k, i+1); j++){
				f[i+1] = f[i+1] - f[i - j + 1] * P[j];
			}
			f[i+1] = f[i+1] * invp;
		}
		Node First = f[n];
		Node Second(0, 0);
		for(int i=1; i<=min(n, k); i++) Second = Second + f[n - i] * P[i];
		Second.y++; Second = Second * power(add(1, Mod - P[0]), Mod - 2);
		ll A = add(Second.y, Mod-First.y), B = add(First.x, Mod - Second.x);
		if(B == 0 && f[p].x != 0){ printf("-1\n"); continue;}
		ll val = mul(A, power(B, Mod - 2));
		printf("%lld\n", add(mul(f[p].x, val), f[p].y));
	} return 0;
}