結果

問題 No.2459 Stampaholic (Hard)
ユーザー srjywrdnprktsrjywrdnprkt
提出日時 2023-08-16 09:49:56
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 415 ms / 4,000 ms
コード長 6,799 bytes
コンパイル時間 3,832 ms
コンパイル使用メモリ 250,236 KB
実行使用メモリ 39,132 KB
最終ジャッジ日時 2024-06-09 14:45:24
合計ジャッジ時間 9,607 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge2
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 16 ms
11,172 KB
testcase_01 AC 415 ms
39,004 KB
testcase_02 AC 105 ms
17,084 KB
testcase_03 AC 15 ms
11,188 KB
testcase_04 AC 14 ms
11,144 KB
testcase_05 AC 15 ms
11,080 KB
testcase_06 AC 15 ms
11,076 KB
testcase_07 AC 15 ms
11,148 KB
testcase_08 AC 200 ms
22,296 KB
testcase_09 AC 107 ms
17,524 KB
testcase_10 AC 406 ms
36,152 KB
testcase_11 AC 206 ms
24,604 KB
testcase_12 AC 413 ms
38,464 KB
testcase_13 AC 403 ms
36,136 KB
testcase_14 AC 108 ms
18,236 KB
testcase_15 AC 414 ms
39,004 KB
testcase_16 AC 413 ms
39,008 KB
testcase_17 AC 415 ms
39,132 KB
testcase_18 AC 411 ms
39,004 KB
testcase_19 AC 409 ms
39,132 KB
testcase_20 AC 15 ms
11,108 KB
testcase_21 AC 399 ms
33,916 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <bits/stdc++.h>
#include <atcoder/modint>
#include <atcoder/convolution>

using namespace std;
using namespace atcoder;
using ll = long long;

using mint = modint998244353;

const ll MAX = 1e6+10;
vector<mint> f, finv;

mint inv(mint x){
    mint ans = 1;
    ll e = 998244351;

    while (e > 0){
        if ((e & 1LL)) ans *= x;
        e = e >> 1LL;
        x *= x;
    }

    return ans;
}

void init(){
    f.resize(MAX+1); finv.resize(MAX+1);
    f[0] = 1;
    for (int i=1; i<=MAX; i++) f[i] = f[i-1]*i;
    finv[MAX] = inv(f[MAX]);
    for (int i=MAX-1; i>=0; i--) finv[i] = finv[i+1] * (i+1);
}

mint C(ll n, ll k){
    if (n < k || k < 0) return 0;

    return f[n] * finv[k] * finv[n-k] ;
}

template<class T>
struct FormalPowerSeries : vector<T> {
    using vector<T>::vector;
    using vector<T>::operator=;
    using F = FormalPowerSeries;

    F operator-() const {
        F res(*this);
        for (auto &e : res) e = -e;
        return res;
    }
    F &operator*=(const T &g) {
        for (auto &e : *this) e *= g;
        return *this;
    }
    F &operator/=(const T &g) {
        assert(g != T(0));
        *this *= g.inv();
        return *this;
    }
    F &operator+=(const F &g) {
        int n = (*this).size(), m = g.size();
        for (int i=0; i<min(n, m); i++) (*this)[i] += g[i];
        return *this;
    }
    F &operator-=(const F &g) {
        int n = (*this).size(), m = g.size();
        for (int i=0; i<min(n, m); i++) (*this)[i] -= g[i];
        return *this;
    }
    F &operator<<=(const int d) {
        int n = (*this).size();
        (*this).insert((*this).begin(), d, 0);
        (*this).resize(n);
        return *this;
    }
    F &operator>>=(const int d) {
        int n = (*this).size();
        (*this).erase((*this).begin(), (*this).begin() + min(n, d));
        (*this).resize(n);
        return *this;
    }
    F inv(int d = -1) const {
        int n = (*this).size();
        assert(n != 0 && (*this)[0] != 0);
        if (d == -1) d = n;
        assert(d > 0);
        F res{(*this)[0].inv()};
        while (res.size() < d) {
            int m = size(res);
            F f(begin(*this), begin(*this) + min(n, 2*m));
            F r(res);
            f.resize(2*m), internal::butterfly(f);
            r.resize(2*m), internal::butterfly(r);
            for (int i=0; i<2*m; i++) f[i] *= r[i];
            internal::butterfly_inv(f);
            f.erase(f.begin(), f.begin() + m);
            f.resize(2*m), internal::butterfly(f);
            for (int i=0; i<2*m; i++) f[i] *= r[i];
            internal::butterfly_inv(f);
            T iz = T(2*m).inv(); iz *= -iz;
            for (int i=0; i<m; i++) f[i] *= iz;
            res.insert(res.end(), f.begin(), f.begin() + m);
        }
        return {res.begin(), res.begin() + d};
    }

