約数畳み込みを使って最大公約数と集合をうまく扱うメモ

書いて置かないと頭に置いておけない気がしたのでメモを残す. 間違ってたらごめん

これについて気になったので

メビウス関数とかを導入するとより形式的に約数とかを扱えるようになるのかなあ
— Niuez (@xiuez) 2020年1月22日

概要

約数畳み込み
メビウス関数
メビウスの反転公式(約数畳み込みの逆操作)
約数畳み込みと逆約数畳み込みのアルゴリズム $O(A \log{\log A})$
最大公約数の扱い
集合の扱い
AGC038C LCMsの解き方

ネタバレあるので気をつけてください

約数畳み込み

関数 $f(n)$ に対する約数畳み込みとは,

$\begin{eqnarray} g(n) = \sum_{ d | n } f(d) \end{eqnarray}$

あとで解説しますが, 方針としてはこの畳み込んだ後の $g(n)$ を問題を解けるように定義してやることでGCDを綺麗に扱うことができます.

メビウス関数

実際に $g(n)$ を定義してみます. 一番有名なのは $g(n) = \delta(n, 1)$ です. $\delta(n, 1)$ はクロネッカーのデルタです. このとき,

$\begin{eqnarray} g(n) = \delta(n, 1) = \sum_{ d | n } f(d) \end{eqnarray}$

を満たす $f(n)$ はメビウス関数と呼ばれ, $\mu(n)$ と書きます.(メビウス関数 - Wikipedia)

メビウスの反転公式(約数畳み込みの逆操作)

上の式のままだと, $f(n)$ を導くのは困難です. ここで登場するのがメビウスの反転公式です. これは, 約数畳み込みの逆操作に当たります.

$$ \begin{eqnarray} g(n) &=& \sum_{ d | n } f(d) \\ f(n) &=& \sum_{ d | n } g(d) \mu(\frac{n}{d}) \end{eqnarray} $$

これで $g(n)$ を定義してから反転公式を適用することで $f(n)$ を導くことができます.

約数畳み込みと逆約数畳み込みのアルゴリズム

約数畳み込みとその逆はnoshi91さんが計算量 $O(A \log{\log A})$ で計算するアルゴリズムの記事を紹介しています.

noshi91.hatenablog.com

逆約数畳み込みの実装例

template <class T>
void inverse_divisor_transform(vector<T> &a) {
    int n = a.size();
    vector<bool> sieve(n, true);
    for (int p = 2; p < n; ++p) {
        if (sieve[p]) {
            for (int k = (n - 1) / p; k > 0; --k) {
                sieve[k * p] = false;
                a[k * p] -= a[k];
            }
        }
    }
}

最大公約数の扱い

自然数 $n, m$ に対して,

$\begin{eqnarray} \sum_{d | n, d | m} f(d) \end{eqnarray}$

を考えると,

$$ \begin{eqnarray} \sum_{d | n, d | m} f(d) &=& \sum_{d | \gcd(n, m)} f(d) \\ &=& g(\gcd(n, m)) \end{eqnarray} $$

となり, $\gcd(n, m)$ に対する操作ができます. 例えば, $f(n) = \mu(n), g(n) = \delta(n, 1)$ とすると, $g(\gcd(n, m))$ は,「 $n, m$ が互いに素であれば $1$ , そうでなければ $0$ 」となり, 互いに素かどうかの判定ができます.

集合の扱い

例えば, $c_m(d) = [d | m$ ]という関数( $d$ が $m$ を割り切るなら $1$ , そうでなければ $0$ )を考えると,

$$ \sum_{d | n, d | m} f(d) = \sum_{ d | n } f(d) c_m(d) $$

と変形できます.

これを応用します. 自然数の集合 $S$ を考え, $c(d) = \sum_{m \in S} c_m(d)$ とすると,

$$ \begin{eqnarray} \sum_{d | n} f(d) c(d) &=& \sum_{m \in S}\sum_{d | n} f(d) c_m(d) \\ &=& \sum_{m \in S} g(\gcd(n, m)) \end{eqnarray} $$

となります.

$f(n) = \mu(n), g(n) = \delta(n, 1)$ を考えてみると, 「集合 $S$ の中に $n$ と互いに素な要素の数」を計算しています.

AGC038C LCMsを解く

atcoder.jp

ネタバレ回避用

$lcm(x, y) = x (\frac{y}{\gcd(x, y)})$ と変形します. 約数畳み込みを使う方針でやると, この $(\frac{y}{\gcd(x, y)})$ が最後に来てほしい気持ちになります. $g(n) = \frac{1}{n}$ と置くと,

$$ \begin{eqnarray} \frac{y}{\gcd(x, y)} &=& y \cdot g(\gcd(x, y)) \\ &=& \sum_{d | gcd(x, y)} f(d) y \\ &=& \sum_{d | x, d | y} f(d) y \\ &=& \sum_{d | x} f(d) s_y(d) \end{eqnarray} $$

ここで $s_y(d)$ を「 $d$ が $y$ を割り切るなら $y$ , そうでなければ $0$ 」としました.
応用して, 自然数の集合 $S$ を考え, $s(d) = \sum_{m \in S} s_m(d)$ とすると,

$$ \begin{eqnarray} \sum_{d | x} f(d) s(d) &=& \sum_{y \in S}\sum_{d | x} f(d) s_y(d) \\ &=& \sum_{y \in S} y \cdot g(\gcd(x, y)) \end{eqnarray} $$

と計算できて, これに $x$ を掛けると「集合 $S$ の中の各要素と $x$ の最大公約数の和」を計算できました.
計算量は, $O(A \log{\log A} + N \sqrt A)$ です.

atcoder.jp

ソースコード

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
#define rep(i,s,e) for(i64 (i) = (s);(i) < (e);(i)++)

/* modint */
/* IO(niu::fin, niu::fout) */

const i64 MOD = 998244353;
using fp = modint<MOD>;


#include <bits/stdc++.h>
using namespace std;
using i64 = long long;

template <class T>
void inverse_divisor_transform(vector<T> &a) {
    int n = a.size();
    vector<bool> sieve(n, true);
    for (int p = 2; p < n; ++p) {
        if (sieve[p]) {
            for (int k = (n - 1) / p; k > 0; --k) {
                sieve[k * p] = false;
                a[k * p] -= a[k];
            }
        }
    }
}

constexpr i64 A = 1e6 + 1;

int main() {
  std::vector<fp> f(A);
  rep(d,1,A) {
    f[d] = fp(d).pow(MOD - 2);
  }
  inverse_divisor_transform(f);

  i64 N;
  niu::fin >> N;
  vector<fp> sum(A);
  fp ans = 0;
  rep(i,0,N) {
    int x;
    niu::fin >> x;
    fp res = 0;
    for(int d = 1; d * d <= x; d++) {
      if(x % d == 0) {
        res += f[d] * sum[d];
        sum[d] += fp(x);
        if(x / d != d) {
          res += f[x / d] * sum[x / d];
          sum[x / d] += fp(x);
        }
      }
    }
    ans += res * fp(x);
  }
  niu::fout << ans.value() << "\n";
}