UOJ #316: 泳池 题解

我觉得这个dp还挺难的
首先计算面积恰好为k的比较难办，不妨容斥，计算面积不大于k的，减去面积不大于k-1的就是恰好为k的
对于一个段，我们关心的显然是最下面的危险区域，因为上面的对最终的矩形划分是没有影响的
考虑dp[i][j]表示考虑宽度为j的区间，前i-1行都没有危险区域，第i行有危险区域，能圈出的最大面积不大于k的概率，这里的概率暂时不考虑1~i-1行都不危险的那部分的q的若干次方
转移我觉得挺棒的，我们考虑枚举当前行的最后一个危险的格子，那么这个位置右边显然是没有危险格子了，转移的时候要从大于i的行转移，左边可以没有也可以有危险的格子，转移的时候要从大于等于i的行转移，把转移方程写出来就是

$$dp[i][j]=\sum_{k\geq i}\sum_{l\ge i}\sum_{t=0}^{j-1}dp[k][t]*dp[l][j-1-t]*(1-q)\\ *q^{t(k-i)}*q^{(j-1-t)*(l-i)}$$
要注意的是dp[i][j]没有处理i下面的行全部不危险的概率，所以转移的时候要乘上最后的q的若干次方
这个看起来是 $O(k^5)$的，但是考虑到 $(i-1)*j\leq k$时才需要转移否则概率一定为0，所以实际上是 $O(k^4)$的，但这个复杂度还是太坏了
我们发现这个dp枚举k和l的部分似乎都可以前缀和优化，我们考虑重新定义dp状态
dp[i][j]表示宽度为j的区域，第一个出现危险区域的行至少为i，能圈出的最大矩形不超过k的概率，同样这里暂时不算第一个危险行下面的全部安全的概率
我们重写一下状态转移方程
$$dp[i][j]=\sum_{t=0}^{j-1}dp[i][t]*dp[i+1][j-1-t]*(1-q)*q^{j-1-t}+q^jdp[i+1][j]$$
注意从i+1行转移来的所有的部分都要乘上第i行全部安全的概率
这个dp是 $O(k^2)$的，十分优秀了，不需要再考虑优化
我们考虑最后的答案的计算，我们容易想到，如果第一行的某个位置是危险的，那泳池相当于被分成了左右两个独立的部分互相不干扰，否则可以利用上面的dp计算
我们考虑DP[i]表示考虑了前i列，当前圈出的最大矩形面积不大于k的概率
可以写出状态转移方程
$$DP[i]=f[2][i]*q^i+\sum_{j=1}^{k+1}DP[i-j]*(1-q)*dp[2][j-1]*q^{1-j}$$
大概是两种转移，要么当前还没有出现过危险格子，否则枚举最后一个危险格子的位置，然后转移
想到这里就能拿70分，我觉得在考场上已经很不容易了
这个dp是 $O(nk)$的，考虑到n的范围是1e9，显然要往矩乘的方面想，矩乘是很容易做的，但复杂度是 $O(k^3logn)$，还不足以通过此题，但可以拿90分，我觉得是考场上的一个非常好的策略
我们注意到DP的转移前半部分只在i很小的情况下有用，而后面的部分，由于dp已经预处理出来，是一个常系数的齐次线性递推，这个可以用特征多项式来优化到 $O(k^2logn)$，感谢出题人k的范围比较小，否则需要用fft来优化多项式取模以达到 $O(klogklogn)$的复杂度

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

#define LL long long
#define LB long double
#define ull unsigned long long
#define x first
#define y second
#define pb push_back
#define pf push_front
#define mp make_pair
#define Pair pair<int,int>
#define pLL pair<LL,LL>
#define pii pair<double,double>
#define LOWBIT(x) x & (-x)

const int INF=2e9;
const LL LINF=2e16;
const int magic=348;
const int MOD=998244353;
const double eps=1e-10;
const double pi=acos(-1);

inline int getint()
{
    bool f;char ch;int res;
    while (!isdigit(ch=getchar()) && ch!='-') {}
    if (ch=='-') f=false,res=0; else f=true,res=ch-'0';
    while (isdigit(ch=getchar())) res=res*10+ch-'0';
    return f?res:-res;
}

const int MAXN=2000;
const int Height=1001;

