Pythonで素数 2

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ミリ秒。