前 6‒05 p ｜ 
もくじ i ｜ 次 6‒07 n 
JavaScript 用に最小構成的な「任意精度整数演算」ライブラリを作ってみた（_26_.js）。実験目的は「2進ベースの計算のパフォーマンス」の確認。「コンピュータだもん、10進より2進ベースの方が高速に決まってる」と思って始めたが、その考えは単純過ぎた。
内部的に約50万桁の計算が可能で、入出力は約3万桁まで対応できる。100桁や1000桁程度の数は軽快に扱えるし、100行のライブラリにしては、そう悪くもない。しかし2進ベースにしたことが逆効果となってしまい、内容的には不満の残るものとなった。
実用上の目的にはそれなりに役立つので「普通の四則計算・メルセンヌ予想の検証・平方根を100桁・円周率を1000桁」などを題材に、インターフェースを紹介したい（といっても四則演算くらいしかなくて、単純だが）。最後の節ではライブラリの構造と問題点について簡単に報告する。
JavaScript の数値型において整数の精度が保証される限界は、絶対値 253 − 1。例えば 253 + 6 = 9007兆1992億5474万0998 について、JavaScript をそのまま使うと正確な計算ができない:
var a = 9007199254740998; Append( a + " + 3 = " + (a + 3) ); Append( a + " - 3 = " + (a - 3) ); Append( a + " * 3 = " + (a * 3) );
よく見ると計算がおかしい。
このように、有効数字が15～16桁では足りない場合、任意精度演算が役立つ。
var a = new TinyBigInt("9007199254740998");
Append( a + " + 3 = " + a.add(3) );
Append( a + " - 3 = " + a.sub(3) );
Append( a + " * 3 = " + a.mul(3) );ライブラリの内容は単純明快で、new TinyBigInt でオブジェクトを作り、add、sub、mul、div でそれぞれ加減乗除を行うことができる。演算結果は新規オブジェクトとして返され、演算に使われたオブジェクトそのものは変化しない。
ライブラリとのデータのやり取り（例えば new するときの引数）では、基本的には文字列型を使う必要がある。絶対値 253 未満の小さな数については数値型でも文字列型でも構わないが、その範囲外の場合、数値型を使うと型エラー例外が発生する。ライブラリによって生成されたオブジェクトを再びライブラリ関数の引数としても構わない。
この記事では、ライブラリからの出力を Append という関数を使って表記する。console.log を使いたければ、それに置き換えてもよく、または
Append = console.log.bind(console); と定義してもいい。console.log を使いたくない場合、次の定義を使えば、任意の HTML ページにソースを埋め込んでテストすることができる（付録参照）。
function Append(str) {
var node = document.createElement("pre");
node.appendChild( document.createTextNode(str) );
document.body.appendChild(node);
}
メルセンヌ数 M67 = 267 − 1 = 147573952589676412927 は合成数。193707721 で割り切れる。メルセンヌ自身は、M67 が素数であると予想した。「こんな大きな数、どうせ誰にも確認できない」と、ハッタリの予想をしたのかもしれない。
原典を確認してみると、メルセンヌの主張は直接的には素数についてではなく、そこから生じる完全数についてのものだったようだ。「M67 から9番目の完全数が作れる」と予想したらしい。
ライブラリで 267 を使いたい場合、次のどちらの書き方でも構わない。
new TinyBigInt(2).pow(67); TinyBigInt.pow(2, 67);
2行目の書き方の方が簡潔だろう。それマイナス1ということなので、次のようになる。
var M67 = TinyBigInt.pow(2, 67).sub(1); Append( "M67 = " + M67 ); Append( "M67 / 193707721 = " + M67.div(193707721) ); Append( "M67 % 193707721 = " + M67.mod(193707721) );
div は整数演算であり、小数点以下を切り捨てた整数商を返す。割り切れたかどうかは mod の結果（割り算の余り）がゼロになるかどうかで判定される。
このくらいの計算 JavaScript 自身でも、できるのでは？ という気もするが…
var M67 = Math.pow(2, 67) - 1; Append( "M67 = " + M67 ); Append( "M67 / 193707721 = " + M67 / 193707721 ); Append( "M67 % 193707721 = " + M67 % 193707721 );
精度不足で M67 の末尾4桁が捨てられてしまい、割り切れるものも割り切れない。
メルセンヌは M257 も素数だと予想した。
var M257 = TinyBigInt.pow(2, 257).sub(1); // 2^257-1 Append( "M257 = " + M257 ); Append( "M257 = " + M257.toString() ); // explicitly Append( "M257 = " + M257.toString(0) ); Append( "M257 = " + M257.toString(-1) );
デフォルトの出力では読みにくい場合、toString() を明示的に呼び出して、引数 0 または -1 を渡すことができる。前者は、10桁ずつ区切って出力を行う（さらに100桁ごとに改行、1000桁ごとに空行が入る）。後者は、桁数が多い場合に、数字の頭と桁数だけを出力する。引数なしの toString() は、デフォルトの（暗黙の）変換と同じ動作となる。
2016年8月7日追記: ライブラリのバージョン0.2以下では、MSIE で toString(-1) が動作しなかった。この問題はバージョン0.3で修正された。
上記に続けて次の2行を実行してみる。
Append( M257.mod("535006138814359").compare(0) === 0 );
Append( M257.mod("1155685395246619182673033").compare(0) === 0 );compare(x) は、x と大小を比較して、大きい・等しい・小さい場合に、それぞれ（通常の数値型の）1, 0, −1 を返す。上記 compare(0) === 0 は「余りが0に等しいか」という判定に当たる。出力はいずれも true であり、要するに割り切れている。
上記の比較を次のように書くこともできる。
M257.mod("535006138814359") == 0;
この場合、左辺値が暗黙にプリミティブな型に変換され（具体的には文字列型に変換され）、それと 0 が比較される。「余りが0か」を調べるエレガントで可読性の高い書き方のようにも思える。しかし、オブジェクトが表している数が例えば1万桁だとどうなるか？ 値が 0 か問い合わせるたびに、いちいち長さ1万の文字列に変換されたらたまらない。だから面倒でも compare を呼び出す方がいい。
もっとも上記の例は、535兆0061億… で割った余りなのだから桁数も高が知れており、どちらでも構わないだろう。
一般に、ライブラリが生成したオブジェクトを普通の数値型のように扱って、通常の数値型の値と等号や不等号で比較することができるが、暗黙の変換コストが無駄に大きくなる場合もある（本物の数値型ではないので === は使えない）。ライブラリが生成したオブジェクトをそのまま JavaScript の数学関数に渡すこともできる。
もちろんオブジェクト同士については、等号や不等号で値を比較することはできない。
商と余りを同時に知りたい場合、次のように第2引数を true にして div を呼び出すと効率が良い（商・余りは内部的にはセットで計算されている）。この場合、最初の要素が商、2番目の要素が余りである配列が返される。
var arr = M257.div("535006138814359", true);
Append( "商 = " + arr[0].toString(0) + "\n余り = " + arr[1] );M257 は、次のように3因数の積に等しい。
Append( new TinyBigInt("53500 6138814359")
.mul("11556 8539524661 9182673033")
.mul("374550598 5018109365 8177663009 6313181393")
.compare(M257) );
ライブラリに渡す数値の文字列には、適当に空白類を挿入しても構わない。
夜中に突然 8 の平方根を100桁、知りたくなって困る。その辺の電卓では、桁数が足りない…。よくあることですね！
そこでご紹介するのが、この逐次近似。平方根を求めたい数 a と、a の平方根の近似値 b があれば、あとは、「a/b と b の平均をあらためて b と置く」を繰り返すだけで、平方根が1000桁でも2000桁でも作れちゃうんです。
「整数ライブラリで平方根？ 小数点以下はどうするの？」
「8 の平方根の代わりに 8×10200 の平方根を求めれば結果は √8×10100 だろ？ ざっと100桁の整数じゃないか」
「なるほどね。でも出発点となる近似値はどうすれば…」
「JavaScript だもん、そんなのは Math.sqrt で一発さ！」
訳の分からない通販番組のノリで、ハイテンションに10回ループしてみましょう。
var a = new TinyBigInt( "8e200" ); // 8 * 10^200
var b = new TinyBigInt( Math.sqrt(8) + "e100" ); // √8 * 10^100 = 2.8284271247461903e100
for( var i = 0; i < 10; ++i ) {
Append( "i = " + i + ": b = " + b );
b = a.div( b ).add( b ).div( 2 ); // (a/b + b)/2
}10回繰り返すまでもなく、「平均値を求める処理」の3回目で既に100桁確定している。
実は、この方法では1回ループするたびに有効桁数がほぼ2倍になる。出発点となる Math.sqrt は、普通は15～16桁有効なので、計算したい桁数が15、30、60、120…桁以内のとき、それぞれ0、1、2、3…回ループすれば足りる。つまり、b の桁数 n について、log2(n/15) 回以上ループすれば全桁有効になる。従って:
function Sqrt10e( x, d ) {
if( isNaN( x ) || isNaN( d ) || x % 1 || d % 1 ) throw new TypeError();
if( x < 0 || x >= 1e15 || d < 0 || d > 5000 ) throw new RangeError();
var a = new TinyBigInt( Number( x ) + "e" + 2*d );
var b = new TinyBigInt( Math.sqrt( x ) + "e" + d );
var n = d + 10; // number of digits of b (overestimated)
var required = Math.log( n / 15 ) / Math.LN2; // Math.log2( n / 15 )
for( var i = 0; i < required; ++i ) {
b = a.div( b ).add( b ).div( 2 );
}
return b;
}
この関数は、15桁以内の任意の（負でない）整数 x と、5000以下の整数 d に対して、√x × 10d の整数部分を返す。実質的に、√x を小数点以下 d 桁（最大5000桁）求めることに当たる。最初の2行は、引数の値・範囲に問題がないかチェックするもので、本筋と関係ない。残りは 例4.1 とほぼ同じで、変数 a と b を指数表記で設定している。
実用的な時間で計算ができる上限桁数は、環境によって異なる。試す場合、いきなり極端な値を入力せず、少しずつ値を増やしながら実験しよう。
JavaScript 自身は 1e+308 くらいの数までしか対応できないが、このライブラリは 9.99e+1000 でも 8e+2000 でも対応できる。ただし、これらの引数は、"8e2000" のように、文字列として扱う必要がある（e の後ろの + は省略可）。仮数部 sig が整数を表している場合、指数部を exp として、
new TinyBigInt(sig + "e" + exp)
よりも
new TinyBigInt(sig).mul(TinyBigInt.pow(10, exp))
または
TinyBigInt.pow(10, exp).mul(sig)
の形を使った方が効率が良いだろうが、どの書き方でも動作する。
ループの必要回数 required は、b の桁数 n に依存する: b の仮数部は「15桁以内の x」の平方根なので8桁以内、指数部は d なのだから、n = d + 10 と見積もれば十分。ES6 より前の JavaScript には Math.log2 がないが、自然対数を Math.LN2 で割れば代用になる。必要回数がゼロ以下になる場合、1度もループする必要がない。
x の大きさに制限のない一般の場合、この逐次近似は、まれに正しい値より 1 大きい値を返すことがある。気になる場合、戻り値を
return b.mul(b).compare(a) <= 0 ? b : b.sub(1) ;
とすれば全桁正確にできる。
厳密に言うと、Math.sqrt の精度は実装依存であり保証されていない。その点まで考慮するなら、ループ部分を次のようにすればいい（速度は低下する）。
for( var i = 0; i < 20; ++i ) {
var old_b = b;
b = a.div( b ).add( b ).div( 2 );
var diff = old_b.sub( b );
if( diff.compare( 1 ) <= 0 && diff.compare( -1 ) >= 0 ) break;
}
最後の行は
if( -1 <= diff && diff <= 1 ) break;
と書いた方が可読性が高いが、[3/6] で記したように、効率の点で問題がある。
diff がゼロになったら（前回の近似値と同じになったら）ループ終了、という実装では問題が生じる。というのは、この逐次近似は必ずしも収束せず、まれに1違いの2個の整数の間を振動するパターンに陥ることがある。
ラマヌジャンによる円周率の近似式のうち、次のものがよく引用される:
![1/π = (√8/9801) ∑ { [(4n)!(1103+26390n)] / [(n!)^4 × 396^(4n)] }; n=0,1,2…](/image/2016/ramanujan78.png)
総和記号の右側の分子を a、分母を b とすると:
var a = TinyBigInt.factorial( 4*n ).mul( 1103 + 26390*n ); var b = TinyBigInt.factorial( n ).pow( 4 ).mul( TinyBigInt.pow( 396, 4*n ) );
ここで TinyBigInt.factorial(x) は、x の階乗を表す。
どんな分数になっているのか、最初の方を出力してみる。
n = 0: a = 1103, b = 1 n = 1: a = 659832, b = 24591257856 n = 2: a = 2172562560, b = 9675679407044507467776 n = 3: a = 38450895436800, b = 19272907305680273946696375374708736 n = 4: a = 2231687537823744000, b = 121329928496325076509677917742087973855188484096 n = 5: a = 323704910893926481920000, b = 1864784723314482690530261346690526668068454009929056911360000
総和記号の部分は a を b で割って足し合わせるだけだが、使っているのが整数ライブラリなので、このままでは、最初の項以外、a/b はゼロになってしまう。円周率を d 桁計算するとして、分子を 10d 倍してから b で割ればいい。D = 10d としておこう。無限級数なので、どこまで足して打ち切るかという問題がある。単純に考えて、商がゼロになったら打ち切ることにしよう。
試しに d = 100 としてやってみると、次のように 9801/(π√8) = 1103.00002683197457346381340888… に向かって収束していくことが分かる。
var d = 100;
var D = TinyBigInt.pow( 10, d );
var sum = new TinyBigInt( 0 );
for( var n = 0; n < 1000; ++n ) {
var a = TinyBigInt.factorial( 4*n ).mul( 1103 + 26390*n );
var b = TinyBigInt.factorial( n ).pow( 4 ).mul( TinyBigInt.pow( 396, 4*n ) );
var c = a.mul( D ).div( b );
sum = sum.add( c );
/// Append( "n = " + n + ": sum = " + sum );
if( c.compare(0) === 0 ) break;
}Append を有効にした場合）この整数に D√8 を掛けて9801で割ると、円周率の逆数に当たる D2/π となる。「D3 ÷ その数」を計算すれば、結果は πD = π10d であり、円周率が d 桁得られる。
要するに、例5.1 の末尾（ループを抜けた後）に次の3行を追加すればいい。例4.2 の Sqrt10e を再利用した。
var inverse_of_pi = sum.mul( Sqrt10e( 8, d ) ).div( 9801 ); var pi = TinyBigInt.pow( D, 3 ).div( inverse_of_pi ); Append( "pi = " + pi.toString(0) );
末尾の真の値は 70680 ではなく 70679。末尾には誤差が入ることがあるので、10桁くらい多めに計算して末尾を捨てるといい。
せっかくなので、d = 1010 として、1000桁くらい計算してみよう。最適化されていない原始的なコードだが、テスト環境では0.2秒台で計算終了した。
動作速度は環境によって異なる。試す場合、少しずつ値を増やしながら実験しよう。
実際には項ごとに毎回階乗やべき乗を1から計算するのは無駄なので、直前の項の構成要素に計算を継ぎ足して次の項を作る方がいい。ラマヌジャン自身が書いたもともとの式も、そうなっている:
![1/(2π√2) = (1103/99^2) + (27493/99^6)(1/2)[(1⋅3)/4^2] + (53883/99^10)[(1⋅3)/(2⋅4)][(1⋅3⋅5⋅7)/(4^2⋅8^2)] + …](/image/2016/ramanujan44.png)
var a1 = new TinyBigInt(1103), a2 = new TinyBigInt(1), a3 = new TinyBigInt(1);
var b1 = new TinyBigInt(9801), b2 = new TinyBigInt(1), b3 = new TinyBigInt(1);
var n = 0;
for(;;) {
var a = a1.mul( a2 ).mul( a3 );
var b = b1.mul( b2 ).mul( b3 );
var c = a.mul( D ).div( b );
if( c.compare(0) === 0 ) break;
sum = sum.add( c );
++n;
a1 = a1.add( 26390 );
a2 = a2.mul( 2*n-1 );
a3 = a3.mul( (4*n-3)*(4*n-1) );
b1 = b1.mul( 96059601 );
b2 = b2.mul( 2*n );
b3 = b3.mul( 16*n*n );
}
var inverse_of_pi = sum.mul( Sqrt10e( 8, d ) );
var pi = TinyBigInt.pow( D, 3 ).div( inverse_of_pi );本当は、スケール D を10のべきではなく、ライブラリのアーキテクチャーに合わせた大きさにして、掛け算の代わりにオブジェクト内部の配列をシフトした方が少し効率がいい。しかし、多少工夫しても、このライブラリでは平方根・円周率とも5000～6000桁くらいが実用上の限界のようだ。
2019年6月23日追記: 例4.2 の Sqrt10e は d > 5000 を受け付けない。5001桁以上を計算したい場合、|| d > 5000 の部分を削除して、この制限を撤廃する必要がある。例5.1 では、ループが最大1000回に制限されている。円周率7900桁くらいまではこれで問題ないが、それを超える場合、ループ回数を増やしてもっと先まで足し算する必要がある。ループ回数の制限をなくすか、または、もともと制限のない例5.3のコード（例5.1 より高速でもある）を使うことができる。
ちゃんとした任意精度ライブラリは既にいろいろある。JavaScript 版はあまりないとしても、言語によっては標準でサポートされている。筆者は、PARI/GP をよく使う。それでも自分で一から組んでみると、いろいろな意味で勉強になった。
ライブラリのコアは98行。そのうち公開用のインターフェースを提供する部分が39行。
40行目から98行目までが実際に加減乗除などを行う実装部。
この他、おまけとして、べき乗と階乗の関数、指数表記された数値文字列を解釈するためのヘルパー関数、出力をフォーマットするためのヘルパー関数などが同梱されている。これらは、無くても動作する（例えば、べき乗は mul を繰り返し呼び出せば同じことができる）。
2016年6月5日追記: バージョン0.2では新機能が追加され、若干行数が変わった。
オブジェクトのプロパティーは、$ という配列と、is_negative というブール値の2個だけ。後者は負なら true、負でなければ false となる。前者には「226進数」で表現された整数の絶対値の各「桁」が、小さい方を先頭に格納されている。例えば、1015 という整数は
13008896×(226)0 + 14901161×(226)1 なので、「226進数」では [13008896, 14901161] となる。
この部分を「232進数」（32ビットベース）にする手もある。比較すると:
ではなぜ26ビットを選択したのか。それは:
他方、10のべきをベースにする方法もある。コンピューター的には2のべきを使った方がシフト演算などが活用できるので有利だが、そうやって計算効率が高くなっても、必ずしもトータルでは速くならない。例えば、ユーザーが1万桁の整数を2個入力して、それを掛け算して答えを出力するとする。そのとき、計算そのものは2のべきベースの方が有利だとしても、「ユーザーが入力した10進数を内部の2進表現に変換する手間」と「計算結果の内部表現を10進数に再変換する手間」を考慮すると…。これらの変換の手間を上回るほど計算が速くならない限り、2進ベースのメリットが生かされない。
2進数と10進数の相互変換はコンピューターの世界では非常に基本的なので、きっと洗練された速い方法があるに違いない。しかしこの小さなライブラリにおいては、その部分が大きなボトルネックとなった。「2進・10進の変換の手間を考えても、トータルでは有利」どころか、むしろ「2進・10進の変換の手間」が計算自体のコストより大きい。これでは本末転倒だ。コードの行数の上でも、演算の実装約60行のうち、約3分の1が2進・10進変換のために使われていて、無駄な感じがする。
速度が低下する典型的な原因は、空間的なものらしい。配列が長くなり過ぎると、ある時点から効率が極端に低下するようだ。それに加えて、2進・10進変換で使われる「割り算」もあまり速くない。
足し算・引き算・掛け算については、ベースが何であれ、基本的なやり方は明らかだろう。「繰り上がり」などを処理しながら、「1桁」ずつ計算するだけだ。
問題は割り算。「2進数の割り算」の基本は1ビットを1桁とする計算であり、筆算の割り算と同様に「上から下が引けるなら商を立てて引き算する」。2進数なので、商が立つときは必ず「1」で、単純に「引けるか引けないか」だけが問題になる。 DIV and MOD: Division and modulo を参考にした。
「ソフトウェアレベルで1ビットずつ割り算なんて…。1回の割り算が何万ステップになるのやら」と心配になるが、やってみると思ったよりは速く、一応実用になることが分かった。
割り算の最初に、分子側の数は、最上位ビットが立つ位置まで左にシフトされる。シフトの大きさを決めるために clz が必要になるが、これについては Find the log base 2 of an integer with a lookup table が参考になり、刺激的でもあった。しかし最初のページにあるような「while ループを使って端から一つずつ見ていく原始的な方法」も実用的だ。分子がランダムに与えられるとすると1/2の確率で最上位ビットはもともと立っているし、3/4の確率でシフト量は1以内、7/8の確率でシフト量は2以内…などとなって、大半のケースでは while ループはすぐ終了し、決して遅くない。ES6で導入された Math.clz32 については、テスト環境において、速度的なメリットが分からなかった。割り算が全体的に遅いため、clz のわずかな速度差はどうでもいいのだろう。
10進ベースの割り算と比べると、2進ベースの割り算は美しい。「シンプルなアルゴリズム一つで、何万ビットの割り算でも対応できる」という点は、気持ちがいい。しかし現実的には、簡潔だから良いとは限らない。汚くても浮動小数点数でガンガン計算する方が得、ということもありうる。
以前10進ベースで任意精度の実装をしたことがあり、その当時から「2進ベースにしたら速くなるだろうな」と思っていた。そこで今回は「2進ベース」の実装を試した。結論としては、2進ベースは（メリットもあるが）思ったほど良くなかった。
特にショックだったのが「1万の階乗」（約3万6000桁）の計算に数秒かかり、書き出しに10秒以上かかること。階乗に特化したアルゴリズムではないし、小さい桁数に最適化されているので直接比較はできないが、昔作った10進ベースのアルゴリズムでは「1万の階乗」は1秒未満だったような…。
入出力の頻度にもよるのだろうが、JavaScript 用のライブラリでは、結局10進ベースが正解かもしれない。
2016年6月5日追記: その後、やはり2進数には強みがあることが分かった。「[JS] メルセンヌ数の分類と分解」参照。
<html>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<head>
<title>Test</title>
<script type="text/javascript" charset="utf-8" src="https://yosei.fi/demo/TinyBigInt/_26_.js"></script>
</head>
<body>
<script type="text/javascript">
//<![CDATA[
try {
var a = TinyBigInt.pow( 2, 10000 );
Append( "2^10000 = " + a.toString(0) );
} catch(e) {
Append(e);
}
//]]>
</script>
</body>
ネットワーク越しに使うと転送の遅延があるので、できれば、_26_.js をローカルに置いて、script タグの src の値をそれに合わせた方が軽快に使える。
[6/6] に記したように、このライブラリには「2進数と10進数の変換がボトルネックになっている」という欠陥がある。しかし実用上、100～1000桁くらいまでの「比較的小さな巨大整数」を少々扱う分には、特に問題もないようだ。パブリックドメインなので、使い道があったら自由に使ってほしい。
n は、値の絶対値が 253 以上の場合、整数として正確であるとは限らない。
それを承知の上で（例えば近似計算において）大きな n をライブラリ関数に渡したい場合には、文字列型に変換してから渡す必要がある。その方法はString(n)"" + nn.toFixed()TypeError 例外が出る。
ユーザーが意識していない部分で精度が失われることを予防するため、意図的にこのようにしてある。絶対値が 253 未満なら数値型でも構わない。"1.23e1000" は "123e998" と同じ）。2016年6月5日追記: 6月4日にバージョン0.2を公開した。上記の最初のバージョンとの主な違いは、16進数の文字列の入出力に対応したことと、2進数の文字列の出力に対応したこと。
2016年8月7日追記: バージョン0.3を公開した。変更点は次の通り。(1) MSIE で toString(-1) が動作しない不具合を修正。1カ所でMSIE未サポートの Math.log10 が使われていたのが原因。 (2) 平方を計算する関数 sqr() を追加した。 (3) 四則演算を多少高速化した。
2016年9月11日追記: バージョン0.4を公開した。new TinyBigInt(true) が 0 を表すオブジェクトを生成していた不具合を修正し、1 を表すオブジェクトを生成するようにした。
メルセンヌ数を素数と合成数に分類する。試行除算と「P − 1」法を使い、メルセンヌ合成数の分解を試みる。1万桁以内のメルセンヌ素数・全26種をリストアップする。
数千万桁のメルセンヌ素数が脚光を浴びるが、その裏では、たった数百桁のメルセンヌ合成数が分解できない。
この記事全体を通じて、p を3以上の素数とし、メルセンヌ数 Mp = 2p − 1 について考える。
次の方法で、メルセンヌ数が素数か合成数か判定できる。
[1] : 4 × 4 − 2 = 14 [2] : 14 × 14 − 2 = 194 [3] : 194 × 194 − 2 = 37634 [4] : 37634 × 37634 − 2 = 1416317954 [5] : 1416317954 × 1416317954 − 2 = 2005956546822746114
最後の数 200京5956兆5468億2274万6114 は 127 で割り切れるので、メルセンヌ数127 は素数と分かる。127 で割った余りだけが分かればいいのだから、計算の途中で現れる数は、全て 127 で割った余りに置き換えて構わない:
[1] : 4 × 4 − 2 = 14 : 14 % 127 = 14
[2] : 14 × 14 − 2 = 194 : 194 % 127 = 67
[3] : 67 × 67 − 2 = 4487 : 4487 % 127 = 42
[4] : 42 × 42 − 2 = 1762 : 1762 % 127 = 111
[5] : 111 × 111 − 2 = 12319 : 12319 % 127 = 0最後の剰余（余り）が 0 なので、素数と判定される。
[1] : 4 × 4 − 2 = 14 : 14 % 2047 = 14
[2] : 14 × 14 − 2 = 194 : 194 % 2047 = 194
[3] : 194 × 194 − 2 = 37634 : 37634 % 2047 = 788
[4] : 788 × 788 − 2 = 620942 : 620942 % 2047 = 701
[5] : 701 × 701 − 2 = 491399 : 491399 % 2047 = 119
[6] : 119 × 119 − 2 = 14159 : 14159 % 2047 = 1877
[7] : 1877 × 1877 − 2 = 3523127 : 3523127 % 2047 = 240
[8] : 240 × 240 − 2 = 57598 : 57598 % 2047 = 282
[9] : 282 × 282 − 2 = 79522 : 79522 % 2047 = 1736最後が 0 にならないので、このメルセンヌ数は合成数と判定される。実際、2047 = 23 × 89。
JavaScript 上で任意精度ライブラリを使い、単純発想で実装してみる:
function Simple_Lucas_Lehmer( Mp, p ) {
var s = new TinyBigInt( 4 );
for( var i = 1; i <= p - 2; ++i )
s = s.mul( s ).sub( 2 ).mod( Mp ); // s := s*s − 2 (mod Mp)
return ( s.compare(0) === 0 );
}この関数は、問題のメルセンヌ数が素数なら true、合成数なら false を返す。
素数の配列 PrimeTable を用意して次のように順々にテストすると、小さいメルセンヌ数が素数か合成数か簡単に判別することができる。
for( var i = 0; i < PrimeTable.length; ++i )
{
var p = PrimeTable[ i ];
if( p > 200 ) break;
var Mp = TinyBigInt.pow( 2, p ).sub( 1 ); // Mp = 2^p − 1, e.g. M3 = 2^3 − 1
var message = "M" + p + " = " + Mp + " is "; // e.g. "M3 = 7 is "
if( p === 2 || Simple_Lucas_Lehmer( Mp, p ) )
Append( message + "prime" );
else
Append( message + "composite" );
}このテストは p = 2 に対しては有効ではない（M2 = 3 は素数だが「素数でない」と判定されてしまう）。この場合については、手動で「素数」と分類している。
Simple_Lucas_Lehmer は、p < 3000 くらいの範囲のメルセンヌ数に対して一応実用になる。テスト環境では、918桁の M3049 が合成数であることを判定するのに約6秒かかった。遅さの主因は剰余演算 mod （本質的には割り算）。掛け算 mul についてもライブラリ側に高速化の余地があるが、割り算に比べると影響はざっと30分の1にすぎず、優先順位は低い。
[5/7] では、実装から割り算をなくして効率を改善し、約25倍の高速化を行う。[6/7] では、それを使って、1桁上の M20000 くらいに挑戦したい。
ところで上記 M3049 は合成数と判定されるものの、具体的にどういう因数を持つのか分かっていない（2016年現在）。
素数の配列を作るには、まず最初の素数 2 を含む長さ1の配列を作り、3 以上の奇数について、それが奇素数であれば配列の末尾に追加すればいい。ここでは 216 未満の素数をデータ化しておくことにする。奇素数であるかの判定は、自分の平方根以下のどれかの奇素数で割り切れれば false、割り切れなければ true。
var PrimeTable = [ 2 ];
for( var odd = 3; odd < 65536; odd += 2 ) {
if( Is_odd_prime( odd ) ) PrimeTable.push( odd );
}
function Is_odd_prime( n ) {
var sqrt = Math.sqrt( n );
/// if( PrimeTable[ PrimeTable.length - 1 ] < sqrt ) throw new RangeError( n ); // *1
for( var i = 1, q = 3; q <= sqrt; q = PrimeTable[ ++i ] ) {
if( n % q === 0 ) return false;
}
return true;
}これで 216 までの素数表ができたので、今後 Is_odd_prime を呼べば、232（約42.9億）くらいまでの数は正しく素数判定できる。
巨大なメルセンヌ数をこの関数に渡して、直接的に素数性を判定させることはできない。「間違って大き過ぎる数を渡したら、エラーになってほしい」という場合は、*1 の行を有効にすればいい。
Lucas–Lehmer テストで合成数と判定された場合、「合成数」という事実だけが与えられる。「合成数だとすれば、具体的に何で割り切れるのか」という点は、別に考える必要がある。
因数を探すのは、試行錯誤によるしかない。それでも多くのショートカットや巧妙なアプローチが存在して、安楽椅子探偵が推理だけで犯人を追い詰めるようなスリルと面白さがある。
でも、そんなのんきなことを言っていられるのは最初のうちだけ。桁数が増えるにつれ、問題は急激に難しくなる。
Mp = 2p − 1 が未知の素数 x により割り切れた（“殺害された”）とする。このとき 2p ≡ 1 (mod x) だが、この p は mod x での 2 の位数（2 の肩に乗せて ≡ 1 (mod x) となる最小の正の数）。
なぜなら 2p ≡ 1 (mod x) を満たす p は 2 の位数の倍数に限られ、言い換えれば 2 の位数は p の約数。約数といっても p は素数なのだから 1 or p 以外の選択肢はない――そのうち、位数 1 という選択肢は不合理（21 ≡ 1 になるわけない）; 消去法で位数は p。
位数が p という事実とフェルマーの小定理 2x − 1 ≡ 1 (mod x) を組み合わせると、x − 1 は p の倍数。従って、ある自然数 a を使って x = ap + 1 と書ける。これは犯人 x についての重要な手掛かりだ。
例えば M67 は 1垓4757京3952兆5896億7641万2927 だが、その最小因数 x = 1億9370万7721 はざっと2億なので、仮に3以上の奇数で順に M67 を割ってそれを見つけようとすると約1億回の試し割りが必要になる。21桁という M67 の大きさを考えると、「割り算1億回」は良いニュースではない。それが原因ではないだろうが、元祖メルセンヌも「M67 は素数」という錯覚に陥ってしまった。
けれど良いニュースがある。上記の考察によれば x = 67a + 1 という形になる。しかもこの a は、偶数でなければならない（奇数だとすると、x が偶数になってしまう）。a = 2k として、x = 2 × 67k + 1 の形の数だけ調べればいい。全部の奇数を試す場合と比べて候補が67分の1に減り、「1億回」から「150万回」まで計算量が激減する。
一般の p, x に対しても、ある自然数 k が存在して x = 2kp + 1 を満たす。従って (x − 1)/2 = kp。 2の肩に乗せると、2(x − 1)/2 = 2kp だが、この左辺は「2 が平方剰余か」についてのオイラーの基準の形になっている。一方、右辺は ≡ 1 (mod x)。なぜなら 2p ≡ 1 (mod x)。オイラーの基準によって「2 は法 x の平方剰余」と分かり、第二補充法則によって x ≡ 1 or 7 (mod 8) と分かる。
最初の結果と合わせると、x の候補は整数 8p 個につき2個に絞られる（以下のコードでは n1 と n2 と呼ばれる）。M67 で言えば割り算回数がさらに50%割引となり「75万回」以内で済む。
「私の推理によれば、M67 を割った犯人は 67a + 1 の形で、a は 8 の倍数か、それに 2 を足したものだろう」
「それはどうして？」
「初等的だよ、ワトソン君。mod 8 において x ≡ 1 なら、明らかに 2kp は 8 の倍数。x ≡ 7 なら 2kp ≡ 6 となり、p = 67 ≡ 3 だから 2k ≡ 2 じゃないか」
初等整数論探偵は、得意げにつぶやくのだった。
この推理を基に a = 8, 16, 24, 32, … と a = 2, 10, 18, 26, 34, … を順に当たると、果たして a = 361395 × 8 のとき M67 が割り切れることが判明した。犯人が見つかったのである。
ここまでを実装してみよう。
function Find_Factor( Mp, p ) {
var t = ( p % 4 === 3 ) ? 2 : 6 ;
var n1 = t * p + 1; // n1 ≡ −1 mod 8
var n2 = 8 * p + 1; // n2 ≡ +1 mod 8
var n_max = 2e+8; // 200,000,000
for( var inc = 8 * p; n1 <= n_max; n1 += inc, n2 += inc ) {
if( Mp.mod( n1 ).compare( 0 ) === 0 ) return n1;
if( Mp.mod( n2 ).compare( 0 ) === 0 ) return n2;
}
return "?";
}この関数は、2億以下の因数を探すように設定されている。検索範囲内では約数が見つからない場合、"?" が返される。
n1 と n2 は、前述の x = ap + 1 に当たる。どちらも ≡ 1 (mod p) であり、かつ、前者は ≡ 7 ≡ −1 (mod 8)、後者は ≡ 1 (mod 8) を満たす（そうなるように、t が選択されている）。a の二つの初期値は、p ≡ 1 (mod 4) なら 6 と 8、p ≡ 3 (mod 4) なら 2 と 8。Mp を n1 と n2 で割ってみて、どちらかで割り切れれば、それが求める因数。割り切れなければ、それぞれ 8*p を足して再試行する。
テスト環境において、M67 の因数 193707721 は2秒未満で発見され、割り算は約72万回だった。
上記の検索には無駄がある。n1 または n2 がかなりの割合で合成数（例えば3の倍数や5の倍数）になっているが、検出されるのは最小因数なのだから、それは必然的に素因数であり、合成数で割り算しても仕方がない（どうせ割り切れない）。割り切れないことが分かっているケースについては、遅い mod の呼び出しをなるべく回避したい。
例えば、if の部分を
if( n1 % 3 && Mp.mod( n1 ).compare( 0 ) === 0 )
のようにすれば、3の倍数による割り算が省略され、割り算回数が3分の2になり、それだけでも30%の高速化が見込まれる。同じことを5の倍数、7の倍数…についてもやるべきだろう。
問題は、それをどこまでやるかだ。
理論上、n1 が素数の場合だけ割り算を実行すれば、割り算回数は最小になる（n2 についても同様）。つまり:
function Find_Factor1( Mp, p ) {
var t = ( p % 4 === 3 ) ? 2 : 6 ;
var n1 = t * p + 1; // n1 ≡ −1 mod 8
var n2 = 8 * p + 1; // n2 ≡ +1 mod 8
var n_max = 2e+8; // 200,000,000
for( var inc = 8 * p; n1 <= n_max; n1 += inc, n2 += inc ) {
if( Is_odd_prime( n1 ) && Mp.mod( n1 ).compare( 0 ) === 0 ) return n1;
if( Is_odd_prime( n2 ) && Mp.mod( n2 ).compare( 0 ) === 0 ) return n2;
}
return "?";
}コード2.1 の Is_odd_prime で予選を行い、Mp.mod の呼び出し回数を節約しているが、Is_odd_prime の呼び出しも無料ではない。トータルでは得だろうか、損だろうか？
チェックを甘くして「小さな素数で割れなければ多分素数」とする手もある。例えば、500未満の素数（95番目の素数499まで）を使うと:
function Find_Factor2( Mp, p ) {
var t = ( p % 4 === 3 ) ? 2 : 6 ;
var n1 = t * p + 1; // n1 ≡ −1 mod 8
var n2 = 8 * p + 1; // n2 ≡ +1 mod 8
var n_max = 2e+8; // 200,000,000
for( var inc = 8 * p; n1 <= n_max; n1 += inc, n2 += inc ) {
if( Maybe_odd_prime( n1 ) && Mp.mod( n1 ).compare( 0 ) === 0 ) return n1;
if( Maybe_odd_prime( n2 ) && Mp.mod( n2 ).compare( 0 ) === 0 ) return n2;
}
return "?";
}
function Maybe_odd_prime( n ) {
if( n <= 249001 ) return Is_odd_prime( n ); // 499*499 = 249001
for( var i = 1; i <= 94; ++i ) { // Prime 499 has index 94
if( n % PrimeTable[ i ] === 0 ) return false;
}
return true;
}Find_Factor: 193707721 (1651.4 ms) 722790 calls to Mp.mod()
Find_Factor1: 193707721 (770.1 ms) 81544 calls to Mp.mod()
Find_Factor2: 193707721 (369.9 ms) 132276 calls to Mp.mod()改良案1も効果があるが、改良案2の方がさらに効果的であることが分かる。改良案2と比べた場合、改良案1では Mp.mod の呼び出し回数は少ないが、呼び出し回数を減らすためのコストがかさんでトータルでは不利になってしまった。
実行時間の計測結果はどれも一例であり、同じ環境でも毎回少し違う。同じコンピューター上でも、JavaScript が遅いブラウザだと何十倍も時間がかかる場合がある。
試行除算の上限が比較的小さい場合（特に約9490万以内の場合）、べき剰余を利用することで計算をさらに高速化することができる。付録ソースコードにある Find_Factor3 がその例で、次のように差は歴然としているが、以下では説明の都合上 Find_Factor2 を使う。
Find_Factor: 13821503 (124.8 ms) Find_Factor1: 13821503 (22.5 ms) Find_Factor2: 13821503 (24.3 ms) Find_Factor3: 13821503 (3.0 ms)
Find_Factor2 を使って試行除算を実施した。
for( var i = 0; i < PrimeTable.length; ++i )
{
var p = PrimeTable[ i ];
if( p > 200 ) break;
var Mp = TinyBigInt.pow( 2, p ).sub( 1 );
var message = "M" + p + " = " + Mp + " is ";
if( p === 2 || Simple_Lucas_Lehmer( Mp, p ) )
Append( message + "prime" );
else
Append( message + "divisible by " + Find_Factor2( Mp, p ) );
}得られた約数はあくまで最小の素因数。「この因数で割れば素因数分解完了」とは限らないが、今は「合成数ならとりあえず一つ素因数を見つける」ことを目標とする。
M101 以降では、? が出るケースがある。合成数と判定されたものの、試行除算では因数が見つからない場合だ。上記の設定（2億までの割り算）では、p < 200 の範囲において次の8個がこれに該当する:
M101 = 2 5353012004 5645880299 3406410751 (31 digits) M103 = 10 1412048018 2583521197 3625643007 (32 digits) M109 = 649 0371073168 5345356631 2041152511 (33 digits) M137 = 17 4224571863 5204932932 4779900506 5324265471 (42 digits) M139 = 69 6898287454 0819731729 9119602026 1297061887 (42 digits) M149 = 71362 3846352979 9405291429 8472474756 8191373311 (45 digits) M157 = 18268770 4666362864 7754606040 8953537745 6991567871 (48 digits) M199 = 8034690221 2949513777 0981046170 5813012611 0149689139 6417650687 (60 digits)
合成数のくせに割り切れないこの連中、どうしてくれよう。
M103, M109, M157 については、試行除算を30億まで進めれば因数が見つかり、上記のアルゴリズムでも一応対応できる。しかし、M101, M139, M199 はやや困難で試行除算では数分～数時間かかると見込まれる。さらに、M137 と M149 の2個は段違いに分解が難しく（推定所要時間・数千年！）、このアルゴリズムでは歯が立たない。もっと強力な新兵器が必要だ。
比較的小さい因数を持つ限り、巨大なメルセンヌ数でも、試行除算による分解は可能だ。例えば 3万0104桁の M100003 が因数 6400193 を持つことは、1秒台で検出された。余因数の表示込みで3秒半。PARI/GP で factor(2^100003-1,6400200) とすると同じ処理に12秒かかるので、試行除算も捨てたものではない（しかも JavaScript で PARI に勝つとは…）。ただし、この場合、数が巨大過ぎて Simple_Lucas_Lehmer は実用にならない。
Find_Factor2 を使うと、p < 200 のメルセンヌ合成数34個のうち26個は割り切れるが、上記のように8個が因数不明となる。試行除算の上限を30億に上げれば、この数は5個まで減るが、かなり時間がかかる（特に Find_Factor1 では時間がかかり過ぎる）。
30億といえば、まだ10桁。たった10桁の因数でも厳しいというのでは（そして因数の長さが1桁増えるごとに困難度が急増するというのでは）、科学の威信にかかわる。10桁の割り算で、もう限界というのでは…。
だがそれが現実だ。10桁が限度ではないけれど、現在の知見では、自由に素因数分解できる数の上限はかなり低い。
さいわい数学の進歩は続いている。1970年代以降、「試行除算では困難な場合でも、他にも道はある」ということが少しずつ分かってきた。
ここでは、英国の数学者 John Pollard が発見した P − 1（ピー・マイナス・ワン）法を試してみたい。
この分解法はメルセンヌ数専用ではなく、一定の条件さえ満たされれば、どんな数に対しても通用する。
メルセンヌ・ウィキの記事（P-1 Factorization Method - Mersennewiki）を参考にして、Step 1（ステージ1）の骨組みをコード化した（アルゴリズムの動作原理については「楕円曲線法の前夜」参照: そちらでは下記の aE が Jm と呼ばれている）。
function P_minus_1_test( N ) {
var B1 = 500;
var a_E = new TinyBigInt( 3 ); // a^E : initially a=3 itself (E=1)
for( var i = 0; i < PrimeTable.length; ++i )
{
var q = PrimeTable[ i ]; // q = 2, 3, 5, 7, 11, ...
if( q > B1 ) break;
var Eq = Math.floor( Math.log( B1 ) / Math.log( q ) ); // floor(log_q(B1)) : max integer Eq such that q^Eq <= B1
var q_Eq = Math.pow( q, Eq ); // q^Eq : small enough for JavaScript since <= B1
// Append( "q=" + q + " : Eq=" + Eq + " : q^Eq=" + q_Eq + " : q^(Eq+1)=" + q_Eq*q ); // *0
a_E = _powmod( a_E, q_Eq, N ); // a_E := (a_E)^(q_Eq) mod N
// Eventually a_E = a^E mod N
// = a^(2^E2 * 3^E3 * 5^E5 * ... * q^Eq) mod N
}
var g = _gcd( a_E.sub( 1 ), N ); // gcd( a^E-1, N )
return g;
}アルゴリズム名「P − 1」の「P」は、見つけようとしている約数を指し、メルセンヌ数 Mp の p とは意味が異なる。
ここでは、この約数を g と呼ぶ。例えば、M67（21桁の整数）の因数 193707721 が発見された場合、p = 67 で g = 193707721 となる。
計算に使われる上限値 B1（コード内では B1）については、暫定的に 500 とした。
Eq は、メルセンヌ・ウィキの E2, E3, E5, … に当たる。Eq の選択については *0 のデバッグ出力を有効にすると、意味が分かりやすい。「q の Eq 乗」がだいたい B1 に等しければいいのだが、ここでは「q の Eq 乗が B1 を超えないような最大の Eq」を選択している。要するに logq B1 の整数部分。
例: B1 = 500 において q = 2 を考える。28 = 256 は B1 以下だが、29 = 512 だと B1 を超えてしまう。この場合、Eq は 8 で q_Eq は 28 = 256。
_powmod(x, y, z) は xy (mod z) を返す。法 z がメルセンヌ数であることを利用した高速化を後から行うが、今のところ法がメルセンヌ数であるという事実は使われていない。つまり、上記のアルゴリズムは、そのままで一般の合成数に対して有効。
_gcd は2数の最大公約数を返す。最大公約数が 1 または N になってしまったら意味がないけれど、それ以外の値になれば、意味のある N の約数が得られたことになる。
コード4.1と4.3には、メルセンヌ・ウィキの例題にある「The number 29 was included as indicated…」に当たる処理が含まれていないが、この処理はメルセンヌ数専用のもので、P − 1 法の基本形（一般バージョン）では必要ない。メルセンヌ数を対象とする場合でも、Mp の p が B1 以下なら、何もしなくても同じことになる。
因数不明の8個のメルセンヌ合成数に対して、上記の P_minus_1_test を試すと、一瞬で
が得られる。発見された約数は5兆6257億台（13桁）であり、試行除算で見つけようとすれば相当時間がかかると見込まれる。
この方法がうまくいくのは、その約数から 1 を引いた数、すなわち g − 1 が
となって、小さな素因数しか持たないから。
基本的には: もし g − 1 のどの素因数も（べき乗の形のものはべき乗込みで）B1 以下なら、このアルゴリズムによってうまく g が得られる。別の例として、コード4.1 で B1 = 5000 とした場合:
var M157 = TinyBigInt.pow( 2, 157 ).sub( 1 ); // 2^157 - 1 var g = P_minus_1_test( M157 ); Append( M157 + " is divisible by " + g );
このとき、発見された約数 g から 1 を引いた数は
となっていて、右辺の積はどの因子も 5000 以下。具体的に、B1 = 4522 では駄目だが、B1 = 4523 以上なら条件が満たされ、うまく g が得られる。
g − 1 がどんな数であれ、それが持つ素因数の大きさは有限なのだから、B1 = 500, 1000, 1500, 2000, ... などと増やしながらアルゴリズムを繰り返せば、いつかは条件が満たされ、約数が見つかるはずだ。
ただし:
B1 が大きくなるにつれ計算時間が長くなるので、現実的には限度がある。g − 1 が非常に大きな素因数を持つ場合、それが壁になる。B1 はこの問題を引き起こしやすい。第1の点は仕方ないとして、第2の点について: このアルゴリズムでは、B1 が増加するにつれて、「N（分解したい数）の約数 g のうち、g − 1 の最大素因数（べき乗込み）が比較的小さいもの」ほど早く条件を満たし、検出の網に掛かる。この順序は g そのものの大小とは必ずしも一致しない。さらに、複数の約数（素数）が同時に条件を満たして、それらの積（合成数）が g として検出されることもある。
検出された約数が合成数だとしても、素因数分解を進める上で前進であることには変わりないし、通常、その約数をあらためて素因数に分解することは容易だろう（単に B1 を下げて、もっと小さい約数を見つければいい）。しかし「P − 1 法で得られた約数は、合成数かもしれない」という点は注意を要する。
コード4.1 で得られた 5兆6257億6724万8687 も、コード4.2 で得られた 8億5213万3201 も、自動的に素因数になるわけではないが、これらは「前の章で2億まで試行除算をやって、因数が見つからなかったメルセンヌ数」の約数なのだから、2億未満の因数を含んでいない。つまり、もし約数が合成数（2個以上の因数の積）なら、最低でも2億の2乗（つまり4京＝17桁）。従って前記の2数は素数と分かる。
最初から大きな B1 を試すと、計算に時間がかかるばかりか、必要以上に g が合成数になりやすい。試行除算で簡単に割れるような小さな数についてあえて P − 1 法を試す場合、一般的に言って、B1 = 500 では既に大き過ぎる。
B1 を増やしながら何度か計算する場合、毎回最初から a_E を計算せず、新しい a_E と古い a_E の差分だけ計算すれば、時間を節約できる。別の問題点として、コード4.1のように真面目に logq を計算するのは遅いし、丸め誤差による予期せぬ挙動も心配なので、整数演算だけで済むようにした方がいいだろう。
function P_minus_1( N ) {
var B1_max = 60000;
var q_Eq_table = [];
// var q_Eq_table = new Array( 6056 + 1 ); // *0 Prime 59999 has index 6056
// var q_Eq_table = new Uint16Array( 6056 + 1 ); // *0 ES6
var a_E = new TinyBigInt( 3 );
var B1, i, q, q_Eq, old_q_Eq, g;
if( B1_max >= PrimeTable[ PrimeTable.length - 1 ] ) // *1
throw new RangeError("PrimeTable is too small for B1_max=" + B1_max);
for( B1 = 500; B1 <= B1_max; B1 += 500 )
{
for( i = 0; ( q = PrimeTable[ i ] ) <= B1; ++i )
{
q_Eq = old_q_Eq = q_Eq_table[ i ] || 1; // old value of q^Eq, e.g. 2^8=256
while( q_Eq * q <= B1 ) {
a_E = _powmod( a_E, q, N ); // general version (slow)
q_Eq *= q;
}
if( q_Eq !== old_q_Eq ) q_Eq_table[ i ] = q_Eq; // new value of q^Eq, e.g. 2^9=512
}
g = _gcd( a_E.sub( 1 ), N );
if( g.compare( 1 ) !== 0 ) break; // if g ≠ 1
}
return g;
}aE の指数部分 E は「素数の累乗 q_Eq」の積。その構成要素である各 q_Eq の現在の値を追跡し、常に q_Eq_table に記録して、次回以降は、それとの違いが生じる場合だけ、その分を q 乗している。初期状態ではテーブルの各要素は未定義値であり、これは 1 を表すものとする（q の 0 乗）。初期状態を別にすれば、配列の各要素は（q の肩の指数ではなく、q をその回数掛け合わせた結果の）q_Eq になる。
例: 仮に B1 = 20 なら q_Eq_table の先頭部は:
[16, 9, 5, 7, 11, 13, 17, 19, undefined]
これは E = 24 × 32 × 5 × 7 × 11 × 13 × 17 × 19 に当たり、24 = 16 や 32 = 9 や 51 = 5 などが、それぞれ一つの q_Eq である（素数 q の Eq 乗。例えば、素数 2 の E2 乗で、このとき E2 = 4）。
B1 = 25 になったとすると:
[16, 9, 25, 7, 11, 13, 17, 19, 23]
最初と比べると、違いは q = 5 と q = 23 が1回ずつ追加されたことだけ。大きくなった B1 に対応する新しい aE を求める場合、一から計算してもいいのだが、古い aE を「追加された q」乗すれば同じことになる。この例で言えば、a_E = _powmod( a_E, 5, N ) と a_E = _powmod( a_E, 23, N ) を実行するだけでいい。
このテーブルについては単純に [] で初期化しても構わないが、やや大きいテーブルであり、使いたい長さが決まっているので、*0 のように要素数を予約しておくといいだろう。ES6 に対応している環境なら、Uint16Array という手もある。
aE は理論上の（観念的な）値で、実際にはどこにも保存されていない。それを N で割った剰余だけが a_E として維持されている。
B1 以下の素数 q を次々と素数表 PrimeTable から読み出すため、PrimeTable はそれに対応する規模にしておく必要がある（B1 の最大値を超える素数が1個は必要）。素数表が小さ過ぎると、*1 の例外が発生する。216 まで素数表を作ってある場合、B1 = 65500 くらいまでは問題ない。
次の4個は試行除算では割れなかったが、コード4.3によって因数が見つかる。発見された因数は、順に10桁・9桁・13桁・9桁。
M103 (32 digits) is divisible (P_minus_1: 101.2 ms) M103 (32 digits) is divisible (P_minus_1_Mp: 26.6 ms) 10141204801825835211973625643007 = 2550183799 × 3976656429941438590393 M109 (33 digits) is divisible (P_minus_1: 548.1 ms) M109 (33 digits) is divisible (P_minus_1_Mp: 108.7 ms) 649037107316853453566312041152511 = 745988807 × 870035986098720987332873 M139 (42 digits) is divisible (P_minus_1: 12.5 ms) M139 (42 digits) is divisible (P_minus_1_Mp: 2.8 ms) 696898287454081973172991196020261297061887 = 5625767248687 × 123876132205208335762278423601 M157 (48 digits) is divisible (P_minus_1: 135.8 ms) M157 (48 digits) is divisible (P_minus_1_Mp: 18.6 ms) 182687704666362864775460604089535377456991567871 = 852133201 × 214388671221792782576324712513502190671
P_minus_1 は、コード4.3そのもの。P_minus_1_Mp はメルセンヌ数に最適化された改造版（次の章のコード5.6）。得られる結果は同じだが、速度に違いがある。
やれば試行除算で割れる場合でも、こちらのアルゴリズムを併用する方がいい。例えば2億まで Find_Factor2 を続けて M67 の約数を見つけると…
M67 is divisible by 193707721 (Find_Factor2: 409.9 ms)
…のようになるが、試行除算の上限を100万に制限して早めに失敗させ、失敗したら P_minus_1[_Mp] を呼び出すようにすると、トータルではずっと速くなる:
M67 is divisible by ? (Find_Factor2: 2.7 ms) M67 is divisible by 193707721 (P_minus_1: 32.0 ms) M67 is divisible by 193707721 (P_minus_1_Mp: 9.8 ms) 147573952589676412927 = 193707721 × 761838257287
一方、p < 200 のメルセンヌ合成数のうち、上記 P_minus_1 をそのまま使っても割り切れない頑強な4個の数は:
M101 = 2 5353012004 5645880299 3406410751 (31 digits) M137 = 17 4224571863 5204932932 4779900506 5324265471 (42 digits) M149 = 71362 3846352979 9405291429 8472474756 8191373311 (45 digits) M199 = 8034690221 2949513777 0981046170 5813012611 0149689139 6417650687 (60 digits)
このうち M101 は、素数表を大きくして B1 = 280000 程度まで進めれば、簡単に7兆台（13桁）の因数が見つかる。このアルゴリズムの延長であるステージ2を併用すれば軽快だ。
今回はステージ2には立ち入らないが、付録のソースコード（mer.js）には、ステージ2の複数の実装例とデモが入っている。
M199 もステージ1で割れないことはないが、ステージ2の方が効果的。1645億台の12桁の因数が見つかる。
M137 はステージ1では難しいが、ステージ2では比較的簡単に割り切れる。20桁の因数 3203京2215兆5964億9643万5569 が見つかる。
M149 だけは、ステージ2でも分解が難しい。これも20桁の因数を持つ（この数は楕円曲線法を使って分解可能）。
既に2億まで試行除算しているので、4京（17桁）以下の因数は素因数と分かる。20桁の因数の素数性については別途検討する必要がある。
結局、p < 200 のメルセンヌ数46個のうち12個が素数、34個が合成数であり、後者のうち30個（ステージ2まで含めれば33個）については、それぞれ具体的な因数を一つ知ることができた。100万まで試行除算を行うとすると、上記33個のうち10個は P − 1 法の成果。
…ということが観察される。
けれど限界もある。実際、M149 には対応できなかった。より一般的に、因数が10桁ならこれでいいとして、100桁や1000桁だとどうなるのか？
「因数が1000桁でも実用的なスピードで対応できるような、うまい分解方法があるか」という問いは、数学上の未解決問題だ。状況証拠から「うまい分解方法は原理的に存在しない」と信じる人もいるが、その状況証拠というのは「ここ50年ほど、いろいろやっても見つからないから」というだけ。逆に言えば「100年に一度の大発見」で簡単にひっくり返り、まだどちらに転ぶか分からない。
メルセンヌ数に限らず一般の合成数に対して P − 1 法を使うことができる。条件さえ満たされれば、約数が得られる。一例として、「1」を59個並べた数を分解してみる:
var R59 = new TinyBigInt( "11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 1111" ); var g = P_minus_1( R59 ); Append( R59 + " = " + g + " × " + R59.div( g ) );
P_minus_1_Mp が JavaScript で簡潔に実装され、この速度で動作するのは、新記録ではないとしても、多分珍記録だろう。要するに「ばかばかしいので、誰もそんなことはしない」というだけだろうが、その気になってやるとしても、高速化の方法は自明ではない。
_powmod は普通の「べき剰余」だが、典型的な繰り返し2乗法とは逆向きに計算している。この向き（左から右）の方が少し得になるようだ（参考リンク: [1], [2]）。
function _powmod( base, power, modulus ) {
var bin = power.toString( 2 );
var res = new TinyBigInt( 1 );
for( var i = 0; i < bin.length; ++i ) {
res = res.mul( res ).mod( modulus );
if( bin.charAt( i ) === "1" ) res = res.mul( base ).mod( modulus );
}
return res;
}[3] には、さらなる高速化のヒントも書かれている。
しかし今回のテーマはその点ではなく、もっと根本的な部分。法がメルセンヌ数専用なら、速度差10倍以上のものすごいショートカットが成り立つ。
function _powmodMp( base, power, Mp, p ) { var bin = power.toString( 2 ); var res = new TinyBigInt( 1 ); for( var i = 0; i < bin.length; ++i ) { res = _modMp( res.mul( res ), Mp, p ); if( bin.charAt(i) === "1" ) res = _modMp( res.mul( base ), Mp, p ); } return res; }
インターフェースが変わって modulus の代わりに Mp と p を受け取るようになった以外、大きな変化はないように見える。核心は、ここから呼び出されている _modMp にある（メルセンヌ数を法とする剰余）。
数値が2進表現で格納されているとする。このとき、0 以上の整数 a を 2p で割った余りは、a の下位 p ビットに等しい。割り算を行うまでもない。
「法がメルセンヌ数 2p − 1 の場合」というのは、この「簡単な場合」と法が1違うだけなので、何かうまいことができないだろうか？
このアイデアは実を結ぶ（参考リンク: [4]）。具体的に:
mod 2p) + (a rshift p) ≡ (a mod Mp)ここで rshift は右論理シフトを表す。
例えば、法 23 = 8 による簡約（余り）と法 M3 = 7 による簡約を比較すると、大筋において: a が [0, 7] の範囲ならどちらの法でも結果は等しく、a が [8, 15] の範囲なら後者の方が1大きくなる。a が [16, 23] の範囲なら後者の方が2大きくなる。以下同様なので、後者を求めるには前者に a/2p の整数部分を足せばいい。この整数部分は a rshift p に等しい。これによって「2p − 1 による面倒な割り算」が「ビットのマスク、ビットのシフト、足し算」に帰着された。足し算の結果が Mp より大きいときは、同じ手続きを繰り返せばいい。Mp 未満になればそれが答え。例外として: 結果が Mp になった場合、Mp を法とする正しい簡約は 0 となる。
実装はオブジェクト内部に踏み込む低水準なものになる。Lucas_Lehmer と P_minus_1_Mp の両方に劇的な高速化をもたらす心臓部なので、参考までにコード例を提示する。詳細は、付録のソースコード参照。
_modMp の実装例（概要）function _modMp( a, Mp, p ) {
var ret = new TinyBigInt( 0 );
var shift = p % 26;
var r_len = ( p - shift ) / 26 + ( shift ? 1 : 0 );
var mask = ( 1 << shift ) - 1;
var tmp$ = a.$;
var q$;
var r$;
while( tmp$.length > r_len || tmp$.length === r_len && tmp$[ r_len - 1 ] > mask )
{
q$ = tmp$.slice( r_len - 1 ); // This will be tmp$ / 2^p
if( shift ) q$ = _26_rshift( q$, shift );
r$ = tmp$.slice( 0, r_len ); // This will be tmp$ % 2^p
if( ! ( r$[ r_len - 1 ] &= mask ) ) __fix_array_len( r$, r_len - 1 );
tmp$ = _26_add( q$, r$ ); // a % Mp
}
if( tmp$.length < r_len || tmp$[ r_len - 1 ] < mask || __compare3( tmp$, Mp.$, r_len - 1 ) ) ret.$ = tmp$;
return ret;
}ライブラリ内部は26ビットベース。32ビットベースの方が効率的だろうが、とにかく2進数なので上記のショートカットが成り立ち、遅くなる最大原因「任意精度の割り算」が回避される。
参考コード5.3では、入力に含まれる a.$ と出力に含まれる ret.$ が同一の配列を共有する場合がある。この記事での _modMp の使用法では、この点は問題にならない。
Simple_Lucas_Lehmer については、*** の行を次のように置き換えればいい。
function Simple_Lucas_Lehmer( Mp, p ) { // Before
var s = new TinyBigInt( 4 );
for( var i = 1; i <= p - 2; ++i ) {
s = s.mul( s ).sub( 2 ).mod( Mp ); // ***
}
return s.compare(0) === 0;
}
function Lucas_Lehmer( Mp, p ) { // After
var s = new TinyBigInt( 4 );
for( var i = 1; i <= p - 2; ++i ) {
s = _modMp( s.mul( s ).sub( 2 ), Mp, p ); // ***
}
return s.compare(0) === 0;
}Simple_Lucas_Lehmer: false (6337.0 ms) Lucas_Lehmer: false (259.4 ms)
約25倍の高速化が実現された。
こうなると今度は mul（掛け算）の遅さが気になるが、それは今後の課題としよう。
一般の合成数に対応できる P_minus_1 についても、次のように改造すれば、メルセンヌ数専用の高速バージョンになる:
function P_minus_1( N ) { // Before
// ...
a_E = _powmod( a_E, q, N );
// ...
g = _gcd( a_E.sub( 1 ), N );
// ...
}
function P_minus_1_Mp( Mp, p ) { // After
// ...
a_E = _powmodMp( a_E, q, Mp, p );
// ...
g = _gcd( a_E.sub( 1 ), Mp );
// ...
}p = 1000 くらいのメルセンヌ数において一般バージョンと比較すると、こちらも約25倍速いようだ（小さいメルセンヌ数では約4～5倍速い）。
_modMp によりパワーアップしたところで、p = 2000 までを一気に判定し、15個のメルセンヌ素数を抜き出してみた。
#1 M2 = 3 (1 digit) is prime #2 M3 = 7 (1 digit) is prime (Lucas_Lehmer: 0.0 ms) #3 M5 = 31 (2 digits) is prime (Lucas_Lehmer: 0.1 ms) #4 M7 = 127 (3 digits) is prime (Lucas_Lehmer: 0.2 ms) #5 M13 = 8191 (4 digits) is prime (Lucas_Lehmer: 0.1 ms) #6 M17 = 131071 (6 digits) is prime (Lucas_Lehmer: 0.1 ms) #7 M19 = 524287 (6 digits) is prime (Lucas_Lehmer: 0.1 ms) #8 M31 = 2147483647 (10 digits) is prime (Lucas_Lehmer: 0.2 ms) #9 M61 = 230584300921369395... (19 digits) is prime (Lucas_Lehmer: 0.6 ms) #10 M89 = 618970019642690137... (27 digits) is prime (Lucas_Lehmer: 0.6 ms) #11 M107 = 162259276829213363... (33 digits) is prime (Lucas_Lehmer: 0.3 ms) #12 M127 = 170141183460469231... (39 digits) is prime (Lucas_Lehmer: 0.3 ms) #13 M521 = 686479766013060971... (157 digits) is prime (Lucas_Lehmer: 2.8 ms) #14 M607 = 531137992816767098... (183 digits) is prime (Lucas_Lehmer: 3.9 ms) #15 M1279 = 104079321946643990... (386 digits) is prime (Lucas_Lehmer: 26.0 ms)
ここまでの所要時間は6秒少々。時間の大部分は、メルセンヌ合成数を「素数でない」と判定するためのもの。探索を続ける:
#16 M2203 = 147597991521418023... (664 digits) is prime (Lucas_Lehmer: 108.8 ms) #17 M2281 = 446087557183758429... (687 digits) is prime (Lucas_Lehmer: 119.6 ms) #18 M3217 = 259117086013202627... (969 digits) is prime (Lucas_Lehmer: 298.2 ms) #19 M4253 = 190797007524439073... (1281 digits) is prime (Lucas_Lehmer: 668.5 ms) #20 M4423 = 285542542228279613... (1332 digits) is prime (Lucas_Lehmer: 696.6 ms) #21 M9689 = 478220278805461202... (2917 digits) is prime (Lucas_Lehmer: 7486.0 ms) #22 M9941 = 346088282490851215... (2993 digits) is prime (Lucas_Lehmer: 8102.4 ms) #23 M11213 = 281411201369737313... (3376 digits) is prime (Lucas_Lehmer: 11.5 sec)
Lucas_Lehmer が1回当たり3秒・4秒…とだんだん遅くなり、ついには1回10秒を超えてしまった。時間節約のため、先に因数検索の関数を呼び、明らかな合成数はテスト前に排除するようにした。#20 から #21 への旅は、間に「合成数砂漠」があり、しかも進むにつれて Lucas_Lehmer の速度が低下するので、ちょっとしんどい。#22 は意外に #21 の近くにあった。
p = 12000 までのメルセンヌ数は全て自力判定できた。これ以降についてもやればできるが、1個1～2分かかる。素数だと知られているものについて Lucas_Lehmer を行い、素数性を追試した:
[#24] M19937 = 431542479738816264... (6002 digits) is prime (Lucas_Lehmer: 68 sec) [#25] M21701 = 448679166119043334... (6533 digits) is prime (Lucas_Lehmer: 86 sec) [#26] M23209 = 402874115778988778... (6987 digits) is prime (Lucas_Lehmer: 109 sec)
テスト環境においてはこの先に壁があり、27番目のメルセンヌ素数 M44497（13395桁）については、現在の実装では判定に時間がかかり過ぎる（推定所要時間＝6時間）。メモリー的な壁、掛け算の壁、またはその両方だろう。もちろん「最先端」のメルセンヌ素数は、もっと先にある。2016年1月には49個目が見つかったという（Mersenne Prime Discovery - 2^74207281-1 is Prime!）。
2016年6月12日追記: その後、掛け算の工夫などによりパフォーマンスが改善され、M44497 については一応、約130秒で判定可能となった。
Lucas–Lehmer テストの結果の数値 s は Lucas–Lehmer 剰余（L-L residue）と呼ばれる。テストでは「それが 0 になるか」だけが問題だが、具体的な数値を見てみよう:
function Lucas_Lehmer1( Mp, p ) {
var s = new TinyBigInt( 4 );
for( var i = 1; i <= p - 2; ++i ) s = _modMp( s.mul( s ).sub( 2 ), Mp, p );
return s;
}
var p = 1277;
var Mp = TinyBigInt.pow( 2, p ).sub( 1 ); // Mp = 2^p - 1
Append( "M" + p + " = " + Mp.toString(0) );
var s = Lucas_Lehmer1( Mp, p );
Append( "L-L residue = " + s.toString(0) );
// Append( "L-L residue = " + s.toString(16) ); // *1
// Append( "L-L residue = " + s.toString(16).slice(-16) ); // *2
// Append( "L-L residue = " + _hex64( s ) ); // *3*1 のように16進数で取り出すこともでき、*2 のようにその末尾16桁（64ビット）を抜き出すと、mersenne.org や mersenne.ca のデータと直接比較できる（*3 は末尾64ビットを取り出すショートカット）。上記の例では 5613a480590e78ba だが、これは Exponent Status に出ている値と一致し、計算に間違いがないことが確認される。
同じページには興味深いデータがいろいろと載っている。B1 を2兆台まで上げて P − 1 法を試したものの、因数が見つからなかったこと。2016年には 263（約922京3372兆）まで試行除算を進めたが、割り切れなかったこと。ECM（楕円曲線法）で猛攻撃を掛けているものの、難攻不落でいまだに一つの因数も出てこないこと…。
こうまで割り切れないところを見ると、ひょっとして M1277 は素数なのだろうか？
もちろん、そうではない。Lucas_Lehmer1 の結果はゼロではないのだから。
「素数ではない」ことは分かっているのだが、「では割ってみせろ」と言われても誰にも割れない。まさに割り切れない気分。
メルセンヌ素数よりメルセンヌ合成数の方が、多くの未知を秘めている。p = 10000 程度のメルセンヌ素数の検出は、この通りブラウザ上の JavaScript でもできるが、同じ範囲のメルセンヌ合成数の素因数分解は非常に困難。
因数が知られているからといって、分解が完了するとも限らない。例えば M1249 の因数をいくつか見つけることは容易だが、2016年現在、余因数は依然として合成数だ（参考リンク: [5], [6]）。
p = 1000 を少し越えたあたりから、こういった「未解決メルセンヌ合成数」がごろごろしている。
メルセンヌ素数探しは確かに魅惑的だけれど、それは「頑張ればできる」というゲームバランスの良いところで遊んでいるようなもの。真の難題はメルセンヌ合成数の分解だ。
例えばの話、宇宙人の先生に「この400桁の数には、この200桁の因数がありますね」と言われたとき、地球人の生徒はすぐそれを検算して「確かにその通りです」と認識できる。ところが「では練習として、今度はこの別の400桁の数を自分で分解してみましょう」と言われたとき、現在の地球の技術では、1万台のコンピューターを使って丸1年ぶっ続けで計算しても「できません」となる。不名誉なことだ。「割り方（アルゴリズム）を教えてください！」と頼みたい。
「大きな数が自由に因数分解できると一部の暗号システムが崩壊するから、できない方がいい」と考える人もいるが、その発想は後ろ向き過ぎる。
[4/7] のアルゴリズムも楕円曲線法も、50年前には存在しなかった。素数性判定のアルゴリズムについては、21世紀に入ってから重大な発見があった。まだまだ良いアルゴリズムが見つかる余地はあるだろう。
2進ベースの任意精度ライブラリを（今回の件とは無関係に）作ってみたものの、最初は2進数と10進数の変換の遅さに不満を感じていた。「こんなことなら10進ベースの方が良かったな」と多少後悔していた。
でもメルセンヌ数を法とする剰余演算 _modMp の高速化は、10進数ではできなかったこと。最大公約数の計算も2進ベースで行っており、今回は2進数がプラスに働いた。こうなるとコロリと気分が変わり、「やっぱり2進数で良かった・何事もやってみるものだな」と感じるようになった。
今度は「掛け算の遅さが壁になっている」と感じられるようになった。同時に「掛け算が遅いことが問題になるなんて、何とぜいたくな悩みだろう」と逆説的な優雅さも感じた。掛け算を速くするアルゴリズムは存在するが、今までは割り算が遅過ぎてそれどころではなかった。
「遅い」という割り算だが、数万桁の数に対する試行除算は意外に速かった。P − 1 法より速い。メルセンヌ数特有の事情と、掛け算の速度の問題で、こうなるのだろう。場合によって PARI/GP より速いこともある。恐らく PARI は与えられた数を「一般の整数」として扱っているのだろう。こちらは「メルセンヌ数」という前提で近道をしている。
既存の汎用ライブラリをそのまま流用したが、メルセンヌ数専用に最初から作れば、違った結果も出せるだろう。ライブラリは比較的小さな数に最適化されている上、インターフェース用の余計なレイヤーを持っている。
探索そのものや、実装のタイムトライアル的な部分より、「いじってみることで、いろいろな事柄が感覚的に分かるようになる」ところが一番楽しい。勘や思い付きで適当にやったことでも、意外と「やってみると気付くこと」が多い。
数学の神秘と永遠は、無常な世にたたずむ私たちに、いささかの慰めと励ましを与えてくれるだろう。
コードのテスト中、M691（209桁）の因数が2個検出された。一つは「2833𥝱6377垓2442京」台の28桁の数で、106秒のステージ1と、6秒のステージ2により計算された。
https://yosei.fi/demo/TinyBigInt/
この記事の初版（2016年6月5日）では以下が使われた。
_26_.js: コア100行のプチ汎用ライブラリ、バージョン0.2（2016年6月4日版）。mer.js: 本文で紹介されたコードの他、記事中では省略された細部、デモ用コード、P − 1 法ステージ2の実装を含む（2016年6月5日版）。mersenne.html: デモ用コードを呼び出す HTMLソースの例。呼び出しの大部分は、初期状態ではコメントアウトされている。試したい部分について1回に1個だけコメントアウトを外し、実行後またコメントアウトするといいだろう。本文もソースもパブリックドメイン。
2016年6月19日追記: mer.js（2016年6月5日版）の _gcd は、2個目の引数が偶数だと正しく動作しなかった（記事中ではこの数は常に奇数なので、発覚していなかった）。このバグは2016年6月13日に修正された。
2016年8月14日追記: 以前の mer.js の P_minus_1_Mp には「B1_max が500の倍数」という制限があった。この制限は2016年8月10日版で除去された。
