[BZOJ4944][NOI2017]泳池-概率与期望-常系数线性递推-多项式取模

泳池

久莲是个爱玩的女孩子。

暑假终于到了，久莲决定请她的朋友们来游泳，她打算先在她家的私人海滩外圈一块长方形的海域作为游泳场。然而大海里有着各种各样的危险，有些地方水太深，有些地方有带毒的水母出没。她想让圈出来的这一块海域都是安全的。

经过初步的分析，她把这块海域抽象成了一个底边长为 N 米，高为 1001 米的长方形网格。其中网格的底边对应着她家的私人海滩，每一个 1米×1米 的小正方形都代表着一个单位海域。她拜托了她爸爸明天去测量每一个小正方形是否安全。在得知了信息之后，她要做的就是圈出她想要的游泳场啦。

她心目中理想的游泳场满足如下三个条件：

必须保证安全性。即游泳场中的每一个单位海域都是安全的。
必须是矩形。即游泳场必须是整个网格中的一个 a×b 的子网格。
必须和海滩相邻。即游泳场的下边界必须紧贴网格的下边界。
例如：当 N=5 时，若测量的结果如下（因为 1001 太大，这儿只画出网格最下面三行的信息，其他部分都是危险的）。

游泳池

那么她可以选取最下面一行的 1×4 的子海域，也可以选择第三列的 3×1 的子海域。注意她不能选取最上面一行的 1×5 的子海域，因为它没有与海滩相邻。

为了让朋友们玩的开心，她想让游泳场的面积尽可能的大。因此她会选取最下面那一行的 1×4 的子海域作为最终方案。

虽然她要明天才能知道每一个单位海域是否安全，但是她现在就想行动起来估计一下她的游泳场面积有多大。经过简单的估计，她假设每一个单位海域都有独立的 q 的概率是安全的，1−q 的概率是不安全的。她想要知道她能选择的最大的游泳场的面积恰好为 K 的概率是多少。

然而久莲对数学并不感兴趣，因此她想让你来帮她计算一下这个数值。

输入格式

从标准输入读入数据。

输入一行四个正整数 N,K,x,y，其中 1≤x< y <998244353。q 的取值为 $\frac{x}{y}$。

输出格式

输出到标准输出。

输出一行一个整数表示答案在模 998244353 意义下的取值。

即设答案化为最简分式后的形式为 ab ，其中 a 和 b 的互质。输出整数 x 使得 $bx≡a\ mod\ 998244353$且$0≤x<998244353$。可以证明这样的整数 xx 是唯一的。

样例一

input

10 5 1 2

output

342025319

样例二

见“样例数据下载”

样例三

见“样例数据下载”

提示

$x^{p−1}≡1\ mod\ p$，其中 $p$ 为质数，$x \in [1,p)$。

限制与约定

$N \leq 10^9 k \leq 1000$
时间限制：3s
空间限制：512Mb

---

好难啊……
研究了好久才算是明白了一点……
菜的不行QAQ

---

思路:
这个最大矩形不好做，考虑转换为最大子矩形$\leq k$的概率减去最大子矩形$\leq k-1$的概率。
于是问题变成了，求最大子矩形面积$\leq k$的概率。

由于需要与海滩相邻，因此可以考虑每一列最靠下的一个不可选位置，这个位置下方均是可选的。
设$f[i][j]$表示，默认下面$i$行是合法的，第$i+1$行存在不合法位置，宽度为$j$的矩形面积合法的概率。

那么有

$$f[i][j]=\sum\limits_{t=0}^{j-1}\sum\limits_{k \geq i}\sum\limits_{l>i} f[k][t]*q^{t*(k-i)}*f[l][j-1-t]*q^{(j-t-1)*(l-i)}*(1-q)$$
意思是，左边和右边不论大小拼起来，中间设置一个不合法的限制最大高度。
由于$f[i][j]$默认下面$i$行合法，因此转移时需要乘上默认的那几行合法的概率。

接着还需要一个统计答案:
设$ans[i]$表示只考虑前$i$列，最大子矩形合法的概率。
于是有

$$ans[i]=f[1][i]*q^{i}+\sum\limits_{j=1}^{min(i,k+1)}ans[i-j]*f[1][j-1]*q^{j-1}*(1-q)$$

此处按连续的一段合法底边转移，前半部分为一整段合法的情况，后半部分为枚举最靠后一段合法底边的情况。

然而这部分是$O(nk)$的，因此需要优化。

