gmpの実装とC++からの利用法，性能について(1) [gmp multi-precision]

Posted on 2 mins read

はじめに

そろそろTwitterのフォロワーが2500人を迎えるそうです． やはりマイナーな世界だけを深堀りするのでなく， もっと大衆受けする活動をしていくべきだと思いました．

そう，つまりマスに訴えかけることでHPC界のヒカキンのような存在になるのです．

では皆が見てくれるような一般的な話題ってなんなのか？

私は考えました．
google analyticsの結果を見るとこのページはPEZYと4倍精度という検索ワードでくる人がほとんどのようですが， そんな一部の研究者しか使わないものでなく，もっとマスに訴えかけるような，そんなテーマを考えました．

選ばれたのはGNUでした

そう．GNUに媚びて生きていこう．
つまり今回のテーマはGMPです．

なんと四則演算やFFTについて触れます．全国民がオオウケ間違いなし!
明日から倍精度を使うのをやめ，楽しい多倍長ライフが始まるはずです！

まぁ真面目に導入すると，2月頃にgccのfloat128の性能について書いたのがそこそこ閲覧数が多かったので， 他の高精度演算ライブラリについてもちょっと書くかー，と思って第2弾としてGMPを取りあげてみるということです．．
実際はもっと前から書いていたのですが，細部が書けずにいて，2500人だしここで公開しようと思い立って今に至ります．

様々な言語で多倍長整数は標準サポートされることも多くなってきましたが， 浮動小数点に関してはほとんど見かけませんし，実装自体も少ないです．

やはり浮動小数点演算は癖が強い(指数部が伸びるのか，仮数部が伸びるのか，どのくらい時間がかかるのか)や， ツールそのものがどういった機能を持っているかが分かりにくいです．

今回はGMPに着目して，機能面や実装面について紹介していきたいと思います．

ただ，GMPにはいろいろな難しい特徴があり，正直私もあまり詳しい事は言えないので，フォロワー2500人記念とかぶち上げておいてあれですが， 今回の記事では私が現状で知っているGMPの実装についてちゃんとまとめ，使い方と性能を簡単に見ることにします． (GMPで論文何件か書いたとは思えないほど内部の実装に自信がないおじさん)

GMPを用いたプログラム

GMPはGNU Multi Precision Arithmetic Libraryのことで，float128などが4倍精度固定であるのに対し， GMPは変数の宣言時，またはどこかでデフォルトの仮数部の精度を指定することで任意の精度型を作ることができます．
ただし，仮数部を自動的に伸ばしたり，自動的に選択するような機能はついていません(そんなことをしたらメモリが大変なことに)．

指数部については固定で，仮数部を指定していくことになります． そのため例題などではよく円周率などが扱われる印象があります．

ここではC++版のgmpxx.hを使ってGMPの多倍長浮動小数点型であるmpf_tのC++ Wrapper, mpf_classを使ってみます． C++版を使えば演算子オーバーロードによって比較的簡単に実装することができます．

すべての変数が同じ精度で良い場合は，mpf_set_default_prec()という関数を使って精度を指定すればよいです． 簡単なプログラムは以下のようになります．

#include<gmpxx.h>
#include<iostream>
int main(){
	mpf_set_default_prec(1024);

	mpf_class a(1.5); // a = 1.5
	mpf_class b = 1.0;
	mpf_class c = 0.0;

	c = a + b;

	std::cout << c << std::endl;
	
	return 0;
}

ね？簡単でしょ？

こんな感じでC++のクラスを使って実装できます．最初に1024 bitを指定しているため，a,b,cは仮数部1024 bitの浮動小数点数になります．

なお，コンストラクタで精度を指定することもでき， mpf_class a(1.5,512)などとするとaだけ512 bitで宣言できます．

数学関数とかもだいたいあります．詳しくは本家を見て下さい．

GMPの実装

さて，ここではGMPの内部実装について触れます．

なぜGMPの内部実装について触れないといけないのか？というのは非常に面倒な問題があって， ここがGMPを使ったプログラムを性能評価する上で一番ややこしいところです． (性能評価しにくい＝計算時間を見積もりにくい)

何が面倒かというと，変数の状態によって内部の実装が大きく切り替わるようになっている点です． そのため変数の状態が計算時間にダイレクトに効いてくるため，小さい系の実行時間を参考にしても大きい系の実行時間を見積もれません．

githubのGMPのミラーリポジトリを例に説明します(cc09e59afe9a23de9e7d0ef10598a395d7e8c850)．

まず，congigure時にgmp.hやgmp-mparam.hが生成される仕組みになっています．
ここで様々なアーキテクチャごとの設定が反映されたヘッダファイルが生成されます．

ヘッダ生成時に行われる大きな切り替えとして，精度によって演算に用いられるアルゴリズムが異なるということがあります．
これはXX_THRESHOLDという変数で決まるので，grepで探してみればわかります(XXはいろんなアルゴリズムの名前が入る)． 詳しい切り替えなどについては次回の記事で触れたいと思います．

何もアーキテクチャを指定しなければgmp_impl.hというヘッダファイルによって決まりますが， tune/ というところでもアーキテクチャごとに色々書いてあるようなので， 実際にどれが呼ばれるのかはconfigure時に判定される．．．のだと思います (要調査)．

