カテゴリー: 多倍長数値計算入門

「多倍長数値計算入門（仮）」執筆ノート

投稿日:

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

[ 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)計算という言葉を使うケースもあったが，無限の精度をコンピュータで扱うことができない以上，過大表現っぽくて著者個人としては好きになれず，使用したことはない。もし見かけた時には本書で扱う多倍長計算と同一のものと考えてよい。

投稿日:

メルセンヌ素数で”FLOPS vs. Memory Bandwidth”を体験する

[mathjax]
　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");
投稿日:

破綻する区間演算の事例

[mathjax]
　精度保証計算（という言葉はあまり好きではないのだけれど）の事例も紹介しようかと，ARB（半径）とMPFI（区間）を使ってロジスティック写像[latex]x_{i+1} := 4x_i(1-x_i)[/latex]の計算をさせてみたところ，桁数が足りないと破綻することが分かった（今更）。つーことで，128bits計算の事例（上）と256bitsの事例（下）をば。

　有効桁数がなくなった途端に区間や半径が急激にデカくなって破綻していることが一目でわかるのはいいとして，精度はともかく計算は最後まで実行したい時には使えん。とりあえず，こういう桁落ちが酷い事例には区間演算だろうとそうでなかろうと桁数増やしてやんなさいという結論になるな。

　あと，実係数2次方程式の事例も作ってみたが，判別式で実数解と複素数解の区別をするif文をどうしようとか，実用に当たってはやはりメンドクサイことを考えないといけないのでその辺をどうかして工夫すると研究になる訳で，こっから先は精度保証の方々にお任せしようっと。まぁ「最初から複素係数でやれ」という結論になるんだろうけど，それはARBの演習問題にしておけばいいか。

　つーことで，使用したプログラム(logistic_arb.c, logisitc_mpfi.c)を最後に掲載しておく。使用時にはARB, MPFIの他，MPFR/GMPが必要。

// logistic 写像
// ARB版
#include <stdio.h>
#include "arb.h"

int main()
{
	int i;
	slong prec;
	arb_t x[102];

	// 初期化
	for(i = 0; i < 102; i++)
		arb_init(x[i]);

	printf("prec(bits) = "); scanf("%ld", &prec);

	// 初期値
	//x[0] = 0.7501;
	arb_set_str(x[0], "0.7501", prec);

	for(i = 0; i <= 90; i++)
	{
		if((i % 10) == 0)
		{
			printf("%5d, ", i);
			arb_printd(x[i], 17);
			printf("\n");
		}

		//x[i + 1] = 4 * x[i] * (1 - x[i]);
		arb_sub_ui(x[i + 1], x[i], 1UL, prec);
		arb_neg(x[i + 1], x[i + 1]);
		arb_mul(x[i + 1], x[i + 1], x[i], prec);
		arb_mul_ui(x[i + 1], x[i + 1], 4UL, prec);

	}

	// 消去
	for(i = 0; i < 102; i++)
		arb_clear(x[i]);

	return 0;
}

// logistic 写像
// MPFR & MPFI版
#include <stdio.h>
#include "mpfr.h"
#include "mpfi.h"

int main()
{
	int i;
	unsigned long prec;
	mpfi_t x[102];
	mpfr_t relerr;

	printf("prec(bits) = "); scanf("%ld", &prec);
	mpfr_set_default_prec(prec);

	// 初期化
	mpfr_init(relerr);
	for(i = 0; i < 102; i++)
		mpfi_init(x[i]);

	// 初期値
	//x[0] = 0.7501;
	mpfi_set_str(x[0], "0.7501", 10);

	for(i = 0; i <= 90; i++)
	{
		if((i % 10) == 0)
		{
			printf("%5d, ", i);
			mpfi_out_str(stdout, 10, 17, x[i]);
			mpfi_diam(relerr, x[i]);
			mpfr_printf("%10.3RNe\n", relerr);
		}

		//x[i + 1] = 4 * x[i] * (1 - x[i]);
		mpfi_ui_sub(x[i + 1], 1UL, x[i]);
		mpfi_mul(x[i + 1], x[i + 1], x[i]);
		mpfi_mul_ui(x[i + 1], x[i + 1], 4UL);

	}

	// 消去
	mpfr_clear(relerr);
	for(i = 0; i < 102; i++)
		mpfi_clear(x[i]);

	return 0;
}
投稿日:

GNU MP(GMP)マニュアル翻訳日記（最終）

　

　「GMP」だと全然関係ない分野の略語と被ることが多いようなので，「GNU MP」も併記して書くことにしようっと。そういえば「gmp_ja」で検索するとこちらの方のページが既にあったことに気づかされたのでご紹介しておく。GNU MPを全然知らん方は目通ししてから拙訳（HTML，PDF，ソース)を読むことをお勧めします。

　さて，翻訳完了から日が経つと忘れそうなことを書きつけておく。

