HSPTV!掲示板

　Wikipediaのミラー-ラビン素数判定法にRubyのソースがあったので
HSPに移植してみました。
といっても、アルゴリズムは理解できないし、Rubyはまるで知らないのですが(^_^;;
　46,690までは正しく判定されるようですが、
46,691以降の素数の多くが合成数と判定されます。

１　何がいけないのでしょう？
　　（私の移植がまずいのか、乱数がまずいのか、Rubyの実装がまずいのか、
　　元々決定的アルゴリズムではないので、アルゴリズムの限界なのか。
　　でも、合成数を素数と判定する可能性はあっても素数を合成数と判定する
　　可能性はなさそうにみえる）

２　正しく判定させるように修正する方法があれば教えてください

#module "primemod"
#defcfunc ModMathPow_ int base_, int power_, int mod
	base=base_
	power=power_
	result=1
	repeat: if power<=0: break
		if power&1: result=result*base\mod
		base=base*base\mod
		power=power>>1
	loop
return result
#defcfunc isprime int n
	if n<0: n=-n
	if n==2: return 1
	if n==1: return 0
	if (n&1)==0: return 0
	d=n-1
	repeat: if d&1: break
		d=d>>1
	loop
	f=1
	repeat 20	;←20.times do
		a=rnd(n-2)+1
		t=d
		y=ModMathPow_(a,t,n)
		repeat: if n-1==t: break
			if y==1: break
			if n-1==y: break
			y=y*y\n
			t=t<<1
		loop
		if n-1!=y: if (t&1)==0: f=0: break
	loop
return f
#global

	ini=43000
	add=(strlen(str(ini))+1)*8
	x=0: y=0
	repeat 4000,ini
		if isprime(cnt) {
			pos x,y
			mes cnt
			y+=16: if y+16>480: y=0: x+=add
		}
	loop
stop

乱数に問題があるとすればrndの乱数発生範囲が32767以内に限られているところですかね

>ミラー-ラビン素数判定法
HSPの整数型までの範囲(〜2147483647)なら、試し割りでも1回10msほどで処理できるので実装する意味が無い。

というのは無粋なので、真面目に考えてみた。
46691をチェックした際に、「素数じゃない」として返す部分を調べると、
どうやら
>if n-1!=y: if (t&1)==0: f=0: break
のところがクサいらしい。ちょっとyを表示させてみると、

あ　る　時　な　ぜ　か　負　数　に　な　っ　て　た

そこで気がついた。HSPのint型の制限は2147483647。つまり、それ以上だと負数になることがままある。だが、
>y=y*y\n
ここで桁あふれが発生しているのではないのか？
もちろん同じことは
>base=base*base\mod
にも言える。別にRubyは使っていないので、(変数の)正しい値を知っているわけではないが

　なるほど、オーバーフローでしたか。
改めて見直してみると、抜けるのは46691からではなく46441からでした(^_^;;
baseが46340(=$7FFFFFFFの平方根)を超えたらダイアログを表示するようにすると、
46441の素数判定時にオーバーフローしていることがわかりました。

　ありがとうございました。