数理解析研究所講究録 1456 巻 2005 年 174-179174

指数対数関数等の超越関数の多倍精度計算
平山 弘 (Hiroshi Hirayama)
神奈川工科大学 *1

森川 敦司 (Atsushi Morikawa)
神奈川工科大学

はじめに級数の高精度計算を有理数を使って計算することによって、計算速度を上げる方法が多くの入によって使われている。 この方法を使って、 1989 年に Chudnovsky 兄弟による円周率 10 億桁の計算が行われた。 しかしながら、 彼らは、 この計算方法の詳細を発表していない。

この方法は 1998 年から 1999 年にかけて、 右田等 [3] [4] や後等 [6] によって、 級数計算への適用の再現が行われている。

この計算法の起源を調べると 1997 年にすでに B.Haible and T.Papanikolaou [1] によって公表されていることがわかる。 この論文では、 トーナメント法とか縮約法とか呼ばれているこれらの計算方法を Binary
Splitting Algorihhm (以降 BSA 法と略す) と呼んでいる。 この論文を読むと BSA 法は、 1976 年に R. P.
1987 年に Borwein 兄弟によって公表され、 Chudnovsky 兄弟による円周率 10 億桁の計算が行われた 1989 年当時でもよく知られた方法であったことがわかる。

$\mathrm{B}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}_{\text{、}}$

BSA 法は、 分割統治法の一種で、 大きな問題を小さな問題に分割し高速化を図るものである。 別の言い方をすれば、 本方式は分割した問題をトーナメント形式に計算することで高速化を図るものである。

高精度数の乗算については高速フーリエ変換 (FFT) を使う高速な計算方法がよく知られている。 FFTは約 1000 桁以上の数値の乗算において有効である。 BSA 法は 1000 桁以下の数値の計算でも高速に計算できる上、 原理が単純なので、様々な計算において容易に利用することが出来る。

本論文では、 指数対数関数の高精度計算として Taylor 展開に BSA 法を使って高速化する方法提案する。
約 1000 桁以下の精度の計算では、 Taylor 展開を使った計算が Sasaki and Kanada[5] によって、様々な計算法を比較して最も高速であることが示されているので、 計算時間が問題となるのは、 1000 桁以上の精度の計算である。 ここで提案した Taylor 展開に BSA 法を適用して高速化した方法と Sasaki and Kanda によって提案された方法を 1000 桁を超えた精度で比較し、 その高速性を示した。
以下の計算では、 高精度計算プログラムとして、 C++言語で自作したものを利用した。 10000 進数表現の浮動小数点演算ルーチンである。 1600 桁を超えると自動的に FFT を利用して計算する。

2

BSA 法
BSA 法は FFT と同様な分割統治法の一種で、大きな問題を小さな問題に分割して計算する方法である。

*hirayarna@sd kanagawa-it.ac$.\mathrm{j}\mathrm{p}$

175
21

簡単な例

B.Haible and T.Papanikolaou の論文にもあるが、 階乗の計算が非常に良い例である。次の計算を行う。
(1)$n!= \prod_{k=1}^{n}k=1\cdot 2\cdot 3\cdots n$が小さい場合、前から順番に計算しても計算時問に大きな違いが生じることはないが、前から順序よく計算する方法は効率的ではない。 が偶数ならば、 2 個つつ組み合わせ$n$$n$が大きい場合、

$n$

(2)

