倍数変換・約数変換
(multiplicative-function/divisor-multiple-transform.hpp)

View this file on GitHub
Last update: 2023-08-30 23:05:07+09:00
Include: #include "multiplicative-function/divisor-multiple-transform.hpp"

倍数変換/約数変換

(約数/倍数)(ゼータ/メビウス)変換を$\mathrm{O}(N\log \log N)$で行うライブラリ。

概要

このライブラリに出てくる4種の変換を理解すると、約数系包除原理やBurnSideの補題が関係する問題をライブラリパンチできる…かもしれない。

以下に簡潔に内容を解説する。

メビウス関数　$\mu(n)$

メビウス関数とは自然数$n$に対して次のように定義される関数である。

\[\mu(n)=\begin{cases}1 & (n=1)\\ 0 & (nが平方因子を持つ) \\ (-1)^k & (nがk個の素因数を持つ)\end{cases}\]

メビウス関数には以下の基本公式が成り立つ。(包除原理で証明できる。)

\[\sum_{d\mid n}\mu (d)=\begin{cases}1 & (n=1) \\ 0 & (n\neq 1)\end{cases}\]

また、メビウス関数が登場する基本的な公式としてメビウスの反転公式が挙げられる。

メビウスの反転公式

関数$f,g$が全ての正の整数$n$に対して

\[g(n) = \sum_{d\mid n}f(d)\]

を満たしているとする。上の式を言い換えると、$g(n)$は全ての$n$の約数$d$について$f(d)$を足し合わせたものである。競技プログラミングの文脈では「$g$は$f$の約数ゼータ変換である」としばしば表現される。(正式な呼称ではない。)

このとき$f$と$g$の間には次の関係式が成り立つ。(メビウスの反転公式)

\[f(n) = \sum_{d\mid n}\mu(d)g\left(\frac{n}{d}\right)=\sum_{d\mid n}\mu\left(\frac{n}{d}\right)g(d)\]

証明は高校数学の美しい物語に詳しい。また、$g$から$f$への変換を約数メビウス変換と呼ぶ。

もう一つのゼータ/メビウス変換

集合関数に対して上位集合へのゼータ/メビウス変換と下位集合へのゼータ/メビウス変換が定義されるように、倍数ゼータ変換/メビウス変換を定義することが出来る。

関数$f,g$が全ての正の整数$n$に対して

\[g(n) = \sum_{n\mid m}f(m)\]

を満たしているとする。上の式を言い換えると、$g(n)$は全ての$n$の倍数$m$について$f(m)$を足し合わせたものである。競技プログラミングの文脈では「$g$は$f$の倍数ゼータ変換である」としばしば表現される。

このとき$f$と$g$の間には次の関係式が成り立つ。(倍数メビウス変換)

\[f(n) = \sum_{n\mid m}\mu\left(\frac{m}{n}\right)g(m)\]

Burnsideの補題

TODO: 書く

(Burnsideの補題と同じ計算が倍数メビウス変換でもできますねという話を書きたかった　Burnsideの補題は前に頑張って理解したが全て忘却した、悲しいね)

(あと関係ない疑問だけどBurnsideの補題は競プロでは「Polyaの定理」という名前で広まっている？)

関連：yosupoさんのブログ　 yukicoder No.125 悪の花弁　

アルゴリズムの解説および正当性

実際に変換を行うアルゴリズムはkazumaさんのブログが初出のようです。参考

$\mathrm{O}(N\log \log N)$の実装についてはnoshiさんが解説しています。参考

ライブラリの使い方

競技プログラミングで倍数・約数変換を使う場面は主に以下の二つが考えられる。

(1) $f(i)\ (1\leq i\leq n)$が与えられた時に$g(j) \ (1\leq j\leq n)$を求める
(2) $f(i)\ (i\mid n$または$n\mid i)$が与えられた時に$g(n)$を求める

以下のライブラリでは、(1)はvector<T>を、(2)はmap<ll, T>を入力する仕様となっている。

関数一覧

divisor_transform::zeta_transform(a): aを約数ゼータ変換する。
divisor_transform::mobius_transform(a): aを約数メビウス変換する。
multiple_transform::zeta_transform(a): aを倍数ゼータ変換する。
multiple_transform::mobius_transform(a): aを倍数メビウス変換する。

計算量は(1)が$\mathrm{O}(n\log \log n)$、(2)が$\mathrm{O}(\sigma(n)^2)$である。

Depends on

prime/prime-enumerate.hpp

Required by

GCD畳み込み (multiplicative-function/gcd-convolution.hpp)

Verified with

Code

#pragma once

#include <map>
#include <vector>
using namespace std;

#include "../prime/prime-enumerate.hpp"

