【NOI2017】 泳池(DP,多项式,常系数齐次线性递推)

89 篇文章 0 订阅
16 篇文章 0 订阅

题面

题目背景

久莲 久条可莲 是个爱玩的女孩子。

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

题目描述

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

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

  • 必须保证安全性。即游泳场中的每一个单位海域都是安全的。
  • 必须是矩形。即游泳场必须是整个网格中的一个 a × b a\times b a×b 的子网格。
  • 必须和海滩相邻。即游泳场的下边界必须紧贴网格的下边界。

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

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

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

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

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

输入格式

输入一行四个正整数 N , K , x , y N,K,x,y N,K,x,y,其中 1 ≤ x < y < 998   244   353 1 \leq x < y < 998\,244\,353 1x<y<998244353 q q q 的取值为 x y \dfrac{x}{y} yx

输出格式

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

即设答案化为最简分式后的形式为 a b \dfrac{a}{b} ba ,其中 a a a b b b 的互质。输出整数 x x x 使得 b x ≡ a m o d    998   244   353 bx \equiv a \mod 998\,244\,353 bxamod998244353 0 ≤ x < 998   244   353 0 \leq x < 998\,244\,353 0x<998244353。可以证明这样的整数 x x x 是唯一的。

样例 #1

样例输入 #1

10 5 1 2

样例输出 #1

342025319

提示

测试点编号 N N N K K K
1,2 = 1 =1 =1 ≤ 1000 \leq 1000 1000
3 ≤ 10 \leq 10 10 ≤ 8 \leq 8 8
4 ≤ 10 \leq 10 10 ≤ 9 \leq 9 9
5 ≤ 10 \leq 10 10 ≤ 10 \leq 10 10
6 ≤ 1000 \leq 1000 1000 ≤ 7 \leq 7 7
7 ≤ 1000 \leq 1000 1000 ≤ 8 \leq 8 8
8 ≤ 1000 \leq 1000 1000 ≤ 9 \leq 9 9
9,10,11 ≤ 1000 \leq 1000 1000 ≤ 100 \leq 100 100
12,13,14 ≤ 1000 \leq 1000 1000 ≤ 1000 \leq 1000 1000
15,16 ≤ 1 0 9 \leq 10^9 109 ≤ 10 \leq 10 10
17,18 ≤ 1 0 9 \leq 10^9 109 ≤ 100 \leq 100 100
19,20 ≤ 1 0 9 \leq 10^9 109 ≤ 1000 \leq 1000 1000

题解

首先转化一下,我们求 ≤ K \leq K K 的概率减去 ≤ K − 1 \leq K-1 K1 的概率。


s t i st_i sti 表示底部长为 i i i 、最大游泳场小于等于 K K K 、且右下角不安全的概率。那么我们要求的就是 s t n + 1 1 − q \frac{st_{n+1}}{1-q} 1qstn+1 ,相当于在原长为 n n n 的泳池右边放个哨兵。

我们考虑怎么递推求出 s t st st

枚举上一个卡到底部的位置,该位置一定距离右边界不超过 K K K ,所以
s t i = ∑ j = 1 min ⁡ ( K , i ) s t i − j ( 1 − q ) P j − 1 令   a K = 1 , a i = − P K − i − 1 ( 1 − q ) ⇒ ∑ j = 0 K a j s t i + j = 0 st_i=\sum_{j=1}^{\min(K,i)} st_{i-j}(1-q)P_{j-1} \\ 令~a_K=1,a_{i}=-P_{K-i-1}(1-q) \\\Rightarrow \sum_{j=0}^{K} a_jst_{i+j}=0 sti=j=1min(K,i)stij(1q)Pj1 aK=1,ai=PKi1(1q)j=0Kajsti+j=0

其中这个 P i P_i Pi 表示底长为 i i i 的泳池、没有触底的危险区、且满足 K K K 的限制的概率。


