とても久しぶりです！
1 年ぶりの投稿となりました、大槻 (通称、けんちょん) です。

去年、『AtCoder 版！マスター・オブ・整数』と題して、プログラミングコンテストで出題される整数問題を解くときに有効な考え方を特集する記事を 2 本書きました！

今回はその続編として、素数を列挙するアルゴリズムであるエラトステネスの篩を特集していきます。なお今回の記事の内容は、競プロへの応用を意識していますが、純粋に数学的興味に沿って読み進めることもできるものになっています。下図は、これから紹介するエラトステネスの篩のイメージ図です。

0. はじめに

エラトステネスの篩は、$1$ 以上 $N$ 以下の素数をすべて列挙する方法です。たとえば $20$ 以下の素数を列挙すると、$2, 3, 5, 7, 11, 13, 17, 19$ となります。

さて、エラトステネスの篩はなんのために必要なのでしょうか？
素数を列挙するのは、$1$ から $N$ までの整数を順に「素数かどうかを判定していく」という愚直な方法でもできそうです。整数 $k$ が素数かどうかの判定については「$2$ 以上 $\sqrt{k}$ 以下の整数の中に $k$ の約数が存在するかどうかを調べる」という方法でできます (詳しくはこの記事で解説しています)。

しかしこのような愚直な方法では、素数判定に $O(\sqrt{N})$ の計算量がかかりますので、全体として $O(N \sqrt{N})$ の計算量となります。

一方、エラトステネスの篩を用いると $O(N \log\log N)$ の計算量で $N$ 以下の素数を列挙できます！！！¹

そして、エラトステネスの篩の効能はこれだけではありません。エラトステネスの篩を用いると、次のようなことが高速にできます！！！！！

$1$ 以上 $N$ 以下の各整数の高速素因数分解²
$1$ 以上 $N$ 以下の各整数の高速約数列挙
$1$ 以上 $N$ 以下の各整数のメビウス関数 (後述) の値を高速に求めること
$1$ 以上 $N$ 以下の各整数のオイラー関数 (後述) の値を高速に求めること
$1$ 以上 $N$ 以下の整数 $n$ に対する関数 $f(n)$ が与えられたときに、$F(k) = \sum_{k|i}f(i)$ によって定義される関数 $F(n)$ を求めること (高速ゼータ変換) ³
$1$ 以上 $N$ 以下の整数 $n$ に対する関数 $F(n)$ が与えられたときに、$F(k) = \sum_{k|i}f(i)$ が成立するような関数 $f(n)$ を求めること (高速メビウス変換)
$1$ 以上 $N$ 以下の整数 $n$ に対する関数 $g(n), h(n)$ が与えられたときに、$f(k) = \sum_{{\rm GCD}(i, j) = k} g(i)h(j)$ で定義される関数 $f(n)$ を求めること (添字 GCD 畳み込み)

ここで、$k|i$ という記法は「$i$ が $k$ の倍数であること」を意味します。さて、1 番目の高速素因数分解は、意外と使える場面が多いです。$1$ 以上 $N$ 以下の各整数を愚直に素因数分解する方法よりも高速です。7 番目の (約数系) 高速メビウス変換は、俗に、約数系包除原理と呼ばれるものでもあります。AtCoder でも、約数系包除原理を用いることで解ける問題は数多く出題されています。

今回は、エラトステネスの篩を理解したうえで、以上の応用的な使用例をまとめて整理していきます。後半は少し難しい話になっていきます。興味関心に応じて読み進めていただけたら嬉しいです。

1. エラトステネスの篩とは

まずはエラトステネスの篩がどんなものなのかを見ていきましょう。$1$ 以上 $N$ 以下の素数をすべて列挙したいのですが、逆に「合成数が残っている限り、それらをすべてふるい落としていく」という方法でやっていきます。具体的には、

まだ生き残っている整数で、しかも "素数ラベル" が貼られていないような最小の整数を $p$ とする (このとき、$p$ は素数になります)
$p$ に "素数ラベル" を貼る
$p$ 以外の $p$ の倍数をすべてふるい落とす

というのを繰り返していきます。最後に生き残ったものが素数です。具体例として、$N = 50$ の場合を考えてみます。下図のように $2$ から $50$ までの整数を並べておきます ($1$ は素数ではないので最初から除外しておきます)。

まず、$2$ に素数ラベルをつけて、それ以外の $2$ の倍数をふるい落としましょう。下図のようになります。

ここで生き残っている整数のうち、素数ラベルのついていない最小の整数は $3$ です。そこで $3$ に素数ラベルをつけて $3$ の倍数をふるい落とします。

ここで生き残っている整数のうち、素数ラベルのついていない最小の整数は $5$ です。そこで $5$ に素数ラベルをつけて $5$ の倍数をふるい落とします。

ここで生き残っている整数のうち、素数ラベルのついていない最小の整数は $7$ です。そこで $7$ に素数ラベルをつけて $7$ の倍数をふるい落とします。

...と、この調子で最後までやると、下図のようになります。

合成数が次々とふるい落とされていき、最後に素数だけが残っている様子が見て取れます。なお余談ですが、実は $7$ の倍数をふるい落とした時点で、合成数が残っていないことに気づかれた方も多いかもしれません。それもそのはずで、$7$ の次の素数は $11$ ですが、$11 \times 11 = 121$ 以下の整数で合成数なものはすべて $2, 3, 5, 7$ のいずれかで割り切れるからです。

2. エラトステネスの篩の実装

それでは、エラトステネスの篩を実装してみましょう。まずサイズ $N + 1$ の配列 isprime を用意します⁴。

isprime[n] ← 整数 $n$ が素数ならば True、そうでないならば False

これを用いて次のように実装できます。ここでは isprime 全体を True で初期化しておいて ($2$ 以上 $N$ 以下の整数すべてに素数ラベルをつけておいて)、合成数は False に変えていく (素数ラベルを剥奪する) という実装をしています。

エラトステネスの篩の実装

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

