試し割りによる素数列挙はエラトステネスの篩か？

背景

前回の記事
IchigoJamでエラトステネスの篩 - Qiita
を受け、ツイート主さんから以下のページのURLが送られてきました。

エラトステネスの篩 [WildTree Wiki]

見ると、それぞれの奇数$T$について、$\sqrt{T}$以下の素数で割っていき、
全て割り切れなければ$T$は素数であると判定するアルゴリズムを使っているようです。

一方、エラトステネスの篩といえば、
確定した最小の素数の倍数を候補から外し、確定していない中で最小の候補を素数として確定させる、
という処理を繰り返すことで、効率よく素数のリストを求めることができるアルゴリズムのはずです。

今回送られたページにあるような、1個ずつ試し割りを行うアルゴリズムは、
果たしてエラトステネスの篩といえるのでしょうか？

そこで、試し割りと一般的なエラトステネスの篩の時間コストを比較してみることにしました。

今回の検証対象

今回は、以下のアルゴリズムを対象に検証を行います。

エラトステネスの篩 (工夫なし)

特に工夫をしないエラトステネスの篩です。
配列全体を走査し、素数$p$が見つかった時は$2p$以上の倍数を候補から外します。

def eratosthenes(num):
    isPrime = [False] * 2 + [True] * (num - 2)
    primes = []
    for i in range(num):
        if isPrime[i]:
            primes.append(i)
            for j in range(i + i, num, i):
                isPrime[j] = False
    return primes

エラトステネスの篩 (工夫あり)

2以外の偶数は素数ではないという性質を利用し、奇数のみを対象に走査します。
さらに、整数$m$について、$1<m<p$のとき、$mp$は既に$m$の倍数として候補から外されているはずなので、
見つかった素数$p$に対して$p^2$以上の倍数のみを候補から外します。

