Nyaan's Library

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

View on GitHub

:heavy_check_mark: 有名な乗法的関数
(multiplicative-function/mf-famous-series.hpp)

有名な乗法的関数

来るべきDGF(ディリクレ母関数)の流行に備えて(本当に流行るのか？)、自分が理解していることをメモしておく…

乗法的関数

$f(n)$が任意の$\gcd(a,b) = 1$である自然数$a,b$に対して$f(ab) = f(a)f(b)$となる時、$f(n)$は乗法的であると呼ぶ。特に任意の$a,b$について$f(ab)=f(a)f(b)$が成り立つ時は完全乗法的であると呼ぶ。

乗法的関数の重要な性質は以下のようなものが挙げられる。

また、乗法的関数に関するアルゴリズムは以下のものが知られている。

有名な乗法的関数

このうち下の3つは次の畳み込みの関係が知られている。

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

メビウス関数に関する説明はここに詳しく書いた。

競技プログラミングにおいてのメビウス関数の利用法は(雑に説明すると)包除原理の$(-1)^k$と似た使い方をすることが多い。実際に具体例を見てみる。

例題1

約数系包除を使えば$\mathrm{O}(\sigma(N)^2)$で解けるがより高速な解法を考えたい。 まずは反転公式を使わずに考察してみる。具体的な$N$についていくつか実験してみると、

\[f(16)=g(16)-g(8)\] \[f(12)=g(12)-g(6)-g(4)+g(2)\] \[f(30)=g(30)-g(15)-g(10)-g(6)+g(5)+g(3)+g(2)-g(1)\]

のようになり、これは包除原理+BIT全探索で解くことが出来る。(tsutajさんの非常にわかりやすい包除PDFに類題の詳しい説明がある。)

一方、反転公式を使うと、例えば$n=12$の時は

\[f(12) = \sum_{d \in \lbrace 1,2,3,4,6,12\rbrace}\mu\left(\frac{n}{d}\right)g(d)=g(12)-g(6)-g(4)+g(2)\]

となり包除原理によって得られる結果と一致する。

下位集合に対するゼータ変換が包除原理で解けることと、約数集合に対するゼータ変換がメビウス関数で解けることは同じような関係にある(？)と言える。

例題2：Cube-loving Numbers　(HackerRank)

$1\leq n\leq N$において自然数$a,b$を用いて$n=a^3\times b$と表せる$n$の個数$g(a)$は$g(a)=\lfloor\frac{N}{a^3}\rfloor$と容易に表せるので、この式をうまく利用して答えを求めたい。(直感的には、$g(n)$は一つの自然数を複数回カウントする関数なのでメビウス変換したいという気持ちになる。)

対象を重複なく数え上げるために、自然数$n$に一対一対応する$(a,b)$を決定したい。具体的には、「$n=A^3\times B$を満たす自然数の組$(A,B)$の中で最も$B$が小さい組」を$n$に対応する組$(a,b)$とおく。そして、$f(a)$を「$(a,\frac{n}{a^3})$と対応している$N$以下の自然数$n$の個数」とおく。

$f$と$g$の関係式を得るために、$g(a)$でカウントされている自然数$n$が$f$のどこでカウントされているかを考える。$n=a^3\times b$としたとき、$b$に一対一対応する整数の組を$(c,d)$とおくと、$n$に対応する組は$(ca,d)$であるから$f(ca)$で数え上げられていることが分かる。逆に$f(ca)$で数え上げられた$n$が$g(a)$で数えられていることも示せる。よって$f(a)$と$g(a)$の間には

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

という倍数変換の関係式が成り立つことがわかる。また、求める答えは$M=\lfloor\sqrt[3]{N}\rfloor$と置いたとき$\sum_{a=2}^Mf(a)$である。 ($M\lt a$のとき$f(a)=g(a)=0$である事実を利用している。)

よって倍数メビウス変換を用いれば$\mathrm{O}(M\log \log M)$で計算できることが示せたが、メビウス関数を用いることでさらなる高速化を図りたい。$\sum_{a=2}^Mf(a)$を$g(m)$の線形和に分解したときの$g(m)$の寄与を考察すると、