// 1 以上 N 以下の整数が素数かどうかを返す
vector<bool> Eratosthenes(int N) {
    // テーブル
    vector<bool> isprime(N+1, true);

    // 1 は予めふるい落としておく
    isprime[1] = false;

    // ふるい
    for (int p = 2; p <= N; ++p) {
        // すでに合成数であるものはスキップする
        if (!isprime[p]) continue;

        // p 以外の p の倍数から素数ラベルを剥奪
        for (int q = p * 2; q <= N; q += p) {
            isprime[q] = false;
        }
    }

    // 1 以上 N 以下の整数が素数かどうか
    return isprime;
}

int main() {
    vector<bool> isprime = Eratosthenes(50);

    for (int v = 2; v <= 50; ++v) {
        cout << v << ": "
             << (isprime[v] ? "prime" : "not") << endl;
    }
}

なおこのコードは、$N = 50$ として、$N$ 以下の整数を素数かどうかを判定するものとなっています。このコードの実行結果は次のようになります。

エラトステネスの篩の実行結果

2: prime
3: prime
4: not
5: prime
6: not
7: prime
8: not

...(中略)...

44: not
45: not
46: not
47: prime
48: not
49: not
50: not

3. エラトステネスの篩の計算量

エラトステネスの篩の計算量が $O(N \log\log N)$ になることを見ていきましょう。

3-1. 劣化エラトステネスの篩の計算量 (調和級数)

まず本物のエラトステネスの篩の計算量を考える前に、下の "劣化エラトステネスの篩" のコードの計算量を考えてみましょう。

これは、エラトステネスの篩において、if (!isprime[p]) continue; を省略したものです。つまり、

$2$ 以外の $2$ の倍数をふるい落とす
$3$ 以外の $3$ の倍数をふるい落とす
$4$ 以外の $4$ の倍数をふるい落とす
$5$ 以外の $5$ の倍数をふるい落とす
$6$ 以外の $6$ の倍数をふるい落とす
...
$N$ 以外の $N$ の倍数をふるい落とす

という挙動になります。本物のエラトステネスの篩は、素数以外の $p$ に対して $p$ の倍数をふるい落とすようなことはしません。たとえば $6$ の倍数をふるい落とそうとするのは明らかに無駄です。なぜなら、$6$ の倍数は $2$ や $3$ の倍数でもあるため、前の時点ですでにふるい落とされているからです。

しかしそれでも、この劣化エラトステネスの篩でも $1$ 以上 $N$ 以下の素数を正しく列挙できますし、計算量も $O(N \log N)$ となります。本物のエラトステネスの篩の $O(N \log\log N)$ には劣りますが、十分高速です。下で、劣化エラトステネスの篩の計算量が $O(N \log N)$ になることを示します。

劣化エラトステネスの篩

vector<bool> Eratosthenes(int N) {
    // テーブル
    vector<bool> isprime(N+1, true);

    // 1 は予めふるい落としておく
    isprime[1] = false;

    // ふるい
    for (int p = 2; p <= N; ++p) {
        // このスキップをやめてみる
        //if (!isprime[p]) continue;

        // p 以外の p の倍数から素数ラベルを剥奪
        for (int q = p * 2; q <= N; q += p) {
            isprime[q] = false;
        }
    }

    // 1 以上 N 以下の整数が素数かどうか
    return isprime;
}

さて、劣化エラトステネスの篩の計算量が $O(N \log N)$ になることを示してみましょう。各 $p$ に対して、2 番目の添字 $q$ は

$q = 2p, 3p, 4p, \dots$

を動きます。よって添字 $q$ の取りうる値の「個数」は $\frac{N}{p}$ に比例しますので、劣化エラトステネスの計算時間 $T(N)$ は、細かい係数を無視して考えると

$T(N) = \frac{N}{2} + \frac{N}{3} + \dots + \frac{N}{N} = N(\frac{1}{2} + \frac{1}{3} + \dots + \frac{1}{N})$

と評価できます。さて、漸近的に

$1 + \frac{1}{2} + \dots + \frac{1}{N} \sim \log N$

となることは有名です ($f(x) = \frac{1}{x}$ の積分が $\log x$ であることから導けます⁵)。このことから

$T(N) \sim N \log N$

と評価できます。以上から、劣化エラトステネスの篩の計算量が $O(N \log N)$ になることが分かりました。

なお、劣化エラトステネスの篩に限らず、下のような二重 for 文を回す計算量は $O(N \log N)$ になります。このことを活用する AtCoder の問題も多数出題されています。いずれも緑色 diff の問題です。

計算量が調和級数になるループ

for (int p = 1; p <= N; ++p) {
    for (int q = p; q <= N; q += p) {

    }
}

3-2. エラトステネスの篩の計算量

劣化エラトステネスの篩の計算時間 $T(N)$ は

$T(N) = N(\frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \frac{1}{5} + \frac{1}{6} + \dots)$

という感じでした。これに対して本物のエラトステネスの篩の計算時間は $P(N)$ は

$P(N) = N(\frac{1}{2} + \frac{1}{3} + \frac{1}{5} + \frac{1}{7} + \frac{1}{11} + \dots)$

というふうになります。素数の逆数和に $N$ をかけたものとなっています。このようになるのは、素数以外の $p$ に対しては「$p$ の倍数をふるい落とす」という処理をスキップするからです。素数の逆数和については、漸近的に $\log\log N$ 程度のスピードで発散することが知られています。より正確には、定数 $b$ を用いて

$\displaystyle \sum_{p \le N, p は素数} \frac{1}{p} = \log\log N + b + o(1)$

という評価が成り立つことが知られています。この定数 $b$ は、Meissel–Mertens 定数と呼ばれています。

以上から、エラトステネスの篩の計算量が $O(N \log\log N)$ となることがわかりました。

3-3. より高速な篩

余談ですが、$1$ 以上 $N$ 以下の素数を列挙するアルゴリズムとして、エラトステネスの篩よりも高速な篩もあります。

アトキンの篩：$O(\frac{N}{\log\log N})$
線形篩：$O(N)$

速度的にはアトキンの篩の方が高速ですが、線形篩は 4 章で解説する「高速素因数分解」にも適用できるメリットがあります。

