GMPの仕組み

前回の記事でGMPについて紹介してみた。
そこでGMPが内部的にどのように任意精度整数を保持しているのか気になったので、調べてみた。

ものぐさな人へ

最後まで目を通すのが面倒な人のために、結論だけ書いておきます。
GMPは内部的には(2^64)進数で値を保存しているようです。

まずはサンプル

前回では触れませんでしたが、C言語でGMPを使うときの基本的な形は次のような感じになります。

C言語でGMPを使うときは、C++用拡張モジュールgmpxxを使用していないので、-lgmpフラグのみでコンパイルできます。
コメントで説明をしておいたので、流れはつかめると思います。

デバッグしてみる

コンパイル後、gdbで中身を覗いてみます。

$ gcc -o gmp gmp.c -lgmp -g
$ gdb gmp
GNU gdb 6.3.50-20050815 
#〜〜〜〜〜（中略）〜〜〜〜〜
done

(gdb) b main
Breakpoint 1 at 0x100000e81: file gmp.c, line 7.
(gdb) r

#〜〜〜〜〜（中略）〜〜〜〜〜

Breakpoint 1, main () at gmp.c:7
7		mpz_init(a);
(gdb) n
Current language:  auto; currently minimal
8		mpz_set_str(a, "48410242354393281104234213124421033", BASE);
(gdb) p a
$1 = {{
    _mp_alloc = 1,
    _mp_size = 0,
    _mp_d = 0x100200000
  }}
(gdb)

まずは、mpz_tという構造体には3つのプロパティがあることが見てとれます。
名前から察するに、_mp_allocがアロケートサイズ、_mp_sizeがこの整数の大きさ、_mp_dがポインタといったところでしょうか。
実体が_mp_dなので、これを覗いてみます。

(gdb) x 0x100200000
0x100200000:	0x00000000
(gdb) x/10b 0x100200000
0x100200000:	0x00	0x00	0x00	0x00	0x00	0x00	0x00	0x00
0x100200008:	0x00	0x00
(gdb) x/10w 0x100200000
0x100200000:	0x00000000	0x00000000	0x00000000	0x00000000
0x100200010:	0x00000000	0x00000000	0x00000000	0x00000000
0x100200020:	0x00000000	0x00000000
(gdb) x/10d 0x100200000
0x100200000:	0	0	0	0
0x100200010:	0	0	0	0
0x100200020:	0	0
(gdb) x/10g 0x100200000
0x100200000:	0	0
0x100200010:	0	0
0x100200020:	0	0
0x100200030:	0	0
0x100200040:	0	0
(gdb)

空ですね。
ちなみに、gdbはxコマンドでメモリの中身を覗くことができます。
x/〜〜のスラッシュの後にフォーマットを指定することができます。
8行目のmpz_set_strを実行して、どうなるか見てみます。

(gdb) n
9		mpz_out_str(stdout, BASE, a);
(gdb) p a
$1 = {{
    _mp_alloc = 3,
    _mp_size = 2,
    _mp_d = 0x100103ad0
  }}
(gdb) x 0x100103ad0
0x100103ad0:	0xe11d35a9
(gdb) x/10b 0x100103ad0
0x100103ad0:	0xa9	0x35	0x1d	0xe1	0xa2	0x55	0x4d	0xd2
0x100103ad8:	0x2c	0x8c
(gdb) x/10w 0x100103ad0
0x100103ad0:	0xe11d35a9	0xd24d55a2	0x29648c2c	0x000952cf
0x100103ae0:	0x00000000	0x00000000	0x00000000	0x00000000
0x100103af0:	0x00000000	0x00000000
(gdb) x/10d 0x100103ad0
0x100103ad0:	-518179415	-766683742	694455340	611023
0x100103ae0:	0	0	0	0
0x100103af0:	0	0
(gdb) x/10g 0x100103ad0
0x100103ad0:	-3292881594488113751	2624324496559148
0x100103ae0:	0	0
0x100103af0:	0	0
0x100103b00:	0	0
0x100103b10:	0	0
(gdb) x/10gu 0x100103ad0
0x100103ad0:	15153862479221437865	2624324496559148
0x100103ae0:	0	0
0x100103af0:	0	0
0x100103b00:	0	0
0x100103b10:	0	0
(gdb)

_mp_dの参照先アドレスがちょっと変わっています。内部的にリアロケートが起こったといったところでしょうか。
確かに何かデータが格納されたようですが、これだけでは何がなんだかわかりません。

ソースコードを読んでみる

仕方ないので、一旦gdbでのデバッグは諦めて、ソースコードを読んでみます。
公式サイトの方からソースコードは落とせますので、読んでみてください。

The GNU MP Bignum Library

解凍したら、自分の環境でconfigureをかけてください。

まずは、configureで生成されたgmp.hを覗いて、mpz_initに辿り着いてみました。

/gmp/mpz/init.c

mpz_ptrというのは、mpz_tのポインタです。また、

