hly1204's library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub hly1204/library

:heavy_check_mark: test/set_power_series/polynomial_composite_set_power_series.0.test.cpp

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/polynomial_composite_set_power_series"

#include "modint.hpp"
#include "sps_in_poly.hpp"
#include <iostream>
#include <vector>

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    using mint = ModInt<998244353>;
    int n, m;
    std::cin >> n >> m;
    std::vector<mint> F(n), G(1 << m);
    for (int i = 0; i < n; ++i) std::cin >> F[i];
    for (int i = 0; i < (1 << m); ++i) std::cin >> G[i];
    const auto FG = sps_in_poly(F, G);
    for (int i = 0; i < (1 << m); ++i) std::cout << FG[i] << ' ';
    return 0;
}
#line 1 "test/set_power_series/polynomial_composite_set_power_series.0.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/polynomial_composite_set_power_series"

#line 2 "modint.hpp"

#include <iostream>
#include <type_traits>

// clang-format off
template<unsigned Mod> class ModInt {
    static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)");
    template<typename Int>
    static std::enable_if_t<std::is_integral_v<Int>, unsigned> safe_mod(Int v) { using D = std::common_type_t<Int, unsigned>; return (v %= (int)Mod) < 0 ? (D)(v + (int)Mod) : (D)v; }
    struct PrivateConstructor {} static inline private_constructor{};
    ModInt(PrivateConstructor, unsigned v) : v_(v) {}
    unsigned v_;

public:
    static unsigned mod() { return Mod; }
    static ModInt from_raw(unsigned v) { return ModInt(private_constructor, v); }
    static ModInt zero() { return from_raw(0); }
    static ModInt one() { return from_raw(1); }
    bool is_zero() const { return v_ == 0; }
    bool is_one() const { return v_ == 1; }
    ModInt() : v_() {}
    template<typename Int, typename std::enable_if_t<std::is_signed_v<Int>, int> = 0> ModInt(Int v) : v_(safe_mod(v)) {}
    template<typename Int, typename std::enable_if_t<std::is_unsigned_v<Int>, int> = 0> ModInt(Int v) : v_(v % Mod) {}
    unsigned val() const { return v_; }
    ModInt operator-() const { return from_raw(v_ == 0 ? v_ : Mod - v_); }
    ModInt pow(long long e) const { if (e < 0) return inv().pow(-e); for (ModInt x(*this), res(from_raw(1));; x *= x) { if (e & 1) res *= x; if ((e >>= 1) == 0) return res; }}
    ModInt inv() const { int x1 = 1, x3 = 0, a = val(), b = Mod; while (b) { const int q = a / b, x1_old = x1, a_old = a; x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q; } return from_raw(x1 < 0 ? x1 + (int)Mod : x1); }
    template<bool Odd = (Mod & 1)> std::enable_if_t<Odd, ModInt> div_by_2() const { if (v_ & 1) return from_raw((v_ + Mod) >> 1); return from_raw(v_ >> 1); }
    ModInt &operator+=(const ModInt &a) { if ((v_ += a.v_) >= Mod) v_ -= Mod; return *this; }
    ModInt &operator-=(const ModInt &a) { if ((v_ += Mod - a.v_) >= Mod) v_ -= Mod; return *this; }
    ModInt &operator*=(const ModInt &a) { v_ = (unsigned long long)v_ * a.v_ % Mod; return *this; }
    ModInt &operator/=(const ModInt &a) { return *this *= a.inv(); }
    ModInt &operator++() { return *this += one(); }
    ModInt operator++(int) { ModInt o(*this); *this += one(); return o; }
    ModInt &operator--() { return *this -= one(); }
    ModInt operator--(int) { ModInt o(*this); *this -= one(); return o; }
    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.v_ == b.v_; }
    friend bool operator!=(const ModInt &a, const ModInt &b) { return a.v_ != b.v_; }
    friend std::istream &operator>>(std::istream &a, ModInt &b) { int v; a >> v; b.v_ = safe_mod(v); return a; }
    friend std::ostream &operator<<(std::ostream &a, const ModInt &b) { return a << b.val(); }
};
// clang-format on
#line 2 "sps_in_poly.hpp"

#line 2 "subset_conv.hpp"

#include <cassert>
#include <vector>

template<typename Tp> inline std::vector<std::vector<Tp>> to_ranked(const std::vector<Tp> &A) {
    const int N = A.size();
    int LogN    = 0;
    while ((1 << LogN) != N) ++LogN;
    std::vector res(LogN + 1, std::vector<Tp>(N));
    for (int i = 0; i < N; ++i) res[__builtin_popcount(i)][i] = A[i];
    return res;
}

template<typename Tp> inline std::vector<Tp> from_ranked(const std::vector<std::vector<Tp>> &A) {
    const int N = A[0].size();
    std::vector<Tp> res(N);
    for (int i = 0; i < N; ++i) res[i] = A[__builtin_popcount(i)][i];
    return res;
}

template<typename Iterator> inline void subset_zeta_n(Iterator a, int n) {
    assert((n & (n - 1)) == 0);
    for (int i = 2; i <= n; i *= 2)
        for (int j = 0; j < n; j += i)
            for (int k = j; k < j + i / 2; ++k) a[k + i / 2] += a[k];
}

template<typename Tp> inline void subset_zeta(std::vector<Tp> &a) {
    subset_zeta_n(a.begin(), a.size());
}

template<typename Iterator> inline void subset_moebius_n(Iterator a, int n) {
    assert((n & (n - 1)) == 0);
    for (int i = 2; i <= n; i *= 2)
        for (int j = 0; j < n; j += i)
            for (int k = j; k < j + i / 2; ++k) a[k + i / 2] -= a[k];
}

