神仙DP题。

ORZRQY!

感谢_rqy的题解。如果以下内容有错误请指教。

题意

给出一个 $n \times 1001$ 的矩形,每个格子有 $q$ 的概率是安全的,要求选出一个与底边相邻且最大的安全矩形区域。求这个最大矩形区域面积恰好为 $k$ 的概率。

思路

建立状态转移方程

直接做根本做不了,考虑差分-前缀和之思想。

类似悬线法,我们把每根悬线加入状态。称下标 $n$ 对应的悬线长为从底部往上第一个危险格子下方的格子数目。

设状态 $f_k(i, j)$ 表示一个大小为 $i \times 1001$ 的泳池,底部所有的悬线长度 $\ge j$ ,且最大矩形面积 $< k$ 的概率。那么答案就是 $f_{k+1}(n, 0) - f_k(n, 0)$.(取 $>k$ 是为了避免负数下标)

决策是什么呢?先解决容易解决的边界情况。

$f(0, *) = 1; f(i, j) = 0 (\text{当 }ij \ge k)$,这是显然的。

再考虑一般情况的决策。我们不妨枚举从右往左第一个悬线长恰好为 $j$ 的位置 $p$ ,则其右边对应 $f(i-p, j+1)$。如图:

Pool

考虑到黄色位置左边仍然可能有位置悬线长恰好为 $j$ ,转移方程为: 直接求,复杂度是 $O(kn)$ 的,可以获得 70 pts.

矩阵快速幂优化

$n \le 10^9, k \le 100$,显然地在提示矩阵快速幂。考虑到我们只需要 $f(*, 0)$,那么就令 $g_i = f(i, 0)$.

对于 $g_1 \sim g_k$ 的部分直接暴力算,有用的状态不超过 $O(n \sqrt{n})$ 个。对于 $i \ge k$ 的情况,显然有 $f(i, j) = 0$ $(j \ne 0)$.

于是:

可以看出这个是常系数齐次线性递推。构造转移矩阵进行优化:

这样做的复杂度是 $O(k^3 \lg n)​$ ,可以获得 90pts.

利用多项式算法

详见利用多项式算法优化常系数齐次线性递推

这里取系数 $a_i = (1-q)f(i-1, 1)$ 即可。

实现

#include <bits/stdc++.h>
using namespace std;

namespace polynomial {
    // 多项式大全略
}

using namespace polynomial;
int n, k, x, y, q, ans, v[MAXN + 10], p[MAXN + 10], c[MAXN + 10];
map<pair<int, int>, int> opt;

int dp(int i, int j) {
    if(opt.count(make_pair(i, j)))
        return opt[make_pair(i, j)];
    else {
        if(i == 0) return opt[make_pair(i, j)] = 1;
        else if(i * j >= k || j > 1000) return opt[make_pair(i, j)] = 0;
        else {
            opt[make_pair(i, j)] = dp(i, j + 1);
            for(int p = 1; p <= i; p++)
                opt[make_pair(i, j)] = pls(opt[make_pair(i, j)], mul(dp(i - p, j + 1), mul(fpow(q, j), mul(dec(1, q), dp(p - 1, j)))));
            return opt[make_pair(i, j)];
        }
    }
}

void solve(int x, int ret[]) {
    static int a[MAXN + 10], tmp[MAXN + 10];
    memset(ret, 0, sizeof(int) * (MAXN + 10));
    ret[0] = 1;
    memset(a, 0, sizeof(a)), memset(tmp, 0, sizeof(tmp));
    a[1] = 1;
    while(x > 0) {
        if(x & 1) {
            conv(k-1, ret, a, ret);
            poly_div(2*k-2, k, ret, p, tmp, ret); // 模k次多项式,余式才为k-1次!!!
        }
        x >>= 1;
        conv(k-1, a, a, a);
        poly_div(2*k-2, k, a, p, tmp, a);
    }
}

void work(int factor) {
    for(int i = 0; i < k; i++) v[i] = dp(i, 0);
    p[k] = 1;
    for(int i = 1; i <= k; i++)
        p[k - i] = mul(dec(q, 1), dp(i - 1, 1));
    
    solve(n, c);
    for(int i = 0; i < k; i++)
        ans = pls(ans, mul(factor, mul(v[i], c[i])));

    opt.clear();
}

int main() {
    cin >> n >> k >> x >> y;
    q = mul(x, mod_inv(y));
    if(n <= 1000 && k <= 1000) {
        k++; ans = pls(ans, dp(n, 0)); opt.clear();
        k--; ans = dec(ans, dp(n, 0)); opt.clear();
    } else {
        k++; work(1);
        k--; work(MOD-1);
    }
    cout << ans;
    return 0;
}