Library/Math/Convolution.hpp

• View this file on GitHub
• Last update: 2024-12-27 17:07:26+09:00
• Include: #include "Library/Math/Convolution.hpp"

Depends on

• Library/Math/Bit.hpp
• Library/Math/Mobius.hpp
• Library/Math/Zeta.hpp

Verified with

• Test/Math/Convolution_and.test.cpp

Code

#pragma once

#include <ranges>
#include <vector>

#include "./Mobius.hpp"
#include "./Zeta.hpp"

namespace mtd::convolution {

  template <class T>
  auto bitwise_and(const std::vector<T>& a, const std::vector<T>& b) {
    auto za = mtd::zeta::bit_supset(a);
    auto zb = mtd::zeta::bit_supset(b);
    auto zab = std::vector<T>();
    for (auto i : std::views::iota(static_cast<size_t>(0), a.size())) {
      zab.emplace_back(za[i] * zb[i]);
    }
    auto zma = mtd::mobius::bit_supset(za);
    return mtd::mobius::bit_supset(zab);
  }

}  // namespace mtd::convolution

#line 2 "Library/Math/Convolution.hpp"

#include <ranges>
#include <vector>

#line 2 "Library/Math/Mobius.hpp"

#line 5 "Library/Math/Mobius.hpp"

#line 2 "Library/Math/Bit.hpp"

namespace mtd {

  constexpr unsigned clz(unsigned int x) {
    int c = 0;
    if (x == 0) { return 0; }
    if (x & 0xffff0000) {
      x &= 0xffff0000;
      c |= 0x10;
    }
    if (x & 0xff00ff00) {
      x &= 0xff00ff00;
      c |= 0x08;
    }
    if (x & 0xf0f0f0f0) {
      x &= 0xf0f0f0f0;
      c |= 0x04;
    }
    if (x & 0xcccccccc) {
      x &= 0xcccccccc;
      c |= 0x02;
    }
    if (x & 0xaaaaaaaa) { c |= 0x01; }
    return c + 1;
  }

  constexpr unsigned ctz(unsigned int n) {
    if (!n) return -1;
    unsigned int c = 32;
    n &= -static_cast<signed int>(n);
    if (n) c--;
    if (n & 0x0000FFFF) c -= 16;
    if (n & 0x00FF00FF) c -= 8;
    if (n & 0x0F0F0F0F) c -= 4;
    if (n & 0x33333333) c -= 2;
    if (n & 0x55555555) c -= 1;
    return c;
  }

  constexpr unsigned long long popcount(unsigned long long x) {
    x = x - ((x >> 1) & 0x5555555555555555);
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333);
    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f;
    x = x + (x >> 8);
    x = x + (x >> 16);
    x = x + (x >> 32);
    return x & 0x0000007f;
  }

}  // namespace mtd
#line 7 "Library/Math/Mobius.hpp"

namespace mtd::mobius {

  template <class T>
  auto n(const std::vector<T>& a) {
    auto ret = a;
    for (auto i : std::views::iota(static_cast<size_t>(1), a.size())) {
      ret[i] = a[i] - a[i - 1];
    }
    return ret;
  }

  template <class T>
  auto bit_subset(const std::vector<T>& a) {
    auto ret = a;
    int size = clz(a.size());
    for (auto b : std::views::iota(0, size)) {
      for (auto bit : std::views::iota(0, 1LL << size)) {
        if (((bit >> b) & 1) && bit < a.size()) {
          ret[bit] -= ret[bit ^ (1LL << b)];
        }
      }
    }
    return ret;
  }

  template <class T>
  auto bit_supset(const std::vector<T>& a) {
    auto ret = a;
    int size = clz(a.size());
    for (auto b : std::views::iota(0, size)) {
      for (auto bit : std::views::iota(0, 1LL << size)) {
        if (((bit >> b) & 1) && bit < a.size()) {
          ret[bit ^ (1LL << b)] -= ret[bit];
        }
      }
    }
    return ret;
  }

}  // namespace mtd::mobius
#line 2 "Library/Math/Zeta.hpp"

#line 5 "Library/Math/Zeta.hpp"

#line 7 "Library/Math/Zeta.hpp"

namespace mtd::zeta {

  template <class T>
  auto n(const std::vector<T>& a) {
    auto ret = a;
    for (auto i : std::views::iota(static_cast<size_t>(1), a.size())) {
      ret[i] += ret[i - 1];
    }
    return ret;
  }

  template <class T>
  auto bit_subset(const std::vector<T>& a) {
    auto ret = a;
    int size = clz(a.size());
    for (auto b : std::views::iota(0, size)) {
      for (auto bit : std::views::iota(0, 1LL << size)) {
        if (((bit >> b) & 1) && bit < a.size()) {
          ret[bit] += ret[bit ^ (1LL << b)];
        }
      }
    }
    return ret;
  }

  template <class T>
  auto bit_supset(const std::vector<T>& a) {
    auto ret = a;
    int size = clz(a.size());
    for (auto b : std::views::iota(0, size)) {
      for (auto bit : std::views::iota(0, 1LL << size)) {
        if (((bit >> b) & 1) && bit < a.size()) {
          ret[bit ^ (1LL << b)] += ret[bit];
        }
      }
    }
    return ret;
  }
}  // namespace mtd::zeta
#line 8 "Library/Math/Convolution.hpp"

namespace mtd::convolution {

  template <class T>
  auto bitwise_and(const std::vector<T>& a, const std::vector<T>& b) {
    auto za = mtd::zeta::bit_supset(a);
    auto zb = mtd::zeta::bit_supset(b);
    auto zab = std::vector<T>();
    for (auto i : std::views::iota(static_cast<size_t>(0), a.size())) {
      zab.emplace_back(za[i] * zb[i]);
    }
    auto zma = mtd::mobius::bit_supset(za);
    return mtd::mobius::bit_supset(zab);
  }

}  // namespace mtd::convolution

Back to top page