HSPTV!掲示板

HSPで1〜1,000,000までの素数を取り出すプログラムを作ったのですが、約2.05秒でした。
(桁数ではありません)

表示に時間がかかるので、あくまで内部で計算して変数に入力するまでの時間です。

全て独学でやった結果なのでそれなりに満足しています。

答え合わせの意味も含め、最も早く素数を取り出せるプログラムの書き方が知りたいです。

よろしくお願いします。

以下のサイトのソースコードを少し改変。
http://monooki.ldblog.jp/archives/33709065.html

#define ctype pek(%1,%2)  ((%1(%2/8)>>(%2\8))&1)	;%1の%2ビット目を抽出
#define ctype pok1(%1,%2) %1(%2/8)|=1<<(%2\8)		;%1の%2ビット目に1を書き込む
#uselib "winmm.dll"
#cfunc timeGetTime "timeGetTime"
n=1000000			;nまで調べる
dim map,n/8+1	;nビット目が1ならnは合成数、0ならnは素数 pek(map,n)
map(0)=3		;0、1を素数判定にするのを回避
開始時間=timegettime()
for i,0,n,1
	;今までに割り切られていないか判定
	if pek(map,i)=0{
		;素数の倍数にフラグ(合成数の目印)を立てる
		for j,i+i,n,i : pok1(map,j) : next
	}
next
title ""+(timegettime()-開始時間)+"ms"

Intel i7-2600で1.25秒です。

ありがとうございます。

ビットはわからないのですがやってることに無駄がないので最善なんでしょうね。

ちなみにですがビットを使わないプログラムで一番早いものは分かりませんか？

自分の2秒くらいのプログラムは相当早い部類に入るのではないかと思っているのですが・・・

環境に依存するので「時間」は意味を持ちません。

個々人のパソコンで、あなたのスクリプトと比較した結果じゃないと
意味ナシです。

ちなみに
○２以上の偶数は素数では無い。（２以上の素数は全て奇数）
○任意の数で最大の素数は、その数の平方根以下。
が原則です。

数の範囲を限定すれば、素数の倍数を除外していく消去方は有効です
が、無制限とした場合は使えません。なので上限の数値により最適な
アルゴリズムは変わるので、条件を明確にしましょう。

Corei7-4770シリーズのPCでソラさんのコードを実行したところ平均910msで処理を完了しました。
素晴らしいですね。しかしながらこのコードは速度を犠牲にしてメモリの消費を最小限に留める配慮が
加えられているようです。HSPで最高速化するように修正すると700msまで高速化しました。
いずれも1秒かかりません。ループ処理も含め根本的に最適化すれば600ms以下にまで高速化し、
当初よりおよそ35％も早く処理を完了します。xxxさんはビットがわからないと仰っていますが、
ソラさんのコード上ではビット演算を行っているものの、アルゴリズムにおいてビットは無関係です。
最適化されたコードをOpenHSPのプロジェクト（Dish等）を利用してネイティブコンパイルした場合の
実行速度は目を見張るものがあることでしょう。
（これをするかしないかによっても最高速となる書き方は変わると思います。）

xxxさんは質問の段階で「答え合わせ」と言っておきながら自分で書いたコードを記載していません。
にもかかわらずソラさんは的確なヒントとサンプルコードを提示してくださいました。
それなのにあなたは意味も理解しないまま「その方法は使わないで一番早い方法はないですか？
　自分のコードは2秒で処理できる。相当早い部類になると思う。」などと滅茶苦茶な返信をしています。
これはないでしょう。
KAさんも仰っていらっしゃいますが、前提となるものを隠蔽して質問をしないでください。
これでは折角回答した人が哀れです。

まず掲示板で質問をするのなら（どこでもそうですが）情報の隠蔽はせず、回答者が的確な返信を
できるよう最大限の努力をしてください。今回の件であれば事のあらましと、あなたの書いた
ソースコードを記載し、これより高速化できますか？と質問するのが確実です。
そうするとこの記事を後から訪れた人にとっても有用な情報となり資産となります。
「自分の自称優れたコードは秘密にしたい。後から来た人の役に立てられてはたまらない。」
などの後ろめたい理由がある場合はここで質問するのはやめてください。

返信ありがとうございます。
一応両方実行した結果で書いてます。
両方の数値もコメントしてたのですが、とある事情でそのコメント削除してしまったのでもう一回書きます。

