Last Updated: 2/6/2024, 5:44:57 AM
Photo by Google

# エラトステネスの篩

エラトステネスでググったら出てきて、さすがにカオス過ぎて笑ってしまいました。 エラトステネスの篩というのは n 以下の整数の素数を求めるアルゴリズムです。

原理自体はとても簡単で、たくさんの数を用意して、次にひたすら割り続けます。 割ることができた数は振い落とします。

以下では 1 から 169 の数を用意して 1, 2, 3, 5, 7, 11, 13 と素数で割り続けて素数を見つけます。

1 ~ 169 の数字をご用意しました。 これをエラトステネスの篩にかけていきます。

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154
  • 155
  • 156
  • 157
  • 158
  • 159
  • 160
  • 161
  • 162
  • 163
  • 164
  • 165
  • 166
  • 167
  • 168
  • 169

４つの書き方を比較検討しました。

項番 項目 時間
1 簡単なもの 9.0777 [msec]
2 高速化 1.6104 [msec]
3 OrderedDict 7.5887 [msec]
4 dict 5.8195 [msec]

計測に使用したコードはこちら。

git clone https://github.com/niconico25/comparison_of_python_code/
cd comparison_of_python_code
python3 9_2_sieve_of_eratosthenes.py

# 1. 簡単なもの

def primes(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, n + 1):
        for j in range(i * 2, n + 1, i):
            is_prime[j] = False
    return [i for i in range(n + 1) if is_prime[i]]

参考にしたサイトは、下記のサイトから抜粋しました。

# 2. 高速版

5 倍くらい早くなる。

