Pythonで素数、Haskellで素数 2の続き。
イテレータとかジェネレータとかあやふやだったんで、素数ジェネレータしてみる。
普通よくみかけるのはこんなやつ。
from itertools import ifilter, count, dropwhile
def sieve():
g = count(2)
while True:
prime = g.next()
yield prime
g = ifilter(lambda x, prime=prime: x % prime, g)
first_prime_GE_sieve = (lambda n , ps = sieve() :
(dropwhile( lambda x : x < n, ps ) ).next()
)
print first_prime_GE_sieve(52951)
# 結果
52951
(結果はどれも paiza.io でtime outにならずに計算できた最大の数。いつも同じではないが速さの目安として。)
なるほど。
素数を出力しつつ、無限長?のイテレータに対してふるいにかけている感じ。
よくできてる。
中味がわかってきたら、前記事でいただいたコメントのコードと同じだといまさら気付く。
関数かオブジェクトかの違いだけ。
あんまり速くはないかも?
Haskell界隈のやりかた。
primes :: Integral a => [a]
primes = map fromIntegral $ [2, 3] ++ primes'
primes' = [5] ++ recursivePrimes 1 7 primes'
recursivePrimes m s (p : ps) =
zonedPrimes m s p ++ recursivePrimes (m * p) (p * p) ps
zonedPrimes m s p =
[n | n <- croppedPrimables s p, gcd m n == 1]
croppedPrimables s p =
[x + y | x <- [s, s + 6 .. p * p - 6], y <- [0, 4]]
main = print . take 1 . dropWhile (<24684009) $ primes
-- 結果
[24684013]
自分でわかりやすいように変形してます。
オリジナルはこちら。
ちょっとずつ計算して積み上げていく感じ。
けっこう速い。
そっちに寄せてみる。
このアルゴリズムを、Pythonでやってみる。
from fractions import gcd
from itertools import chain, tee, dropwhile
cropped_primables = (lambda s, p :
[x + y for x in xrange(s,p * p - 5, 6) for y in [0, 4]]
)
zoned_primes = (lambda m, s, p :
[n for n in cropped_primables(s, p) if gcd(m, n) == 1]
)
def primes() :
primes_list = [2, 3, 5]
last_prime = primes_list[-1]
primes_iter = iter(primes_list)
seeds_iter = iter([])
m = 1
s = 7
p = 5
while True :
prime = primes_iter.next()
yield prime
if prime == last_prime :
primes_list = zoned_primes(m, s, p)
last_prime = primes_list[-1]
it1, it2 = tee(primes_list)
primes_iter = it1
seeds_iter = chain(seeds_iter, it2)
m = m * p
s = p * p
p = seeds_iter.next()
first_prime_GE = (lambda n , ps = primes() :
(dropwhile( lambda x : x < n, ps ) ).next()
)
print first_prime_GE(3003079)
# 結果
3003079
Haskellの方は遅延評価を使っているので、Pythonでそこをどうするか?
でも、アルゴリズムとしては次の素数がわかればいいだけなので、イテレータを、値を返す用と、内部で素数を計算する用とに分けて工夫してみた。
ちなみに簡易的だが、IDLE上で引数100000で一回実行して時間を計ってみたら
>>> measure_t(first_prime_GE,100000)
2.5987625122070312e-05
約0.026 ミリ秒。
冒頭の方は
>>> measure_t( first_prime_GE_sieve,100000)
16.222843885421753
約16秒。
けっこう速いんじゃない?
追記
上のは実行時エラーを起さないようにやってみたけど普通はこうなのかな?
# coding: utf-8
from fractions import gcd
from itertools import chain, tee, dropwhile
cropped_primables = (lambda s, p :
(x + y for x in xrange(s,p * p - 5, 6) for y in [0, 4])
)
zoned_primes = (lambda m, s, p :
(n for n in cropped_primables(s, p) if gcd(m, n) == 1)
)
def primes() :
primes_iter = iter([2, 3, 5])
seeds_iter = iter([])
m = 1
s = 7
p = 5
while True :
try :
prime = primes_iter.next()
except StopIteration :
primes_base_iter = zoned_primes(m, s, p)
it1,it2 = tee(primes_base_iter)
primes_iter = it1
seeds_iter = chain(seeds_iter, it2)
m = m * p
s = p * p
p = seeds_iter.next()
else :
yield prime
first_prime_GE = (lambda n , ps = primes() :
(dropwhile( lambda x : x < n, ps ) ).next()
)
print first_prime_GE(2728100)
##結果
2728129
StopIterationをキャッチして処理するように変更。
配列じゃなくてよくなったものはイテレータに変更。
見た目はすっきり。けど、少し遅くなったよう。
同様に実行時間を測ると
>>> measure_t(first_prime_GE,100000)
3.409385681152344e-05
0.034ミリ秒。