相加相乗平均による初等超越関数の計算

(1)

>研究論文@

相加相乗平均による初等超越関数の計算

平山弘・加藤俊二

自動車システム開発工学科

Calculation of Elementary Transcendental Function

by Arithmetic-Geometric Mean

Hiroshi HIRAYAMA and Shunji KATOH

Abstract

In 1976, Brent showed that the elementary transcendental functions (exp, log, atan, sin, cosh, etc) can be calculated in O(𝑛𝑛(log 𝑛𝑛)2_{log log 𝑛𝑛) operations with relative error O(2}−𝑛𝑛_{) as 𝑛𝑛 → ∞. This algorithms} depends on the theory of elliptic integrals, using the arithmetic-geometric mean method.

It's learned that this method is effective when calculating by the high precision, but It isn't known in what kind of area this way is effective. In this paper, Elementary transcendental function was calculated and its validity was checked using algorithm of this Brent.

Keywords:multi-precision, elementary transcendental functions, the arithmetic-geometric mean はじめにこれまで多くの高精度計算プログラムが開発され来たが、相加相乗平均法である高速な %UHQW のアルゴリズム _{を使ったプログラムは開発されて来なか} った。このアルゴリズムは極限において非常に効率的なものであるが、桁とか桁程度の計算では、あまり有効ではなかったためではないかと思われる。最近計算機が速くなり記憶装置も十分な大きさを持つようになったことから、%UHQW のアルゴリズムの有効性を調べた。 %UHQW の文献にもあるが、6DODPLQ が提案した計算法に相加相乗平均を使う円周率の高速な計算法がある。この計算法は一時期円周率の計算によく使われた計算法である。その当時、それと類似した方法で対数関数の計算 _{の研究などもある。最近は、} 分割法 _{によって円周率の計算が行われるように} なって来ている。関数の計算も将来的には、その方法で計算が行われると思われる。本論文では、高精度の計算でしか有効でないと言われた相加相乗平均法を使って、初等超越関数を計算し、どの程度有効かを調べ、その有効範囲を決定する。以下の計算には、計算機として ,QWHO L. *+] を使用し、コンパイラーとしては、 0LFURVRIW9LVXDO6WXGLR&を使用した。 ))7 による高精度数の乗算高精度数の演算で問題になるのは計算時間である。このため高速アルゴリズムは不可欠である。その中で、高精度数の乗算は、高速フーリエ変換（))7）を使うと高速に行うことができることが知られている。

n

桁の数値の乗算には通常

n

2のオーダーの計算時間が必要であるが、))7 を使うと

n

log

n

のオーダーで計算できる。相加相乗平均法で高速に計算するためにも、))7 を利用した乗算が不可欠である。

(2)

.DUDWVXED による高速乗算法 _{もある。低精度の} 乗算では通常の乗算法が有効で、高い精度計算では ))7 を利用した乗算方法が有効であることが知られているので、その中間の精度で .DUDWVXED の乗算方法が有効ではないかと推定される。この方法で高精度計算プログラムを作成してみたが、自作のプログラムでは有効な範囲が見つからなかった。このため、自作プログラムには、.DUDWVXED の乗算方法は、複素数の乗算には使っているが、多倍長精度の数値の乗算には使っていない。

n

が十分大きな数の場合すなわち高い精度の場合、 ))7 を使った計算法は非常に効率的である。多くの計算機で、進数で約桁程度を越えれば、))7 アルゴリズムを使った計算法が、通常の計算法より高速に計算できる。 ))7 を使って、高精度数の積を計算するには、通常の二つの計算方法が考えられる。有限体を使って計算する方法と浮動小数点を使って計算する方法である。有限体を使った計算方法の場合、整数演算だけを使うので、誤差の入らない厳密な計算結果が得られる。一方、浮動小数点演算を使う方法は、三角関数などの近似値を扱うために厳密な計算ができない。しかし、最終結果は整数であることが分かっているので、誤差が十分に小さいならば、丸め処理によって厳密な計算が可能である。これまでの多くの計算機では、最も精度の高い整数は、ビットであり、浮動小数点数は、ビットである。大型計算機では、ビットの浮動小数点数もある。このため、計算精度の高い浮動小数点数を使った計算方法が多く使われてきた。最近のマイクロプロセッサでは、ビット精度の整数を扱うことが出来るので有限体を使った計算法も行われ始めている。

2.1. 実数用 FFT による高精度数の乗算法

高精度整数

x

を

b

進数

m

桁で表現する。すなわち

