Pythonで素数日を求める

パクリ元

Pythonで「素数日」を求めたら処理が長すぎて日が暮れそうになった
素数日っていうのを初めて聞いて面白かったので自分でもやってみる。

コード

素数を求める方法はいろいろあるけど、僕が知っている中で一番早そうなエラトステネスの篩を使った。
$N$までの数値を素数判定したいのなら、$\sqrt{N}$までの素数を一回数え上げれば十分である。

prime.py

from math import sqrt, ceil
from itertools import count, takewhile
import datetime
from datetime import date

# maxの素数判定に必要な素数を数え上げる
def sieve(max):
    c = ceil(sqrt(max))
    s = [2]
    yield 2
    for n in range(3, c, 2):
        for m in s:
            if n % m == 0:
                break
            if n < m * m:
                s.append(n)
                yield n
                break
        else:
            s.append(n)
            yield n

# 素数のリストsieveを使ってnが素数か判定する
def is_prime(n, sieve):
    for m in sieve:
        if n % m == 0:
            return False
    else:
        return True

# [begin, end)の範囲の日付の数値表現のリストを求める
def daynum(begin, end):
    days = (begin + datetime.timedelta(days=i) for i in count())
    days = takewhile(lambda d: d < end, days)
    return [int(d.strftime('%Y%m%d')) for d in days]

def prime_days(begin, end):
    ns = daynum(begin, end) 
    si = list(sieve(ns[-1]))
    return (d for d in ns if is_prime(d, si))

def main():
    from datetime import date
    for d in prime_days(date(2019, 1, 1), date(2119, 1, 1)): # 2019/01/01 ～ 2119/01/01
        print(d)

if __name__ == "__main__":
    from timeit import timeit
    import sys
    print(timeit(main, number=1), file=sys.stderr)

結果

100年分計算しても1秒かからなかった。

$ python ./prim.py > out.txt
0.3709903

素数日はこんな感じ。100年で2205日だった。
次の素数日は2019/07/19らしい。

out.txt

20190221
20190227
20190301
20190319
20190323

……

21181003
21181007
21181103
21181109
21181117

追記 (2019/06/27)

なんか改めて見るとワチャワチャしててもっと簡単に書けそうだったので書き直した。
要するにフィボナッチ数列と同じで同じ計算を二回しなければいいんだろ、という発想でlru_cacheを使ってメモ化した。

memo.py

from math import sqrt, ceil
from itertools import count, takewhile
import datetime
from functools import lru_cache

def daynum(begin, end):
    days = (begin + datetime.timedelta(days=i) for i in count())
    days = takewhile(lambda d: d < end, days)
    return [int(d.strftime('%Y%m%d')) for d in days]

def prime_days(begin, end):
    ns = daynum(begin, end) 
    return (d for d in ns if is_prime_memo(d))

@lru_cache(maxsize=None)
def is_prime_memo(n):
    if n == 2:
        True
    if n % 2 == 0:
        return False
    for m in range(3, ceil(sqrt(n)), 2):
        if not is_prime_memo(m):
            continue
        if n % m == 0:
            return False
    else:
        return True

def main():
    from datetime import date
    for d in prime_days(date(2019, 1, 1), date(2119, 1, 1)):
        print(d)

if __name__ == "__main__":
    from timeit import timeit
    import sys
    print(timeit(main, number=1), file=sys.stderr)

思った以上に遅くなってしまった。
割る数を二乗してるループの終了判定が足をひっぱてるんじゃないかと思う。
まあでもこっちのほうがきれいに書けてるし一般的だよね。

$ python ./prim.py > out.txt
2.1147849