上位集合に対する高速ゼータ変換・高速メビウス変換
(library/transform/supset.hpp)

• View this file on GitHub
• Last update: 2022-03-21 02:24:20+09:00
• Include: #include "library/transform/supset.hpp"

• シグネチャ

  template <typename T, auto add = default_operator::add<T>>
  void zeta(std::vector<T> &a)
  template <typename T, auto sub = default_operator::sub<T>>
  void mobius(std::vector<T> &a)
  

• 概要

  長さ $N=2^L$ の列 $(A_0=0,A_1,\ldots,A_{N-1})$ に上位集合版のゼータ変換を施す関数 zeta およびその逆変換 (メビウス変換) を施す関数 mobius を提供します．各変換は inplace に行われ，引数として渡した列は書き換えられます．

  ゼータ変換では，引数として与えた列 $A$ に対して，以下を満たす列 $B$ を計算します．

  \[B_i = \sum_{i\subseteq j} A_j\]

  メビウス変換では，引数として与えた列 $A$ に対して，以下を満たす列 $B$ を計算します．

  \[A_i = \sum_{i\subseteq j} B_j\]

• テンプレート引数

  • T: 列の要素の型．
  • add: 二項演算 (加算)．デフォルトでは operator+ が呼ばれるようになっています．
  • sub: 二項演算 (減算)．デフォルトでは operator- が呼ばれるようになっています．

• 制約

  • ある非負整数 $L$ が存在して，$\vert A\vert = 2 ^ L$
  • $0\leq L\leq 20$

• 時間計算量

  $\Theta(N\log N)$

Depends on

• クロネッカー冪による線形変換 (仮称) (library/transform/kronecker_power.hpp)
• Default Operator (library/util/default_operator.hpp)

Required by

• Bitwise And Convolution (library/convolution/and_convolution.hpp)

Verified with

• test/src/convolution/and_convolution/and_convolution.test.cpp

Code

#ifndef SUISEN_SUPSET_TRANSFORM
#define SUISEN_SUPSET_TRANSFORM

#include "library/transform/kronecker_power.hpp"

namespace suisen::supset_transform {
    namespace internal {
        template <typename T, auto add = default_operator::add<T>>
        void zeta_unit_transform(T &x0, T &x1) {
            x0 = add(x0, x1);   // 1, 1
                                // 0, 1
        }
        template <typename T, auto sub = default_operator::sub<T>>
        void mobius_unit_transform(T &x0, T &x1) {
            x0 = sub(x0, x1);   // 1, -1
                                // 0,  1
        }
    } // namespace internal

    using kronecker_power_transform::kronecker_power_transform;

    template <typename T, auto add = default_operator::add<T>>
    void zeta(std::vector<T> &a) {
        kronecker_power_transform<T, 2, internal::zeta_unit_transform<T, add>>(a);
    }
    template <typename T, auto sub = default_operator::sub<T>>
    void mobius(std::vector<T> &a) {
        kronecker_power_transform<T, 2, internal::mobius_unit_transform<T, sub>>(a);
    }
} // namespace suisen::supset_transform

#endif // SUISEN_SUPSET_TRANSFORM

#line 1 "library/transform/supset.hpp"



#line 1 "library/transform/kronecker_power.hpp"



#include <cassert>
#include <utility>
#include <vector>

#line 1 "library/util/default_operator.hpp"



namespace suisen {
    namespace default_operator {
        template <typename T>
        auto zero() -> decltype(T { 0 }) { return T { 0 }; }
        template <typename T>
        auto one()  -> decltype(T { 1 }) { return T { 1 }; }
        template <typename T>
        auto add(const T &x, const T &y) -> decltype(x + y) { return x + y; }
        template <typename T>
        auto sub(const T &x, const T &y) -> decltype(x - y) { return x - y; }
        template <typename T>
        auto mul(const T &x, const T &y) -> decltype(x * y) { return x * y; }
        template <typename T>
        auto div(const T &x, const T &y) -> decltype(x / y) { return x / y; }
        template <typename T>
        auto mod(const T &x, const T &y) -> decltype(x % y) { return x % y; }
        template <typename T>
        auto neg(const T &x) -> decltype(-x) { return -x; }
        template <typename T>
        auto inv(const T &x) -> decltype(one<T>() / x)  { return one<T>() / x; }
    } // default_operator
    namespace default_operator_noref {
        template <typename T>
        auto zero() -> decltype(T { 0 }) { return T { 0 }; }
        template <typename T>
        auto one()  -> decltype(T { 1 }) { return T { 1 }; }
        template <typename T>
        auto add(T x, T y) -> decltype(x + y) { return x + y; }
        template <typename T>
        auto sub(T x, T y) -> decltype(x - y) { return x - y; }
        template <typename T>
        auto mul(T x, T y) -> decltype(x * y) { return x * y; }
        template <typename T>
        auto div(T x, T y) -> decltype(x / y) { return x / y; }
        template <typename T>
        auto mod(T x, T y) -> decltype(x % y) { return x % y; }
        template <typename T>
        auto neg(T x) -> decltype(-x) { return -x; }
        template <typename T>
        auto inv(T x) -> decltype(one<T>() / x)  { return one<T>() / x; }
    } // default_operator
} // namespace suisen


