質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

Q&A

解決済

5回答

2501閲覧

1億以下の素数判定が遅い(区間篩 javaから書き換え)

opyon

総合スコア1009

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

0グッド

0クリップ

投稿2018/10/08 16:38

編集2018/10/09 06:49

#結論:JavaやCと比べるとPythonで外部ライブラリを使わずに高速化することは厳しい。
PythonはJavaやC ++に比べて非常に遅いのに、なぜ機械学習のような高速アルゴリズム処理にPythonが使われるのですか?

#追記
回答やコメントでご指摘くださった方ありがとうございました。

ここまでの回答やご指摘を取り込み修正した現状です。
修正後実行結果10^8
コードは長いのでGitHubへ

実用的かどうかは別ですが、まだ改善の余地があるのではないかと試行錯誤しておりますので引き続きヒントや気づいたことがあればご教示頂ければと思います。

素数判定に使われていると思われるアルゴリズム?が何なのか知りたい
以前頂いた回答からエラトステネスの篩を使った素数判定の実装は出来ました。

その後偶然見つけた記事から、エラトステネスの区間篩なるアルゴリズムを知り実装してみました。
元がjavaで書かれていたものをなるべくそのままpythonに書き換えたつもりなのですが旧来ののみに比べて倍近く遅くなりました。
そのままでは動かない部分を自分なりに改修したのが影響してるのかなと思います。

1億以下の正の整数を素数判定した実行速度比較
javaコード :1秒未満
旧来の篩  :11秒以上
今回の区間篩:24秒以上

####知りたいこと
書き換えでどこか間違っているのか
改善出来そうなことがあれば何かしらヒントやご教示頂けたら助かります。

import timeit #旧:エラトステネスの篩 def eratosthenes_list(n): primes = [True] * n result = [2] for prime in range(3, n, 2): if primes[prime]: result.append(prime) for i in range(prime * prime, n, prime * 2): primes[i] = False return result # 旧 nn = 100000000 print(timeit.timeit(lambda: (eratosthenes_list(nn)), number = 1)) # nn = 100000000 # 11.725847288999997
import timeit import math def primesSieve(n): sieve = [True] * n sieve[0] = sieve[1] = False for i in range(2, n): for j in range(i * 2, n, i): sieve[j] = False return sieve #新:区間篩 def primesList(n): SIEVE_MAX = int(math.sqrt(n)) primes = [] sieve = primesSieve(SIEVE_MAX) for i in range(len(sieve)): if sieve[i]: primes.append(i) for i in range(1, int(n / SIEVE_MAX)): sieve = [True] * len(sieve) start_index = SIEVE_MAX * i for p in primes: if (p * p > start_index + SIEVE_MAX): break jj = start_index + (p - (start_index % p)) % p # print(jj, start_index + SIEVE_MAX, p) for j in range(jj, start_index + SIEVE_MAX, p): sieve[j - start_index] = False for j in range(0, len(sieve)): if sieve[j]: primes.append(start_index + j) return primes # 新 nn = 100000000 print(timeit.timeit(lambda: (primesList(nn)), number = 1)) # print(timeit.timeit(lambda: print(primesList(nnn)), number = 1)) # nn = 100000000 # 24.349327380000002

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答5

0

区間篩にした方はふるい落としを素数の2倍から始めているせいでは?

python

1# 旧篩 2for i in range(prime * prime, n, prime * 2) 3 4# 区間篩 5for j in range(i * 2, n, i)

おまけに区間篩の方では素数でない数をスキップしていませんから、余計に時間を食うと思います。

投稿2018/10/09 01:41

swordone

総合スコア20649

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

hayataka2049

2018/10/09 01:53

回答のコード書きながらプロファイル取ってましたが、そこはたかだかsqrt(n)までの篩なので、大勢に影響はなかったです。確かに改善点ですが
opyon

2018/10/09 02:05

ご指摘ありがとうございます。 >区間篩にした方はふるい落としを素数の2倍から始めているせいでは? 仰る通りです。 コードを書き換える時に上手くいかなかったので元のjavaの形に合わせたままでした。 まだ改善の余地がありそうですね。 現状@hayataka2049さんの回答を参考に修正がほぼ終わりテストしているところなので落ち着いたらそこを重点的に見てみます。
opyon

2018/10/09 02:42

>区間篩にした方はふるい落としを素数の2倍から始めているせいでは? 10^8までの比較テストが終わったので、早速その部分のみ差し替え合わせて回答に追記リンク貼りました。 >おまけに区間篩の方では素数でない数をスキップしていませんから、余計に時間を食うと思います。 こちらは具体的にこれから探ってみますが引き続きヒントなどあればご教示お願いします。
opyon

