#1012632 (C++17) No.2883 K-powered Sum of Fibonacci

• 問題: No.2883 K-powered Sum of Fibonacci
• ユーザー: nonon
• 提出日時: 2024-09-21 12:55:55
• 言語: C++17 (gcc 12.3.0 + boost 1.83.0)
• 結果: AC
• 実行時間: 4 ms / 3,000 ms
• コード長: 15,288 bytes
• コンパイル時間: 3,359 ms
• コンパイル使用メモリ: 225,632 KB
• 実行使用メモリ: 6,944 KB
• 最終ジャッジ日時: 2024-09-21 12:56:00
• 合計ジャッジ時間: 4,080 ms
• ジャッジサーバーID （参考情報）: judge1 / judge3

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	2 ms; 6,816 KB
testcase_01	AC	2 ms; 6,944 KB
testcase_02	AC	1 ms; 6,940 KB
testcase_03	AC	2 ms; 6,944 KB
testcase_04	AC	2 ms; 6,944 KB
testcase_05	AC	2 ms; 6,944 KB
testcase_06	AC	2 ms; 6,940 KB
testcase_07	AC	3 ms; 6,944 KB
testcase_08	AC	2 ms; 6,940 KB
testcase_09	AC	2 ms; 6,944 KB
testcase_10	AC	2 ms; 6,940 KB
testcase_11	AC	3 ms; 6,944 KB
testcase_12	AC	2 ms; 6,940 KB
testcase_13	AC	2 ms; 6,940 KB
testcase_14	AC	4 ms; 6,940 KB
testcase_15	AC	3 ms; 6,944 KB
testcase_16	AC	1 ms; 6,944 KB
testcase_17	AC	2 ms; 6,940 KB
testcase_18	AC	2 ms; 6,944 KB
testcase_19	AC	3 ms; 6,940 KB
testcase_20	AC	3 ms; 6,940 KB
testcase_21	AC	3 ms; 6,940 KB
testcase_22	AC	3 ms; 6,940 KB
testcase_23	AC	3 ms; 6,944 KB
testcase_24	AC	3 ms; 6,940 KB
testcase_25	AC	3 ms; 6,940 KB
testcase_26	AC	4 ms; 6,944 KB
testcase_27	AC	3 ms; 6,940 KB
testcase_28	AC	3 ms; 6,944 KB
testcase_29	AC	3 ms; 6,940 KB
testcase_30	AC	2 ms; 6,940 KB
testcase_31	AC	2 ms; 6,944 KB
testcase_32	AC	2 ms; 6,944 KB
testcase_33	AC	2 ms; 6,944 KB
testcase_34	AC	2 ms; 6,944 KB
testcase_35	AC	2 ms; 6,940 KB
testcase_36	AC	2 ms; 6,940 KB
testcase_37	AC	2 ms; 6,944 KB
testcase_38	AC	2 ms; 6,944 KB
testcase_39	AC	1 ms; 6,940 KB
testcase_40	AC	1 ms; 6,944 KB
testcase_41	AC	1 ms; 6,944 KB
testcase_42	AC	4 ms; 6,940 KB

ソースコード

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

