题解

（请在此补充题目的中文题解与思路描述。）

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

const int MOD = 998244353;
using int64 = long long;

int addmod(int a, int b){ int s = a+b; if(s>=MOD) s-=MOD; return s; }
int submod(int a, int b){ int s = a-b; if(s<0) s+=MOD; return s; }
int mulmod(int64 a, int64 b){ return int(a*b%MOD); }
int powmod(int a, long long e){
    int r=1;
    while(e){
        if(e&1) r=mulmod(r,a);
        a=mulmod(a,a);
        e>>=1;
    }
    return r;
}
int invmod(int a){ return powmod(a, MOD-2); }

/* --------------------------------------------------------------- */
/*  compute F(T) = probability that maximal rectangle area ≤ T      */

int compute_F(long long N, int T, int q){          // q = x * inv(y) % MOD
    if (T==0){
        // all heights must be 0 -> each column height 0 with prob (1-q)
        int p0 = submod(1, q);
        return powmod(p0, N);
    }

    // pre‑compute powers of q up to T+1
    vector<int> powq(T+2);
    powq[0]=1;
    for(int i=1;i<=T+1;i++) powq[i]=mulmod(powq[i-1], q);
    int q_Tplus1 = powq[T+1];               // q^{T+1}

    // arrays for a = 1..T
    vector<int> pLow(T+2), pEq(T+2), pGeTrunc(T+2), L(T+2);
    for(int a=1;a<=T;a++){
        pLow[a] = submod(1, powq[a]);               // 1 - q^a
        pEq[a]  = mulmod(powq[a], submod(1, q));    // q^a * (1-q)
        pGeTrunc[a] = submod(powq[a], q_Tplus1);    // q^a - q^{T+1}
        L[a] = T / a;
    }
    // L[T+1] = 0 (not used)
    L[T+1]=0;

    // ----------  table R[a][len]  ----------
    vector<vector<int>> R(T+2);               // R[a][len] for 0 ≤ len ≤ L[a]
    R[T+1].assign(1,1);                       // R[T+1][0]=1

    for(int a=T;a>=1;a--){
        int maxlen = L[a];
        R[a].assign(maxlen+1,0);
        R[a][0]=1;
        for(int n=1;n<=maxlen;n++){
            int64 cur = mulmod(pEq[a], R[a][n-1]);   // first column = a
            int maxEll = min( (a+1<=T? L[a+1]:0), n);
            for(int ell=1; ell<=maxEll; ++ell){
                int high = R[a+1][ell];               // probability of a high block
                if (ell < n){
                    cur += (int64)high * pEq[a] % MOD * R[a][n-ell-1] % MOD;
                }else{
                    cur += high;
                }
                if(cur>= (int64)MOD*MOD) cur %= MOD;
            }
            R[a][n] = int(cur % MOD);
        }
    }

    // ----------  initial values of U[1] ----------
    int d = T+1;                     // order of the recurrence
    vector<int> init(d,0);           // U[0..T]
    init[0]=1;
    for(int n=1;n<=T;n++){
        int64 cur = mulmod(pLow[1], init[n-1]);
        int maxEll = min(L[1], n);
        for(int ell=1; ell<=maxEll; ++ell){
            int r = R[1][ell];
            if (ell < n){
                cur += (int64)r * pLow[1] % MOD * init[n-ell-1] % MOD;
            }else{
                cur += r;
            }
            if(cur>= (int64)MOD*MOD) cur %= MOD;
        }
        init[n] = int(cur % MOD);
    }

    if (N <= T) return init[N];

    // ----------  recurrence coefficients ----------
    // U[n] = Σ_{i=0}^{d-1} coeff[i] * U[n-i-1]
    vector<int> coeff(d,0);
    coeff[0] = pLow[1];
    for(int k=1;k<=T;k++){
        coeff[k] = mulmod(pLow[1], R[1][k]);
    }

    // ----------  Kitamasa :  compute x^N (mod characteristic) ----------
    auto combine = [&](const vector<int>& A, const vector<int>& B){
        vector<int64> tmp(2*d-1,0);
        for(int i=0;i<d;i++) if(A[i])
            for(int j=0;j<d;j++) if(B[j]){
                tmp[i+j] += (int64)A[i]*B[j];
                if(tmp[i+j] >= (int64)MOD*MOD) tmp[i+j] %= MOD;
            }
        // reduction
        for(int i=2*d-2;i>=d;i--){
            if(tmp[i]==0) continue;
            int64 v = tmp[i] % MOD;
            for(int j=0;j<d;j++){
                tmp[i-1-j] += v * coeff[j];
                if(tmp[i-1-j] >= (int64)MOD*MOD) tmp[i-1-j] %= MOD;
            }
        }
        vector<int> res(d,0);
        for(int i=0;i<d;i++) res[i] = int(tmp[i]%MOD);
        return res;
    };

    // binary exponentiation of polynomial x
    vector<int> base(d,0), result(d,0);
    base[1]=1;               // x
    result[0]=1;             // 1
    long long e = N;
    while(e){
        if(e&1) result = combine(result, base);
        base = combine(base, base);
        e >>= 1;
    }
    // result now holds coefficients of x^N  (degree < d)
    int64 ans = 0;
    for(int i=0;i<d;i++){
        ans += (int64)result[i] * init[i];
        if(ans >= (int64)MOD*MOD) ans %= MOD;
    }
    return int(ans % MOD);
}

/* --------------------------------------------------------------- */

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    long long N;
    int K, x, y;
    if(!(cin>>N>>K>>x>>y)) return 0;
    int q = mulmod(x, invmod(y));          // q = x / y  (mod MOD)

    int FK = compute_F(N, K, q);
    int FK1 = (K==0 ? 0 : compute_F(N, K-1, q));
    int answer = submod(FK, FK1);
    cout << answer << "\n";
    return 0;
}