倍数についてのゼータ変換/メビウス変換
(convolution/multiple-zeta-moebius-transform.hpp)

View this file on GitHub
Last update: 2024-08-21 21:41:37+09:00
Include: #include "convolution/multiple-zeta-moebius-transform.hpp"

TODO: mapで渡すタイプのゼータ/メビウス変換について書く

概要

計算量は $ O(N \log \log N) $ にできるけど、よくわからないので妥協。

倍数についてのゼータ変換

vector<T> zeta_transform_naive<T, op>(const vector<T>& A)

可換な二項演算(大抵の場合、数の和) $\oplus$ があるとき、長さ$N$の数列$A$が与えられたときに、 $ B_i = \bigoplus_{i \mid j} A_j $ なる数列 $B$ を求める。

opは、例えばT=intで二項演算が和なら

int op(int a, int b){
  return a + b;
}

などとして定義された関数を渡す。

計算量

$O(N \log N) $

倍数についてのメビウス変換

可換な二項演算(大抵の場合、数の和) $\oplus$ に逆元があるとき、長さ$N$の数列$B$が与えられたときに、$B_i = \bigoplus_{i \mid j} A_j $ なる数列$A$ を求める。

vector<T> moebius_transform_naive<T, invop>(const vector<T>& B)

invopは、例えばT=intで、二項演算が和なら

int invop(int a, int b) {
  return a - b;
}

などとして定義された関数を渡す。

計算量

$O(N \log N) $

mapでやるやつ

map<I,T> zeta_transform<I,T,op>(const map<I,T>& mp)

map<I,T> moebius_transform<I,T,invop>(const map<I, T>& mp)

どちらもゼータ変換/メビウス変換を行うが、vectorで変換する場合と違い、mapのkeyだけを添え字として扱う。

メビウス関数を使うと、

$ g(n) = \bigoplus_{n \mid i} f(i) \Leftrightarrow f(n) = \bigoplus_{ n \mid i} \mu (\frac i n) g(i) $ と表せるらしい。

計算量

mapのサイズを$N$として

$O(N^2 \log N) $

何が嬉しいのか

数列$A$と$B$が与えられたとして、倍数についてのゼータ変換をすることで $ A^{\prime}$と $B^{\prime}$ を得る。 $ A^{\prime}i B^{\prime}_i = \Sigma{i \mid j , i \mid k} A_j B_k $ であり、

$ {i \mid j , i \mid k} $ とは $ i \mid gcd(j,k) $ と言い換えられるから、$A^{\prime} B^{\prime} $ を倍数についてメビウス変換することで、gcd畳み込みができる、みたいな感じ。(書いてて自信なくなってきたな。間違ってるかも)

Verified with

Code

#ifndef HARUILIB_CONVOLUTION_MULTIPLE_ZETA_MOEBIUS_TRANSFORM_HPP
#define HARUILIB_CONVOLUTION_MULTIPLE_ZETA_MOEBIUS_TRANSFORM_HPP

#include <vector>
#include <map>

namespace multiple {

  // 倍数についてのゼータ変換。 g_n = \Sigma_{n|m} f_m なる g を求める。
  // n|mというのは、m%n==0という意味。
  // O(N log N) (調和級数)
  // うまくやるとO(Nlog(log(N)))にできることがよく知られているが、難しいしlogは定数なので妥協する。
  template <typename T, T (*op)(T, T)>
  std::vector<T> zeta_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = 1; i <= N; i++) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = op(g[i], f[j]);
      }
    }

    return g;
  }

  // 倍数についてのメビウス変換
  // f_n = \Sigma_{n|m} g_m なる g を求める。
  // O(N log N) (調和級数)
  template <typename T, T(*invop)(T, T)>
  std::vector<T> moebius_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = N; i >= 1; i--) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = invop(g[i], g[j]);
      }
    }
    return g;
  }


  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> zeta_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (std::pair<I, T> pit : ret) {
      for (auto p2itr = ret.rbegin(); (*p2itr).first != pit.first; p2itr++) {
        if ((*p2itr).first % pit.first == 0) {
          ret[pit.first] = op(ret[pit.first], (*p2itr).second);
        }
      }
    }

    return ret;
  }


  template <typename I, typename T, T (*invop)(T, T)>
  std::map<I, T> moebius_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (auto p1itr = ret.rbegin(); p1itr != ret.rend(); p1itr++) {
      for (auto p2itr = ret.rbegin(); p2itr != p1itr; p2itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          (*p1itr).second = invop((*p1itr).second, (*p2itr).second);
        }
      }
    }

    return ret;
  }

} // namespace multiple

#endif // HARUILIB_CONVOLUTION_MULTIPLE_ZETA_MOEBIUS_TRANSFORM_HPP

#line 1 "convolution/multiple-zeta-moebius-transform.hpp"



#include <vector>
#include <map>

namespace multiple {

  // 倍数についてのゼータ変換。 g_n = \Sigma_{n|m} f_m なる g を求める。
  // n|mというのは、m%n==0という意味。
  // O(N log N) (調和級数)
  // うまくやるとO(Nlog(log(N)))にできることがよく知られているが、難しいしlogは定数なので妥協する。
  template <typename T, T (*op)(T, T)>
  std::vector<T> zeta_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = 1; i <= N; i++) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = op(g[i], f[j]);
      }
    }

    return g;
  }

  // 倍数についてのメビウス変換
  // f_n = \Sigma_{n|m} g_m なる g を求める。
  // O(N log N) (調和級数)
  template <typename T, T(*invop)(T, T)>
  std::vector<T> moebius_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = N; i >= 1; i--) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = invop(g[i], g[j]);
      }
    }
    return g;
  }


  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> zeta_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (std::pair<I, T> pit : ret) {
      for (auto p2itr = ret.rbegin(); (*p2itr).first != pit.first; p2itr++) {
        if ((*p2itr).first % pit.first == 0) {
          ret[pit.first] = op(ret[pit.first], (*p2itr).second);
        }
      }
    }

    return ret;
  }


  template <typename I, typename T, T (*invop)(T, T)>
  std::map<I, T> moebius_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (auto p1itr = ret.rbegin(); p1itr != ret.rend(); p1itr++) {
      for (auto p2itr = ret.rbegin(); p2itr != p1itr; p2itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          (*p1itr).second = invop((*p1itr).second, (*p2itr).second);
        }
      }
    }

    return ret;
  }

} // namespace multiple

倍数についてのゼータ変換/メビウス変換 (convolution/multiple-zeta-moebius-transform.hpp)

概要