4. 活用例 (1): 高速素因数分解・高速約数列挙

ここから先は、エラトステネスの篩のさまざまな活用例を見ていくことにしましょう。

4-1. 高速素因数分解

まずは、エラトステネスの篩を活用して、$1$ 以上 $N$ 以下の整数をすべて素因数分解する処理を高速に実施してみましょう。これも

愚直な方法 (それぞれ素因数分解)：$O(N \sqrt{N})$
エラトステネスの篩を活用：$O(N \log N)$

というように、愚直な方法よりずっと高速な方法となっています⁶。より正確には、エラトステネスの篩を用いた前処理 ($O(N \log\log N)$) のあと、$1$ 以上 $N$ 以下の各整数 $i$ をそれぞれ $O(\log N)$ で素因数分解できるということです。アイデアも単純です。エラトステネスの篩を実行中に、次の配列 minfactor の値を同時に求めていきます。

minfactor[n] ← 整数 $n$ を割り切る最小の素数

たとえば minfactor[12] = 2, minfactor[49] = 7 というふうになります。次のコードのように、エラトステネスの篩に少し手を加えるだけで実現できます。

minfactorを求める

vector<bool> Eratosthenes(int N) {
    // テーブル
    vector<bool> isprime(N+1, true);
    
    // 整数 i を割り切る最小の素数
    vector<int> minfactor(N+1, -1);
    
    // 1 は予めふるい落としておく
    isprime[1] = false;
    minfactor[1] = 1;

    // ふるい
    for (int p = 2; p <= N; ++p) {
        // すでに合成数であるものはスキップする
        if (!isprime[p]) continue;

        // p についての情報更新
        minfactor[p] = p;

        // p 以外の p の倍数から素数ラベルを剥奪
        for (int q = p * 2; q <= N; q += p) {
            // q は合成数なのでふるい落とす
            isprime[q] = false;

            // q は p で割り切れる旨を更新
            if (minfactor[q] == -1) minfactor[q] = p;
        }
    }
}