struct divisor_transform {
  template <typename T>
  static void zeta_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = 1; k * p <= N; ++k) a[k * p] += a[k];
  }
  template <typename T>
  static void mobius_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = N / p; k > 0; --k) a[k * p] -= a[k];
  }

  template <typename I, typename T>
  static void zeta_transform(map<I, T> &a) {
    for (auto p = rbegin(a); p != rend(a); p++)
      for (auto &x : a) {
        if (p->first == x.first) break;
        if (p->first % x.first == 0) p->second += x.second;
      }
  }
  template <typename I, typename T>
  static void mobius_transform(map<I, T> &a) {
    for (auto &x : a) {
      for (auto p = rbegin(a); p != rend(a); p++) {
        if (x.first == p->first) break;
        if (p->first % x.first == 0) p->second -= x.second;
      }
    }
  }
};

struct multiple_transform {
  template <typename T>
  static void zeta_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = N / p; k > 0; --k) a[k] += a[k * p];
  }
  template <typename T>
  static void mobius_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = 1; k * p <= N; ++k) a[k] -= a[k * p];
  }

  template <typename I, typename T>
  static void zeta_transform(map<I, T> &a) {
    for (auto &x : a)
      for (auto p = rbegin(a); p->first != x.first; p++)
        if (p->first % x.first == 0) x.second += p->second;
  }
  template <typename I, typename T>
  static void mobius_transform(map<I, T> &a) {
    for (auto p1 = rbegin(a); p1 != rend(a); p1++)
      for (auto p2 = rbegin(a); p2 != p1; p2++)
        if (p2->first % p1->first == 0) p1->second -= p2->second;
  }
};

/**
 * @brief 倍数変換・約数変換
 * @docs docs/multiplicative-function/divisor-multiple-transform.md
 */

#line 2 "multiplicative-function/divisor-multiple-transform.hpp"

#include <map>
#include <vector>
using namespace std;

#line 2 "prime/prime-enumerate.hpp"

// Prime Sieve {2, 3, 5, 7, 11, 13, 17, ...}
vector<int> prime_enumerate(int N) {
  vector<bool> sieve(N / 3 + 1, 1);
  for (int p = 5, d = 4, i = 1, sqn = sqrt(N); p <= sqn; p += d = 6 - d, i++) {
    if (!sieve[i]) continue;
    for (int q = p * p / 3, r = d * p / 3 + (d * p % 3 == 2), s = 2 * p,
             qe = sieve.size();
         q < qe; q += r = s - r)
      sieve[q] = 0;
  }
  vector<int> ret{2, 3};
  for (int p = 5, d = 4, i = 1; p <= N; p += d = 6 - d, i++)
    if (sieve[i]) ret.push_back(p);
  while (!ret.empty() && ret.back() > N) ret.pop_back();
  return ret;
}
#line 8 "multiplicative-function/divisor-multiple-transform.hpp"

struct divisor_transform {
  template <typename T>
  static void zeta_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = 1; k * p <= N; ++k) a[k * p] += a[k];
  }
  template <typename T>
  static void mobius_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = N / p; k > 0; --k) a[k * p] -= a[k];
  }

  template <typename I, typename T>
  static void zeta_transform(map<I, T> &a) {
    for (auto p = rbegin(a); p != rend(a); p++)
      for (auto &x : a) {
        if (p->first == x.first) break;
        if (p->first % x.first == 0) p->second += x.second;
      }
  }
  template <typename I, typename T>
  static void mobius_transform(map<I, T> &a) {
    for (auto &x : a) {
      for (auto p = rbegin(a); p != rend(a); p++) {
        if (x.first == p->first) break;
        if (p->first % x.first == 0) p->second -= x.second;
      }
    }
  }
};

struct multiple_transform {
  template <typename T>
  static void zeta_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = N / p; k > 0; --k) a[k] += a[k * p];
  }
  template <typename T>
  static void mobius_transform(vector<T> &a) {
    int N = a.size() - 1;
    auto sieve = prime_enumerate(N);
    for (auto &p : sieve)
      for (int k = 1; k * p <= N; ++k) a[k] -= a[k * p];
  }

  template <typename I, typename T>
  static void zeta_transform(map<I, T> &a) {
    for (auto &x : a)
      for (auto p = rbegin(a); p->first != x.first; p++)
        if (p->first % x.first == 0) x.second += p->second;
  }
  template <typename I, typename T>
  static void mobius_transform(map<I, T> &a) {
    for (auto p1 = rbegin(a); p1 != rend(a); p1++)
      for (auto p2 = rbegin(a); p2 != p1; p2++)
        if (p2->first % p1->first == 0) p1->second -= p2->second;
  }
};

/**
 * @brief 倍数変換・約数変換
 * @docs docs/multiplicative-function/divisor-multiple-transform.md
 */

Nyaan's Library

This documentation is automatically generated by online-judge-tools/verification-helper

倍数変換・約数変換
(multiplicative-function/divisor-multiple-transform.hpp)