| title | emoji | type | topics | published | ||
|---|---|---|---|---|---|---|
128ビット符号付き整数の最大値は素数 - Rustで任意精度整数演算 |
🦒 |
tech |
|
true |
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には組み込みの整数型として usize, isize はここでは除外]がそれぞれ符号付き・符号なしで備わっています^[Integer Types - The Rust Programming Language]。そのうち符号付き整数は、他の多くの言語と同様、2の補数によって負の数が表現されます。したがって、ビット数
各々の型で表現可能な最大値を実際に見てみましょう。型の定数として値は（128ビットの場合） i128::MAXから取ることができます。
:::details 一応コード
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);
}:::
| 桁数(十進) | ||
|---|---|---|
表には十進法表記を示しました。語呂合わせでも作ってちょっと覚えてみてもいいかもしれない、と思わせるような数字の並びですが、二進法と十六進法なら自明です：
$$ 2^n-1 = (\underbrace{111\cdots 1}{n,個}){2} =(7\underbrace{{\rm f,f,f}\cdots{\rm f}}{m-1,個}){16} $$
ただし、$n=4m$ です。二進法では全ての桁の数が1であり、基数2におけるレピュニット数になっています。
表には素因数分解も示しました。$n=16,,64$ の等号の右辺が素因数分解です。一方、等号を含まない行、$n=8,,32,,128$ ではこれらの値は素数になっています。
つまり、
確かめてみましょう。
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世紀フランスのカトリックの司祭・数学者メルセンヌ^[Marin Mersenne(1588-1648) https://en.wikipedia.org/wiki/Marin_Mersenne]に由来します。彼は
の11個^[A109461 Mersenne's original list of "Mersenne" exponents.] に限られると予想しました。後の時代になって彼のリストには
- 3つの抜け：$n=61, 89, 107$ に対しても
$M_n$ は素数 - 2つの誤り：$n=67, 257$ に対して
$M_n$ は実際には合成数
があることが示されましたが、彼に敬意を表して今日でもその名が残っています。なお、$n\leq 257$ の範囲で
を得る方法は後で見ることになります。
さて、$M_n$ が素数であるためには
という分解^[$2^{ab}-1$ は
しかし逆は成り立ちません。すなわち、素数
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
);
}
}
}
https://github.com/roiban1344/mersenne_primes/blob/main/packages/mersenne_factor_naive/src/main.rs
出力：
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列目: 素因数たち（昇順）
を意味します。この範囲で既に
:::details メルセンヌ数の素因数と、素数が無限に存在することについて
特に
特筆すべきことに、この事実は素数が無限に存在することの証明にもなります。最大の素数
参考：M.アイグナー，G.M.ツィーグラー， 訳 蟹江幸博『天書の証明』，丸善出版，2002, 第1章. :::
メルセンヌ素数は完全数^[完全数 - Wikipedia。A000396] との密接な繋がりから関心が向けられてきました。
完全数とは、自身を除く約数の和が自身と一致する非負整数のことです。小さいほうから
と、確かに
:::details 証明 偶数の完全数とメルセンヌ素数が一対一対応すること
正整数
したがって
一方、$n$ が完全数であることから、$\sigma(n)=2n=2^{a+1}b$。2つの式を合わせて、
左辺は整数だから、右辺第2項
よって
オイラーに帰せられる後者の証明は非常に面白いです。
参考：高木貞治『初頭整数論講義 第2版』，共立出版，1971． :::
この事実により、古代ギリシャから知られていた偶完全数の探索はメルセンヌ数探索と表裏一体になりました。一方、奇数の完全数は未だ一つも知られていないうえ存在しないとも証明されておらず、数論における有名な未解決問題になっています。
本題に入ります。1876年、リュカ^[Édouard Luca(1842-1891) https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%89%E3%82%A5%E3%82%A2%E3%83%BC%E3%83%AB%E3%83%BB%E3%83%AA%E3%83%A5%E3%82%AB]は以下の判定法を利用し、$M_{127}=2^{127}-1$ が素数であることを示しました。
(リュカテスト)
確かに
一方
確かに
リュカ・テストが正しいことの証明は、その核心でフィボナッチ数（特に整除性と周期性）が関わり、初等的でありながら全く自明でない結果が導かれる、非常に興味深いものです。
:::details 証明 リュカテスト 証明の大筋は参考文献[中村]p.82によるものです。
フィボナッチ数列を
リュカ数列を
によってそれぞれ再帰的に定義する。なお、$L_{i}=F_{i+1}+F_{i-1}$。このとき、
とおくと、
証明：帰納法による。
証明：帰納法。
証明：
ほぼ同様にして
組み合わせると、
を得る。$\square$
正整数
証明: まず、定義から
が成り立つ。この関係式を繰り返し適用して、
一般に、
であることから、$F_{e-1}$ は
により、$F_{r}\equiv 0 \mod q$ が成り立つことになり、$e$ の最小性に反する。$\square$
準備ができたので以下本証明に入ります。
$$ M_p=2^{4n+3}-1=8\cdot 16^n-1\equiv 2 \mod 5. $$
補題(2), (3)から、
一般に奇素数
のように振る舞うため
このとき
に
補題(3)から
:::
確率的な素数判定法が実用に耐えるのは、判定対象の整数の 対数に関する 多項式時間のアルゴリズムであるからです。一方、例えば試し割りによる方法だと、判定対象の整数それ自体に関する多項式時間です。
計算量の観点から見たリュカテストの著しい性質は、決定的であるにもかかわらず確率的な素数判定法と同等かそれ以上に高速であることです。
この非常に効率的な判定法を得たリュカは、1876年に
とはいえ実際に手計算してみようとすると分かりますが、$M_{19}$ くらいでも相当骨が折れます^[本気でやるならカラツバ法を使うのはいいアイデアだと思います。実演してみませんか？]。歴史に名を遺すとはどういうことかよく分かります。
リュカテストを実装します。$R_{i}$ はおよそ前の項の2倍の桁数を取り急速に増加しますが、$M_p$ で剰余を取ればよいので桁数は有限に留まります。それでも、$M_{127}$ の場合には i128 で扱うと2乗を取る過程でオーバーフローが起こります。
そこで任意桁数の整数演算を導入します。これを実現するRustのクレートとそれぞれの特徴を以下に挙げます^[num-bigintの"Alternatives"の項の表より。また、やや古いがRedditの投稿。]
- num-bigint: 純Rust製。使いやすい。
- rug: GMPライブラリ等のバンドル。使いやすい。高速。
- rust-gmp: 同上（のようだがあまりドキュメントが充実しておらず詳細不明）。
- ramp: Nightlyビルド限定。
- apint: 純Rust製、未完。
今回はこの中から rug を採用することにしました。
:::details num-bigintについて
num-bigintは直感的に扱えて使いやすいですが、どうもバグらしき挙動を示します。$2^{2112}-1$ の自乗を取ろうとすると panic! します。
https://github.com/roiban1344/mersenne_primes/blob/main/packages/num_bigint/product_panic/src/main.rs
よりによって
もっとも性能ではrugが上回っているため、num-bigintを使う必要はあまりなさそうです。
:::
rugはGNUに属す3つの任意精度演算用ライブラリ(GMP, MPFR, MPC)のラッパーです。任意精度整数はこの中のGMPによって提供されます。Mathematicaの内部でも使われるなど^[内部実装について—Wolfram言語ドキュメント]、歴史と信頼性のあるライブラリです。任意精度整数演算を行うなら最上の選択肢の一つであることは間違いないでしょう。
rug は低レベルのFFIを提供するクレートである gmp-mpfr-sys に依存するため、ビルドにはCargo.tomlにdependenciesを追加する以外にも操作が必要になることがあります。詳細は gmp-mpfr-sys のドキュメントを参照してください^[筆者はWindowsのホスト上での構築に挫折し、WSL2に移行しました。Rust公式のDockerイメージ(rust:latest)を使っても特別な操作が不要であることを確認しています。]
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);
}
}
}
}
https://github.com/roiban1344/mersenne_primes/blob/main/packages/rug/lucas_test/src/main.rs
primesはただの篩です。lucas_testがリュカテストの実装です。
初期化や演算にInteger::fromを介していることを除けばプリミティブ型で書くのとそう変わりません。rugの使い方とは別の面で実装上の注意が2つだけあります。
まず、メルセンヌ数に関する剰余演算は汎用の剰余演算子である % より高速に実行できます。
から、$a=0$ でない限り
これは二進法版の九去法のようなものです。実装上は
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 はとりません。また、$M_p$ に関して
さて、実行すると
3
7
19
31
107
127
607
1279
2203
4423
たった10個！ この区間には
リュカテストの明らかな弱点は
(リュカ・レーマーテスト)
初項を
リュカ・レーマーテストもリュカテストと同様、フィボナッチ・リュカ数列に対応する2階線形回帰数列を考えることで証明が可能ですが、より簡潔な証明がBruceによって与えられています^[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.、]。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);
}
}
}
}https://github.com/roiban1344/mersenne_primes/blob/main/packages/rug/lucas_lehmer_test/src/main.rs
判定対象の上限を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 - Getting Started、またGIMPSの始め方 - Qiitaを参照してください。
GIMPSの探索プログラムには前処理や種々の最適化、検算が含まれますが、核心部分ではやはりリュカ・レーマーテストが利用されています^[2018年からはフェルマーテストも採用されているようです。Is there a better method than Lucas-Lehmer?]。先述の確率的素数判定法や素因数分解アルゴリズムがコンピューター時代以後に洗練されてきたことを考えると、1世紀半に渡って最前線で利用されているリュカ・レーマーテストの強力さは驚嘆に値します。
リュカ・レーマーテストは素因数に関する情報を与えてくれません。素因数分解と素数判定が質的に異なる問題であることを示す実例の一つです。
メルセンヌが誤って素数だと予想した2つの合成数のうち、$M_{67}$ は試し割りでも素因数分解可能です。コンピュータ時代以前の偉業として、コール^[Frank Nelson Cole(1861-1926) Frank Nelson Cole] が下記の十進で9桁と12桁の素数への分解を行っています。
一方
https://github.com/roiban1344/mersenne_primes/blob/main/packages/rug/mersenne_factor_rho/src/main.rs
数学計算ソフトSageMathはオンライン版SageMathCellが利用可能です。ここに単に
factor(2^257-1)
と打ち込んで"Evaluate"を押すと、10秒程度で素因数分解が得られます。$M_{67}$ と比較しましょう。
今日の素因数分解アルゴリズムとしては一般数体篩法が最高効率です。類似のアイデアに基づくもう少し簡単なアルゴリズムには二次篩法もあります。これらの実装によるメルセンヌ数の素因数分解は今後の課題です。
メルセンヌ数の歴史と基本的な性質を概観し、その素数判定法であるリュカ・レーマーテストを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]はまだ
メルセンヌ数に関する最新のデータや歴史はこれらが詳しいです。本記事の数学的な内容や歴史はほとんどここにも載っています。