Loj #2304. 「NOI2017」泳池
题目描述
久莲是个爱玩的女孩子。
暑假终于到了，久莲决定请她的朋友们来游泳，她打算先在她家的私人海滩外圈一块长方形的海域作为游泳场。然而大海里有着各种各样的危险，有些地方水太深，有些地方有带毒的水母出没。她想让圈出来的这一块海域都是安全的。
经过初步分析，这块海域可视为一个底边长为 (N) 米，高为 (1001) 米的长方形网格。其中网格的底边对应着她家的私人海滩，每一个 (1: extrm{m} imes1: extrm{m}) 的小正方形都代表着一个单位海域。她拜托了她爸爸明天去测量每一个小正方形是否安全。在得知了信息之后，她要做的就是圈出她想要的游泳场啦。
她心目中理想的游泳场满足如下三个条件：
* 必须保证安全性。即游泳场中的每一个单位海域都是安全的。
* 必须是矩形。即游泳场必须是整个网格中的一个 (a imes b) 的子网格。
* 必须和海滩相邻。即游泳场的下边界必须紧贴网格的下边界。
例如：当 (N = 5) 时，若测量的结果如下（因为 (1001) 太大，这儿只画出网格最下面三行的信息，其他部分都是危险的）。
那么她可以选取最下面一行的 (1 imes4) 的子海域，也可以选择第三列的 (3 imes1) 的子海域。注意她不能选取最上面一行的 (1 imes5) 的子海域，因为它没有与海滩相邻。
为了让朋友们玩的开心，她想让游泳场的面积尽可能的大。因此她会选取最下面那一行的 (1 imes4) 的子海域作为最终方案。
虽然她要明天才能知道每一个单位海域是否安全，但是她现在就想行动起来估计一下她的游泳场面积有多大。经过简单的估计，她假设每一个单位海域都有独立的 (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) 是唯一的。
数据范围与提示
(Nleq 10^9,Kleq 1000)
神仙题。
我们可以用求最大面积(leq k)的概率减去最大面积(leq k+1)的概率得到答案。
先考虑(n)比较小的情况。设(f_{i,j})表示有(i)列，(1)到(j)行都安全时最大矩阵(leq k)的概率。对于(i*j>k)的状态，(f_{i,j}=0)。
转移时，可以考虑第(j+1)行是否全部没有危险，否则枚举从左往右第一个有危险的地方。于是得到转移方程：
由于有效状态只有(k)个，转移复杂度也是(O(k))，所以这个(DP)的复杂度(O(k^2))。答案就是(f_{n,0})。
考虑(n)很大的情况，因为不可能出现连续的(k+1)列至少第一行是安全的情况，所以一个合法的矩阵一定是若干段矩阵，中间用第一行的危险的格子连接。
设(F_i)表示(i)列的概率，则：
可以发现这就是一个(k)阶其次线性递推。所以我们可以将转移写成
其中(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)，不可能直接取模。
由：
于是就可以用类似于快速幂的方法来求(x^n\%A(x))。
代码：
#include<bits/stdc++.h>
#define ll long long
#define K 1005
using 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++)
for(int j=0;j<len;j++)
(tem[i+j]+=x[i]*y[j])%=mod;
for(int i=2*len-1;i>=len;i--) {
if(tem[i]) {
for(int j=0;j<len;j++) (tem[i-len+j]+=a[j]*tem[i])%=mod;
}
}
for(int i=0;i<len;i++) ans[i]=tem[i];
}
ll solve(int lim) {
memset(f,0,sizeof(f));
memset(F,0,sizeof(F));
memset(r,0,sizeof(r));
memset(b,0,sizeof(b));
for(int i=0;i<=maxx;i++) f[0][i]=1;
for(int i=1;i<=1000;i++) {
for(int j=1000;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;
}