我们可以通过讨论最低危险区的高度来设 DP 进行转移计算 P P P,令 d p i , j dp_{i,j} dpi,j 表示底长为 i i i ,最低危险区高度为 j + 1 j+1 j+1 的方案数。可以用笛卡尔树的思路转移:
∀ i j > K , d p i , j = 0 d p i , j = ( 1 − q ) q j ∑ t = 1 i ( ∑ d = j + 1 ∞ d p [ t − 1 ] [ d ] ) ( ∑ d = j ∞ d p [ i − t ] [ d ] ) \forall ij>K,dp_{i,j}=0\\ dp_{i,j} =(1-q)q^j\sum_{t=1}^i\left(\sum_{d=j+1}^{\infty}dp[t-1][d]\right)\left(\sum_{d=j}^{\infty}dp[i-t][d]\right) ij>K,dpi,j=0dpi,j=(1q)qjt=1id=j+1dp[t1][d]d=jdp[it][d]

后缀和优化得
d p i , j = ( 1 − q ) q j ∑ t = 1 i s u f t − 1 , j + 1 ⋅ s u f i − t , j dp_{i,j} =(1-q)q^j\sum_{t=1}^isuf_{t-1,j+1}\cdot suf_{i-t,j} dpi,j=(1q)qjt=1isuft1,j+1sufit,j

由于有效的 DP 状态只有 O ( K ln ⁡ K ) O(K\ln K) O(KlnK) 个,所以可以暴力转移。

根据定义可知 P i = s u f i , 1 P_i=suf_{i,1} Pi=sufi,1


然后我们求出 s t st st 的前 0 ∼ K 0\sim K 0K 项( s t 0 = 1 st_0=1 st0=1),就可以进行常系数齐次线性递推了!

我知道很多人死活学不会,我也是,所以我直接把具体要干什么写下来了。


参考

设状态转移矩阵为 A A A ,初始矩阵为 S = [ s t 0 , s t 1 , . . . , s t K ] S=[st_0,st_1,...,st_{K}] S=[st0,st1,...,stK],我们需要的是
( S ⋅ A n + 1 ) 1 , 1 (S\cdot A^{n+1})_{1,1} (SAn+1)1,1

而事实上存在一个神奇的多项式 C ( A ) C(A) C(A) 满足
A n + 1 = C ( A ) = ∑ i = 0 K c i A i A^{n+1}=C(A)=\sum_{i=0}^K c_iA^i An+1=C(A)=i=0KciAi

我们只要求出了 C ( A ) C(A) C(A) ,就可以直接计算目标答案:
( S ⋅ A n + 1 ) 1 , 1 = ∑ i = 0 K c i ⋅ ( S ⋅ A i ) 1 , 1 = ∑ i = 0 K c i s t i (S\cdot A^{n+1})_{1,1}=\sum_{i=0}^K c_i\cdot (S\cdot A^i)_{1,1}=\sum_{i=0}^Kc_ist_i (SAn+1)1,1=i=0Kci(SAi)1,1=i=0Kcisti

怎么求出 C C C 呢?我们考虑构造一个多项式 F F F ,满足
F ( A ) = ∑ i = 0 K + 1 f i A i = [ 0 ] F(A)=\sum_{i=0}^{K+1} f_iA^i=[0] F(A)=i=0K+1fiAi=[0]

那么 x n + 1 m o d    F = C x^{n+1}\mod F=C xn+1modF=C ,可以用快速幂套多项式取模求出 C C C ,时间复杂度 O ( K log ⁡ n log ⁡ K ) O(K\log n\log K) O(KlognlogK)

F F F 怎么构造呢?常系数齐次线性递推告诉我们一个真理:
f i = a i f_{i}=a_i fi=ai

就是上面满足 ∑ j = 0 K a j s t i + j = 0 \sum_{j=0}^{K} a_jst_{i+j}=0 j=0Kajsti+j=0 的那个 a i a_i ai

CODE

不会厂除法 😢

⚠长代码警告⚠

#include<map>
#include<set>
#include<cmath>
#include<ctime>
#include<queue>
#include<stack>
#include<random>
#include<bitset>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
// #pragma GCC optimize(2)
// #pragma GCC optimize("Ofast")
using namespace std;
#define MAXN 2005
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
int xchar() {
	static const int maxn = 1000000;
	static char b[maxn];
	static int pos = 0,len = 0;
	if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
	if(pos == len) return -1;
	return b[pos ++];
}
// #define getchar() xchar()
LL read() {
	LL f = 1,x = 0;int s = getchar();
	while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
	while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
	return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar((x%10)^48);}
