「素数」を編集中

'''素数'''とは、1とその数でしか割り切れない2以上の[[整数]]であり、[[東工大生]]の大好物である。

==素数判定==

=== [[F Sharp|F#]]による記述例 ===
<source lang="fsharp">
let isPrime (number : bigint) =
    match number with
        | _ -> seq { bigint(2) .. bigint(sqrt(float number))}
        |> Seq.exists (fun x -> if (number % x = bigint(0)) then true else false)
        |> not

let primes =
    Seq.initInfinite (fun i -> i + 2) //need to skip 0 and 1 for isPrime
    |> Seq.map (fun i -> bigint(i))
    |> Seq.filter isPrime

printfn "%A" primes;;
</source>


===愚直な方法===
<source lang="c">
int isPrime(int n) {
	int i;
	if(n<=1)return 0;
	for(i=2;i<n;i++) {
		if(n%i==0)return 0;
	}
	return 1;
}
</source>
直感的だが、かなり遅い。

===少し改良した方法===
<source lang="c">
int isPrime(int n) {
	int i;
	if(n<=1)return 0;
	for(i=2;i*i<=n;i++) {
		if(n%i==0)return 0;
	}
	return 1;
}
</source>
前のコードを数バイト書き換えただけだが、[[枝刈り]]により比較的速く判定できる。

===高速な方法===
<source lang="c">
/* suppose that a<c */
unsigned long long repeat_add(
		unsigned long long a,unsigned long long b,unsigned long long c) {
	unsigned long long result=0;
	while(b>0) {
		if(b&1) {
			result+=a;
			if(result>=c)result-=c;
		}
		a<<=1;
		if(a>=c)a-=c;
		b>>=1;
	}
	return result;
}

/* suppose that n<(1ll<<63) */
int isPrime(unsigned long long n) {
	const int smallPrimes[30]={
		 2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
		31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
		73, 79, 83, 89, 97,101,103,107,109,113
	};
	int i,ok;
	unsigned long long d,a;
	unsigned long long pownow,powcurrent,powresult;
	int r,s;

	if(n<2)return 0;
	for(i=0;i<30;i++) {
		if(n==smallPrimes[i])return 1;
		if(n%smallPrimes[i]==0)return 0;
	}

	s=0;d=n-1;
	while((d&1)==0){s++;d>>=1;}
	for(i=0;i<12;i++) {
		powresult=1;
		pownow=smallPrimes[i];
		powcurrent=d;
		while(powcurrent>0) {
			if(powcurrent&1)powresult=repeat_add(powresult,pownow,n);
			pownow=repeat_add(pownow,pownow,n);
			powcurrent>>=1;
		}
		if(powresult==1)continue;
		ok=1;
		for(r=0;r<s;r++) {
			if(powresult==n-1){ok=0;break;}
			if(powresult==1)break;
			powresult=repeat_add(powresult,powresult,n);
		}
		if(ok)return 0;
	}
	return 1;
}
</source>
処理は複雑だが、ある程度大きな数でもかなり速く判定してくれる。

===巨大素数の確率的判定===
大きな数が素数であるかを決定的に判定することは難しいが、[[ミラー・ラビン法]]により確率的に判定することができる。
<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>

===優秀なソフトウェアに丸投げする方法===
[[Mathematica]]が使えるなら、<nowiki>PrimeQ[判定したい数]</nowiki>という関数で素数判定ができる。

しかし、調子に乗って<nowiki>PrimeQ[(2^57885161)-1]</nowiki>を計算しようとすると、かなり時間がかかりそうだ。

==素数リストの作成==
ある範囲の素数の一覧が欲しい場合、[[エラトステネスの篩]]を使うと効率よく一覧が作れる。
<source lang="c">
#define PRIME_MAX 1000000

int primeNum;
int primeList[PRIME_MAX];
char primeFlag[PRIME_MAX+1];

void makePrimeList(void) {
	int i,j;
	primeNum=0;
	for(i=0;i<=PRIME_MAX;i++)primeFlag[i]=1;
	primeFlag[0]=primeFlag[1]=0;
	for(i=2;i<=PRIME_MAX;i++) {
		if(primeFlag[i]) {
			primeList[primeNum++]=i;
			for(j=i*2;j<=PRIME_MAX;j+=i)primeFlag[j]=0;
		}
	}
}
</source>
このコードは、1からPRIME_MAXまでの素数のリストを求めるプログラムである。<br />
makePrimeList関数を実行すると、primeNumに素数の数が、primeListに素数のリストが格納される。<br />
また、primeFlagにその数が素数かどうかが格納されるので、O(1)で素数判定ができて便利である。

==素数を使った問題==
[[競技プログラミング]]においては、道順などの求めたい数を100,000,007などの素数で割った数を出力させる問題が出題されることがある。
素数で割った余りを扱う際には、[[逆元]](割る数をMとしたとき、a*xをMで割った余りが1になるようなxをaの逆元という)が簡単に計算できるので、都合がいい。

逆元を求めるのには拡張[[ユークリッドの互除法]]を使うのがいいとされているが、
そんな複雑なことをしなくても累乗の計算をするだけでとりあえず求めることができる。
<source lang="c">
#define MOD_BY 100000007

int getGyakugen(int a) {
	int zyou=MOD_BY-2;
	int result=1;
	while(zyou>0) {
		if(zyou&1)result=(int)((long long)result*a%MOD_BY);
		zyou>>=1;
		a=(int)((long long)a*a%MOD_BY);
	}
	return result;
}
</source>

==関連項目==
* [[東工大生]]
* [[PrimeGrid]]
* [[LLR]]
* [[ミラー・ラビン法]]
* [[エラトステネスの篩]]
* [[AKS素数判定法]]
* [[合成数]]
* [[素因数分解]]