· 開発

Project Euler 7: 10001番目の素数を、大きさを2倍ずつ広げるエラトステネスの篩で求める

n番目の素数を篩で求めるとき、どこまで篩えばよいかが前もって分からない問題を、篩を2倍ずつ広げる方法で解きます。試し割りからエラトステネスの篩を導き、2倍に広げても仕事が定数倍で済む理由と、C++17・Python 3 の実装まで扱います。

Project Eulerアルゴリズム数学素数

ネタバレ注意 · この記事には解き方と答えがそのまま書かれています。先に自分で解いてみたい方は、Project Euler の問題 7からどうぞ。

Project Euler シリーズの7問目です。素数を順に数える問題なので、3問目の試し割りをそのまま使っても答えは出ます。この記事はもう一歩進んでエラトステネスの篩を導き、篩を使うたびに出てくる「どこまで篩えばよいのか」という問いに、根拠のない定数ではなく篩を2倍ずつ広げる方法で答えます。篩は10問目など後の問題でも再び使う道具です。

ネタバレ注意: この記事には解法の全体と答えが含まれます。先に自分で解きたい場合は、元の問題から始めてください。このシリーズの記事には広告を載せません。問題文のライセンス(CC BY-NC-SA 4.0)に非営利の条件があるためです。

問題

素数を小さい順に並べると、最初の6個は 2, 3, 5, 7, 11, 13 で、6番目の素数は13です。同じように数えたときの10001番目の素数を求めるのが問題です。

入力はなく、n は10001に固定されています。求めるのは整数1つです。この記事では n 番目の素数を pₙ と書きます(p₁ = 2、p₆ = 13)。

前提知識

素数と合成数。 1より大きい整数のうち、約数が1と自分自身だけのものが素数、それ以外の1より大きい整数が合成数です。1はどちらでもありません。そのため数え始めは p₁ = 2 です。

√m まで割ってみれば十分な理由。 合成数 m = a × b(1 < a ≤ b)なら a × a ≤ a × b = m なので a ≤ √m です。つまり合成数には必ず √m 以下の約数があります。逆に、2から √m までのどの数でも割り切れない m ≥ 2 は素数です。コードでは d × d ≤ m の代わりに同じ意味の d ≤ m / d(整数除算)を使い、掛け算があふれる心配をなくします。正の整数では d ≤ ⌊m / d⌋ と d × d ≤ m はまったく同じ条件です。

最小の約数は素数。 m ≥ 2 の1より大きい約数のうち最小のものを q とすると、q は素数です。q が合成数なら q より小さい1より大きい約数があり、それは m の約数でもあるので、q が最小であることに矛盾します。下の2つの証明で使います。

素数は尽きない(ユークリッド)。 素数が p₁, …, pₖ だけだと仮定します。N = p₁ × … × pₖ + 1 の1より大きい最小の約数は素数ですが、N をどの pᵢ で割っても余りは1なので、その素数はリストにありません。矛盾するので素数は無限にあります。この事実が「篩を広げ続ければいつか n 個集まる」こと、つまり下のアルゴリズムが終わることの根拠です。

エラトステネスの篩。 0から L までの数を並べ、2から順に「まだ消されていない数」を素数として確定し、その倍数を消していく方法です。終わったとき、消されずに残った2以上の数がちょうど L 以下の素数です。「解法の導出」で一歩ずつ組み立てます。

素朴な方法

2から数を1つずつ上げながら試し割りで素数かどうかを調べ、n 番目の素数が出たら止めます。

def is_prime(m):
    if m < 2:
        return False
    d = 2
    while d <= m // d:
        if m % d == 0:
            return False
        d += 1
    return True

count = 0
m = 1
while count < 10001:
    m += 1
    if is_prime(m):
        count += 1
# これで m が10001番目の素数

この方法は正しく、前もって上限を知る必要もありません。費用は数 m 1つにつき割り算が多くても √m 回なので、pₙ までのすべての数を調べると全体で多くても pₙ × √pₙ 程度、つまり O(pₙ√pₙ) です。合成数は小さな約数で早く止まりますが、n 個の素数は毎回 √m まで割り切る必要があり、費用の大部分を占めます。n = 10001 ならこれでも十分速いものの、n が大きくなると √pₙ 倍がすぐに効いてきます。1つ1つの数に「素数か」と尋ねる代わりに、素数から出発して倍数をまとめて消せば、この √pₙ 倍をほとんどなくせます。それが篩です。

