1 基本的なアルゴリズム

1.1 多倍長演算とは

CPUが直接扱える数の範囲には限りがあり、昔からのプログラミング言語ではその制限が残っていることが多い。C言語を例にとってみると、int型の変数は通常CPUのレジスタのサイズとなるので、32bitのCPUでは $2^{32}=4294967296$種類の数しか表現することができないことが多い。それに対して、Lisp系の言語や最近のスクリプト言語では任意の桁数の数を扱えるものがよくある。これらの言語は、限られた桁数の数のみによる任意の桁数の数の計算、つまり多倍長演算(Multiple Precision Arithmetic)を行うプログラムを作ってCPUの制限を解消しているのである。

1.2 多倍長数の表現と加減乗算

多倍長演算の対象となる数、つまりCPUが直接扱える範囲を超えるような数を多倍長数(Multiple Precision Number)と呼ぶことにする。多倍長演算を行うにはまず多倍長数を表現しなければならないが、大きな数を日常的に3141592と一桁の数の並びで表すように、一桁ごとに変数を割り当てるのが最も単純な方法である。つまり

$$u_6=3,u_5=1,u_4=4, \ldots ,u_0=2$$

のようにすればよいのである。実装する際には10進数よりも$2^{32}$進数などで表現した方が、桁上げがビット演算で行えるのでよい。

多倍長数をこのように表現しておけば、加減乗算は筆算(Schoolbook Method)を使うことで簡単に実現できるので説明は必要ないだろう。$n$桁同士の加減乗算に必要な時間は、それぞれ $\Theta(n),\Theta(n),\Theta(n^2)$である。除算については次の節で詳しく見ることにする。

1.3 除算

実際にやってみれば分かるように、筆算では$(n+1)$桁と$n$桁の除算が重要となってくる。そこで$u,v$

$$u=(u_{n}u_{n-1}u_{n-2} \ldots u_0)_b,v=(v_{n-1}v_{n-2} \ldots v_0)_b$$

に対して$u=q_0v+r$ $(q_0<b,r<v)$となる場合の除算を考える。$q_0$ の満たす必要条件は

$$q_0v \le u \nonumber \Longrightarrow q_0v_{n-1}b^{n-1} \le u... ...\le \biggl\lfloor \frac{u_{n}b+u_{n-1}}{v_{n-1}} \biggr\rfloor$$

となる。よって、$q_0 \le b-1$ であることを考えて

$$q_0=\min\left(\biggl\lfloor \frac{u_{n}b+u_{n-1}}{v_{n-1}} \biggr\rfloor,b-1 \right),r=u-qv$$ (1.1)

とおき、$r \ge 0$ となるまで

$$q_0:=q_0-1,r:=r+v$$ (1.2)

を繰り返せば $q_0,r$ が求まる。

この繰り返しの回数は、

$$q_0'=\min\left( \biggl\lfloor \frac{u_{n}b+u_{n-1}}{v_{n-1}} \biggr\rfloor,b-1 \right)$$

とおけば$q_0'-q_0$で与えられる。ここで

$$q_0'\le\frac{u_{n}b+u_{n-1}}{v_{n-1}}\le\frac{u}{v_{n-1}b^{n-1}},q_0>\frac{u}{v}-1$$

であるから、

$$q_0'-q_0 < \frac{u}{v_{n-1}b^{n-1}}-\frac{u}{v}+1 =\frac{u}{v} \cdot \frac{v-v_{n-1}b^{n-1}}{v_{n-1}b^{n-1}}+1$$

$$\le \frac{u}{v} \cdot \frac{1}{v_{n-1}}+1$$

$$< \frac{b}{v_{n-1}}+1$$ (1.3)

となる。よって、 $v_{n-1}\ge b/2$ならば繰り返しは必ず2回以下となるのである。この条件を満たすためには適当な数$k$を$v$に掛ければ良い。$u$にも$k$を掛けて除算を行えば、商は$q$、余りは$kr$となるので、$kr$を$k$で割れば求める結果が得られる。$b$が2の冪の場合は$k$も2の冪となるので、シフト演算を用いることでより高速に行える。$2n$桁と$n$桁の除算の実行時間は$\Theta(n^2)$である。

1.4 平方根

平方根の計算方法を調べるために、 $u=(u_{2n-1}u_{2n-2} \ldots u_0)_b$の平方根 $s=(s_{n-1}s_{n-2} \ldots s_0)_b$を求めることを考える。まず、$s_{n-1}$は

$$(s_{n-1}00 \ldots 0)_b^2 \le u \Longleftrightarrow s_{n-1}^2 \le u_{2n-1}b+u_{2n-2}$$

を満たす最大の整数であるから、

$$s_{n-1}=\Bigl\lfloor \sqrt{u_{2n-1}b+u_{2n-2}} \Bigr\rfloor$$ (1.4)

と計算できる。ここで誤差(余り) $r_{n-1}$を

$$r_{n-1}=(u_{2n-1}b+u_{2n-2})-s_{n-1}^2$$ (1.5)

とおく。次に、$s_{n-2}$は

$$\begin{eqnarray*} && (s_{n-1}s_{n-2}0 \ldots 0)_b^2 \le u \\ &\Longleftrightarr... ...w& (2s_{n-1}b+s_{n-2})s_{n-2} \le r_{n-1}b^2+u_{2n-3}b+u_{2n-4} \end{eqnarray*}$$