ここまでは $O(N \log\log N)$ でできます。次に、$1$ 以上 $N$ 以下の整数 $n$ を $O(\log N)$ で素因数分解する方法を考えてみましょう。まず、$p$ = minfactor[n] とすると、$n$ は素数 $p$ を素因数にもつことがわかっていますので、$n$ を $p$ で割れるだけ割ります。そうしてできる整数を $n'$ としましょう。このとき、ふたたび $p'$ = minfactor[n'] とすると $n'$ は素数 $p'$ を素因数にもつことがわかりますので、$n'$ を $p'$ で割れるだけ割ります。

以上の処理を、$n' = 1$ となるまで繰り返していくことで素因数分解できます。具体的な実装を下のコードに示します。なおここでは、エラトステネスの篩全体を構造体として表現してあげて、整数 $n$ を引数として受け取って素因数分解するメンバ関数 factorize(int n) を実装する形をとっています

最後に、関数 factorize() の計算量を評価します。$n$ の素因数分解の結果を

$n = p_1^{e_1} p_2^{e_2} \dots p_K^{e^K}　(p_1 < p_2 < \dots < p_K)$

と表すと、割り算の回数は $e_1 + e_2 + \dots + e_K$ となります。これは少なくとも $\log_{p_1} N$ 以下の値となります。よって計算量は $O(\log N)$ と評価できます。

高速素因数分解メソッドを備えたエラトステネス

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

// エラトステネスの篩
struct Eratosthenes {
    // テーブル
    vector<bool> isprime;
    
    // 整数 i を割り切る最小の素数
    vector<int> minfactor;

    // コンストラクタで篩を回す
    Eratosthenes(int N) : isprime(N+1, true),
                          minfactor(N+1, -1) {
        // 1 は予めふるい落としておく
        isprime[1] = false;
        minfactor[1] = 1;

        // 篩
        for (int p = 2; p <= N; ++p) {
            // すでに合成数であるものはスキップする
            if (!isprime[p]) continue;

            // p についての情報更新
            minfactor[p] = p;
            
            // p 以外の p の倍数から素数ラベルを剥奪
            for (int q = p * 2; q <= N; q += p) {
                // q は合成数なのでふるい落とす
                isprime[q] = false;
                
                // q は p で割り切れる旨を更新
                if (minfactor[q] == -1) minfactor[q] = p;
            }
        }
    }

    // 高速素因数分解
    // pair (素因子, 指数) の vector を返す
    vector<pair<int,int>> factorize(int n) {
        vector<pair<int,int>> res;
        while (n > 1) {
            int p = minfactor[n];
            int exp = 0;

            // n で割り切れる限り割る
            while (minfactor[n] == p) {
                n /= p;
                ++exp;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }  
};

int main() {
    // エラトステネスの篩
    Eratosthenes er(50);

    // 結果表示
    for (int n = 2; n <= 50; ++n) {
        auto pf = er.factorize(n);
        cout << n << ": ";
        for (int i = 0; i < pf.size(); ++i) {
            if (i > 0) cout << " * ";
            cout << pf[i].first << "^" << pf[i].second;
        }
        cout << endl;
    }
}

上記コードの結果は次のようになります。

2: 2^1
3: 3^1
4: 2^2
5: 5^1
6: 2^1 * 3^1
7: 7^1

...(中略)...

45: 3^2 * 5^1
46: 2^1 * 23^1
47: 47^1
48: 2^4 * 3^1
49: 7^2
50: 2^1 * 5^2

4-2. 高速約数列挙

整数 $n$ の素因数分解と同様、$n$ の約数列挙も、配列 minfactor を求めておくと高速に実現できます。上のエラトステネス構造体の新たなメンバ関数 divisors(int n) として実装しておきます。

計算量は $n$ の約数の個数を $\sigma(n)$ として、$O(\sigma(n))$ となります。$n \le 10^9$ の範囲では、$\sigma(n) \le 1344$ ($n = 735134400$ で最大) であることが知られていますので、十分高速です。

高速約数列挙

struct Eratosthenes {
    // テーブル
    vector<bool> isprime;
    
    // 整数 i を割り切る最小の素数
    vector<int> minfactor;

    // コンストラクタと factorize は省略

    // 高速約数列挙
    vector<int> divisors(int n) {
        vector<int> res({1});

        // n を素因数分解 (メンバ関数使用)
        auto pf = factorize(n);

        // 約数列挙
        for (auto p : pf) {
            int s = (int)res.size();
            for (int i = 0; i < s; ++i) {
                int v = 1;
                for (int j = 0; j < p.second; ++j) {
                    v *= p.first;
                    res.push_back(res[i] * v);
                }
            }
        }
        return res;
    }
};

5. 活用例 (2): メビウス関数値の列挙

メビウス関数 $\mu(n)$ とは、正の整数 $n$ に対して $-1, 0, 1$ のいずれかの値を返す関数であり、次のように定義されます。

$\mu(1) = 1$
$n$ がある素数 $p$ で $2$ 回以上割り切れるとき、$\mu(n) = 0$
$n = p_1 \times p_2 \times \dots p_K$ と素因数分解できるとき、$\mu(n) = (-1)^K$

小さな $n$ に対するメビウス関数 $\mu$ の値は次のようになります。

$\mu(1) = 1$
$\mu(2) = -1$
$\mu(3) = -1$
$\mu(4) = 0　(4 = 2^2)$
$\mu(5) = -1$
$\mu(6) = 1　(6 = 2 \times 3)$
$\mu(7) = -1$
$\mu(8) = 0　(8 = 2^3)$
$\mu(9) = 0　(9 = 3^2)$
$\mu(10) = 1　(10 = 2 \times 5)$
$\mu(11) = -1$
$\mu(12) = 0　(12 = 2^2 \times 3)$
$\mu(13) = -1$
$\mu(14) = 1　(14 = 2 \times 7)$
$\mu(15) = 1　(15 = 3 \times 5)$
$\mu(16) = 0　(16 = 2^4)$
...
$\mu(30) = -1　(30 = 2 \times 3 \times 5)$
...

つまり、素因数分解したときの指数に $2$ 以上となる箇所がある場合は $0$ となり、指数がすべて $1$ である場合には素数の個数の偶奇によって $\pm 1$ のいずれかとなります。

さて、$1$ 以上 $N$ 以下の各整数 $n$ に対するメビウス関数の値も、エラトステネスの篩を回すのと同時に求めることができます。次の配列 mobius を用意しましょう。

mobius[n] ← 整数 $n$ に対するメビウス関数の値 $\mu(n)$

エラトステネスの篩のコンストラクタを下のコードのように変更することで、配列 mobius が求められます。素数 $p$ の倍数 $n$ について篩にかけるとき、

$n$ が $p$ で 2 回以上割り切れるならば、mobius[n] = 0
そうでないならば、mobius[n] = -mobius[n]

というようにします。ところで、メビウス関数なんて何に使うんだろうと思った方も多いかもしれません。後述するように、約数系包除原理で大活躍します。

エラトステネスの篩でメビウス関数を求める

struct Eratosthenes {
    // テーブル
    vector<bool> isprime;
    
    // 整数 n を割り切る最小の素数
    vector<int> minfactor;

    // メビウス関数値
    vector<int> mobius;

    // コンストラクタで篩を回す
    Eratosthenes(int N) : isprime(N+1, true),
                          minfactor(N+1, -1),
                          mobius(N+1, 1) {
        // 1 は予めふるい落としておく
        isprime[1] = false;
        minfactor[1] = 1;

        // 篩
        for (int p = 2; p <= N; ++p) {
            // すでに合成数であるものはスキップする
            if (!isprime[p]) continue;

            // p についての情報更新
            minfactor[p] = p;
            mobius[p] = -1;
            
            // p 以外の p の倍数から素数ラベルを剥奪
            for (int q = p * 2; q <= N; q += p) {
                // q は合成数なのでふるい落とす
                isprime[q] = false;
                
                // q は p で割り切れる旨を更新
                if (minfactor[q] == -1) minfactor[q] = p;
                if ((q / p) % p == 0) mobius[q] = 0;
                else mobius[q] = -mobius[q];
            }
        }
    }
    
    // 高速素因数分解、高速約数列挙は省略
};

6. 活用例 (3): 高速ゼータ変換

いよいよ、エラトステネスの篩の本格的な活用例を見てみましょう。正の整数 $n$ に対する関数 $f(n)$ があって、$f(1), f(2), \dots, f(N)$ が与えられている状況を考えます。ただし $n > N$ のとき $f(n) = 0$ であるものと仮定します。このとき、

$F(n) = \sum_{n | i} f(i)$

で定義される関数 $F(n)$ を考えます⁷。定義から $n > N$ のときは $F(n) = 0$ となります。$F(1), F(2), \dots, F(N)$ の値を高速に求める処理は高速ゼータ変換と呼ばれます。なお、ここでいう高速ゼータ変換は、正確には約数系の高速ゼータ変換と呼ぶべきものです。単に高速ゼータ変換といった場合には、有限集合 $V$ の部分集合の集合 $2^{V}$ 上で定義された関数 $f$ に対して、$F(S) = \sum_{S \in T} f(T)$ と定義される関数 $F$ を求めるアルゴリズムを指すのが普通です⁸。

さて、上記の関数 $F$ は、$O(N \log\log N)$ の計算量で求められます。たとえば $N = 12$ のとき、次のようになります。

$F(1) = f(1) + f(2) + f(3) + \dots + f(12)$
$F(2) = f(2) + f(4) + f(6) + f(8) + f(10) + f(12)$
$F(3) = f(3) + f(6) + f(9) + f(12)$
$F(4) = f(4) + f(8) + f(12)$
$F(5) = f(5) + f(10)$
$F(6) = f(6) + f(12)$
$F(7) = f(7)$
$F(8) = f(8)$
$F(9) = f(9)$
$F(10) = f(10)$
$F(11) = f(11)$
$F(12) = f(12)$

これらの表式は何を表しているのでしょうか。まず一般に $N$ 以下の素数をすべて列挙して $(p_1, p_2, \dots, p_K)$ としてみましょう。今回は $N = 12$ ですので、$(2, 3, 5, 7, 11)$　となります。そして、各整数 $n$ に対して

$n = 2^{e_1} \times 3^{e_2} \times 5^{e_3} \times 7^{e_4} \times 11^{e_5}$

と表したときの指数からなるベクトル $(e_1, e_2, e_3, e_4, e_5)$ を考えます。そして表式 $F(n) = \sum_{n | i} f(i)$ について、$i$ を同じく指数ベクトル $(d_1, d_2, d_3, d_4, d_5)$ のように表したとしましょう。このとき $n | i$ であることは、$d_1 \ge e_1$, $d_2 \ge e_2$, $d_3 \ge e_3$, $d_4 \ge e_4$, $d_5 \ge e_5$ が成り立つことと同値です。よって $\sum_{n | i}$ という表式は、多次元空間上で超直方体領域の総和 (いわゆる多次元累積和) をとっていることを意味します。なお、累積和については、次の記事で解説しています。

累積和を何も考えずに書けるようにする！

簡単のため、$n$ が $2, 3$ 以外の素因数をもつときには $f(n) = 0$ であることとして、二次元の指数ベクトル $(e_1, e_2)$ で考えてみましょう。このとき、$F(1)$, $F(2)$, $F(3)$, $F(4)$, $F(6)$, $F(12)$ はそれぞれ下図のように表せます。

これらは右上を始点とした二次元累積和をとっているものとみなせます。一般の場合には、多次元累積和をとります。具体的には下のコードのように実装できます。各素数 $p$ に対して、座標の大きい方 (座標値 k * p ) から座標の小さい方 (座標値 k ) へと、足し込むように処理しています。

計算量を評価してみましょう。前処理としてのエラトステネスの篩自体の計算量も $O(N \log\log N)$ ですが、高速ゼータ変換のループに要する計算量も、エラトステネスの篩と全く同様であることがわかります。よって高速ゼータ変換の計算量は $O(N \log\log N)$ となります。

高速ゼータ変換

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

// エラトステネスの篩
vector<bool> Eratosthenes(int N) {
    vector<bool> isprime(N+1, true);
    isprime[1] = false;
    for (int p = 2; p <= N; ++p) {
        if (!isprime[p]) continue;
        for (int q = p * 2; q <= N; q += p) {
            isprime[q] = false;
        }
    }
    return isprime;
}

// 高速ゼータ変換
// 入力 f が in-place に更新されて、F になる
template<class T> void fast_zeta(vector<T> &f) {
    int N = f.size();
    
    // エラトステネスの篩を用いて素数を列挙
    vector<bool> isprime = Eratosthenes(N);
    
    // 各素数 p 軸に対して
    // 大きい座標 (k * p) から小さい座標 (k) へと足し込む
    for (int p = 2; p < N; ++p) {
        if (!isprime[p]) continue;
        
        // 座標が大きい方を起点として累積和をとる
        for (int k = (N - 1) / p; k >= 1; --k) {
            f[k] += f[k * p];
        }
    }
}

int main() {
    // f[0] は関係ないので適当な値を入れておく
    vector<int> f = {-1000, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
    fast_zeta(f);

    for (int i = 1; i < f.size(); ++i) {
        cout << i << ": " << f[i] << endl;
    }
}

実行結果は次のようになります。

高速ゼータ変換の実行結果

7. 活用例 (4): 高速メビウス変換と約数系包除原理

次はいよいよ、競技プログラミングでも実用性の高い高速メビウス変換です。

7-1. 高速メビウス変換

高速メビウス変換とは、高速ゼータ変換の逆変換です。高速ゼータ変換のときと同様に、$n > N$ に対しては $f(n) = 0$ であるような関数 $f(n)$ があって、

$F(n) = \sum_{n | i} f(i)$

で定義される関数 $F(n)$ を考えます。今度は逆に $F(1), F(2), \dots, F(N)$ が与えられたときに、$f(1), f(2), \dots, f(N)$ を求めることを考えます。これが高速に求める処理は高速メビウス変換と呼ばれます ⁹。実はメビウス関数 $\mu(n)$ を用いて、次のように求められることが知られています。

$f(n) = \displaystyle \sum_{n|i} \mu(\frac{i}{n}) F(i)$

この表式をメビウスの反転公式と呼びます¹⁰。具体的に $N = 12$ の場合に書き出してみると、次のようになります。

$f(1) = F(1) - F(2) - F(3) - F(5) - F(7) - F(11) + F(6) + F(10)$
$f(2) = F(2) - F(4) - F(6) - F(10) + F(12)$
$f(3) = F(3) - F(6) - F(9)$
$f(4) = F(4) - F(8) - F(12)$
$f(5) = F(5) - F(10)$
$f(6) = F(6) - F(12)$
$f(7) = F(7)$
$f(8) = F(8)$
$f(9) = F(9)$
$f(10) = F(10)$
$f(11) = F(11)$
$f(12) = F(12)$

これらの表式の意味については、やはり高速ゼータ変換のときと同様に、多次元累積和を考えると分かりやすいです。次のような意味合いを持っています。

ゼータ変換は、累積和をとる
メビウス変換は、差分をとる (累積和の逆操作)

簡単のため、$n$ が $2, 3$ 以外の素因数をもつときには $f(n) = 0$ であるとしましょう。このとき、$n$ が $2, 3$ 以外の素因数をもつときは $F(n) = 0$ も成立しますので、たとえば $f(1)$ は

$f(1) = F(1) - F(2) - F(3) + F(6)$

と表せます。これは下図のような演算を意味しています。累積和を足し引きすること (包除原理) により、確かに $f(1)$ が求められていることがわかります。

さて、高速メビウス変換の実装自体は、メビウス関数を陽に求めなくてもできます。高速ゼータ変換の実装を逆方向に回すイメージで、下のコードのように実装できます。計算量は $O(N \log\log N)$ です。

高速メビウス変換

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

// エラトステネスの篩
vector<bool> Eratosthenes(int N) {
    vector<bool> isprime(N+1, true);
    isprime[1] = false;
    for (int p = 2; p <= N; ++p) {
        if (!isprime[p]) continue;
        for (int q = p * 2; q <= N; q += p) {
            isprime[q] = false;
        }
    }
    return isprime;
}

// 高速ゼータ変換
// 入力 f が in-place に更新されて、F になる
template<class T> void fast_zeta(vector<T> &f) {
    int N = f.size();
    
    // エラトステネスの篩を用いて素数を列挙
    vector<bool> isprime = Eratosthenes(N);
    
    // 各素数 p 軸に対して
    // 大きい座標 (k * p) から小さい座標 (k) へと足し込む
    for (int p = 2; p < N; ++p) {
        if (!isprime[p]) continue;
        
        // 座標が大きい方を起点として累積和をとる
        for (int k = (N - 1) / p; k >= 1; --k) {
            f[k] += f[k * p];
        }
    }
}

// 高速メビウス変換
// 入力 F が in-place に更新されて、f になる
template<class T> void fast_mobius(vector<T> &F) {
    int N = F.size();
        
    // エラトステネスの篩を用いて素数を列挙
    vector<bool> isprime = Eratosthenes(N);
    
    // 各素数 p 軸に対して
    // 小さい座標 (k) から大きい座標 (k * p) を引いていく
    for (int p = 2; p < N; ++p) {
        if (!isprime[p]) continue;
        
        // 座標が小さい方を起点として差分をとる
        for (int k = 1; k * p < N; ++k) {
            F[k] -= F[k * p];
        }
    }
}

int main() {
    // f[0] は関係ないので適当な値を入れておく
    vector<int> f = {-1000, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};

    // 高速ゼータ変換
    fast_zeta(f);

    // 高速メビウス変換で元に戻す
    fast_mobius(f);

    for (int i = 1; i < f.size(); ++i) {
        cout << i << ": " << f[i] << endl;
    }
}

実行結果は次のようになります。たしかに $f$ が元に戻ることが見てとれます。

7-2. 約数系包除原理

ここでは、関数 $F(n)$ から、$f(1), f(2), \dots, f(N)$ の値を求めました。しかし競技プログラミングでは、これら $N$ 種類の値をすべて求める必要はないことも多いです。$f(1)$ のみを求めると、次のようになります。

$f(1) = \sum_{i=1}^N \mu(i) F(i)$

このような表式で $f(1)$ を求めることを、競技プログラミング界隈では、とくに約数系包除原理と呼ぶことがあります。累積和を足し引きするという「包除原理」に基づく表式となっています。

$f(1)$ を直接求めることは難しいが
各 $n$ に対して、$F(n)$ は容易に求められる

という状況は、意外と多くあります。このような場面で約数系包除原理を用いることで、$f(1)$ が求められます。なお高速メビウス変換の実装ではメビウス関数値を陽に求める必要はありませんでしたが、約数系包除原理の実装では、あらかじめ $1$ 以上 $N$ 以下の各整数 $n$ に対するメビウス関数値 $\mu(n)$ を求めておくのが便利でしょう。

例題 (1): オイラー関数

ここで、約数系包除原理の考え方を体験するための具体例として、オイラー関数 $\varphi(N)$ の表式を求めてみましょう。オイラー関数を求める実装はこのページでジャッジできます。オイラー関数とは次のようなものです。

$\varphi(N)$ = $1$ 以上 $N$ 以下の整数のうち、$N$ と互いに素であるものの個数

たとえば小さな $N$ に対しては、

$\varphi(1) = 1$
$\varphi(2) = 1$　($1$ の $1$ 個)
$\varphi(3) = 2$　($1, 2$ の $2$ 個)
$\varphi(4) = 2$　($1, 3$ の $2$ 個)
$\varphi(5) = 4$　($1, 2, 3, 4$ の $4$ 個)
$\varphi(6) = 2$　($1, 5$ の $2$ 個)
$\varphi(7) = 6$　($1, 2, 3, 4, 5, 6$ の $6$ 個)
$\varphi(8) = 4$　($1, 3, 5, 7$ の $4$ 個)
$\varphi(9) = 6$　($1, 2, 4, 5, 7, 8$ の $6$ 個)
$\varphi(10) = 4$　($1, 3, 7, 9$ の $4$ 個)
$\varphi(11) = 10$　($1, 2, 3, 4, 5, 6, 7, 8, 9, 10$ の $10$ 個)
$\varphi(12) = 4$　($1, 5, 7, 11$ の $4$ 個)

という感じです。さて、$\varphi(N)$ を具体的に求めていきたいのですが、「$N$ と互いに素である」という条件は扱いづらそうなので余事象を考えてみましょう。たとえば $N$ が $2$ の倍数であったならば、$2$ の倍数は条件を満たさないので除外したいです。$N$ が $3$ の倍数でもあったならば、さらに $3$ の倍数も除外したいです。しかしこのとき、$6$ の倍数が重複して除外されてしまうことに注意しなければなりません。

ここで感触を掴むため、$N$ が $2, 3$ を素因数に持つがそれ以外の素因数を持たない場合を考えてみましょう。この場合は、$1$ 以上 $N$ 以下の整数のうち、$2$ の倍数でも $3$ の倍数でもないものの個数を求めればよいことになります。下図のようにして、

$\varphi(N) = N - \frac{N}{2} - \frac{N}{3} + \frac{N}{6}$

と求められます。

一般に $N = p_1^{e_1} p_2^{e_2} \dots p_K^{e_K}$ と素因数分解できる場合には、$1$ 以上 $N$ 以下の整数のうち、$p_1, p_2, \dots, p_K$ のいずれの倍数でもないものの個数を求めればよいことになります。このようなややこしい状況は「約数系包除原理」を用いることで解決できることが多々あります (直接、通常の包除原理を適用してもよいです)。具体的には、

$f(n)$ = $1$ 以上 $N$ 以下の整数のうち、$N$ との最大公約数が $n$ であるようなものの個数

として、$f(1)$ を求めたいです。なお、ここで $n$ は $N$ の約数に限ってよいです。なぜなら、$n$ が $N$ の約数でない限り、$n$ と $N$ との最大公約数が $n$ になることはありえないからです。言い換えれば、$n$ が $N$ の約数でないときは $f(n) = 0$ となります。

さて、$f(n)$ のように「$N$ との最大公約数がちょうど $n$ になる」ものの個数を求めることは困難ですが、「$N$ との最大公約数が $n$ の倍数になる」ものの個数 (これは $F(n)$ です) を求めることは簡単です。「$N$ との最大公約数が $n$ の倍数である」とは、「$n$ の倍数である」ということに他なりません。たとえば $N = 12$, $n = 2$ のとき、$2$ の倍数 $2, 4, 6, 8, 10, 12$ はすべて $N = 12$ との最大公約数が $2$ の倍数となります。これらのうちの $6$ などは $N = 12$ との最大公約数は $6$ になりますが、それでも $2$ の倍数であることには変わりありません。以上より、

$F(n) = \sum_{n | i} f(i)$

とおいたときに、

$F(n) = \frac{N}{n}$　($n$ が $N$ の約数であるとき)

となることがわかりました。よって、約数系包除原理を適用することで、オイラー関数は次のように求められます。

$\varphi(N) = f(1) = \sum_{n | N} \mu(n) F(n) = \sum_{n | N} \mu(n) \frac{N}{n}$

さらに一歩

このままでも $O(\sqrt{N})$ の計算量で求められますが、より明快な表式を得てみましょう。$N = p_1^{e_1} p_2^{e_2} \dots p_K^{e_K}$ と素因数分解できるとします。簡単のため $K = 2$ の場合の式を因数分解すると、次のようになります。

$\varphi(N) = N - \frac{N}{p_1} - \frac{N}{p_2} + \frac{N}{p_1 p_2}$

$= N(1 - \frac{1}{p_1})(1 - \frac{1}{p_2})$

より一般には、

$\varphi(N) = N(1 - \frac{1}{p_1})(1 - \frac{1}{p_2}) \dots (1 - \frac{1}{p_K})$

と表せます。とてもシンプルな表式が得られました。

オイラー関数の直感的な意味

オイラー関数の表式は、直感的にも納得しやすいものです。

$1$ 以上 $N$ 以下の整数のうち、$p_1$ の倍数であるものの割合は $\frac{1}{p_1}$ですので、まず $p_1$ の倍数であるものを除外すると、$N(1 - \frac{1}{p_1})$ 個となります
次に残った整数のうち、$p_2$ の倍数であるものの割合は $\frac{1}{p_2}$ ですので、$p_2$ の倍数であるものを除外すると、$N(1 - \frac{1}{p_1})(1 - \frac{1}{p_2})$ 個となります
...
最後に残った整数のうち、$p_K$ の倍数であるものの割合は $\frac{1}{p_K}$ ですので、$p_K$ の倍数であるものを除外すると、$N(1 - \frac{1}{p_1})(1 - \frac{1}{p_2})\dots (1 - \frac{1}{p_K})$ 個となります

8. 活用例 (5)：添字 GCD 畳み込み

最後に、高速ゼータ変換と高速メビウス変換を鮮やかに活用して、添字 GCD 畳み込みを高速に計算するアルゴリズムを求めてみましょう。

$1$ 以上 $N$ 以下の整数 $n$ に対して定義された関数 $f(n), g(n)$ が与えられます。$n > N$ のときは $f(n) = g(n) = 0$ であるとします。このとき

$h(n) = \displaystyle \sum_{{\rm GCD}(i, j) = n} f(i)g(j)$

として、関数 $h(n)$ を求めたいとします。ここで ${\rm GCD}(i, j)$ とは、$i, j$ の最大公約数を表します。添字 $i, j$ を $i + j = n$ に変更した場合は通常の畳み込み演算を表しており、高速フーリエ変換 (FFT) によって高速に求められます。今回は添字の取り方が異なるので一見難しそうですが、実は FFT と同様の発想で解くことができます。$f(n), g(n), h(n)$ をそれぞれ高速ゼータ変換して得られる関数を

$F(n) = \sum_{n | i} f(i)$
$G(n) = \sum_{n | i} g(i)$
$H(n) = \sum_{n|i}h(i)$

としましょう。このとき実は、任意の整数 $n$ に対して次のことが成立するのです！！！

$H(n) = F(n)G(n)$

$N = 12$, $n = 4$ の場合を具体的に見てましょう。まず、

$F(4) = f(4) + f(8) + f(12)$
$G(4) = g(4) + g(8) + g(12)$

となっています。ここで $H(4)$ を計算してみると、

$H(4) = h(4) + h(8) + h(12)$
$= (f(4)g(4) + f(4)f(8) + f(4)f(12) + f(8)g(4) + f(12)g(4)) + (f(12)g(12) + f(8)g(12) + f(12)g(8)) + f(12)g(12)$
$= (f(4) + f(8) + f(12))(g(4) + g(8) + g(12))$
$=F(4)G(4)$

となります。確かに $H(4) = F(4)G(4)$ が成立しています。より一般には、

$H(n)$：$i, j$ の最大公約数が $n$ の倍数であるような $(i, j)$ についての $f(i)g(j)$ の総和
$F(n)G(n)$：$i, j$ がともに $n$ の倍数であるような $(i, j)$ についての $f(i)g(j)$ の総和

と解釈できます。「$i, j$ の最大公約数が $n$ の倍数であること」と「$i, j$ がともに $n$ の倍数であること」は同値ですので、$H(n) = F(n)G(n)$ が成立します。さて、この定理が成立することから、次の方法によって関数 $h(n)$ が求められることがわかります。

$f(n), g(n)$ をそれぞれ高速ゼータ変換して、$F(n), G(n)$ とする
$1$ 以上 $N$ 以下の各整数 $n$ に対して、$H(n) = F(n)G(n)$ を計算する
$H(n)$ を高速メビウス変換して、$h(n)$ を得る

全体として、$O(N \log\log N)$ の計算量で求められます。具体的には次のコードのように実装できます。

添字GCD畳み込み

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

// エラトステネスの篩
vector<bool> Eratosthenes(int N) {
    vector<bool> isprime(N+1, true);
    isprime[1] = false;
    for (int p = 2; p <= N; ++p) {
        if (!isprime[p]) continue;
        for (int q = p * 2; q <= N; q += p) {
            isprime[q] = false;
        }
    }
    return isprime;
}

// 高速ゼータ変換
template<class T> void fast_zeta(vector<T> &f) {
    int N = f.size();
    vector<bool> isprime = Eratosthenes(N);
    for (int p = 2; p < N; ++p) {
        if (!isprime[p]) continue;
        for (int k = (N - 1) / p; k >= 1; --k) {
            f[k] += f[k * p];
        }
    }
}

// 高速メビウス変換
template<class T> void fast_mobius(vector<T> &F) {
    int N = F.size();
    vector<bool> isprime = Eratosthenes(N);
    for (int p = 2; p < N; ++p) {
        if (!isprime[p]) continue;
        for (int k = 1; k * p < N; ++k) {
            F[k] -= F[k * p];
        }
    }
}

// 添字 GCD 畳み込み
template<class T> vector<T> gcd_conv(const vector<T> &f, const vector<T> &g) {
    int N = max(f.size(), g.size());
    vector<T> F(N, 0), G(N, 0), H(N);
    for (int i = 0; i < f.size(); ++i) F[i] = f[i];
    for (int i = 0; i < g.size(); ++i) G[i] = g[i];

    // 高速ゼータ変換
    fast_zeta(F);
    fast_zeta(G);

    // H を求める
    for (int i = 1; i < N; ++i) H[i] = F[i] * G[i];

    // 高速メビウス変換
    fast_mobius(H);

    return H;
}

int main() {
    // f[0], g[0] は関係ないので適当な値を入れておく
    vector<int> f = {-1000, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
    vector<int> g = {-1000, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23};

    // 畳み込み
    auto h = gcd_conv(f, g);

    for (int i = 1; i < h.size(); ++i) {
        cout << i << ": " << h[i] << endl;
    }
}

9. 問題例

ここまででの知見を生かして、競技プログラミングの問題を具体的に解いてみましょう。

例題 (2): AtCoder ABC 206 E - Divide Both

問題概要

整数 $L, R$ ($L \le R$) が与えられます。

以下の条件をすべて満たす整数組 $(x, y)$ の個数を求めてください。

$L \le x, y \le R$
$g = {\rm GCD}(x, y)$ としたとき $g \neq 1$ かつ $\frac{x}{g} \neq 1$ かつ $\frac{y}{g} \neq 1$ が成立する

制約

$1 \le L \le R \le 10^6$
入力はすべて整数

条件の言い換え

まずは問題文の条件を適切に言い換えてみましょう。$\frac{x}{g} \neq 1$ という条件は分かりづらいので、代わりに $\frac{x}{g} = 1$ という条件を考察してみます。

$\frac{x}{g} = 1$

⇔ $x = {\rm GCD}(x, y)$

⇔ $x | y$

ということから、$\frac{x}{g} \neq 1$ とは「$x$ が $y$ の約数ではないこと」と同値です。同様に $\frac{y}{g} \neq 1$ は「$y$ が $x$ の約数でないこと」と同値ですので、問題の条件は次のように言い換えられます。

$g \neq 1$
$x$ は $y$ の約数でない
$y$ は $x$ の約数でない

だいぶ分かりやすくなりました。さらにわかりやすくするために、$x < y$ としましょう。ここで、$x = y$ のときは $g = 1$ より、$x = y = 1$ となるしかなく、$\frac{x}{g} \neq 1$ を違反することに注意しましょう。よって $x < y$ として数え上げて、最後に 2 倍することにします。$x < y$ とすると、条件は

$g \neq 1$
$y$ は $x$ の約数でない

とさらに簡単になります。

数え上げの方針

次の方針で数えることにしましょう。

$A$ ： $x < y$ であって、$g \neq 1$ であるものの個数
$B$ ：$1 < x < y$ であって、$x | y$ であるものの個数

このとき求める答えは $2(A - B)$ となります。$B$ については簡単ですので、先に求めておきましょう。まず $x$ を固定して考えます。$x = L, L+1, \dots, R$ に対して、

$x | y$
$x < y \le R$　($L \le x$ なので、$L \le y$ という条件は無視してよいことに注意)

を満たす整数 $y$ の個数を求めればよいです。これは R / x - 1 と求められます (1 を引いているのは $y = x$ の場合を除外しています)。

約数系包除原理へ

それでは最後に、

$L \le x < y \le R$
$g \neq 1$

という条件を満たす $(x, y)$ を数え上げましょう。このように最大公約数が絡むときは「最大公約数を固定して考えてみる」というのが定石です。つまり、各 $g$ に対して、$L \le x < y \le R$ であって $x, y$ の最大公約数が $g$ となるものの個数を数えようと考えるのです。しかし、「$x, y$ の最大公約数がちょうど $g$」という条件を扱うのは得てして非常に大変です。

そこで定石は「$x, y$ の最大公約数が $g$ の倍数」という条件を満たすものを数え上げて、約数系包除原理に持ち込むことです。$x, y$ の最大公約数が $g$ の倍数であることは、「$x, y$ がともに $g$ の倍数であること」と同値です。こうなるととても扱いやすい条件となります！！！

$f(n)$： $L \le x < y \le R$ であって、$x, y$ の最大公約数が $n$ であるような $(x, y)$ の個数
$F(n)$： $L \le x < y \le R$ であって、$x, y$ がともに $n$ の倍数であるような $(x, y)$ の個数

としましょう。このとき、$F(1) - f(1)$ を求めたいわけですが、約数系包除原理より

$F(1) - f(1) = F(1) - \displaystyle \sum_{n=1}^R \mu(n) F(n) = - \sum_{n = 2}^R \mu(n) F(n)$

という求められます。あとは $F(n)$ を評価できれば OK です。$L$ 以上 $R$ 以下の $n$ の倍数の個数を $c$ とすると、$F(n)$ は「$c$ 個のものから $2$ 個選ぶ場合の数」に一致します。よって

$F(n) = \frac{c (c - 1)}{2}$