ALLOC (x) => x->_mp_alloc
PTR (x) =>x->_mp_d
SIZ (x) => x->_mp_size

とマクロ定義されています。
__gmp_allocate_funcはmallocのラッパでした。
つまり、mpz_initは本当に単なるmallocのラッパと構造体の初期化処理を行っているということになります。

まずはmpz_tという型がなんなのか見てみます。

mpz_tはこの__mpz_structをtypedefしています。
だいたいgdbで見た通りの中身です。
ここで、mp_limb_tとはなんなのか気になります。実際に探してみると、単にunsigned longのtypedefでした。
というわけで、GMPはunsigned long単位で保持していることになります。
が、「/* Pointer to the limbs. */」という少し気になるコメントが。
公式ドキュメントをひっくり返して有益な情報がないかと調べたところ…

Nomenclature and Types - GNU MP 5.1.3

A limb means the part of a multi-precision number that fits in a single machine word. (We chose this word because a limb of the human body is analogous to a digit, only larger, and containing several digits.) Normally a limb is 32 or 64 bits. The C data type for a limb is mp_limb_t.

（limbはひとつの機械語と対応する、任意精度整数の一部分だよ！）

わけわかんねえ…

しかたないのでソースコードをおとなしく読みます。

では続いて、mpz_set_strを探してみます。

/gmp/mpz/set_str.c

この関数、いろいろ長々と処理を行っているようで、実は大したことをしてないです。
だいたいは空白スキップとか、マイナスを処理したりとか、基数を判定したりとか。
で、最終的に97行目でmpn_set_strという関数に処理を投げています。
引数を確認しておくと、

PTR(x)：x->_mp_d
begs：渡された任意精度整数の文字列strのアドレス
str_size：strlen(str)
base：基数

です。では、mpn_set_strを見ます。

この関数は基数に応じて処理を分岐させています。
詳しく読めばわかるんですが、今は簡単に言うと基数が2のベキ乗のときはmpn_set_str_bits、それ以外のときはmpn_set_str_otherを呼びます。
とりあえず、今は10進数を考えるとして、後者の方を参照します。

これが文字列を実際に変換している部分になります。
この関数のポイントは3つ。

11行目のfirst
13〜16行目の処理
28〜31行目の処理

最初に全体的な流れを説明しておきたいと思います。
この関数でやっていることは、任意精度整数を上の桁から数桁単位で切り取って、変換しています。

f:id:levelfour:20140122215651p:plain

切り取り方ですが、基本的には19桁単位で切り取ります。
なぜ19桁かというと、

$10^{19} < 2^{64} < 10^{20}$

すなわち、unsigned longに入る最大の10のベキが10^19だからです。
そして、下の位から19桁ずつ区切り、残った端数桁の部分がfirstになります。
このfirstの桁数を11行目で計算しています。

次に、13〜16行目です。これは、firstをまず変換しています。
spには各桁の数字が入っていますので、結局何をしているのかというと

$\cdots ((a_n\times 10+a_{n-1})\times 10 + a_{n-2})\times 10 +\cdots \\= a_n\times 10^{n-1} + a_{n-1}\times 10^{n-2}+\cdots+a_1\cdots 10^0 \\= \sum_{k=1}^{n}a_k\times 10^{k-1}$

（ただしa_n=firstの下から数えてn桁目とする）

ですね。各桁の数字の配列を元の数字に復元しています。適当な数で試してみればわかると思います。

なぜわざわざいったん各桁の数字の配列を間にはさんでいるのかというと、10進数以外にも対応するためです。

20行目からのforループがちょっと面倒なのですが、結局はfirstよりあとの部分（上の図ではexp1、exp2、…）をfirstと同様に処理しています。
24〜26行目で（13〜16行目と同様に）exp[k]を元の数字に復元します。
28、29行目ではさらに外部の関数を呼んでいますが、名前の通りmpn_mul_1はかけ算（rp*=info->bb）、mpn_add_1は足し算（rp+=w）です。
興味がある人は読んでみるといいと思います。両方とも、人間が手で筆算をするのと同様のアルゴリズムで計算されています。キャリーオーバー（繰り上がり）が綺麗に処理されててなかなか面白いです。

結局、この部分はrp = rp * info->bb + wとなります。つまり、16行目や26行目と同じ計算です。
これを繰り返して任意精度整数のいちばん下の位まで計算すると、2^64を基数とした2^64進数に変換できています。

f:id:levelfour:20140122223405p:plain

たとえば、最初に挙げたサンプルで代入した数字48410242354393281104234213124421033について、gdbでguフォーマット（ようするにunsigned longです）でダンプしてみると15153862479221437865と2624324496559148という2つのデータが見えました。
これは
15153862479221437865 * (2^64)^0 + 2624324496559148 * (2^64)^1 == 48410242354393281104234213124421033
ということです！（だいたいのスクリプト言語は任意精度整数をサポートしているので計算して確かめてみるといいかと思います）