128ビット符号付き整数の最大値は素数 - Rustで任意精度整数演算
概要
rug を利用して実装します。
実行環境
- CPU: Intel Core i7 1.8GHz
- メモリ: 16GB
- OS(ホスト): Windows 10 Home 21H1
- WSL2: Ubuntu 20.04.3
- rustc: Ver. 1.55.0
- cargo: Ver. 1.55.0
符号付き整数型の範囲について
Rustには組み込みの整数型として
各々の型で表現可能な最大値を実際に見てみましょう。型の定数として値は（128ビットの場合） i128::MAXから取ることができます。
一応コード
fn main() {
println!("{} {}", i8::BITS, i8::MAX);
println!("{} {}", i16::BITS, i16::MAX);
println!("{} {}", i32::BITS, i32::MAX);
println!("{} {}", i64::BITS, i64::MAX);
println!("{} {}", i128::BITS, i128::MAX);
}
| 桁数(十進) | ||
|---|---|---|
表には十進法表記を示しました。語呂合わせでも作ってちょっと覚えてみてもいいかもしれない、と思わせるような数字の並びですが、二進法と十六進法なら自明です：
ただし、
表には素因数分解も示しました。
つまり、
確かめてみましょう。
fn main() {
let nums = vec![
i8::MAX as i128,
i16::MAX as i128,
i32::MAX as i128,
i64::MAX as i128,
i128::MAX,
];
for n in nums {
for m in 2.. {
let (q, r) = (n / m, n % m);
if r == 0 {
println!("{} is composite", n);
break;
}
if q <= m {
println!("{} is prime", n);
break;
}
}
}
}
mersenne_primes/main.rs at main · roiban1344/mersenne_primes
出力（中断）：
127 is prime
32767 is composite
2147483647 is prime
9223372036854775807 is composite
…………
終わらない！ 当然です。170141183460469231731687303715884105727 is prime が出力されるまでには long longのループを回した苦い経験が蘇ります。
ところがこの
メルセンヌ数とは
歴史的背景について。
非負整数
「メルセンヌ数」の呼び名は、16世紀フランスのカトリックの司祭・数学者メルセンヌ[3]に由来します。彼は
の11個[4] に限られると予想しました。後の時代になって彼のリストには
- 3つの抜け：
n=61, 89, 107 M_n - 2つの誤り：
n=67, 257 M_n
があることが示されましたが、彼に敬意を表して今日でもその名が残っています。なお、
を得る方法は後で見ることになります。
さて、
という分解[5]から、
しかし逆は成り立ちません。すなわち、素数
64ビット整数の範囲の
const PRIMES: [i32; 18] = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
];
fn main() {
for p in PRIMES {
let m: u64 = (1 << p) - 1;
{
let mut m = m;
let mut factors = vec![];
for d in 2.. {
while m % d == 0 {
m /= d;
factors.push(d);
}
if m <= d * d {
if m != 1 {
factors.push(m);
}
break;
}
}
factors.sort();
println!(
"{:>2} {} {:?}",
p,
if factors.len() == 1 { "p" } else { "c" },
factors
);
}
}
}
出力：
2 p [3]
3 p [7]
5 p [31]
7 p [127]
11 c [23, 89]
13 p [8191]
17 p [131071]
19 p [524287]
23 c [47, 178481]
29 c [233, 1103, 2089]
31 p [2147483647]
37 c [223, 616318177]
41 c [13367, 164511353]
43 c [431, 9719, 2099863]
47 c [2351, 4513, 13264529]
53 c [6361, 69431, 20394401]
59 c [179951, 3203431780337]
61 p [2305843009213693951]
出力の各行は
- 1列目:
p - 2列目:
PならM_p Cなら合成数 - 3列目: 素因数たち（昇順）
を意味します。この範囲で既に
メルセンヌ数の素因数と、素数が無限に存在することについて
特に
特筆すべきことに、この事実は素数が無限に存在することの証明にもなります。最大の素数
参考：M.アイグナー，G.M.ツィーグラー， 訳 蟹江幸博『天書の証明』，丸善出版，2002, 第1章.
メルセンヌ数と完全数
メルセンヌ素数は完全数[6] との密接な繋がりから関心が向けられてきました。
完全数とは、自身を除く約数の和が自身と一致する非負整数のことです。小さいほうから
と、確かに
証明 偶数の完全数とメルセンヌ素数が一対一対応すること
M_p n=2^{p-1}M_p
正整数
したがって
n M_p n=2^{p-1}M_p
一方、
左辺は整数だから、右辺第2項
よって
オイラーに帰せられる後者の証明は非常に面白いです。
参考：高木貞治『初頭整数論講義 第2版』，共立出版，1971．
この事実により、古代ギリシャから知られていた偶完全数の探索はメルセンヌ数探索と表裏一体になりました。一方、奇数の完全数は未だ一つも知られていないうえ存在しないとも証明されておらず、数論における有名な未解決問題になっています。
リュカテスト
本題に入ります。1876年、リュカ[7]は以下の判定法を利用し、
(リュカテスト)
確かに
一方
確かに
リュカ・テストが正しいことの証明は、その核心でフィボナッチ数（特に整除性と周期性）が関わり、初等的でありながら全く自明でない結果が導かれる、非常に興味深いものです。
証明 リュカテスト
証明の大筋は参考文献[中村]p.82によるものです。
補題(1)
フィボナッチ数列を
リュカ数列を
によってそれぞれ再帰的に定義する。なお、
とおくと、
証明：帰納法による。
補題(2)
証明：帰納法。
補題(3)
証明：
ほぼ同様にして
組み合わせると、
を得る。
補題(4)
正整数
証明:
まず、定義から
が成り立つ。この関係式を繰り返し適用して、
一般に、
であることから、
により、
準備ができたので以下本証明に入ります。
リュカテストが正しいことの証明
(A) M_p R_{p-1} \equiv 0 \mod M_p
補題(2), (3)から、
一般に奇素数
のように振る舞うため
(B) R_{p-1} \equiv 0 \mod M_p M_p
このとき
に
補題(3)から
確率的な素数判定法が実用に耐えるのは、判定対象の整数の 対数に関する 多項式時間のアルゴリズムであるからです。一方、例えば試し割りによる方法だと、判定対象の整数それ自体に関する多項式時間です。
計算量の観点から見たリュカテストの著しい性質は、決定的であるにもかかわらず確率的な素数判定法と同等かそれ以上に高速であることです。
この非常に効率的な判定法を得たリュカは、1876年に
とはいえ実際に手計算してみようとすると分かりますが、
任意精度演算用クレート rug
リュカテストを実装します。i128 で扱うと2乗を取る過程でオーバーフローが起こります。
そこで任意桁数の整数演算を導入します。これを実現するRustのクレートとそれぞれの特徴を以下に挙げます[10]
- num-bigint: 純Rust製。使いやすい。
- rug: GMPライブラリ等のバンドル。使いやすい。高速。
- rust-gmp: 同上（のようだがあまりドキュメントが充実しておらず詳細不明）。
- ramp: Nightlyビルド限定。
- apint: 純Rust製、未完。
今回はこの中から rug を採用することにしました。
num-bigintについて
num-bigintは直感的に扱えて使いやすいですが、どうもバグらしき挙動を示します。panic! します。
よりによって
もっとも性能ではrugが上回っているため、num-bigintを使う必要はあまりなさそうです。
rugはGNUに属す3つの任意精度演算用ライブラリ(GMP, MPFR, MPC)のラッパーです。任意精度整数はこの中のGMPによって提供されます。Mathematicaの内部でも使われるなど[11]、歴史と信頼性のあるライブラリです。任意精度整数演算を行うなら最上の選択肢の一つであることは間違いないでしょう。
rug は低レベルのFFIを提供するクレートである gmp-mpfr-sys に依存するため、ビルドにはCargo.tomlにdependenciesを追加する以外にも操作が必要になることがあります。詳細は gmp-mpfr-sys のドキュメントを参照してください[12]
rug の任意精度整数演算は Integer 構造体が担っています。
-
Copyトレイトを備えていないこと。 -
Incomplete-computationという機構があること。参照同士の演算の結果は直ちに反映されず、
fromやassignを介する必要があります。
に注意すれば利用は非常に簡単です。基本的な算術演算・ビット操作は演算子オーバーロードされており、プリミティブ型との演算(+ 12 や == 60)も可能です。
リュカテストの実装
use rug::ops::Pow;
use rug::Integer;
//Returns a list of prime numbers less than n.
fn primes(n: u32) -> Vec<u32> {
let mut is_prime = vec![true; n as usize];
let mut primes = vec![];
for i in 2..n {
if is_prime[i as usize] {
primes.push(i);
for j in 2.. {
let k = i * j;
if !(k < n) {
break;
}
is_prime[k as usize] = false;
}
}
}
primes
}
fn mersenne_number(n: u32) -> Integer {
(Integer::from(1) << n) - 1
}
//Executes the Lucas test for p-th Mersenne number.
//p must be a prime number of the form 4n + 3.
fn lucas_test(p: u32) -> bool {
let m = mersenne_number(p);
let mut s = Integer::from(3);
for _ in 2..=p - 1 {
s = s.pow(2) - 2;
while s >= m {
s = Integer::from(&s >> p) + (s & &m);
if s == m {
s = Integer::from(0);
break;
}
}
}
s == 0
}
fn main() {
let primes = primes(10000);
for p in primes {
if p % 4 == 3 {
if lucas_test(p) {
println!("{}", p);
}
}
}
}
primesはただの篩です。lucas_testがリュカテストの実装です。
初期化や演算にInteger::fromを介していることを除けばプリミティブ型で書くのとそう変わりません。rugの使い方とは別の面で実装上の注意が2つだけあります。
まず、メルセンヌ数に関する剰余演算は汎用の剰余演算子である % より高速に実行できます。
から、
これは二進法版の九去法のようなものです。実装上は
while s >= m {
s = Integer::from(&s >> p) + (s & &m);
if s == m {
s = Integer::from(0);
break;
}
}
もう一つ、減算を行う箇所ではコード上 s が負値をとりえます:
s = s.pow(2) - 2;
ところが、上記リュカテストの証明中で述べた理由から、右辺の s はループ中 0 はとりません。また、
さて、実行すると
3
7
19
31
107
127
607
1279
2203
4423
たった10個！ この区間には
リュカ・レーマーテスト
リュカテストの明らかな弱点は
(リュカ・レーマーテスト)
初項を
リュカ・レーマーテストもリュカテストと同様、フィボナッチ・リュカ数列に対応する2階線形回帰数列を考えることで証明が可能ですが、より簡潔な証明がBruceによって与えられています[13]。Wikipediaに掲載されているものがこれに基づくため、詳細はそちらに譲ります: Proof of correctness
先のリュカテストに数行の変更を加えて、リュカ・レーマーテストの実装は以下のようになります。
use rug::ops::Pow;
use rug::Integer;
//Returns a list of prime numbers less than n.
fn primes(n: u32) -> Vec<u32> {
let mut is_prime = vec![true; n as usize];
let mut primes = vec![];
for i in 2..n {
if is_prime[i as usize] {
primes.push(i);
for j in 2.. {
let k = i * j;
if !(k < n) {
break;
}
is_prime[k as usize] = false;
}
}
}
primes
}
fn mersenne_number(n: u32) -> Integer {
(Integer::from(1) << n) - 1
}
-//Executes the Lucas test for p-th Mersenne number.
-//p must be a prime number of the form 4n + 3.
-fn lucas_test(p: u32) -> bool {
+//Executes the Lucas-Lehmer test for p-th Mersenne number.
+//p must be an odd prime number.
+fn lucas_lehmer_test(p: u32) -> bool {
let m = mersenne_number(p);
- let mut s = Integer::from(3);
+ let mut s = Integer::from(4);
for _ in 2..=p - 1 {
s = s.pow(2) - 2;
while s >= m {
s = Integer::from(&s >> p) + (s & &m);
if s == m {
s = Integer::from(0);
break;
}
}
}
s == 0
}
fn main() {
let primes = primes(10000);
for p in primes {
- if p % 4 == 3 {
- if lucas_test(p) {
+ if p != 2 {
+ if lucas_lehmer_test(p) {
println!("{}", p);
}
}
}
}
判定対象の上限を30000にまで上げて実行：
3
5
7
13
17
19
31
61
89
107
127
521
607
1279
2203
2281
3217
4253
4423
9689
9941
11213
19937
21701
23209
ものの数分の実行で26番目のメルセンヌ素数
メルセンヌ素数探索の歴史が表にまとめられています。いま見つかった中で最大の
ちなみに24番目の
現在の世界記録は2018年に発見された
51個のうち半分までの発見は家庭用PCで瞬時に再現できる一方、残り半分は今日までの人類のメルセンヌ素数探索の歴史全てがかけられてようやく見つかるのです。
GIMPS
既知最大のメルセンヌ素数
このプロジェクトには誰でもアカウントを作って貢献することができます。詳しくは公式によるGIMPS - Getting Started、またGIMPSの始め方 - Qiitaを参照してください。
GIMPSの探索プログラムには前処理や種々の最適化、検算が含まれますが、核心部分ではやはりリュカ・レーマーテストが利用されています[16]。先述の確率的素数判定法や素因数分解アルゴリズムがコンピューター時代以後に洗練されてきたことを考えると、1世紀半に渡って最前線で利用されているリュカ・レーマーテストの強力さは驚嘆に値します。
素因数分解
リュカ・レーマーテストは素因数に関する情報を与えてくれません。素因数分解と素数判定が質的に異なる問題であることを示す実例の一つです。
メルセンヌが誤って素数だと予想した2つの合成数のうち、
一方
数学計算ソフトSageMathはオンライン版SageMathCellが利用可能です。ここに単に
factor(2^257-1)
と打ち込んで"Evaluate"を押すと、10秒程度で素因数分解が得られます。
今日の素因数分解アルゴリズムとしては一般数体篩法が最高効率です。類似のアイデアに基づくもう少し簡単なアルゴリズムには二次篩法もあります。これらの実装によるメルセンヌ数の素因数分解は今後の課題です。
まとめ
メルセンヌ数の歴史と基本的な性質を概観し、その素数判定法であるリュカ・レーマーテストをRustで実装しました。アルゴリズム的にはごく単純であるため、特に最適化を目指さなかったこの範囲ではプログラミング的には少々味気なくはありますが、それもrugとGMPライブラリのおかげです。単純とはいえ、任意精度整数演算が本質的に役に立つ例としても良い題材であるように思います。
128ビット符号付き整数の最大値にまた出会う機会があれば、その背景に拡がる数論の世界に思いを馳せましょう。
参考
書籍
- [中村]中村滋『［改訂版］フィボナッチ数の小宇宙／フィボナッチ数・リュカ数・黄金分割』，日本評論社，2008. 特に第10章『リュカ・テスト』．
- [Dickson]Dickson, Leonard, E, History of the Theory of Numbers Volume 1: DIvisibility and Primality, Dover, 2005. (原版:Carnegie Institution of Washington, 1919. https://archive.org/details/historyoftheoryo01dick) 特にChap. I "Perfect, Multipy Perfect, and Amicable Numbers." 及び Chap. XVII "Recurring Series ; Lucas'
u_n,\,v_n
本記事の内容は[中村]に大きく依っています。[Dickson]はまだ
ウェブサイト
メルセンヌ数に関する最新のデータや歴史はこれらが詳しいです。本記事の数学的な内容や歴史はほとんどここにも載っています。
-
アーキテクチャ依存の
usize,isizeはここでは除外 ↩︎ -
Marin Mersenne(1588-1648) https://en.wikipedia.org/wiki/Marin_Mersenne ↩︎
-
A109461 Mersenne's original list of "Mersenne" exponents. ↩︎
-
2^{ab}-1 2^a-1 2^b-1 2^{2\cdot 4}-1=255 2^2-1=3 2^4-1=15 3\cdot 15=45 2^n-1 -
Édouard Luca(1842-1891) https://ja.wikipedia.org/wiki/エドゥアール・リュカ ↩︎
-
M_{127} M_{61}, M_{89}, M_{107} -
num-bigintの"Alternatives"の項の表より。また、やや古いがRedditの投稿。 ↩︎
-
筆者はWindowsのホスト上での構築に挫折し、WSL2に移行しました。Rust公式のDockerイメージ(
rust:latest)を使っても特別な操作が不要であることを確認しています。 ↩︎ -
Bruce, James W. "A really trivial proof of the Lucas-Lehmer primality test." The American Mathematical Monthly 100.4 (1993): 370-371. https://research.edgehill.ac.uk/en/publications/a-really-trivial-proof-of-the-lucas-lehmer-primality-test-2
Fermat's Library | A really trivial proof of the Lucas-Lehmer test annotated/explained version.、 ↩︎ -
Tuckerman, Bryant. "The 24th Mersenne prime." Proceedings of the National Academy of Sciences 68.10 (1971): 2319-2320. https://doi.org/10.1073/pnas.68.10.2319 ↩︎
-
GIMPSによるプレスリリース：Mersenne Prime Discovery - 2^82589933-1 is Prime! ↩︎
-
2018年からはフェルマーテストも採用されているようです。Is there a better method than Lucas-Lehmer? ↩︎
-
Frank Nelson Cole(1861-1926) Frank Nelson Cole ↩︎
Discussion