2018/10/09 06:07

>区間篩の方では素数でない数をスキップしていませんから これは具体的にどこのことを指しておられますか? 区間毎に素数リストを使ったループで素数分スキップしてるはずですが・・・ 他に見落としてるようでしたらご指摘ください。 for j in range(jj, SIEVE_MAX, p): sieve[j] = False
guest

0

ベストアンサー

省メモリ化を意図しているアルゴリズムなので、そんなに速くならないのでは? という気がします。


とりあえず、挙動が怪しかったので、アルゴリズムを元のjavaのコードにある程度忠実にして、エラトステネスの篩と同じ結果になるようにしておきました。

python

1import timeit 2import math 3 4#旧:エラトステネスの篩 5def eratosthenes_list(n): 6 primes = [True] * n 7 result = [2] 8 for prime in range(3, n, 2): 9 if primes[prime]: 10 result.append(prime) 11 for i in range(prime * prime, n, prime * 2): 12 primes[i] = False 13 return result 14 15#新:区間篩 16def primesSieve(n): 17 sieve = [True] * (n+1) 18 sieve[0] = sieve[1] = False 19 for i in range(2, n+1): 20 for j in range(i * 2, n+1, i): 21 sieve[j] = False 22 return sieve 23 24def primesList(n): 25 n -= 1 # なぜか必要だった 26 SIEVE_MAX = int(math.sqrt(n)) + 1 27 primes = [] 28 sieve = primesSieve(SIEVE_MAX-1) 29 30 for i, x in enumerate(sieve): 31 if x: 32 primes.append(i) 33 34 for i in range(1, (n // (SIEVE_MAX))+1): 35 sieve = [True] * len(sieve) 36 start_index = SIEVE_MAX * i 37 38 for p in primes: 39 if (p * p > start_index + SIEVE_MAX): 40 break 41 42 j = start_index + (p - (start_index % p)) % p 43 while j - start_index < SIEVE_MAX: 44 sieve[j - start_index] = False 45 j += p 46 47 48 j = 0 49 while j < SIEVE_MAX and start_index + j <= n: 50 if sieve[j]: 51 primes.append(start_index + j) 52 j += 1 53 54 return primes 55 56for nn in range(10, 10**4): 57 a = eratosthenes_list(nn) 58 b = primesList(nn) 59 if not a == b: 60 print(nn) 61 print("miss") 62 print(len(a), len(b)) 63 print(a) 64 print(b) 65 break 66else: 67 print("all ok") 68 69nn = 10**6 70print(timeit.timeit(lambda: (eratosthenes_list(nn)), number = 1)) 71print(timeit.timeit(lambda: (primesList(nn)), number = 1)) 72

実はこの状態だと質問文の状況より単純なエラトステネスの篩と差が開くのですが……

一番内側のwhileは明らかに冗長なので、とりあえず単純化します。

python

1 j = (p - (start_index % p)) % p 2 while j < SIEVE_MAX: 3 sieve[j] = False 4 j += p

改めてforにします。whileより余計な処理が減らせれば速くなるはずです。

python

1 jj = (p - (start_index % p)) % p 2 for j in range(jj, SIEVE_MAX, p): 3 sieve[j] = False

同様に最後のwhileもforにしておきます。

python

1 for j in range(min([n - start_index + 1, SIEVE_MAX])): 2 if sieve[j]: 3 primes.append(start_index + j)

ここまでのコード。

python