#line 9 "library/transform/kronecker_power.hpp"

namespace suisen {
    namespace kronecker_power_transform {
        namespace internal {
            template <typename UnitTransform, typename ReferenceGetter, std::size_t... Seq>
            void unit_transform(UnitTransform transform, ReferenceGetter ref_getter, std::index_sequence<Seq...>) {
                transform(ref_getter(Seq)...);
            }
        }

        template <typename T, std::size_t D, auto unit_transform>
        void kronecker_power_transform(std::vector<T> &x) {
            const std::size_t n = x.size();
            for (std::size_t block = 1; block < n; block *= D) {
                for (std::size_t l = 0; l < n; l += D * block) {
                    for (std::size_t offset = l; offset < l + block; ++offset) {
                        const auto ref_getter = [&](std::size_t i) -> T& { return x[offset + i * block]; };
                        internal::unit_transform(unit_transform, ref_getter, std::make_index_sequence<D>());
                    }
                }
            }
        }

        template <typename T, typename UnitTransform>
        void kronecker_power_transform(std::vector<T> &x, const std::size_t D, UnitTransform unit_transform) {
            const std::size_t n = x.size();
            std::vector<T> work(D);
            for (std::size_t block = 1; block < n; block *= D) {
                for (std::size_t l = 0; l < n; l += D * block) {
                    for (std::size_t offset = l; offset < l + block; ++offset) {
                        for (std::size_t i = 0; i < D; ++i) work[i] = x[offset + i * block];
                        unit_transform(work);
                        for (std::size_t i = 0; i < D; ++i) x[offset + i * block] = work[i];
                    }
                }
            }
        }

        template <typename T, auto e = default_operator::zero<T>, auto add = default_operator::add<T>, auto mul = default_operator::mul<T>>
        auto kronecker_power_transform(std::vector<T> &x, const std::vector<std::vector<T>> &A) -> decltype(e(), add(std::declval<T>(), std::declval<T>()), mul(std::declval<T>(), std::declval<T>()), void()) {
            const std::size_t D = A.size();
            assert(D == A[0].size());
            auto unit_transform = [&](std::vector<T> &x) {
                std::vector<T> y(D, e());
                for (std::size_t i = 0; i < D; ++i) for (std::size_t j = 0; j < D; ++j) {
                    y[i] = add(y[i], mul(A[i][j], x[j]));
                }
                x.swap(y);
            };
            kronecker_power_transform<T>(x, D, unit_transform);
        }
    }
} // namespace suisen



#line 5 "library/transform/supset.hpp"

namespace suisen::supset_transform {
    namespace internal {
        template <typename T, auto add = default_operator::add<T>>
        void zeta_unit_transform(T &x0, T &x1) {
            x0 = add(x0, x1);   // 1, 1
                                // 0, 1
        }
        template <typename T, auto sub = default_operator::sub<T>>
        void mobius_unit_transform(T &x0, T &x1) {
            x0 = sub(x0, x1);   // 1, -1
                                // 0,  1
        }
    } // namespace internal

    using kronecker_power_transform::kronecker_power_transform;

    template <typename T, auto add = default_operator::add<T>>
    void zeta(std::vector<T> &a) {
        kronecker_power_transform<T, 2, internal::zeta_unit_transform<T, add>>(a);
    }
    template <typename T, auto sub = default_operator::sub<T>>
    void mobius(std::vector<T> &a) {
        kronecker_power_transform<T, 2, internal::mobius_unit_transform<T, sub>>(a);
    }
} // namespace suisen::supset_transform

Back to top page