void putnum(LL x) {
	if(!x) {putchar('0');return ;}
	if(x<0) putchar('-'),x = -x;
	return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}

const int MOD = 998244353;
int n,m,s,o,k;
int xm[MAXN<<2],om,rev[MAXN<<2];// tool
int qkpow(int a,int b) { // It's necessary
	int res = 1; while(b > 0) {
		if(b & 1) res = res *1ll* a % MOD;
		a = a *1ll* a % MOD; b >>= 1;
	}return res;
}
const int RM = 3;     // Prime Root
// Let's Start
inline void NTT(int *s,int n) {
	for(int i = 1;i < n;i ++) {
		rev[i] = ((rev[i>>1] >> 1) | ((i&1) ? (n>>1):0));
		if(rev[i] < i) swap(s[rev[i]],s[i]);
	}
	om = qkpow(RM,(MOD-1)/n); xm[0] = 1;
	for(int k = 1;k < n;k ++) xm[k] = xm[k-1]*1ll*om % MOD;
	for(int k = 2,t = (n>>1);k <= n;k <<= 1,t >>= 1)
		for(int j = 0;j < n;j += k)
			for(int i = j,l = 0;i < j+(k>>1);i ++,l += t) {
				int A = s[i],B = s[i+(k>>1)];
				s[i] = (A + xm[l]*1ll*B%MOD) % MOD, s[i+(k>>1)] = (A +MOD- xm[l]*1ll*B%MOD) % MOD;
			}
	return ;
}
inline void INTT(int *s,int n) {
	for(int i = 1;i < n;i ++) {
		rev[i] = ((rev[i>>1] >> 1) | ((i&1) ? (n>>1):0));
		if(rev[i] < i) swap(s[rev[i]],s[i]);
	}
	om = qkpow(qkpow(RM,(MOD-1)/n),MOD-2); xm[0] = 1;
	for(int k = 1;k < n;k ++) xm[k] = xm[k-1]*1ll*om % MOD;
	for(int k = 2,t = (n>>1);k <= n;k <<= 1,t >>= 1)
		for(int j = 0;j < n;j += k)
			for(int i = j,l = 0;i < j+(k>>1);i ++,l += t) {
				int A = s[i],B = s[i+(k>>1)];
				s[i] = (A + xm[l]*1ll*B%MOD) % MOD, s[i+(k>>1)] = (A +MOD- xm[l]*1ll*B%MOD) % MOD;
			}
	int invn = qkpow(n,MOD-2);
	for(int i = 0;i <= n;i ++) s[i] = s[i] *1ll* invn % MOD;
	return ;
}
//------------------------NTT End
inline void polymul(int *A,int *B,int n) { // A *= B mod x^(n+1), and clear B
	int le = 1;
	while(le <= n*2) le <<= 1;
	for(int i = n+1;i <= le;i ++) A[i] = B[i] = 0;
	NTT(A,le);NTT(B,le);
	for(int i = 0;i <= le;i ++) A[i] = A[i]*1ll*B[i] % MOD,B[i]=0;
	INTT(A,le);for(int i = n+1;i <= le;i ++) A[i]=0;
}
int a_[MAXN<<2],b_[MAXN<<2];
inline void polyinv(int *s,int n) { // mod x^n
	for(int i = 0;i <= n;i ++) a_[i] = s[i],s[i]=0;
	s[0] = qkpow(a_[0],MOD-2);
	int le = 1;
	while(le < n) {
		le <<= 1;
		for(int i = 0;i < le;i ++) b_[i]=a_[i];
		NTT(b_,le<<1),NTT(s,le<<1);
		for(int i = 0;i <= (le<<1);i ++)
			s[i] = (2ll+MOD-b_[i]*1ll*s[i] % MOD) % MOD *1ll* s[i] % MOD,b_[i]=0;
		INTT(s,le<<1);
		for(int i = le;i <= (le<<1);i ++) s[i] = 0;
	}
	for(int i = n;i <= le;i ++) s[i] = 0;
	for(int i = 0;i <= n;i ++) a_[i] = 0; // tool cleared
	return ;
}
inline void polyrev(int *s,int n) { // 0~n
	for(int i = 0;(i<<1) <= n;i ++) swap(s[i],s[n-i]); return ;
}
int s1_[MAXN<<2];
inline void polydiv(int *f,int *g,int *q,int *r,int n,int m) { // f=q*g+r , q,r is wanted
	for(int i = 0;i <= n;i ++) q[i] = f[i];polyrev(q,n);
	for(int i = n-m+1;i <= n;i ++) q[i] = 0;
	for(int i = 0;i <= m;i ++) s1_[i] = g[i];polyrev(s1_,m);
	for(int i = n-m+1;i <= m;i ++) s1_[i] = 0;
	polyinv(s1_,n-m+1);
	polymul(q,s1_,n-m); polyrev(q,n-m);
	for(int i = 0;i <= n;i ++) s1_[i] = g[i],r[i] = q[i];
	polymul(r,s1_,n);
	for(int i = 0;i <= n;i ++) r[i] = (f[i]+MOD-r[i]) % MOD;
	for(int i = m;i <= n;i ++) r[i] = 0;
	return ;
}
//-----------------------END----------------
int q;
int dp[1005][1005],sdp[1005][1005],tl[MAXN<<2],tl2[MAXN<<2];
int A[MAXN<<2],B[MAXN<<2],F[MAXN<<2],f[MAXN];
int solve(int K) {
	if(K < 0) return 0;
	if(K == 0) return qkpow((1+MOD-q)%MOD,n);
	memset(dp,0,sizeof(dp));
	memset(sdp,0,sizeof(sdp));
	for(int i = 0;i <= K+1;i ++) sdp[0][i] = 1;
	for(int i = 1;i <= K;i ++) {
		for(int j = K/i;j > 0;j --) {
			int p = qkpow(q,j) *1ll* (MOD+1-q) % MOD;
			for(int k = 1;k <= i;k ++) {
				(dp[i][j] += sdp[k-1][j+1] *1ll* sdp[i-k][j] % MOD) %= MOD;
			}
			dp[i][j] = dp[i][j] *1ll* p % MOD;
			sdp[i][j] = (sdp[i][j+1] + dp[i][j]) % MOD;
		}
	}
	f[0] = 1; K ++;
	for(int i = 1;i < K;i ++) {
		f[i] = 0;
		for(int j = 1;j <= i;j ++) (f[i] += f[i-j] *1ll* sdp[j-1][1] % MOD * (1+MOD-q) % MOD) %= MOD;
	}
	for(int i = 1;i <= K;i ++) {
		F[K-i] = sdp[i-1][1] *1ll* (q+MOD-1) % MOD;
	} F[K] = 1;
	int nn = n+1;
	int le = 1; while(le <= K+K)le<<=1;
	for(int i = 0;i < le;i ++) A[i] = B[i] = 0;
	A[0] = 1; B[1] = 1;
	while(nn > 0) {
		if(nn & 1) {
			NTT(A,le); NTT(B,le);
			for(int i = 0;i < le;i ++) A[i] = A[i] *1ll* B[i] % MOD;
			INTT(A,le); INTT(B,le);
			polydiv(A,F,tl2,tl,K<<1,K);
			for(int i = 0;i < K;i ++) A[i] = tl[i];
			for(int i = K;i < le;i ++) A[i] = 0;
		}
		NTT(B,le);
		for(int i = 0;i < le;i ++) B[i] = B[i] *1ll* B[i] % MOD;
		INTT(B,le);
		polydiv(B,F,tl2,tl,K<<1,K);
		for(int i = 0;i < K;i ++) B[i] = tl[i];
		for(int i = K;i < le;i ++) B[i] = 0;
		nn >>= 1;
	}
	int ans = 0;
	for(int i = 0;i < K;i ++) {
		(ans += A[i] *1ll* f[i] % MOD) %= MOD;
	}
	ans = ans *1ll* qkpow((1+MOD-q)%MOD,MOD-2) % MOD;
	return ans;
}
int main() {
	n = read();m = read();
	q = read(); q = q *1ll* qkpow(read(),MOD-2) % MOD;
	int ans = (solve(m) +MOD- solve(m-1)) % MOD;
	AIput(ans,'\n');
	return 0;
}

哦,终于学会了,厂除法真简单。

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值