Fast Zeta Transform
(src/algebra/convolution/zeta.hpp)

View this file on GitHub
Last update: 2021-05-31 22:43:54+09:00
Include: #include "src/algebra/convolution/zeta.hpp"

Required by

Verified with

test/library-checker/bitwise_and_convolution.test.cpp

Code

#pragma once

/**
 * @file zeta.hpp
 * @brief Fast Zeta Transform
 */

#include <cstddef>

#include "lib/bit"

namespace workspace {

/**
 * @brief Fast zeta transform on a lattice.
 * @param f Value list of the function
 * @param n Number of join-irreducible elements, bit-width of indices
 * @return Value list of resulting function.
 */
template <class A> A fast_zeta(A f, size_t n) {
  while (n--)
    for (size_t s = 0; s != std::size(f); ++s)
      if (s >> n & 1) f[s] += f[s ^ 1 << n];
  return f;
}

/**
 * @brief Fast zeta transform on a lattice.
 * @param f Value list of the function
 * @return Value list of resulting function.
 */
template <class A> A fast_zeta(A const &f) {
  return fast_zeta(f, std::__countr_zero(std::__bit_floor(std::size(f))));
}

/**
 * @brief Inverse of fast zeta transform.
 * @param f Value list of the function
 * @param n Number of join-irreducible elements, bit-width of indices
 * @return Value list of resulting function.
 */
template <class A> A fast_mobius(A f, size_t n) {
  while (n--)
    for (size_t s = 0; s != std::size(f); ++s)
      if (s >> n & 1) f[s] -= f[s ^ 1 << n];
  return f;
}

/**
 * @brief Inverse of fast zeta transform.
 * @param f Value list of the function
 * @return Value list of resulting function.
 */
template <class A> A fast_mobius(A const &f) {
  return fast_mobius(f, std::__countr_zero(std::__bit_floor(std::size(f))));
}

}  // namespace workspace

#line 2 "src/algebra/convolution/zeta.hpp"

/**
 * @file zeta.hpp
 * @brief Fast Zeta Transform
 */

#include <cstddef>

#line 2 "lib/bit"

#if __cplusplus > 201703L

#include <bit>

#elif __cplusplus > 201402L

#ifndef _GLIBCXX_BIT
#define _GLIBCXX_BIT 1

#include <limits>
#include <type_traits>

namespace std {

template <typename _Tp> constexpr int __countl_zero(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;

  if (__x == 0) return _Nd;

  constexpr auto _Nd_ull = numeric_limits<unsigned long long>::digits;
  constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits;
  constexpr auto _Nd_u = numeric_limits<unsigned>::digits;

  if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u) {
    constexpr int __diff = _Nd_u - _Nd;
    return __builtin_clz(__x) - __diff;
  } else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul) {
    constexpr int __diff = _Nd_ul - _Nd;
    return __builtin_clzl(__x) - __diff;
  } else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull) {
    constexpr int __diff = _Nd_ull - _Nd;
    return __builtin_clzll(__x) - __diff;
  } else  // (_Nd > _Nd_ull)
  {
    static_assert(_Nd <= (2 * _Nd_ull),
                  "Maximum supported integer size is 128-bit");

    unsigned long long __high = __x >> _Nd_ull;
    if (__high != 0) {
      constexpr int __diff = (2 * _Nd_ull) - _Nd;
      return __builtin_clzll(__high) - __diff;
    }
    constexpr auto __max_ull = numeric_limits<unsigned long long>::max();
    unsigned long long __low = __x & __max_ull;
    return (_Nd - _Nd_ull) + __builtin_clzll(__low);
  }
}

template <typename _Tp> constexpr int __countr_zero(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;

  if (__x == 0) return _Nd;

  constexpr auto _Nd_ull = numeric_limits<unsigned long long>::digits;
  constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits;
  constexpr auto _Nd_u = numeric_limits<unsigned>::digits;

  if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u)
    return __builtin_ctz(__x);
  else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul)
    return __builtin_ctzl(__x);
  else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull)
    return __builtin_ctzll(__x);
  else  // (_Nd > _Nd_ull)
  {
    static_assert(_Nd <= (2 * _Nd_ull),
                  "Maximum supported integer size is 128-bit");

    constexpr auto __max_ull = numeric_limits<unsigned long long>::max();
    unsigned long long __low = __x & __max_ull;
    if (__low != 0) return __builtin_ctzll(__low);
    unsigned long long __high = __x >> _Nd_ull;
    return __builtin_ctzll(__high) + _Nd_ull;
  }
}