ソラさんの約1125ms　　自分の約2050ms

数の上限を決めないでそれなりにスピードが出るものを作りたいです。


ソラさんのプログラムより劣っているソースコードを書くことに意味があるのかと思うのですが、ご希望とのことですのでソースコードを書きます。

時間取得をgettimeからソラさんが使っていたtimeGetTimeに変えてあります。

申し訳ないのですが、こちらは自分が作った一つ前のバージョンになります。
自分のパソコンでは約2400msかかります。

最新版はこれとはまた少し違うアイデアを取り入れており、現在そのアイデアをさらに進化させる方法を思いついたので、試してから公開します。
これは別に隠したいとかではなく、ここに書くと簡単に最適化されたやり方を示されると思うので、それを見る前に自分でやるだけやっておきたいからです。

ソースコードの見にくさはご容赦を。
ループの内容はもっと最適化できそうですね。

#uselib "winmm.dll"
#cfunc timeGetTime "timeGetTime"

dim s,16777215:dim r,16777215　　;sは素数用　rは素数の除算数用
s.0=1:s.1=2　　　　　;素数1と2を代入
r.0=0:r.1=1　　　　　;素数の除算上限回数を代入

e=1000000　;eまで調べる
n=3　　　　;nから開始
x=1　　　　;素数の除算可能回数　nに対するx
y=1　　　　;除算実行回数

開始時間=timegettime()

*main
if n>e:goto *mesb
r.n=x:i=2          ;0、1を素数判定にするのを回避

z=sqrt (n):z=int (z):if ( z \ 2 ) = 0:z-　　;平方根を調べる
y=r.z:if n<10:y=x-1    ;0で除算するのを回避

if ( n \ 3) = 0 and ( n ! 3 ):n=n+2:r.n=x:i+:y-　;3割切で次の数へ移行し、5の素数の除算から始める

repeat y
if ( n \s.i ) = 0:j=1:break　　　;素数判定
i+:loop

n=n+2:if j=1:j=0:goto *main

x+:s.x=n-2:goto *main

*mesb
	title ""+(timegettime()-開始時間)+"ms"+x

>隠したいとかではなく、ここに書くと簡単に最適化されたやり方を示されると思うので、それを見る前に自分でやりたい

誰が返信すんのこれ。
なら勝手にやってくださいって感じです。

ソースコードも書かずにって言われたから書いたらこの仕打ち。
何がしてほしいか書かなかった自分が悪いですね。

ソースコード見ればわかると思いますが、自分が何を目的としていたかを説明します。


素数抽出プログラムで一番時間が掛かるのが「割り切れるか調べる」部分です。

たとえば83が素数かどうかを調べるためには83未満の素数で割り切れるか調べるわけですが、
「2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79」の22回も調べなければなりません。

実際は2や√83以上の素数は調べなくてもいいので、突き詰めていけば「3,5,7」の素数で割り切れなければ素数判定ができます。
22回が3回で収まれば大幅な短縮になるのは明らかです。
しかしループの中に√83以上の素数になったらループを抜ける計算式を書くと余計な時間が掛かります。

この83という数を過不足なく3回の素数の除算で収めたいと思って書いたのが上記のものです。

「素数で割る回数を過不足なく必要最低限の回数に収める」

83という数値から3を出せる計算式など自分では思いつきませんでしたので余計な処理をしています。
ここが改善できれば間違いなく今より早くなりますのでアドバイスをお願いします。

誤解されないように追記しておきます。

自分が時間短縮に成功した部分とさらに改良を加えたい部分は今説明したところとは関係がありません。

「素数で割る回数を過不足なく必要最低限の回数に収める」部分においては
知識不足でノータッチですので、ぜひアドバイスをお願いします。

エラトステネスのふるい　というアルゴリズムがあってですね
数値ひとつを　変数なりビット１つで表現して
一回素数じゃないと判定された要素は　変数なりビットなりを　素数じゃないと表現すれば
割り算の回数は減ります

詳しくは「エルネストのふるい」でググってくださいな

昔、自分が素数のプログラムを組んだ時は
2n+1確実に2で割り切れない数字なので
2+1 4+1 6+1...と素数っぽいと判定します
割り算じゃなく掛け算と足し算なので　処理が軽くなります
6n+1と6n+5で表現したらもっと要素が減りますが
どこまでそういう要素を用意するかは
書き手次第　みたいな