template<typename Tp> inline void subset_moebius(std::vector<Tp> &a) {
    subset_moebius_n(a.begin(), a.size());
}

template<typename Tp>
inline std::vector<Tp> subset_convolution(const std::vector<Tp> &A, const std::vector<Tp> &B) {
    assert(A.size() == B.size());
    const int N = A.size();
    int LogN    = 0;
    while ((1 << LogN) != N) ++LogN;
    auto rankedA = to_ranked(A);
    auto rankedB = to_ranked(B);

    for (int i = 0; i <= LogN; ++i) {
        subset_zeta(rankedA[i]);
        subset_zeta(rankedB[i]);
    }

    // see: https://codeforces.com/blog/entry/126418
    // see: https://oeis.org/A025480
    std::vector<int> map(LogN + 1);
    for (int i = 0; i <= LogN; ++i) map[i] = (i & 1) ? map[i / 2] : i / 2;

    std::vector rankedAB(LogN / 2 + 1, std::vector<Tp>(N));
    for (int i = 0; i <= LogN; ++i)
        for (int j = 0; i + j <= LogN; ++j)
            for (int k = (1 << j) - 1; k < N; ++k)
                rankedAB[map[i + j]][k] += rankedA[i][k] * rankedB[j][k];

    for (int i = 0; i <= LogN / 2; ++i) subset_moebius(rankedAB[i]);

    std::vector<Tp> res(N);
    for (int i = 0; i < N; ++i) res[i] = rankedAB[map[__builtin_popcount(i)]][i];
    return res;
}
#line 4 "sps_in_poly.hpp"
#include <algorithm>
#line 7 "sps_in_poly.hpp"

// returns F(G)
// requires deg(F)<=n, G(0)=0
// see:
// [1]: Elegia. Optimal Algorithm on Polynomial Composite Set Power Series.
//      https://codeforces.com/blog/entry/92183
template<typename Tp>
inline std::vector<Tp> sps_in_egf(const std::vector<Tp> &F, const std::vector<Tp> &G) {
    const int N = (int)F.size() - 1;
    assert((int)G.size() == (1 << N));
    assert(G[0] == 0);

    auto conv_ranked = [](const auto &rankedA, const auto &rankedB, int LogN, auto res) {
        const int N = (1 << LogN);
        std::vector<int> map(LogN + 1);
        for (int i = 0; i <= LogN; ++i) map[i] = (i & 1) ? map[i / 2] : i / 2;
        std::vector rankedAB(LogN / 2 + 1, std::vector<Tp>(N));
        for (int i = 0; i <= LogN; ++i)
            for (int j = 0; i + j <= LogN; ++j)
                for (int k = (1 << j) - 1; k < N; ++k)
                    rankedAB[map[i + j]][k] += rankedA[i][k] * rankedB[j][k];
        for (int i = 0; i <= LogN / 2; ++i) subset_moebius(rankedAB[i]);
        for (int i = 0; i < N; ++i) res[i] = rankedAB[map[__builtin_popcount(i)]][i];
    };

    std::vector<std::vector<std::vector<Tp>>> rankedG;
    std::vector res = {F[N]};
    for (int i = 0; i < N; ++i) {
        auto &&rankedGi = rankedG.emplace_back(
            to_ranked(std::vector(G.begin() + (1 << i), G.begin() + (2 << i))));
        auto rankedRes = to_ranked(res);
        for (int j = 0; j <= i; ++j) {
            subset_zeta(rankedGi[j]);
            subset_zeta(rankedRes[j]);
        }
        res.resize(1 << (i + 1));
        res[0] = F[N - (i + 1)];
        for (int j = 0; j <= i; ++j) conv_ranked(rankedG[j], rankedRes, j, res.begin() + (1 << j));
    }

    return res;
}

template<typename Tp> inline std::vector<Tp> sps_in_poly(std::vector<Tp> F, std::vector<Tp> G) {
    const int N = G.size();
    int LogN    = 0;
    while ((1 << LogN) != N) ++LogN;

    if (G[0] != 0) {
        std::vector<Tp> bin(LogN + 1), pw(F.size() + 1), FF(LogN + 1);
        pw[0] = 1;
        for (int i = 1; i < (int)pw.size(); ++i) pw[i] = pw[i - 1] * G[0];
        G[0]   = 0;
        bin[0] = 1;
        for (int i = 0; i < (int)F.size(); ++i) {
            for (int j = 0; j <= std::min(LogN, i); ++j) FF[j] += F[i] * bin[j] * pw[i - j];
            for (int j = LogN; j > 0; --j) bin[j] += bin[j - 1];
        }
        FF.swap(F);
    }

    F.resize(LogN + 1);
    Tp c = 1;                                       // factorial
    for (int i = 1; i <= LogN; ++i) F[i] *= c *= i; // to EGF
    return sps_in_egf(F, G);
}
#line 7 "test/set_power_series/polynomial_composite_set_power_series.0.test.cpp"

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    using mint = ModInt<998244353>;
    int n, m;
    std::cin >> n >> m;
    std::vector<mint> F(n), G(1 << m);
    for (int i = 0; i < n; ++i) std::cin >> F[i];
    for (int i = 0; i < (1 << m); ++i) std::cin >> G[i];
    const auto FG = sps_in_poly(F, G);
    for (int i = 0; i < (1 << m); ++i) std::cout << FG[i] << ' ';
    return 0;
}
Back to top page