Subset Convolution
(src/convolution/subset_convolution.hpp)

• View this file on GitHub
• Last update: 2024-06-24 04:02:51+09:00
• Include: #include "src/convolution/subset_convolution.hpp"

入力

$\lbrack n \rbrack := \lbrace 0, 1, \dots , n - 1 \rbrace$，$\mathbb{K}$ を環とする． このもとでの集合関数 $f \colon 2 ^ {\lbrack n \rbrack} \to \mathbb{K},\ g \colon 2 ^ {\lbrack n \rbrack} \to \mathbb{K}$．

出力

関数 $f, g$ の subset convolution． すなわち，

\[(f \ast g) \colon 2 ^ {\lbrack n \rbrack} \ni S \to \sum_{T \subseteq S} f(T) g(S \setminus T) \in \mathbb{K}\]

で定義される集合関数 $(f \ast g)$．

計算量

時間計算量 $\mathrm{O}(n ^ 2 2 ^ n)$

Links

• Xmas Contest 2020 H: Hierarchical Phylogeny 解説 - hos.lyric’s blog
• Subset Convolutionのアルゴリズム – 37zigenのHP
• Optimal Algorithm on Polynomial Composite Set Power Series - Codeforces
• 集合冪級数 (Set Power Series) のダブリングまとめ - suisen のブログ
• 集合べき級数解説 | maspyのHP

出題例

Required by

• Exp of Set Power Series (src/math/exp_of_set_power_series.hpp)

Verified with

• test/yosupo/exp_of_set_power_series.test.cpp
• test/yosupo/subset_convolution.test.cpp

Code

#pragma once
#include <array>
#include <cassert>
#include <vector>

namespace internal {

template <typename T, int LIM> std::vector<std::array<T, LIM + 1>> ranked_zeta(const std::vector<T>& a) {
    int n = a.size(), len = __builtin_ctz(n);
    assert((n & (n - 1)) == 0);
    std::vector<std::array<T, LIM + 1>> res(n);
    for (int i = 0; i < n; i++) res[i][__builtin_popcount(i)] = a[i];
    for (int step = 1; step < n; step <<= 1) {
        for (int start = 0; start < n; start += 2 * step) {
            for (int i = start; i < start + step; i++) {
                for (int j = 0; j <= len; j++) {
                    res[i | step][j] += res[i][j];
                }
            }
        }
    }
    return res;
}

template <typename T, int LIM> std::vector<T> ranked_mobius(std::vector<std::array<T, LIM + 1>>& ranked) {
    int n = ranked.size(), len = __builtin_ctz(n);
    assert((n & (n - 1)) == 0);
    for (int step = 1; step < n; step <<= 1) {
        for (int start = 0; start < n; start += 2 * step) {
            for (int i = start; i < start + step; i++) {
                for (int j = 0; j <= len; j++) {
                    ranked[i | step][j] -= ranked[i][j];
                }
            }
        }
    }
    std::vector<T> res(n);
    for (int i = 0; i < n; i++) res[i] = ranked[i][__builtin_popcount(i)];
    return res;
}

}  // namespace internal

template <typename T, int LIM = 20>
std::vector<T> subset_convolution(const std::vector<T>& a, const std::vector<T>& b) {
    auto ra = internal::ranked_zeta<T, LIM>(a);
    auto rb = internal::ranked_zeta<T, LIM>(b);
    int n = ra.size(), len = __builtin_ctz(n);
    for (int i = 0; i < n; i++) {
        auto &f = ra[i], &g = rb[i];
        for (int j = len; j >= 0; j--) {
            T sum = 0;
            for (int k = 0; k <= j; k++) sum += f[k] * g[j - k];
            f[j] = sum;
        }
    }
    return internal::ranked_mobius<T, LIM>(ra);
}

#line 2 "src/convolution/subset_convolution.hpp"
#include <array>
#include <cassert>
#include <vector>

namespace internal {

template <typename T, int LIM> std::vector<std::array<T, LIM + 1>> ranked_zeta(const std::vector<T>& a) {
    int n = a.size(), len = __builtin_ctz(n);
    assert((n & (n - 1)) == 0);
    std::vector<std::array<T, LIM + 1>> res(n);
    for (int i = 0; i < n; i++) res[i][__builtin_popcount(i)] = a[i];
    for (int step = 1; step < n; step <<= 1) {
        for (int start = 0; start < n; start += 2 * step) {
            for (int i = start; i < start + step; i++) {
                for (int j = 0; j <= len; j++) {
                    res[i | step][j] += res[i][j];
                }
            }
        }
    }
    return res;
}

template <typename T, int LIM> std::vector<T> ranked_mobius(std::vector<std::array<T, LIM + 1>>& ranked) {
    int n = ranked.size(), len = __builtin_ctz(n);
    assert((n & (n - 1)) == 0);
    for (int step = 1; step < n; step <<= 1) {
        for (int start = 0; start < n; start += 2 * step) {
            for (int i = start; i < start + step; i++) {
                for (int j = 0; j <= len; j++) {
                    ranked[i | step][j] -= ranked[i][j];
                }
            }
        }
    }
    std::vector<T> res(n);
    for (int i = 0; i < n; i++) res[i] = ranked[i][__builtin_popcount(i)];
    return res;
}

}  // namespace internal

template <typename T, int LIM = 20>
std::vector<T> subset_convolution(const std::vector<T>& a, const std::vector<T>& b) {
    auto ra = internal::ranked_zeta<T, LIM>(a);
    auto rb = internal::ranked_zeta<T, LIM>(b);
    int n = ra.size(), len = __builtin_ctz(n);
    for (int i = 0; i < n; i++) {
        auto &f = ra[i], &g = rb[i];
        for (int j = len; j >= 0; j--) {
            T sum = 0;
            for (int k = 0; k <= j; k++) sum += f[k] * g[j - k];
            f[j] = sum;
        }
    }
    return internal::ranked_mobius<T, LIM>(ra);
}

Back to top page