template<int MOD>
struct modint {
    modint() : x(0) {}
    modint(long long v) {
        long long y = v % m();
        if (y < 0) y += m();
        x = (unsigned int)(y);
    }
    static modint raw(int v) {
        modint a;
        a.x = v;
        return a;
    }
    static constexpr int mod() { return m(); }
    unsigned int val() const { return x; }
    modint& operator++() {
        x++;
        if (x == m()) x = 0;
        return *this;
    }
    modint& operator--() {
        if (x == 0) x = m();
        x--;
        return *this;
    }
    modint operator++(int) {
        modint res = *this;
        ++*this;
        return res;
    }
    modint operator--(int) {
        modint res = *this;
        --*this;
        return res;
    }
    modint& operator+=(const modint &r) {
        x += r.x;
        if (x >= m()) x -= m();
        return *this;
    }
    modint& operator-=(const modint &r) {
        x -= r.x;
        if (x >= m()) x += m();
        return *this;
    }
    modint& operator*=(const modint &r) {
        unsigned long long y = x;
        y *= r.x;
        x = (unsigned int)(y % m());
        return *this;
    }
    modint &operator/=(const modint &r) {
        return *this = *this * r.inv();
    }
    friend modint operator+(const modint &a, const modint &b) {
        return modint(a) += b;
    }
    friend modint operator-(const modint &a, const modint &b) {
        return modint(a) -= b;
    }
    friend modint operator*(const modint &a, const modint &b) {
        return modint(a) *= b;
    }
    friend modint operator/(const modint &a, const modint &b) {
        return modint(a) /= b;
    }
    friend bool operator==(const modint &a, const modint &b) {
        return a.x == b.x;
    }
    friend bool operator!=(const modint &a, const modint &b) {
        return a.x != b.x;
    }
    modint operator+() const { return *this; }
    modint operator-() const { return modint() - *this; }
    modint pow(long long k) const {
        assert(k >= 0);
        modint a = *this;
        modint res = 1;
        while (k > 0) {
            if (k & 1) res *= a;
            a *= a;
            k >>= 1;
        }
        return res;
    }
    modint inv() const {
        long long a = x, b = m(), u = 1, v = 0;
        while (b > 0) {
            long long t = a / b;
            a -= t * b;
            swap(a, b);
            u -= t * v;
            swap(u, v);
        }
        return modint(u);
    }
private:
    unsigned int x;
    static constexpr int m() { return MOD; }
};

template<typename mint>
struct Number_Theoretic_Transform {
    static vector<mint> dw, dw_inv;
    static int log;
    static mint root;
    static void ntt(vector<mint> &f) {
        init();
        int n = f.size();
        for (int m = n; m >>= 1;) {
            mint w = 1;
            for (int s = 0, k = 0; s < n; s += (m << 1)) {
                for (int i = s, j = s + m; i < s + m; i++, j++) {
                    mint a = f[i], b = f[j] * w;
                    f[i] = a + b;
                    f[j] = a - b;
                }
                w *= dw[__builtin_ctz(++k)];
            }
        }
    }
    static void intt(vector<mint> &f) {
        init();
        int n = f.size();
        for (int m = 1; m < n; m <<= 1) {
            mint w = 1;
            for (int s = 0, k = 0; s < n; s += (m << 1)) {
                for (int i = s, j = s + m; i < s + m; i++, j++) {
                    mint a = f[i], b = f[j];
                    f[i] = a + b;
                    f[j] = (a - b) * w;
                }
                w *= dw_inv[__builtin_ctz(++k)];
            }
        }
        mint invn = mint(n).inv();
        for (mint &x : f) x *= invn;
    }
private:
    Number_Theoretic_Transform() = default;
    static void init() {
        if (log > 0) return;
        int mod = mint::mod();
        int tmp = mod - 1;
        log = 1;
        while (tmp % 2 == 0) {
            tmp >>= 1;
            log++;
        }
        dw.resize(log);
        dw_inv.resize(log);
        for (int i = 0; i < log; i++) {
            dw[i] = -root.pow((mod - 1) >> (i + 2));
            dw_inv[i] = dw[i].inv();
        }
    }
};
template<typename mint>
vector<mint>Number_Theoretic_Transform<mint>::dw = vector<mint>();
template<typename mint>
vector<mint>Number_Theoretic_Transform<mint>::dw_inv = vector<mint>();
template<typename mint>
int Number_Theoretic_Transform<mint>::log = 0;
template<typename mint>
mint Number_Theoretic_Transform<mint>::root = 3;

