题解 P3824 【[NOI2017]泳池】

· · 题解

首先，恰好等于k的概率可以通过差分转化为小于等于k的概率减去小于等于k-1的概率。

考虑怎么求小于等于k的概率。设f_n为底部宽为n时的答案，那么我们可以枚举最后连续的一段白色的长度，即\displaystyle f_n=\sum_{i=0}^{k}f_{n-i-1}\cdot\sum_{j=2}^{\lfloor{k/i}\rfloor+1}dp_{i,j}。其中dp_{i,j}表示一个宽为i的矩形，高度为1至j-1都全部是白色，高度为j的一行开始出现障碍，并且包含的所有贴住底部的白色矩形(包括宽不一定为i和高超过j的那些)大小都不超过k的概率。记g_{i,j}=\sum_{l\geq j}dp_{i,l}。

如果我们可以求出dp和g数组，那么转移方程就变成了齐次常系数线性递推的形式：\displaystyle f_n=\sum_{i=0}^{k}f_{n-i-1}\cdot g_{i,2}。注意到k比较大，传统的矩阵乘法为O(k^3\log n)无法通过，可以用多项式取模的方法来做。

先考虑怎么求出dp和g。枚举高度为j处的第一个障碍，

\displaystyle dp_{i,j}=[i\cdot(j-1)\leq k](1-p)p^{j-1}\sum_{l=1}^{i}\left(\sum_{q>j}dp_{l-1,q}\right)\left(\sum_{q\geq j}dp_{i-l,q}\right)

前缀和优化一下，就可以做到O(k^2\log k)了，因为合法的状态中要求i\cdot(j-1)\leq k，所以一共有\sum_{i=1}^{k}\lfloor{k/i}\rfloor=O(k\log k)个非零状态，每次转移为O(k)，已经可以通过此题。实际上还可以用FFT优化到k\log^2k。

~

下面是用多项式取模优化齐次常系数线性递推的方法：

对于递推f_n=\sum_{i=1}^ka_kf_{n-k}

设矩阵乘法中的转移矩阵为A，设初始状态矩阵为B[f_k,f_{k-1},...,f_1]，我们现在要求BA^n，即求A^n。

考虑A的特征多项式\lambda(x)=det(xI-A)，I为单位矩阵。由于转移矩阵的特殊性，在此处我们可以手动求一下特征多项式，\lambda(x)=x^k-\sum_{i=1}^ka_ix^{k-i}，是一个k次多项式，并且满足\lambda(A)=0。

如果我们把A^n表示成若干倍的\lambda(A)，即A^n=\lambda(A)g(A)+r(A)。由于\lambda(A)=0，所以A^n=r(A)。而r(A)=A^n\mod\lambda(A)，我们只需要求出x^n模特征多项式后的多项式，然后将A带入即可求出A^n。

~

多项式取模可以使用Picks的方法，用FFT做到k\log k，加上多项式快速幂复杂度就是O(k\log k\log n)，边乘边取模。当然在这里没有必要，直接暴力取模即可。

~

如何暴力取模?假设要求f(x)\mod g(x)，设f(x)=g(x)h(x)+r(x)，我们就是要求r(x)，并且要求r(x)的次数小于g(x)的次数。不妨将g(x)看成0，那么此时f(x)就等于r(x)。设g(x)=\sum_{i=0}^kb_ix^i=0，那么我们可以得出\displaystyle x^k=-\sum_{i=0}^{k-1}\frac{b_i}{b_k}x^i，我们只需要把f(x)中系数过大的一部分迭代，即可将系数一步步压到比k小。

两个k次多项式相乘，超出的部分系数是O(k)的，所以暴力取模的复杂度为O(k^2)，套上快速幂就是O(k^2\log n)，可以通过此题。

~

然后我们考虑怎么统计答案。我们已经用一个关于A的k-1次多项式A^n=R(A)=\sum_{i=0}^{k-1}r_iA^i表示出了A^n。

根据矩阵乘法的分配律，B\cdot R(A)=\sum_{i=0}^{k-1}r_i\cdot B\cdot A^i。我们要的结果是矩阵的第一个元素，而B\cdot A^i的第一个元素就是f_{k+i}，所以我们要先求出f_k,f_{k+1},...,f_{2k}，然后直接扫一遍统计答案即可。