x

x b

_k k k m



 



0 1 とする。このとき高精度整数

y

との積

z

は、



 



j k k j k j

x

y

z

0 となる。この計算は畳み込み演算と呼ばれ高速フーリエ変換によって効率よく計算できることが知られている_{。この場合、倍精度実数を利用した ))7 を} 利用する計算方法がよく使われる。))7 にも多くの計算法があり、今回利用したプログラムには、実数用で基数がおよびの ))7 プログラムである。 ))7 には、複素数用と実数用があるがここで扱うデータが実数であるので、実数用を利用した。実数用は、同じデータ数であるとき複素数用の約倍高速である。基数２の ))7 は、プログラムは短く、よく使われる計算法である。%HUJODQG_{って発表された} 基数のプログラムを使えば、さらに高速に計算できるので、))7 には基数のものを使用した。実数用 ))7 の計算では、計算の途中で打切り誤差が入る。この誤差は、最後の丸め処理によって厳密な値になるためには、次のような関係式を満たさなければならない。この式は +HQULFL１）_{によって導か} れたものである。

b

進数

m

桁の数値を厳密に乗算できるには、計算精度の相対誤差を



マシン・イプシロンとすると 2 2 2

₍

₂

_log

₇

₎

192

1 b

m







を満たさなければならない。この式から基数

b

が大きいほど要求精度が高くなることがわかる。桁数

m

も大ききなると要求精度が高くなる。この式は、相対誤差の２乗以上の高次の項を省略する方法で、誤差を評価しているので、ほぼ十分条件になる。現実にはもっと低い精度でも計算可能である。たとえば、上の式で

b



10000

、





2 .

22 

10

16（,((( 方式の倍精度浮動小数点）の場合、この式では計算可能桁数は桁となる。実際には万桁の数も計算可能である。式より精密な評価を得るために、区間演算などを試みたが、若干適用範囲が増加したが、実用的な範囲にはならなかった。式を誤差の評価式と見たとき、誤差は、基数

b

の２乗、対数部分を省略すると桁数

m

の約２乗に比例することがわかる。この場合の基数

b

とは、１語の中に入る最大の整数という意味になるので、

(3)

誤差が最大になるのは、各語に

b 1

の数値を入れたとき最大の誤差になることがわかる。

b  10000

とすれば、各語にを入れたとき最大の誤差が生じることになる。

b



2

nのときは、各語に

2

n



1

を入れたとき最大の誤差が生じる。計算結果は、整数であることが知られているので、誤差がより小さいならば、四捨五入の計算によって、厳密な計算が可能である。誤差は、最大の数値を入れて計算することによって計算が可能なので、計算できる限界も容易にわかる。

2.2. 最大誤差の計算

前節で述べたことを確かめるために、 ,((( 方式の倍精度浮動小数点をもつ計算機を使って誤差を計算した。

b  100

と固定して、桁数

m

を増やす。このときの誤差は、表１のようになる。表１桁数 P を増やしたときの誤差 P HUURU ( ( ( ( ( ( このときの計算方法は、基数２の ))7 である。桁数

m

が２倍になると、誤差は約４倍となり、上の公式が成り立つことがわかる。次に、桁数P と固定して、基数

b

を変化させる。このとき、誤差は表２のようになる。このときの計算法は、基数２の ))7 である。表２基数 E のビット幅を増やした時の誤差（P ） E (UURU ( ( ( ( ( ( この計算結果から、誤差は、基数

b

のかなり正確に２乗に比例することがわかる。この誤差が、より小さければ、実数を使った ))7 による高精度の数値の乗算が厳密に行うことができる。このようにして、限界を求めると次のようになる。 ))7 として基数８の %HUJODQG のプログラム_を利用した。基数

b

として、２のべき乗とのべき乗を使った。その結果は表３のようである。E 場合の計算は、使用したプログラムの限界で、求められなかった。ここで示した E の結果は、 2RXUD によって作成したプログラム _{による結果で} ある。表．基数の ))7 の場合基数

b

誤差計算可能桁数

m

計算可能桁数進換算 ! !

(4)

2.3. FFT を使った計算時間

高速フーリエ変換))7を使うと高精度の数値を高速に乗算を行うことができる。同じ桁の数値の掛け算を行うとき、進数で約桁（進数で約桁）以上のとき、))7 を使った計算法を使い、それ以下では通常の乗算法を使っている。いろいろな桁数の計算時間を表に示す。計算精度を上げていくと、途中再び通常の計算法が一旦速くなることもあるが、計算時間にそれほど大きな違いにならないので、その精度の時でも、))7 を使った計算法を使用している。通常の計算法と ))7 を利用した計算法の境界の桁数は、高速計算機ほど、大きくなる傾向がある。最初にこのプログラム作成した計算機 ,QWHO 社 Lでは、乗算の限界が進数で約桁であった。相加相乗平均法を計算相加相乗の平均の計算とは、二つの数列

a ,

_n

b

_nに対して、次のような計算を次々と行うことである。 n n n n n n

a

b

a

b

a

_₁





,

_₁



2

この計算は、２乗収束することが知られている。この計算法を利用した関数計算が %UHQW によって提案されている。%UHQW によると、関数

U

(

m

),

T

(

m

)

を、表のように定義する。このとき

))