template<typename mint>
struct Formal_Power_Series : vector<mint> {
    using FPS = Formal_Power_Series;
    using vector<mint>::vector;
    using NTT = Number_Theoretic_Transform<mint>;
    void ntt() { NTT::ntt(*this); }
    void intt() { NTT::intt(*this); }
    FPS &operator+=(const mint &r) {
        if (this->empty()) this->resize(1);
        (*this)[0] += r;
        return *this;
    }
    FPS &operator-=(const mint &r) {
        if (this->empty()) this->resize(1);
        (*this)[0] -= r;
        return *this;
    }
    FPS &operator*=(const mint &r) {
        for (mint &x : *this) x *= r;
        return *this;
    }
    FPS &operator/=(const mint &r) {
        mint invr = r.inv();
        return *this *= invr;
    }
    FPS &operator+=(const FPS &f) {
        int n = this->size(), m = f.size();
        if (n < m) this->resize(m);
        for (int i = 0; i < m; i++) (*this)[i] += f[i];
        return *this;
    }
    FPS &operator-=(const FPS &f) {
        int n = this->size(), m = f.size();
        if (n < m) this->resize(m);
        for (int i = 0; i < m; i++) (*this)[i] -= f[i];
        return *this;
    }
    FPS &operator*=(const FPS &f) {
        *this = convolution(*this, f);
        return *this;
    }
    FPS &operator/=(const FPS &f) {
        return *this *= f.inv();
    }
    FPS &operator%=(const FPS &f) {
        *this -= div(f) * f;
        this->shrink();
        return *this;
    }
    FPS div(const FPS &f) const {
        if (this->size() < f.size()) return FPS{};
        int n = this->size() - f.size() + 1;
        return (rev().pre(n) * f.rev().inv(n)).pre(n).rev(n);
    }
    FPS operator+(const mint &r) const { return FPS(*this) += r; }
    FPS operator-(const mint &r) const { return FPS(*this) -= r; }
    FPS operator*(const mint &r) const { return FPS(*this) *= r; }
    FPS operator/(const mint &r) const { return FPS(*this) /= r; }
    FPS operator+(const FPS &f) const { return FPS(*this) += f; }
    FPS operator-(const FPS &f) const { return FPS(*this) -= f; }
    FPS operator*(const FPS &f) const { return FPS(*this) *= f; }
    FPS operator/(const FPS &f) const { return FPS(*this) /= f; }
    FPS operator%(const FPS &f) const { return FPS(*this) %= f; }

