原文链接www.cnblogs.com/zhouzhendong/p/UOJ316.html

题解

首先，我们将答案转化成最大矩形大小 \(\leq k\) 的概率 减去 \(\leq k-1\) 的概率。

然后我们考虑 DP。

设 \(dp[i][j]\) 表示矩形宽度为 \(j\) ，当前已知最底下 \(i\) 行是安全的，在这个情况下，最大安全区域 \(\leq k\) 或 \(\leq k-1\) 的概率。

状态的转移分两种：一种是第 \(i+1\) 层全部安全，一种是枚举第 \(i+1\) 层的第一个不安全点。

\[ dp[i][j] = q ^ {j} dp[i+1][j] + \sum_{k=1}^j (1-q)q^{k-1} dp[i+1][k-1] \cdot dp[i][j-k] \]

至此，我们得到了 70 分的做法。

接下来进入鬼畜时间：

通过打(bai)表(du)，我们可以得知，数列 \(dp[0][1],dp[0][2],\cdots ,dp[0][n]\) 的最短线性递推式很短，几千就够了。

大力上一波 BM 求个递推式，然后大力套用 Cayley-Hamilton定理 求出第 \(n\) 项即可。

代码

#include <bits/stdc++.h>
#define clr(x) memset(x,0,sizeof x)
#define For(i,a,b) for (int i=(a);i<=(b);i++)
#define Fod(i,b,a) for (int i=(b);i>=(a);i--)
#define fi first
#define se second
#define pb(x) push_back(x)
#define mp(x,y) make_pair(x,y)
#define outval(x) printf(#x" = %d\n",x)
#define outtag(x) puts("---------------"#x"---------------")
#define outarr(a,L,R) printf(#a"[%d..%d] = ",L,R);                        For(_x,L,R)printf("%d ",a[_x]);puts("")
using namespace std;
typedef long long LL;
LL read(){
    LL x=0,f=0;
    char ch=getchar();
    while (!isdigit(ch))
        f|=ch=='-',ch=getchar();
    while (isdigit(ch))
        x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return f?-x:x;
}
const int N=2005*2,mod=998244353;
int Pow(int x,int y){
    int ans=1;
    for (;y;y>>=1,x=(LL)x*x%mod)
        if (y&1)
            ans=(LL)ans*x%mod;
    return ans;
}
void Add(int &x,int y){
    if ((x+=y)>=mod)
        x-=mod;
}
void Del(int &x,int y){
    if ((x-=y)<0)
        x+=mod;
}
int n,k,p,ip;
int getp(){
    int x=read(),y=read();
    return (LL)x*Pow(y,mod-2)%mod;
}
int pwp[N];
int dp[N][N];
vector <int> BM(int *a,int n){
    static vector <int> R[N];
    static int Fail[N],d[N];
    clr(Fail),clr(d);
    R[0].clear();
    int c=0;
    For(i,1,n){
        d[i]=a[i];
        For(j,0,(int)R[c].size()-1)
            Del(d[i],(LL)R[c][j]*a[i-j-1]%mod);
        if (d[i]==0)
            continue;
        Fail[c]=i;
        if (c==0){
            R[++c].clear();
            R[c].resize(i,0);
            continue;
        }
        int id=c-1;
        For(j,0,c-1)
            if (i-Fail[j]+(int)R[j].size()<i-Fail[id]+(int)R[id].size())
                id=j;
        int tmp=(LL)d[i]*Pow(d[Fail[id]],mod-2)%mod;
        R[c+1]=R[c];
        R[c+1].resize(max((int)R[c].size(),i-Fail[id]+(int)R[id].size()),0);
        Add(R[c+1][i-Fail[id]-1],tmp);
        For(j,0,(int)R[id].size()-1)
            Del(R[c+1][i-Fail[id]+j],(LL)tmp*R[id][j]%mod);
        c++;
    }
    return R[c];
}
void Mul(int *x,int *y,int *a,int n){
    static int z[N];
    clr(z);
    For(i,0,n-1)
        For(j,0,n-1)
            Add(z[i+j],(LL)x[i]*y[j]%mod);
    Fod(i,n*2-2,n)
        if (z[i])
            For(j,1,n)
                Add(z[i-j],(LL)a[j]*z[i]%mod);
    For(i,0,n-1)
        x[i]=z[i];
}
int CH(int *a,int *b,int n,int k){
    static int c[N],x[N];
    clr(c),clr(x),c[0]=x[1]=1;
    if (n==1)
        x[1]=0,x[0]=a[1];
    int y=k;
    for (;y;y>>=1,Mul(x,x,a,n))
        if (y&1)
            Mul(c,x,a,n);
    int ans=0;
    For(i,0,n-1)
        Add(ans,(LL)b[i]*c[i]%mod);
    return ans;
}
int calc(int *a,int n,int k){
    vector <int> vb=BM(a,n);
    static int b[N];
    clr(b);
    n=vb.size();
    For(i,0,n-1)
        b[i+1]=vb[i];
    return CH(b,a,n,k);
}
int Solve(int k){
    int lim=min(n,2000);
    clr(dp);
    For(i,0,k+1)
        dp[i][0]=1;
    Fod(i,k,0)
        for (int j=1;j<=lim&&j*i<=k;j++){
            Add(dp[i][j],(LL)dp[i+1][j]*pwp[j]%mod);
            For(t,1,j)
                Add(dp[i][j],(LL)dp[i+1][t-1]*pwp[t-1]%mod*ip%mod*dp[i][j-t]%mod);
        }
    if (n<=lim)
        return dp[0][n];
    return calc(dp[0],lim,n);
}
int main(){
    n=read(),k=read(),p=getp(),ip=(mod+1-p)%mod;
    pwp[0]=1;
    For(i,1,min(n,4000))
        pwp[i]=(LL)pwp[i-1]*p%mod;
    printf("%d\n",(Solve(k)-Solve(k-1)+mod)%mod);
    return 0;
}

UOJ#316. 【NOI2017】泳池 动态规划,Berlekamp-Massey,Cayley-Hamilton定理

原文：https://www.cnblogs.com/zhouzhendong/p/UOJ316.html