を満たす最大の整数であり、誤差$r_{n-2}$は右辺から左辺を引いた

$$r_{n-2}=r_{n-1}b^2+u_{2n-3}b+u_{2n-4}-(2s_{n-1}b+s_{n-2})s_{n-2}$$

となる。同様にして、一般の$s_k$は

$$(2(s_{n-1}s_{n-2} \ldots s_{k+1})_bb+s_k)s_k \le r_{k+1}b^2+u_{2k+1}b+u_{2k}$$ (1.6)

を満たす最大の整数で、誤差$r_k$は

$$r_k=r_{k+1}b^2+u_{2k+1}b+u_{2k}-(2(s_{n-1}s_{n-2} \ldots s_{k+1})_bb+s_k)s_k$$ (1.7)

となることが分かる。

$s_k,r_k$を求める方法を考えよう。(1.7)から

$$(2(s_{n-1}s_{n-2} \ldots s_{k+1})_bb+s_k)s_k \le r_{k+1}b^2+u_{2k+1}b+u_{2k}$$

$$\Longrightarrow 2(s_{n-1}s_{n-2} \ldots s_{k+1})_bbs_k \le r_{k+1}b^2+u_{2k+1}b+u_{2k}$$

$$\Longleftrightarrow 2(s_{n-1}s_{n-2} \ldots s_{k+1})_bs_k \le r_{k+1}b+u_{2k+1}$$

$$\Longleftrightarrow s_k\le\frac{r_{k+1}b+u_{2k+1}}{2(s_{n-1}s_{n-2} \ldots s_{k+1})_b}$$ (1.8)

となる。最後の除算の商を$s_k$、余りを$r_k'$とおけば、$r_k$は

$$r_k = r_{k+1}b^2+u_{2k+1}b+u_{2k}-(2(s_{n-1}s_{n-2} \ldots s_{k+1})_bb+s_k)s_k$$

$$= r_k'b+u_{2k}-s_k^2$$ (1.9)

となる。$s_k$が真値よりも大きければ$r_k$は負となるので、$r_k$の値が正となるまで

$$r_k:=r_k+2s_k-1,s_k:=s_k-1$$ (1.10)

を繰り返せば良い。

さて、除算のときと同じように繰り返しの回数を検討しよう。$k=n-1$ の場合は単純に平方根をとるだけであるから、$k \le n-2$ の場合を考える。

$$s_k'=\biggl\lfloor\frac{r_{k+1}b+u_{2k+1}}{2(s_{n-1}s_{n-2} \ldots s_{k+1})_b} \biggr\rfloor$$

とおけば、繰り返しの回数は $s_k^{\prime}-s_k$で与えられる。

$$\begin{eqnarray*} s_k' &\le& \frac{r_{k+1}b+u_{2k+1}}{2(s_{n-1}s_{n-2} \ldots s_... ...2+u_{2k+1}b+u_{2k}}{2(s_{n-1}s_{n-2} \ldots s_{k+1}0)_b+s_k+1}-1 \end{eqnarray*}$$

であるから、

$$s_k^{\prime}-s_k < \frac{r_{k+1}b^2+u_{2k+1}b+u_{2k}}{2(s_{n-1}s_{n-2} \ldots s_{k+1... ...frac{r_{k+1}b^2+u_{2k+1}b+u_{2k}}{2(s_{n-1}s_{n-2} \ldots s_{k+1}0)_b+s_k+1} +1$$

$$= \frac{r_{k+1}b^2+u_{2k+1}b+u_{2k}}{2(s_{n-1}s_{n-2} \ldots s_{k+1}0)_b+s_k+1} \cdot \frac{s_k+1}{2(s_{n-1}s_{n-2} \ldots s_{k+1}0)_b} +1$$

$$< \frac{b^2}{2(s_{n-1}s_{n-2} \ldots s_{k+1}0)_b}+1$$

$$\le \frac{b^{k+3-n}}{2s_{n-1}}+1$$ (1.11)

となる。この式から、繰り返しは$k<n-2$の場合常に1回以下であり、また $s_{n-1} \ge b/2$ならば$k=n-2$の場合も1回以下となることが分かる。よって $s_{n-1} \ge b/2$とするために、$u$に適当な数$k^2$を掛けて $u_{2n-1} \ge b/4$とすると、平方根が$ks$、余りが$k^2r$となるので、それぞれ$k,k^2$で割って補正すれば良い。実行時間は$\Theta(n^2)$である。

1.5 基数変換

基数を10の冪以外にしている場合には、数を表示する際に基数変換が必要となることが多い。

自然数$n$を基数$b$で表すには、$n$を$b$で割った商と余りを求め、商を新たに $n$ と置く、という操作を繰り返せば良い。出てきた余りの列は、$n$を$b$ 進表記したときの下位の桁から順に並べたものになる。

次に整数部が$0$の正の実数$r$の基数変換を考える。今度は$r$を$b$倍して整数部と小数部に分け、小数部を新たに$r$と置く、という操作を繰り返せば良い。出てきた整数部の列は、$r$の小数部を上位の桁から順に並べたものになる。

一般の実数は整数部と小数部に分け、それぞれを上に述べた方法で基数変換すれば良い。実行時間は$\Theta(n^2)$である。