素数列の生成

Python で Project Euler（http://projecteuler.net/ ）（日本語訳は http://odz.sakura.ne.jp/projecteuler/）をとき始めた。

現在 3 問目。素因数、素数、エラトステネスの篩、と調べていたら、次のたとえ話にたどり着いた。

雲梯を進む怠惰な猿たち

雲梯を進む怠惰な猿たち

3以上の奇数の番号が振られている横棒が半無限に並んでいる雲梯。

そこに不思議な管理人が現れる。管理人は 3 という番号の棒に何もないことを見つけると、2段飛ばしで進んでいく猿を 9 という番号が振られている棒にぶら下げ、3 という番号を伝書鳩に託して休憩する。猿は怠け者なので進めと言われないと進まない。

管理人の元へ、子どもがやってきて何やら言葉を交わし、頭を下げると帰っていく。管理人は腰を上げ、5 に猿がいないことを確認すると、4段飛ばしで進んでいく猿を 25 という番号の棒にぶら下げる。やはり伝書鳩に今度は 5 という番号を託して休憩する。

もう一度同じようなことがあって 7 という番号を付けた鳩が飛んでいった後、次に子どもがやってきた時、管理人のすることが少し変わる。 2段飛ばしで進んでいく猿が 9 番の棒にぶら下がっているのを見ると、管理人は猿に進むように促す。猿は眠たそうに 15 の棒に移る。管理人はそれを見届けると、11 番に何もいないことを確認し、10段飛ばしで進んでいく猿を 121 番の棒にぶら下げる。 11 という番号を付けた鳩が飛んでいく。

眠たくなるような景色は 37 の鳩が飛んでいった後、ちょっとした変化を迎える。いつものように管理人の元へ子どもがやってきて去っていく。管理人は立ち上がると 39 の横棒にぶら下がる猿に進むように促す。猿は2段飛ばしで 45 の棒に行くがそこには既に4段飛ばしの猿がいる。すると、あの怠惰な4段飛ばしの猿が歯を剥き威嚇する。 2段飛ばしの猿は逃げるように、でもいつも通り2段飛ばしで 51 の棒に移り、そこに落ち着く。一つの横棒に2匹がぶら下がることはできずに、後から来た猿が更に移動していくようだ。管理人は何事も無かったかのように、40段飛ばしの猿を遠くの方にぶら下げ、鳩に 41 の番号を託して放すと、いつものように休憩する。

http://lowlife.jp/mft/weblog/2006/04/21.html

なるほど、こうやってもエラトステネスの篩を実装できるのか。

つまり、

• 雲梯に振られた番号、3以上の奇数列、すなわち初項 3 公差 2 の等差数列、イテレータ
• n - 1 段飛ばしの怠惰な猿、初項 n ** 2、公差 n * 2 の等差数列、イテレータ
• 雲梯の状態、つまりぶら下がっている猿たち、今猿がぶら下がっている番号をキー、猿を値とした辞書

となる、と。とりあえず等差数列ジェネレータを作れば脇役はそろう。残る主役、鳩使いの不思議な管理人は一言ではあらわせないけれど、そこまで難しくはない。というわけで、書いてみた。

（たとえ話だと、すすんでいる最中の猿は別の猿がすでにぶら下がっている箇所を飛ばして進むとなっているが、このコードだと管理人がすぐ次に行くように指示をだしているように見える…。まぁ、たとえ話の忠実再現が目的ではないし、これでいいか。）

# coding: utf-8

u"""素数列の生成 - 遅延評価によるエラトステネスの篩
"""

def _seq_with_common_difference(first, diff):
    u"""初項 first 公差 diff の等差数列"""

    while True:
        yield first
        first += diff

def seq_of_prime():
    u"""素数列"""

    yield 2
    nums = _seq_with_common_difference(3, 2)
    sequences = {}
    for num in nums:
        seq = sequences.pop(num, None)
        if seq is None:
            yield num
            seq = _seq_with_common_difference(num ** 2, num * 2)
            i = next(seq)
        else:
            for i in seq:
                if i not in sequences:
                    break
        sequences[i] = seq

def test():
    from itertools import takewhile
    for i, prime in enumerate(
            takewhile(lambda x: x < 100, seq_of_prime())):
        print i, prime

if __name__ == '__main__':
    test()

エピローグ

この雲梯はいつまでも使えるように見えますが、猿が一定以上に増えると重みで壊れてしまうのです。

http://lowlife.jp/mft/weblog/2006/04/21.html

seq_of_prime は無限長のイテレータなので、使用の際には終了条件をつけるのを忘れずに。とはいえ、生成した素数の数だけ sequences 辞書にジェネレータが蓄えられていくので、有限のメモリ上で動かすといつか止まってしまう。

>>> from prime import seq_of_prime
>>> for i in seq_of_prime():
...   pass
...
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "prime.py", line 36, in seq_of_prime
    sequences[i] = seq
MemoryError