    FPS operator-() const {
        return FPS{} - *this;
    }
    FPS operator<<(int n) const {
        FPS res(*this);
        res.insert(res.begin(), n, mint());
        return res;
    }
    FPS operator>>(int n) const {
        if (int(this->size()) <= n) return FPS{};
        FPS res(*this);
        res.erase(res.begin(), res.begin() + n);
        return res;
    }
    FPS &operator<<=(int n) {
        return *this = (*this) << n;
    }
    FPS &operator>>=(int n) {
        return *this = (*this) >> n;
    }
    FPS pre(int n) const {
        n = min(n, int(this->size()));
        return FPS(this->begin(), this->begin() + n);
    }
    FPS rev(int deg = -1) const {
        FPS res(*this);
        if (deg != -1) res.resize(deg, 0);
        reverse(res.begin(), res.end());
        return res;
    }
    FPS dot(const FPS &f) const {
        int n = min(this->size(), f.size());
        FPS res(n);
        for (int i = 0; i < n; i++) res[i] = (*this)[i] * f[i];
        return res;
    }
    void shrink() {
        while (this->size() && this->back() == 0) {
            this->pop_back();
        }
    }
    mint operator()(const mint &x) const {
        mint res = 0, powx = 1;
        for (const mint &a : *this) {
            res += a * powx;
            powx *= x;
        }
        return res;
    }
    FPS diff() const {
        int n = this->size();
        if (n == 0) return FPS{};
        FPS res(n - 1);
        for (int i = 1; i < n; i++) {
            res[i - 1] = i * (*this)[i];
        }
        return res;
    }
    FPS integral() const {
        int n = this->size();
        FPS res(n + 1);
        res[0] = 0;
        for (int i = 0; i < n; i++) {
            res[i + 1] = (*this)[i] / (i + 1);
        }
        return res;
    }
    FPS inv(int deg = -1) const {
        int n = this->size();
        assert(n > 0);
        mint c = (*this)[0];
        assert(c != 0);
        if (deg == -1) deg = n;
        FPS res(deg);
        res[0] = c.inv();
        for (int d = 1; d < deg; d <<= 1) {
            FPS f(d << 1), g(d << 1);
            for (int i = 0; i < n && i < d << 1; i++) f[i] = (*this)[i];
            for (int i = 0; i < d; i++) g[i] = res[i];
            f.ntt();
            g.ntt();
            f = f.dot(g);
            f.intt();
            for (int i = 0; i < d; i++) f[i] = 0;
            f.ntt();
            f = f.dot(g);
            f.intt();
            for (int i = d; i < deg && i < d << 1; i++) res[i] -= f[i];
        }
        return res;
    }
    FPS exp(int deg = -1) const {
        int n = this->size();
        if (deg == -1) deg = n;
        if (n == 0) {
            FPS res(deg);
            res[0] = 1;
            return res;
        }
        assert((*this)[0] == 0);
        vector<mint>inv;
        inv.reserve(deg + 1);
        inv.push_back(0);
        inv.push_back(1);
        auto inplace_diff = [](FPS &f) -> void {
            if (f.empty()) return;
            f.erase(f.begin());
            for (int i = 0; i < int(f.size()); i++) f[i] *= i + 1;
        };
        auto inplace_integral = [&](FPS &f) -> void {
            constexpr int mod = mint::mod();
            while (inv.size() <= f.size()) {
                int i = inv.size();
                inv.push_back(-inv[mod % i] * (mod / i));
            }
            f.insert(f.begin(), 0);
            for (int i = 1; i < int(f.size()); i++) f[i] *= inv[i];
        };
        FPS b = {1, 1 < n ? (*this)[1] : 0};
        FPS c = {1}, z1, z2 = {1, 1};
        for (int d = 2; d < deg; d <<= 1) {
            FPS  y = b;
            y.resize(d << 1);
            y.ntt();
            z1 = z2;
            FPS z = y.dot(z1);
            z.intt();
            fill(z.begin(), z.begin() + (d >> 1), 0);
            z.ntt();
            z = z.dot(-z1);
            z.intt();
            c.insert(c.end(), z.begin() + (d >> 1), z.end());
            z2 = c;
            z2.resize(d << 1);
            z2.ntt();
            FPS x(this->begin(), this->begin() + min(n, d));
            inplace_diff(x);
            x.push_back(0);
            x.ntt();
            x = x.dot(y);
            x.intt();
            x -= b.diff();
            x.resize(d << 1);
            for (int i = 0; i < d - 1; i++) x[i + d] = x[i], x[i] = 0;
            x.ntt();
            x = x.dot(z2);
            x.intt();
            x.pop_back();
            inplace_integral(x);
            for (int i = d; i < min(n, d << 1); i++) x[i] += (*this)[i];
            fill(x.begin(), x.begin() + d, 0);
            x.ntt();
            x = x.dot(y);
            x.intt();
            b.insert(b.end(), x.begin() + d, x.end());
        }
        return FPS(b.begin(), b.begin() + deg);
    }
    FPS log(int deg = -1) const {
        assert((*this)[0] == 1);
        if (deg == -1) deg = this->size();
        return (diff() * inv()).integral().pre(deg);
    }
    FPS pow(long long k, int deg = -1) const {
        if (deg == -1) deg = this->size();
        if (k == 0) {
            FPS res(deg);
            res[0] = 1;
            return res;
        }
        FPS res(*this);
        int p = 0;
        while (p < int(res.size()) && res[p] == 0) p++;
        if (p > (deg - 1) / k) return FPS(deg);
        res >>= p;
        deg -= p * k;
        mint c = res[0];
        res = ((res / c).log(deg) * k).exp(deg) * c.pow(k);
        res <<= p * k;
        return res;
    }
private:
    FPS convolution(FPS f, FPS g, int deg = -1) {
        int n = f.size(), m = g.size();
        if (n == 0 || m == 0) return FPS{};
        int sz = 1;
        while (sz < n + m - 1) sz <<= 1;
        f.resize(sz);
        f.ntt();
        g.resize(sz);
        g.ntt();
        f = f.dot(g);
        f.intt();
        if (deg == -1) deg = n + m - 1;
        f.resize(deg);
        return f;
    }
};

