Subset convolution (tools/subset_convolution.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2026-02-21 18:33:47+09:00
• Include: #include "tools/subset_convolution.hpp"

(1)
template <std::ranges::input_range R1, std::ranges::input_range R2>
std::vector<std::common_type_t<std::ranges::range_value_t<R1>, std::ranges::range_value_t<R2>>> subset_convolution(R1&& a, R2&& b);

(2)
template <tools::ring R, std::ranges::input_range R1, std::ranges::input_range R2>
requires std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R1>>
      && std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R2>>
std::vector<typename R::add::T> subset_convolution(R1&& a, R2&& b);

Given size $2^N$ sequences $(a_0, a_1, \ldots, a_{2^N - 1})$ and $(b_0, b_1, \ldots, b_{2^N - 1})$, it returns the sequence $(c_0, c_1, \ldots, c_{2^N - 1})$ where $\displaystyle c_k = \sum_{\substack{(i~\mathrm{AND}~j) = 0\\ (i~\mathrm{OR}~j) = k}} a_i b_j$. In (2), addition and multiplication are defined by R.

Constraints

• $|a|$ is a power of $2$.
• $|b|$ is a power of $2$.
• $|a| = |b|$

Time Complexity

• $O(2^N N^2)$

License

• CC0

Author

• anqooqie

Depends on

• Concept for arithmetic types including tools::int128_t and tools::uint128_t (tools/arithmetic.hpp)
• Concept for commutative group (tools/commutative_group.hpp)
• Concept for commutative monoid (tools/commutative_monoid.hpp)
• std::gcd(m, n) extended for my library (tools/gcd.hpp)
• Concept for group (tools/group.hpp)
• Typical groups (tools/groups.hpp)
• Concept for integral types including tools::int128_t and tools::uint128_t (tools/integral.hpp)
• std::is_integral extended for tools::int128_t and tools::uint128_t (tools/is_integral.hpp)
• Concept for monoid (tools/monoid.hpp)
• Typical monoids (tools/monoids.hpp)
• tools::integral except for bool (tools/non_bool_integral.hpp)
• Concept for ring (tools/ring.hpp)
• Typical rings (tools/rings.hpp)
• Concept for semiring (tools/semiring.hpp)
• Typical semirings (tools/semirings.hpp)

Verified with

• tests/subset_convolution.test.cpp

Code

#ifndef TOOLS_SUBSET_CONVOLUTION_HPP
#define TOOLS_SUBSET_CONVOLUTION_HPP

#include <algorithm>
#include <bit>
#include <cassert>
#include <concepts>
#include <iterator>
#include <ranges>
#include <type_traits>
#include <utility>
#include <vector>
#include "tools/ring.hpp"
#include "tools/rings.hpp"

namespace tools {
  template <tools::ring R, std::ranges::input_range R1, std::ranges::input_range R2>
  requires std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R1>>
        && std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R2>>
  auto subset_convolution(R1&& a, R2&& b) {
    if constexpr (std::ranges::sized_range<R1> && std::ranges::sized_range<R2>) {

      using Add = typename R::add;
      using Mul = typename R::mul;
      using T = typename Add::T;

      assert(std::has_single_bit(std::ranges::size(a)));
      assert(std::has_single_bit(std::ranges::size(b)));
      assert(std::ranges::size(a) == std::ranges::size(b));

      const int N = std::countr_zero(std::ranges::size(a));
      const int pow2_N = 1 << N;

      auto za = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      {
        int i = 0;
        for (auto&& a_i : a) {
          za[i][std::popcount<unsigned int>(i)] = a_i;
          ++i;
        }
      }

      auto zb = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      {
        int i = 0;
        for (auto&& b_i : b) {
          zb[i][std::popcount<unsigned int>(i)] = b_i;
          ++i;
        }
      }

      for (int w = 0; w < N; ++w) {
        for (int i = 0; i < pow2_N; ++i) {
          if (i & (1 << w)) {
            for (int j = 0; j <= N; ++j) {
              za[i][j] = Add::op(za[i][j], za[i ^ (1 << w)][j]);
              zb[i][j] = Add::op(zb[i][j], zb[i ^ (1 << w)][j]);
            }
          }
        }
      }

      auto zc = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      for (int i = 0; i < pow2_N; ++i) {
        for (int j = 0; j <= N; ++j) {
          for (int k = 0; k <= j; ++k) {
            zc[i][j] = Add::op(zc[i][j], Mul::op(za[i][k], zb[i][j - k]));
          }
        }
      }

      for (int w = 0; w < N; ++w) {
        for (int i = 0; i < pow2_N; ++i) {
          if (i & (1 << w)) {
            for (int j = 0; j <= N; ++j) {
              zc[i][j] = Add::op(zc[i][j], Add::inv(zc[i ^ (1 << w)][j]));
            }
          }
        }
      }

      std::vector<T> c;
      c.reserve(pow2_N);
      for (int i = 0; i < pow2_N; ++i) {
        c.push_back(zc[i][std::popcount<unsigned int>(i)]);
      }

      return c;

    } else {
      return tools::subset_convolution<R>(
        std::forward<R1>(a) | std::ranges::to<std::vector<std::ranges::range_value_t<R1>>>(),
        std::forward<R2>(b) | std::ranges::to<std::vector<std::ranges::range_value_t<R2>>>()
      );
    }
  }

  template <std::ranges::input_range R1, std::ranges::input_range R2>
  auto subset_convolution(R1&& a, R2&& b) {
    using T = std::common_type_t<std::ranges::range_value_t<R1>, std::ranges::range_value_t<R2>>;
    return tools::subset_convolution<tools::rings::plus_multiplies<T>, R1, R2>(std::forward<R1>(a), std::forward<R2>(b));
  }
}

#endif

#line 1 "tools/subset_convolution.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <concepts>
#include <iterator>
#include <ranges>
#include <type_traits>
#include <utility>
#include <vector>
#line 1 "tools/ring.hpp"



#line 1 "tools/commutative_group.hpp"



#line 1 "tools/commutative_monoid.hpp"



#line 1 "tools/monoid.hpp"



#line 5 "tools/monoid.hpp"

namespace tools {
  template <typename M>
  concept monoid = requires(typename M::T x, typename M::T y) {
    { M::op(x, y) } -> std::same_as<typename M::T>;
    { M::e() } -> std::same_as<typename M::T>;
  };
}


#line 5 "tools/commutative_monoid.hpp"

namespace tools {
  template <typename M>
  concept commutative_monoid = tools::monoid<M>;
}


#line 1 "tools/group.hpp"



#line 6 "tools/group.hpp"

namespace tools {
  template <typename G>
  concept group = tools::monoid<G> && requires(typename G::T x) {
    { G::inv(x) } -> std::same_as<typename G::T>;
  };
}


#line 6 "tools/commutative_group.hpp"

namespace tools {
  template <typename G>
  concept commutative_group = tools::group<G> && tools::commutative_monoid<G>;
}


#line 1 "tools/semiring.hpp"



#line 6 "tools/semiring.hpp"

namespace tools {
  template <typename R>
  concept semiring = tools::commutative_monoid<typename R::add> && tools::monoid<typename R::mul> && std::same_as<typename R::add::T, typename R::mul::T>;
}


#line 6 "tools/ring.hpp"

namespace tools {
  template <typename R>
  concept ring = tools::semiring<R> && tools::commutative_group<typename R::add>;
}


#line 1 "tools/rings.hpp"



#line 1 "tools/groups.hpp"



#include <cstddef>
#line 1 "tools/arithmetic.hpp"



#line 1 "tools/integral.hpp"



#line 1 "tools/is_integral.hpp"



#line 5 "tools/is_integral.hpp"

namespace tools {
  template <typename T>
  struct is_integral : std::is_integral<T> {};

  template <typename T>
  inline constexpr bool is_integral_v = tools::is_integral<T>::value;
}


#line 5 "tools/integral.hpp"

namespace tools {
  template <typename T>
  concept integral = tools::is_integral_v<T>;
}


#line 6 "tools/arithmetic.hpp"

namespace tools {
  template <typename T>
  concept arithmetic = tools::integral<T> || std::floating_point<T>;
}


#line 7 "tools/groups.hpp"

namespace tools {
  namespace groups {
    template <typename G>
    struct bit_xor {
      using T = G;
      static T op(const T& x, const T& y) {
        return x ^ y;
      }
      static T e() {
        return T(0);
      }
      static T inv(const T& x) {
        return x;
      }
    };

    template <typename G>
    struct multiplies {
      using T = G;
      static T op(const T& x, const T& y) {
        return x * y;
      }
      static T e() {
        return T(1);
      }
      static T inv(const T& x) {
        return e() / x;
      }
    };

    template <typename G>
    struct plus {
      using T = G;
      static T op(const T& x, const T& y) {
        return x + y;
      }
      static T e() {
        return T(0);
      }
      static T inv(const T& x) {
        return -x;
      }
    };
  }
}


#line 1 "tools/monoids.hpp"



#line 8 "tools/monoids.hpp"
#include <limits>
#line 1 "tools/gcd.hpp"



#include <numeric>
#line 7 "tools/gcd.hpp"

namespace tools {
  namespace detail::gcd {
    template <typename M, typename N>
    struct impl {
      constexpr decltype(auto) operator()(const M m, const N n) const noexcept(noexcept(std::gcd(m, n))) {
        return std::gcd(m, n);
      }
    };
  }

  template <typename M, typename N>
  constexpr decltype(auto) gcd(M&& m, N&& n) noexcept(noexcept(tools::detail::gcd::impl<std::remove_cvref_t<M>, std::remove_cvref_t<N>>{}(std::forward<M>(m), std::forward<N>(n)))) {
    return tools::detail::gcd::impl<std::remove_cvref_t<M>, std::remove_cvref_t<N>>{}(std::forward<M>(m), std::forward<N>(n));
  }
}


#line 1 "tools/non_bool_integral.hpp"



#line 7 "tools/non_bool_integral.hpp"

namespace tools {
  template <typename T>
  concept non_bool_integral = tools::integral<T> && !std::same_as<std::remove_cv_t<T>, bool>;
}


#line 14 "tools/monoids.hpp"

namespace tools {
  namespace monoids {
    template <typename M>
    struct bit_and {
      using T = M;
      static T op(const T& x, const T& y) {
        return x & y;
      }
      static T e() {
        return std::numeric_limits<T>::max();
      }
    };

    template <typename M>
    struct bit_or {
      using T = M;
      static T op(const T& x, const T& y) {
        return x | y;
      }
      static T e() {
        return T(0);
      }
    };

    template <typename M>
    requires requires (M x, M y) {
      {tools::gcd(x, y)} -> std::convertible_to<M>;
    }
    struct gcd {
      using T = M;
      static T op(const T& x, const T& y) {
        return tools::gcd(x, y);
      }
      static T e() {
        return T(0);
      }
    };

    template <typename M, M ...dummy>
    struct max;

    template <tools::arithmetic M>
    struct max<M> {
      using T = M;
      static T op(const T& x, const T& y) {
        return std::max(x, y);
      }
      static T e() {
        if constexpr (tools::integral<M>) {
          return std::numeric_limits<M>::min();
        } else {
          return -std::numeric_limits<M>::infinity();
        }
      }
    };

    template <std::totally_ordered M, M E>
    struct max<M, E> {
      using T = M;
      static T op(const T& x, const T& y) {
        assert(E <= x);
        assert(E <= y);
        return std::max(x, y);
      }
      static T e() {
        return E;
      }
    };

    template <typename M, M ...dummy>
    struct min;

    template <tools::arithmetic M>
    struct min<M> {
      using T = M;
      static T op(const T& x, const T& y) {
        return std::min(x, y);
      }
      static T e() {
        if constexpr (tools::integral<M>) {
          return std::numeric_limits<M>::max();
        } else {
          return std::numeric_limits<M>::infinity();
        }
      }
    };

    template <std::totally_ordered M, M E>
    struct min<M, E> {
      using T = M;
      static T op(const T& x, const T& y) {
        assert(x <= E);
        assert(y <= E);
        return std::min(x, y);
      }
      static T e() {
        return E;
      }
    };

    template <typename M>
    struct multiplies {
      using T = M;
      static T op(const T& x, const T& y) {
        return x * y;
      }
      static T e() {
        return T(1);
      }
    };

    template <>
    struct multiplies<bool> {
      using T = bool;
      static T op(const bool x, const bool y) {
        return x && y;
      }
      static T e() {
        return true;
      }
    };

    template <typename M, M E>
    struct update {
      using T = M;
      static T op(const T& x, const T& y) {
        return x == E ? y : x;
      }
      static T e() {
        return E;
      }
    };
  }
}


#line 1 "tools/semirings.hpp"



#line 8 "tools/semirings.hpp"

namespace tools {
  namespace semirings {
    template <tools::commutative_monoid A, tools::monoid M>
    struct of {
      using add = A;
      using mul = M;
    };

    template <typename R>
    using min_plus = tools::semirings::of<tools::monoids::min<R>, tools::groups::plus<R>>;

    template <typename R>
    using max_plus = tools::semirings::of<tools::monoids::max<R>, tools::groups::plus<R>>;

    template <typename R>
    using min_max = tools::semirings::of<tools::monoids::min<R>, tools::monoids::max<R>>;

    template <typename R>
    using max_min = tools::semirings::of<tools::monoids::max<R>, tools::monoids::min<R>>;
  }
}


#line 9 "tools/rings.hpp"

namespace tools {
  namespace rings {
    template <tools::commutative_group A, tools::monoid M>
    using of = tools::semirings::of<A, M>;

    template <typename R>
    using plus_multiplies = tools::rings::of<tools::groups::plus<R>, tools::monoids::multiplies<R>>;

    template <typename R>
    using xor_and = tools::rings::of<tools::groups::bit_xor<R>, tools::monoids::bit_and<R>>;
  }
}


#line 15 "tools/subset_convolution.hpp"

namespace tools {
  template <tools::ring R, std::ranges::input_range R1, std::ranges::input_range R2>
  requires std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R1>>
        && std::assignable_from<typename R::add::T&, std::ranges::range_value_t<R2>>
  auto subset_convolution(R1&& a, R2&& b) {
    if constexpr (std::ranges::sized_range<R1> && std::ranges::sized_range<R2>) {

      using Add = typename R::add;
      using Mul = typename R::mul;
      using T = typename Add::T;

      assert(std::has_single_bit(std::ranges::size(a)));
      assert(std::has_single_bit(std::ranges::size(b)));
      assert(std::ranges::size(a) == std::ranges::size(b));

      const int N = std::countr_zero(std::ranges::size(a));
      const int pow2_N = 1 << N;

      auto za = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      {
        int i = 0;
        for (auto&& a_i : a) {
          za[i][std::popcount<unsigned int>(i)] = a_i;
          ++i;
        }
      }

      auto zb = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      {
        int i = 0;
        for (auto&& b_i : b) {
          zb[i][std::popcount<unsigned int>(i)] = b_i;
          ++i;
        }
      }

      for (int w = 0; w < N; ++w) {
        for (int i = 0; i < pow2_N; ++i) {
          if (i & (1 << w)) {
            for (int j = 0; j <= N; ++j) {
              za[i][j] = Add::op(za[i][j], za[i ^ (1 << w)][j]);
              zb[i][j] = Add::op(zb[i][j], zb[i ^ (1 << w)][j]);
            }
          }
        }
      }

      auto zc = std::vector(pow2_N, std::vector(N + 1, Add::e()));
      for (int i = 0; i < pow2_N; ++i) {
        for (int j = 0; j <= N; ++j) {
          for (int k = 0; k <= j; ++k) {
            zc[i][j] = Add::op(zc[i][j], Mul::op(za[i][k], zb[i][j - k]));
          }
        }
      }

      for (int w = 0; w < N; ++w) {
        for (int i = 0; i < pow2_N; ++i) {
          if (i & (1 << w)) {
            for (int j = 0; j <= N; ++j) {
              zc[i][j] = Add::op(zc[i][j], Add::inv(zc[i ^ (1 << w)][j]));
            }
          }
        }
      }

      std::vector<T> c;
      c.reserve(pow2_N);
      for (int i = 0; i < pow2_N; ++i) {
        c.push_back(zc[i][std::popcount<unsigned int>(i)]);
      }

      return c;

    } else {
      return tools::subset_convolution<R>(
        std::forward<R1>(a) | std::ranges::to<std::vector<std::ranges::range_value_t<R1>>>(),
        std::forward<R2>(b) | std::ranges::to<std::vector<std::ranges::range_value_t<R2>>>()
      );
    }
  }

  template <std::ranges::input_range R1, std::ranges::input_range R2>
  auto subset_convolution(R1&& a, R2&& b) {
    using T = std::common_type_t<std::ranges::range_value_t<R1>, std::ranges::range_value_t<R2>>;
    return tools::subset_convolution<tools::rings::plus_multiplies<T>, R1, R2>(std::forward<R1>(a), std::forward<R2>(b));
  }
}

Back to top page