概要

$2^n-1$ 型の数はメルセンヌ数と呼ばれ、更に素数である場合にメルセンヌ素数といいます。本記事では、メルセンヌ数に対する高速な素数判定法であるリュカ・レーマーテストを、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には組み込みの整数型として $8,\,16,\,32,\,64,\,128$ ビット整数^[1]がそれぞれ符号付き・符号なしで備わっています^[2]。そのうち符号付き整数は、他の多くの言語と同様、2の補数によって負の数が表現されます。したがって、ビット数 $n = 8, 16, 32, 64, 128$ に対して、表現可能な値の範囲は $-2^{n-1}$ から $2^{n-1}-1$ （両端を含む）に渡ります。

各々の型で表現可能な最大値を実際に見てみましょう。型の定数として値は（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);
}

$n$	$2^{n-1}-1$	桁数(十進)
$8$	$127$	$3$
$16$	$32767=7 \cdot 31 \cdot 151$	$5$
$32$	$2147483647$	$10$
$64$	$9223372036854775807=7^2\cdot 73\cdot 127\cdot 337\cdot 92737\cdot 649657$	$19$
$128$	$170141183460469231731687303715884105727$	$39$

表には十進法表記を示しました。語呂合わせでも作ってちょっと覚えてみてもいいかもしれない、と思わせるような数字の並びですが、二進法と十六進法なら自明です：

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$ ではこれらの値は素数になっています。

つまり、 $8,\,32,\,128$ ビット符号付き整数の最大値は素数なのです。

確かめてみましょう。

/packages/naive_factorization/src/main.rs

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 が出力されるまでには $\sqrt{2^{127}-1} \simeq 2^{63}$ 回のループ、つまり64ビット整数全てに渡るのとほぼ同じ回数のループを回すことになるためです。long longのループを回した苦い経験が蘇ります。

ところがこの $2^{127}-1$ という数は、今から遡ること実に145年前、西暦1876年に素数であることが証明されているのです。当然、電子計算機が発明される以前の時代です。

メルセンヌ数とは

歴史的背景について。

非負整数 $n$ に対して $2^{n}-1$ 型の整数はメルセンヌ数（Mersenne number）と呼ばれ、記号 $M_n=2^n-1$ で表されます。特に素数のメルセンヌ数はメルセンヌ素数（Mersenne prime）といいます。すなわち、 $8, 32, 128$ ビット符号付き整数の最大値はそれぞれ $M_{7}, M_{31}, M_{127}$ で表されるメルセンヌ素数です。

「メルセンヌ数」の呼び名は、16世紀フランスのカトリックの司祭・数学者メルセンヌ^[3]に由来します。彼は $2^n-1$ が素数となるのは $n\leq 257$ の範囲で

n=2,3,5,7,13,17,19,31,67,127,257

の11個^[4] に限られると予想しました。後の時代になって彼のリストには

3つの抜け： $n=61, 89, 107$ に対しても $M_n$ は素数
2つの誤り： $n=67, 257$ に対して $M_n$ は実際には合成数

があることが示されましたが、彼に敬意を表して今日でもその名が残っています。なお、 $n\leq 257$ の範囲で $M_n$ が素数になる $n$ たちの正しいリスト：

n=2,3,5,7,13,17,19,31,61,89,107,127

を得る方法は後で見ることになります。

さて、 $M_n$ が素数であるためには $n$ が素数であることが必要です。なぜなら、 $n$ が合成数であるとき2つの非負整数 $1<a,b<n$ によって $n=ab$ と表され、

M_n = 2^{ab}-1 = (2^a-1)(2^{a(b-1)}+2^{a(b-2)+\cdots + 2 + 1})

という分解^[5]から、 $M_n$ は $1 < 2^a-1 < M_n$ を約数に持つためです。

しかし逆は成り立ちません。すなわち、素数 $p$ に対して $M_p$ が素数になるとは限りません。最小の反例が $p=11$ の場合で、 $M_{11}=2047=23\cdot 89$ と分解されます。以降、本記事では記号 $M_p$ を「素数 $p$ に対するメルセンヌ数（ $M_p$ 自身が素数とは限らない）」の意味で用いることにします。

