「素数」を編集中 (節単位)

===巨大素数の確率的判定===
大きな数が素数であるかを決定的に判定することは難しいが、[[ミラー・ラビン法]]により確率的に判定することができる。
<source lang="perl">
#!/usr/bin/perl

use strict;
use Math::BigInt;
use Math::BigInt::Random qw/ random_bigint /;

my($result);
if(@ARGV==0 || @ARGV>2) {
	die("Usage: miraarabinntesuto.pl <number to check> [loop count]\n");
} elsif(@ARGV==2) {
	$result=&miraarabinntesuto($ARGV[0],$ARGV[1]);
} else {
	$result=&miraarabinntesuto($ARGV[0]);
}
print $result?"Prime!\n":"Not prime...\n";
exit(0);

# http://d.hatena.ne.jp/pashango_p/20090704/1246692091

# if the number is prime, return 1. otherwise, return 0.
sub miraarabinntesuto {
	my($q,$loopnum)=@_;
	my($s,$d,$a,$r,$i,$j,$ok);
	my($powresult,$pownow,$powcurrent);
	if(!$loopnum){$loopnum=50;}
	$q=Math::BigInt->new($q);
	if($q<=1){return 0;}
	if($q==2){return 1;}
	if($q%2==0){return 0;}
	# q-1 = 2^s * d
	$s=0;$d=$q-1;
	while($d%2==0) {
		$s++;$d/=2;
	}
	for($i=0;$i<$loopnum;$i++) {
		# 1～q-1までの間の数を適当に選びaとし
		$a=random_bigint( min => "1", max => $q-1);

		# 「pow(a,d,q) != 1」がTrueだった場合
		$powresult=Math::BigInt->new("1");
		$pownow=$a;
		$powcurrent=$d;
		while($powcurrent>0) {
			if($powcurrent & 1){$powresult=($powresult*$pownow)%$q;}
			$pownow=($pownow*$pownow)%$q;
			$powcurrent>>=1;
		}
		if($powresult==1) {next;}

		# 0～s-1までのカウンタをrとし、すべての「pow(a, 2**(r*d), d)!=-1」がTrue
		# もしかして：pow(a, d*(2**r), q)!=-1
		$ok=1;
		# powresultは前の結果を利用
		for($r=0;$r<=$s-1;$r++) {
			if($powresult==$q-1){$ok=0;last;}
			if($powresult==1){last;}
			$powresult=($powresult*$powresult)%$q;
		}
		if($ok){return 0;}
	}
	return 1;
}
</source>