メルセンヌ素数と仮想通貨のパラダイム

Visual StudioにはNuGetというパッケージマネージャがあって、様々なパッケージをうまくいけば一瞬で導入できます。うまくいけばインストールバトルが発生しないので便利です。（うまくいかなかった場合は、Microsoftのやる気なさげな企画が発する独特の雰囲気を体験できます）例えばC++ BoostとかOpenCVとかはそれでいけます。Boostには任意精度演算機能があるので、ちょっとC++で多倍長整数を扱いたいときとかに

#include <boost/multiprecision/cpp_int.hpp> 
typedef boost::multiprecision::cpp_int Bigint;

とかでできます。任意精度演算ライブラリたちの中では比較的遅いらしいのですが手軽です。

有名な任意精度演算ライブラリのGMPというのがあって、それをWindows向けにしたmpirというライブラリをNuGetで入れられると最近知ったので試してみました。今回私が使ったのはmpir-vc140-x64という名前のパッケージです。（vc140と書いてある通り、Visual Studio 2015のコンパイラがインストールされてないと動きません。そういう雰囲気です）

とりあえずメルセンヌ数が素数かどうか判定するやつ(Lucas-Lehmer test)を書いてみました。

#include <iostream>
#include <chrono>

#include <gmp.h>

void LLtest(const uint32_t p) {

	mpz_t integer_mp;
	mpz_init(integer_mp);
	mpz_set_ui(integer_mp, 1UL);
	mpz_mul_2exp(integer_mp, integer_mp, p);
	mpz_sub_ui(integer_mp, integer_mp, 1UL);

	mpz_t integer_a1;
	mpz_init(integer_a1);
	mpz_set_ui(integer_a1, 4UL);

	mpz_t integer_lo;
	mpz_t integer_hi;
	mpz_init(integer_lo);
	mpz_init(integer_hi);

	for (uint32_t i = 1; i <= p - 2; ++i) {

		const auto start = std::chrono::system_clock::now();

		//↓プロファイリング結果これが97%以上を占めている。
		mpz_mul(integer_a1, integer_a1, integer_a1);

		mpz_add(integer_a1, integer_a1, integer_mp);
		mpz_sub_ui(integer_a1, integer_a1, 2);

		while (true) {
			const int c = mpz_cmp(integer_mp, integer_a1);
			if (c > 0)break;
			if (c == 0) {
				mpz_set_ui(integer_a1, 0UL);
				break;
			}

			mpz_tdiv_r_2exp(integer_lo, integer_a1, p);
			mpz_tdiv_q_2exp(integer_hi, integer_a1, p);

			mpz_add(integer_a1, integer_lo, integer_hi);
		}

		const auto end = std::chrono::system_clock::now();
		const double elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();

		std::cout << "iteration " << i << " / " << p - 2 << ", time = " << elapsed << " ms" << std::endl;
	}

	if (mpz_sgn(integer_a1) == 0 || mpz_cmp(integer_a1, integer_mp) == 0) {
		std::cout << "2^" << p << "-1 is a prime!" << std::endl;
	}
	else {
		std::cout << "2^" << p << "-1 is not a prime." << std::endl;
	}

	mpz_clear(integer_mp);
	mpz_clear(integer_a1);
	mpz_clear(integer_lo);
	mpz_clear(integer_hi);
}

int main() {

	const uint32_t p = 412876283;

	std::cout << "start: Lucas-Lehmer test, p = " << p << std::endl;

	LLtest(p);

	return 0;
}

手元のCore i7 4700 (Haswell, 3.4GHz)で走らせると、p = 412876283のとき1イテレーションあたり6.2秒くらいかかりました。（ちなみにLL testは入力pが素数であることを前提としていて、pが合成数のとき出力はナンセンスです。412876283が素数であることはこのコード内では確認していませんが、予め確認済みです。）

有名な話ですが、GIMPS（ https://www.mersenne.org/ ）は巨大なメルセンヌ素数に懸賞金をかけています。というか、正確に言うと、電子フロンティア財団という別組織が巨大な素数（メルセンヌ素数でなくてもよい）の発見に懸賞金をかけていて、GIMPSの参加者がGIMPSのソフト(prime95)を使ってメルセンヌ素数を発見すると、電子フロンティア財団からGIMPS運営に与えられている懸賞金の一部を運営から発見者に渡すという感じです。

世界中の参加者が計算した結果がGIMPSのサイトで公開されています。（トップページからCurrent Progress → Recent Results とか）計算時間も載っていたので、最近のやつを適当に集めてグラフにしてみました。

f:id:eukaryo:20200711094045p:plain

グラフの横軸は試行されたpで、縦軸はCPUコア数*クロック周波数*時間です。Lucas-Lehmer testではp桁の2乗をp-2回やるので、乗算にSchönhage–Strassen algorithmを使うとすると、全体の時間計算量は $O(n^{2}log(n)log(log(n)))$ です。とりあえずエクセルの二次関数のやつで適当にフィッティングしました。これで雑に外挿すると、pが4億くらいだと約 6000 (GHz*day)くらいかかると予想されます。

すごく気になるのは、prime95は前述のコードよりも異常に速いということです。前述のコードは乗算1回あたり約 3.4*6.2 (GHz*sec) だったので、単純に412876283倍してテスト全体の所要時間を求めると約100000 (GHz*day) になります。15倍以上の差があります。HaswellはAVX2演算器が128bit幅だから～とかそういうレベルではなくないですか？　GIMPSのサイトからprime95のソースコードをダウンロードできますが、ダウンロードページの下のほうには「めっちゃ読みにくいぞ」みたいなことが長々と書いてありました。

ちなみに超雑な計算ですが、 $2^{32}$ 進N桁同士の乗算に $N log_{2}(N) log_{2}( log_{2}( N ))$ 回の演算が必要だとして、1クロックサイクルあたり1回の演算ができるとして、その乗算を412876283回やるとすると、N=412876283/32のとき、約6600 (GHz*day)になります。そう考えるとprime95の速度はギリギリ実現不可能ではない気がしませんか？　（実現してるわけですが）　もちろん他にも諸説あって、gmpとmpirに何らかの性能差がある説、NuGetのmpirが本来の性能を出せてない説、上記の私のコードがタコい説とかあるでしょう。

でもよく考えると、prime95が他と比べて異様に速いという事態は、多くの人々がprime95に参加している状況そのものから予想できたことです。もしgmpとかで短く書いたコードが同等に速いなら、懸賞金狙いの人は誰も初手ではprime95を使わずに手元で計算するはずです。なぜなら予備的な実験結果やネガティブリザルトをライバルに知らせるのは不利にしかならないからです。もっと言うと、GIMPSがそういう状況を防ぎたいならば、prime95が素数を発見した場合に懸賞金を出すだけでなく、同じprime95で合成数を発見した場合にもその労力に応じて支払うべきです。これは突き詰めると仮想通貨のプールマイニングと同じような図式に帰着するでしょう。