前の記事で「素数検索」「素数判定」「素因数分解」の64bit版のプログラムを作成しましたが、ここでは 128bit で計算し、数字39桁 (安全に扱えるのは38桁) まで完全に利用可能なコードを考えてみぬくました。一般によく使われているPC等 (OSや言語) は、64bitまでの整数の扱いが限界ですが、Cなどの処理系では128bit整数の計算が可能です。
Perl (64bit版)でも、Math::BigInt パッケージで64bitを超える整数を扱えるのですが、処理が遅くなります。そこで、Inline::C を使った実装によって、Cの符号なし128bit整数演算を使い、モンゴメリ乗算と組み合わせれば、完全な128bit整数の扱いが可能になります。あくまでMPUモジュールの優秀さに対抗して、興味の探求から始まったプログラム習作ですが、128bitまでフルに使えれば実務的にも有用かなと思います。
C側で128bit計算はできますが、Perl側は整数型としては認識出来ない大きさですので、Cへの引数や結果の素数値は、Perl側では文字列として扱います。そのため、C側で数字列を128bit整数に変換し、結果の128bit整数も数字列に変換してPerlの配列構造にセットします。この入り口と出口のデータ型のパース処理だけやれば、後はC内部の128bit計算をオーバーフローに気をつけて実現すれば良いことになります。128bit版で作成したプログラム (関数) は以下の3種類です。「何に使うの？」と聞かれれば、研究者でもないので特に用途はないのですが、単純に「アルゴリズム的にはどうなの？」と思っただけなんですけどね😅。
- 素数の判定
64bit版の素数判定を128bit整数版に拡張しました。64bit未満の数値はミラー・ラビン法 (確定)。64bit以上の数値は、BPSW (ミラー・ラビン基数2＋強Lucasテスト) を使った高速判定です。 - 素数の検索
64bit版では、mod30 および mod210 の車輪篩 (wheel sieve) を使い、広いレンジの検索では高速でしたが、メモリの制限もあり100億程度が限界でした (100億のフルレンジ検索で約15GBのメモリを使用)。そこで、128bit版では狭い範囲の検索に特化し、窓篩＋個別判定の方式に変えました。この方が大きな数字でピンポイントに素数を検索できますので、より実務的と判断しました。 - 素因数分解
これも64bit版からの拡張です。論理も同じく「ポラード・ロー素因数分解法」(ブレント版) を使っています。素数判定ロジックも自前のものを組み込んでいます。また実験的ですが、合成数1つあたりの分解の反復回数 (上限) を指定できるようにしました。
※本稿で掲載の自作のプログラムコード、および動作確認の結果等については、一切の保証はありません。ご了承の上でご確認お願いいたします。《かなり真面目な記事になりました。5万字を超えています！😲》
1．素数の判定
使い方は64bit版同様に簡単です。調べたい数字列を Inline::C の関数 (is_prime128) に引数で渡すだけです。結果は [ 1: 素数、0: 合成数 ] で返します。64bit未満の数値 (264-1) までは、ミラー・ラビン法で決定的・数学的に確定した判定になります。64bit以上の大きな数値は、BPSW (Baillie-PSW) のミラー・ラビン基数2＋強Lucasテストで判定しています。
Lucasテストについては、64bitまでなら不要です。ミラー・ラビン法の7個の (基数) witness (2, 325, 9375, 28178, 450775, 9780504, 1795265022) は、264未満の全ての整数で判定が正しいことが、計算機による全数検証で確かめられていると理解しています (反例が見つかっていないのではなく、検証が済んでいる範囲です)。BPSW が必要になるのは、264 を超える範囲を扱う場合です。
264 超では、確定できる固定 witness 集合は知られていません。標準はBPSWで反例は見つかっていませんが、証明はされていません。64bitでは Lucasテストは不要ですが、128bitでは必要になります。その結果は「ほぼ確実」になります。なお、128bitの乗算は256bit積が必要になります。そのため、128bit版の各素数プログラムでは、オーバーフロー対策で統一してモンゴメリ乗算を使っています。また素数判定ロジックには、小さな素数での試し割りと平方数の検査を含めています。
プログラムコードの検証では、(AI) Claude の協力で以下のテストにパスしています (誤判定なし)。
- ミラー・ラビン基数2：
65bit〜128bitの2万件。 - 強 Lucas テスト：
65bit〜128bitの2万件。小さい奇数10万件 (101〜120万、擬素数が多い領域)。既知の強Lucas擬素数12個は正しく判定しました。素数3000件で偽陰性なし。 - 関数全体の試験：
18579件 (素数1901個)。2128 直前の素数3000個全部。64bit境界 (264 ±300)、2127 ±300、メルセンヌ素数 (289−1、2107−1、2127−1)、素数の平方、半素数を含みます。
この素数判定関数は、素数検索や素因数分解でも利用しますので、速度よりも正確な判定が必須なのです。
※実行速度の検証は使用マシンの関係で、Windows 11 環境＋Strawberry Perl を使っています。GMP (GNU Multiple Precision Arithmetic Library) なしのフォールバックなどの関係で、処理プロセスに一定のコストがあるかもしれません。Linux等のOSでGMPありの環境で計測すれば、MPUとの速度差が無くなる可能性があります。
※実行環境 (使用マシンの性能やOSの種類等) により実行時間は異なります。時間はあくまで目安です。
【 判定速度の比較 (素数) 】時間単位：ミリ秒 (ms)
| No. | 判定する数 (素数) | Inline::C (128bit) |
MPU is_prime |
| 1 | 18446744073709551557 264 未満の最大素数 |
0.005 ms | 0.008 ms |
| 2 | 18446744073709551629 264 超の最小素数 |
0.007 ms | 0.048 ms |
| 3 | 618970019642690137449562111 メルセンヌ素数 289 -1 |
0.008 ms | 0.023 ms |
| 4 | 162259276829213363391578010288127 メルセンヌ素数 2107 -1 |
0.009 ms | 0.023 ms |
| 5 | 170141183460469231731687303715884105727 メルセンヌ素数 2127 -1 |
0.009 ms | 0.024 ms |
| 6 | 340282366920938463463374607431768211297 128bit未満で最大の素数 2128 -159 |
0.013 ms | 0.094 ms |
| 7 | 963892709107619255203 素数 70bit |
0.008 ms | 0.049 ms |
| 8 | 837929339126532632949089 素数 80bit |
0.010 ms | 0.058 ms |
| 9 | 761410555838264820800620043 素数 90bit |
0.010 ms | 0.142 ms |
| 10 | 709333186601603974736322613229 素数 100bit |
0.011 ms | 0.124 ms |
| 11 | 837856659879287811671387066351609 素数 110bit |
0.011 ms | 0.085 ms |
| 12 | 665437864705649490934581250675891361 素数 120bit |
0.012 ms | 0.074 ms |
| 13 | 89894305334003814839566387884490424443 素数 127bit |
0.012 ms | 0.122 ms |
【 判定速度の比較 (合成数) 】時間単位：ミリ秒 (ms)
| No. | 判定する数 (合成数) | Inline::C (128bit) |
MPU is_prime |
| 1 | 18446744073709551615 264 -1 |
0.003 ms | 0.005 ms |
| 2 | 18446744073709551617 264 +1 |
0.007 ms | 0.026 ms |
| 3 | 170141183460469231731687303715884105729 2127 +1 |
0.004 ms | 0.007 ms |
| 4 | 340282366920938463463374607431768211455 2128 -1 |
0.003 ms | 0.007 ms |
| 5 | 1114830017192323777514569 素数の平方 (p 約40bit) |
0.005 ms | 0.028 ms |
| 6 | 38494287782063039473137139965219900721 素数の平方 (p 約63bit) |
0.006 ms | 0.029 ms |
| 7 | 35102443323622819405329210077365945517 半素数 63bit x 63bit |
0.006 ms | 0.029 ms |
| 8 | 194190127723507801860413143590896545919 半素数 64bit x 64bit |
0.006 ms | 0.028 ms |
| 9 | 26726860543977851007246381436349 素数 x 37 (試し割りで除外) |
0.004 ms | 0.017 ms |
| 10 | 26832559556697902539588317647777 素数 x 41 (試し割りをすり抜ける) |
0.005 ms | 0.018 ms |
| 11 | 1296001987165015643369032371289 カーマイケル数 (Chernick k=1000000511) |
0.006 ms | 0.029 ms |
| 12 | 1296002356525428293844563788009 カーマイケル数 (Chernick k=1000000606) |
0.011 ms | 0.044 ms |
| 13 | 1296003188558614941390839227921 カーマイケル数 (Chernick k=1000000820) |
0.006 ms | 0.028 ms |
| 14 | 2417851664969925135785653 合成なのに基数2の強確率素数 (p*(2p-1)) |
0.009 ms | 0.041 ms |
| 15 | 2417851671619771500544621 合成なのに基数2の強確率素数 (p*(2p-1)) |
0.009 ms | 0.042 ms |
【 判定速度について 】
全ての速度比較に於いて、MPU::is_prime を超える高速な素数判定になっています。素数・合成数いずれも13マイクロ秒以下で判定しています。ただしMPUとの差が純粋なアルゴリズムの差とは限りません。MPUは任意の桁を扱えるので、その分の汎用性が時間に出ている可能性があります。
素数判定は桁が増える程、少しずつ遅くなります。 0.005 → 0.013 ms と 64bit から 128bit にかけて、ゆるやかに増えています。判定1回の中の乗算が、64bit → 128bit で重くなる分です。合成数は、ほとんどが 0.003〜0.009 ms です。 試し割りや、ミラー・ラビン基数2だけで弾けるため、素数よりずっと軽く済みます。強Lucasテストが効いたのは、最後の「基数2の強確率素数」2件です。 ミラー・ラビン基数2では素数と判定されますが、強Lucasテストで合成数と判定されて、0.009 ms かかっています。ほかの合成数より少し重いのはそのためです。
【 プログラムソース 】
このプログラムのソースは以下のとおりです。Windows 環境でも「Strawberry Perl for Windows」で動作確認済みです (Inline::C は CPAN からインストールする必要があります)。
# =====================================================================
# (Perl＋Inline::C) 128bit版 素数判定
#
# N < 2^64 : ミラー・ラビン法 (決定的・数学的に確定)
# N >= 2^64 : BPSW (ミラー・ラビン基数2＋強Lucasテスト)
#
# 数の最大値: 340282366920938463463374607431768211455 (2^128-1,39桁)
#
# 関数の引数: SV* is_prime128(char *s) ※判定値は数字列で渡す
# 関数戻り値: 1: 素数、0: 合成数 (不正な数字列はcroakで中断する)
#
# This Script is in the public domain, No rights reserved.
# Script written by N.O, Updated Last on 2026/09/30
# Claude (AI) assisted with the script optimization.
# =====================================================================
use strict;
use Inline (
C => 'DATA',
DIRECTORY => '_Inline',
);
# ---------------------------------------------
# テストサンプル (コマンドラインで数値指定)
# ---------------------------------------------
use Math::Prime::Util qw(is_prime);
use Time::HiRes qw(gettimeofday);
my ($n, $xs) = @ARGV;
$n =~ s/[,_]*//go;
if ($n eq "" || $n eq "0" || $n =~ /\D/o) {
print "\n[HELP] >perl $0 判定数 [MPU指定]\n\n";
print "判定数は 340282366920938463463374607431768211455 以内の正の整数です。\n\n";
print "MPU指定:'0'以外を入力なら MPU::is_prime を使用します。\n";
exit;
}
my($prime, $start, $finish);
if ($xs) {
$start = gettimeofday();
$prime = is_prime($n);
$finish = gettimeofday();
print "\n*** 素数判定(MPU::is_prime)の結果 ***\n\n";
} else {
$start = gettimeofday();
$prime = is_prime128($n);
$finish = gettimeofday();
print "\n*** 素数判定 (128bit版) の結果 ***\n\n";
}
my $elapse = int(($finish - $start) * 1_000_000) / 1000;
if ($prime) {
print "判定数: $n は素数です\n\n";
} else {
print "判定数: $n は素数ではありません\n\n";
}
print "実行時間: ", sprintf("%.3f", $elapse), " ミリ秒(ms)\n";
exit;
# -----------------------------------
# ここからInline::C(言語)での記述
# -----------------------------------
__DATA__
__C__
/* ==============================================================
* 128bit版 素数判定: BPSW (ミラー・ラビン基数2＋強Lucasテスト)
* SV* normalize128(char *s) 数字列 <=> 128bit整数変換の確認用
* ============================================================== */
#include <stdint.h>
#include <math.h>
typedef uint64_t u64;
typedef unsigned __int128 u128;
/* 数字列 <-> 128bit整数の変換。成功:0, 失敗:-1 */
static int parse_u128(const char *s, u128 *out) {
const u128 MAXV = ~(u128)0;
u128 v = 0;
int nd = 0;
while (*s == ' ' || *s == '\t') s++;
for (; *s >= '0' && *s <= '9'; s++) {
unsigned dg = (unsigned)(*s - '0');
if (v > (MAXV - dg) / 10) return -1; /* v*10+dg が 2^128 を超える */
v = v * 10 + dg;
nd++;
}
while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') s++;
if (nd == 0 || *s != '\0') return -1;
*out = v;
return 0;
}
/* buf は 40 バイト以上 (39桁＋NUL) */
static void u128_to_str(u128 v, char *buf) {
char tmp[40];
int i = 0, j = 0;
if (v == 0) {
buf[0] = '0'; buf[1] = '\0'; return;
}
while (v) {
tmp[i++] = (char)('0' + (int)(v % 10)); v /= 10;
}
while (i) buf[j++] = tmp[--i];
buf[j] = '\0';
}
/* 128bit 補助関数 */
static int ctz128(u128 x) { /* x != 0 */
u64 lo = (u64)x;
return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll((u64)(x >> 64));
}
static int bitlen128(u128 x) { /* x != 0 */
u64 hi = (u64)(x >> 64);
return hi ? 128 - __builtin_clzll(hi) : 64 - __builtin_clzll((u64)x);
}
static u128 isqrt128(u128 n) { /* floor(sqrt(n)) */
if (n < 2) return n;
u128 x = (u128)1 << ((bitlen128(n) + 1) / 2); /* sqrt(n) 以上の初期値 */
for (;;) {
u128 y = (x + n / x) >> 1;
if (y >= x) return x;
x = y;
}
}
/* ヤコビ記号 (a/n) n は奇数の正の数。戻り値 -1, 0, 1 */
static int jacobi128(u128 a, u128 n) {
int result = 1;
a %= n;
while (a) {
while ((a & 1) == 0) {
a >>= 1;
unsigned r = (unsigned)(n & 7);
if (r == 3 || r == 5) result = -result;
}
u128 t = a; a = n; n = t;
if ((a & 3) == 3 && (n & 3) == 3) result = -result;
a %= n;
}
return n == 1 ? result : 0;
}
/* 小さな符号付き整数 x の n による剰余 [0, n) */
static u128 small_mod(long x, u128 n) {
u128 m = (u128)(x < 0 ? -(long long)x : x) % n;
return x < 0 ? (m ? n - m : 0) : m;
}
/* 64bit: モンゴメリ乗算＋ミラー・ラビン (n < 2^64 で決定的) */
static u64 mont_inv64(u64 n) {
u64 x = n;
for (int i = 0; i < 5; i++) x *= 2 - n * x;
return x;
}
static u64 mont_mul64(u64 a, u64 b, u64 n, u64 ninv) {
u128 t = (u128)a * b;
u64 q = (u64)t * ninv;
u64 hi = (u64)(t >> 64);
u64 mh = (u64)(((u128)q * n) >> 64);
return hi >= mh ? hi - mh : hi - mh + n;
}
static const u64 SMALL_PRIMES[] = {2,3,5,7,11,13,17,19,23,29,31,37};
static int is_prime_u64(u64 n) {
if (n < 2) return 0;
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return n == SMALL_PRIMES[i];
}
if (n < 37ULL * 37) return 1;
static const u64 b2[] = {31, 73};
static const u64 b3[] = {2, 7, 61};
static const u64 b4[] = {2, 13, 23, 1662803};
static const u64 b7[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
const u64 *bases; int nb;
if (n < 9080191ULL) { bases = b2; nb = 2; }
else if (n < 4759123141ULL) { bases = b3; nb = 3; }
else if (n < 1122004669633ULL) { bases = b4; nb = 4; }
else { bases = b7; nb = 7; }
u64 ninv = mont_inv64(n);
u64 one = (0 - n) % n;
u64 mone = n - one;
u64 r2 = (u64)(((u128)one * one) % n);
u64 d = n - 1;
int s = __builtin_ctzll(d);
d >>= s;
for (int i = 0; i < nb; i++) {
u64 a = bases[i] % n;
if (a == 0) continue;
u64 b = mont_mul64(a, r2, n, ninv);
u64 x = one;
for (u64 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul64(x, b, n, ninv);
b = mont_mul64(b, b, n, ninv);
}
if (x == one || x == mone) continue;
int composite = 1;
for (int r = 1; r < s; r++) {
x = mont_mul64(x, x, n, ninv);
if (x == mone) {
composite = 0;
break;
}
}
if (composite) return 0;
}
return 1;
}
/* 128bit: モンゴメリ乗算 (R = 2^128, n は奇数) */
/* 128bit×128bit → 256bit (hi:lo) */
static void mul_128x128(u128 a, u128 b, u128 *hi, u128 *lo) {
u64 a0 = (u64)a, a1 = (u64)(a >> 64), b0 = (u64)b, b1 = (u64)(b >> 64);
u128 p00 = (u128)a0 * b0, p01 = (u128)a0 * b1;
u128 p10 = (u128)a1 * b0, p11 = (u128)a1 * b1;
u128 mid = (p00 >> 64) + (u64)p01 + (u64)p10; /* 3*2^64 未満なので溢れない */
*lo = ((u128)(u64)mid << 64) | (u64)p00;
*hi = p11 + (p01 >> 64) + (p10 >> 64) + (mid >> 64);
}
/* n^-1 mod 2^128 (ニュートン法) 有効bit: 3→6→12→24→48→96→192 */
static u128 mont_inv128(u128 n) {
u128 x = n;
for (int i = 0; i < 6; i++) x *= 2 - n * x;
return x;
}
/* a*b*R^-1 mod n (a, b < n) */
static u128 mont_mul128(u128 a, u128 b, u128 n, u128 ninv) {
u128 hi, lo, mh, dummy;
mul_128x128(a, b, &hi, &lo);
u128 q = lo * ninv; /* mod 2^128 で巻き戻る */
mul_128x128(q, n, &mh, &dummy); /* q*n の上位128bit */
return hi >= mh ? hi - mh : hi - mh + n;
}
static u128 add_mod128(u128 a, u128 b, u128 n) {
return a >= n - b ? a - (n - b) : a + b;
}
static u128 sub_mod128(u128 a, u128 b, u128 n) {
return a >= b ? a - b : n - (b - a);
}
/* x / 2 mod n (n は奇数) */
static u128 half_mod128(u128 x, u128 n) {
return (x & 1) == 0 ? x >> 1 : (x >> 1) + (n >> 1) + 1;
}
/* 強確率素数テスト: n を基数 a で判定。ninv/one/r2 はモンゴメリ準備値 */
static int mr128(u128 n, u128 a, u128 ninv, u128 one, u128 r2) {
u128 nm1 = n - 1;
int s = ctz128(nm1);
u128 d = nm1 >> s;
a %= n;
if (a == 0) return 1;
u128 mone = n - one;
u128 b = mont_mul128(a, r2, n, ninv);
u128 x = one;
for (u128 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul128(x, b, n, ninv);
b = mont_mul128(b, b, n, ninv);
}
if (x == one || x == mone) return 1;
for (int r = 1; r < s; r++) {
x = mont_mul128(x, x, n, ninv);
if (x == mone) return 1;
}
return 0;
}
/* 強LucasTest (Selfridgeの方法A: D = 5, -7, 9, -11,...P = 1, Q = (1-D)/4) */
/* n: 奇数。n+1 が溢れないこと (n != 2^128-1) は呼び出し側の試し割りで保証 */
static int strong_lucas128(u128 n, u128 ninv, u128 one, u128 r2) {
u128 rt = isqrt128(n);
if (rt * rt == n) return 0; /* 平方数は D が見つからない */
long D = 5;
for (int guard = 0; guard < 100000; guard++) {
u128 dm = small_mod(D, n);
int j = jacobi128(dm, n);
if (j == -1) break;
if (j == 0) return n == (u128)(D < 0 ? -D : D); /* 共通因数あり(nが|D|と等しい場合だけ素数) */
D = D > 0 ? -(D + 2) : -(D - 2);
}
long Q = (1 - D) / 4;
u128 Dm = mont_mul128(small_mod(D, n), r2, n, ninv);
u128 Qm = mont_mul128(small_mod(Q, n), r2, n, ninv);
u128 np1 = n + 1;
int s = ctz128(np1);
u128 d = np1 >> s;
u128 U = one, V = one, Qk = Qm; /* k = 1: U=1, V=P=1, Q^1 */
for (int i = bitlen128(d) - 2; i >= 0; i--) {
U = mont_mul128(U, V, n, ninv); /* U_2k = U_k V_k */
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n); /* V_2k = V_k^2 - 2Q^k */
Qk = mont_mul128(Qk, Qk, n, ninv);
if ((d >> i) & 1) { /* k -> k+1 */
u128 U2 = half_mod128(add_mod128(U, V, n), n);
u128 V2 = half_mod128(add_mod128(mont_mul128(Dm, U, n, ninv), V, n), n);
U = U2; V = V2;
Qk = mont_mul128(Qk, Qm, n, ninv);
}
}
if (U == 0 || V == 0) return 1;
for (int r = 1; r < s; r++) {
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n);
Qk = mont_mul128(Qk, Qk, n, ninv);
if (V == 0) return 1;
}
return 0;
}
/* 素数判定 (128bit) */
static int is_prime_u128(u128 n) {
if ((n >> 64) == 0) return is_prime_u64((u64)n); /* 64bit以下は確定判定 */
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return 0; /* n > 2^64 なので n 自身が素数ではない */
}
u128 ninv = mont_inv128(n);
u128 one = (0 - n) % n; /* R mod n */
u128 r2 = one; /* R^2 mod n : R を128回倍加 (256bit除算を避ける) */
for (int i = 0; i < 128; i++) r2 = add_mod128(r2, r2, n);
if (!mr128(n, 2, ninv, one, r2)) return 0;
return strong_lucas128(n, ninv, one, r2);
}
/* Perlへの公開関数 */
int is_prime128(char *s) {
u128 n;
if (parse_u128(s, &n) != 0)
croak("is_prime128: 不正な数字列です (数字のみ最大 340282366920938463463374607431768211455)");
return is_prime_u128(n);
}
/* 数字列を128bit整数に変換して数字列に戻す(変換動作確認用) */
SV* normalize128(char *s) {
u128 n;
char buf[40];
if (parse_u128(s, &n) != 0)
croak("normalize128: 不正な数字列です (数字のみ最大 340282366920938463463374607431768211455)");
u128_to_str(n, buf);
return newSVpv(buf, 0);
}
2．素数の検索
128bit版では、大きな数を扱えること、メモリの制限に掛からず素速く結果を返せることを設計方針にしました。64bit版ではワイドレンジで素数一覧を出すのに向いたアルゴリズム (mod30とmod210による車輪篩 wheel sieve) を使いましたが、128bit版は「窓篩＋個別素数判定」の方式にしています。
検索の仕方は、最小数～最大数の指定と、要求数 (欲しい素数の個数) を指定すると、1個または最大100個までの素数を検索して返す仕様です。
- 最小数 (min) は指定必須
0以上 2128 -1 (39桁) までの数。ただし、2128 -159 が128bitで最大の素数ですので、これを超える数を指定しても検索結果は 0個 になります。 - 最大数 (max) は指定または省略可能
指定の場合：min～max までの素数を検索します。
省略の場合：min 以上の素数を要求数まで検索します。 - 要求数 (count) は指定または省略可能
指定の場合：min～max 条件の検索から要求数までを返します。
省略の場合：
⇨ max が省略の場合：1個
⇨ max が指定の場合：最大100個 (ただし検索条件の結果内)
最大数と要求数は 0 を指定すると、省略と同じ意味になります。要求数の指定は最大100です。最小数～最大数の指定は、それ自身の判定も含みます。P の次の素数が欲しい場合は、最小数に P＋1 を指定します。
このプログラムは広範囲の素数一覧出しには向いていませんが、128bitの範囲で最大100個までの素数を高速に検索できます。そのための安全策として、最大100個までの検索結果の限定と、最小数からの探索範囲をクラメール予想の100倍を限度として制限しています (ただし最小探索範囲は 100,000 です)。
窓篩 (window sieve) と個別判定の仕組み
- 最小数の次の奇数から、固定サイズの窓 (奇数の候補を最大16384個) を作ります。
- 小さい素数 (窓の大きさの16倍まで最大65536) の倍数を窓の中で消します。
- 生き残った候補だけを素数判定します。
- 個数が揃うか範囲の末尾に達したら終わりです。
窓の中の処理は、数が大きくなっても同じ作業量です。メモリも窓の分だけで、範囲の広さに依存しません。篩う素数の上限は、窓の大きさの2〜512倍で比べました。大きくしすぎると、準備のほうが重くなって遅くなります。8〜32倍程度が最速だったので、中間の16倍にしています。実質は「篩で候補を絞ってから、残りを個別に判定する」方式です。2128 付近の素数を篩だけで確定させるには、√(1038) ≒ 1019 の素数が必要になって不可能です。篩の部分と個別判定の部分を、それぞれ数えてみたのが以下の表です。
実測した内訳（要求数100を指定）
| 大きさ | 窓の候補数 | 生き残った数 (判定した数) |
素数1個あたり 判定回数 |
うち篩の部分の 処理時間の割合 |
| 40bit | 約1800 | 約150 | 1.5回 | 35% |
| 64bit | 約2900 | 約230 | 2.3回 | 22% |
| 100bit | 約4500 | 約350 | 3.5回 | 8% |
| 127bit | 約5800 | 約445 | 4.5回 | 6.5% |
篩で約9割が消えます。 篩を使わずに全部判定すると、2127 付近では素数1個につき約44個の奇数を調べることになります。それが約4.5回で済んでいます (約10分の1)。処理の大半は個別判定です。127bitでは全体の94%が判定 (BPSW)、篩は6.5%程度です。また、大きい数ほど素数1個あたりの判定が増えます。 素数の密度が下がるためです。
【 窓篩の仕組みを簡単に図解 】
候補が篩で絞られて、残りが判定にかかる流れは以下のようになります。
1．小さい素数で倍数を消す (30個→12個)
2．生き残りだけを1個ずつ判定 (12回)
図の読み方
- 上段：3, 5, 7, 11, 13 の倍数を窓から消した結果です。灰色のマスの下の ÷3 などは、最初に消した素数です。30個のうち18個が消えて、白いマス12個が残ります。
- 下段：残った12個だけを、素数判定にかけた結果です。9個が素数 (緑)、3個が合成数 (赤) でした。赤のマスの 17×59 などは、素因数分解の中身です。
- 赤い3個の理由：篩に使った素数が小さいので (13まで) 、17 や 19 を因数に持つ合成数が生き残ったためです。篩で消し残した分を、素数判定が最終的に取り除いています。
実際のコードとの違い
図は仕組みを見やすくするための例です。実際の動作とは次の点が違います。
- 窓の大きさ：実際は数千個 (最大16384個) の奇数を並べます。
- 篩に使う素数：実際は 3 から 65536 までの素数 (最大 65521) を使います (窓の大きさの16倍まで)。そのため、生き残りの大半が素数になります。2127付近でも、生き残りの約10%を判定にかけるだけで済みます。
- 素数判定：生き残りが、実際に篩で使った最大の素数の平方根の範囲に収まる小さな数では、判定せずに篩を通っただけで素数と確定します。
【 プログラムコードの検証 】
(AI) Claude の協力で、独立した正解 (gmpy2) を使って以下のことを確認しました。
| テスト項目 | 件数 | 誤り件数 |
| 小さい範囲 (min = 0〜399、区間) | 2600 | 0 |
| 小〜中のランダム (個数1〜150) | 1700 | 0 |
| 境界 (264、2127、2128 直前) | 120 | 0 |
| ランダム 65bit〜128bit | 1500 | 0 |
| ランダム区間 (最大100個) | 500 | 0 |
| max / min が素数ちょうど | 1500 | 0 |
| 篩だけで確定する境界の前後 | 140 | 0 |
| 固定サンプルテスト | 20 | 0 |
【 素数検索の速度 】時間単位：ミリ秒 (ms)
※MPU (primes) は、min, max を必ず指定での検索結果のため参考値。
| No. | 最小数 最大数 |
要求数 | 摘要・結果 | Inline::C (128bit) |
MPU primes |
| 1 | 0 0 |
5 | 最初の素数5個 2, 3, 5, 7, 11 |
0.021 ms | 0.007 ms |
| 2 | 2 0 |
0 | 次の素数 (min自身) 1個の結果 |
0.021 ms | 0.005 ms |
| 3 | 14 0 |
0 | 14の次の素数 17 |
0.027 ms | 0.005 ms |
| 4 | 1000 0 |
0 | 1000の次の素数 1009 |
0.028 ms | 0.006 ms |
| 5 | 1000 1100 |
5 | 1000〜1100 5個の結果 |
0.028 ms | 0.007 ms |
| 6 | 1000 1100 |
0 | 1000〜1100 16個の結果 |
0.035 ms | 0.007 ms |
| 7 | 90 96 |
0 | 素数がない区間 0 個 |
0.028 ms | 0.005 ms |
| 8 | 4294967291 4294967311 |
0 | 232 付近 2個の結果 |
0.046 ms | 0.009 ms |
| 9 | 18446744073709551557 0 |
0 | 264 未満の最大素数 (min自身) 1個の結果 |
0.030 ms | 0.008 ms |
| 10 | 18446744073709551558 0 |
0 | 264 を越える最初の素数 1個の結果 |
0.038 ms | 104.255 ms |
| 11 | 18446744073709551616 0 |
3 | 264 から3個 3個の結果 |
0.052 ms | 105.227 ms |
| 12 | 1000000000000000000 0 |
5 | 1029 から5個 5個の結果 |
0.045 ms | 0.014 ms |
| 13 | 170141183460469231731687303715884105727 0 |
3 | 2127 -1 (素数) から3個 3個の結果 |
0.063 ms | 105.011 ms |
| 14 | 340282366920938463463374607431768211297 0 |
0 | 2128 未満の最大の素数 (min自身) 1個の結果 |
0.042 ms | 103.461 ms |
| 15 | 340282366920938463463374607431768211298 0 |
0 | それより上は素数なし 0個の結果 |
0.058 ms |
103.435 ms |
| 16 | 340282366920938463463374607431768211455 0 |
0 | 2128 -1 (素数なし) 0個の結果 |
0.023 ms | 105.197 ms |
| 17 | 340282366920938463463374607431768210955 340282366920938463463374607431768211455 |
0 | 2128 直前の区間 6個の結果 |
0.226 ms | 106.031 ms |
| 18 | 18446744073709551557 340282366920938463463374607431768211455 |
0 | 264 未満の最大素数～ 2128 -1 までの100個 |
0.662 ms | 応答なし |
※MPU (primes) は、64bit未満では高速ですが、64bit以上の数では極端に悪くなります。これはテスト環境 (OS) のためかもしれません。GMPありのLinux等で計測すれば、本来の実行時間になる可能性があります。
MPU (primes) がクラッシュする件：
上記 No.18 の範囲指定では、MPU (primes) が Windows環境では「応答なし」になり、Claude と原因分析しました。GMPありのLinux環境で、primes 単体で実験してみても、約13秒後に「Segmentation fault」でクラッシュします。範囲が広くなっても、落ちるまでの時間は一定 (約13秒) です。つまり「巨大すぎる範囲」だけの問題ではなく、指定上限が約1021 ( 270 ) で既に落ちるようです (No.18 の指定上限は 2128 -1 )。原因までは不明です。実行中のメモリ(RSS)は、約13秒間13MBのまま増えませんでした。メモリ不足が原因ではないようです。なお、264 をまたぐ幅が143と約4.8万の範囲では正常に終了しました。検索結果も、次の1個と1068個を正しく応答しました。
今回作成した Inline::C (128bit版) では、検索結果の100個制限や検索区間 (幅) をクラメール予想の100倍を限度としているため、巨大な範囲を指定しても安全に短時間で応答できます。
【 プログラムソース 】
このプログラムのソースは以下のとおりです。Perlからの引数では「最大数」「要求数」を省略可能にするために、Inline Stack から直に引数を取り出し、返り値 (素数の配列リファレンス) も Inline Stack にプッシュする方法にしています。プログラム間のリンケージ規約に直接触るので、行儀の悪い書き方です。公開関数はプログラム間のインターフェイスのみで、本来の関数の実体 (主エントリ) は primes128_core です。
# ======================================================================
# (Perl＋Inline::C) 128bit版 素数検索
#
# 素数の検索: 窓篩＋素数判定
# 素数判定等: (BPSW)ミラー・ラビン法＋強Lucasテスト＋モンゴメリ乗算
#
# 数の最大値: 340282366920938463463374607431768211455 (2^128-1,39桁)
#
# 関数の引数: SV* primes128(char *min, char *max, UV count)
# 関数戻り値: 昇順の素因数"文字列"の配列リファレンスを返す
#
# min : min以上の素数を順に探索し数字列の配列リファレンスで返す
# max : 探索の上限 (その値を含む)。"0"なら上限なし
# count : 要求数 (欲しい個数：最大100。101以上の結果は100に丸める)
# 0 なら既定値 (max 指定あり→100個、max が"0"→1個)。
#
# max が"0"のとき、探索は min から (log min)^2 の100倍(最小10万)の
# 距離で打ち切る (クラメール予想の100倍)。打ち切りまでに見つかった
# 素数を返す。2^128未満の素数は 2^128-159 が最大。それを超える範囲
# では、結果が 0個 の場合がある。
#
# This Script is in the public domain, No rights reserved.
# Script written by N.O, Updated Last on 2026/10/05
# Claude (AI) assisted with the script optimization.
# ======================================================================
use strict;
use Inline (
C => 'DATA',
DIRECTORY => '_Inline',
);
# ---------------------------------------------
# テストサンプル (コマンドラインで数値指定)
# ---------------------------------------------
use Math::Prime::Util qw(primes);
use Time::HiRes qw(gettimeofday);
my ($min, $max, $cnt) = @ARGV;
$min =~ s/[,_]*//go;
$max =~ s/[,_]*//go;
$max ||= '0'; $cnt ||= 0;
if ($min eq "" || $min =~ /\D/o || $max =~ /\D/o || $cnt < 0 || $cnt > 100) {
print "\n[HELP] >perl $0 最小数 [最大数 [要求数]]\n\n";
print "最小数は39桁以内の正の整数。最大数と要求数は'0'または省略できます。\n";
exit;
}
my $start = gettimeofday();
my $primes = primes128($min, $max, $cnt);
my $finish = gettimeofday();
my $elapse = int(($finish - $start) * 1_000_000) / 1000;
my $pieces = @$primes;
print "\n*** 素数の検索 (128bit版) の結果 ***\n\n";
print "最小数: ", $min, "\n最大数: ", $max, "\n要求数: ", $cnt, "\n\n";
print "見つかった素数: ", $pieces, " 個\n\n";
if ($pieces > 0) {
if ($pieces <= 10) {
print join(",", @$primes), "\n\n";
} else {
print "最初の5個: ", join(",", @{$primes}[0..5]), "\n";
print "最後の5個: ", join(",", @{$primes}[-5..-1]), "\n\n";
}
}
print "実行時間: ", sprintf("%.3f", $elapse), " ミリ秒(ms)\n";
exit;
# -----------------------------------
# ここからInline::C(言語)での記述
# -----------------------------------
__DATA__
__C__
/* ======================================================== *
* 128bit版 素数検索: 窓篩＋素数判定(ミラー・ラビン/BPSW) *
* ======================================================== */
#include <stdint.h>
#include <math.h>
typedef uint64_t u64;
typedef unsigned __int128 u128;
/* 数字列 <-> 128bit整数の変換。成功: 0, 失敗: -1 */
static int parse_u128(const char *s, u128 *out) {
const u128 MAXV = ~(u128)0;
u128 v = 0;
int nd = 0;
while (*s == ' ' || *s == '\t') s++;
for (; *s >= '0' && *s <= '9'; s++) {
unsigned dg = (unsigned)(*s - '0');
if (v > (MAXV - dg) / 10) return -1; /* v*10+dg が 2^128 を超える */
v = v * 10 + dg;
nd++;
}
while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') s++;
if (nd == 0 || *s != '\0') return -1;
*out = v;
return 0;
}
/* buf は 40 バイト以上 (39桁 + NUL) */
static void u128_to_str(u128 v, char *buf) {
char tmp[40];
int i = 0, j = 0;
if (v == 0) { buf[0] = '0'; buf[1] = '\0'; return; }
while (v) { tmp[i++] = (char)('0' + (int)(v % 10)); v /= 10; }
while (i) buf[j++] = tmp[--i];
buf[j] = '\0';
}
/* 128bit 補助関数 */
static int ctz128(u128 x) { /* x != 0 */
u64 lo = (u64)x;
return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll((u64)(x >> 64));
}
static int bitlen128(u128 x) { /* x != 0 */
u64 hi = (u64)(x >> 64);
return hi ? 128 - __builtin_clzll(hi) : 64 - __builtin_clzll((u64)x);
}
static u128 isqrt128(u128 n) { /* floor(sqrt(n)) */
if (n < 2) return n;
u128 x = (u128)1 << ((bitlen128(n) + 1) / 2); /* sqrt(n) 以上の初期値 */
for (;;) {
u128 y = (x + n / x) >> 1;
if (y >= x) return x;
x = y;
}
}
/* ヤコビ記号 (a/n) n は奇数の正の数。戻り値 -1, 0, 1 */
static int jacobi128(u128 a, u128 n) {
int result = 1;
a %= n;
while (a) {
while ((a & 1) == 0) {
a >>= 1;
unsigned r = (unsigned)(n & 7);
if (r == 3 || r == 5) result = -result;
}
u128 t = a; a = n; n = t;
if ((a & 3) == 3 && (n & 3) == 3) result = -result;
a %= n;
}
return n == 1 ? result : 0;
}
/* 小さな符号付き整数 x の n による剰余 [0, n) */
static u128 small_mod(long x, u128 n) {
u128 m = (u128)(x < 0 ? -(long long)x : x) % n;
return x < 0 ? (m ? n - m : 0) : m;
}
/* 64bit: モンゴメリ乗算 + ミラー・ラビン (n < 2^64 で決定的) */
static u64 mont_inv64(u64 n) {
u64 x = n;
for (int i = 0; i < 5; i++) x *= 2 - n * x;
return x;
}
static u64 mont_mul64(u64 a, u64 b, u64 n, u64 ninv) {
u128 t = (u128)a * b;
u64 q = (u64)t * ninv;
u64 hi = (u64)(t >> 64);
u64 mh = (u64)(((u128)q * n) >> 64);
return hi >= mh ? hi - mh : hi - mh + n;
}
static const u64 SMALL_PRIMES[] = {2,3,5,7,11,13,17,19,23,29,31,37};
static int is_prime_u64(u64 n) {
if (n < 2) return 0;
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return n == SMALL_PRIMES[i];
}
if (n < 37ULL * 37) return 1;
static const u64 b2[] = {31, 73};
static const u64 b3[] = {2, 7, 61};
static const u64 b4[] = {2, 13, 23, 1662803};
static const u64 b7[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
const u64 *bases; int nb;
if (n < 9080191ULL) { bases = b2; nb = 2; }
else if (n < 4759123141ULL) { bases = b3; nb = 3; }
else if (n < 1122004669633ULL) { bases = b4; nb = 4; }
else { bases = b7; nb = 7; }
u64 ninv = mont_inv64(n);
u64 one = (0 - n) % n;
u64 mone = n - one;
u64 r2 = (u64)(((u128)one * one) % n);
u64 d = n - 1;
int s = __builtin_ctzll(d);
d >>= s;
for (int i = 0; i < nb; i++) {
u64 a = bases[i] % n;
if (a == 0) continue;
u64 b = mont_mul64(a, r2, n, ninv);
u64 x = one;
for (u64 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul64(x, b, n, ninv);
b = mont_mul64(b, b, n, ninv);
}
if (x == one || x == mone) continue;
int composite = 1;
for (int r = 1; r < s; r++) {
x = mont_mul64(x, x, n, ninv);
if (x == mone) { composite = 0; break; }
}
if (composite) return 0;
}
return 1;
}
/* 128bit: モンゴメリ乗算 (R = 2^128, n は奇数) */
/* 128bit×128bit → 256bit (hi:lo) */
static void mul_128x128(u128 a, u128 b, u128 *hi, u128 *lo) {
u64 a0 = (u64)a, a1 = (u64)(a >> 64), b0 = (u64)b, b1 = (u64)(b >> 64);
u128 p00 = (u128)a0 * b0, p01 = (u128)a0 * b1;
u128 p10 = (u128)a1 * b0, p11 = (u128)a1 * b1;
u128 mid = (p00 >> 64) + (u64)p01 + (u64)p10; /* 3*2^64 未満なので溢れない */
*lo = ((u128)(u64)mid << 64) | (u64)p00;
*hi = p11 + (p01 >> 64) + (p10 >> 64) + (mid >> 64);
}
/* n^-1 mod 2^128 (ニュートン法: 有効bit 3→6→12→24→48→96→192) */
static u128 mont_inv128(u128 n) {
u128 x = n;
for (int i = 0; i < 6; i++) x *= 2 - n * x;
return x;
}
/* a*b*R^-1 mod n (a, b < n) */
static u128 mont_mul128(u128 a, u128 b, u128 n, u128 ninv) {
u128 hi, lo, mh, dummy;
mul_128x128(a, b, &hi, &lo);
u128 q = lo * ninv; /* mod 2^128 で巻き戻る */
mul_128x128(q, n, &mh, &dummy); /* q*n の上位128bit */
return hi >= mh ? hi - mh : hi - mh + n;
}
static u128 add_mod128(u128 a, u128 b, u128 n) {
return a >= n - b ? a - (n - b) : a + b;
}
static u128 sub_mod128(u128 a, u128 b, u128 n) {
return a >= b ? a - b : n - (b - a);
}
/* x / 2 mod n (n は奇数) */
static u128 half_mod128(u128 x, u128 n) {
return (x & 1) == 0 ? x >> 1 : (x >> 1) + (n >> 1) + 1;
}
/* 強確率素数テスト: n を基数 a で判定。ninv/one/r2 はモンゴメリ準備値 */
static int mr128(u128 n, u128 a, u128 ninv, u128 one, u128 r2) {
u128 nm1 = n - 1;
int s = ctz128(nm1);
u128 d = nm1 >> s;
a %= n;
if (a == 0) return 1;
u128 mone = n - one;
u128 b = mont_mul128(a, r2, n, ninv);
u128 x = one;
for (u128 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul128(x, b, n, ninv);
b = mont_mul128(b, b, n, ninv);
}
if (x == one || x == mone) return 1;
for (int r = 1; r < s; r++) {
x = mont_mul128(x, x, n, ninv);
if (x == mone) return 1;
}
return 0;
}
/* 強Lucasテスト(Selfridge の方法A: D = 5, -7, 9, -11,... P = 1, Q = (1-D)/4) */
/* n: 奇数。n+1 が溢れないこと (n != 2^128-1) は呼び出し側の試し割りで保証 */
static int strong_lucas128(u128 n, u128 ninv, u128 one, u128 r2) {
u128 rt = isqrt128(n);
if (rt * rt == n) return 0; /* 平方数は D が見つからない */
long D = 5;
for (int guard = 0; guard < 100000; guard++) {
u128 dm = small_mod(D, n);
int j = jacobi128(dm, n);
if (j == -1) break;
if (j == 0) return n == (u128)(D < 0 ? -D : D); /* 共通因数あり(n が|D|と等しい場合だけ素数) */
D = D > 0 ? -(D + 2) : -(D - 2);
}
long Q = (1 - D) / 4;
u128 Dm = mont_mul128(small_mod(D, n), r2, n, ninv);
u128 Qm = mont_mul128(small_mod(Q, n), r2, n, ninv);
u128 np1 = n + 1;
int s = ctz128(np1);
u128 d = np1 >> s;
u128 U = one, V = one, Qk = Qm; /* k = 1: U=1, V=P=1, Q^1 */
for (int i = bitlen128(d) - 2; i >= 0; i--) {
U = mont_mul128(U, V, n, ninv); /* U_2k = U_k V_k */
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n); /* V_2k = V_k^2 - 2Q^k */
Qk = mont_mul128(Qk, Qk, n, ninv);
if ((d >> i) & 1) { /* k -> k+1 */
u128 U2 = half_mod128(add_mod128(U, V, n), n);
u128 V2 = half_mod128(add_mod128(mont_mul128(Dm, U, n, ninv), V, n), n);
U = U2; V = V2;
Qk = mont_mul128(Qk, Qm, n, ninv);
}
}
if (U == 0 || V == 0) return 1;
for (int r = 1; r < s; r++) {
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n);
Qk = mont_mul128(Qk, Qk, n, ninv);
if (V == 0) return 1;
}
return 0;
}
/* 素数判定 (128bit) */
static int is_prime_u128(u128 n) {
if ((n >> 64) == 0) return is_prime_u64((u64)n); /* 64bit以下は確定判定 */
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return 0; /* n > 2^64 なので n 自身が素数ではない */
}
u128 ninv = mont_inv128(n);
u128 one = (0 - n) % n; /* R mod n */
u128 r2 = one; /* R^2 mod n : R を128回倍加 (256bit除算を避ける) */
for (int i = 0; i < 128; i++) r2 = add_mod128(r2, r2, n);
if (!mr128(n, 2, ninv, one, r2)) return 0;
return strong_lucas128(n, ninv, one, r2);
}
/* 素数検索 (窓篩) */
#define PL_BMAX 65536 /* 篩に使う素数の上限 (準備用の表) */
#define PL_WMAX 16384 /* 1つの窓に入れる奇数の候補の最大数 */
#define PL_BFACT 16 /* 篩う素数の上限 = 窓の候補数×この値(PL_BMAXで頭打ち) */
static unsigned int plist[6600]; /* 3 以上 PL_BMAX 以下の奇素数 */
static int plist_n = 0;
/* lim 以下の奇素数を out[] に昇順で入れ個数を返す (out は 6600個必要) */
static int odd_primes_upto(unsigned lim, unsigned *out) {
unsigned char comp[PL_BMAX / 2 + 1]; /* comp[i]=1 : 2i+1 は合成数 */
int n = 0;
memset(comp, 0, lim / 2 + 1);
for (unsigned i = 3; i * i <= lim; i += 2)
if (!comp[i >> 1]) for (unsigned j = i * i; j <= lim; j += 2 * i) comp[j >> 1] = 1;
for (unsigned i = 3; i <= lim; i += 2)
if (!comp[i >> 1]) out[n++] = i;
return n;
}
/* 文字列に変換してPerl配列にpush */
static void push_str(AV *av, u128 v) {
char buf[40];
u128_to_str(v, buf);
av_push(av, newSVpv(buf, 0));
}
/* 公開関数の実体 */
static SV* primes128_core(const char *smin, const char *smax, UV count_in) {
const u128 MAXV = ~(u128)0;
u128 mn, mx;
if (parse_u128(smin, &mn) != 0 || parse_u128(smax, &mx) != 0)
croak("primes128_c: 不正な数字列です (数字のみ最大 340282366920938463463374607431768211455)");
int limited = (mx != 0);
if (limited && mx < mn) croak("primes128_c: max が min より小さいです");
unsigned count = count_in ? (count_in > 100 ? 100u : (unsigned)count_in) : (limited ? 100u : 1u);
u128 end = mx; /* 探索範囲の末尾 (この値を含む) */
if (!limited) { /* クラメールの予想の100倍で打ち切る */
double ln = mn > 3 ? log((double)mn) : 1.0;
double dist = 100.0 * ln * ln;
if (dist < 100000.0) dist = 100000.0;
u128 cap = (u128)dist;
end = (MAXV - mn < cap) ? MAXV : mn + cap;
}
AV *av = newAV();
unsigned found = 0;
if (mn <= 2 && end >= 2) {
push_str(av, 2); found++;
}
if (found < count) {
u128 lo = (mn <= 3) ? 3 : ((mn & 1) ? mn : mn + 1); /* 最初の奇数の候補 */
double lnlo = lo > 3 ? log((double)lo) : 1.0;
int w0 = (int)(count * lnlo * 0.65 + 64); /* 窓の候補数: 必要な範囲の約1.3倍 */
if (w0 > PL_WMAX) w0 = PL_WMAX;
/* 篩に使う奇素数の表 */
unsigned beff_max = (unsigned)((unsigned long long)w0 * PL_BFACT > PL_BMAX ? PL_BMAX : (unsigned long long)w0 * PL_BFACT);
unsigned int plist[6600];
int plist_n = odd_primes_upto(beff_max, plist);
unsigned char mark[PL_WMAX]; /* mark[i]=1 : lo+2i は合成数 */
while (found < count && lo <= end) {
u128 avail = (end - lo) / 2 + 1; /* lo, lo+2,...end 以下の個数 */
int w = (avail < (u128)w0) ? (int)avail : w0;
memset(mark, 0, (size_t)w);
unsigned long long bl = (unsigned long long)w * PL_BFACT;
unsigned beff = bl > PL_BMAX ? PL_BMAX : (unsigned)bl;
unsigned pmax = 1; /* 実際に使った最大の素数 */
for (int k = 0; k < plist_n; k++) {
unsigned p = plist[k];
if (p > beff) break;
pmax = p;
unsigned r = (unsigned)(lo % p);
/* lo+2*i0 が p の倍数になる最小の i0 */
unsigned i0 = (unsigned)((((u64)(p - r)) % p) * ((p + 1) / 2) % p);
if (i0 >= (unsigned)w) continue;
if (lo + 2 * (u128)i0 == p) i0 += p; /* p 自身は消さない */
for (unsigned i = i0; i < (unsigned)w; i += p) mark[i] = 1;
}
/* x < sq で篩を通れば素数が確定 */
u128 sq = (u128)(pmax + 1) * (pmax + 1);
for (int i = 0; i < w; i++) {
if (mark[i]) continue;
u128 x = lo + 2 * (u128)i;
if (x < sq || is_prime_u128(x)) {
push_str(av, x);
if (++found >= count) break;
}
}
if ((u128)w == avail) break; /* 探索範囲の末尾に達した */
lo += 2 * (u128)w;
}
}
/* 配列リファレンスで返す */
return newRV_noinc((SV*)av);
}
/* Perl への公開関数 (maxとcountは省略可) */
void primes128(...) {
Inline_Stack_Vars;
const char *smin;
const char *smax = "0";
UV count = 0;
if (Inline_Stack_Items < 1)
croak("Usage: primes128(min [, max [, count]])");
smin = SvPV_nolen(Inline_Stack_Item(0));
if (Inline_Stack_Items > 1 && SvOK(Inline_Stack_Item(1))) smax = SvPV_nolen(Inline_Stack_Item(1));
if (Inline_Stack_Items > 2 && SvOK(Inline_Stack_Item(2))) count = SvUV(Inline_Stack_Item(2));
SV *rv = primes128_core(smin, smax, count);
Inline_Stack_Reset;
Inline_Stack_Push(sv_2mortal(rv));
Inline_Stack_Done;
}
3．素因数分解
この関数も使い方は64bit版同様に簡単です。調べたい数字列を Inline::C の関数 (factor128) に引数で渡すだけです。結果は、素因数の配列リファレンスで返します。引数には、合成数1つあたりの分解の反復回数 (上限) も指定できます。指定がない場合には既定の反復回数になります。
素因数分解には、ポラード・ロー素因数分解法 (Pollard's rho algorithm) を使っています。ただし、素因数分解は素数判定ほど簡単には128bit化できません。素数判定は BPSW で十分に実用的ですが、素因数分解 (ポラード・ロー法) は、小さい方の素因数 に対して約 回の反復が必要です。38桁の半素数で両方の因数が大きい場合は、現実的な時間では終わらない可能性があります。
そこで、安全のために反復回数に上限を付けて、見つからなければ「未分解の合成数」として返す方式にしました。関数の引数に反復回数を指定できるようにしたので、分解できなかった場合には、人手で任意の回数を指定して再試行できます。現在の反復回数 (上限の既定値) は、134,217,728 回 (約1億3千万回。最大で4～6秒程度？) です。
素数判定は、先に作成の自前の関数を組み込んでいます。264 以下の部分は、64bitモンゴメリ＋ブレント法をそのまま使い、必ず完全に分解できます。264 超の部分は、128bitモンゴメリ＋ブレント法＋GCDまとめ取りです。GCDは128bit対応のバイナリGCDです。前処理として、1000未満の素因数は先に試し割りで取り除いています。平方数は平方根を取って2回分解します。また、完全冪の検出も入れてあります (3, 5, 7, 11乗の検出)。
今回作成の関数は、ポラード・ロー法のみ適用ですが、その結果、下記の速度比較 No.25 のようなケースでは未分解で終わります。このアルゴリズムでは分解の反復回数が足りないのが理由ですが、MPUは処理時間のバラつきはありますが、必ず分解しています。これが不思議だったので、MPUのドキュメントで調べたところ、MPUの factor は、複数のアルゴリズムを試しながら素因数分解しているとのことでした。
- Math::Prime::Util (64bit以下を担当)
少しの試し割り、完全冪の検出。その後で、ポラード・ロー、SQUFOF、p−1法の組み合わせ。見つけた非素数の因子ごとに、この組み合わせを適用します。 - Math::Prime::Util::GMP (大きな数を担当)
試し割り、完全冪の検出、ポラード・ロー、p−1法、HOLF (フェルマー法の一種)、ECM (楕円曲線法)、QS (二次篩) を試しながら分解します。
ドキュメントには、26桁の半素数なら100ms未満、36桁の半素数なら1秒未満で分解できると書かれています。今回の計測で時間のブレが大きかったのは、乱数を使う ECM (楕円曲線法) で分解していた可能性があります。ECMは乱数で曲線を選ぶ手法で、当たりの曲線を引くまでの時間に運が絡み、同じ数でも実行毎に時間が変わります。また、複数の手法を順に試すので、どの手法で見つかるかで時間が変わります。
ポラード・ロー法での反復回数は、小さい方の素因数 の平方根に比例します。 が約264 の No.25 では、約232 回かかります (反復回数に43億回指定で解けました！実行時間は約2分！😲)。ECMの計算量は、 の大きさに対して平方根より緩やかに増えます。そのため、約264 の素因数でも、数百ミリ秒で見つかる場合があります。No.25 のケースでは、この手法が効いたものだと思います (どの手法で分解したかは確認していません)。二次篩はさらに大きい数向けです。
「単一のアルゴリズムでは、得意な形と苦手な形がある」というのは、素因数分解の本質的な性質でもあります。ポラード・ロー法は、小さい素因数に強く、2つの大きい素因数に弱い手法です。MPUのような実用ライブラリが、複数の手法を組み合わせているのは、そのためだと思います。
【 反復回数の検証 】
ポラード・ロー法での素因数の分解における反復回数は、どの程度が適当なのかを実験してみました。2種類の桁数で測った結果、小さい方の素因数 に対して の2〜4倍にしておけば、ほぼ確実に分解できるようでした。4倍より増やしても、成功する割合は変わりません。ランダムな半素数で分解できた割合は以下のとおりです。
| 反復回数の上限 | 小さい素因数が 約234 (300件) |
小さい素因数が 約240 (150件) |
|---|---|---|
| の 0.5倍 | 29% | 31% |
| の 1倍 | 78% | 69% |
| の 2倍 | 99% | 97% |
| の 4倍 | 100% | 100% |
| の 8倍・16倍 | 100% | 100% |
【 プログラムコードの検証 】
(AI) Claude の協力で、独立した正解 (gmpy2) を使って以下のことを確認しました。積が元の数に一致すること、すべての因数が素数であること、昇順であること、未分解と返したものが本当に合成数であること。《徹底的に試験してくれますな》
| テスト項目 | 件数 | 誤り件数 |
| 特殊形（2128−1、2127+1、2127など） | 28 | 0 |
| 素数の冪（2〜5乗） | 671 | 0 |
| 中くらいの素数の積 × 大きい素数 | 1,500 | 0 |
| 純ランダム（2128 未満、反復上限は 220） | 3,000 | 0（うち93件は未分解） |
| 半素数（小さい方が 220〜238） | 200 | 0 |
| カーマイケル数 | 5 | 0 |
| サンプルデータテスト | 26 | 0 |
【 実行速度 】
この結果だけからは、64bit以下の数では、Inline::C (64bit版) が最も速いようです。MPUの factor も小さめな数では優勢ですが、64bit以上の数値では総じて Inline::C (128bit版) が優勢で、大差になっているケースも多いです。最後 No.25 のケースは、64bit素数×64bit素数ですが、 Inline::C (128bit版) の反復数の上限が小さく未分解で終わっています (既定の反復回数では解けません)。MPUの factor では分解しますが、計測時間が大きくブレる (160～2000 ms) ので参考値です。No.3～No.7 は、64bit版でも試している数 (15桁～20桁) ですので比較できます。(括弧内が64bit版の計測値)
※実行環境 (使用マシンの性能やOSの種類等) により実行時間は異なります。64bit以上の数における factor の平均的な 100ms コストは、実行環境による可能性があります。GMPありのLinux環境等では、本来の高速な結果になるかもしれません。
【 素因数分解の速度比較 】時間単位：ミリ秒 (ms)
| No. | 素因数分解する数 （下段：素因数） |
Inline::C (128bit) |
MPU factor |
| 1 | 1000 2, 2, 2, 5, 5, 5 |
0.009 ms (0.009 ms) |
0.006 ms |
| 2 | 1018081 1009, 1009 |
0.010 ms (0.016 ms) |
0.006 ms |
| 3 | 564281496252099 3, 3, 2639401, 23754611 |
0.022 ms (0.014 ms) |
0.016 ms |
| 4 | 9000000000000037 7, 223, 22739, 253552703 |
0.013 ms (0.007 ms) |
0.013 ms |
| 5 | 19473684210526317 3, 6491228070175439 |
0.011 ms (0.004 ms) |
0.010 ms |
| 6 | 999999866000004473 999999929, 999999937 |
0.109 ms (0.101 ms) |
0.123 ms |
| 7 | 18446744030759878681 4294967291, 4294967291 |
0.012 ms (0.005 ms) |
0.010 ms |
| 8 | 340282366920938463463374607431768211455 3, 5, 17, 257, 641, 65537, 274177, 6700417, 67280421310721 |
0.240 ms | 105.000 ms |
| 9 | 170141183460469231731687303715884105729 3, 56713727820156410577229101238628035243 |
0.025 ms |
102.607 ms |
| 10 | 170141183460469231731687303715884105728 2127 |
0.026 ms | 103.648 ms |
| 11 | 170141183460469231731687303715884105727 170141183460469231731687303715884105727 |
0.023 ms | 102.996 ms |
| 12 | 340282366920938463463374607431768211297 340282366920938463463374607431768211297 |
0.036 ms | 102.825 ms |
| 13 | 18446744073709551615 3, 5, 17, 257, 641, 65537, 6700417 |
0.014 ms |
0.016 ms |
| 14 | 147573952589676412927 193707721, 761838257287 |
0.251 ms | 103.472 ms |
| 15 | 147808829414345923316083210206383297601 380 |
0.021 ms | 102.162 ms |
| 16 | 74853500292876717928978827574247424 2100, 310 |
0.023 ms | 102.625 ms |
| 17 | 14282391973023475368482165658371922649 3779205203878650907, 3779205203878650907 |
0.022 ms | 104.303 ms |
| 18 | 73321016560309734457617669197591862911 4185456435071, 4185456435071, 4185456435071 |
0.023 ms | 103.010 ms |
| 19 | 1493998239615372206582229290524177943 17173943, 17173943, 17173943, 17173943, 17173943 |
0.023 ms | 104.429 ms |
| 20 | 48263651268012733650346617808893744017 10858187, 4444908829440194173331755827091 |
0.089 ms | 114.223 ms |
| 21 | 38240069382460403919538506226994224933 776620883, 49239043424564214196566365351 |
0.514 ms | 113.261 ms |
| 22 | 12300651168316729166140348424007860951 787513, 830080577, 18210378217, 1033311238903 |
12.345 ms | 105.428 ms |
| 23 | 821764905721378443954584372014114243 856441, 959511403262312808418308292123 |
0.056 ms | 102.422 m |
| 24 | 1296001987165015643369032371289 6000003067, 12000006133, 18000009199 |
4.805 ms | 108.574 ms |
| 25 | 274377929787655441099400654922386830051 15019950122537348993, 18267565973867844707 |
114970.653 ms (反復43億回指定) |
238.792 ms (計測不安定) |
検証データの説明：
11.はメルセンヌ素数 (2127-1)。12.は2128未満で最大の素数。13.は264-1。17.は素数の平方。18.は素数の3乗。19.は素数の5乗。20.は半素数 224 × 2102。21.は半素数 230 × 296。22.は4個の素数の積。23.は素数(100bit) × 素数(20bit)。24.はカーマイケル数。
【 プログラムソース 】
このプログラムのソースは以下のとおりです。「反復回数」を省略可能にするために、素数検索と同様に Inline Stack から直に引数を取り出し、返り値 (素因数の配列リファレンス) も Inline Stack にプッシュする方法にしています。公開関数はプログラム間のインターフェイスのみで、本来の関数の実体 (主エントリ) は factor128_core です。
# ======================================================================
# (Perl＋Inline::C) 128bit版
#
# 素因数分解: ポラード・ロー素因数分解法 (ブレント版)
# 素数判定等: (BPSW)ミラー・ラビン法＋強Lucasテスト＋モンゴメリ乗算
#
# 数の最大値: 340282366920938463463374607431768211455 (2^128-1,39桁)
#
# 関数の引数: SV* factor128(char *s, UV max_steps)
# 関数戻り値: 昇順の素因数"数字列"の配列リファレンスを返す
#
# ※max_steps: 合成数1個あたりの反復回数の上限 (0なら既定値 134217728)
# ※上限内に分解できなかった合成数は"C:数字列"という形でそのまま返す
#
# This Script is in the public domain, No rights reserved.
# Script written by N.O, Updated Last on 2026/10/03
# Claude (AI) assisted with the script optimization.
# ======================================================================
use strict;
use Inline (
C => 'DATA',
DIRECTORY => '_Inline',
);
# ---------------------------------------------
# テストサンプル (コマンドラインで数値指定)
# ---------------------------------------------
use Math::Prime::Util qw(factor);
use Time::HiRes qw(gettimeofday);
my ($n, $s, $xs) = @ARGV;
$n =~ s/[,_]*//go;
if ($n eq "" || $n eq "0" || $n =~ /\D/o || $s =~ /\D/o) {
print "\n[HELP] >perl $0 判定数 [反復回数 [MPU指定]] \n\n";
print "判定数は 340282366920938463463374607431768211455 以内の正の整数です。\n\n";
print "反復回数:'0'または省略で既定値。MPU指定:'0'以外ならMPU::factorを使用。\n";
exit;
}
$s ||= 0;
my(@primes, $start, $finish);
if ($xs) {
$start = gettimeofday();
@primes = factor($n);
$finish = gettimeofday();
print "\n*** 素因数分解(MPU::factor)の結果 ***\n\n";
} else {
$start = gettimeofday();
my $ref = factor128($n, $s);
$finish = gettimeofday();
@primes = @$ref if ref $ref;
print "\n*** 素因数分解 (128bit版) の結果 ***\n\n";
}
my $elapse = int(($finish - $start) * 1_000_000) / 1000;
print "判定数: $n\n";
print "素因数: ", join(", ", @primes), "\n\n";
print "実行時間: ", sprintf("%.3f", $elapse), " ミリ秒(ms)\n";
exit;
# -----------------------------------
# ここからInline::C(言語)での記述
# -----------------------------------
__DATA__
__C__
/* ==================================================== *
* 128bit版 素因数分解: ポラード・ロー法 (ブレント版) *
* ==================================================== */
#include <stdint.h>
#include <math.h>
typedef uint64_t u64;
typedef unsigned __int128 u128;
/* 数字列 <-> 128bit整数の変換。成功: 0, 失敗: -1 */
static int parse_u128(const char *s, u128 *out) {
const u128 MAXV = ~(u128)0;
u128 v = 0;
int nd = 0;
while (*s == ' ' || *s == '\t') s++;
for (; *s >= '0' && *s <= '9'; s++) {
unsigned dg = (unsigned)(*s - '0');
if (v > (MAXV - dg) / 10) return -1; /* v*10+dg が 2^128 を超える */
v = v * 10 + dg;
nd++;
}
while (*s == ' ' || *s == '\t' || *s == '\n' || *s == '\r') s++;
if (nd == 0 || *s != '\0') return -1;
*out = v;
return 0;
}
/* buf は 40 バイト以上 (39桁 + NUL) */
static void u128_to_str(u128 v, char *buf) {
char tmp[40];
int i = 0, j = 0;
if (v == 0) {
buf[0] = '0'; buf[1] = '\0'; return;
}
while (v) {
tmp[i++] = (char)('0' + (int)(v % 10)); v /= 10;
}
while (i) buf[j++] = tmp[--i];
buf[j] = '\0';
}
/* 128bit 補助関数 */
static int ctz128(u128 x) { /* x != 0 */
u64 lo = (u64)x;
return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll((u64)(x >> 64));
}
static int bitlen128(u128 x) { /* x != 0 */
u64 hi = (u64)(x >> 64);
return hi ? 128 - __builtin_clzll(hi) : 64 - __builtin_clzll((u64)x);
}
static u128 isqrt128(u128 n) { /* floor(sqrt(n)) */
if (n < 2) return n;
u128 x = (u128)1 << ((bitlen128(n) + 1) / 2); /* sqrt(n) 以上の初期値 */
for (;;) {
u128 y = (x + n / x) >> 1;
if (y >= x) return x;
x = y;
}
}
/* ヤコビ記号 (a/n) n は奇数の正の数。戻り値 -1, 0, 1 */
static int jacobi128(u128 a, u128 n) {
int result = 1;
a %= n;
while (a) {
while ((a & 1) == 0) {
a >>= 1;
unsigned r = (unsigned)(n & 7);
if (r == 3 || r == 5) result = -result;
}
u128 t = a; a = n; n = t;
if ((a & 3) == 3 && (n & 3) == 3) result = -result;
a %= n;
}
return n == 1 ? result : 0;
}
/* 小さな符号付き整数 x の n による剰余 */
static u128 small_mod(long x, u128 n) {
u128 m = (u128)(x < 0 ? -(long long)x : x) % n;
return x < 0 ? (m ? n - m : 0) : m;
}
/* 64bit: モンゴメリ乗算 + ミラー・ラビン(n < 2^64 で決定的) */
static u64 mont_inv64(u64 n) {
u64 x = n;
for (int i = 0; i < 5; i++) x *= 2 - n * x;
return x;
}
static u64 mont_mul64(u64 a, u64 b, u64 n, u64 ninv) {
u128 t = (u128)a * b;
u64 q = (u64)t * ninv;
u64 hi = (u64)(t >> 64);
u64 mh = (u64)(((u128)q * n) >> 64);
return hi >= mh ? hi - mh : hi - mh + n;
}
static const u64 SMALL_PRIMES[] = {2,3,5,7,11,13,17,19,23,29,31,37};
static int is_prime_u64(u64 n) {
if (n < 2) return 0;
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return n == SMALL_PRIMES[i];
}
if (n < 37ULL * 37) return 1;
static const u64 b2[] = {31, 73};
static const u64 b3[] = {2, 7, 61};
static const u64 b4[] = {2, 13, 23, 1662803};
static const u64 b7[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
const u64 *bases; int nb;
if (n < 9080191ULL) { bases = b2; nb = 2; }
else if (n < 4759123141ULL) { bases = b3; nb = 3; }
else if (n < 1122004669633ULL) { bases = b4; nb = 4; }
else { bases = b7; nb = 7; }
u64 ninv = mont_inv64(n);
u64 one = (0 - n) % n;
u64 mone = n - one;
u64 r2 = (u64)(((u128)one * one) % n);
u64 d = n - 1;
int s = __builtin_ctzll(d);
d >>= s;
for (int i = 0; i < nb; i++) {
u64 a = bases[i] % n;
if (a == 0) continue;
u64 b = mont_mul64(a, r2, n, ninv);
u64 x = one;
for (u64 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul64(x, b, n, ninv);
b = mont_mul64(b, b, n, ninv);
}
if (x == one || x == mone) continue;
int composite = 1;
for (int r = 1; r < s; r++) {
x = mont_mul64(x, x, n, ninv);
if (x == mone) { composite = 0; break; }
}
if (composite) return 0;
}
return 1;
}
/* 128bit: モンゴメリ乗算 (R = 2^128, n は奇数) */
/* 128bit × 128bit → 256bit (hi:lo) */
static void mul_128x128(u128 a, u128 b, u128 *hi, u128 *lo) {
u64 a0 = (u64)a, a1 = (u64)(a >> 64), b0 = (u64)b, b1 = (u64)(b >> 64);
u128 p00 = (u128)a0 * b0, p01 = (u128)a0 * b1;
u128 p10 = (u128)a1 * b0, p11 = (u128)a1 * b1;
u128 mid = (p00 >> 64) + (u64)p01 + (u64)p10; /* 3*2^64 未満なので溢れない */
*lo = ((u128)(u64)mid << 64) | (u64)p00;
*hi = p11 + (p01 >> 64) + (p10 >> 64) + (mid >> 64);
}
/* n^-1 mod 2^128 (ニュートン法: 有効bit 3→6→12→24→48→96→192) */
static u128 mont_inv128(u128 n) {
u128 x = n;
for (int i = 0; i < 6; i++) x *= 2 - n * x;
return x;
}
/* a*b*R^-1 mod n (a, b < n) */
static u128 mont_mul128(u128 a, u128 b, u128 n, u128 ninv) {
u128 hi, lo, mh, dummy;
mul_128x128(a, b, &hi, &lo);
u128 q = lo * ninv; /* mod 2^128 で巻き戻る */
mul_128x128(q, n, &mh, &dummy); /* q*n の上位128bit */
return hi >= mh ? hi - mh : hi - mh + n;
}
static u128 add_mod128(u128 a, u128 b, u128 n) {
return a >= n - b ? a - (n - b) : a + b;
}
static u128 sub_mod128(u128 a, u128 b, u128 n) {
return a >= b ? a - b : n - (b - a);
}
/* x / 2 mod n (n は奇数) */
static u128 half_mod128(u128 x, u128 n) {
return (x & 1) == 0 ? x >> 1 : (x >> 1) + (n >> 1) + 1;
}
/* 強確率素数テスト: n を基数 a で判定。ninv/one/r2 はモンゴメリ準備値 */
static int mr128(u128 n, u128 a, u128 ninv, u128 one, u128 r2) {
u128 nm1 = n - 1;
int s = ctz128(nm1);
u128 d = nm1 >> s;
a %= n;
if (a == 0) return 1;
u128 mone = n - one;
u128 b = mont_mul128(a, r2, n, ninv);
u128 x = one;
for (u128 e = d; e; e >>= 1) {
if (e & 1) x = mont_mul128(x, b, n, ninv);
b = mont_mul128(b, b, n, ninv);
}
if (x == one || x == mone) return 1;
for (int r = 1; r < s; r++) {
x = mont_mul128(x, x, n, ninv);
if (x == mone) return 1;
}
return 0;
}
/* 強Lucasテスト (Selfridge の方法A: D = 5, -7, 9, -11, ... P = 1, Q = (1-D)/4) */
/* n: 奇数。n+1 が溢れないこと (n != 2^128-1) は呼び出し側の試し割りで保証 */
static int strong_lucas128(u128 n, u128 ninv, u128 one, u128 r2) {
u128 rt = isqrt128(n);
if (rt * rt == n) return 0; /* 平方数は D が見つからない */
long D = 5;
for (int guard = 0; guard < 100000; guard++) {
u128 dm = small_mod(D, n);
int j = jacobi128(dm, n);
if (j == -1) break;
if (j == 0) return n == (u128)(D < 0 ? -D : D); /* 共通因数あり (nが|D|と等しい場合だけ素数) */
D = D > 0 ? -(D + 2) : -(D - 2);
}
long Q = (1 - D) / 4;
u128 Dm = mont_mul128(small_mod(D, n), r2, n, ninv);
u128 Qm = mont_mul128(small_mod(Q, n), r2, n, ninv);
u128 np1 = n + 1;
int s = ctz128(np1);
u128 d = np1 >> s;
u128 U = one, V = one, Qk = Qm; /* k = 1: U=1, V=P=1, Q^1 */
for (int i = bitlen128(d) - 2; i >= 0; i--) {
U = mont_mul128(U, V, n, ninv); /* U_2k = U_k V_k */
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n); /* V_2k = V_k^2 - 2Q^k */
Qk = mont_mul128(Qk, Qk, n, ninv);
if ((d >> i) & 1) { /* k -> k+1 */
u128 U2 = half_mod128(add_mod128(U, V, n), n);
u128 V2 = half_mod128(add_mod128(mont_mul128(Dm, U, n, ninv), V, n), n);
U = U2; V = V2;
Qk = mont_mul128(Qk, Qm, n, ninv);
}
}
if (U == 0 || V == 0) return 1;
for (int r = 1; r < s; r++) {
V = sub_mod128(mont_mul128(V, V, n, ninv), add_mod128(Qk, Qk, n), n);
Qk = mont_mul128(Qk, Qk, n, ninv);
if (V == 0) return 1;
}
return 0;
}
/* 素数判定 (128bit) */
static int is_prime_u128(u128 n) {
if ((n >> 64) == 0) return is_prime_u64((u64)n); /* 64bit以下は確定判定 */
for (int i = 0; i < 12; i++) {
if (n % SMALL_PRIMES[i] == 0) return 0; /* n > 2^64 なので n 自身が素数ではない */
}
u128 ninv = mont_inv128(n);
u128 one = (0 - n) % n; /* R mod n */
u128 r2 = one; /* R^2 mod n : R を128回倍加 (256bit除算を避ける) */
for (int i = 0; i < 128; i++) r2 = add_mod128(r2, r2, n);
if (!mr128(n, 2, ninv, one, r2)) return 0;
return strong_lucas128(n, ninv, one, r2);
}
/* GCD (バイナリ法) / 完全平方判定 */
static u64 gcd_u64(u64 a, u64 b) {
if (a == 0) return b;
if (b == 0) return a;
int sh = __builtin_ctzll(a | b);
a >>= __builtin_ctzll(a);
do {
b >>= __builtin_ctzll(b);
if (a > b) {
u64 t = a; a = b; b = t;
}
b -= a;
} while (b);
return a << sh;
}
static u128 gcd_u128(u128 a, u128 b) {
if (a == 0) return b;
if (b == 0) return a;
int sh = ctz128(a | b);
a >>= ctz128(a);
do {
b >>= ctz128(b);
if (a > b) {
u128 t = a; a = b; b = t;
}
b -= a;
} while (b);
return a << sh;
}
/* n が完全平方数なら平方根、そうでなければ 0 */
static u64 exact_sqrt64(u64 n) {
u64 r = (u64)sqrt((double)n);
if (r > 4294967295ULL) r = 4294967295ULL;
while ((u128)r * r > n) r--;
while ((u128)(r + 1) * (r + 1) <= n) r++;
return (u128)r * r == n ? r : 0;
}
/* r < 2^64 なので r*r は溢れない */
static u128 exact_sqrt128(u128 n) {
u128 r = isqrt128(n);
return r * r == n ? r : 0;
}
/* r^k == n か (途中で n を超えるか128bitを超えるなら 0 を返す) */
static int pow_equals(u128 r, int k, u128 n) {
u128 v = 1;
for (int i = 0; i < k; i++) {
if (__builtin_mul_overflow(v, r, &v)) return 0;
if (v > n) return 0;
}
return v == n;
}
/* n が k乗数 (k >= 3) なら k乗根、そうでなければ 0 (浮動小数で見当をつけて前後を厳密に確かめる) */
static u128 exact_root128(u128 n, int k) {
u128 r = (u128)(pow((double)n, 1.0 / k) + 0.5);
for (u128 c = (r > 2 ? r - 1 : 2); c <= r + 1; c++) {
if (pow_equals(c, k, n)) return c;
}
return 0;
}
/* ポラード・ロー法 (ブレント版・GCDまとめ取り) */
#define BATCH 128
/* 64bit: n は奇数の合成数 (平方数・小さい素因数は排除済み)。必ず真の因数を返す */
static u64 rho_brent64(u64 n) {
u64 ninv = mont_inv64(n);
for (u64 c = 1; ; c++) {
u64 nc = n - c;
u64 x = 2, y = 2, ys = 2, q = 1, g = 1;
for (u64 r = 1; g == 1; r <<= 1) {
x = y;
for (u64 i = 0; i < r; i++) {
y = mont_mul64(y, y, n, ninv); y = y >= nc ? y - nc : y + c;
}
for (u64 k = 0; k < r && g == 1; k += BATCH) {
ys = y;
u64 steps = (r - k < BATCH) ? r - k : BATCH;
for (u64 i = 0; i < steps; i++) {
y = mont_mul64(y, y, n, ninv); y = y >= nc ? y - nc : y + c;
q = mont_mul64(q, x > y ? x - y : y - x, n, ninv);
}
g = gcd_u64(q, n);
}
}
if (g == n) {
do {
ys = mont_mul64(ys, ys, n, ninv); ys = ys >= nc ? ys - nc : ys + c;
g = gcd_u64(x > ys ? x - ys : ys - x, n);
} while (g == 1);
}
if (g != n) return g;
}
}
/* 128bit: n は奇数の合成数。反復回数が max_steps を超えたら 0 を返す (未分解) */
static u128 rho_brent128(u128 n, u64 max_steps) {
u128 ninv = mont_inv128(n);
u64 total = 0;
/* 失敗したら c を変えてリトライ */
for (u64 c = 1; ; c++) {
u128 nc = n - c;
u128 x = 2, y = 2, ys = 2, q = 1, g = 1;
for (u64 r = 1; g == 1; r <<= 1) {
x = y;
for (u64 i = 0; i < r; i++) {
y = mont_mul128(y, y, n, ninv); y = y >= nc ? y - nc : y + c;
}
total += r;
for (u64 k = 0; k < r && g == 1; k += BATCH) {
ys = y;
u64 steps = (r - k < BATCH) ? r - k : BATCH;
for (u64 i = 0; i < steps; i++) {
y = mont_mul128(y, y, n, ninv); y = y >= nc ? y - nc : y + c;
q = mont_mul128(q, x > y ? x - y : y - x, n, ninv);
}
total += steps;
g = gcd_u128(q, n);
}
if (g == 1 && total > max_steps) return 0;
}
/* まとめ過ぎたらやり直す */
if (g == n) {
do {
ys = mont_mul128(ys, ys, n, ninv); ys = ys >= nc ? ys - nc : ys + c;
g = gcd_u128(x > ys ? x - ys : ys - x, n);
} while (g == 1);
}
if (g != n) return g;
}
}
/* 素因数分解 */
#define MAXF 130
static void push_f(u128 *out, char *fl, int *cnt, u128 v, int composite) {
out[*cnt] = v; fl[*cnt] = (char)composite; (*cnt)++;
}
/* n > 1 で、素因数はすべて 1000 より大きい */
static void factorize_rec(u128 n, u128 *out, char *fl, int *cnt, u64 max_steps) {
if (n == 1) return;
if ((n >> 64) == 0) { /* 64bit以下: 必ず完全に分解できる */
u64 m = (u64)n;
if (is_prime_u64(m)) {
push_f(out, fl, cnt, n, 0);
return;
}
u64 r = exact_sqrt64(m);
if (r) {
factorize_rec(r, out, fl, cnt, max_steps);
factorize_rec(r, out, fl, cnt, max_steps);
return;
}
u64 f = rho_brent64(m);
factorize_rec(f, out, fl, cnt, max_steps);
factorize_rec(m / f, out, fl, cnt, max_steps);
return;
}
if (is_prime_u128(n)) {
push_f(out, fl, cnt, n, 0);
return;
}
u128 r = exact_sqrt128(n);
if (r) {
factorize_rec(r, out, fl, cnt, max_steps);
factorize_rec(r, out, fl, cnt, max_steps);
return;
}
/* 完全冪の検出: 3乗・5乗・7乗・11乗 (素因数は全て1009以上なので 2^128 未満なら12乗以上はない) */
{
static const int PK[] = {3, 5, 7, 11};
for (int i = 0; i < 4; i++) {
u128 rt = exact_root128(n, PK[i]);
if (rt) {
for (int j = 0; j < PK[i]; j++) factorize_rec(rt, out, fl, cnt, max_steps);
return;
}
}
}
u128 f = rho_brent128(n, max_steps);
/* 上限内に分解できなかった合成数 */
if (f == 0) {
push_f(out, fl, cnt, n, 1);
return;
}
factorize_rec(f, out, fl, cnt, max_steps);
factorize_rec(n / f, out, fl, cnt, max_steps);
}
/* 公開関数の実体 */
static SV* factor128_core(const char *s, UV max_steps_in) {
u128 n;
if (parse_u128(s, &n) != 0)
croak("factor128: 不正な数字列です (数字のみ最大 340282366920938463463374607431768211455)");
u64 max_steps = max_steps_in ? (u64)max_steps_in : 134217728ULL;
u128 f[MAXF];
char fl[MAXF];
int cnt = 0;
/* 1000未満の素因数を先に除く */
if (n >= 2) {
for (u64 d = 2; d < 1000; d += (d == 2 ? 1 : 2)) {
while (n % d == 0) {
push_f(f, fl, &cnt, d, 0); n /= d;
}
}
/* 素因数が全て 1009 以上 → n は素数 */
if (n > 1) {
if (n < 1009ULL * 1009) {
push_f(f, fl, &cnt, n, 0);
} else {
factorize_rec(n, f, fl, &cnt, max_steps);
}
}
}
/* 昇順に整列 (挿入ソート、フラグも連動) */
for (int i = 1; i < cnt; i++) {
u128 v = f[i]; char c = fl[i]; int j = i - 1;
while (j >= 0 && f[j] > v) {
f[j + 1] = f[j]; fl[j + 1] = fl[j]; j--;
}
f[j + 1] = v; fl[j + 1] = c;
}
AV *av = newAV();
for (int i = 0; i < cnt; i++) {
char buf[44];
char *p = buf;
if (fl[i]) {
*p++ = 'C'; *p++ = ':';
}
u128_to_str(f[i], p);
av_push(av, newSVpv(buf, 0));
}
return newRV_noinc((SV*)av);
}
/* Perlへの公開関数 (反復回数を省略可能) */
void factor128(...) {
Inline_Stack_Vars;
const char *s;
UV max_steps = 0;
if (Inline_Stack_Items < 1)
croak("Usage: factor128(n [, max_steps])");
s = SvPV_nolen(Inline_Stack_Item(0));
if (Inline_Stack_Items > 1 && SvOK(Inline_Stack_Item(1))) max_steps = SvUV(Inline_Stack_Item(1));
SV *rv = factor128_core(s, max_steps);
Inline_Stack_Reset;
Inline_Stack_Push(sv_2mortal(rv));
Inline_Stack_Done;
}