ただ，ここでも更に問題があり，「精度によって切り替わる」というのが変数の初期化時に指定した精度ではなく， 実際に値が入っている桁数によって動的に切り替わるという点です． これはGMPの浮動小数点型であるmpf_t型の宣言を見るのがわかりやすいです．

typedef struct
{
	int _mp_prec;         /* Max precision, in number of `mp_limb_t's. Set by mpf_init and modified by mpf_set_prec.  The area pointed to by the _mp_d field contains `prec' + 1 limbs.  */
	int _mp_size;         /* abs(_mp_size) is the number of limbs the last field points to.  If _mp_size is negative this is a negative number.  */
	mp_exp_t _mp_exp;     /* Exponent, in the base of `mp_limb_t'.  */
	mp_limb_t *_mp_d;     /* Pointer to the limbs.  */
} __mpf_struct;
typedef __mpf_struct mpf_t[1];

これがmpf_t型の実装です．英語でも説明がありますが簡単に説明すると， mp_expは整数型(デフォルトだとint型)で，指数部です． *_mp_dは基本的にdouble型で，仮数部です．

つまりGMPはint型の指数部，およびdouble型の配列から構成される仮数部によって多倍長浮動小数点数を実装しています．

_mp_precはユーザが指定した精度です．ただしmpf_set_default_prec()には仮数部のbit数を入力するのに対し， _mp_precには_mp_dの配列の長さ+1が入ります．
つまり指数部64bitなら1+1=2, 128bitなら2+1=3, 512bitなら8+1=9です．

_mp_sizeは_mp_dの配列の何番目まで値が入っているかです．

簡単な構造体なので，データサイズなどを計算するのは簡単かと思います．

問題になるのは，_mp_sizeです． アルゴリズムを切り替えるTHRESHOLDは_mp_sizeを見るものがほとんどです． つまりどれだけ大きい数で初期化したとしても値が1であれば最も軽いアルゴリズムで，_mp_dの先頭要素だけが使われるので 演算としては一瞬で終わります．

つまり_mp_prec (=指定した精度)は，「これ以上はallocateされませんよ」という意味であって， 実際には_mp_size (=入っている値)によって挙動が変化するということです．

ここがGMPの性能が評価しにくい点です． 内部の値が仮数部の配列をどれだけ埋めているかによって演算の種類や参照する配列の範囲が動的に異なることにより， 単純なforループでも回しているうちに計算がどんどん重くなるような挙動が起こることがよくあります．

また，演算する際の精度が変数によって異なる場合にどうなるのかなど，分岐がややこしく，私もどういう条件になっているのか完全にはわかっていません．

ただし，このように値によって挙動が変化するということを覚えておかないと，演算時間がスケールしなかったり，入力によって性能が安定しないかのような印象を受けるでしょう．

(もしかするとプログラムによっては演算順序を入れ替えることによって値が変わるため性能が向上する例もあるかもしれません．特におもしろい例は思いついていませんが．)

性能評価

試しに内積を実装して時間を測ってみた．
コードはここにおきました．

gccの4倍精度と違い，ompのreductionが共通化できなかったため， templateを用いて共通化することはできませんでした．

gcpで16コアのHaswellマシンを借りて実行してみました．
gccの4倍精度のときと同じマシンです．あんまり参考にはならないが Intel(R) Xeon(R) CPU @ 2.30GHzとのこと．

gccはgcc version 8.2.1 20180905 (Red Hat 8.2.1-3)
OSはCentOS Linux release 8.0.1905 (Core)

-O3 -fopenmp をつけ，gmpはyumから落としてきた． 参考までに自分でオプションを色々変えてgmpを色々ビルドしてみたが，ほとんど性能は変わらなかった．

値については初期化後に1/3を入れることで必ず_mp_dに値がすべて詰まる状態にした． サイズはgccの4倍精度では10^9までやったが，あまりにも遅くて待てなかったので10^3と10^6だけにした．

表 各ベクトルサイズにおける倍精度と4倍精度の実行時間 [ms] (1 thread)

double gmp(1024) gmp(2048) gmp (4096)
10^3 0.013 0.33 0.90 2.63
10^6 2.14 340.3 914.6 2648.2

表 各ベクトルサイズにおける倍精度と4倍精度の実行時間 [ms] (16 threads)

double gmp(1024) gmp(2048) gmp (4096)
10^3 0.003 0.07 0.13 0.37
10^6 1.28 39.3 107.5 305.4

マルチスレッド化の効果はGMPの精度に関わらずどれもでています． 一方で精度ごとの関係は微妙で，16 threadsのときで見比べると，精度が1024->2048->4096と2倍になっても， 計算時間は2.7->2.85倍と変化しており，精度の増加とは異なる傾向があるように見えました．

今後はもう少し実験データを増やすことと，それぞれの精度で何のアルゴリズムが使われているかを調査していきたいと思います． それについては次回の課題としたいと思います．

今後の課題

  • GMPのTHRESHOLDを変更し，特定のアルゴリズムに固定した上でそれぞれの性能を評価したい
  • GMPのFFTは何を使っているんだろう．．FFTWやIntelのライブラリに置き換えたら速いだろうか？
  • QDライブラリについてまとめたい