\[\sum_{a=2}^M f(a)=\sum_{2\leq a\leq M, a\mid m} \left( \mu\left(\frac{m}{a}\right)g(m)\right)\] \[=\sum_{2\leq m\leq M} g(m)\left(\sum_{a\mid m,a\neq 1}\mu\left(\frac{m}{a}\right)\right)\] \[=\sum_{2\leq m\leq M} g(m)(-\mu(m)+\sum_{a\mid m}\mu(a))=-\sum_{2\leq m\leq M} g(m)\mu(m)\]

と非常にきれいな式になる。$\mu(m)$および$g(m)$は線形で列挙できるため、求める答えも線形で列挙できる。以上よりこの問題を$\mathrm{O}(\sqrt[3]{N})$で解くことが出来た。

Depends on

Required by

Verified with

Code

#pragma once

#include "enumerate-multiplicative-function.hpp"

namespace multiplicative_function {
template <typename T>
T moebius(int, int, int c) {
  return c == 0 ? 1 : c == 1 ? -1 : 0;
}
template <typename T>
T sigma0(int, int, int c) {
  return c + 1;
}
template <typename T>
T sigma1(int n, int p, int) {
  return (n - 1) / (p - 1) + n;
}
template <typename T>
T totient(int n, int p, int) {
  return n - n / p;
}
}  // namespace multiplicative_function

template <typename T>
static constexpr vector<T> mobius_function(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::moebius<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> sigma0(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::sigma0<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> sigma1(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::sigma1<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> totient(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::totient<T>> em(
      n);
  return em.run();
}

/**
 * @brief 有名な乗法的関数
 */
#line 2 "multiplicative-function/mf-famous-series.hpp"

#line 2 "multiplicative-function/enumerate-multiplicative-function.hpp"

#include <vector>
using namespace std;

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

#include <cmath>
#line 5 "prime/prime-enumerate.hpp"
using namespace std;

// 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 7 "multiplicative-function/enumerate-multiplicative-function.hpp"

// f(n, p, c) : n = pow(p, c), f is multiplicative function

template <typename T, T (*f)(int, int, int)>
struct enumerate_multiplicative_function {
  enumerate_multiplicative_function(int _n)
      : ps(prime_enumerate(_n)), a(_n + 1, T()), n(_n), p(ps.size()) {}

  vector<T> run() {
    a[1] = 1;
    dfs(-1, 1, 1);
    return a;
  }

 private:
  vector<int> ps;
  vector<T> a;
  int n, p;
  void dfs(int i, long long x, T y) {
    a[x] = y;
    if (y == T()) return;
    for (int j = i + 1; j < p; j++) {
      long long nx = x * ps[j];
      long long dx = ps[j];
      if (nx > n) break;
      for (int c = 1; nx <= n; nx *= ps[j], dx *= ps[j], ++c) {
        dfs(j, nx, y * f(dx, ps[j], c));
      }
    }
  }
};

template <typename T, T (*f)(int, int, int)>
using enamurate_multiplicative_function =
    enumerate_multiplicative_function<T, f>;

/**
 * @brief 乗法的関数の列挙
 */
#line 4 "multiplicative-function/mf-famous-series.hpp"

namespace multiplicative_function {
template <typename T>
T moebius(int, int, int c) {
  return c == 0 ? 1 : c == 1 ? -1 : 0;
}
template <typename T>
T sigma0(int, int, int c) {
  return c + 1;
}
template <typename T>
T sigma1(int n, int p, int) {
  return (n - 1) / (p - 1) + n;
}
template <typename T>
T totient(int n, int p, int) {
  return n - n / p;
}
}  // namespace multiplicative_function

template <typename T>
static constexpr vector<T> mobius_function(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::moebius<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> sigma0(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::sigma0<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> sigma1(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::sigma1<T>> em(
      n);
  return em.run();
}

template <typename T>
static constexpr vector<T> totient(int n) {
  enumerate_multiplicative_function<T, multiplicative_function::totient<T>> em(
      n);
  return em.run();
}

/**
 * @brief 有名な乗法的関数
 */
Back to top page