首先设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;
}