解法の導出

ステップ1: 問題文の例を手で確かめる。 2から素数を数えると 2(1番目)、3(2番目)、5(3番目)、7(4番目)、11(5番目)、13(6番目)です。間にある 4, 6, 8, 9, 10, 12 はそれぞれ2か3で割り切れます。p₆ = 13 は問題文の例と一致します。

ステップ2: 尋ねる代わりに消す。 試し割りは m ごとに「どれか小さい数が m を割り切るか」を尋ねます。逆に、素数 p を1つ見つけるたびに p の倍数をすべて「合成数」と印を付けておけば、ある数に着いたとき印がないことだけで素数だと分かります。L = 16 でやってみると次のとおりです。

i = 2: 印なし → 素数。4, 6, 8, 10, 12, 14, 16 を消す
i = 3: 印なし → 素数。9, 12, 15 を消す (6 はもう消えている)
i = 4: 消えている → 飛ばす
i = 5: 印なし → 素数。5 × 5 = 25 > 16 なので消すものはない
i = 6 … 16: 7, 11, 13 だけ印なし → 素数
結果: 2, 3, 5, 7, 11, 13

6個出てきて、6番目は13です。ステップ1と一致します。

ステップ3: 消すのは i × i から、i × i > L なら省く。 素数 i の倍数 k × i のうち k < i のものは、k の素因数(i より小さい素数)の倍数なので、その素数を処理したときにもう消えています。だから i × i から消しても漏れはありません。上の i = 3 で6をもう一度消さなかったのがその例です。同じ理由で i × i > L なら消すものはなく、そういう i は印を確かめて素数のリストに入れるだけです。

ステップ4: 篩には上限が要るのに、pₙ の大きさは分からない。 篩は「L 以下の素数すべて」を与える道具で、問題は「n 番目の素数」を尋ねています。よくあるのは適当に大きな L を決めておくことですが、その L で足りる根拠がなければ、答えが正しい保証もありません。この記事では L を前もって決めません。

  • L = 16 で篩います。
  • 得られた素数が n 個以上なら、そのうち n 番目を返します。
  • 足りなければ L を2倍にして、最初から篩い直します。

n 個以上あれば n 番目が答えである理由。 篩は L 以下の素数を漏れなく、小さい順に返します。ですからリストの n 番目は「n 番目に小さい素数」、つまり pₙ です。リストが L 以下しか含まないことは問題になりません。pₙ より小さい素数はすべて pₙ ≤ L よりも小さいので、全部リストに入っています。

終わる理由。 素数は無限にあるので、L を2倍にし続ければいつか L ≥ pₙ になり、そのときリストには n 個以上が入っています。

小さな例でたどってみます。n = 6 なら L = 16 の最初の篩で6個集まるので、すぐに13を返します。n = 7 なら L = 16 では6個しかなく足りないので、L = 32 で篩い直すと 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 の11個が出て、7番目の17を返します。

ステップ5: 2倍に広げても無駄は大きくない。 最初の L = 16 で終わらなかったとします。最後の篩の大きさを L とすると、直前の L / 2 の篩では素数が n 個未満だったので pₙ > L / 2、最後の篩で終わったので pₙ ≤ L です。

L / 2 < pₙ ≤ L,   つまり   L < 2pₙ

最後の篩は必要な大きさの2倍未満です。それより前の篩を全部足しても 16 + 32 + … + L < 2L < 4pₙ で、最初から正確な上限が分かっていて1回だけ篩った場合と比べて定数倍しか増えません。これが、証明済みの上限なしでも使える理由です。

最終的なアルゴリズム

  1. primes_up_to(L): 大きさ L + 1 の印の配列を作り、i = 2 から L まで、印がなければ素数のリストに加え、i ≤ L / i なら i × i, i × i + i, … を消します。
  2. nth_prime(n): L = 16 から始め、リストが n 個以上になるまで L を2倍にしながら1を繰り返し、リストの n 番目(0から数えて n − 1 番の位置)を返します。
  3. 答えを出力する前に確認します。p₁ = 2、p₆ = 13 であること、そして試し割りで見つけた最初の100個の素数と nth_prime(1) から nth_prime(100) までがすべて一致することです。1つでも違えば0以外の終了コードで終わります。