    // fast: FMT-friendly modulus only
    F &operator*=(const F &g) {
        int n = (*this).size();
        *this = convolution(*this, g);
        (*this).resize(n);
        return *this;
    }
    F &operator/=(const F &g) {
        int n = (*this).size();
        *this = convolution(*this, g.inv(n));
        (*this).resize(n);
        return *this;
    }

    // sparse
    F &operator*=(vector<pair<int, T>> g) {
        int n = (*this).size();
        auto [d, c] = g.front();
        if (d == 0) g.erase(g.begin());
        else c = 0;
        for (int i=n-1; i>=0; i--){
            (*this)[i] *= c;
            for (auto &[j, b] : g) {
                if (j > i) break;
                (*this)[i] += (*this)[i-j] * b;
            }
        }
        return *this;
    }
    F &operator/=(vector<pair<int, T>> g) {
        int n = (*this).size();
        auto [d, c] = g.front();
        assert(d == 0 && c != T(0));
        T ic = c.inv();
        g.erase(g.begin());
        for (int i=0; i<n; i++){
            for (auto &[j, b] : g) {
            if (j > i) break;
            (*this)[i] -= (*this)[i-j] * b;
            }
            (*this)[i] *= ic;
        }
        return *this;
    }

    // multiply and divide (1 + cz^d)
    void multiply(const int d, const T c) {
        int n = (*this).size();
        if (c == T(1)) for (int i=n-d-1; i>=0; i--) (*this)[i+d] += (*this)[i];
        else if (c == T(-1)) for (int i=n-d-1; i>=0; i--) (*this)[i+d] -= (*this)[i];
        else for (int i=n-d-1; i>=0; i--) (*this)[i+d] += (*this)[i] * c;
    }
    void divide(const int d, const T c) {
        int n = (*this).size();
        if (c == T(1)) for (int i=0; i<n-d; i++) (*this)[i+d] -= (*this)[i];
        else if (c == T(-1)) for (int i=0; i<n-d; i++) (*this)[i+d] += (*this)[i];
        else for (int i=0; i<n-d; i++) (*this)[i+d] -= (*this)[i] * c;
    }

    T eval(const T &a) const {
        T x(1), res(0);
        for (auto e : *this) res += e * x, x *= a;
        return res;
    }

    F operator*(const T &g) const { return F(*this) *= g; }
    F operator/(const T &g) const { return F(*this) /= g; }
    F operator+(const F &g) const { return F(*this) += g; }
    F operator-(const F &g) const { return F(*this) -= g; }
    F operator<<(const int d) const { return F(*this) <<= d; }
    F operator>>(const int d) const { return F(*this) >>= d; }
    F operator*(const F &g) const { return F(*this) *= g; }
    F operator/(const F &g) const { return F(*this) /= g; }
    F operator*(vector<pair<int, T>> g) const { return F(*this) *= g; }
    F operator/(vector<pair<int, T>> g) const { return F(*this) /= g; }
};

using mint = modint998244353;
using fps = FormalPowerSeries<mint>;

int main(){

    init();

    ll H, W, N, K, hx, wx, hy, wy;
    mint c, ans, M;
    cin >> H >> W >> N >> K;
    vector<mint> hps(N+1), wps(N+1), hs(N+1), ws(N+1);
    fps a(N+1), b(N+1);

    if (K*2 <= H){
        hx = K-1; hy = H-K*2+2;
    } 
    else{
        hx = H-K; hy = H-(H-K)*2;
    }
    if (K*2 <= W){
        wx = K-1; wy = W-K*2+2;
    }   
    else{
        wx = W-K; wy = W-(W-K)*2;
    } 

    c = 1;
    for (int i=0; i<=N; i++){
        b[i] = finv[i+1];
        c *= hx+1;
        a[i] = c * finv[i+1];
    }
    a *= b.inv();
    for (int i=0; i<=N; i++) hps[i] = a[i] * f[i];
    hps[0] -= 1;

    c = 1;
    for (int i=0; i<=N; i++){
        b[i] = finv[i+1];
        c *= wx+1;
        a[i] = c * finv[i+1];
    }
    a *= b.inv();
    for (int i=0; i<=N; i++) wps[i] = a[i] * f[i];
    wps[0] -= 1;

    ans = (hx * 2 + hy) * (wx * 2 + wy); 

    c = 1;
    for (int i=0; i<=N; i++){
        hs[i] = hps[i] * 2 + c * hy;
        c *= hx+1;
    }

    c = 1;
    for (int i=0; i<=N; i++){
        ws[i] = wps[i] * 2 + c * wy;
        c *= wx+1;
    }

    M = -mint((H-K+1) * (W-K+1)).inv();
    c = 1;
    for (int i=0; i<=N; i++){
        ans -= c * C(N, i) * hs[i] * ws[i];
        c *= M;
    }

    cout << ans.val() << endl;

    return 0;
}
0