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

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

kronecker_power_transform

• シグネチャ

  template <typename T, std::size_t D, auto unit_transform>
  void kronecker_power_transform(std::vector<T> &x) // (1)
  
  template <typename T, typename UnitTransform>
  void kronecker_power_transform(std::vector<T> &x, const std::size_t D, UnitTransform unit_transform) // (2)
  
  template <typename T, auto e = default_operator::zero<T>, auto add = default_operator::add<T>, auto mul = default_operator::mul<T>>
  void kronecker_power_transform(std::vector<T> &x, const std::vector<std::vector<T>> &A) // (3)
  

• 概要

  $D \times D$ 行列 $A$ のクロネッカー積における $N$ 乗 $A ^ {\otimes N}$ と $D ^ N$ 次元ベクトル $x$ の積を inplace に $\Theta ( N \cdot D ^ N )$ 時間で計算する。

  ゼータ変換・メビウス変換・アダマール変換らはこの具体例であり、以下のように対応付けられる。

  • 上位集合版のゼータ変換

    $D=2,A=\begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}$

  • 上位集合版のメビウス変換

    $D=2,A=\begin{pmatrix} 1 & -1 \\ 0 & 1 \end{pmatrix}$

  • 下位集合版のゼータ変換

    $D=2,A=\begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}$

  • 下位集合版のメビウス変換

    $D=2,A=\begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}$

  • アダマール変換

    $D=2,A=\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$

  これらは全て $D=2$ の例であるが、例えば $3$ 進ゼータ変換は $D=3,A=\begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{pmatrix}$、メビウス変換は $D=3,A=\begin{pmatrix} 1 & 0 & 0 \\ -1 & 1 & 0 \\ 0 & -1 & 1 \end{pmatrix}$ となる (はず)。

• 使い方

  1. $D$ の値がコンパイル時に分かっている場合はこれを用いるのがパフォーマンス的に最適 (なはず)。テンプレート引数 unit_transform は引数に $x _ 1,\ldots, x _ D$ を $A x$ に置き換える関数である。具体的には、$D$ 個の参照 T &x1, ..., T &xD を引数に取る関数で、各 $i = 1, \ldots, D$ に対して $x_i \leftarrow (Ax)_i$ という代入が行われなければならない。返り値は捨てられる。
  2. $D$ の値がコンパイル時に分からない場合であって、行列積をそのまま計算すると無駄が多い場合はこれを用いる。unit_transform はサイズ D の vector<T> & を受け取り、1 と同様に中身を $Ax$ で置き換える必要がある。
  3. $D$ の値がコンパイル時に分からない場合であって、行列をそのまま与えてしまいたい場合はこれを用いる。加算や乗算、加法単位元をテンプレート引数に渡すことができる。デフォルトでは加算は operator+、乗算は operator*、加法単位元は T{0} として扱われる。

• 返り値

  1, 2, 3 ともに無し。変換は inplace に行われる。

• 制約

  • 引数 x のサイズ $L$ に対してある非負整数 $N$ が存在して $L=D^N$ を満たす。

• 時間計算量

  $\Theta(N\cdot D ^ N)$

Depends on

• Default Operator (library/util/default_operator.hpp)

Required by

• Bitwise And Convolution (library/convolution/and_convolution.hpp)
• Array Subset Convolution (library/convolution/array_subset_convolution.hpp)
• Bitwise Or Convolution (library/convolution/or_convolution.hpp)
• Subset Convolution (library/convolution/subset_convolution.hpp)
• Bitwise Xor Convolution (library/convolution/xor_convolution.hpp)
• Hafnian (完全マッチングの数え上げ) (library/linear_algebra/hafnian.hpp)
• Array Sps (library/math/array_set_power_series.hpp)
• Set Power Series (library/math/set_power_series.hpp)
• 下位集合に対する高速ゼータ変換・高速メビウス変換 (library/transform/subset.hpp)
• 上位集合に対する高速ゼータ変換・高速メビウス変換 (library/transform/supset.hpp)
• Walsh Hadamard 変換 (library/transform/walsh_hadamard.hpp)

Verified with

• test/src/convolution/and_convolution/and_convolution.test.cpp
• test/src/convolution/polynomial_eval/nim_counting.test.cpp
• test/src/convolution/polynomial_eval_multipoint_eval/nim_counting.test.cpp
• test/src/convolution/subset_convolution/subset_convolution.test.cpp
• test/src/convolution/xor_convolution/xor_convolution.test.cpp
• test/src/linear_algebra/hafnian/hafnian_of_matrix.test.cpp
• test/src/math/array_set_power_series/abc213_g.test.cpp
• test/src/math/array_set_power_series/abc236_h.test.cpp
• test/src/math/array_set_power_series/abc253_h.test.cpp
• test/src/math/array_set_power_series/abc253_h_2.test.cpp
• test/src/math/array_set_power_series/arc105_f.test.cpp
• test/src/math/set_power_series/abc213_g.test.cpp
• test/src/math/set_power_series/abc236_h.test.cpp
• test/src/math/set_power_series/abc253_h.test.cpp
• test/src/math/set_power_series/abc253_h_2.test.cpp
• test/src/math/set_power_series/arc105_f.test.cpp
• test/src/math/set_power_series/exp_of_set_power_series.test.cpp
• test/src/math/set_power_series/polynomial_composite_set_power_series.test.cpp
• test/src/transform/kronecker_power/agc044_c.test.cpp
• test/src/transform/kronecker_power/arc132_f.test.cpp

Code

#ifndef SUISEN_KRONECKER_POWER
#define SUISEN_KRONECKER_POWER

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

#include "library/util/default_operator.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


#endif // SUISEN_KRONECKER_POWER

#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

Back to top page