1def primesList(n): 2 n -= 1 3 SIEVE_MAX = int(math.sqrt(n)) + 1 4 primes = [] 5 sieve = primesSieve(SIEVE_MAX-1) 6 7 for i, x in enumerate(sieve): 8 if x: 9 primes.append(i) 10 11 for i in range(1, (n // (SIEVE_MAX))+1): 12 sieve = [True] * len(sieve) 13 start_index = SIEVE_MAX * i 14 15 for p in primes: 16 if (p * p > start_index + SIEVE_MAX): 17 break 18 19 jj = (p - (start_index % p)) % p 20 for j in range(jj, SIEVE_MAX, p): 21 sieve[j] = False 22 23 for j in range(min([n - start_index + 1, SIEVE_MAX])): 24 if sieve[j]: 25 primes.append(start_index + j) 26 27 return primes

この状態で10^8までだと、手元環境ではエラトステネスの篩の倍弱の実行時間になります。

軽くプロファイリングした感じだと、一番内側のforループで大半の時間を食っているようです。

追記

最終的に、numbaのちからで区間篩の方が速くなるところまで持っていきました。numbaのjitコンパイルを通すためにコードをまただいぶ書き換えてあります。

python

1import timeit 2import math 3 4import numba 5 6#旧:エラトステネスの篩 7@numba.jit(nopython=True) 8def eratosthenes_list(n): 9 primes = [True] * n 10 result = [2] 11 for prime in range(3, n, 2): 12 if primes[prime]: 13 result.append(prime) 14 for i in range(prime * prime, n, prime * 2): 15 primes[i] = False 16 return result 17 18#新:区間篩 19@numba.jit(nopython=True) 20def primesList(n): 21 n -= 1 22 SIEVE_MAX = int(math.sqrt(n)) + 1 23 24 sieve = [True] * (SIEVE_MAX) 25 sieve[0] = sieve[1] = False 26 for i in range(2, SIEVE_MAX): 27 for j in range(i * 2, SIEVE_MAX, i): 28 sieve[j] = False 29 30 primes = [] 31 for i, x in enumerate(sieve): 32 if x: 33 primes.append(i) 34 35 for i in range(1, (n // (SIEVE_MAX))+1): 36 sieve = [True] * len(sieve) 37 start_index = SIEVE_MAX * i 38 39 for p in primes: 40 if (p * p > start_index + SIEVE_MAX): 41 break 42 43 jj = (p - (start_index % p)) % p 44 for j in range(jj, SIEVE_MAX, p): 45 sieve[j] = False 46 47 tmp = n - start_index + 1 48 if tmp < SIEVE_MAX: 49 x = tmp 50 else: 51 x = SIEVE_MAX 52 for j in range(x): 53 if sieve[j]: 54 primes.append(start_index + j) 55 return primes 56 57if True: 58 for nn in range(10, 10**3): 59 a = eratosthenes_list(nn) 60 b = primesList(nn) 61 if not a == b: 62 print(nn) 63 print("miss") 64 print(len(a), len(b)) 65 print(a) 66 print(b) 67 break 68 else: 69 print("all ok") 70 71nn = 10**8 72print(timeit.timeit(lambda: (eratosthenes_list(nn)), number = 1)) 73print(timeit.timeit(lambda: (primesList(nn)), number = 1)) 74""" => 75all ok 761.2398602159810252 771.059988280001562 78"""

投稿2018/10/08 21:42

編集2018/10/08 23:47
hayataka2049

総合スコア30933

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

退会済みユーザー

退会済みユーザー

2018/10/08 22:14

スクリプトの実行速度だと原因を探すのが難しいですよね。計算量的には区間篩のほうが確かに速くなるようです。 hayataka2049さんに質問なのですが、このレベルのコードでメモリ不足のため実行時間が遅くなるといったことはあり得るのでしょうか。 もしメモリ不足が原因がとしたら、3重のfor文を手書きイテレータに書き換えてしまうと良いのかなと思いました。
hayataka2049

2018/10/08 22:50 編集

>ALMIさん あまりこのアルゴリズムに詳しくないのと、探してもそんなに情報が出てこないので(とりあえず日本語でググっただけだけど)あまり何も言えないのですが、そんなに時間計算量変わりますか? 空間計算量がケチれるという話は見当たりますが (逆に、空間計算量をケチってるのに時間計算量の面でメリットがあるのかという気も)
hayataka2049

2018/10/08 22:58 編集

メモリ不足はたぶんないです。10**8のboolaen配列が、今実測しましたが800MBくらい。一要素あたり8byteで作れます
hayataka2049

2018/10/08 23:33 編集

>opyonさん ループの条件とかを元のコードに忠実にしただけで、そんなにいじっていないです。あと、エラトステネスの篩の方とこっちのコードで「以下」と「未満」の違いがあるみたいだったので、nをデクリメントしております
hayataka2049

2018/10/08 23:32

どうでも良いですが、エラトステネスの篩の方にnumba.jitつけたら10倍速くなった・・・
opyon

2018/10/08 23:35

10倍以上なら劇的高速化といえますね^^
opyon

2018/10/08 23:39

>元のコードに忠実にしただけで、 そうなんですよね。 最初は自分流に書いてたのですが上手く動かせなくて最後になるべく忠実に書き直したつもりだったのですが、javaのfor文の条件式を上手く再現出来てなかったようです。 @hayataka2049さんの書き直しを見てその辺りを工夫されていて流石だなと感銘を受けております。
hayataka2049

2018/10/08 23:48

jitコンパイルでちゃんと速くなりました。コンパイルを通すためにだいぶ苦労しましたが・・・
opyon

2018/10/08 23:51

凄いですね。 jitコンパイルなるものも初耳なので修正が終わったら調べてみます。
opyon

2018/10/09 00:25

>この状態で10^8までだと、手元環境ではエラトステネスの篩の倍弱の実行時間になります。 (jitで10倍速出ているので)今更なのですが、初期回答時点では篩よりも区間篩の方が遅かったということで良いですか? 上下逆に見て勘違いしていました。 回答に実行結果画像貼っています。
hayataka2049

2018/10/09 01:28

そうですね、倍の実行時間はそのままの意味です
opyon

2018/10/09 01:33

ありがとうございます。 アルゴリズムが同じでも言語によるということなのでしょうか・・・ どちらにしても前回同様pythonで他言語と比較するような速度を求めるのは良くないという教訓を得た気がします。
hayataka2049

2018/10/09 02:12

まあ、どこのオーバーヘッドが大きいかとかは言語によるので、計算量が同じでも優劣が覆ることはありえますね 割とレアケースだとは思います
opyon

2018/10/09 06:18

レアケースですか。jitなどを使わないとするとこの辺りが限界なのかもしれませんね。
guest

0

そもそもPythonはコンパイル言語ではないので同じアルゴリズムでも実行結果が遅くなります。
高速化しないのであればCythonなどをつかってみてはいかがですか?
参考

投稿2018/10/08 21:08

Nippun

総合スコア1147

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

opyon

2018/10/08 23:25

ありがとうございます。 修正が落ち着いたら参考にしてみます。
guest

0

旧来の方は、偶数を考慮したものですが、今回の方は、考慮されていないと思われます。

追記
クリップ後に追加された回答をみてなかったです。すでにかれていました、すいません。

投稿2018/10/09 03:26

編集2018/10/09 03:30
tmp

総合スコア277

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

Zuishin

2018/10/09 03:27

間違った出力がされるということですか?
opyon

2018/10/09 03:35

ありがとうございます。 @swordoneさんのご指摘から、初期素数リストの取得で旧来の篩を使ったものに差し替えています。 コメントに書くと長すぎるので自己回答の中に実行結果とコードへのリンク貼っています。 >追記 入れ違いでしたね、大丈夫です。
Zuishin

2018/10/09 03:37

回答ではなく質問を編集して追記してください。
opyon

2018/10/09 03:40

すみません質問に追記して良いのか分からずに回答に書いてました。 質問の追記としても編集しておきます。(書いてしまった回答はどうしましょ・・・)
Zuishin

2018/10/09 03:45

自己解決の回答は見る人が混乱しないように一行目に注釈を入れておいてください。
opyon

2018/10/09 03:53

ご指摘ありがとうございます。質問と回答を編集しました。
Zuishin

2018/10/09 03:56

そうではなく、質問の補足は質問を編集して追記してください。 あちこちに情報が分散されていると見る人が混乱します。 回答の方は、回答ではないこととそこに書いてしまった経緯を書いて質問に誘導してください。
Zuishin

2018/10/09 03:58

質問の補足は基本的にコメントや回答ではなく、質問を編集することでしてください。 コメントはあくまでも追記・編集依頼とそれに対する返答で、質問は質問文だけで完結するようにしてください。 新しく質問を見る回答者をあちこち振り回してはいけません。
opyon

2018/10/09 06:11

分かりづらくてすみません。 書き直したつもりだったのですが、まだ至らない点があるようでしたら編集します。
guest

0

#※コメント欄にコードを返信すると見にくいので回答を使っています※
これは@hayataka2049さんから頂いた回答へやコメントへの返信です。

(追記2は質問を編集したので削除)

追記1
実行結果

@hayataka2049さん回答より差異がある箇所を抜粋

これを参考に自分のコードも修正してみます。
(コメント欄に書くと見にくいので回答を使っています。)

#差異がある箇所 ① def primesSieve(n): sieve = [True] * (n+1) ② def primesList(n): n -= 1 SIEVE_MAX = int(math.sqrt(n)) + 1 ③ for i in range(1, (n // (SIEVE_MAX))+1): ④ jj = (p - (start_index % p)) % p for j in range(jj, SIEVE_MAX, p): sieve[j] = False ⑤ for j in range(min([n - start_index + 1, SIEVE_MAX])): if sieve[j]: primes.append(start_index + j)

私の質問にあるコード

#差異がある箇所 ① def primesSieve(n): sieve = [True] * n ② def primesList(n): SIEVE_MAX = int(math.sqrt(n)) ③ for i in range(1, int(n / SIEVE_MAX)): ④ jj = start_index + (p - (start_index % p)) % p for j in range(jj, start_index + SIEVE_MAX, p): sieve[j - start_index] = False ⑤ for j in range(0, len(sieve)): if sieve[j]:

投稿2018/10/08 23:24

編集2018/10/09 04:22
opyon

総合スコア1009

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問