$p\leq 19$ では $p=11$ を除いて全て $M_p$ が素数となるため、「それなりの」頻度で素数になりそうに感じられますが、 $p$ が大きくなるほどメルセンヌ素数は急速に疎らになっていきます。

64ビット整数の範囲の $M_p$ を素因数分解してみます。この範囲ならナイーブな実装による $M_p$ の素因数分解も数秒で完了します。

/packages/mersenne_factor_naive/src/main.rs

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列目: 素因数たち（昇順）

を意味します。この範囲で既に $p=37, 41, 43, 47, 53, 59$ の連続する6つに対して $M_p$ が合成数になっています。

メルセンヌ数の素因数と、素数が無限に存在することについて

$M_p$ の素因数を眺めると、素因数として「小さい」素数 ( $3, 7, 13, \cdots$ 等）がほとんど現れないことに気付きます。実際、 $M_p$ の素因数は ( $p$ の倍数) $+ 1$ の形に限られる（したがって $p+1$ より大きい）ことが以下のように示されます：

$M_p$ の最小の素因数を $q$ とすると、 $2^p \equiv 1 \mod q$ 。これは $\mathbb{Z}/q\mathbb{Z}$ の乗法群 $(\mathbb{Z}/q\mathbb{Z})^*=\mathbb{Z}/q\mathbb{Z}-\{0\}$ において $2$ の位数が $p$ であることを意味する。したがって $p$ は $(\mathbb{Z}/q\mathbb{Z})^*$ の位数 $q-1$ を割り切り、正整数 $k$ によって $q=kp+1$ と表される。 $\square$

特に $p$ が奇素数の場合、ありうる最小の素因数は $2p+1$ で、 $M_{11}=2047,\, M_{23}=8388607$ がそれぞれ実際に $2\cdot 11+1=23, 2\cdot23+1=47$ を素因数に持ちます。

特筆すべきことに、この事実は素数が無限に存在することの証明にもなります。最大の素数 $p$ が存在するとすれば、 $M_p$ が $p+1$ 以上の素因数を持って矛盾するためです！

参考：M.アイグナー，G.M.ツィーグラー，訳蟹江幸博『天書の証明』，丸善出版，2002, 第1章.

メルセンヌ数と完全数

メルセンヌ素数は完全数^[6] との密接な繋がりから関心が向けられてきました。

完全数とは、自身を除く約数の和が自身と一致する非負整数のことです。小さいほうから $1,2,3$ 番目の完全数の実例はそれぞれ $n=6,28,496$ で、

\begin{align*} 6&=1+2+3,\\ 28&=1+2+4+7+14,\\ 496&=1+2+4+8+16+31+62+124+248 \end{align*}

と、確かに $n=(n$ 自身を除く $n$ の約数の和)となっています。実は、 $M_p$ が素数なら $2^{p-1}M_p$ は完全数であり、逆に偶数の完全数はこの形に限られることが知られています。

証明　偶数の完全数とメルセンヌ素数が一対一対応すること

$M_p$ が素数なら $n=2^{p-1}M_p$ は完全数であること

正整数 $m$ に対して $m$ の約数の和を $\sigma(m)$ とする（約数関数）。 $\sigma(m)=2m$ なら $m$ は完全数である。仮定から $M_p$ は3以上の素数だから、約数関数の性質から、

\sigma(n) = \frac{2^{p}-1}{2-1}\frac{M_p^2-1}{M_p-1} = M_p\cdot 2^p = 2n.

したがって $n$ は完全数。 $\square$

$n$ が偶数の完全数ならメルセンヌ素数 $M_p$ が存在して $n=2^{p-1}M_p$ が成り立つこと

$n$ は偶数だから、正整数 $a$ と正の奇数 $b$ によって $n=2^a b$ と（一意に）表される。約数関数の性質から、

\sigma(n) = \sigma(2^a)\sigma(b) = (2^{a+1}-1)\sigma(b).

一方、 $n$ が完全数であることから、 $\sigma(n)=2n=2^{a+1}b$ 。2つの式を合わせて、

\sigma(b) = b + \frac{b}{2^{a+1}-1}.

左辺は整数だから、右辺第2項 $b/(2^{a+1}-1)=c$ も整数でなくてはならない。したがって分母の $2^{a+1}-1=d$ は $b$ の約数。さらに $a\geq 1$ だから、 $d>1$ と $c<b$ が成り立つ。 $\sigma(b)=b+c$ は「『 $b$ の約数の和』は『 $b$ 自身と $b$ より小さい $b$ の約数 $c$ の和』に等しい」ことを意味する。

よって $b$ は $b$ 自身と $c$ のふたつの約数しか持たないため素数であり、 $c=1$ 。さらに、 $b=2^{a+1}-1$ が素数であるためには $a+1$ は素数でなくてはならず（既述）、これを $p$ とおくと $n=2^{p-1}M_p$ と表される。 $\square$

オイラーに帰せられる後者の証明は非常に面白いです。

参考：高木貞治『初頭整数論講義　第2版』，共立出版，1971．

この事実により、古代ギリシャから知られていた偶完全数の探索はメルセンヌ数探索と表裏一体になりました。一方、奇数の完全数は未だ一つも知られていないうえ存在しないとも証明されておらず、数論における有名な未解決問題になっています。

リュカテスト

本題に入ります。1876年、リュカ^[7]は以下の判定法を利用し、 $M_{127}=2^{127}-1$ が素数であることを示しました。

(リュカテスト)
$p$ を $4n+3$ 型の素数とする。 $R_1=3,\, R_i=R_{i-1}^2-2\,(i\geq 2)$ で再帰的に定まる整数列 $R_i$ について、 $R_{p-1}\equiv 0 \mod M_p$ なら $M_p$ は素数、そうでなければ合成数。

$M_{p} \mod \square$ ではなく $\square \mod M_p$ であることに注意してください。 $M_p$ に対して一度も割り算を施さず、非自明な約数を一つも得ることなく素数か合成数か判定できるのです！

$M_p$ が実際に素数になるケース、 $p=7,\,M_7=127$ で試します。以下の式で $\mod M_p$ は省略します。

\begin{align*} R_1 &\equiv 3,\\ R_2 &\equiv 3^2-2\equiv 7,\\ R_3 &\equiv 7^2-2\equiv 47,\\ R_4 &\equiv 47^2-2\equiv 48,\\ R_5 &\equiv 48^2-2\equiv 16,\\ R_6 &\equiv 16^2-2\equiv 0. \end{align*}

確かに $R_{p-1} \equiv 0 \mod M_p$ になりました。

一方 $M_p$ が合成数になるケース、 $p=11,\,M_{11}=2047=23\cdot 89$ で試します。

\begin{align*} R_1 &\equiv 3,\\ R_{2} &\equiv 3^2-2 \equiv 7,\\ R_{3} &\equiv 7^2-2 \equiv 47,\\ R_{4} &\equiv 47^2-2 \equiv 160,\\ R_{5} &\equiv 160^2-2 \equiv 1034,\\ R_{6} &\equiv 1034^2-2 \equiv 620,\\ R_{7} &\equiv 620^2-2 \equiv 1609,\\ R_{8} &\equiv 1609^2-2 \equiv 1471,\\ R_{9} &\equiv 1471^2-2 \equiv 160,\\ R_{10} &\equiv 160^2-2 \equiv 1034. \end{align*}

確かに $R_{p-1}\not\equiv 0 \mod M_p$ となりました。特に $R_{4}\equiv R_{9}\equiv 160$ でループに入っています。

リュカ・テストが正しいことの証明は、その核心でフィボナッチ数（特に整除性と周期性）が関わり、初等的でありながら全く自明でない結果が導かれる、非常に興味深いものです。

証明　リュカテスト

証明の大筋は参考文献[中村]p.82によるものです。

補題(1)

フィボナッチ数列を

F_0=F_1=1,\,F_{i}=F_{i-1}+F_{i-2},

リュカ数列を

L_0=2,\, L_1=1,\, L_{i}=L_{i-1}+L_{i-2}

によってそれぞれ再帰的に定義する。なお、 $L_{i}=F_{i+1}+F_{i-1}$ 。このとき、

\alpha = \frac{1+\sqrt{5}}{2},\, \beta = \frac{1-\sqrt{5}}{2}

とおくと、

F_i = \frac{\alpha^i-\beta^i}{\sqrt{5}},\, L_i = \alpha^i + \beta^i.

証明：帰納法による。

補題(2)

$i\geq 1$ について、

R_{i} = \alpha^{2^i} + \beta^{2^i} = L_{2^i}.

証明：帰納法。

補題(3)

$q\equiv 1\mod 5$ であるような奇素数 $q$ について、 $F_{q-1} \equiv 0 \mod q$ .

$q\equiv 2\mod 5$ であるような奇素数 $q$ について、 $F_{q+1} \equiv 0,\, L_{q+1} \equiv -2\mod q$ .

証明：

$5$ ではない奇素数 $q$ について、

\begin{align*} F_{q} &= \frac{1}{2^q\sqrt{5}}\sum_{k=0}^q \binom{q}{k} (\sqrt{5})^k - (-\sqrt{5})^k\\ &=\frac{1}{2^q\sqrt{5}}\sum_{kは奇数,\, 1\leq k \leq q} \binom{q}{k} 2\cdots 5^{\frac{k-1}{2}},\\ \therefore 2^{q-1} F_{q} &= \sum_{kは奇数,\, 1\leq k \leq q} \binom{q}{k} \, 5^{\frac{k-1}{2}} \end{align*}

$\mod q$ を取ると、 $1\leq k\leq q-1$ で $\binom{q}{k}\equiv 0\mod q$ であることと $2^{q-1}\equiv 1 \mod q$ （ $\because$ フェルマーの小定理）から和は $k=q$ の項のみが残り、法 $5$ の平方剰余に関するオイラーの規準から、

F_{q} \equiv 5^{\frac{q-1}{2}} \equiv \left(\frac{q}{5}\right) \equiv \left\{\begin{array}{cc} 1 & {\rm if}\, q\equiv \pm 1 \mod 5\\ -1 & {\rm if}\, q\equiv \pm 2 \mod 5 \end{array}\right. \mod q.

ほぼ同様にして

\begin{align*} F_{q+1} \equiv \left\{\begin{array}{cc} 1 & {\rm if}\, q\equiv \pm 1 \mod 5\\ 0 & {\rm if}\, q\equiv \pm 2 \mod 5 \end{array}\right. \mod q. \end{align*}

組み合わせると、

\begin{align*} F_{q-1} = F_{q+1}-F_{q} &\equiv \left\{\begin{array}{cc} 0 & {\rm if}\, q\equiv \pm 1 \mod 5\\ 1 & {\rm if}\, q\equiv \pm 2 \mod 5 \end{array}\right. \mod q, \\ L_{q+1} = 2F_{q} + F_{q+1} &\equiv \left\{\begin{array}{cc} 3 & {\rm if}\, q\equiv \pm 1 \mod 5\\ -2 & {\rm if}\, q\equiv \pm 2 \mod 5 \end{array}\right. \mod q. \end{align*}

を得る。 $\square$

補題(4)

正整数 $n$ に対して $F_{i}\equiv 0 \mod n$ となる最小の正整数 $i$ が存在して $e$ であるとする。正整数 $j$ が $F_{j} \equiv 0 \mod n$ を満たすなら、 $e$ は $j$ を割り切る。

証明:
まず、定義から

\begin{align*} \begin{pmatrix} F_{i} & F_{i-1}\\ F_{i-1} & F_{i-2} \end{pmatrix} = {\bm A} \begin{pmatrix} F_{i-1} & F_{i-2}\\ F_{i-2} & F_{i-3} \end{pmatrix},\, {\bm A} := \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix} \end{align*}

が成り立つ。この関係式を繰り返し適用して、

\begin{align*} \begin{pmatrix} F_{i+1} & F_{i}\\ F_{i} & F_{i-1} \end{pmatrix} = {\bm A}^{i}. \end{align*}

$F_{j}\equiv 0 \mod q$ を満たす正整数 $j$ が $e$ で割り切れないと仮定して矛盾を導く。仮定から、整数 $k$ と $1\leq r\leq e-1$ を満たす整数 $r$ が存在して $j=ke+r$ が成り立つ。この $j$ に対して、

\begin{align*} \begin{pmatrix} F_{j+1} & F_{j}\\ F_{j} & F_{j-1} \end{pmatrix} &={\bm A}^{ke+r}\\ &=\left(({\bm A}^e)\right)^k {\bm A}^r\\ &= \begin{pmatrix} F_{e+1} & F_{e}\\ F_{e} & F_{e-1} \end{pmatrix}^k \begin{pmatrix} F_{r+1} & F_{r}\\ F_{r} & F_{r-1} \end{pmatrix}. \end{align*}

$\mod n$ を取ると、

\begin{pmatrix} F_{j-1} & 0\\ 0 & F_{j-1} \end{pmatrix} \equiv \begin{pmatrix} F_{e-1} & 0\\ 0 & F_{e-1} \end{pmatrix}^k \begin{pmatrix} F_{r+1} & F_{r}\\ F_{r} & F_{r-1} \end{pmatrix}\mod n.

一般に、

{\rm gcd}(F_{i}, F_{i+1}) = {\rm gcd}(F_{i}, F_{i+1}-F_{i}) = {\rm gcd}(F_{i}, F_{i-1}) = \cdots = {\rm gcd}(F_{0}, F_{1}) = 1

であることから、 $F_{e-1}$ は $F_{e}$ とは互いに素で、 $\mathbb{Z}/m\mathbb{Z}$ の乗法群における逆元が存在する。したがって

\begin{pmatrix} F_{r+1} & F_{r}\\ F_{r} & F_{r-1} \end{pmatrix} \equiv F_{e-1}^{-k}F_{j} \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \mod n

により、 $F_{r}\equiv 0 \mod q$ が成り立つことになり、 $e$ の最小性に反する。 $\square$

準備ができたので以下本証明に入ります。

リュカテストが正しいことの証明

(A) $M_p$ が素数なら $R_{p-1} \equiv 0 \mod M_p$ であること

$p=4n+3$ と表せることから、

M_p=2^{4n+3}-1=8\cdot 16^n-1\equiv 2 \mod 5.

補題(2), (3)から、

L_{M_p+1} = L_{2^p} = R_{p} \equiv -2 \mod M_p.

一般に奇素数 $q$ に対して、 $R_{k}\equiv -2 \mod q$ となる $k$ が存在するなら、その前後で

\begin{align*} R_{k-1} &\equiv 0 \mod q,\\ R_{k} = R_{k-1}^2-2&\equiv -2 \mod q,\\ R_{k+1} = R_{k}^2-2 &\equiv 2 \mod q,\\ R_{k+2} = R_{k+1}^2-2 &\equiv 2 \mod q,\\ &\vdots \end{align*}

のように振る舞うため $R_{k-1}$ のみ唯一 $q$ で割り切れる。この奇素数 $q$ を $M_p$ とみなすと、 $R_{p-1}$ が $M_p$ で割り切れて、しかも他の $R_i$ は $M_p$ では割り切れないことが分かる。

(B) $R_{p-1} \equiv 0 \mod M_p$ なら $M_p$ は素数であること

$M_p$ の素因数の1つをとって $q$ とする。 $q$ は $2,5$ のどちらでもない。

このとき $R_{p-1} \equiv 0 \mod q$ だが、(A)で行った観察から、 $R_{k} \equiv 0 \mod q$ を満たす唯一の $k$ が $k=p-1$ になる。ところで、

\begin{align*} F_{2^{k}} &= F_{2^{k-1}}L_{2^{k-1}} & (\because 補題(1)を利用)\\ &=F_1L_1L_2L_4\cdots L_{2^{k-1}}\\ &=R_{1}R_{2}\cdots R_{k-1} \end{align*}

に $k=p$ を代入すると、仮定から $F_{2^p} \equiv 0 \mod q$ であり、補題(4)から、 $F_{e} \equiv 0 \mod q$ を満たす最小の $e$ は $2^p$ の約数。ところが上述の理由から $e=2^p$ 。

補題(3)から $F_{q+1}\equiv 0$ または $F_{q-1}\equiv 0$ であり、再び補題(4)から、 $q+1$ か $q-1$ のどちらか一方は $2^p$ で割り切れる。ところが $q$ は $M_p$ の約数だから、 $q \leq M_p = 2^p -1$ であり、 $q+1=2^p$ しかありえない。すなわち $M_p$ の素因数は $q=2^p-1=M_p$ 自身に限られ、 $M_p$ は素数である。 $\square$

$M_{127}$ より遥かに大きな数千ビットの素数はRSA暗号などのキーとして日常的に用いられていますが、そこで生成されるランダムな整数の素数性の判定には、確率的なアルゴリズムが採用されます。ミラー・ラビンテストは実用的な確率的素数判定法の1つです。ここでいう「確率」は「ハッシュ値の衝突確率は現実的には無視できる」という文脈における確率と似ていて、「実際には合成数であるのに素数であると間違って判定してしまう（擬素数）」確率は十分無視できる精度まで任意に小さくすることができます。

確率的な素数判定法が実用に耐えるのは、判定対象の整数の 対数に関する 多項式時間のアルゴリズムであるからです。一方、例えば試し割りによる方法だと、判定対象の整数それ自体に関する多項式時間です。

計算量の観点から見たリュカテストの著しい性質は、決定的であるにもかかわらず確率的な素数判定法と同等かそれ以上に高速であることです。

この非常に効率的な判定法を得たリュカは、1876年に $M_{127}$ が素数であることを手計算により示しました。その後電子計算が登場して次に小さな $M_{521}$ が（後述のリュカ・レーマーテストにより）素数であることが示されるまでの約75年間^[8]、 $M_{127}$ は人類に知られた最大の素数として君臨しました。128ビット符号付き整数の最大値はそんな歴史を持つ数です。

とはいえ実際に手計算してみようとすると分かりますが、 $M_{19}$ くらいでも相当骨が折れます^[9]。歴史に名を遺すとはどういうことかよく分かります。

任意精度演算用クレート `rug`

リュカテストを実装します。 $R_{i}$ はおよそ前の項の2倍の桁数を取り急速に増加しますが、 $M_p$ で剰余を取ればよいので桁数は有限に留まります。それでも、 $M_{127}$ の場合には i128 で扱うと2乗を取る過程でオーバーフローが起こります。

そこで任意桁数の整数演算を導入します。これを実現するRustのクレートとそれぞれの特徴を以下に挙げます^[10]

num-bigint: 純Rust製。使いやすい。
rug: GMPライブラリ等のバンドル。使いやすい。高速。
rust-gmp: 同上（のようだがあまりドキュメントが充実しておらず詳細不明）。
ramp: Nightlyビルド限定。
apint: 純Rust製、未完。

今回はこの中から rug を採用することにしました。

num-bigintについて

num-bigintは直感的に扱えて使いやすいですが、どうもバグらしき挙動を示します。 $2^{2112}-1$ の自乗を取ろうとすると panic! します。

よりによって $M_{21701}$ のリュカ・レーマーテストの最後のステップでこけた事象が最初の遭遇でした。そこから何か原因が見つけられるかと思いましたが、まだ何も分かっていません。

もっとも性能では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)も可能です。

リュカテストの実装

/packages/rug/lucas_test/main.rs

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つだけあります。

まず、メルセンヌ数に関する剰余演算は汎用の剰余演算子である % より高速に実行できます。

$M_n$ で剰余を取る対象の正整数 $m$ を $m=2^n a + b$ と表します。 $0 \leq b \leq 2^n-1$ に制限すると $a, b$ は一意に決まり、 $b$ が下位 $n$ ビット、 $2^n a$ がその上位ビットの値です。このとき

\begin{align*} m &= 2^n a + b\\ &= (2^n-1)a + a + b\\ &\equiv a + b \mod M_n \end{align*}

から、 $a=0$ でない限り $m > a+b$ であるので、 $a=0$ 　になるまで $m$ へ $a+b$ の再代入を繰り返せば $M_n$ に関する剰余が決まります。ただしちょうど $M_n$ に一致する場合のみ $0$ への置き換えが必要です。

これは二進法版の九去法のようなものです。実装上は $a$ が右シフト、 $b$ が $M_n$ との AND 演算を行うことでそれぞれ計算できます：

        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$ が非平方剰余であることから、1 にもなりません。したがってループ中では $s$ が負値に転じることは起こらず、負値に関する剰余の振る舞いを考慮する必要はありません。

さて、実行すると $10^4$ 以下の $4n+3$ 型の素数 $p$ で $M_p$ が素数になるものが以下のようにリストアップされます。

たった10個！　この区間には $619$ 個の $4n+3$ 型の素数があるため、かなり疎らだと言えます。

リュカ・レーマーテスト

リュカテストの明らかな弱点は $4n+1$ 型の素数に関する判定が行えないことです。仮に $M_5=31$ に対して適用すると、合成数と誤判定されます。

$4n+1$ 型素数を判定対象に含めるには、リュカテストにほんの少しの改良を加えるだけで済みます。それがリュカ・レーマーテストです：

(リュカ・レーマーテスト)
$S_1=4, S_{i}=S_{i-1}^2-2 (i\leq 2)$ によって数列 $S_i$ を再帰的に定義する。奇素数 $p$ に対して $S_{p-1}$ が $M_p$ で割り切れれば $M_p$ は素数、割り切れなければ合成数。

初項を $4$ に置き換えるだけです！

リュカ・レーマーテストもリュカテストと同様、フィボナッチ・リュカ数列に対応する2階線形回帰数列を考えることで証明が可能ですが、より簡潔な証明がBruceによって与えられています^[13]。Wikipediaに掲載されているものがこれに基づくため、詳細はそちらに譲ります: Proof of correctness

先のリュカテストに数行の変更を加えて、リュカ・レーマーテストの実装は以下のようになります。

/packages/mersenne_primes/packages/rug/lucas_test/src/main.rs

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にまで上げて実行：

ものの数分の実行で26番目のメルセンヌ素数 $M_{23209}$ まで見つかりました。特にメルセンヌが述べた $p\leq 257$ の範囲はこれでカバーされました。

メルセンヌ素数探索の歴史が表にまとめられています。いま見つかった中で最大の $M_{23209}$ は、歴史的には1979年の記録です。

ちなみに24番目の $M_{19937}$ はちょうど50年前、1971年の発見です。論文^[14]には、IBM 360/91に実装したプログラムがこの値に対して最初にリュカ・レーマーテストをパスしたのが3月4日の夕方であることまで書かれています。更には検算に30分かかっていたことも。

$M_{19937}$ のもう一つの注目すべき点は、質の高い擬似乱数発生法であるメルセンヌ・ツイスタ（の実装MT19937）で広く利用されている値であることです。ある意味では $M_{127}$ と同じくらい身近な数かもしれません。

現在の世界記録は2018年に発見された $M_{82589933}$ です。小さいほうから数えて何番目のメルセンヌ素数であるかは未確定ですが、人類に知られたメルセンヌ素数としては51番目です。既知の最大の素数でもあります。 $82589933 < 10^8 < 2^{32}$ ですから、まだ $2^{2^{32}}$ は超えていないことになります。

51個のうち半分までの発見は家庭用PCで瞬時に再現できる一方、残り半分は今日までの人類のメルセンヌ素数探索の歴史全てがかけられてようやく見つかるのです。

GIMPS

既知最大のメルセンヌ素数 $M_{82589933}$ の発見^[15]を含め過去四半世紀に渡り、メルセンヌ素数発見の世界記録はGIMPS（Great Internet Mersenne Prime Search）という、分散コンピューティングプロジェクトが保持しています。

GIMPSトップページ

このプロジェクトには誰でもアカウントを作って貢献することができます。詳しくは公式によるGIMPS - Getting Started、またGIMPSの始め方 - Qiitaを参照してください。

GIMPSの探索プログラムには前処理や種々の最適化、検算が含まれますが、核心部分ではやはりリュカ・レーマーテストが利用されています^[16]。先述の確率的素数判定法や素因数分解アルゴリズムがコンピューター時代以後に洗練されてきたことを考えると、1世紀半に渡って最前線で利用されているリュカ・レーマーテストの強力さは驚嘆に値します。

素因数分解

リュカ・レーマーテストは素因数に関する情報を与えてくれません。素因数分解と素数判定が質的に異なる問題であることを示す実例の一つです。

メルセンヌが誤って素数だと予想した2つの合成数のうち、 $M_{67}$ は試し割りでも素因数分解可能です。コンピュータ時代以前の偉業として、コール^[17] が下記の十進で9桁と12桁の素数への分解を行っています。

一方 $M_{257}$ は、試し割りよりは強力ながら実装が容易なポラードのρ法でも10分程度の実行では素因数を見つけられませんでした。

数学計算ソフトSageMathはオンライン版SageMathCellが利用可能です。ここに単に

factor(2^257-1)

と打ち込んで"Evaluate"を押すと、10秒程度で素因数分解が得られます。 $M_{67}$ と比較しましょう。

\begin{align*} M_{67} &= 193707721 \cdot 761838257287,\\ M_{257} &= 535006138814359 \cdot 1155685395246619182673033 \cdot 374550598501810936581776630096313181393 \end{align*}

$M_{67}$ の最小素因数が十進9桁であるのに対して、 $M_{257}$ は15桁です。ここに壁があります。

今日の素因数分解アルゴリズムとしては一般数体篩法が最高効率です。類似のアイデアに基づくもう少し簡単なアルゴリズムには二次篩法もあります。これらの実装によるメルセンヌ数の素因数分解は今後の課題です。

まとめ

メルセンヌ数の歴史と基本的な性質を概観し、その素数判定法であるリュカ・レーマーテストを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]はまだ $M_{127}$ が既知最大の素数だった時代に書かれた、数論の歴史をまとめた本です。

ウェブサイト

メルセンヌ数に関する最新のデータや歴史はこれらが詳しいです。本記事の数学的な内容や歴史はほとんどここにも載っています。

脚注

アーキテクチャ依存の usize, isize はここでは除外 ↩︎
Integer Types - The Rust Programming Language ↩︎
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^{2\cdot 4}-1=255$ は $2^2-1=3$ と $2^4-1=15$ で割り切れますが、その積 $3\cdot 15=45$ では割り切れません。 $2^n-1$ の「因数分解」による整数の積への分解は、少なくとも円分多項式までは行えます。「円分多項式に2を代入した値」の数列がA019320として登録されています。合成数に対する「円分多項式に2を代入した値」が合成数となる最小の例は18番目の57です。 ↩︎
完全数 - Wikipedia。A000396 ↩︎
Édouard Luca(1842-1891) https://ja.wikipedia.org/wiki/エドゥアール・リュカ ↩︎
$M_{127}$ より小さな $M_{61}, M_{89}, M_{107}$ が素数であることが示されたのはその間の期間。Table of Known Mersenne Primes ↩︎
本気でやるならカラツバ法を使うのはいいアイデアだと思います。実演してみませんか？ ↩︎
num-bigintの"Alternatives"の項の表より。また、やや古いがRedditの投稿。 ↩︎
内部実装について—Wolfram言語ドキュメント ↩︎
筆者は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