(

log(

)

(

m

T

m

U



e

U(m)



T

(

m

)

この式を使って指数関数

e

xを計算するには、まず次の方程式を解く。

U

(

m

)



x

方程式を解き解

m

₀を求める。上の関係式から次の式が成り立つ。

e

U(m0)



T

(

m

₀

)

0

m

が求まれば、

T

(

m

₀

)

を計算すれば、

e

xの値が計算できる。の方程式を解くために、次のように 1HZWRQ 法を使う。

)

('

)

(

1 n n n n

m

U

_U

m

_m

x

m

_





（）得られた

m

を使って

T

(m

)

を計算することによって

e

xを求めることができる。逆に

T

(

m

)



x

）を解き、求められた

m

を使って

U

(m

)

を計算することによって

log

x

を求めることができる。自動微分法による高速化 %UHQW>@は 1HZWRQ 法で使う微分係数を計算するために数値微分を使うことを提案している。このような数値微分法では、（）の右辺を計算するた

表４乗算の実行時間（単位 PVHF）

進の桁数 ))7 通常表３ %UHQW の

U

(m

)

と

T

(m

)

の関数 $OJRULWKPIRU

U

(m

)

A



1 B



1 

m

ZKLOH

A



B



2

n/2GR EHJLQ

C



(

A



B

)

/

2 B 

AB

A 

C

HQG

A





/(

A



B

)

S 

m

ZKLOH

1 

S



2

n/2GR EHJLQ

A



A

( S

1 

)

/

2 S



2 S

/(

1 

S

)

HQG UHWXUQ

A 

( S

1 )

/

2

$OJRULWKPIRU

T

(m

)

1 

V

S 

m

ZKLOH

1 

S



2

nGR EHJLQ

W



2 SV

/(

1 

V

2

)

W



W

/(

1 

1 

W

2

)

1 /(

)

(

V

W

VW

W







)

1