首先显然这个转移可以写成矩阵。
复杂度降至$O(k^3logn)$，$k=1000$ 时依旧过不去。

把转移写成$ans[i]=\sum\limits_{j=1}^{k+1}a_jans[i-j]$的形式。
可以发现这个转移是个常系数线性递推式，考虑化简。

考虑矩阵乘法最初的使用方法，现在有一个向量$T$，在乘以$A$转移矩阵$n$次后，得到$[TA^n]_0$即为答案。
$T$很好求，瓶颈在于$A^n$的求法。

构造$A$的特征多项式。
一个矩阵$A$的特征多项式为$g(\lambda)=det(\lambda E-A)$，其中$E$为单位矩阵。
可以发现当$\lambda==A$时$g(A)=0$。
同时可以发现$g(\lambda)=\lambda^{|A|}-\sum\limits_{i=1}^{|A|} a_i \lambda^{|A|-i}$

那么将$A$看成变量，则令多项式$B$满足$A^n \equiv B\ mod\ g(A)$。
由于已知$g(A)=0$，因此只需要关注$B$，即$A^n$等价于$B$。
于是运用多项式取模等操作得到$B$的各项系数。

于是就有
$A^n=\sum\limits_{i=0}^{k} b_i*A^i$
由于咱们只关心$[TA^n]_0$，考虑到每个$[TA^i]_0$表示$ans[i+1]$的值，那么答案便是

$$ans[n]=\sum\limits_{i=0}^{k} b_i*ans[i]$$

直接暴力多项式快速幂+取模即可!

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

typedef long long ll;
const int K=1009;
const ll md=998244353;

ll n,k,p,q,x,y;
ll f[K][K],g[K],powq[K],ans[K];
ll a[K],b[K];

inline ll qpow(ll a,ll b)
{
    ll ret=1;
    while(b)
    {
        if(b&1)
            ret=ret*a%md;
        a=a*a%md;
        b>>=1;
    }
    return ret;
}

inline void mul(ll *a,ll *b,int k)
{
    static ll c[K*2];
    memset(c,0,sizeof(c));
    for(int i=0;i<=k;i++)
        if(a[i])
            for(int j=0;j<=k;j++)
                (c[i+j]+=a[i]*b[j]%md)%=md;
    for(int i=k*2;i>=k+1;i--)
        if(c[i])
            for(int j=1;j<=k+1;j++)
                (c[i-j]+=c[i]*g[j]%md)%=md;
    memcpy(a,c,sizeof(a[0])*(k+1));
}

inline ll solve(int k)
{
    if(!k)return qpow(1-q,n);
    memset(f,0,sizeof(f));
    f[k+1][0]=1;

    for(int j=k;j>=1;j--)
    {
        f[j][0]=1;
        int mi=min((int)n,k/j);
        int ml=min(mi-1,k/(j+1));
        for(int l=0;l<=ml;l++)
            g[l]=f[j+1][l]*powq[l]%md*p%md;
        for(int i=1;i<=mi;i++)
        {
            int ml=min(i-1,k/(j+1));
            for(int l=0;l<=ml;l++)
                (f[j][i]+=g[l]*f[j][i-1-l]%md)%=md;
            (f[j][i]+=f[j+1][i]*powq[i])%=md;
        }
    }   

    ans[0]=1;
    for(int i=1;i<=k+1;i++)
        g[i]=f[1][i-1]*powq[i-1]%md*p%md;
    for(int i=1;i<=k;i++)
    {
        ans[i]=f[1][i]*powq[i]%md;
        for(int j=1;j<=i;j++)
            (ans[i]+=ans[i-j]*g[j])%=md;
    }

    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    a[1]=1;b[0]=1;

    int pcnt=n;
    while(pcnt)
    {
        if(pcnt&1)
            mul(b,a,k);
        mul(a,a,k);
        pcnt>>=1;
    }

    ll ret=0;
    for(int i=0;i<=k;i++)
        (ret+=b[i]*ans[i]%md)%=md;

    return ret;
}

int main()
{
    scanf("%lld%lld%lld%lld",&n,&k,&x,&y);
    x%=md;y%=md;
    q=x*qpow(y,md-2)%md;
    p=(y-x)*qpow(y,md-2)%md;
    powq[0]=1;
    for(int i=1;i<=k;i++)
        powq[i]=powq[i-1]*q%md;
    printf("%lld\n",(solve(k)-solve(k-1)+md)%md);
    return 0;
}