Loj #2304. 「NOI2017」泳池

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)行是否全部没有危险，否则枚举从左往右第一个有危险的地方。于是得到转移方程：

[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)^2equiv g(x)^2pmod{A(x)}\ ]

于是就可以用类似于快速幂的方法来求(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;
}