素数判定 (primality test) とは、自然数 (natural number) nが素数 (prime number) であるか合成数 (composite number) であるかを判定する問題です。素数判定を行うアルゴリズムを素数判定法といい、nが「合成数である」または「素数ではないかと思われる数＝確率的素数 (PRP: probable prime) 」と判定する確率的素数判定法 (probabilistic primality test) と、nが「素数である」か否かを判定する決定的素数判定法 (deterministic primality test) に大別されます。
ここでは、幾つかの素数判定法とその応用例を集めてみました。素因数分解については "PrimeFactor!" をご覧ください。
ソースコードでは関数リテラルとしていますが、関数宣言(function文)でも同様に動作します。
[ 戻る ]
試し割り法 (trial division) とは、最も単純な決定的素数判定法で、自然数nを小さい順に割ってみて、割り切れるかどうかを調べていく手法です。割り切れる (divisible) ときは「合成数」、最後まで割り切れない (indivisible) ときは「素数」と確実に判定できますが、nが大きいときは時間がかかってしまいます。
単純に2からn−1までの数で割ったときの時間計算量 (time complexity) はΟ(n)となります。nがpで割り切れたときの商 (quotient) をqとすると、n=pq (2≤p≤√n≤q<n) の関係が成り立ちますので、√nまでの数で確かめれば十分です。これにより計算量をΟ(√n)に減少させることができます。
さらに、既に確かめた数の倍数 (multiple) については確かめる必要はないため、√nまでの素数で割るようにすると計算量の係数を小さくすることができます。このため、2と3以外の素数は6の倍数の前後 (6i±1) にしか存在しないことを利用すると計算量を1/3に減らすことができます。
var trial_division = function(n) {
/* 試し割り法(trial division)による素数判定 */
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
for (var divisor = 5; divisor * divisor <= n; divisor += 6) {
if (n % divisor == 0) return false;
if (n % (divisor + 2) == 0) return false;
}
return true;
};2以上9007199254740992以下の自然数nに対応しています。
[ 戻る ]
エラトステネスの篩 (Sieve of Eratosthenes) は自然数n以下の全ての素数を列挙するための決定的素数判定法の1つで、時間計算量はΟ(n log log n)になります。基本的なアルゴリズムは以下の通りです。
このとき、倍数の探索に除算を使わず、配列の添字 (subscript / index) に素数を順次加算していくことで探索を行い、消去フラグを立てています。
試し割り法と同様に2と3以外は 6i±1 (i=1, 2, …) だけを調べるようにすると効率がよくなります。そこで、配列pを作成し、p[0]に2、p[1]に3を代入し、p[2]以降に素数の候補として 6i±1 の数をnに達するまで代入していきます。
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| 00 | 2 | 3 | 5 | 7 | 11 | 13 | 17 | 19 | 23 | 25 |
| 10 | 29 | 31 | 35 | 37 | 41 | 43 | 47 | 49 | 53 | 55 |
| 20 | 59 | 61 | 65 | 67 | 71 | 73 | 77 | 79 | 83 | 85 |
| 30 | 89 | 91 | 95 | 97 | 101 | 103 | 107 | 109 | 113 | 115 |
| 40 | 119 | 121 | 125 | 127 | 131 | 133 | 137 | 139 | 143 | 145 |
上表より配列p上で、5、7、11、…の倍数の出現する奇数回目と偶数回目の間隔とその和を求めてみます。
| m | i | 奇数回目 | 偶数回目 | 間隔の和 |
|---|---|---|---|---|
| 5 | 2 | 7 | 3 | 10 |
| 7 | 3 | 9 | 5 | 14 |
| 11 | 4 | 15 | 7 | 22 |
| 13 | 5 | 17 | 9 | 26 |
| 17 | 6 | 23 | 11 | 34 |
これより2回分の間隔の和は 2m、偶数回目の間隔は2i−1になることがわかります。このことを利用して、i=2から順に最小の数miを残してその倍数になる要素p[j]に0を代入していきます。これをmiが ⌊√n⌋に達するまで繰り返します。
| i | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| 00 | 2 | 3 | 5 | 7 | 11 | 13 | 17 | 19 | 23 | 0 |
| 10 | 29 | 31 | 0 | 37 | 41 | 43 | 47 | 0 | 53 | 0 |
| 20 | 59 | 61 | 0 | 67 | 71 | 73 | 0 | 79 | 83 | 0 |
| 30 | 89 | 0 | 0 | 97 | 101 | 103 | 107 | 109 | 113 | 0 |
| 40 | 0 | 0 | 0 | 127 | 131 | 0 | 137 | 139 | 0 | 0 |
filter()メソッドを使って配列pから0の要素を削除した配列を新たに作成し、それを戻り値とします。
var sieve = function(n) {
/* Eratosthenesの篩(sieve of Eratisthenes) */
var p = [];
var i, j, len;
var i0, i1; // 篩内のiの増分算出用
var m = 5, m1 = 2;
// p[i]に値を代入
if (n > 1) p[0] = 2;
if (n > 2) p[1] = 3;
if (n < 5) return p;
for (i = 2; m <= n; i++) {
// p[i]にm=6i±1を代入
p[i] = m;
m += m1;
m1 = 6 - m1; // mの増分は交互に2と4
}
// 篩
len = p.length;
m1 = int_root(p[len - 1]); // 計算の上限値
for (i = 2; i < len; i++) {
m = p[i];
if (m > m1) break;
if (m == 0) continue;
i0 = m + m; // 増分の基準値
i1 = i0 - i - i + 1; // iの増分
for (j = i + i1; j < len; j += i1) {
p[j] = 0; // mの倍数(合成数)ならばp[i]に0を代入
i1 = i0 - i1;
}
}
return p.filter(Boolean); // 0の要素を削除
};134217728 (=227) 以下の自然数nに対応しており、指定されたnまでの素数の個数［素数計数関数 (prime counting function)］Π(n) と最大1万個の素数を列挙しています。素数の個数が1万個以上のときは小さい素数は表示されません。nが31以上のときはスクロールバーを操作して結果をご確認ください。nに大きい値を指定すると処理に時間がかかります。
[ 戻る ]
フェルマーテスト (Fermat primality test) とは、自然数nが確率的素数であるかどうかの判定を行うアルゴリズムで、フェルマーの小定理 (Fermat's little theorem) の対偶 (contraposition) を利用しています。
フェルマーの小定理とは、pを素数、aを任意の自然数としたとき ap≡a (mod p)、特にaとpが互いに素 (coprime) であるとき ap−1≡1 (mod p) が成り立つという定理です。
後者の対偶を取ると、自然数nと互いに素な自然数aが an−1≢1 (mod n) を満たすとき、nは合成数になるという十分条件 (sufficient condition) が得られます。即ち an−1≡1 (mod n) である場合、nは素数である可能性があります。このことを用いて判定を行います。
素数ではないにもかかわらず「確率的素数」と判定されることもあり、このような合成数を擬素数 (psp: pseudoprime) といいます。「確率的素数」と判定されたときは、異なるaを用いて複数回のテストを行うことで素数である可能性を高めることができます。それでも、nと互いに素である全てのaでフェルマーテストを通過してしまう擬素数が存在し、それらはカーマイケル数 (Carmichael number, OEIS A002997) と呼ばれます。
var fermat_test = function(n, a) {
/* フェルマーテストによる確率的素数判定 */
if (n < 2) return false;
if (a < 2 || a == n) return '--'; // 判定不能
if (euclidean(n, a) > 1) return false;
return mod_pow(a, n - 1, n) == 1;
};94906265(226.5)<n≤134217728(=227)では冪剰余 (modular exponentiation) 演算で偶数に丸められてしまう丸め誤差 (round-off error / rounding error) により、素数がうまく判定できないこともあります。
[ 戻る ]
ミラー・ラビン素数判定法 (Miller–Rabin primality test) は、フェルマーテストを改良した確率的素数判定法の1つです。
自然数nが奇素数 (odd prime) であるとき、n-1は偶数 (even number) なので、2s×tと表すことができます。ここで、sは自然数、tは奇数 (odd number) です。このとき、2≤a≤n−1を満たすaについて at≡1 (mod n)、または0≤i≤s−1に対して a2i·t≡−1 (mod n) が成立します。このときのnは底 (base) aの強い確率的素数 (a-SPRP: strong probable prime) になります。これを用いた判定の流れは:
奇数nが合成数であると正しく判別できる基数 (base) aは合成性の証人 (witness of compositeness) であり、1〜n−1の中に少なくとも75%は存在します。ところが、合成数なのに確率的素数と判定されることがあり、このときのnは底aについての強擬素数 (a-spsp: strong pseudoprime) で、aは証人にならない強い嘘つき (strong liar) となります。このため独立に選んだaで判定を繰り返します。
判定対象となるnの範囲を限定し、基底aをとして小さい素数を組み合わせることで決定的素数判定法として用いることができます。ここでのnの上限値は各々の全てのaに対する最小の強擬素数 (a-spsp, OEIS A014233) になります。
演算にIEEE754型倍精度浮動小数点数 (double precision floating point number) を使用していて、冪剰余演算内で偶数に丸められても一の位がかろうじて保持される値である227=134217728をnの上限としてaを選びました。
var miller_rabin = function(n) {
/* ミラー・ラビン素数判定法 */
if (n % 2 == 0 || n < 2) return n == 2; // 偶数 or 2未満で素数は2のみ
var s = 0, t = n - 1;
var a = [2, 3, 5, 7];
var len = a.length;
while (t % 2 == 0) {
// n-1=s*2^t に分解
s += 1;
t /= 2;
}
var witness = function(ai) {
// 合成性の証人(true:合成数)
var r = mod_pow(ai, t, n);
if (r == 1 || r == n - 1) return false;
for (var i = 0; i < s; i++) {
r = r * r % n;
if (r == n - 1) return false;
}
return true;
};
for (var i = 0; i < len; i++) {
if (a[i] >= n) break;
if (witness(a[i])) return false;
}
return true;
};94906265(226.5)<n≤134217728(=227)では冪剰余演算で偶数に丸められてしまう丸め誤差の影響により、素数がうまく判定きないこともあります。
[ 戻る ]
エラトステネスの篩の考え方を応用すると、自然数n1からn2までの約数の個数を列挙することができます。
n1<n2であることを確認してから、長さn2−n1、先頭要素d[0]がn1となる配列dを作成します。開始値となるn1以上で最小の i の倍数を ⌊(n1+i−1)÷i⌋×iで求め、配列d内の i の倍数になる要素に個数を加算していきます。このとき、自然数が i で割り切れたときの商も約数になりますので+2を、自然数が i の平方数のときは+1を加算します。これを1から√n2までの自然数 i で繰り返します。
var divisors = function(n1, n2) {
/* n1からn2までの約数の個数 */
var d = [];
var j, s;
if (n1 > n2) {
s = n1;
n1 = n2;
n2 = s;
}
for (var i = 0; i <= n2 - n1; i++){
d[i] = 0;
}
for (i = 1; i * i <= n2; i++){
s = Math.floor((n1 + i - 1) / i) * i;
for (j = s; j <= n2; j += i){
d[j - n1] += i * i < j
? 2
: i * i == j ? 1 : 0;;
}
}
return d;
};9007199254740992 (=253) までの自然数に対応しています。大きな値を指定すると処理に時間がかかり、値の範囲を広くすると"Out of Memory"エラーが生じることがあります。
[ 戻る ]
エラトステネスの篩では素数であるかどうかを配列に格納しますが、2からnまでの自然数の最小素因数 (smallest prime factor: SPF / least prime factor: lpf, OEIS A020639 ) を格納した配列を得ることもできます。ここで得られた最小素因数を用いて素因数分解をすることができます。
最小素因数を格納する配列spの初期値としてその数自身にしておきます。次に2の倍数（偶数）については2を代入します。これらを1つに纏めるとループの回数を減らすことができます。
3以上√n以下の奇数 i については、すでに i 未満の数が格納されていたら無視し、そうでなければ全ての i の倍数に最小素因数として i を代入します。このとき、i2未満の i の倍数にはすでに i 未満の数が記録されていますので、i2から2i 毎に i を代入していきます。
var spf = function(n) {
/* 2からnまでの自然数の最小素因数 */
var sp = [0];
var j
for (var i = 1; i <= n; i++){
sp[i] = i % 2 ? i : 2;
}
for (i = 3; i * i <= n; i += 2){
if (sp[i] < i) continue;
for (j = i * i; j <= n; j += i + i){
if (sp[j] < j) continue;
sp[j] = i;
}
}
return sp;
};ここで得られた最小素因数sp[]を用いて素因数分解をすることができます。基本的には現在の数を最小素因数で割ることを繰り返し、最後に1になったときに終了します。whileループ中の(prime[sp[n]]||0)+1は論理和演算子 (logical OR operator) の左辺がfalseと評価されるときのみ右辺が評価される短絡評価 (short-circuit evaluation) を利用し、prime[sp[n]]がundefinedのときは1を、それ以外のときは現在の値に1を加算した値を返しています。
var prime_factor_spf = function(n, sp) {
/* 最小素因数sp[]を用いた素因数分解 */
var prime = {};
while (sp[n] > 1) {
prime[sp[n]] = (prime[sp[n]] || 0) + 1;
n /= sp[n];
}
return prime;
};2以上1000000 (=106) 以下の自然数に対応しています。
最小素因数を用いた素因数分解は最小素因数sp[]を保持するためメモリ消費量が多いことと前処理の篩処理に要する時間により、大きな数に対して適切な方法とは言えません。一度に多く (≈104個) の小さい数 (≲106) の素因数分解を処理するのに適しています。
[ 戻る ]
素数日 (prime day) とは日付を年月日で並べたときにできる5桁以上の数が素数となる日のことで、娯楽数学 (entertainment math) の1つになります。ここでは日付を表す文字列ymdから区切り文字を除去して作成した数値を上述のミラー・ラビン法を用いて素数判定をしています。
var ymd = '2022-04-17';
var IsPrimeDay = miller_rabin(ymd.replace(/[-/.年月日]/g,''));9491年〜9999年では判定漏れが起こることがあります。
[ 戻る ]
指定した年yearの素数日pd[]を列挙します。素数判定には上述のミラー・ラビン法を用いています。
var PrimeDay = function(year) {
/* 指定した年の素数日の列挙 */
var pd = []; // 年間素数日
var i = 0; // 配列pdの添字
var yymm;
var dim; // 月の日数(day in month)
for (var month = 1; month <= 12; month++) {
yymm = year * 10000 + month * 100;
dim = new Date(year, month, 0).getDate(); // 月の日数の取得
for (var dd = 1; dd <= dim; dd += 2) { // 日が偶数の場合を除く
if (miller_rabin(yymm + dd)) pd[i++] = yymm + dd;
}
}
return pd;
};9491年〜9999年では判定漏れが起こることがあります。
[ 戻る ]
[ 戻る ]
[ 戻る ]