template <typename _Tp> constexpr int __popcount(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;

  if (__x == 0) return 0;

  constexpr auto _Nd_ull = numeric_limits<unsigned long long>::digits;
  constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits;
  constexpr auto _Nd_u = numeric_limits<unsigned>::digits;

  if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u)
    return __builtin_popcount(__x);
  else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul)
    return __builtin_popcountl(__x);
  else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull)
    return __builtin_popcountll(__x);
  else  // (_Nd > _Nd_ull)
  {
    static_assert(_Nd <= (2 * _Nd_ull),
                  "Maximum supported integer size is 128-bit");

    constexpr auto __max_ull = numeric_limits<unsigned long long>::max();
    unsigned long long __low = __x & __max_ull;
    unsigned long long __high = __x >> _Nd_ull;
    return __builtin_popcountll(__low) + __builtin_popcountll(__high);
  }
}

template <typename _Tp> constexpr _Tp __bit_ceil(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;
  if (__x == 0 || __x == 1) return 1;
  auto __shift_exponent = _Nd - __countl_zero((_Tp)(__x - 1u));
#ifdef _GLIBCXX_HAVE_BUILTIN_IS_CONSTANT_EVALUATED
  if (!__builtin_is_constant_evaluated()) {
    __glibcxx_assert(__shift_exponent != numeric_limits<_Tp>::digits);
  }
#endif
  using __promoted_type = decltype(__x << 1);
  if _GLIBCXX17_CONSTEXPR (!is_same<__promoted_type, _Tp>::value) {
    const int __extra_exp = sizeof(__promoted_type) / sizeof(_Tp) / 2;
    __shift_exponent |= (__shift_exponent & _Nd) << __extra_exp;
  }
  return (_Tp)1u << __shift_exponent;
}

template <typename _Tp> constexpr _Tp __bit_floor(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;
  if (__x == 0) return 0;
  return (_Tp)1u << (_Nd - __countl_zero((_Tp)(__x >> 1)));
}

template <typename _Tp> constexpr _Tp __bit_width(_Tp __x) noexcept {
  constexpr auto _Nd = numeric_limits<_Tp>::digits;
  return _Nd - __countl_zero(__x);
}

}  // namespace std

#endif

#endif
#line 11 "src/algebra/convolution/zeta.hpp"

namespace workspace {

/**
 * @brief Fast zeta transform on a lattice.
 * @param f Value list of the function
 * @param n Number of join-irreducible elements, bit-width of indices
 * @return Value list of resulting function.
 */
template <class A> A fast_zeta(A f, size_t n) {
  while (n--)
    for (size_t s = 0; s != std::size(f); ++s)
      if (s >> n & 1) f[s] += f[s ^ 1 << n];
  return f;
}

/**
 * @brief Fast zeta transform on a lattice.
 * @param f Value list of the function
 * @return Value list of resulting function.
 */
template <class A> A fast_zeta(A const &f) {
  return fast_zeta(f, std::__countr_zero(std::__bit_floor(std::size(f))));
}

/**
 * @brief Inverse of fast zeta transform.
 * @param f Value list of the function
 * @param n Number of join-irreducible elements, bit-width of indices
 * @return Value list of resulting function.
 */
template <class A> A fast_mobius(A f, size_t n) {
  while (n--)
    for (size_t s = 0; s != std::size(f); ++s)
      if (s >> n & 1) f[s] -= f[s ^ 1 << n];
  return f;
}

/**
 * @brief Inverse of fast zeta transform.
 * @param f Value list of the function
 * @return Value list of resulting function.
 */
template <class A> A fast_mobius(A const &f) {
  return fast_mobius(f, std::__countr_zero(std::__bit_floor(std::size(f))));
}

}  // namespace workspace

Fast Zeta Transform (src/algebra/convolution/zeta.hpp)

Required by

Verified with

Code

Fast Zeta Transform
(src/algebra/convolution/zeta.hpp)