Luogu3824 [NOI2017]泳池 【多项式取模】【递推】【矩阵快速幂】

题目分析：

用数论分块的思想，就会发现其实就是连续一段的长度$i$的高度不能超过$lfloor frac{k}{i} floor$,然后我们会发现最长的非$0$一段不会超过$k$，所以我们可以弄一个长度为$i$的非$0$段的个数称为"元"，然后用
"元"去递推。

这个"元"的求法用DP：令数论分块之后第$i$段的长度为$g[i]$

$$f[i][j] = f[i-1][j] + f[i-1][k]*f[i][j-k-1]*g[i]$$

$$f[1][1] = 1$$

然后接下来就是一个很简单的矩阵快速幂了。

对于$k$很大的情况用多项式取模优化就行了。

代码：

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

const int mod = 998244353;

int n,k,x,y;
int maxx[1060],pos[1060],cas[150],num;
int g[120][1060],f[1060];

int FF[1060],ans[1060],dt[1060],zeta[2060];

int fast_pow(int now,int pw){
    int ans = 1,dt = now,bit = 1;
    while(bit <= pw){
    if(bit & pw) ans = 1ll*ans*dt%mod;
    bit<<=1; dt = 1ll*dt*dt%mod;
    }
    return ans;
}

void multi(int *A,int *B,int len){
    memset(zeta,0,sizeof(zeta));
    for(int i=0;i<len;i++){
    for(int j=0;j<len;j++) zeta[i+j]+=1ll*A[i]*B[j]%mod,zeta[i+j]%=mod;
    }
    for(int i=len*2;i>=len;i--){
    if(zeta[i] == 0) continue;
    int pp = zeta[i],dr = (FF[len]==1?1:mod-1);
    for(int j=i,kk=len;kk>=0;kk--,j--){
        zeta[j] -= 1ll*pp*FF[kk]%mod*dr%mod;
        if(zeta[j] <0) zeta[j] += mod;
    }
    }
    for(int i=0;i<len;i++) A[i] = zeta[i];
}

void fastpow(int len,int pw){
    memset(ans,0,sizeof(ans)); ans[0] = 1;
    memset(dt,0,sizeof(dt)); dt[1] = 1;
    int bit = 1;
    while(bit<=pw){
    if(bit & pw) multi(ans,dt,len);
    multi(dt,dt,len); bit<<=1;
    }
}

void getbase(int now){
    memset(g,0,sizeof(g));
    g[0][0] = 1;
    for(int i=1;i<=num;i++){
    g[i][0] = 1;
    for(int j=1;j<=now/cas[i];j++){
        g[i][j] += g[i-1][j]; g[i][j] %= mod;
        for(int l=0;l<j;l++){
        for(int ec = cas[i+1]+1;ec<=cas[i];ec++){
            g[i][j]+=1ll*g[i-1][l]%mod*g[i][j-l-1]%mod*pos[ec]%mod;
            g[i][j] %= mod;
        }
        }
    }
    }
}

int work(int now){
    memset(cas,0,sizeof(cas)); memset(maxx,0,sizeof(maxx));
    memset(f,0,sizeof(f));
    num = 0;
    for(int i=1;i<=now;i++) {
    maxx[i] = now/i;
    if(maxx[i] != maxx[i-1]) cas[++num] = maxx[i];
    }
    getbase(now);
    f[0] = 1;
    for(int i=1;i<=1000;i++){
    for(int j=0;j<=now;j++){
        if(j > i) break;
        if(j == i){f[i] += g[num][i]; f[i] %= mod;}
        else{f[i]+=1ll*f[i-j-1]*pos[0]%mod*g[num][j]%mod;f[i]%=mod;}
    }
    }
    if(n <= 1000) return f[n];
    else{
    memset(FF,0,sizeof(FF));
    FF[now+1] = 1;
    for(int i=1;i<=now+1;i++)
        FF[now+1-i]=(mod-1ll*g[num][i-1]*pos[0]%mod)%mod;
    if(!(n&1)){for(int i=0;i<=now+1;i++) FF[i] = (mod-FF[i])%mod;}
    fastpow(now+1,n);
    int res = 0;
    for(int i=0;i<=now;i++){res += 1ll*ans[i]*f[i]%mod; res %= mod;}
    return res;
    }
}

int main(){
    scanf("%d%d%d%d",&n,&k,&x,&y);
    x = 1ll*x*fast_pow(y,mod-2)%mod;
    for(int i=0;i<=k;i++){pos[i] = 1ll*fast_pow(x,i)*(mod+1-x)%mod;}
    if(n == 1){printf("%d
",pos[k]);return 0;}
    int ans = work(k)-work(k-1);
    if(ans < 0) ans += mod;
    printf("%d
",ans);
    return 0;
}