「NOI2017」泳池

DP式子比后面的東西難推多了

LOJ2304 

Luogu P3824

UOJ #316

---

題意

給定一個長度為$ n$高為$ \infty$的矩形

每個點有$ 1-P$的概率不可被選擇

求最大的和底邊重合的不包含不可選點的矩形的面積為$ K$的概率

$ n \leq 10^9 k \leq 10^3$

---

題解

K可以出到50000的

首先考慮DP

面積恰好為$ K$的概率可以差分為不高於$ K$的概率減去不高於$ K-1$的概率

設$ f[i][j]$表示長度為$ i$的矩形，從底邊起$ j$行都可選，最大面積不大於$ K$的概率

邊界為$ f[0][j]=1,f[i][j]=0 當且僅當i*j>k$

考慮轉移，要么第$ j+1$行也都可選，要么第$ j+1$行有不可選的位置

對於第二種情況我們枚舉從左到右第一個不可選的位置

有轉移方程式

$$f[i][j]=f[i][j+1]*P^i+\sum_{k=1}^iP^{k-1}(1-P)f[k-1][j+1]·f[i-k][j]$$

我們要求的是$ f[n][0]$

由於$ i*j \leq K$因此復雜度大致是$ n·k \log k$的

可以得$ 70$分

 

容易發現當$ n$遠大於$ k$的時候，每連續$ k$列必然有一列最低端有不可選點

令$ F[i]$表示當$ i>k$時，長度為$ i$的矩形的答案

枚舉從右往左第一個不可選點，有轉移方程式

$$ F[i]=f[i-k]*(1-P)*f[k-1][1]*P^{k-1}$$

這是一個線性遞推的標准形式，可以用特征多項式優化到$ O(k^2 \log n)$甚至$ O(k \log k \log n)$

如果采用后一種的話復雜度的瓶頸在於前面的$ O(k^2)DP$，這部分可以用分治$ NTT$優化

然而我並沒有寫

---

代碼

#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#define p 998244353
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
    ll x=0;char zf=1;char ch=getchar();
    while(ch!='-'&&!isdigit(ch))ch=getchar();
    if(ch=='-')zf=-1,ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
