M.Hiroi's Home Page

JavaScript Programming

JavaScript Junk Scripts


Copyright (C) 2025 Makoto Hiroi
All rights reserved.

素数

「素数 (prime number)」は正の約数が 1 と自分自身しかない自然数のことです。1 は素数に含めません。1 でない正整数で、素数ではないものを「合成数 (composite number)」といいます。たとえば、100 以下の素数は 25 個あります。

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

今回は素数にまつわるプログラムを JavaScript で作ってみましょう。詳しい説明は拙作のページをお読みください。

●仕様

●プログラムリスト

//
// primes.js : 素数
//
// Copyright (c) 2025 Makoto Hiroi
//
// Released under the MIT license
// https://opensource.org/license/mit/
//

// べき剰余
function modpow(x, y, n) {
    let r = 1;
    while (y > 0) {
        if ((y & 1) > 0) r = r * x % n;
        x = x * x % n;
        y >>= 1;
    }
    return r;
}

// エラトステネスの篩
function sieve(n) {
    let ps = [2];
    let x = 3;
    let primes = new Uint8Array(Math.floor(n / 2));
    while (x * x <= n) {
        let y = Math.floor((x - 3) / 2);
        if (!primes[y]) {
            ps.push(x);
            y += x;
            while (y < primes.length) {
                primes[y] = 1;
                y += x;
            }
        }
        x += 2;
    }
    while (x <= n) {
        let y = Math.floor((x - 3) / 2);
        if (!primes[y]) ps.push(x);
        x += 2;
    }
    return ps;
}

// 区間篩
function segmented_sieve(low, high) {
    if (low <= 2) return sieve(high)
    let xs = new Uint8Array(high - low + 1);
    for (let p of sieve(Math.floor(Math.sqrt(high)))) {
        let s = Math.ceil(low / p) * p - low;
        if (low <= p) s += p;
        for (let i = s; i < xs.length; i += p) xs[i] = 1;
    }
    let ps = [];
    for (let i = 0; i < xs.length; i++) {
        if (!xs[i]) ps.push(i + low);
    }
    return ps;
}

// 素数の判定 (試し割り)
function is_prime(n) {
    let wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6];
    let [i, s] = [0, 3];
    let p = 2;
    while (p * p <= n) {
        if (n % p == 0) return false;
        p += wheel[i];
        if (++i >= wheel.length) i = s;
    }
    return true;
}


function big_random_sub(n, d) {
    let a = 0n;
    while (d-- > 0) {
        a = (a << 32n) + BigInt(Math.floor(0x1_0000_0000 * Math.random()));
    }
    return a % n;
}

// s 以上 e 未満の乱数を BigInt で生成する
function big_random(s, e) {
    let n = e - s;
    let d = 0;
    for (let m = n; m > 0n; m >>= 32n) d++;
    return s + big_random_sub(n, d);
}

// k 個の乱数を生成
function big_randoms(s, e, k) {
    let n = e - s;
    let xs = [];
    let d = 0;
    for (let m = n; m > 0n; m >>= 32n) d++;
    while (k-- > 0) {
        xs.push(s + big_random_sub(n, d));
    }
    return xs;
}

// ミラーラビン素数判定法

// n を (2 ** c) * d に分解する
function factor2(n) {
    let c = 0n;
    while (n % 2n == 0n) {
        n /= 2n;
        c++;
    }
    return [c, n];
}

// べき剰余 (BigInt 版)
function modpow_bigint(x, y, n) {
    let r = 1n;
    while (y > 0n) {
        if ((y & 1n) > 0n) r = r * x % n;
        x = x * x % n;
        y >>= 1n;
    }
    return r;
}

function miller_rabin_sub(n, s, d, ks) {
    for (let a of ks) {
        if (a >= n) break;
        if (modpow_bigint(a, d, n) != 1n) {
            let r = 0n;
            while (r < s) {
                if (modpow_bigint(a, (2n ** r) * d, n) == n - 1n) break;
                r++;
            }
            if (r >= s) return false;
        }
    }
    return true;
}

function miller_rabin_test(n, k = 20) {
    if (n < 2n) return false;
    if (n == 2n) return true;
    if (n % 2n == 0) return false;
    let [s, d] = factor2(n - 1n);
    let ks;
    if (n < 4759123141n) {
        ks = [2n, 7n, 61n];
    } else if (n <= 0x1_0000_0000_0000_0000n) {
        ks = [2n, 325n, 9375n, 28178n, 450775n, 9780504n, 1795265022n];
    } else {
        ks = big_randoms(1n, n, k);
    }
    return miller_rabin_sub(n, s, d, ks);
}


// リュカ-レーマー・テスト
function lucas_lehmer_test(p) {
    if (p % 2 == 0) return false;
    let m = 2n ** BigInt(p) - 1n;
    let x = 4n;
    for (let i = 0; i < p - 2; i++) {
        x = (x * x - 2n) % m
    }
    return x % m == 0;
}

// 双子素数
function twin_primes(n) {
    let xs = sieve(n + 2);
    let ys = [];
    for (let i = 0; i < xs.length - 1; i++) {
        if (xs[i] + 2 == xs[i + 1]) ys.push([xs[i], xs[i+1]]);
    }
    return ys;
}