また2n+1で表現する場合
3で割り切れる数字は3*3から始めたり
5なら5*5から始めるといいでしょう
5*3は3で割り切れる的な　チェック済みみたいな
大きい数字まで適用すると　幸せになれますよー　的な

自分がやった時は　インラインアセンブラを使ってみたら
やたら早くなってびっくりした思い出があります

ただ　自分の場合
どう正しく求めたか　チェックの仕方がわからなくて　やめましたｗ

横から失礼します。Dripです。
xxxさん、こんにちは。

>数の上限を決めない

当初の質問と言っていることが大きく変わってきている気がします。
そもそも上限がないならint型で処理をしている時点で天井が見えていると思います。
64bit版HSPを使えば32bitアプリメモリ上限の壁を突破できそうですが、
配列管理がintなのでやっぱり天井が見えています。
多倍長整数の拡張機能を独自に作るか外から持って来るかして、
sqrt関数等も独自に整備する必要があると思います。
配列管理も改造して擬似的に多倍長整数に対応してください。
これらは動作速度に大きく影響を与える部分になるでしょう。

現在のxxxさんのスクリプトはメモリを食いまくるので
素数が21億に達するよりはるか手前でエラーする気がします。
（xxxさんの方法ですと調査範囲２百万あたりにつき10MB程度消費するようなので
　32bitアプリメモリ上限の縛りを受け調査範囲（〜8億程度）が求められる素数の限界？）
ソラさんのスクリプトのようにメモリを節約する工夫もすべきかと思います。
（ソラさんの方法ですとint型の限界値スレスレの調査範囲（〜21億程度）が求められる素数の限界？）

数の上限を決めずに素数を取得し続けるプログラムを書くとしたら、
求めた素数はHDDにダンプしてメモリを食い続ける構造を何とかするのが先だと思います。
メモリを食い続ける構造が何とかなったら多倍長整数を管理するためのシステムを整備してみてください。
範囲はintの範囲まででいい…ということであれば、もう答えは１レス目で出てしまった気がします。
というかY_.repeatさんのエラトステネスのふるいそのものです＾＾；

兎に角お手軽・高速に！ということであれば、以下が最速じゃないでしょうか。
（ソラさんよりメモリ食いますが、理論上21億まで行けると思います。）
私の環境では一千万までの素数もわずか6秒で取得してしまいました。

#uselib "winmm.dll"
#cfunc time "timeGetTime"
n=1000000:sdim map,n:t=time()	//とりあえず100万まで取得してみる。（n=取得範囲）
repeat n-2,2
	if peek(map,cnt)=0:for p,cnt,n,cnt:poke map,p,1:next:poke map,cnt,2	//2は素数フラグ
loop
//★以下表示部分（表示できる範囲で表示します）
color ,,255:mes ""+n+"個の素数取得にかかった所要時間："+(time()-t)+"ms":color
y=20
repeat n
	if peek(map,cnt)=2:mes cnt:if ginfo_cy>ginfo_winy-20:pos ginfo_cx+50,20:if ginfo_cx>ginfo_winx:break
loop

返答ありがとうございます。

Dripさんのが早すぎて満足しました。
もう素数を抽出するプログラムを組むことはないと思います。
色々書き込んで頂きましたが、そもそもの知識量が足りないので正直半分も理解できていません。



最後に、初めに書いた自分が2秒で計算できるプログラムに使用した考え方は

∔2ではなく、+4,+2,+4,+2,+4,+6,+2,+6　以下ループです。
∔2ループだと数えなくていい数が被りまくるので必要のないものを飛ばします。
これを更に改良してもう少し短縮しようと考えていました。

素数の世界では常識？
現行プログラムでは必要ないから無視されてるのかな。

それではありがとうございました。

再度すみません。

+4,+2,+4,+2,+4,+6,+2,+6ではなく∔2ループが使われている理由は

わざわざ加算数値を変更する手間のほうが時間が掛かるからでしょうか？
それともそもそも違う方法で計算しなくてもいい数値を飛ばしているのでしょうか？

最後にそれだけ知りたいので教えてくれると助かります。