namespace poly{
    vector<int>R;
    int ksm(int x,int y=p-2){
        int ans=1;
        for(rt i=y;i;i>>=1,x=1ll*x*x%p)if(i&1)ans=1ll*ans*x%p;
        return ans;
    }
    void NTT(int n,vector<int>&A,int fla){
        A.resize(n);
        for(rt i=0;i<n;i++)if(i>R[i])swap(A[i],A[R[i]]);
        for(rt i=1;i<n;i<<=1){
            int w=ksm(3,(p-1)/2/i);
            for(rt j=0;j<n;j+=i<<1){
                int K=1;
                for(rt k=0;k<i;k++,K=1ll*K*w%p){
                    int x=A[j+k],y=1ll*K*A[i+j+k]%p;
                    A[j+k]=(x+y)%p,A[i+j+k]=(x-y)%p;
                }
            }
        }
        if(fla==-1){
            reverse(A.begin()+1,A.end());
            int invn=ksm(n);
            for(rt i=0;i<n;i++)A[i]=1ll*A[i]*invn%p;
        }
    }
    vector<int>Mul(vector<int>x,vector<int>y){
        int lim=1,sz=x.size()+y.size()-1;
        while(lim<=sz)lim<<=1;R.resize(lim);
        for(rt i=0;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
        NTT(lim,x,1);NTT(lim,y,1);
        for(rt i=0;i<lim;i++)x[i]=1ll*x[i]*y[i]%p;
        NTT(lim,x,-1);x.resize(sz);
        return x;
    }
    vector<int>sqr(vector<int>x){
        int lim=1,sz=x.size()*2-1;
        while(lim<=sz)lim<<=1;R.resize(lim);
        for(rt i=0;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
        NTT(lim,x,1);for(rt i=0;i<lim;i++)x[i]=1ll*x[i]*x[i]%p;
        NTT(lim,x,-1);x.resize(sz);
        return x;
    }
    vector<int>Inv(vector<int>A,int n=-1){
        if(n==-1)n=A.size();
        if(n==1)return vector<int>(1,ksm(A[0]));
        vector<int>b=Inv(A,(n+1)/2);
        int lim=1;while(lim<=n+n)lim<<=1;R.resize(lim);
        for(rt i=0;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
        A.resize(n);NTT(lim,A,1);NTT(lim,b,1);
        for(rt i=0;i<lim;i++)A[i]=1ll*b[i]*(2ll-1ll*A[i]*b[i]%p)%p;
        NTT(lim,A,-1);A.resize(n);
        return A;
    }
    vector<int>Div(vector<int>A,vector<int>B){
        int n=A.size(),m=B.size();
        reverse(A.begin(),A.end());
        reverse(B.begin(),B.end());
        A.resize(n-m+1),B.resize(n-m+1);
        int lim=1;while(lim<=2*(n-m+1))lim<<=1;R.resize(lim);
        for(rt i=0;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
        vector<int>ans=Mul(A,Inv(B));ans.resize(n-m+1);
        reverse(ans.begin(),ans.end());
        return ans;
    }
    vector<int>add(vector<int>A,vector<int>B){
        int len=max(A.size(),B.size());A.resize(len+1);
        for(rt i=0;i<=len;i++)(A[i]+=B[i])%=p;
        return A;
    }
    vector<int>sub(vector<int>A,vector<int>B){
        int len=max(A.size(),B.size());A.resize(len+1);
        for(rt i=0;i<=len;i++)(A[i]-=B[i])%=p;
        return A;
    }
    vector<int>Mod(vector<int>x,vector<int>y){
        if(x.size()<=y.size())return x;
          vector<int>ans=Div(x,y);
        ans=sub(x,Mul(y,ans));
        while(!ans[ans.size()-1])ans.pop_back();
        if(ans.size()>y.size())ans.resize(y.size());
        return ans;
    }
}
using namespace poly;
int a[100010];
vector<int>fmo;
vector<int>ksm(vector<int>x,int y){
    if(y==1)return x;
    vector<int>ans=Mod(sqr(ksm(x,y>>1)),fmo);
    if(y&1){
        ans.push_back(0);
        for(rt i=ans.size()-2;i>=0;i--)ans[i+1]=ans[i],ans[i]=0;
    }
    return ans;
}
using namespace poly;
int k,m,n,x,y,z,cnt,ans,K,P;
int f[1010][1010],mi[1010];
//f[i][j]長度為i合法高度最低至少為j的合法概率 
void calc(int n,int K,int fla){
    memset(f,0,sizeof(f));
    for(rt i=0;i<=K+10;i++)f[0][i]=1;
    for(rt i=1;i<=K+2;i++)
    for(rt j=K/i;j>=0;j--){
        f[i][j]=1ll*f[i][j+1]*mi[i]%p;
        for(rt k=1;k<=i;k++)(f[i][j]+=1ll*f[k-1][j+1]*mi[k-1]%p*(1+p-P)%p*f[i-k][j]%p)%=p;
    }
    fmo.resize(K+2);
    for(rt j=1;j<=K+1;j++)fmo[K+1-j]=-1ll*mi[j-1]*f[j-1][1]%p*(p+1-P)%p;
    fmo[K+1]=1;int ret=0;
    if(n<=K+1)ret=f[n][0];else {
        vector<int>x;x.push_back(0);x.push_back(1);
        x=ksm(x,n);
        for(rt i=0;i<=K+1;i++)(ret+=1ll*f[i][0]*x[i]%p)%=p;
    }
    (ans+=ret*fla)%=p;
}
int main(){
//    file("pool");
    n=read();K=read();x=read();y=read();P=1ll*x*ksm(y)%p;
    mi[0]=1;for(rt i=1;i<=K+2;i++)mi[i]=1ll*mi[i-1]*P%p;
    calc(n,K,1);calc(n,K-1,-1);cout<<(ans+p)%p;
    return 0;
}