Bitwise AND convolution (tools/and_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/and_convolution.hpp"

template <typename InputIterator, typename RandomAccessIterator>
void and_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{AND}~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 \log N + M \log M + K)$

License

Author

anqooqie

Depends on

Verified with

Code

#ifndef TOOLS_AND_CONVOLUTION_HPP
#define TOOLS_AND_CONVOLUTION_HPP

#include <iterator>
#include <vector>
#include <algorithm>
#include "tools/superset_zeta.hpp"
#include "tools/superset_moebius.hpp"

namespace tools {
  template <typename InputIterator, typename RandomAccessIterator>
  void and_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 N = a.size();
    const int M = b.size();
    const int K = ::std::distance(c_begin, c_end);

    ::tools::superset_zeta(a.begin(), a.end());
    ::tools::superset_zeta(b.begin(), b.end());

    if (::std::min(N, M) <= K) {
      for (int i = 0; i < ::std::min(N, M); ++i) {
        c_begin[i] = a[i] * b[i];
      }
      ::tools::superset_moebius(c_begin, c_begin + ::std::min(N, M));
      ::std::fill(c_begin + ::std::min(N, M), c_end, T(0));
    } else {
      ::std::vector<T> c;
      c.reserve(::std::min(N, M));
      for (int i = 0; i < ::std::min(N, M); ++i) {
        c.push_back(a[i] * b[i]);
      }
      ::tools::superset_moebius(c.begin(), c.end());
      ::std::move(c.begin(), c.begin() + K, c_begin);
    }
  }
}

#endif

#line 1 "tools/and_convolution.hpp"



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



#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/superset_zeta.hpp"

namespace tools {
  template <typename RandomAccessIterator>
  void superset_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] += begin[i + ::tools::pow2(w)];
        }
      }
    NEXT_W:
      ;
    }
  }

  template <typename InputIterator, typename OutputIterator>
  void superset_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::superset_zeta(a.begin(), a.end());
    ::std::move(a.begin(), a.end(), result);
  }
}


#line 1 "tools/superset_moebius.hpp"



#line 8 "tools/superset_moebius.hpp"

namespace tools {
  template <typename RandomAccessIterator>
  void superset_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] -= begin[i + ::tools::pow2(w)];
        }
      }
    NEXT_W:
      ;
    }
  }

  template <typename InputIterator, typename OutputIterator>
  void superset_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::superset_moebius(b.begin(), b.end());
    ::std::move(b.begin(), b.end(), result);
  }
}


#line 9 "tools/and_convolution.hpp"

namespace tools {
  template <typename InputIterator, typename RandomAccessIterator>
  void and_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 N = a.size();
    const int M = b.size();
    const int K = ::std::distance(c_begin, c_end);

    ::tools::superset_zeta(a.begin(), a.end());
    ::tools::superset_zeta(b.begin(), b.end());

    if (::std::min(N, M) <= K) {
      for (int i = 0; i < ::std::min(N, M); ++i) {
        c_begin[i] = a[i] * b[i];
      }
      ::tools::superset_moebius(c_begin, c_begin + ::std::min(N, M));
      ::std::fill(c_begin + ::std::min(N, M), c_end, T(0));
    } else {
      ::std::vector<T> c;
      c.reserve(::std::min(N, M));
      for (int i = 0; i < ::std::min(N, M); ++i) {
        c.push_back(a[i] * b[i]);
      }
      ::tools::superset_moebius(c.begin(), c.end());
      ::std::move(c.begin(), c.begin() + K, c_begin);
    }
  }
}