~

复杂度O(k^2(\log k+\log n))，其实可以做到O(k\log k(\log k+\log n))。

我代码中的dp部分和上面写的转移不太一样。

#include<bits/stdc++.h>

#define For(i,_beg,_end) for(int i=(_beg),i##end=(_end);i<=i##end;++i)
#define Rep(i,_beg,_end) for(int i=(_beg),i##end=(_end);i>=i##end;--i)

template<typename T>T Max(const T &x,const T &y){return x<y?y:x;}
template<typename T>T Min(const T &x,const T &y){return x<y?x:y;}
template<typename T>int chkmax(T &x,const T &y){return x<y?(x=y,1):0;}
template<typename T>int chkmin(T &x,const T &y){return x>y?(x=y,1):0;}
template<typename T>void read(T &x){
    T f=1;char ch=getchar();
    for(;ch<'0'||ch>'9';ch=getchar())if(ch=='-')f=-1;
    for(x=0;ch>='0'&&ch<='9';ch=getchar())x=x*10+ch-'0';
    x*=f;
}

typedef long long LL;
const int N=1010,mod=998244353;
int n,m;
LL a,b,p,q,dp[N][N],g[N][N],pw[N];
LL A[N],f[N<<1];

LL power(LL,LL);
LL Solve(int);

int main(){
    read(n);read(m);read(a);read(b);
    p=a*power(b,mod-2)%mod;q=(mod+1-p)%mod;
    pw[0]=1;
    For(i,1,m) pw[i]=pw[i-1]*p%mod;

    printf("%lld\n",(Solve(m)-Solve(m-1)+mod)%mod);
    return 0;
}

LL Solve(int k){
    memset(dp,0,sizeof dp);
    memset(g,0,sizeof g);
    For(i,1,k+2) g[0][i]=dp[0][i]=1;
    For(i,1,k) Rep(j,k/i+1,2){
        For(l,1,i) dp[i][j]=(dp[i][j]+g[l-1][j+1]*g[i-l][j]%mod*pw[l-1]%mod*q)%mod;
        g[i][j]=(g[i][j+1]*pw[i]+dp[i][j])%mod;
    }
    memset(A,0,sizeof A);
    For(i,0,k) A[i+1]=q*g[i][2]%mod*pw[i]%mod;
    memset(f,0,sizeof f);
    f[0]=1;
    For(i,1,k){
        f[i]=g[i][2]*pw[i]%mod;
        For(j,1,i) f[i]=(f[i]+A[j]*f[i-j])%mod;
    }
    k++;
    For(i,k,k<<1) For(j,1,k) f[i]=(f[i]+A[j]*f[i-j])%mod;
    if(n<=k)return f[n];
    int y=n-k,len=1,L=0;
    LL res[N<<2],tmp[N<<2],x[N<<2];
    memset(res,0,sizeof res);
    memset(x,0,sizeof x);
    res[0]=1;x[1]=1;
    for(;y;y>>=1){
        if(y&1){
            For(i,0,L+len) tmp[i]=0;
            For(i,0,L) For(j,0,len) tmp[i+j]=(tmp[i+j]+res[i]*x[j])%mod;
            L+=len;
            Rep(i,L,k) For(j,1,k) tmp[i-j]=(tmp[i-j]+tmp[i]*A[j])%mod;
            chkmin(L,k-1);
            For(i,0,L) res[i]=tmp[i];
        }
        For(i,0,len+len) tmp[i]=0;
        For(i,0,len) For(j,0,len) tmp[i+j]=(tmp[i+j]+x[i]*x[j])%mod;
        len<<=1;
        Rep(i,len,k) For(j,1,k) tmp[i-j]=(tmp[i-j]+tmp[i]*A[j])%mod;
        chkmin(len,k-1);
        For(i,0,len) x[i]=tmp[i];
    }
    LL ans=0;
    For(i,0,k-1) ans=(ans+res[i]*f[i+k])%mod;
    return ans;
}
LL power(LL x,LL y){
    LL res=1;
    for(;y;y>>=1,x=x*x%mod) if(y&1) res=res*x%mod;
    return res;
}