Subset Convolution
(math/subset_convolution.hpp)

• View this file on GitHub
• Last update: 2026-07-13 05:19:17+09:00
• Include: #include "math/subset_convolution.hpp"

Overview

Subset convolution combines two arrays indexed by bit masks. For every mask $S$, it computes

\[c_S = \sum_{T \subseteq S} a_T b_{S \setminus T}.\]

Equivalently, every term partitions S into two disjoint masks. This operation appears in subset dynamic programming when two independently solved parts must cover the current set exactly.

#include "math/subset_convolution.hpp"

The function is in m1une::math.

Interface

template <typename T>
std::vector<T> subset_convolution(
    std::vector<T> first,
    std::vector<T> second
);

Function	Description	Complexity
subset_convolution(first, second)	Returns the subset convolution of two equally sized mask arrays.	$O(K^2 2^K)$ time and $O(K 2^K)$ additional memory

Both input vectors must have the same size. Their size must be zero or a power of two; a nonempty vector of size $2^K$ represents masks on K bits. Empty inputs produce an empty result. The result has the same size as each input.

The arguments are passed by value and used as workspace, so the caller’s vectors are unchanged unless explicitly moved into the function. The element type T must be default-constructible and support multiplication, operator+=, and operator-=. Modular integers and ordinary integer types satisfy these requirements.

The implementation separates coefficients by mask popcount, applies ranked subset-zeta transforms, multiplies the rank polynomials pointwise, and applies the ranked Mobius transform. It uses two flat ranked work buffers rather than allocating a separate vector for every mask.

Example

#include "math/subset_convolution.hpp"

#include <iostream>
#include <vector>

int main() {
    std::vector<long long> a{1, 2, 3, 4};
    std::vector<long long> b{5, 6, 7, 8};

    std::vector<long long> c =
        m1une::math::subset_convolution(a, b);

    // c[3] = a[0] * b[3] + a[1] * b[2]
    //      + a[2] * b[1] + a[3] * b[0] = 60.
    std::cout << c[3] << '\n';
}

Required by

• Math All (math/all.hpp)
• Set Power Series (math/set_power_series.hpp)

Verified with

• verify/math/math_algorithms.test.cpp
• verify/math/set_power_series_exp.test.cpp
• verify/math/set_power_series_log.test.cpp
• verify/math/subset_convolution.test.cpp

Code

#ifndef M1UNE_MATH_SUBSET_CONVOLUTION_HPP
#define M1UNE_MATH_SUBSET_CONVOLUTION_HPP 1

#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

namespace m1une {
namespace math {

template <typename T>
std::vector<T> subset_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    assert(first.size() == second.size());
    if (first.empty()) return {};
    assert((first.size() & (first.size() - 1)) == 0);

    const std::size_t size = first.size();
    std::size_t bit_count = 0;
    while ((std::size_t(1) << bit_count) < size) ++bit_count;
    const std::size_t rank_count = bit_count + 1;

    std::vector<T> first_ranked(size * rank_count);
    std::vector<T> second_ranked(size * rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        const std::size_t rank = std::popcount(mask);
        first_ranked[mask * rank_count + rank] = std::move(first[mask]);
        second_ranked[mask * rank_count + rank] = std::move(second[mask]);
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] +=
                    first_ranked[source + rank];
                second_ranked[destination + rank] +=
                    second_ranked[source + rank];
            }
        }
    }

    std::vector<T> product(rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        for (T& value : product) value = T{};
        const std::size_t offset = mask * rank_count;
        const std::size_t rank_limit = std::popcount(mask);
        for (std::size_t left = 0; left <= rank_limit; ++left) {
            const std::size_t right_limit =
                std::min(rank_limit, bit_count - left);
            for (std::size_t right = 0; right <= right_limit; ++right) {
                product[left + right] +=
                    first_ranked[offset + left] *
                    second_ranked[offset + right];
            }
        }
        for (std::size_t rank = 0; rank < rank_count; ++rank) {
            first_ranked[offset + rank] = std::move(product[rank]);
        }
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] -=
                    first_ranked[source + rank];
            }
        }
    }

    std::vector<T> result(size);
    for (std::size_t mask = 0; mask < size; ++mask) {
        result[mask] = std::move(
            first_ranked[mask * rank_count + std::popcount(mask)]
        );
    }
    return result;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_SUBSET_CONVOLUTION_HPP

#line 1 "math/subset_convolution.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

namespace m1une {
namespace math {

template <typename T>
std::vector<T> subset_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    assert(first.size() == second.size());
    if (first.empty()) return {};
    assert((first.size() & (first.size() - 1)) == 0);

    const std::size_t size = first.size();
    std::size_t bit_count = 0;
    while ((std::size_t(1) << bit_count) < size) ++bit_count;
    const std::size_t rank_count = bit_count + 1;

    std::vector<T> first_ranked(size * rank_count);
    std::vector<T> second_ranked(size * rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        const std::size_t rank = std::popcount(mask);
        first_ranked[mask * rank_count + rank] = std::move(first[mask]);
        second_ranked[mask * rank_count + rank] = std::move(second[mask]);
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] +=
                    first_ranked[source + rank];
                second_ranked[destination + rank] +=
                    second_ranked[source + rank];
            }
        }
    }

    std::vector<T> product(rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        for (T& value : product) value = T{};
        const std::size_t offset = mask * rank_count;
        const std::size_t rank_limit = std::popcount(mask);
        for (std::size_t left = 0; left <= rank_limit; ++left) {
            const std::size_t right_limit =
                std::min(rank_limit, bit_count - left);
            for (std::size_t right = 0; right <= right_limit; ++right) {
                product[left + right] +=
                    first_ranked[offset + left] *
                    second_ranked[offset + right];
            }
        }
        for (std::size_t rank = 0; rank < rank_count; ++rank) {
            first_ranked[offset + rank] = std::move(product[rank]);
        }
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] -=
                    first_ranked[source + rank];
            }
        }
    }

    std::vector<T> result(size);
    for (std::size_t mask = 0; mask < size; ++mask) {
        result[mask] = std::move(
            first_ranked[mask * rank_count + std::popcount(mask)]
        );
    }
    return result;
}

}  // namespace math
}  // namespace m1une

Back to top page