50番目のメルセンヌ素数が見つかった，つーか，素数であることの検証が済んだ，という記事に呼応して下記のようなTweetをした。
これって，”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");