正しさの証明

  • 素数は消されない。 消される数は i × i + k × i = i × (i + k)(i ≥ 2、k ≥ 0)の形で、1より大きい2数の積、つまり合成数です。
  • L 以下の合成数はすべて消される。 合成数 m ≤ L の最小の約数 q は素数で、前提知識のとおり q ≤ √m ≤ √L です。ループが i = q に来たとき、q は(素数なので)消されておらず、q ≤ L / q なので消す処理をします。m = q × (m / q) で m / q ≥ q なので、m は q × q から始まる消す範囲に入っています。そして q < m なので、これは m に着く前に済んでいます。
  • したがって i に着いたとき印がないことは i が素数であることと同値で、i を2から L まで順に見るので、リストは L 以下の素数すべてを小さい順に含みます。
  • ステップ4で示したとおり、リストが n 個以上なら n 番目は pₙ で、素数は無限にあるので繰り返しは終わります。

計算量

大きさ L の篩1回の仕事を3つに分けて考えます。

  • 篩の構築(消す処理)。 素数 p ≤ √L ごとに多くても L / p 個を消します。すべての i ≤ √L について L / i を足しても L × (1 + 1/2 + …) = O(L log L) なので、これが手軽な上限です。素数の逆数の和が log log 程度にしか増えないという既知の結果(メルテンスの定理)を使えば、より正確には O(L log log L) です。
  • 走査(リストの作成)。 外側のループは i = 2 から L まで印を1回ずつ確かめるので O(L) です。
  • 選択。 リストの n − 1 番を読むのは O(1) です。

2倍に広げる繰り返し全体では、ステップ5のとおり篩の大きさの合計は 2L 未満で、L < 2pₙ(最初の篩で終われば L = 16)なので、全体の時間は O(pₙ log log pₙ)、空間は L + 1 個の印の配列と素数のリストを合わせて O(pₙ) です。素数定理によれば pₙ はおおよそ n ln n なので、n で書けば O(n log n) より少し大きい程度です。素朴な方法の O(pₙ√pₙ) と比べると、√pₙ 倍が log log pₙ 倍に縮みました。

この記事は、素数定理や pₙ の明示的な上限をアルゴリズムの根拠には使っていません。上のおおよその大きさは費用の見積もりのためだけで、正しさは2倍に広げることだけで保証されます。

整数と添字の範囲。

  • 保存される値はすべて最後の L 以下です。印の配列の添字は 0 から L までで、リストに入る素数もすべて L 以下です。L を超えうるのは一時的なループの値だけです。最後の j += i の後の j は j + i < 2L まで進み、配列に渡す大きさ L + 1 も L より 1 大きくなります。2倍に広げる limit *= 2 は篩が足りないときにしか実行されないので、その結果の最大値は最後の L そのものです。n = 10001 では最後の L は 131072 なので(「答えと検証」を参照)、一時的な値も含めてすべての値が 2 × 131072 = 262144 未満で、符号なし32ビットの範囲に十分収まります。一般的な条件は3つ目の項目にあります。
  • C++ では値も添字も std::size_t にしています。印の配列が実際に L + 1 個の要素を占めるので、その配列を作れる環境なら L は std::size_t で表せます。
  • 掛け算 i × i は i ≤ L / i を確かめてからしか計算しないので i × i ≤ L で、あふれません。消すループの j += i は j ≤ L、i ≤ L から出発するので多くても 2L で、次の篩の大きさ L × 2 と同じく、L が std::size_t の最大値の半分以下である限りあふれません。それほど大きな配列はメモリが先に足りなくなるので、現実の制約はメモリです。
  • nth_prime は n ≥ 1 を前提にしています。n − 1 番を読むからです(下の「境界のケース」)。
  • std::vector<bool> は規格で空間を節約する特殊化と定められていて、一般的な実装は1要素を1ビットで持ちます。Python の bytearray は1要素に1バイトを使います。Python の int は精度に制限がないので、あふれることはありません。

C++17 による実装

#include <cstddef>
#include <iostream>
#include <vector>

