Loj #2304. 「NOI2017」泳池

题目描述

久莲是个爱玩的女孩子。

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

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

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

* 必须保证安全性。即游泳场中的每一个单位海域都是安全的。

* 必须是矩形。即游泳场必须是整个网格中的一个 \(a\times b\) 的子网格。

* 必须和海滩相邻。即游泳场的下边界必须紧贴网格的下边界。

例如：当 \(N = 5\) 时，若测量的结果如下（因为 \(1001\) 太大，这儿只画出网格最下面三行的信息，其他部分都是危险的）。

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

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

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

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

输入格式

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

输出格式

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

即设答案化为最简分式后的形式为 \(\frac{a}{b}\) ，其中 \(a\) 和 \(b\) 的互质。输出整数 \(x\) 使得 \(bx \equiv a \mod 998244353\) 且 \(0 \leq x < 998244353\)。可以证明这样的整数 \(x\) 是唯一的。

数据范围与提示

\(N\leq 10^9,K\leq 1000\)

神仙题。

我们可以用求最大面积\(\leq k\)的概率减去最大面积\(\leq k+1\)的概率得到答案。

先考虑\(n\)比较小的情况。设\(f_{i,j}\)表示有\(i\)列，\(1\)到\(j\)行都安全时最大矩阵\(\leq k\)的概率。对于\(i*j>k\)的状态，\(f_{i,j}=0\)。

转移时，可以考虑第\(j+1\)行是否全部没有危险，否则枚举从左往右第一个有危险的地方。于是得到转移方程：

\[ f_{i,j}=f_{i,j+1}*q^i+\sum_{k=1}^if_{k-1,j+1}*q^{k-1}*(1-q)*f_{i-k,j} \]

由于有效状态只有\(k\)个，转移复杂度也是\(O(k)\)，所以这个\(DP\)的复杂度\(O(k^2)\)。答案就是\(f_{n,0}\)。

考虑\(n\)很大的情况，因为不可能出现连续的\(k+1\)列至少第一行是安全的情况，所以一个合法的矩阵一定是若干段矩阵，中间用第一行的危险的格子连接。

设\(F_i\)表示\(i\)列的概率，则：

\[ F_i=\sum_{j=i-k}^iF_{j-1}*(1-q)*f_{i-j,1}*q^{i-j} \]

可以发现这就是一个\(k\)阶其次线性递推。所以我们可以将转移写成

\[ F_i=\sum_{j=0}^{k}F_{i-k-1+j}*a_j \]

其中\(a_i=(1-q)*f_{k-i,1}*q^{k-i}\)。

我们发现这就是个\(k\)阶齐次线性递推。

设\(A(x)=x^{k+1}-\sum_{i=0}^ka_ix^i，B(x)=x^n\)，然后我们要求出\(C(x)=B(x)\%A(x)\)，答案就是\(\sum_{i=0}^k c_i*F_i\)，其中\(F_i=f_{i,0}\)。

但是\(C(x)\)最高次可以达到\(10^9\)，不可能直接取模。

由：

\[ f(x)\equiv g(x)\pmod{A(x)}\\ f(x)^2\equiv g(x)^2\pmod{A(x)}\\ \]

于是就可以用类似于快速幂的方法来求\(x^n\%A(x)\)。

代码：

#include
    
     #define ll long long#define K 1005using namespace std;inline int Get() {int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}const ll mod=998244353;ll ksm(ll t,ll x) {    ll ans=1;    for(;x;x>>=1,t=t*t%mod)        if(x&1) ans=ans*t%mod;    return ans;}const int maxx=1001;int n,k;ll q,p;ll f[K][K];ll F[K];ll G[K],a[K];ll r[K],b[K];void mul(ll *x,ll *y,ll *ans,int len) {    static ll tem[K<<1];    for(int i=0;i<2*len;i++) tem[i]=0;    for(int i=0;i
     
      =len;i--) {        if(tem[i]) {            for(int j=0;j
      
       =0;j--) {            if(i*j>lim) continue ;            f[i][j]=f[i][j+1]*ksm(q,i)%mod;            for(int k=1;k<=i;k++) {                (f[i][j]+=f[k-1][j+1]*ksm(q,k-1)%mod*p%mod*f[i-k][j])%=mod;            }        }    }    for(int i=0;i<=lim;i++) a[lim-i]=f[i][1]*ksm(q,i)%mod*p%mod;    r[0]=b[1]=1;    int now=n;    for(;now;now>>=1,mul(b,b,b,lim+1)) {        if(now&1) mul(r,b,r,lim+1);    }    ll ans=0;    for(int i=0;i<=lim;i++) (ans+=r[i]*f[i][0])%=mod;    return ans;}int main() {    n=Get(),k=Get();    int x=Get(),y=Get();    q=x*ksm(y,mod-2)%mod;    p=(1-q+mod)%mod;    cout<<(solve(k)-solve(k-1)+mod)%mod;    return 0;}