素数
素数判定編集
F#による記述例編集
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;;
愚直な方法編集
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;
}
直感的だが、かなり遅い。
少し改良した方法編集
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;
}
前のコードを数バイト書き換えただけだが、枝刈りにより比較的速く判定できる。
高速な方法編集
/* 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;
}
処理は複雑だが、ある程度大きな数でもかなり速く判定してくれる。
巨大素数の確率的判定編集
大きな数が素数であるかを決定的に判定することは難しいが、ミラー・ラビン法により確率的に判定することができる。
#!/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;
}
優秀なソフトウェアに丸投げする方法編集
Mathematicaが使えるなら、PrimeQ[判定したい数]という関数で素数判定ができる。
しかし、調子に乗ってPrimeQ[(2^57885161)-1]を計算しようとすると、かなり時間がかかりそうだ。
素数リストの作成編集
ある範囲の素数の一覧が欲しい場合、エラトステネスの篩を使うと効率よく一覧が作れる。
#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;
}
}
}
このコードは、1からPRIME_MAXまでの素数のリストを求めるプログラムである。
makePrimeList関数を実行すると、primeNumに素数の数が、primeListに素数のリストが格納される。
また、primeFlagにその数が素数かどうかが格納されるので、O(1)で素数判定ができて便利である。
素数を使った問題編集
競技プログラミングにおいては、道順などの求めたい数を100,000,007などの素数で割った数を出力させる問題が出題されることがある。 素数で割った余りを扱う際には、逆元(割る数をMとしたとき、a*xをMで割った余りが1になるようなxをaの逆元という)が簡単に計算できるので、都合がいい。
逆元を求めるのには拡張ユークリッドの互除法を使うのがいいとされているが、 そんな複雑なことをしなくても累乗の計算をするだけでとりあえず求めることができる。
#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;
}