// Sieve of Eratosthenes over 0..limit. Returns every prime <= limit in increasing order.
std::vector<std::size_t> primes_up_to(std::size_t limit) {
    std::vector<bool> composite(limit + 1, false);
    std::vector<std::size_t> primes;
    for (std::size_t i = 2; i <= limit; ++i) {
        if (composite[i]) {
            continue;
        }
        primes.push_back(i);
        // i <= limit / i is i * i <= limit without computing i * i.
        if (i <= limit / i) {
            for (std::size_t j = i * i; j <= limit; j += i) {
                composite[j] = true;
            }
        }
    }
    return primes;
}

// The n-th prime, for n >= 1. Doubles the sieve limit until the sieve holds at
// least n primes, so no bound on the n-th prime has to be known in advance.
std::size_t nth_prime(std::size_t n) {
    std::size_t limit = 16;
    while (true) {
        const std::vector<std::size_t> primes = primes_up_to(limit);
        if (primes.size() >= n) {
            return primes[n - 1];
        }
        limit *= 2;
    }
}

// Trial division. Used only to cross-check the sieve.
bool is_prime_by_trial(std::size_t m) {
    if (m < 2) {
        return false;
    }
    for (std::size_t d = 2; d <= m / d; ++d) {
        if (m % d == 0) {
            return false;
        }
    }
    return true;
}

int main() {
    if (nth_prime(1) != 2 || nth_prime(6) != 13) {
        return 1;
    }
    std::size_t count = 0;
    for (std::size_t m = 2; count < 100; ++m) {
        if (is_prime_by_trial(m)) {
            ++count;
            if (nth_prime(count) != m) {
                return 1;
            }
        }
    }
    std::cout << nth_prime(10001) << '\n';
    return 0;
}

is_prime_by_trial は素朴な方法の試し割りをそのまま移した関数で、最初の100個の素数で篩と同じ答えを出すかを照合するためだけに使います。篩の消す範囲や2倍に広げる繰り返しに誤りがあれば、この照合で引っかかります。小さな例の確認は assert ではなく if で行いました。assert は NDEBUG が定義されたビルドでは何もしませんが、if の比較と return 1 はプログラムの観測可能な動作なので、-O2 の最適化でも消えません。

Python 3 による実装