１．GMPの開発チーム
　って表紙(PDF)に書いてあるんだもんな，誰だか知りたくなるってもんであるが，これについてはAppendix Aにまとめてある。ズラズラ時系列的に並べてあるけど，Tröbjorn以外の主要人物を抜き出すと，Paul Zimmermann(仏，MPFR開発陣の一人でもある），Kevin Ryde(豪，Perl ModuleやEmacsスクリプト著作多数)，Niels Möller（瑞），Marco Bodrato（伊，Toom-Cookアルゴリズムに詳しい)，Marc Glisse（仏，最近ではmini-gmpの提唱と実装)あたりが主要開発者ということになる。今でもGMPのMLでお見かけするメンツで，保守管理もマメにやり，中核人物の狷介さを和らげるクッションにもなっているわけで，実力もさることながら人間的も優れている（MLを見る限りだけど）貴重な人材である。

２．イマイチ理解できてない部分
　やっぱり最後の「アルゴリズム」の章が難物で，ここについては数論アルゴリズムと各CPUのアーキテクチャに通じている方（つーたら日本だとかなり稀少な部類な絶滅危惧種）にリライトか再翻訳をお願いしたい。特にアセンブラで書かなきゃいけなくなった理由である「桁上がり・桁借り上げ処理」（素直に「キャリー操作」の方が良かったかな・・・），「ループ展開」，「アセンブラコードの書き方」なんて，実際にアセンブラ職人になってシコシコとマイクロ秒単位のチューニングを体験しないことにはスッキリ訳せないよ！　GMPマニュアルの翻訳は十数年来挫折しまくってきたけど，原因はこの辺りの解読に自信が持てないからだったりする。とはいえ，もういい加減イイ歳なので恥だろうが何だろがやっておかないと後悔するとばかりにえいやっと直訳してみたのである。とゆーことでリライトをやってくれる貴重な人材求む！

３．Texinfoとmakeinfo
　TexinfoについてはTrueRoadさんのパッケージを一部利用させて頂いたが，W32Tex環境では通らない部分があって，結局，Ubuntu 16.04にtexlive2016をフルで突っ込んでxetex使ったら何とかPDFは生成できた。この辺，最近のTeX環境に疎いのでよー分からん。HTML化もUbuntuのmakeinfo –htmlで生成。こちらはスムーズにできた。今時デカい文書をTexinfo使う方が珍しいのかしらん？
　次はMPFR 4.0のTexinfoを処理する予定なんだが，まぁ当面は上記のように処理することになりそうである。TexnicianじゃないのであまりTexinfo辺で苦労したくないんで，本来ゴシックにすべき部分も明朝になってたりするけど，この辺でご勘弁を。

４．翻訳の意義
　今時は「ない」と断言してもいいんじゃないのかな。既に古いMPFRやQDの翻訳やってみたけど，ワシの場合は完全に自己満足，オナニー的なものである。Web公開しているのはGoogleという世界最高の検索ツールが使いたいからであって，備忘録以上の意味はない。ワシの業績リストにも一切入れてないし（何人かの方からお褒めの言葉は頂いたけど）。ただ，「一通り読んでみましたよ」ぐらいの意思表示はできるかなと。多倍長数値計算の基本ライブラリであるMPFR,QD,GMPぐらいはマニュアルぐらい一通り読んでおかないとこの辺の専門家でございなんて言えないし，この辺の後追いをやるにしても，先人の業績を参考文献に使うぐらいの節度は欲しいし，参考文献に入れるなら「全部読んだ」ぐらいのことは言いたいしねぇ。
　GMPについては，もう当分これを超える多倍長自然数演算カーネルは出てこないだろうなと思うし，GPUではカーネルでmalloc/freeを繰り返すと覿面に遅くなるんで，MPNカーネルそのままの移植は向いてないだろうと感じる。CUMPは今んとこbasecase演算にとどまっているし，テンポラリメモリは固定で取ってるので，これそのままだとKaratsubaやToom-Cookの移植は大変そうである（割算だけ移植して挫折した奴がここにいる）。何にせよ，多倍長精度を素直にサポートするライブラリとしてはエベレスト的存在がGMPなので，当面この方向での進化もさしてなさそうに思える。
　GMP自体，Version 4以降はアーキテクチャの進化に追随しているだけで（それだけでも大したもんなんだけど），アルゴリズム的な大幅改善はない。Version 5はあっという間にVersion 6になっちゃったけど，メジャーバージョンアップの割には新規性はかなり薄い（だもんで不満も多少でてた）。どん詰まりっぽい状況だから，マニュアル翻訳も今やっておけばさほど古びずに済みそうという予想があってやっている訳である。オナニーなら性欲のあるうちじゃないと詰まらんしねぇ。

　とりあえず，こんなもんかしら？　後はどなたかに（翻訳する人がいれば）お任せする。