template<typename mint>
vector<mint> Berlekamp_Massey(const vector<mint> &a) {
    int n = a.size();
    vector<mint> b = {1}, c = {1};
    b.reserve(n + 1);
    c.reserve(n + 1);
    mint y = 1;
    for (int d = 1; d <= n; d ++) {
        int k = b.size(), l = c.size();
        mint x = 0;
        for (int i = 0; i < l; i++) {
            x += c[i] * a[d - l + i];
        }
        b.push_back(0);
        k++;
        if (x == 0) continue;
        mint buf = x / y;
        if (l < k) {
            auto tmp = c;
            c.insert(c.begin(), k - l, 0);
            for (int i = 0; i < k; i++) {
                c[k - 1 - i] -= buf * b[k - 1 - i];
            }
            b = tmp;
            y = x;
        } else {
            for (int i = 0; i < k; i++) {
                c[l - 1 - i] -= buf * b[k - 1 - i];
            }
        }
    }
    reverse(c.begin(), c.end());
    for (mint &x : c) x = -x;
    return c;
}

template<typename FPS>
typename FPS::value_type Bostan_Mori(FPS p, FPS q, long long k) {
    using mint = typename FPS::value_type;
    mint res = 0;
    if (p.size() >= q.size()) {
        FPS r = p.div(q);
        p -= q * r;
        p.shrink();
        if (k < int(r.size())) res += r[k];
    }
    if (p.empty()) return res;
    p.resize(q.size() - 1);
    auto sub = [&](const FPS &f, bool odd = false) -> FPS {
        int n = (f.size() + !odd) / 2;
        FPS g(n);
        for (int i = odd; i < int(f.size()); i += 2) {
            g[i / 2] = f[i];
        }
        return g;
    };
    while (k > 0) {
        FPS q2 = q;
        for (int i = 1; i < int(q2.size()); i += 2) {
            q2[i] = -q2[i];
        }
        p = sub(p * q2, k & 1);
        q = sub(q * q2);
        k >>= 1;
    }
    return res + p[0];
}

template<typename FPS>
typename FPS::value_type linear_recurrence(FPS a, FPS c, long long k) {
    assert(a.size() == c.size());
    c = FPS{1} - (c << 1);
    return Bostan_Mori((a * c).pre(a.size()), c, k);
}

template<typename mint>
mint BMBM(const vector<mint> &x, long long k) {
    using FPS = Formal_Power_Series<mint>;
    auto tmp = Berlekamp_Massey(x);
    int n = tmp.size() - 1;
    FPS a(n), c(n);
    for (int i = 0; i < n; i++) {
        a[i] = x[i];
        c[i] = tmp[i + 1];
    }
    return linear_recurrence(a, c, k);
}

using mint = modint<998244353>;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    long long N;
    int K;
    cin >> N >> K;
    int M = 4 * K;
    vector<mint> F(M), A(M);
    F[0] = F[1] = 1;
    A[0] = 1, A[1] = 2;
    for (int i = 2; i < M; i++) {
        F[i] = F[i - 1] + F[i - 2];
        A[i] = A[i - 1] + F[i].pow(K);
    }
    cout << BMBM(A, N - 1).val() << endl;
}