Luogu3824 NOI2017 泳池

题面：luogu3824

题意：宽为$n$米，高为$1001$米的长方形网格，每个格子有$q$的概率是安全的。问以最下面一行为底边的最大的安全子矩阵的大小恰好为$K$的概率。$n\le 10^9$，$K\le 1000$

题解：首先转化问题为求矩阵大小小于等于$K$的概率。答案即为$\le K$减去$\le K-1$

设$dp[i][j]$表示宽度为$i$，高度至少为$j$（从下往上数$j$行），最大子矩阵大小小于等于$K$的概率。

考虑从$j+1$行转移：一种是第$j+1$行都可选；一种是第$j+1$行有不可选的位置，那么枚举第一个不可选的位置。

$dp[i][j]=dp[i][j+1]\times q^i+{\sum }_{k=1}^i{q}^{k-1}(1-q)\times dp[k-1][j+1]\times dp[i-k][j]$

边界：$dp[0][i]=1$，$dp[i][j]=0,i\times j>K$

要求的即为$dp[n][0]$

因为$i\times j\le K$，复杂度为$O(n\cdot k\log k)$，70pts。

考虑当$n$很大的时候，最下面一行每隔$K$个格子一定会有不可选的位置。

设$f[i]$表示当$i>K$时矩阵大小小于等于$K$的概率。

枚举第一个不可选的位置转移

$f[i]={\sum }_{j=1}^{K+1}f[i-j]\times q^{j-1}(1-q)\times dp[j-1][1]$，$i>K$

可以发现这是一个常系数齐次递推式，可以用特征多项式优化至$O(k^2\cdot \log n)$。

(或丧心病狂地用分治NTT + 多项式取模做到$O(k\log k\cdot \log n)$)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll ;
const int mod = 998244353 ;
inline ll power (ll x, ll y) {
    ll res = 1 ;
    while (y) {
        if (y & 1) res = res * x % mod ;
        x = x * x % mod; y >>= 1 ;
    }
    return res ;
}
int dp[1010][1010] ;
int f[2010], a[2010], res[2010], now[2010] ;
inline void mul (int x[], int y[], int k) {
    static int tmp[2010] ;
    memset (tmp, 0, sizeof tmp) ;
    for (int i = 0; i < k; i ++)
        for (int j = 0; j < k; j ++)
            tmp[i + j] = (tmp[i + j] + 1ll * x[i] * y[j] % mod) % mod ;
    for (int i = 2 * k - 2; i >= k; i --) {
        for (int j = 1; j <= k; j ++)
            tmp[i - j] = (tmp[i - j] + 1ll * tmp[i] * a[j] % mod) % mod ;
        tmp[i] = 0 ;
    }
    for (int i = 0; i < k; i ++) x[i] = tmp[i] ;
}
inline int cal (int n, int K, int p) {
    memset (dp, 0, sizeof dp) ;
    memset (f, 0, sizeof f) ;
    for (int i = 0; i <= 1000; i ++) dp[0][i] = 1 ;
    for (int i = 1; i <= 1000; i ++)
        for (int j = 1000; j >= 0; j --) {
            if (i * j > K) {
                dp[i][j] = 0; continue ;
            }
            dp[i][j] = 1ll * dp[i][j + 1] * power (p, i) % mod ;
            for (int k = 1; k <= i; k ++)
                dp[i][j] = (dp[i][j] + 1ll * power (p, k - 1) * (1 + mod - p) % mod * dp[k - 1][j + 1] % mod * dp[i - k][j] % mod) % mod ;
        }
    if (n <= K) return dp[n][0] ;
    K ++ ;
    for (int i = 1; i <= K; i ++)
        a[i] = 1ll * (1 + mod - p) * power (p, i - 1) % mod * dp[i - 1][1] % mod ;
    for (int i = 0; i < K; i ++)
        f[i] = dp[i + 1][0] ;
    memset (res, 0, sizeof res) ;
    memset (now, 0, sizeof now) ;
    n -= K; res[0] = now[1] = 1 ;
    while (n) {
        if (n & 1) mul (res, now, K) ;
        mul (now, now, K); n >>= 1 ;
    }
    for (int i = K; i <= 2 * K - 2; i ++)
        for (int j = 1; j <= K; j ++)
            f[i] = (f[i] + 1ll * f[i - j] * a[j] % mod) % mod ;
    int ans = 0 ;
    for (int i = 0; i < K; i ++)
        ans = (ans + 1ll * f[K - 1 + i] * res[i] % mod) % mod ;
    return ans ;
}
int main() {
    int n, k, x, y ;
    scanf("%d%d%d%d", &n, &k, &x, &y) ;
    int p = 1ll * x * power (y, mod - 2) % mod ;
    ll ans = cal (n, k, p) ;
    ans = (ans - cal (n, k - 1, p) + mod) % mod ;
    printf("%lld\n", ans) ;
    return 0 ;
}