$n!=(1\cdot 2)\cdot(3\cdot 4)\cdots((n-1)\cdot n)$$n$が 4 の倍数ならば、 さらに 2 個づっ組み合わせ$n!=\{(1\cdot 2)\cdot(3\cdot 4)\}\cdots\{[(n-3)\cdot(n-2)_{\rfloor}\urcorner.

$(3\rangle$[(n-1)\cdot n_{\rfloor}^{\rceil}\}$計算する。 このように計算することで、 通常の計算方法よりも計算桁数が小さいままで計算を行うことが出来るので、計算速度の向上が期待できる。 これは 2 のべき乗の場合に限らず、 一般の の場合においても$\mathrm{n}$同様のことがいえる。

211

階乗計算例

10000! の計算を行う。 この計算では、 BSA 法を使うだけでなく、 1600 桁以上の数値に対しては FFT を利用して乗算を行っている。計算結果は$10000!=2.\mathrm{S}4625968091705451890641321211$ $9\mathrm{S}68890\mathrm{x}10^{35659}$(4)35660 桁の数値となる。 この計算には厳密に計算するため 5 万桁の浮動少数点演算) レーチンを使って計算しで行うと以下のようになる。 このとき使用したコンパイラー {よ Borland
た。計算時間は Pentium 4$3.06\mathrm{G}\mathrm{H}\mathrm{z}$$\mathrm{C}++6.0$である。
$\equiv-\{’\supset\ovalbox{\tt\small REJECT} \mathrm{p}_{\mathrm{V}}\ovalbox{\tt\small REJECT} g$$\equiv-\mathrm{p}\mathrm{t}\ovalbox{\tt\small REJECT}\hslash\backslash \mathit{1}\xi$BSA$357\mathrm{S}$$\not\in J\backslash \ovalbox{\tt\small REJECT} \mathit{0}2\mathrm{i}F\backslash l\yen$

このように、 大きな

2L2$n$(msec)47に対する階乗計算は非常に速く計算できる。

級数の計算

次の形で表現できる級数は、$S= \frac{1}{C_{0}}(A_{0}+\frac{B_{0}}{c_{1}}(A_{1}+\frac{B_{1}}{c_{2}}(A_{2}+\frac{B_{2}}{C_{3}}(A_{3}+\frac{B_{3}}{c_{4}}(A_{4}+\cdots$(5)

次のように変形することができる。
$S= \frac{1}{C_{0}}(A_{0}+\frac{B_{0}}{C_{1}C_{2}}(A_{1}C_{2}+A_{2}B_{1}+\frac{B_{1}B_{\underline{9}}}{c_{3}}\langle A_{3}+\frac{B_{3}}{c_{4}}(A_{4}+\cdots$

これらの関係を行列を使って表すと、 (5) 式の第$n$項までの値$S_{n}$(6)

は.
$B_{1}$$(\begin{array}{ll}Q_{n} P_{n}0 R_{n}\end{array})=(\begin{array}{ll}B_{0} A_{0}0 C_{0}\end{array})($$A_{1}C_{1})$0$(\begin{array}{ll}B_{2} A_{2}0 c_{2}\end{array})\cdots$$(\begin{array}{ll}B_{n} A_{n}0 c_{n}\end{array})$(7)176となる。 ここで $S_{n}=P_{n}/R_{n}$ である。 (6) 式は、右辺第 2 と第 3 の行列をかけ算したものに相当する。 この計算も、 階乗の計算と同じように、 乗算を行うことによって高速に計算することができる。
これらの計算の各数値は短い桁数の数値であるから、 非常に高速に計算できる。桁数が大きくなった場合、 高速フーリエ変換を利用した乗算法を使うと、 $O(n(\log n)^{3})$ の計算量で計算できる。

213

指数関数の計算

指数関数の Taylor 展開式から、 (5) 式に相当する式は$e^{x}= \frac{1}{1}(1+\frac{x}{1}(1+\frac{x}{2}(1+.\frac{x}{3}(1+\frac{x}{4}(1+\cdots(1+\frac{x}{n}(1+\cdots$(S)となる。 これを行列の積で表現すると$(\begin{array}{ll}Q_{n} P_{n}0 R_{n}\end{array})=(\begin{array}{ll}x 10 1\end{array})(\begin{array}{ll}x 10 1\end{array})(\begin{array}{ll}x 10 2\end{array})(\begin{array}{ll}x 10 3\end{array})\cdots(\begin{array}{ll}x 10 n\end{array})$(9)となる。

2.14

対数関数の計算

対数関数の場合も同様で Taylor 展開式から、 (5) 式に相当する式は$\frac{1_{0_{\mathrm{t}}^{\sigma}},(1+x)}{x}=\frac{1}{1}(1$$-$ $\frac{x}{2}(1-\frac{2x}{3}(1-\frac{3x}{4}(1-\cdots(1-\frac{(n-1)x}{n}(1-\cdots$(10)

行列の乗算形式に変形すると$(\begin{array}{ll}Q_{n} P_{n}0 R_{n}\end{array})=(\begin{array}{ll}-x 10 1\end{array})(\begin{array}{ll}-2x 1\mathrm{o} 2\end{array})(\begin{array}{ll}-3x 10 3\end{array})\cdots(\begin{array}{ll}-nx 10 n\end{array})$(11)

このほか、 三角関数、逆三角関数、双曲線関数など簡単な規則で各項の係数が表現でき、 多くの関数がこの行列の乗算形式に変形できます。 Taylor 展開の係数が簡単な規則で表現できない $\tan x$ が例外的に表現できないだけである。

2.2

連分数の計算法

級数と同様に、 連分数の計算も BSA 法によって高速化することが出来る。 (有限項で打ち切った連分数を次のように定義する。

.

$\frac{P_{n}}{Q_{\Gamma\iota}}.=b_{0}+\frac{a_{1}}{b_{\grave{[perp]}}+\frac{a}{b_{2}+\underline{\mathrm{s}}}}\underline,=b_{0}+\frac{a_{1}}{b_{1}}+\frac{a_{2}}{b_{2}}+\cdot$..

$+ \frac{a_{n}}{b_{n}}$

(12)

+器

このとき、 よく知られているように以下の漸化式が成り立つ。
$P_{n}=b_{n}P_{n-1}+a_{n}P_{n-2}$

Qn=bnQ ユー l+anQn-2$(n=1,2, \cdots)$

(13)

177
だだし、 $Q\mathrm{o}=1_{7}P_{0}=b_{0},$ $P_{-1}=1,$ $Q_{-1}=0$ とする。 この連分数の漸化式は次のような行列の乗算の形式で書くことができる。

(14)$(\begin{array}{ll}P_{n} P_{n-1}Q_{n} Q_{n-\mathrm{l}}\end{array})=(\begin{array}{ll}P_{n-\mathrm{l}} P_{n-2}Q_{n-1} Q_{n-2}\end{array})(\begin{array}{ll}b_{n} 1a_{n} 0\end{array})$

この関係式を利用すると、 (12) 式は、( QPnn 二 l )$Q_{n}P_{n}$$=(\begin{array}{ll}b_{0} 1a\mathrm{o} 0\end{array})(\begin{array}{ll}b_{1} 1a_{\mathrm{l}} 0\end{array})(\begin{array}{ll}b_{2} \mathrm{l}a_{2} 0\end{array})\cdots(\begin{array}{lll}b_{n} 1 a_{n} 0 \prime\end{array})$(15)となる。 ここで、 $a_{0}=1$ である。 この式を使って、 Pn/Q、を計算すれば、 連分数の値を計算できる。

221

逆正接関数の計算法

逆正接関数 $\tan x^{-1}$ は次のように連分数展開できる。

$\tan^{-1}x=\frac{x}{1}+\frac{x^{2}}{3}+\frac{4x^{2}}{5}+\frac{9x^{2}}{7}+\cdot$..

$+ \frac{n^{2}x^{2}}{2n+1}+\cdot$

..

(16)

これを行列の乗算形式になおすと$(\begin{array}{ll}P_{n} P_{n-1}Q_{n} Q_{n-1}\end{array})=(\begin{array}{ll}0 11 0\end{array})(\begin{array}{ll}1 1x 0\end{array})(\begin{array}{ll}3 1x^{2} 0\end{array})(\begin{array}{ll}5 14x^{2} 0\end{array})\cdots(\begin{array}{ll}2n+1 1n^{2}x^{2} 0\end{array})$(17)となる。 $\pi=4\tan^{-1}1$ となるので $x=1$ を代入することで計算できる。 この式を使って 10 万桁計算するのに$\mathrm{A}\mathrm{t}\mathrm{h}\mathrm{l}\mathrm{o}\mathrm{n}800\mathrm{M}\mathrm{H}\mathrm{z}$t.a$\mathrm{n}^{-1}1$において 45 秒で計算できた。

を計算するために、 $\tan^{-1}x$ の Taylor 展開式に $x=1$ を代入することは理論的には収束するが、非常に遅 \prec 、実用的には収束しない級数のため、 このような方法でこの式の値を計算することには無理がある。 しかし、 $\tan^{-1}x$ に対応する連分数においては、 実用に耐えうる速度で収束するため、 このように計算すること可能である。

23

指数関数の高精度計算指数関数は次の方法で計算した。
(18)$e^{x}=(e^{\frac{x}{\prime \mathrm{A}}}\underline’.\cdot)^{2^{n}}$この式は、$t= \frac{x}{2^{\mathrm{n}}}$と置くと、次の 2 式に分割できる。
$e^{t}=1+t+ \frac{t^{2}}{2!}+\frac{t^{3}}{3!}+\cdots$(19)$e^{x}=(e^{t})^{2^{n}}$$(20\rangle$とに分割できる。 (19) と (20) の計算量の総和が小さくなるように を選ぶ。 このとき、 は計算精度を 10 進での桁数 $P$ とした時 $n=2.25\sqrt{P}$ となる。 ただし 225 の数値は数値実験によって決定した値である。
$n$$n$

178
23I

BSA 法の利用一を高速に計算するために、 次のように変形する。

(21)$e^{t}=e^{t_{0}+t_{1}}=e^{t_{0}}e^{t_{1}}$

ここで$t_{0}$に小さい桁の数をとり、いため BSA 法で高速化できる。

は to の桁数が少なにはその先の数値にとる。 このように取ると、がきわめて小さいため急速に収束する。 このことから数値をの値は$e^{t_{\mathrm{O}}}$$t_{1}$$e^{t_{1}}$$l_{1}$

少ない上位桁と残りの桁に分割することは計算する上で非常に有効な方法となる。
を計算した。
この方法で、 5000 桁の$e^{\pi}$$e^{\pi}=23.1406926327792690057290\mathrm{S}6367948547380266106$

これまで24$530\mathrm{m}\sec$かかったものを、 BSA 法を使うことによって、 $293\mathrm{m}\sec$ に短縮できた。

対数関数の高精度計算

対数関数は次の方法で計算した。

.

(22)

$x= \log t=x0+\log(1+\frac{t-t_{0}}{t_{0}}),x_{0}=\log t_{0}$$-rightarrow-$

ここで$t_{0}$

として小さい桁の近{$\mathrm{L}^{\backslash }\mathrm{J}\not\in\llcorner$をとる。

$\log t_{0}$は’$\backslash$$t_{0}$

to の桁数力$>^{\grave{\neg}}/\mathrm{J}\grave{\prime}$$\vee$$\backslash$$\cdot$$\int$

ないため BSA 法速化できる。

$\text{て^{}\backslash }\backslash B\cap-$$\log(1+\frac{t-t_{0}}{t_{0}})$がきわめて小さいためすぐに収束する。の値は上の式を適用する前に、 BSA 法の計算を速く収束させるために、$\frac{t-}{t}\mathrm{p}t0$(23)$t= \frac{x}{2^{n}}$として、 (22) で使う値を 1 以下になるように$n$を選んだ。 このようにすると、 log2 の定数値の高精度の値が必要となる。 これを次の式を利用して、計算した。
$\log(2)=3\log(81/80)+5\log(25/24)+7\log(16/15)$$\log(81/80)_{\text{、}}1\mathrm{o}_{\mathrm{t}\supset}^{\sigma}(25/24)_{\text{、}}10_{\epsilon}^{\sigma},(16/15)$

は、 BSA 法を使って、 次の式$\log\frac{x+1}{x-1}=\log\frac{1+\frac{1}{x}}{1-\frac{1}{x}}=2(\frac{1}{x}+\frac{1}{3x^{3}}+\frac{1}{5x^{5}}+\frac{1}{7x^{7}}+\cdots$(24)を使って計算した。
この方法で、 5000 桁の$10_{\epsilon\supset}^{\sigma}\pi$を計算した。

$\log\pi=1.1447298\mathrm{S}584940017414342735135305\mathrm{S}7116472948$に短縮できた。 これまでの従来の方法で 431msec かかったものを、 BSA 法を使うことによって、対数の高速高精度計算法として、 Sasaki and Kanada による計算法が知られている。 この方法で上の計算であった。 ここで使った Sasaki and Kanada の方法とは、 元の Sasaki and Kanada の方を行うと$375\mathrm{m}\sec$$385\mathrm{m}\sec$

法では、 高精度数の乗算には FTT を使用していないがここででは、 FFT を利用した乗算を使った。 さらに引数を 1000 分の 1 以下にする条件を 10000 分の 1 にした。

Sasaki and Kanada の方法では、 高精度で計算された log2 と円周率の定数が必要であるがこれは事前に計算されているとしている。 このため上の計算時問にはこれら定数計算の時問を含まない。
BSA 法による計算でも高精度の log2 の定数が必要であるが、 毎回計算を行っている。 もし log2 の計算が必要ないならば、 上の問題の計算時聞はさらに速くなり$250\mathrm{m}\sec$となる。

178
3まとめ

指数関数や対数関数の Taylor 展開に BSA 法を適用することによって、 BSA を使わない従来の方法に比べ 40 %程度の高速化ができた。対数関数に対しては、 5000 桁程度の精度で最も高速な計算方法として知られた Sasaki and Kanada の方法を超えることを示した。

参考文献[1] Haible B., Papanikolaou T.: Fast multiprecision evaluation of series of rational numbers,
Technical Report No. TI-7/97, Darmstadt University of Technology, http://www.informatik.tudarmst$\mathrm{a}\mathrm{d}\mathrm{t}.\mathrm{d}\mathrm{e}/\mathrm{T}\mathrm{I}/\mathrm{M}\mathrm{i}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{b}\mathrm{e}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}/\mathrm{p}\mathrm{a}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{k}/(1997)$[2] 平山、 連分数の多倍精度高速計算法、情報処理学会論文誌 Vo1.41, No$\mathrm{S}\mathrm{I}\mathrm{G}8$, pp. 85-91(2000)[3] 右田剛史、 天野晃、 浅田尚紀、藤野清次: 級数の集約による多倍長数の計算法と報処理学会研究報告、 vo1.98, No. 74, pp. 31-36(1998)$\pi$の計算への応用、 情[4] 右田岡史、天野晃、浅田尚紀、 藤野清次: 級数の再帰的集約による多倍長数の計算法と用、 情報処理学会論文誌、 vo1.40, No. 12, pp. $4193-4200(1999)$[5] Sasaki T., Kanada Y.: Practically Fast Multiple-Precision Evaluation of. $247-250(19\mathrm{S}2)$
mation Processing,$\mathrm{V}\mathrm{o}\mathrm{l}.4,$$\mathrm{L}\mathrm{O}\mathrm{G}(\mathrm{X}),$$\pi$の計算への応

Journal of Infor-$[perp]\backslash ^{\gamma}0.4,$$\mathrm{p}\mathrm{p}$[6] 後保範、金田康正、高橋大介、 無限級数に基づく多数桁計算の演算量劇減を実現する分割有理数化法 ‘. $60-71(1999)$
京都大学数理解析研究所講究録$\mathrm{V}\mathrm{o}\mathrm{l}.1084,$$\mathrm{p}\mathrm{p}$[7] 後保範、 金田康正、 高橋大介、 級数に基づく多数桁計算の演算量削減を実現する分割有理数化法、 情報処理学会論文誌 Vo1.41, No. 6, PP 1811-1819(2000)