日: 2018年1月5日

[mathjax]
　50番目のメルセンヌ素数が見つかった，つーか，素数であることの検証が済んだ，という記事に呼応して下記のようなTweetをした。

見間違えたので再TW。古いXeonマシンでもGNU MP(GMP)のmpz_ui_pow_ui & mpz_sub_ui関数で計算して30秒足らず。ファイルDLするより計算した方が速い。素数検証の方がずっと大変。＞「史上最大の素数」約2年ぶりに更新 最長だったものより約100万桁大 #ldnews https://t.co/lO9BwReyng pic.twitter.com/jEzQbfQhZu

— 幸谷智紀＠駿府(+髭) (@tkouya) 2018年1月5日

　これって，”FLOPS vs. Memory Bandwidth”に通じる話でもあるなと思ったので，ここでもう少し詳細にメモっておく。上記のTweetより少し手直ししたプログラムは下記の通り。

// Mersenne Prime Number Calculator by GNU MP
// Copyright (c) 2018 Tomonori Kouya
//
// Compile: cc mpz_mersennne.c -lgmp
//
#include <stdio.h>
#include <stdlib.h>

#include "gmp.h"

/* Mersenne Prime Number */
/* 2^p - 1

/* cf.) http://primes.utm.edu/mersenne/ */

int max_index = 50; // 2018-01-05
unsigned long known_mersenne_prime_p[] = {
	       2,        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,    44497,    86243,   110503,   132049, \
	  216091,   756839,   859433,  1257787,  1398269,  2976221,  3021377,  6972593, 13466917, 20996011, \
	24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917 \
};

int main(void)
{
	mpz_t a;
	int index;

	mpz_init(a);

	// input
	fprintf(stderr, "2^p[index] - 1 = ? ");
	fprintf(stderr, "index(1 - %d) = ", max_index);
	scanf("%d", &index);
	if((index > max_index) || (index <= 0))
	{
		fprintf(stderr, "ERROR: index = %d is illegal!\n", index);
		return EXIT_FAILURE;
	}

	index--;
	printf("p[%d] = %ld\n", index + 1, known_mersenne_prime_p[index]);

	// a := 2^p - 1
	mpz_ui_pow_ui(a, 2UL, known_mersenne_prime_p[index]);
	//gmp_printf("2^%ld     = %Zd\n", known_mersenne_prime_p[index], a);
	mpz_sub_ui(a, a, 1UL);
	gmp_printf("2^%ld - 1 = %Zd\n", known_mersenne_prime_p[index], a);

	mpz_clear(a);

	return 0;
}

　Tweetに挙げたのは\(2^p\)(mpz_ui_pow_ui関数で計算)と\(2^{p – 1}\)(mpz_sub_ui関数で計算)の両方をファイルに書き出す奴だったので，後者のみ書き出すようにするとほぼ半分の時間で済む。

　10進数テキストで書き出すと大体45MBぐらいのファイルになる訳だが，ネット越しだとダウンロードする方が時間がかかる。100BASE-TXの環境でもイーブン。大体，計算そのものにかかる時間よりファイルに書き出す時間の方が問題で，更に言うと，メインメモリ内部のデータ移動の法に時間が費やされている。128MBぐらいのL3キャッシュがあれば，FLOPS(この場合は整数演算だからMIPSだけど)上げるよりもずっと高速に計算終了している筈で，さらに言えば，ファイルに書き出しせずにオンメモリで済むならもっと高速。つまりCPU内部の演算効率よりメモリ帯域(memory bandwidth)の方がよっぽど問題，というお話になるわけで。

　という自分用メモ。

[2018-01-26(Fri) 追記] 「多倍長数値計算入門（仮）」執筆開始につき，上記のCプログラムをJuliaスクリプトで書き直してみた。なーるほど，これは使いやすいや。

#  Mersenne Prime Number Calculator by GNU MP
#  Copyright (c) 2018 Tomonori Kouya
# 
#  Run: julia mersennne.jl 
# 

# Mersenne Prime Number
# 2^p - 1

# cf.) http:# primes.utm.edu/mersenne/

max_index = 50; #  2018-01-05
known_mersenne_prime_p  = [
	       2,        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,    44497,    86243,   110503,   132049,
	  216091,   756839,   859433,  1257787,  1398269,  2976221,  3021377,  6972593, 13466917, 20996011,
	24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917
];

#  input
print("2^p[index] - 1 = ? ");
print("index(1 - $max_index) = ");
index = parse(Int, chomp(readline(STDIN)));
if (index > max_index) || (index <= 0)
	print("ERROR: index = $index is illegal!\n");
	return EXIT_FAILURE;
end

println("p[$index] = ", known_mersenne_prime_p[index]);

#  a := 2^p - 1
a = BigInt(2)^known_mersenne_prime_p[index] - 1;
println("2^", known_mersenne_prime_p[index], " - 1 = $a");