[NOI2017]泳池

Description
　　有一个$1001×n$的的网格，最下面一行连着海滩，每个格子有$q$的概率是安全的，$1−q$的概率是危险的。
　　定义一个矩形是合法的当且仅当：

　　- 这个矩形中每个格子都是安全的
　　- 必须紧贴网格的下边界

　　问你最大的合法子矩形大小为$k$的概率是多少。

　　$n≤10^9,k≤1000$

Solution
　　“快喊救护车来，这里有一棵白菜被题目虐死了！！！”
　　这题目做了一天。。（我好菜菜啊

　　首先拆开询问，等价于是用小于等于$k$的答案减去小于等于$k - 1$的答案。
　　我们来一个 根本想不出的 $DP$，令$f[i][j]$表示考虑一共有$j$列，现在从最靠近海滩的那一行（称为第$1$行）到第$i$行的情况不考虑，默认为一定安全。然后我们在第$i + 1$行放至少一个危险格子，保证当前最大合法子矩形合法的概率。同时注意保证$i * j \le k$，即就算下面的$i$行全是安全的仍然满足条件。
　　那么有：
　　

$$f[i][j] = \sum_{t = 0}^{j - 1} \sum_{u > i} \sum_{v \ge i}(f[u][t] * q^{t*(u-i)})(1 - q)(f[v][j - t - 1] * q^{(j - t - 1)(v - i)})$$
　　意思是在第$i + 1$行枚举一个危险格子，它的坐标是$(i + 1, t + 1)$，它的左边和右边我们分别枚举，并中间有多少空的就应该乘$q$的那么多次方。
　　（为什么一个是等于一个是大于等于呢？因为若存在几个连续高度为一的危险格子，那么就会把方案算重，这么写就保证只在连续的那一段最后一个高度为一的地方计入答案）
　　前缀和优化，我们记$F[i][j] = \sum_{u \ge i}f[u][j]$，枚举$i$和$j$的时间为$k$，所以我们复杂度可以做到$O(k^2)$

　　定义$h[i]$为有$i$列的答案是多少，于是有
　　

$$h[i] = f[1][i] * q^{i} + \sum_{j = 1}^{i} f[1][j - 1] * q^{j - 1} * (1 - q) * h[i - j] ---i \le k$$
　　

$$h[i] = \sum_{j = 1}^{k + 1} f[1][j - 1] * q^{j - 1} * (1 - q) * h[i - j] ---i > k$$
　　那么求解出 $h[n]$即为答案。

---

　　然而并没有介素，还过不了呐筒子们啊！矩阵乘法优化直接炸飞起～
　　
　　来来来，先看一下这个戳我戳我:)，学习一下常系数齐次线性递推优化矩阵快速幂，找了一上午资料最后看见这个感觉整棵菜都好起来了。然后这个会写那个链接所补充的模板题就好。
　　我们可以先推出$i \le k$的$h[i]$们，之后那个式子，$h[i - j]$的前面都是常数对吧，直接就是个其次线性方程了，就用那个方法做就行。暴力多项式取模是$O(k^2)$的，总复杂度$O(k^2logn)$可以通过此题了。
　　
　　好像还有$FFT$优化多项式取模？！留坑。

Source

//2018-4-28
//miaowey
//
#include <bits/stdc++.h>
using namespace std;

#define Set(a, v) memset(a, v, sizeof(a))
#define For(i, a, b) for(int i = (a); i <= (int)(b); ++i)
#define Forr(i, a, b) for(int i = (a); i >= (int)(b); --i)

#define N (1000 + 5)
const int P = 998244353;

inline int Mod(int a){
    if(a < 0) return a + P;
    return a >= P? a - P: a;
}
inline int Mul(int a, int b){
    return (long long)a * b % P;
}

inline int Pow(int a, int anum){
    int ret = 1;
    while(anum){
        if(anum & 1) ret = Mul(ret, a);
        a = Mul(a, a); anum >>= 1;
    }
    return ret;
}

int n, x, y, p, q, powq[N], f[N][N], g[N], h[N], b[N], c[N], tmp[N << 1];

void ExpMul(int *A, int *B, int k){
    For(i, 0, k) For(j, 0, k) tmp[i + j] = Mod(tmp[i + j] + Mul(A[i], B[j]));

    Forr(i, k << 1, k){
        For(j, 0, k - 1) tmp[i + j - k] = Mod(tmp[i + j - k] + Mul(tmp[i], g[k - j]));
        tmp[i] = 0;
    }

    For(i, 0, k - 1) A[i] = tmp[i], tmp[i] = 0;
}

void ExpPow(int *A, int anum, int *ret, int k){
    while(anum){
        if(anum & 1) ExpMul(ret, A, k);
        ExpMul(A, A, k); anum >>= 1;
    }
}

int Solve(int k){
    if(!k) return Pow(p, n);

    Set(f, 0); f[k + 1][0] = 1;

    Forr(i, k, 1){
        int m = min(n, k / i);

        f[i][0] = 1;
        For(j, 0, m) g[j] = Mul(Mul(f[i + 1][j], powq[j]), p);

        For(j, 1, m){
            For(t, 0, j - 1) f[i][j] = Mod(f[i][j] + Mul(g[t], f[i][j - t - 1]));
            f[i][j] = Mod(f[i][j] + Mul(f[i + 1][j], powq[j]));
        }
    }

    ++k;
    For(i, 1, k) g[i] = Mul(Mul(f[1][i - 1], powq[i - 1]), p);

    h[0] = 1;
    For(i, 1, k - 1){
        h[i] = Mul(powq[i], f[1][i]);
        For(j, 1, i) h[i] = Mod(h[i] + Mul(h[i - j], g[j]));
    }

    if(n < k) return h[n];

    int ret = 0;

    Set(b, 0); Set(c, 0); b[0] = c[1] = 1;
    ExpPow(c, n, b, k);

    For(i, 0, k - 1) ret = Mod(ret + Mul(b[i], h[i]));

    return ret;
}

int main(){
#ifndef ONLINE_JUDGE
    freopen("pool.in", "r", stdin);
    freopen("pool.out", "w", stdout);
#endif

    int k;

    scanf("%d%d%d%d", &n, &k, &x, &y);  
    q = Mul(x, Pow(y, P - 2)); p = Mod(1 - q);

    powq[0] = 1;
    For(i, 1, 1000) powq[i] = Mul(powq[i - 1], q);

    printf("%d\n", Mod(Solve(k) - Solve(k - 1)));

    return 0;
}