投稿者: T.Kouya

[ 2019-10-25追記 ] いよいよ「多倍長精度数値計算」（森北出版）が2019年11月22日に発売になります。この部分も手直ししたものが収録されますんで・・・(売れそうもないので)お高いですが・・・よろしくです。

第1章の最初の記述がまとまらんので，ざっと書き出してみる。こんな感じでどうっすかね？＞Fさん

—-
1.1 そもそも「多倍長」ってどういう意味？

コンピュータ内部におけるデータは全て2進数(binary number)，即ち，0もしくは1の並びで表現される「数値」(numerical value)である。データを格納するメモリ(memory)の大きさは制限されているので，全ての数値の長さは有限でなければならない。標準的に使用される数値の形式（format），即ちデータ型(data type)は規定されており，大別して整数(integer)型と実数型(real number)型に分類される。C/C++の標準的な整数型と実数型の変数，即ち，数値を格納するために確保されるメモリ領域を使用するために宣言するデータ型は

• 整数型
  • 符号付き整数型(signed integer)
    • • int型
      • long int型
  • 符号なし整数型(unsigned integer)
    • • unsigned int型
      • unsigned long int型
• 実数型
  • • float型
    • double型

となる。現在主流の64bit CPUではint型，unsigned int型，float型が32bit長，long int型，unsigned long int型, double型が64bit長である。これらの標準的なデータ型はCPU内部のハードウェア演算回路を使用し，1命令(instruction)での高速な演算処理(arithmetic)が可能である。
本書では，これらの標準的なデータ型を複数組み合わせ，64bitより長い数値を扱う計算全般を広義の多倍長(multiple precision)計算と呼ぶことにする。もっとも主題は多倍長の実数型を使用する数値計算(numerical computation)なので，特に断らない限り，double型より長い実数型を扱う数値計算を多倍長計算と呼ぶ。これらの多倍長計算は標準的なfloat型，double型を使用した数値計算に比べて多大な計算時間を要するのが普通である。

本書で扱うのは下記の多倍長整数型と多倍長浮動小数点型である。

• 多倍長整数型(mpz_t, mpz_class)・・・整数型を可変長の配列(固定長にする場合もある)にして長い桁をサポートする
• 多倍長浮動小数点型
  • 多数桁方式(mpf_t, mpf_class, mpfr_t, mpfr::mpreal)・・・固定長の整数型の配列を仮数部(小数部)とする
  • マルチコンポーネント方式(dd_real, qd_real)・・・固定長の実数型(浮動小数点型)の配列を丸ごと利用する

多倍長有理数型(mpq_t, mpq_class)と多倍長複素数型(mpc_t, complex)については，前者は分子と分母にそれぞれ多倍長整数型を，後者は実数部と虚数部にそれぞれ多倍長浮動小数点型を使用して実装される。

ところで，本書では「多倍長」をかように定義するが，複数のデータ型を混合して計算する場合もmultiple(複数） precision（精度）と呼称しているケースもある。後述するように，現在の計算環境は標準データ型演算速度の高速化が限界に達しており，計算効率を上げるために，例えば，高速なfloat型と比較的低速なdouble型を混合させる方法も盛んである。本書で取り上げる長い実数型でも，精度は様々なものを混在させて利用することが前提となっており，その意味では「より長い精度の実数型を取り混ぜで使用する数値計算」という意味合いも「多倍長計算」に込められていると解釈してもらっても構わない。

しかし，これでは「より長い精度桁数」を利用しているかどうかは判然としないこともあるので，そこを強調したい向きには，可変精度(variable precision)計算，任意精度(arbitrary precision)計算という言葉を使用すると良い。無限精度(infinite precision)計算という言葉を使うケースもあったが，無限の精度をコンピュータで扱うことができない以上，過大表現っぽくて著者個人としては好きになれず，使用したことはない。もし見かけた時には本書で扱う多倍長計算と同一のものと考えてよい。

　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");