約数系ゼータ変換・メビウス変換
(library/transform/divisor.hpp)

• View this file on GitHub
• Last update: 2022-01-31 13:34:34+09:00
• Include: #include "library/transform/divisor.hpp"

• シグネチャ

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

• 概要

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

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

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

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

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

  ここで，$j\mid i$ は $j$ が $i$ を割り切ること，言い換えれば $j$ が $i$ の約数であることを表します．

• テンプレート引数

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

• 制約

  • $1\leq N \leq 10^7$

• 時間計算量

  $\Theta(N\log\log N)$

Depends on

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

Required by

• LCM Convolution (library/convolution/lcm_convolution.hpp)

Verified with

• test/src/convolution/lcm_convolution/lcm_convolution.test.cpp
• test/src/transform/multiple/divide_both.test.cpp

Code

#ifndef SUISEN_DIVISOR_TRANSFORM
#define SUISEN_DIVISOR_TRANSFORM

#include <vector>
#include "library/util/default_operator.hpp"

namespace suisen::divisor_transform {
    // Calculates `g` s.t. g(n) = Sum_{d | n} f(d) inplace.
    template <typename T, auto add = default_operator::add<T>>
    void zeta(std::vector<T> &f) {
        const int n = f.size();
        std::vector<char> is_prime(n, true);
        auto cum = [&](const int p) {
            for (int q = 1, pq = p; pq < n; ++q, pq += p) {
                f[pq] = add(f[pq], f[q]);
                is_prime[pq] = false;
            }
        };
        for (int p = 2; p < n; ++p) if (is_prime[p]) cum(p);
    }
    // Calculates `f` s.t. g(n) = Sum_{d | n} f(d) inplace.
    template <typename T, auto sub = default_operator::sub<T>>
    void mobius(std::vector<T> &f) {
        const int n = f.size();
        std::vector<char> is_prime(n, true);
        auto diff = [&](const int p) {
            const int qmax = (n - 1) / p, rmax = qmax * p;
            for (int q = qmax, pq = rmax; q >= 1; --q, pq -= p) {
                f[pq] = sub(f[pq], f[q]);
                is_prime[pq] = false;
            }
        };
        for (int p = 2; p < n; ++p) if (is_prime[p]) diff(p);
    }
} // namespace suisen::divisor_transform

#endif // SUISEN_DIVISOR_TRANSFORM

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



#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 6 "library/transform/divisor.hpp"

namespace suisen::divisor_transform {
    // Calculates `g` s.t. g(n) = Sum_{d | n} f(d) inplace.
    template <typename T, auto add = default_operator::add<T>>
    void zeta(std::vector<T> &f) {
        const int n = f.size();
        std::vector<char> is_prime(n, true);
        auto cum = [&](const int p) {
            for (int q = 1, pq = p; pq < n; ++q, pq += p) {
                f[pq] = add(f[pq], f[q]);
                is_prime[pq] = false;
            }
        };
        for (int p = 2; p < n; ++p) if (is_prime[p]) cum(p);
    }
    // Calculates `f` s.t. g(n) = Sum_{d | n} f(d) inplace.
    template <typename T, auto sub = default_operator::sub<T>>
    void mobius(std::vector<T> &f) {
        const int n = f.size();
        std::vector<char> is_prime(n, true);
        auto diff = [&](const int p) {
            const int qmax = (n - 1) / p, rmax = qmax * p;
            for (int q = qmax, pq = rmax; q >= 1; --q, pq -= p) {
                f[pq] = sub(f[pq], f[q]);
                is_prime[pq] = false;
            }
        };
        for (int p = 2; p < n; ++p) if (is_prime[p]) diff(p);
    }
} // namespace suisen::divisor_transform

Back to top page