1 /(

_W

2

W

V





S



2 S

/(

1 

S

)

HQG UHWXUQ

(

1 

V

)

/(

1 

V

)

(5)

めに最低でも関数

U

(m

)

を２回計算しなければならない。ところが上の式を見ればわかるように、除算や平方根の計算が含まれている。この場合、除算や平方根を含む式を計算する時間より少ない時間で微分係数を計算することができる。次のような行の単純なプログラムで [ について微分することを考える。 S [AT VTUWS S、ｑの [ についての微分を GS、GT とすると GS [ G[GT GS VTUWS として計算される。関数 VTUWSは、元々の関数ですでに計算されているものであるから、微分係数を計算するには平方根の計算は不要になり、高速に計算_{することができる。実際の計算では、平方根を} 計算するには、平方根の逆数を計算し、その値から平方根を計算されているため、上の計算では平方根による除算も不要となり次のように高速化ができる。 S [AW UHFLSBVTUWST W S UHFLSBVTUW[は平方根の逆数関数 GS [ G[GT GS W と計算できる。高精度計算では、短い桁数の整数によるわり算は高速に行えるので、２によるわり算は問題とならず、全体としてかなり高速化できる。わり算を含む式の微分係数の計算は、乗算の回数は増えるものの除算が不要になるのでかなり高速化できる。指数関数

e

xの微分係数は、この関数が微分しても同じ値となるため、この性質を使うことで高速化できる。対数関数

log

x

の微分係数は、この関数の微分が

x

1

_{と単純な関数になることを利用し} て高速化できる。三角関数（

sin

x

、

cos

x

）の計算は、これらの関数が同時に計算すると高速計算することができることから、高速化できる。この方法でプログラムを微分すると

U

(m

)

と

U

(m

)

が同時に計算できるプログラムを作ることができる。

4.1 微分係数を使った効果

指数関数の計算例として、

e

2



4.1132503



を桁と桁の精度で計算した。計算は（）の方程式を倍精度数を使って 1HZWRQ 法で計算し、高精度計算を始めるための初期値を計算した。その値を使い高精度で 1HZWRQ 法を使ってその計算に必要な精度を確保しながら計算した。1HZWRQ 法では、微分係数は計算精度の半分で十分なので、微分係数の計算式を計算精度を半分にして評価した。桁桁差分を利用した方法秒秒微分法を利用した方法秒秒微分係数を差分近似で計算する方法に比べ、約パーセント高速になった。微分係数の計算が差分法に比べ、４倍程度高速になったことになる。相加相乗平均法が有効な範囲相加相乗平均法（$*0 法）の有効範囲を求めて計算精度を上げて調べた。$*0 法として微分係数を利用した高速化した $*0 法を使った。比較の対象は自作の多倍長プログラム033DFNと比較した。ここでは、の問題を使って計算した。 033DFN では、細かい点を除けば、

e

xの

x

を

x

_n

2

として、7D\ORU 展開式から計算し、その値を

n

回乗計算する方法で計算する。

n

は計算する 7D\ORU 展開式の項数と２乗計算の回数がになりように選んでいる。その計算結果を表に示す。$*0 法は万桁以下では通常の多倍長ルーチン033DFNより遅いことがわかった。万桁の計算でようやく速くなった。この多倍長ルーチンは、その当時の計算機が遅かったこともあり、万桁程度の数値想定して作成したものである。このため、さらに改良が行われると思われるが、現段階では $*0 法は万桁を超えたところでは、効率的な計算法と言える。表５ $*0 法と多倍長ルーチン（033DFN）の比較計算精度高速化 $*0 法 033DFN

(6)

まとめ相加相乗平均法（$*0 法）による初等超越関数の計算法を解析的な微分係数の計算（自動微分法）を導入することによって、高速化を行った。高速化された $*0 法で初等超越関数の計算を行い、その計算の有効範囲を調べた。その有効範囲は万桁以上の範囲とわかった。この範囲では個の関数値を計算するのに１～２分以上かかり、あまり実用的とは思えない範囲であった。円周率の計算では、最近 $*0 法を使わない傾向があり、その方法を使えば、通常の計算法がさらに高速化される可能性があり、さらに $*0 法による計算の有効範囲がさらに高精度に追いやられる可能性はある。参考文献 >@ %UHQW 5 3 )DVW PXOWLSOHSUHFLVLRQ HYDOXDWLRQRIHOHPHQWDU\IXQFWLRQV-$VVRF &RPSXW0DFK

>@ 6DODPLQ ( &RPSXWDWLRQ RI π XVLQJ DULWKPHWLFJHRPHWULF PHDQ 0DWK &RPSXW

>@6DVDNL7DQG.DQDGD<3UDFWLFDOO\)DVW

0XOWLSOH3UHFLVLRQ (YDOXDWLRQ RI /RJ[ - ,QIR3URF

>@ +HQULFL 3 $SSOLHG DQG &RPSXWDWLRQDO &RPSOH[ $QDO\VLV 9RO &KDS -RKQ :LOH\ 6RQV1HZ<RUN

>@ %HUJODQG * ' $ 5DGL[(LJKW )DVW )RXULHU7UDQVIRUP6XEURXWLQHIRU5HDO9DOXHG 6HULHV,(((7UDQV$($8SS >@.DUDWVXED$DQG2IPDQ<0XOWLSOLFDWLRQ RI PXWLGLJLW QXPEHUV RQ DXWRPDWD 'RNODG\ $NDG1DXN66659ROSS >@ 後保範金田康正高橋大介：級数に基づく多数桁計算の演算量削減を実現する分割有理数化法，情報処理学会論文誌， >@平山弘：連分数の多倍長精度高速計算法情報処理学会論文誌、

>@ 大浦拓哉 2RXUD’V 0DWKHPDWLFDO 6RIWZDUH 3DFNDJHV WWSZZZNXULPVN\RWRXDFMSRRXUDLQGH[ MKWPO >@ 5DOO/ % $XWRPDWLF 'LIIHUHQWLDWLRQ 7HFKQLTXHDQG$SSOLFDWLRQV/HFWXUH1RWHVLQ &RPSXWHU6FLHQFH9RO6SULQJHU9HUODJ %HUOLQ+HLGHOEHUJ1HZ<RUN