def primes(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if not is_prime[i]:
            continue
        for j in range(i * 2, n + 1, i):
            is_prime[j] = False
    return [i for i in range(n + 1) if is_prime[i]]

for 文の方が while 文よりも速い気配があるので、for 文で書きました。 その検証は下のサイトでやりました。for 文の方が綺麗に書けますしね。 素因数分解するときは while 文の方が早かったです。なぜかは、わかりません。

# 3. OrderedDict

# ◯ 背景

枝狩りをした 2. 高速版のアルゴリズムでも、 既に素数でないと判定されたものにも重複して、素数でないと判定しています。 同じ j に対して複数回 is_prime[j] = False を実行しています。

例えば 3 * 4, 3 * 6, 3 * 8 のような数でそのような重複がなされています （素数 * 既に素数と判定された数の倍数）。

そこで順序付き集合を用いて、 同じ処理を繰り返し行わないようにしたいと考えました。

# ◯ 速度

しかしそこまで速くならず、簡単なものより 1.5 倍速いくらい。 アルゴリズム的には一番速いはずだけど。

Python は、四則演算やらなんやらが遅いので、 何が原因かはいまいちよくわかりません。

# ◯ 本当に重複してないの？

素数の判定を重複して行いません。１度だけしか行いません。 本当にそうなの？みたいな感じですが、重複して削除していた場合 del number_dict[non_prime] で例外を投げ返されます。

# ◯ non_prime_list

その時は for 文の中では del number_dict[non_prime] ができません。 できいないので上のコードでは non_prime_list を用いて 一時的に素数ではない数をリストとして保存している。

# ◯ ポイント

ポイントは以下の箇所です。 ここの動作は上のアニメーションと照らし合わせると、 追いやすいかなと思います。

def primes_ordered_dict(n):

    ...

        # 次の約数 = prime X prime
        #
        # なので、もし prime X prime が n より大きくなったら break
        if not (prime * prime <= n):
            break
        non_prime_list.append(prime * prime)

        # 約数 = prime X k
        # k は prime より大きいまだ篩い落されていない数
        for k in number_dict:
            if prime * k <= n:
                non_prime_list.append(prime * k)
            else:
                break

# ◯ コード

Python 3.6 以前の set は Ordered, 順序付きではないので OrderedDict で代用しました。

import collections


def primes_ordered_dict(n):
    #
    prime_list = []
    non_prime_list = []
    number_dict = collections.OrderedDict(
        ((i, None) for i in range(n + 1)))
    #
    del number_dict[0]
    del number_dict[1]
    while True:
        # get prime
        prime, _ = number_dict.popitem(False)
        prime_list.append(prime)
        # break
        if not (prime * prime <= n):
            break
        # get non primes
        non_prime_list.append(prime * prime)
        for k in number_dict:
            if prime * k <= n:
                non_prime_list.append(prime * k)
            else:
                break
        # delete non primes
        for non_prime in non_prime_list:
            del number_dict[non_prime]
        non_prime_list.clear()

    # add rest of numbers as prime
    prime_list.extend(
        [prime for prime, _ in number_dict.items()])

    return prime_list

# 4. dict

Python 3.6 以前では辞書は順序が保たれていなかったのですが、 Python 3.7 から順序付きになりました。

What’s New In Python 3.7 (opens new window)
Python data model improvements:

書いてみたのですが、そこまで速くできませんでした。 while 文を多用せざるを得ず、結果的に Python のコードが長くなってしまったためです。

Python は遅いので C で実装されているコードに処理を投げられれば大抵速くなります。 逆に言えば Python で書くコードが長くなると遅くなってしまいます。

なぜ while 文で書いていたかというと、 for 文だと正順 (FIFO) で要素を取り出すのに dict.popitem だと逆順(LIFO)で要素を取り出すためです。

>>> reversed({})
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'dict' object is not reversible
>>>

popitem() (opens new window)
任意の (key, value) 対を辞書から消去して返します。 対は LIFO の順序で返却されます。

以下が実装したコードです。

def primes_dict(n):
    #
    prime_list = []
    non_prime_list = []
    k_list = []
    number_dict = {
        i: None for i in reversed(range(n + 1))}
    #
    del number_dict[0]
    del number_dict[1]
    while True:
        # get prime
        prime, _ = number_dict.popitem()
        prime_list.append(prime)
        if prime * prime > n:
            break

        # get non primes
        non_prime_list.append(prime * prime)
        while True:
            k, _ = number_dict.popitem()
            k_list.append(k)
            if prime * k > n:
                break
            non_prime_list.append(prime * k)

        while k_list:
            number_dict[k_list.pop()] = None

        # delete non prime
        while non_prime_list:
            del number_dict[non_prime_list.pop()]

    while number_dict:
        prime, _ = number_dict.popitem()
        prime_list.append(prime)

    return prime_list

ちなみにこの更新は、日本人の方がコミットされています。 初めて知ったときは驚きました。

# 5. ワンライナー

これはエラトステネスの篩では、ありません。速度的には 1 より少し遅いくらいでした。 自分のパソコンだと i=3 を代入すると計算が終わらない。しかし、こんなやり方、思いつくなんて凄すぎる...。 Twitter の "だよ〜" ってコメントなんだよ、軽すぎるよ (´；ω；`)ﾌﾞﾜｯ

# 素数リスト 2**i 桁以下
primes = lambda i: primes(i - 1) + [k for k in range(10**2**(i-1), 10**2**i) if all(k % p for p in primes(i - 1))] if i else [2, 3, 5, 7]

# 素数リスト 2**i 桁以下
def primes(i):
    if i:
        # 素数リスト 2**(i-1) 桁以下
        L = primes(i - 1)

        # 素数リスト 2**(i-1)+1 桁以上 2**i 桁以下
        n = range(10**2**(i-1), 10**2**i)
        M = [k for k in n if all(k % p for p in L)]

        # 素数リスト 2**i 桁以下
        return L + M
    else:
        # 素数リスト 2**0 桁以下
        return [2, 3, 5, 7]

# ◯ ポイント

def primes(i):
        ...
        # 素数リスト 2**(i-1)+1 桁以上 2**i 桁以下
        n = range(10**2**(i-1), 10**2**i)
        M = [k for k in n if all(k % p for p in L)]
        ...

2**(i-1)+1 桁以上 2**i 桁以下の合成数（素数でない数）は、 2**(i-1) 桁以下の素数を素因数に持ちます。逆に言えば、2**(i-1)+1 桁以上 2**i 桁以下の合成数（素数でない数）は、2**(i-1)+1 桁以上の素数を素因数に持ちません。

2 桁の合成数は、1 桁以下の素因数を持ちます。3 桁以上4 桁以下の合成数は、2桁以下の素因数を持ちます。5 桁以上 8 桁以下の合成数は、4 桁以下の素因数を持ちます。

例えば 2 桁の合成数で 1 桁の素因数を持たない数を考えて見てください。存在しません。逆に言えば、例えば 2 桁の素数を使って 2 桁の合成数を作って見てください。作ることはできません。

これは、任意の 2 桁の素数 p, q について考えると、必ず p * q >= 100 となるからです。

Python でエラトステネスの篩を使って素数のリストを作る
← 約数 はじめに →