def primes_up_to(limit: int) -> list[int]:
    """Sieve of Eratosthenes over 0..limit. Returns every prime <= limit in increasing order."""
    composite = bytearray(limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if composite[i]:
            continue
        primes.append(i)
        # The range is empty when i * i > limit.
        for j in range(i * i, limit + 1, i):
            composite[j] = 1
    return primes


def nth_prime(n: int) -> int:
    """The n-th prime, for n >= 1.

    Doubles the sieve limit until the sieve holds at least n primes, so no bound
    on the n-th prime has to be known in advance.
    """
    limit = 16
    while True:
        primes = primes_up_to(limit)
        if len(primes) >= n:
            return primes[n - 1]
        limit *= 2


def is_prime_by_trial(m: int) -> bool:
    """Trial division. Used only to cross-check the sieve."""
    if m < 2:
        return False
    d = 2
    while d <= m // d:
        if m % d == 0:
            return False
        d += 1
    return True


def main() -> None:
    if nth_prime(1) != 2 or nth_prime(6) != 13:
        raise SystemExit(1)
    count = 0
    m = 2
    while count < 100:
        if is_prime_by_trial(m):
            count += 1
            if nth_prime(count) != m:
                raise SystemExit(1)
        m += 1
    print(nth_prime(10001))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ確認です。C++ の if (i <= limit / i) に当たる条件は別に書いていません。i × i > L なら range(i * i, limit + 1, i) が空になり、同じ働きをするからです。bytearray(limit + 1) はすべての要素が0(印なし)で始まります。小さな例が失敗すれば raise SystemExit(1) で終了コード1を返します。assert は python -O で実行すると取り除かれるので使っていません。

答えと検証

10001番目の素数は 104743 です。

手でたどれる部分は次のとおりです。

  • 篩の大きさ。 篩は 16, 32, 64, … と大きくなり、2¹⁶ = 65536 < p₁₀₀₀₁ ≤ 2¹⁷ = 131072 なので、L = 65536 までは素数が足りず、L = 131072 で止まります。ステップ5の L / 2 < pₙ ≤ L と一致します。
  • 素数かどうか。 上の値は奇数で、各桁の和が19なので3の倍数ではなく、末尾が0でも5でもないので5の倍数でもありません。素数であることの完全な確認には平方根以下(323以下)のすべての数での試し割りが必要で、それは手ではなくプログラムに任せます。

検証は三重です。

  1. 小さな例と最初の100個の照合。 2つのプログラムはどちらも答えを出力する前に p₁ = 2、p₆ = 13 を確認し、試し割りで見つけた最初の100個の素数と、篩が出した1番目から100番目までがすべて一致するかを確認します。1つでも違えば0以外の終了コードで終わります。p₆ = 13 と p₇ = 17 は「解法の導出」で手でもたどりました。
  2. 2つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。2つの出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則上、この検査を通ったコードだけを記事に載せます。
  3. 証明。 篩が L 以下の素数を正確に出すこと、2倍に広げる繰り返しが pₙ を返して終わることは「解法の導出」と「正しさの証明」で示しました。

境界のケース

  • n = 1: 最初の篩(L = 16)の最初の要素2を返します。1は素数ではないのでリストに入りません(外側のループが i = 2 から始まる)。小さな例に p₁ = 2 が入っています。
  • n = 6: 最初の篩にちょうど6個あるので、2倍にせずに終わります。比較が >= だからで、誤って > と書くと無駄にもう1回篩います(答えは同じですが)。
  • n = 7: 初めて篩を広げる必要があるケースで、L = 32 で17を返します。
  • n = 0: 0番目の素数というものはなく、関数は n ≥ 1 を前提にしています。0で呼ぶと、C++ では n − 1 が符号なし整数の最大値に回り込んで範囲外アクセスになり、Python では primes[-1] が最後の要素を黙って返します。このプログラムは n ≥ 1 でしか呼びませんが、別の場所で再利用するなら先に入力を確かめる必要があります。
  • 最初の大きさ: 16は問題文の例が最初の篩に収まるように選んだだけで、1以上ならどの値でも正しく動きます。0から始めると、2倍にしても0のままなので終わりません。
  • i × i から消すこと: 2 × i から消しても答えは同じですが、もう消した数を消し直す仕事が増えます。逆に i × i より大きい数から始めると合成数が残り、間違った答えになります。

発展

  • 1回だけ篩う。 pₙ の証明済みの明示的な上限が分かれば、篩は1回で済みます。その場合、その上限が実際に成り立つ n の範囲も合わせて確かめる必要があります。2倍に広げる方法は、そうした確認なしに、多くても定数倍の仕事を足すだけで安全に終わります。
  • 奇数だけを篩う。 2を別に扱い、配列に奇数だけを入れると、メモリと仕事がおよそ半分になります。添字と数の対応(k 番の位置が 2k + 1)を自分で合わせてみてください。
  • 篩い直さない。 今は L を2倍にするたびに最初から篩い直しています。すでに見つけた素数を使って、新しく増えた区間 (L, 2L] だけを篩うように変えられるでしょうか。区間ごとにどこから消し始めるかが鍵です(区間篩、segmented sieve)。
  • 関連する問題。 このシリーズの3問目では試し割りで素因数を見つけました。10問目は決まった上限より下の素数をすべて扱うので、上限が分かっている場合の篩がそのまま使えます。

出典とライセンス

  • Project Euler, Problem 7 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、原文と同じ CC BY-NC-SA 4.0 で共有します。問題の著作者表示は Project Euler です。
  • Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません(ads: false)。
  • std::vector<bool> は C++ 規格の作業草案 [vector.bool]、assert と NDEBUG は [assertions.assert] に従いました。Python の assert と -O は The assert statement、bytearray は Bytearray Objects に従いました。
  • メルテンスの定理と素数定理は、費用の大きさを見積もるためだけに引いたよく知られた整数論の結果で、この記事では証明していません。アルゴリズムの正しさはこの2つに頼っていません。
  • 解法、証明、コードは自分で書いたものです。Project Euler の解説(overview)、フォーラム、会員限定の資料や他の解答サイトは参照していません。Motitory は Project Euler と提携・後援の関係になく、Project Euler がこの記事を確認したこともありません。

ブログへ