Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
4
Help us understand the problem. What is going on with this article?
@ttatsf

Pythonで素数 2

More than 3 years have passed since last update.

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

4
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
4
Help us understand the problem. What is going on with this article?