基底変換
多倍長計算をする上で重要だけれど意外と忘れがちなのが基底変換である。 計算そのものを 10 進数で行うのであれば特に必要がないが、 超高速を実現しようとなるとコンピュータが使う 2 進数で計算するのが必然となり 人間が見る形にするため 10 進数へ変換する必要がある。 ここでは 2 進数 (正確には $2^{32}$ 進数や $2^{64}$ 進数) で計算がなされる環境を前提とし、 具体的な対象として 10 進数 (正確には $10^9$ 進数や $10^{19}$ 進数）を例に挙げながら 基底変換のアルゴリズムを紹介する。
整数の変換
まずは基本的に変換を正確に行える(非負)整数の変換を考える。 問題を数学的にもう少し厳密に定義しよう。 つまり今我々が求めたいのは整数 $x$ と基数 $b$ に対して
\[ x = s_{n-1}b^{n-1} + s_{n-2}b^{n-2} + \cdots + s_1 b + s_0 \]を満たす整数列 $\{s_i | 0 \leqq s_i \lt b\}$ を求めることである。 今回は $n$ が分かっている前提で話を進めるが、 未知の場合は $n=\lceil\log_{10} x\rceil$ で求めることができる。
基本アルゴリズム
上の式を見て分かる単純なアルゴリズムは、$x$ を $b$ で割っていく方法である。
\[ x = (s_{n-1}b^{n-2} + s_{n-2}b^{n-3} + \cdots + s_2b + s_1) b + s_0 \]となっているので $x \bmod {b}$ が $s_0$、$\lfloor x/b\rfloor$ が 2 桁目以上の数を表すということが分かる。 これを $x := \lfloor x / b \rfloor$ として $n$ 桁分繰り返せば変換後の数字列 $\{s_i\}$ が求まる。
この方法では下の桁から求まるため、通常表記と同様に Big Endian で処理する場合は 最後に文字列の反転を行うか、処理の途中でもどの位置に出力するかを把握しつつ 処理していく必要がある。
- function ConvertIntegerBase($x$, $n$)
- $s$ = {サイズ $n$}
- for $i$ = 0 $\cdots n-1$:
- $s_i$ = $x \bmod b$
- $x$ = $\lfloor x/b \rfloor$
- return $s$
- End of function
Karatsuba 法 [名称要確認]
Binary splitting や DRM 法 と同様に 同程度の桁数で分割統治する方法を考えよう。つまり
\[ \begin{eqnarray} x &=& s_{n-1}b^{n-1} + s_{n-2}b^{n-2} + \cdots + s_1 b + s_0\\ &=& (s_{n-1}b^{n_1-1} + s_{n-2}b^{n_1-2} + \cdots + s_{n_1})b^{n_1} + (s_{n_1-1}b^{n_1-1} + \cdots + s_0)\\ &=& X_1 B_1 + X_0 \end{eqnarray} \]という風に分割することを考える。 このとき $B_1=b^{n_1}$ は 2 の冪乗ではないので $B_1$、$X_1$、$X_0$ は全て計算して求める必要がある。
\[ \begin{cases} B_1 = b^{n_1}\\ X_1 = \lfloor x / B_1 \rfloor\\ X_0 = x \bmod B_1 \end{cases} \]そして$n_k$ が適当な大きさになるまでこの方法を $X_0$, $X_1$ に再帰的に適用し、 それ以降は基本アルゴリズムを適用していくことで $b$ 進数の文字列にすることができる。 計算効率や再利用性の観点から、$B_i$ は $b$ の 2 の冪乗乗、つまり $B_i=b^{2^i}$ とする方法が知られているが、負荷分散を考えると分割統治方式で定めることも良いかもしれない。[要実験]
- function ConvertIntegerKaratsuba($x$, $n$)
- $N$ = ($n$ 以下で最大の 2 の冪乗)
- return KaratsubaCore($x$, $N$)
- End of function
- .
- function KaratsubaCore($x$, $n$)
- if $n \lt n_t$:
- return ConvertIntegerBase($x$, $n$)
- $B$ = $b^N$ (先に作ってキャッシュしておく)
- $x_0$ = $x \bmod B$
- $x_1$ = $\lfloor x/B \rfloor$
- return KaratsubaCore($x_1$, $N/2$) + KaratsubaCore($x_0$, $N/2$)
- if $n \lt n_t$:
- End of function
Scaled Remainder Tree 法
Karatsuba 法のアルゴリズムを見てみると、かなりの数の多倍長除算が必要なことがわかる。 これを減らしたアルゴリズムが Scaled Remainder Tree 法である。 このアルゴリズムの要点は浮動小数の変換を行う形で考える方が分かりやすいので そちらで述べる。 整数の変換では $x$ が $b$ 進数で高々 $n$ 桁であることを前提に
\[ y = x / b^n \]として $0\leqq y \lt 1$ を変換していくことになる。
- function ConvertIntegerSRT($x$, $n$)
- $s$ = SRTCore($2^mx/b^n$, $n$, $m$) # $m$ は $x$ のビット数
- return $s$
- End of function
浮動小数点数の変換
今度は浮動小数点数の基底を変換する方法を紹介しよう。ここに言及している資料はあまり見られないのが残念である。
前提条件として、16 進数で管理されている多倍長浮動小数点数 $x$ は $0\leqq x \lt 1$ とする。 $\pi$ のように $x \geqq 1$ である数を変換したい場合は、(1) 整数部を引く (2)整数部と小数部を別処理にする (3)$b$ の冪乗で割って $x \geqq 1$ に落とす などの方法で条件に合うようにすることができる。 また、入力となる $x$ は必ずしも $b$ 進数でキリの良い数であるとは限らない。そのため変換後の数 $y$ を
\[ \begin{eqnarray} y &=& s_{n-1}b^{-1} + s_{n_2}b^{-2} + \cdots + s_{0}b^{-n}\\ &=& (s_{n-1}b^{n-1} + s_{n_2}b^{n-2} + \cdots + s_{0})b^{-n}\\ \end{eqnarray} \]とするとき、$|x-y| \lt b^{-n}$ となる $y$ ならびに $\{s_i\}$ を求めることとしよう。
整数化
多倍長浮動小数点数の基底変換を一番素直に行う方法は一旦整数化して基底変換し、 小数点を移動させる方法である。 つまり $x$ を $\lfloor b^kx\rfloor$ にして変換し、先頭に "0." を置くことになる。 この方法では整数でのアルゴリズムを流用することができるため実装を手軽に行うことができる。
- function ConvertFloatByInt$(x, n)$:
- $x$ = $\lfloor b^nx \rfloor$
- $s$ = ConvertIntegerKaratsuba$(x, n)$
- return "0." + $s$
- End of function
基本変換
浮動小数点数の変換では $b$ をかけて整数部分を取り出すことを繰り返すと 整数での変換とは逆に上位の桁からの変換を簡単に行うことができる。 [FB03]
- function ConvertFloatBase($x$, $n$):
- $s$ = FloatCore($x$, $n$)
- return "0." + $s$
- End of function
- .
- function FloatCore($x$, $n$):
- for $i$ in $(n, 0]$:
- $x = bx$
- $s_{n-1-i} = \lfloor x \rfloor$
- $x = x \bmod 1$
- while end
- return $s$
- for $i$ in $(n, 0]$:
- End of function
Scaled Remainder Tree 法
この方法も同様にを半分ずつの桁数で分割統治することはできるだろうか。 具体的な数字を見ながら考えていこう。
\[ x = 0.243f6a8885a308d31319_{[16]} \]を 10 進数に変えていきたいとする。 ちなみにこの数は $\pi$ の小数点以下 20 桁分であり、10 進数では $20 \log_{10}{16} = 24.08$ から 24 桁分に相当するので 24 桁の 10 進数に変換しよう。 まずはこれを 10 進数で半分ずつの桁数に分け、同時に有効桁数も半分にする。 正確には変換後の桁数で数えて $n_0+n_1$ 桁を上位 $n_1$ 桁(実質 $n_1-1$ 桁) と下位 $n_0+1$ 桁に分割する。
\[ \begin{cases} x &=& x_1 + 10^{-12}x_0 \\ x_1 &=& 0.243f6a8885_{[16]} &(= 0.14159265358)\\ x_0 &=& 0.fab4f7b3e5e_{[16]} &(= 0.9793238462642) \end{cases} \]上半分($x_1$)を求める計算手順としては、単純に有効桁数を落とせばよいだけである。 今回は 12-12 桁に分けるので 10 進数 12 桁相当が残る 16 進数 10 桁を残す。 そして下半分($x_0$) は 10 進数 12 桁目からを担当するので $10^{11}x \bmod 1$ を計算する。その際に元々有効桁数が 16 進数 20 桁であったが 10 進数 13 桁分をカバーする 16 進数 11 桁に丸める。
同様に次は (16 進数, 10進数) で (5,6),(6,7),(6,6),(6,7) 桁ごとに分ける。
\[ \begin{cases} x_{11} &=& 0.243f6_{[16]} &(=0.14159)\\ x_{10} &=& 0.43ee8f_{[16]} &(=0.265358)\\ x_{01} &=& 0.fab4f7_{[16]} &(=0.97932)\\ x_{00} &=& 0.6276e0_{[16]} &(=0.3846263) \end{cases} \]ある程度有効桁数が小さくなったらあとは基本変換を適用することで各桁の数字を得ることができる。
この考え方をもとに整数演算で計算するアルゴリズムが [FT15] に以下のようにまとめられている。 この中では 16 進数ではなく 2 進数への変換としている他、有効桁数の管理を容易にするため演算自体は整数に変換して行われている。
- function ConvertFloatSRT$(x, n)$:
- $N =$ 基本変換を使う桁数の閾値(固定)
- $g = \max(\lceil \log_2n\rceil + 1, N)$
- $k = \lceil \log_2 (4gb^n) \rceil$ (2進数で使うビット数)
- if $x$ is an integer:
- $x$ = $\lfloor (x+1)2^k/b^n \rfloor - 1$
- return RecursiveSRT$(x, n, k, g)$
- End of function
- .
- function RecursiveSRT$(x, n, k, g)$:
- if $n \leqq N$
- return FloatCore$(x, n)$
- $n_1$ = $\lfloor (n + 1) / 2 \rfloor$
- $n_0$ = $n - n_1 + 1$
- $k_1$ = $\lceil \log_2 (4gb^{n_1}) \rceil$
- $k_0$ = $\lceil \log_2 (4gb^{n_0}) \rceil$
- $x_1$ = $\lfloor x/2^{k-k_1}\rfloor$
- $x_0$ = $\lfloor (b^{n-n_0}x \bmod 2^k) / 2^{k-k_0}\rfloor$
- $s$ = RecursiveSRT$(x_1, n_1, k_1, g)$
- $t$ = RecursiveSRT$(x_0, n_0, k_0, g)$
- if $s$ の一番下の数字が $b-1$ && $t$ の最上位の数字が 0:
- $s = s + 1 \bmod b^{n_1}$
- if $s$ の一番下の数字が 0 && $t$ の最上位の数字が $b - 1$:
- $t = 000 \cdots 000$ ($n_0$ 桁の 0)
- return $\lfloor s/b \rfloor + t$ (上位の最下位桁を除いて、下位をくっつける)
- if $n \leqq N$
- End of function