P4457 [BJOI2018]治疗之雨
题意描述比较迷惑，剩下 \(m\) 个数是无穷的意思是任意操作对它们不产生影响，操作影响的，或者说影响选哪个数进行操作的，只有第一个数，所以当前的状态可以用第一个数的值来表示。
设 \(dp_i\) 表示第一个数为 \(i\) 时把它变为 \(0\) 的期望步数，可得：
若 \(dp_p\) 无解则减不到 \(0\)，否则 \(dp_p\) 即为答案。
对于 \(dp_0\)，显然不需要任何操作。
对于 \(dp_n\) 的转移，枚举第一个数一共减了 \(j\) 次以及具体在哪几次操作中减了第一个数，对于每一个减了 \(j\) 次的方案，有 \((\frac{1}{m+1})^j \times (\frac{m}{m+1})^{k-j}\) 的概率出现，又因一共有 \(\binom{k}{j}\) 种方案，所以把 \(n\) 减为 \(n-j\) 的概率为 \(\binom{k}{j} \times (\frac{1}{m+1})^j \times (\frac{m}{m+1})^{k-j}\)。
对于 \(dp_i\) 的转移，类似于 \(dp_n\)。第一个数 \(i\) 有 \(\frac{1}{m+1}\) 的概率变为 \(i+1\)，有 \(\frac{m}{m+1}\) 的概率不变，分类以后两种状态的转移与 \(dp_n\) 的转移相同。
需注意，由于第一个数减为 \(0\) 后不能再减小，设当前第一个数为 \(i\)，那么最多进行 \(min(i, k)\) 次减 \(1\) 操作。
\(\text{dp}\) 状态和转移设计好了，接下来思考怎么转移。发现转移有后效性，再看一眼数据范围，\(1 \leq p \leq n \leq 1500\)，似乎可以用高斯消元？
但正常高斯消元复杂度是 \(O(n^3)\) 的，\(O(T \cdot n^3)\) 过不了。/kk
但观察转移方程或打印一下增广矩阵发现矩阵都是形如这样的：
\( \begin{bmatrix} k_{1,1} & k_{1,2} & ... & 0 & 0 & k_{1,n+1}\\ k_{2,1} & k_{2,2} & k_{2,3} & ... & 0 & k_{2,n+1}\\ ... & ... & ... & ... & ... & ...\\ 0 & k_{n-1,n-1-min(n-1,k)} & k_{n-1,n-min(n-1,k)} & ... & k_{n-1,n} & k_{n-1,n+1}\\ 0 & 0 & k_{n,n-min(n,k)} & ... & k_{n,n} & k_{n,n+1} \end{bmatrix} \)
看得出矩阵的左下和右上会存在一堆 \(0\)。当 \(k \ge n-1\) 时，共有右上角 \(\frac{(n-1)(n-2)}{2}\) 个 \(0\)。从上往下消元，当前行最多只有 \(2\) 个数，于是高斯消元的复杂度可以降到 \(O(n^2)\)。
具体实现时，由于 \(k\) 很大，但是是固定的，可以递推预处理出所有的概率系数。
设：
展开得：
\(f(x)\) 递推式为：
起始项：
还有一些小细节，\(k=0\) 但 \(p \ne 0\) 时无解，\(m=0\) 时直接计算，高斯消元第二层函数里不能有求逆元操作，否则复杂度多一个 \(\text{log}\)。
整体复杂度 \(O(T \cdot n^2)\)
点击查看代码
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll SIZE = 1505;
const ll mod = 1000000007;
ll a[SIZE][SIZE];
ll f[SIZE];
ll n, st, m, K;
ll T;
ll power(ll x, ll y){
ll jl = 1;
while(y){
if(y & 1) jl = (jl * x) % mod;
x = (x * x) % mod;
y >>= 1;
}
return jl;
}
int main(){
scanf("%lld", &T);
while(T--){
scanf("%lld%lld%lld%lld", &n, &st, &m, &K);
//特判
if(K == 0){
printf("-1\n");
continue;
}
if(m == 0){
if(K <= 1) printf("-1\n");
else{
ll cnt = 0;
while(st > 0){
if(st < n) st++;
st -= K;
cnt++;
}
printf("%lld\n", cnt);
}
continue;
}
//构造增广矩阵
ll t = power(m+1, mod-2), tt = ((1-t)%mod+mod)%mod;
for(ll i = 1; i <= n; i++){
for(ll j = 1; j <= n+1; j++){
a[i][j] = 0;
}
f[i] = 0;
}
ll xx;
f[0] = power((m * t) % mod, K);
for(ll i = 1; i <= n; i++){
f[i] = ((((K-i+1+mod)%mod) * f[i-1])%mod * power((m*i)%mod, mod-2)) % mod;
}
for(ll i = 1; i < n; i++){
a[i][i] = (-1+mod)%mod;
for(ll j = 0; j <= min(i+1, K); j++){
xx = (f[j] * t) % mod;
a[i][i+1-j] = (a[i][i+1-j] + xx) % mod;
}
for(ll j = 0; j <= min(i, K); j++){
xx = (f[j] * tt) % mod;
a[i][i-j] = (a[i][i-j] + xx) % mod;
}
a[i][n+1] = (-1+mod)%mod;
}
a[n][n] = mod-1;
for(ll j = 0; j <= min(n, K); j++){
a[n][n-j] = (a[n][n-j] + f[j]) % mod;
}
a[n][n+1] = mod-1;
//高斯消元
for(ll i = 1; i <= n; i++){
ll uu = power(a[i][i], mod-2);
for(ll j = i+1; j <= n; j++){
xx = (a[j][i] * uu) % mod;
for(ll k = i; k <= n; k++){
if(a[i][k] == 0) break;
a[j][k] = (a[j][k] - ((xx * a[i][k]) % mod) + mod) % mod;
}
a[j][n+1] = (a[j][n+1] - ((xx * a[i][n+1]) % mod) + mod) % mod;
}
}
for(ll i = n-1; i >= 1; i--){
xx = (a[i][i+1] * power(a[i+1][i+1], mod-2)) % mod;
a[i][i+1] = 0;
a[i][n+1] = (a[i][n+1] - ((xx*a[i+1][n+1])%mod) + mod) % mod;
}
//判无解
bool ff = 1;
for(int i = 1; i <= n; i++){
if(a[i][i] == 0 && a[i][n+1] != 0){
ff = 0;
break;
}
}
if(!ff) printf("-1\n");
else{
ll ans = (a[st][n+1] * power(a[st][st], mod-2)) % mod;
printf("%lld\n", ans);
}
}
return 0;
}
浙公网安备 33010602011771号