Bitwise OR convolution (tools/or_convolution.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2024-10-08 23:47:47+09:00
• Include: #include "tools/or_convolution.hpp"

template <typename InputIterator, typename RandomAccessIterator>
void or_convolution(InputIterator a_begin, InputIterator a_end, InputIterator b_begin, InputIterator b_end, RandomAccessIterator c_begin, RandomAccessIterator c_end);

Given two infinite sequences $(a_0, a_1, \ldots, a_{N - 1}, 0, 0, \ldots)$ and $(b_0, b_1, \ldots, b_{M - 1}, 0, 0, \ldots)$, it returns the first $K$ terms of the infinite sequence $(c_0, c_1, \ldots)$ where

\[\begin{align*} N &= \mathrm{a\_end} - \mathrm{a\_begin}\\ M &= \mathrm{b\_end} - \mathrm{b\_begin}\\ K &= \mathrm{c\_end} - \mathrm{c\_begin}\\ c_k &= \sum_{(i~\mathrm{OR}~j) = k} a_i b_j \end{align*}\]

Constraints

• InputIterator is an input iterator type.
• RandomAccessIterator is a random access iterator type.

Time Complexity

• $O(N + M + K \log K)$

License

• CC0

Author

• anqooqie

Depends on

• $2^x$ (tools/pow2.hpp)
• Subset Möbius transform (tools/subset_moebius.hpp)
• Subset Zeta transform (tools/subset_zeta.hpp)

Verified with

• tests/or_convolution/different_lengths.test.cpp
• tests/or_convolution/regular.test.cpp

Code

#ifndef TOOLS_OR_CONVOLUTION_HPP
#define TOOLS_OR_CONVOLUTION_HPP

#include <iterator>
#include <vector>
#include "tools/subset_zeta.hpp"
#include "tools/subset_moebius.hpp"

namespace tools {
  template <typename InputIterator, typename RandomAccessIterator>
  void or_convolution(InputIterator a_begin, InputIterator a_end, InputIterator b_begin, InputIterator b_end, RandomAccessIterator c_begin, RandomAccessIterator c_end) {
    using T = typename ::std::iterator_traits<InputIterator>::value_type;
    ::std::vector<T> a(a_begin, a_end);
    ::std::vector<T> b(b_begin, b_end);
    const int K = ::std::distance(c_begin, c_end);

    a.resize(K, T(0));
    ::tools::subset_zeta(a.begin(), a.end());
    b.resize(K, T(0));
    ::tools::subset_zeta(b.begin(), b.end());

    for (int i = 0; i < K; ++i) {
      c_begin[i] = a[i] * b[i];
    }
    ::tools::subset_moebius(c_begin, c_end);
  }
}

#endif

#line 1 "tools/or_convolution.hpp"



#include <iterator>
#include <vector>
#line 1 "tools/subset_zeta.hpp"



#line 6 "tools/subset_zeta.hpp"
#include <algorithm>
#line 1 "tools/pow2.hpp"



#include <type_traits>
#include <cstddef>

namespace tools {

  template <typename T, typename ::std::enable_if<::std::is_unsigned<T>::value, ::std::nullptr_t>::type = nullptr>
  constexpr T pow2(const T x) {
    return static_cast<T>(1) << x;
  }

  template <typename T, typename ::std::enable_if<::std::is_signed<T>::value, ::std::nullptr_t>::type = nullptr>
  constexpr T pow2(const T x) {
    return static_cast<T>(static_cast<typename ::std::make_unsigned<T>::type>(1) << static_cast<typename ::std::make_unsigned<T>::type>(x));
  }
}


#line 8 "tools/subset_zeta.hpp"

namespace tools {
  template <typename RandomAccessIterator>
  void subset_zeta(const RandomAccessIterator begin, const RandomAccessIterator end) {
    const int N = end - begin;

    for (int w = 0; ::tools::pow2(w) < N; ++w) {
      for (int i = 0; ; i += ::tools::pow2(w)) {
        for (; !((i >> w) & 1); ++i) {
          if (!(i + ::tools::pow2(w) < N)) goto NEXT_W;
          begin[i + ::tools::pow2(w)] += begin[i];
        }
      }
    NEXT_W:
      ;
    }
  }

  template <typename InputIterator, typename OutputIterator>
  void subset_zeta(const InputIterator begin, const InputIterator end, const OutputIterator result) {
    using T = typename ::std::iterator_traits<InputIterator>::value_type;
    ::std::vector<T> a(begin, end);
    ::tools::subset_zeta(a.begin(), a.end());
    ::std::move(a.begin(), a.end(), result);
  }
}


#line 1 "tools/subset_moebius.hpp"



#line 8 "tools/subset_moebius.hpp"

namespace tools {
  template <typename RandomAccessIterator>
  void subset_moebius(const RandomAccessIterator begin, const RandomAccessIterator end) {
    const int N = end - begin;

    for (int w = 0; ::tools::pow2(w) < N; ++w) {
      for (int i = 0; ; i += ::tools::pow2(w)) {
        for (; !((i >> w) & 1); ++i) {
          if (!(i + ::tools::pow2(w) < N)) goto NEXT_W;
          begin[i + ::tools::pow2(w)] -= begin[i];
        }
      }
    NEXT_W:
      ;
    }
  }

  template <typename InputIterator, typename OutputIterator>
  void subset_moebius(const InputIterator begin, const InputIterator end, const OutputIterator result) {
    using T = typename ::std::iterator_traits<InputIterator>::value_type;
    ::std::vector<T> b(begin, end);
    ::tools::subset_moebius(b.begin(), b.end());
    ::std::move(b.begin(), b.end(), result);
  }
}


#line 8 "tools/or_convolution.hpp"

namespace tools {
  template <typename InputIterator, typename RandomAccessIterator>
  void or_convolution(InputIterator a_begin, InputIterator a_end, InputIterator b_begin, InputIterator b_end, RandomAccessIterator c_begin, RandomAccessIterator c_end) {
    using T = typename ::std::iterator_traits<InputIterator>::value_type;
    ::std::vector<T> a(a_begin, a_end);
    ::std::vector<T> b(b_begin, b_end);
    const int K = ::std::distance(c_begin, c_end);

    a.resize(K, T(0));
    ::tools::subset_zeta(a.begin(), a.end());
    b.resize(K, T(0));
    ::tools::subset_zeta(b.begin(), b.end());

    for (int i = 0; i < K; ++i) {
      c_begin[i] = a[i] * b[i];
    }
    ::tools::subset_moebius(c_begin, c_end);
  }
}

Back to top page