export {
    modpow, sieve, segmented_sieve, is_prime, big_random, big_randoms,
    miller_rabin_test, lucas_lehmer_test, twin_primes
}

●簡単なテスト

//
// test_primes.js : 素数のテスト
//
// Copyright (c) 2025 Makoto Hiroi
//
// Released under the MIT license
// https://opensource.org/license/mit/
//
import * as P from './primes.js';

console.log('---- modepow -----');
console.log(P.modpow(7, 654321, 10));   // 7
console.log(P.modpow(7, 6543210, 10));  // 9
console.log('---- sieve, segmented_sieve -----');
console.log(P.sieve(100));         // 2,3,5,7,...
console.log(P.sieve(100).length);  // 25
console.log(P.sieve(100000).length);  // 9592
console.log(P.sieve(1000000).length);  // 78498
console.log(P.sieve(10000000).length);  // 664579
console.log(P.segmented_sieve(100, 200)); // 101 ... 199
console.log(P.segmented_sieve(100, 200).length); // 21
console.log(P.segmented_sieve(1000000, 1000100)); // [1000003, 1000033, 1000037, 1000039, 1000081, 1000099]
console.log('---- is_prime -----');
console.log(P.is_prime(2));  // true
console.log(P.is_prime(97));  // true
console.log(P.is_prime(111));  // false
console.log(P.is_prime(2 * 31 - 1));  // true
console.log(P.is_prime(2 * 32 - 1));  // false
console.log('---- miller_rabin_test -----');
console.log(P.miller_rabin_test(2n));  // true
console.log(P.miller_rabin_test(97n));  // true
console.log(P.miller_rabin_test(111n));  // false
console.log(P.miller_rabin_test(2n ** (2n ** 0n) + 1n));  // true
console.log(P.miller_rabin_test(2n ** (2n ** 1n) + 1n));  // true
console.log(P.miller_rabin_test(2n ** (2n ** 2n) + 1n));  // true
console.log(P.miller_rabin_test(2n ** (2n ** 3n) + 1n));  // true
console.log(P.miller_rabin_test(2n ** (2n ** 4n) + 1n));  // true
console.log(P.miller_rabin_test(2n ** (2n ** 5n) + 1n));  // false
console.log(P.miller_rabin_test(2n ** 31n - 1n));  // true
console.log(P.miller_rabin_test(2n ** 32n - 1n));  // false
console.log(P.miller_rabin_test(2n ** 89n - 1n));  // true
console.log(P.miller_rabin_test(2n ** 90n - 1n));  // false
console.log(P.miller_rabin_test(2n ** 607n - 1n));  // true
console.log(P.miller_rabin_test(2n ** 608n - 1n));  // false
console.log('---- lucas_lehmer_test -----');
for (let i = 3; i < 1000; i += 2) {
    if (P.lucas_lehmer_test(i)) console.log(i);
}
console.log('---- twin_primes -----');
console.log(P.twin_primes(100));
console.log(P.twin_primes(3000000).length); // 20932
console.log('---- big_random -----');
console.log(2n**32n);
for (let i = 0; i < 4; i++) console.log(P.big_random(0n, 2n ** 32n));
console.log(P.big_randoms(2n**16n, 2n ** 32n, 4));
console.log(2n**64n);
for (let i = 0; i < 4; i++) console.log(P.big_random(0n, 2n ** 64n));
console.log(P.big_randoms(2n**48n, 2n ** 64n, 4));
console.log(2n**96n);
for (let i = 0; i < 4; i++) console.log(P.big_random(0n, 2n ** 96n));
console.log(P.big_randoms(2n**80n, 2n ** 96n, 4));

●実行例

$ node test_primes.js
---- modepow -----
7
9
---- sieve, segmented_sieve -----
[
   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
]
25
9592
78498
664579
[
  101, 103, 107, 109, 113,
  127, 131, 137, 139, 149,
  151, 157, 163, 167, 173,
  179, 181, 191, 193, 197,
  199
]
21
[ 1000003, 1000033, 1000037, 1000039, 1000081, 1000099 ]
---- is_prime -----
true
true
false
true
false
---- miller_rabin_test -----
true
true
false
true
true
true
true
true
false
true
false
true
false
true
false
---- lucas_lehmer_test -----
3
5
7
13
17
19
31
61
89
107
127
521
607
---- twin_primes -----
[
  [ 3, 5 ],   [ 5, 7 ],
  [ 11, 13 ], [ 17, 19 ],
  [ 29, 31 ], [ 41, 43 ],
  [ 59, 61 ], [ 71, 73 ]
]
20932
---- big_random -----
4294967296n
445165036n
67746803n
3111766799n
317258130n
[ 1243407988n, 368579738n, 2794824690n, 489971548n ]
18446744073709551616n
14754816802008137100n
9315702261375939981n
12849699641133499619n
6410765238664596166n
[
  7159733700449148631n,
  14753923102134523062n,
  9021503864372089985n,
  495178360824506896n
]
79228162514264337593543950336n
45897100895909178701239571339n
30936128865905323611882010292n
26503787494083486129373503926n
15771361232900844706876833616n
[
  52136448438308779376695304250n,
  60460452991954675642761860585n,
  48552450636253818932392439895n,
  35381167006374857514585795697n
]

初版 2025 年 2 月 15 日