inline int add(int x) {if (x>=MOD) x-=MOD;return x;}
inline int sub(int x) {if (x<0) x+=MOD;return x;}

inline int quick_pow(int x,int y)
{
    int res=1;
    while (y)
    {
        if (y&1) res=(1ll*res*x)%MOD,y--;
        x=(1ll*x*x)%MOD;y>>=1;
    }
    return res;
}

int n,k,q;
int dp[MAXN+48][MAXN+48],DP[MAXN+48];
int pw[MAXN+48];
int a[MAXN+48],g[MAXN+48],inv,K;

inline void construct_g()
{
    int i;
    for (i=1;i<=K;i++) a[i]=((1ll*sub(1-q)*pw[i-1])%MOD*dp[2][i-1])%MOD;
    g[K]=1;for (i=K;i>=1;i--) g[K-i]=sub(-a[i]);
    if (K&1) for (i=0;i<=K;i++) g[i]=sub(-g[i]);
    inv=quick_pow(g[K],MOD-2);
}

struct poly
{
    int a[MAXN*4+48],len;
    inline void clear() {for (register int i=0;i<len;i++) a[i]=0;}
    inline poly operator * (poly other)
    {
        poly res;res.len=len+other.len-1;res.clear();
        for (register int i=0;i<len;i++)
            for (register int j=0;j<other.len;j++)
                res.a[i+j]=add(res.a[i+j]+(1ll*a[i]*other.a[j])%MOD);
        for (register int i=res.len-1;i>=K;i--)
        {
            int tmp=(1ll*res.a[i]*inv)%MOD;
            for (register int j=i;j>=i-K;j--)
                res.a[j]=sub(res.a[j]-(1ll*g[j-i+K]*tmp)%MOD);
        }
        if (res.len>K) res.len=K;
        return res;
    }
};

inline poly quick_pow(poly x,int y)
{
    poly res;res.len=1;res.a[0]=1;
    while (y)
    {
        if (y&1) res=res*x,y--;
        x=x*x;y>>=1;
    }
    return res;
}

inline int solve(int limarea)
{
    int i,j,t;
    for (i=0;i<=Height+1;i++)
        for (j=0;j<=limarea;j++)
            dp[i][j]=0;
    dp[Height+1][0]=1;
    for (i=Height;i>=2;i--)
    {
        dp[i][0]=1;
        for (j=1;(i-1)*j<=limarea;j++)
        {
            for (t=0;t<=j-1;t++)
                dp[i][j]=add(dp[i][j]+(((1ll*dp[i][t]*sub(1-q))%MOD*dp[i+1][j-1-t])%MOD*pw[j-1-t])%MOD);
            dp[i][j]=add(dp[i][j]+(1ll*dp[i+1][j]*pw[j])%MOD);
        }
    }
    /*
    if (n<=MAXN)
    {
        DP[0]=1;
        for (i=1;i<=n;i++)
        {
            if (i<=limarea) DP[i]=(1ll*dp[2][i]*pw[i])%MOD; else DP[i]=0;
            for (j=0;j<=min(limarea,i-1);j++) DP[i]=add(DP[i]+(((1ll*DP[i-j-1]*sub(1-q))%MOD*dp[2][j])%MOD*pw[j])%MOD);
        }
        return DP[n];
    }
    */
    DP[0]=1;
    for (i=1;i<=limarea;i++)
    {
        DP[i]=(1ll*dp[2][i]*pw[i])%MOD;
        for (j=0;j<=min(limarea,i-1);j++) DP[i]=add(DP[i]+(((1ll*DP[i-j-1]*sub(1-q))%MOD*dp[2][j])%MOD*pw[j])%MOD);
    }
    K=limarea+1;construct_g();
    poly base;base.len=2;base.clear();base.a[1]=1;
    poly res=quick_pow(base,n);int ans=0;
    for (i=0;i<=K-1;i++) ans=add(ans+(1ll*res.a[i]*DP[i])%MOD);
    return ans;
}

int main ()
{
    int i,x,y;
    n=getint();k=getint();x=getint();y=getint();
    q=(1ll*x*quick_pow(y,MOD-2))%MOD;
    pw[0]=1;for (i=1;i<=MAXN;i++) pw[i]=(1ll*pw[i-1]*q)%MOD;
    printf("%d\n",sub(solve(k)-solve(k-1)));
    return 0;
}