def eratosthenes2(num):
    isPrime = [False, False, True, True] + [False, True] * ((num - 4) // 2)
    primes = [2]
    for i in range(3, num, 2):
        if isPrime[i]:
            primes.append(i)
            for j in range(i * i, num, i):
                isPrime[j] = False
    return primes

エラトステネスの篩もどき

「倍数を候補から外す」処理を、素数以外についても行います。
暗黒通信団の『素数表150000個』に載っていたアルゴリズムです。

def pseudoEratosthenes(num):
    flags = [True] * num
    primes = []
    for i in range(2, num):
        if flags[i]:
            primes.append(i)
        for j in range(i, num, i):
            flags[j] = False
    return primes

試し割り (素数限定なし)

3以上の奇数について、3以上その数の平方根以下の奇数で割っていくことで素数判定を行います。

def trialDivision(num):
    primes = [2]
    for i in range(3, num, 2):
        isPrime = True
        for j in range(3, i, 2):
            if j * j > num:
                break
            if i % j == 0:
                isPrime = False
                break
        if isPrime:
            primes.append(i)
    return primes

試し割り (素数限定あり)

3以上の奇数について、3以上その数の平方根以下の素数で割っていくことで素数判定を行います。
今回送られたページに載っていたコードのアルゴリズムです。

def trialDivision2(num):
    primes = [2]
    for i in range(3, num, 2):
        isPrime = True
        for j in primes[1:]:
            if j * j > i:
                break
            if i % j == 0:
                isPrime = False
                break
        if isPrime:
            primes.append(i)
    return primes

時間コストの比較

今回は、計算コストとして「ループを回る回数」を用います。
また、エラトステネスの篩系のアルゴリズムについては、配列を初期化するコストも加えます。

costs.py

def getCosts(num):
    costs = []

    # eratosthenes
    count = num
    isPrime = [False] * 2 + [True] * (num - 2)
    primes = []
    for i in range(num):
        count += 1
        if isPrime[i]:
            primes.append(i)
            for j in range(i + i, num, i):
                count += 1
                isPrime[j] = False
    costs.append(count)

    # eratosthenes2
    count = num
    isPrime = [False, False, True, True] + [False, True] * ((num - 4) // 2)
    primes = [2]
    for i in range(3, num, 2):
        count += 1
        if isPrime[i]:
            primes.append(i)
            for j in range(i * i, num, i):
                count += 1
                isPrime[j] = False
    costs.append(count)

    # pseudoEratosthenes
    count = num
    flags = [True] * num
    primes = []
    for i in range(2, num):
        count += 1
        if flags[i]:
            primes.append(i)
        for j in range(i, num, i):
            count += 1
            flags[j] = False
    costs.append(count)

    # trialDivision
    count = 0
    primes = [2]
    for i in range(3, num, 2):
        count += 1
        isPrime = True
        for j in range(3, i, 2):
            count += 1
            if j * j > num:
                break
            if i % j == 0:
                isPrime = False
                break
        if isPrime:
            primes.append(i)
    costs.append(count)

    # trialDivision2
    count = 0
    primes = [2]
    for i in range(3, num, 2):
        count += 1
        isPrime = True
        for j in primes[1:]:
            count += 1
            if j * j > i:
                break
            if i % j == 0:
                isPrime = False
                break
        if isPrime:
            primes.append(i)
    costs.append(count)

    return costs

delta = 10000
for num in range(delta, delta * 100 + 1, delta):
    print("\t".join(map(str, [num] + getCosts(num))))

実行結果をグラフにすると、以下のようになりました。

試し割り (素数限定なし) の計算コストがぶっちぎりで大きくなっています。
次に、試し割り (素数限定あり) とエラトステネスの篩もどきの計算コストが同程度になっています。
そして、これらに比べエラトステネスの篩はかなり効率がいいようです。

試し割りよりエラトステネスの篩の方が計算コストが少なそうなことがわかりましたが、
これはオーダーレベルで改善しているのでしょうか？それともただの定数倍でしょうか？
これを確かめるため、各アルゴリズムで$5\times {10}^5$未満の素数を求めた時のループ回数を$1$とし、
ループ回数の変化の仕方を比較してみます。

試し割り (素数限定なし) のループ回数の増え方が一番速く、
試し割り (素数限定あり) もエラトステネスの篩よりは試し割り (素数限定なし) に近い増え方をしています。
エラトステネスの篩もどきは、試し割りよりかなりエラトステネスの篩に近い増え方をしています。
そして、エラトステネスの篩2種類は、かなり近い増え方をしています。

他の記事

素因数分解とかエラトステネスの篩(ふるい)とかのメモ - 唯物是真 @Scaled_Wurm

試し割り (素数限定なし) の時間計算量は$O(n\sqrt{n})$
試し割り (素数限定あり) の時間計算量は$O(n\sqrt{n} / (\log n)^2)$
エラトステネスの篩の時間計算量は$O(n \log \log n)$

のようです。
エラトステネスの篩もどきの時間計算量は載っていないようでした。

エラトステネスの篩の計算量を実験で確認する - Qiita

こちらの記事には実際の計算時間のグラフが載っています。
また、エラトステネスの篩 (工夫なし) とエラトステネスの篩 (工夫あり) の計算量のオーダーは変わらないらしいです。

結論

エラトステネスの篩 [WildTree Wiki]
に載っている、奇数を1個ずつ素数で割ることで素数かどうかを判定するアルゴリズムは、
エラトステネスの篩よりオーダーレベルで効率が悪く、エラトステネスの篩ではないといえるでしょう。

なお、このページでは、エラトステネスの篩は先頭にリンクが載っているだけで、
「エラトステネスの篩を使った」とは書かれていないので、間違いとはいえないでしょう。
しかし、このページには以下のツイートと同じ秒数が書かれており、
以下のツイートのための実験に使われたコードが載っていると推測できます。
そして、ツイートにははっきりと「エラトステネスの篩で」「処理をさせた」と書かれているので、
このツイートは間違っているといえるでしょう。

エラトステネスの篩で、32768未満の素数をリストアップする処理をさせたところ、LPC1114版では968秒ほどかかったのが、RISC-V版では202秒と、約4.8倍高速になってた。> IchigoDake

— 荒木 宏之 (@h_araki) January 6, 2022