円周率を計算するのに使われるChudnovskyの公式やそのプログラムについて色々書いていますが、多倍長計算ライブラリGMPが配布しているプログラム gmp-chudnovsky.c が爆速な割りに解説されているのを見たことがないなと思って書いてみました。
概要
まずは計算に使われている公式はChudnovskyの公式
\frac{1}{\pi}=\frac{12}{\sqrt{C}^3}\sum_{k=0}^\infty
\frac{(-1)^k(6k)!(A+Bk)}{(3k)!(k!)^3C^{3k}}
(A,B,C) =(13591409, 545140134,640320)
です。これに Binary Splitting 法を適用して計算していくことになります。
Binary Splitting 法
Binary Splitting 法 (BS、BSM と略す) は私の過去の記事では類似のDRMの説明で置き換えていましたが、ここではプログラム内の変数説明が楽なので元の定義に沿って説明を書いておきます。
\frac{q(0,n)}{p(0,n)} = \sum_{k=0}^{n-1}\frac{a_k}{b_k}
\prod_{j=0}^k\frac{ q_j}{p_j}
という形の計算式において、適切に $p(i,j)$、$q(i,j)$、$g(i,j)$ を定義すると $a \leq m \leq b$ に対して
\begin{cases}
g(a,b) = g(a,m)g(m,b)\\
p(a,b) = p(a,m)p(m,b)\\
q(a,b) = q(a,m)p(m,b) + q(m,b)g(a,m)
\end{cases}
が成り立ち、分割統治で計算することができます。しかも桁数の大きな数同士の乗算についてはFFTを利用することで高速に計算できるので、$m$をバランスよく取れれば全体として $O(n(\log n)^3)$ の計算量で計算することができます。
Chudnovsky の公式を利用する場合の具体的な値は
\begin{cases}
g(k-1,k) = (6k-5)(2k-1)(6k-1)\\
p(k-1,k) = k^3 C^3 / 24\\
q(k-1,k) = (-1)^kg(k-1,k)(A+Bk)
\end{cases}
となっています。(Code)
正確には $p(0,0)$ が本来なら特殊な値となるためプログラム的には別途処理する必要がありますが、gmp-chudnovsky.c では $m=1$ で最初に分割して例外化処理を省き、最後に
\frac{q(0,n)}{p(0,n)} = \frac{q(1,n)+Ap(1,n)}{p(1,n)}
という形で計算しています。(Code)
高速化の工夫
さて、BS法を使うことで計算量としては最適になりましたが、gmp-chudnovsky.cでは他の工夫を施すことで実時間での高速化を図っているので紹介します。ここからが本番。
多倍長数のメモリ領域のプール
BS 法は分割統治方式なので、ほぼ必然的に再帰呼び出しを行うことになります。しかも2つに分けた領域からそれぞれ多倍長整数を3つずつもらう形になります。この6個の多倍長整数のそれぞれに対して毎回 mpz_init() をかけてメモリ割り当てを行っていては遅くなってしまいます。
そこで gmp-chudnovsky.c では pstack、qstack、gstack というスタックを作り、その上2つずつを使ってメモリ領域の使い回しをしています。スタックの高さは top という変数を使って管理し、スタックの参照を簡単にするために p1, p2, q1, q2 といった alias を定義しています。
static mpz_t *pstack, *qstack, *gstack;
static fac_t *fpstack, *fgstack;
static long int top = 0;
#define p1 (pstack[top])
#define q1 (qstack[top])
#define g1 (gstack[top])
#define fp1 (fpstack[top])
#define fg1 (fgstack[top])
#define p2 (pstack[top+1])
#define q2 (qstack[top+1])
#define g2 (gstack[top+1])
#define fp2 (fpstack[top+1])
#define fg2 (fgstack[top+1])
不要な演算の省略
再帰計算の最後では $p(0,n)$、$q(0,n)$ だけ必要になるので $g(0,n)$ は使いません。するとその再帰計算を一段分見ると $g(m,n)$ の計算が必要ないということが分かります。つまり再帰呼び出しにおいて必要のない $g(m,b)$ の計算を省いて高速化することができます。(Code)
if (gflag) {
mpz_mul(g1, g1, g2);
fac_mul(fg1, fg2);
}
約分と因数分解
BS の数式だけではわかりにくいですが、分割統治の途中で $p(m,b)$ と $g(a,m)$ は約分できます。数式で追っていきたい場合は約分についてコメントしているわけではないですが こちらの解説を御覧ください。
そして約分するにあたって、計算対象になる $p$ と $g$ は因数分解した形 fac_t (Code) でも値を保持しています。
typedef struct {
unsigned long max_facs;
unsigned long num_facs;
unsigned long *fac;
unsigned long *pow;
} fac_t[1];
この max_facs は保存できる素因数の個数の最大値、num_facs が実際に保存している素因数の個数を示しています。この 2 つがバラバラに管理されているのはデータ更新時にメモリアロケーション(実は結構重い処理)を頻発させないためです。
2 つの fac_t 型変数を約分をするにあたっては、ざっくり以下のようなことをしています。(Code)
- それぞれが持つ素因数を小さい方から順に見ていきます
- もしどちらかの素因数が全て見終わっていれば終わります
- もし両方で見ている素因数が同じであれば
- 指数が大きい方は小さい方の指数だけ小さくします
- 指数が小さい方はその素因数を消します
- 見ている素因数が違う場合は、素因数が小さい方だけ次の素因数を見ます
for (i=j=k=0; i<fp->num_facs && j<fg->num_facs; ) {
if (fp->fac[i] == fg->fac[j]) {
c = min(fp->pow[i], fg->pow[j]);
fp->pow[i] -= c;
fg->pow[j] -= c;
fmul->fac[k] = fp->fac[i];
fmul->pow[k] = c;
i++; j++; k++;
} else if (fp->fac[i] < fg->fac[j]) {
i++;
} else {
j++;
}
}
プログラム内では同時に mpz_t 型でも約分を行う必要があるのですが、またこの手順も少し凝った形になっています。まず上記の fac_t 型の約分と同時並行で fac_t 型の最大公約数を求め、分割統治方式で mpz_t 型に変換します。そして mpz_divexact 関数を使って最大公約数を除きます。mpz_divexact 関数は割り切れることが分かっている前提で高速に割り算ができる関数です。
if (fmul->num_facs) {
bs_mul(gcd, 0, fmul->num_facs);
mpz_divexact(p, p, gcd);
mpz_divexact(g, g, gcd);
}
篩
次に気になるのは素因数分解の方法ですが、これは一部界隈で osa_k法 と呼ばれる エラトステネスの篩を拡張した方法、にさらに改良を加えた方法です。
まずは基本となるエラトステネスの篩をざっとおさらいします。Chudnovskyの公式での利用を考えて偶数の対処は考えていません。
- $1..n$ の数 $i$ について $s[i] = $True とします
- $i$ を 3 から $n$ まで小さい順に回します
- もし $s[i] = $False であれば次の $i$ に移ります
- $k$ を $i^2$ から $n$ まで $2i$ ずつ増やして
- $s[i] = $False とします
これでそれぞれの数が「素数なのかそうでないのか」がわかるので、素因数する際には素数を小さい方から割れるかどうか確認していくことになります。エラトステネスの篩は素数を見つけていく手順としては高速ですが、その結果を残していたとしても素因数分解の効率化は見込めません。
次に同じような篩の処理を「どの素数で割れるのか」という情報を残しつつ行うことになる osa_k 法を紹介します。
- $1..n$ の数 $i$ について $s[i] = i$ とします
- $i$ を 3 から $n$ まで小さい順に回します
- もし $s[i] \neq i$ であれば次の $i$ に移ります
- $k$ を $i^2$ から $n$ まで $2i$ ずつ増やして
- $s[k] = k$ であれば $s[k] = i$ とします
上と見比べて分かる通り、処理としては保存する値を True/False から素因数に変更しただけです。しかし、この新しい保存データが有ることで、範囲内の任意の奇数 $m$ が $s[m]$ を(最小の)素因数として持つことがわかるので $m$ が 1 になるまで $s_m$ で割っていく形で簡単に素因数分解を行うことができます。ちなみに最後の $s[k]=k$ の確認をしないで $s[k]$ を上書きする場合、記録された数は最大の素因数になります。
そして gmp-chudnovsky.c ではさらに指数も保持することで素因数分解の方の処理も高速化します。(Code)
なお $p[i]$ が Code 内の s[i/2].fac に対応し、$i$ の最小素因数を表しています。あとは $e[i]$ が s[i/2].pow、$f[i]$ が s[i/2].nxt に対応しています。ループ変数については $i$ が素数、$j$ は $i$ を素因数に持つ合成数、$k$ は $i/j$ の値、に相当します。
- $1..n$ の数 $i$ について $(p[i], e[i], f[i]) = (0,0,1)$ とします
- $i$ を 3 から $n$ まで小さい順に回します
- もし $p[i] \neq 0$ であれば次の $i$ に移ります
- $(p[i], e[i]) = (i, 1)$ とします
- $j$ を $i^2$ から $n$ まで $2i$ ずつ、$k$ を $i$ から 2 ずつ増やして
- もし $p[j] \neq 0$ であれば次の $(j,k)$ に移ります
- もし $p[k] = i$ であれば
- $(p[j], e[j], f[j]) = (i, e[k]+1, e[k])$ とします
- そうでなければ ($f[k] \neq i$ であれば)
- $(p[j], e[j], f[j]) = (i, 1, k)$ とします
static void build_sieve(long int n, sieve_t *s) {
long int m = (long int)sqrt(n);
memset(s, 0, sizeof(sieve_t) * n / 2);
s[1/2].fac = 1;
s[1/2].pow = 1;
for (i=3; i<=n; i+=2) {
if (s[i/2].fac == 0) {
s[i/2].fac = i;
s[i/2].pow = 1;
if (i<=m) {
for (j=i*i, k=i/2; j<=n; j+=i+i, k++) {
if (s[j/2].fac==0) {
s[j/2].fac = i;
if (s[k].fac == i) {
s[j/2].pow = s[k].pow + 1;
s[j/2].nxt = s[k].nxt;
} else {
s[j/2].pow = 1;
s[j/2].nxt = k;
}
}
}
}
}
}
}
すると整数 $m$ について $m = p[m]^{e[m]} \times f[m]$ となる上に $f[m]$ は $p[m]$ を因数に持たないので、$f[m]$ を再帰的に $m$ にしていくことで素因数分解を高速に行うことができます。プログラム内では $k^3$ といった指数方式での計算もあるので、それに合わせた形で実装されています。(Code)
/* f = base^pow */
static void fac_set_bp(fac_t f, unsigned long base, long int pow) {
long int i;
for (i=0, base/=2; base>0; i++, base = sieve[base].nxt) {
f[0].fac[i] = sieve[base].fac;
f[0].pow[i] = sieve[base].pow * pow;
}
f[0].num_facs = i;
}
計測
せっかくなので、ここで解説した高速化技術を敢えて取り除いた場合にどの程度必要メモリや時間が増えるか計測してみました。1億桁をターゲットにして1回ずつの計測なのでブレはあるでしょうがとりあえず。
| オリジナル | 遅いバージョン | |
|---|---|---|
| sieve | 0.359 | - |
| bs | 45.489 | 58.662 |
| gcd | 0.000 | - |
| div | 4.518 | 4.412 |
| sqrt | 2.503 | 2.461 |
| mul | 1.635 | 1.653 |
| total | 54.539 | 67.211 |
| P size | 145605885 | 248778666 |
| Q size | 145605879 | 248778659 |
| Memory | 1.21 GB | 728 MB |
メモリ計測は CHECK_MEMUSAGE を定義して出てきた値の最大値を使ってます。この結果によるとメモリは篩のために倍増しますが、計算時間は20％程削減されているようです。
また、2009年に当時世界記録となる2.7兆桁の円周率計算を行ったBellardのレポートによると、gmp-chudnovsky.c の実装では篩に使うメモリは桁数に比例する $O(N)$ となるが、これは工夫することで $O(N^{1/2})$ で済むので(何兆桁という世界記録スケールになると)実質無視できる規模となるらしいです。