HSPによる配列参照は重いのです。
そのほかにもrepeatとforではrepeatの方がずっと速いという特性を生かし、
repeatにはcntという自動で増える変数があるにもかかわらず、
あえてforも使わずrepeat内でcnt_を手動で+2ずつしてrepeatのループ数を減らすと…

#uselib "winmm.dll"
#cfunc time "timeGetTime"
n=1000000:sdim map,n+1:t=time()	//とりあえず100万まで取得してみる。（n=取得範囲）
poke map,2,2:cnt_=3				//偶数は最初の２だけ
repeat (n-cnt_)/2+1,cnt_	//以降の偶数は全て無視で探す
	if peek(map,cnt_)=0:for p,cnt_,n,cnt_:poke map,p,1:next:poke map,cnt_,2	//2は素数フラグ
	cnt_+=2
loop
//★以下表示部分（表示できる範囲で表示します）
color ,,255:mes ""+n+"個の素数取得にかかった所要時間："+(time()-t)+"ms":color
y=20
repeat n+1
	if peek(map,cnt)=2:mes cnt:if ginfo_cy>ginfo_winy-20:pos ginfo_cx+50,20:if ginfo_cx>ginfo_winx:break
loop

おやおや、もっと速くなったりします＾＾；；；おもしろいですね。
処理速度に関しては自分でベンチを取って色々調べてみるのがよいかと思います。
以下に５百万回の処理を５回テストして平均を比較するベンチマークを貼ります。

#uselib "winmm.dll"
#cfunc time "timeGetTime"
workCnt=5000000	//処理のループ回数
testCnt=5		//同じ試験を何回繰り返して平均を取るか
screen 0,(testCnt+2)*110,480

font msgothic,12:objsize 64,32,24	//テスト用紙の作成
repeat testCnt+1:a=cnt:repeat ginfo_winy/24
	hsvcolor 0,0,250-cnt\2*20-a\2*5:boxf a*110,cnt*24,a*110+110,cnt*24+24
loop:loop

color:dim ANS,100,testCnt	//ANS(処理項目,処理回数)にかかった時間を代入していく。
repeat testCnt
	pos cnt*110+5,5:ptn=0:a=0:b=0:mes "READY...":wait 10
	
	t=time():repeat workCnt:	a+	//この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a+:         "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	a++	//この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a++:        "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	a+=1	//この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a+=1:       "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	a=a+1	//この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a=a+1:      "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	b=a*2	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "b=a*2:      "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	b=a/2	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "b=a/2:      "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	if b=a	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "if b=a:     "+ANS(ptn,cnt):ptn++:wait 10

	t=time():repeat workCnt:	if b=a(0)	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "if b=a(0):  "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	b(0)=16:	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "b(0)=16:    "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	poke b,0,16:	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "poke b,0,16:"+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	a=b(0)	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a=b(0):     "+ANS(ptn,cnt):ptn++:wait 10
	
	t=time():repeat workCnt:	a=peek(b,0)	;この式と↓の項目名を合わせること
	loop:ANS(ptn,cnt)=(time()-t):mes "a=peek(b,0):"+ANS(ptn,cnt):ptn++:wait 10
loop

color ,,255:pos testCnt*110+5,5:mes "平均時間（小さい程速い）"
repeat ptn
	ptn=cnt:a=0
	repeat testCnt
		a+=ANS(ptn,cnt)
	loop
	mes "----    "+(a/testCnt)
loop

b(0)=16　と　poke b,0,16では　poke b,0,16が速いのに
a=b(0)　と　a=peek(b,0)では　a=peek(b,0)が遅い。
今回の処理ではpokeの方がずっと多いしpokeはこの４つの中で一番速いので、
peekは重いけどフラグに使う変数は断然dimよりsdimでpeek/pokeを使った方が速い…
なんて比較しながらプログラムしてたら日が暮れそうですが、
頭の片隅に置いておいてプログラムするとHSPの速いスクリプトが作成できるかもしれません。

追伸、私が最初に投稿したスクリプトはループ回数が１回足りなかったです。
２回目に書いたようにn+1として要素数とループ回数を直せば
ループ最後のひとつが素数だった場合の処理が正しくなります。
お恥ずかしい…

丁寧に解説ありがとうございます。
処理速度がそれぞれ違うのを比較対象して一番早いものでやらないとダメということですね。
色々試してみます。