多倍長計算手法

目次
1. ４倍精度データ形式
2 . DD形式に関する注意事項
3 . Rump’s の例題
4 . ファイマンループ積分のmasslessのケース
5. 条件数の大きい線形計算
5．1 非対称Toeplitz行列
5. 2 ヒルベルト行列
6．３次元ポアソン方程式の反復解法
6.1 bcg法
6.2 cgs法

多倍長計算の手法
最近では扱う問題規模が大きくなり,演算量が増えた事により,倍精度演算のみでは,十分な精度の結果が得られなかったりする事もあり,一部の処理部分をより高精度な演算を使用する事が行われています。
また, から を引いて有限値を求める（いわゆるくりこみと言われる）計算は計算機から見れば桁落ちの激しい計算となります。

ここでは拡張倍精度(符号部１ビット,指数部15ビット,仮数部64ビット）より多くの仮数部を持つものを多倍長精度と記述します。多倍長精度演算では計算機のアーキテクチャーに依存し,またコンパイラの最適化手法とも係りがあるため,1つの項目にまとめました。

1 . 4倍精度データ形式
4倍精度には２つのデータ形式があります。
(ア) ４倍精度変数Ｑを２つの倍精度変数A,Bの和で表現する
方式 Q=A+B
更に精度を拡張する場合は倍精度変数の数を増やします。 (３つで６倍精度,4つで８倍精度）(イ) IEEE754-2008に準拠したもの
符号部１ビット,指数部15ビット,仮数部112ビット
更に精度を拡張する場合は仮数部のビット数を増やします。 (仮数部176ビットで６倍精度,240ビット８倍精度）
このあとの記述では,(ア)の形式をＤＤ形式,(イ)の形式をIEEE形式とします。

この事に関連した精度に関する２つの簡単な例題をあとの部分で記載しています。

2. DD形式に関する注意事項
DD形式での多倍長演算の基本は倍精度変数の加減算, 乗除算の結果を２つの倍精度変数の和で表す事です。
a  b, a  b, a  b  c  da, b, c, d倍精度変数;
加算 : c  a  bt cad  (a  (c  t ))  (b  t )減算 : c  a  bt cad  (a  (c  t ))  (b  t )(加算, 減算ともに最適化オプションで括弧をはずした演算順序にならない様に注意が必要です。
)
乗算(乗加算命令あり)
c  a *bd  a *b  c（最適化オプションによりd  0とならないように注意が必要です。
)

乗算(乗加算命令なし)
r  134217729.0d 0c  a *bt1  a * r
a1  t1  (t1  a )( 2 27  1)a 2  a  a1t2  b * r
b1  t 2  (t 2  b)b 2  b  b1d  ((a1* b1  c)  a1* b 2  a 2 * b1)  a 2 * b 2

コンパイラの最適化機能はd=0にならない様にする事と関連します。
アーキテクチャーは乗加算命令がサポート（かつコンパイラのサポート)と関連します。

3. Rump’sの例題
以下の非常に簡単な例題です。

a  77617.0, b  33096.0
f  333.75b 6  a 2 (11a 2b 2  b 6  121b 4  2)  5.5b8 a 2  5.5b 2  1からf  2 a2ba54767
 0.827396059946...
2b66192

計算機で実行しますと,以下の様に全くでたらめな結果がでてきます。

倍精度演算結果 : f  0.118059162071741130 10 224倍精度演算結果 : f  1.172603940053......
これは以下のA+Bの計算時に121ビットの桁落ちが発生することが原因です。

A  333.75b 6  a 2 (11a 2 b 2  b 6  121b 4  2)=-7917111340668961361101134701542942850.0

B= 5.5b 8  7917111340668961361101134701542942848.0
このため,正しい結果を得るには有効桁数が１０進３７桁必要になります。実際,IEEE形式５倍精度(有効桁数１０進４３桁),6倍精度(有効桁数１０進５３桁),DD形式６倍精度（有効桁数１０進４８桁),拡張倍精度変数を２つつなげた（有効桁数１０進３９桁）演算を用いると,f=-0.827396059946…..となります。
(注)拡張倍精度変数をつなげるときは乗算でr  2 27  1(倍精度)をr  233  1に変更します。

4.ファイマンループ積分のmassless のケース
これも非常に簡単な二重積分の計算の場合です。
1 1 x1I  dydx
1
( xy )
0 0
これは, x  0, y  0で特異点(端点)を持つため二重指数関数型積分を使用して計算します。
求める解は以下の様になります。
1  11
1 ( )(  1) 1  2 ( )I   x (1  x)  B( ,  1) 0
 (2  1)
2 (2 )1(注)二重指数関数型積分法

x   (t ), a   ( ), b   ( )baI   f ( x ) dx  g (t )dtg (t )  f ( (t )) ' (t )(b  a )( a  b)r ,c とすると
22

 (t )  r tanh( ' (t ) sinh(t ))  c2
r cosh(t )2 cosh 2 (2sinh(t ))で数値積分を行います。このようなケースは
a=0,b=1で変数値と重みのテーブルを作成します。



1での実測例を示します。
100

解析値=9998.378865794316870325682352395113
DD形式４倍精度 =9328.2209999549855615709671950420315
倍精度 =9328.22099995494682

DD形式の４倍精度は表現できる数値の範囲が倍精度と同じです。
このため,この積分計算は演算精度ではなく,表現できる数値の範囲に依存することがわかります。この事からDD形式はさらに6倍精度,8倍精度演算を使用しても正しい積分値が得られない事になります。
解析値
IEEE４倍精度
拡張 倍精度=9998.378865794316870325682352395113=9998.37886579431687032567867936706=9998.37886579431643419

演算精度の問題か,表現できる数値範囲(アーキテクチャー)の問題かの判断を誤ると,問題解決までに長い時間がかかります。

5. 条件数の大きい線形計算
条件数の大きい線形計算では小さな次元の行列でもその演算精度を挙げる必要が生じます。その典型的な例２つを記載しました。
5.1 非対称Toeplitz行列
Ax=b A:行列,x,bはベクトルＡを定め,解がすべて１となる様にbを設定しています。
その条件数及び演算精度、最大誤差は次ページの結果となっています。

A=20r0120r0
12000
1
2

....
....
....
....

00
000000

0
:
00

0:
00

r
:
000:
00

....
:
....
....

0:
20

0
:
123b=

33 r
3 r
3 r
:
3 r
2r

次元数 N=1000

Toeplitz 行列の条件数と最大誤差一覧表r

条件数
倍精度 ４倍精度
2.5 1.04E+07 5.47E-09 5.55E-28
10 1.21E+257 9.17E+239 9.62E+221
25 2.27E+411 ーーーー 2.67E+376r

８倍精度 １６倍精度 ３２倍精度
2.5 4.87E-66 4.59E-143 4.06E-297
10 1.70E+183 1.09E+106 2.34E-49
25 9.70E+336 4.68E+260 2.28E+106

5.2 ヒルベルト行列
1H (i, j ) (1  i, j  n)i  j 1(n  i  1)!(n  j  1)!
1i jH j ,i  (1)(i  j  1)[(i  1)!( j  1)!]2 (n  i )!(n  j )!
条件数 || H || || H 1 ||n|| H ||  max(1  i  n) | hi ,k |
k 1

N 5

条件数0.943 106

N  10

条件数0.353 1014

行列の次元数N ( 1), 条件数10 mm ≒ 1. 5  N
Hx  b, x(i )  1(i  1,2,...N )となるようにbを設定して連立一次方程式を解きました。
次元数,演算精度、最大誤差は次ページの結果となっています。

ヒルベルト行列最大誤差一覧表
次元数

倍精度 ４倍精度 ８倍精度
2051.616 1.44E-07 2.13E-47
40 182.107
62.034 8.72E-16
60 1971.569
62.804 251.132
80 642.074
300.734 252.653
100 201.355
378.428 242.672

次元数

16倍精度 ３２倍精度
20 1.47E-99 1.44E-276
40 4.09E-86 4.30E-243
60 2.66E-60 1.36E-210
80 3.65E-31 1.08E-178
100 2.58E+00 1.41E-150

6. ３次元ポアソン方程式の反復解法
下の３次元ポアソン方程式を反復解法bcg法とcgs法を使用して解いた結果です。反復回数と演算精度により実行時間がもっとも短くて済む、演算精度がきまります。

6.1 ｂｃｇ 法
演算精度と反復回数は以下の様になっています。
これより,４倍精度演算が良い事がわかります。
またそれぞれの演算精度の誤差状況を次ページ以降に記載しました。
bcg 法

反復回数
3500
3000

2911

2500
2000
反復回数

1500

1000
500

1075

346

173

163

060 68 76 81 86 94 102113177241

ﾋﾞｯﾄ数BCG法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12-log10(残差)倍精度
1917505
513
1001
1009
1017
2005
2013
2790
2798
2806

141210864
20-2-4-6-8

反復回数

BCG法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12-log10(残差)1510556ﾋﾞｯﾄ
80ﾋﾞｯﾄ-5181522293643503203451091
1098

0-10
反復回数

BCG法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12-log10(残差)151051173349658197113
129145161
177

0-5-10
反復回数

4倍精度
5倍精度
6倍精度
7倍精度
8倍精度

6.2 cgs 法
演算精度と反復回数は以下の様になっています。
これより,６倍精度演算が良い事がわかります。
またそれぞれの演算精度の誤差状況を次ページ以降に記載しました。

cgs 法

反復回数
3500
3000

2911

2500
2000

反復回数

1500
1000
587

500

364
151

12809799

102 106 145 177 212

ﾋﾞｯﾄ数CGS法の収束状況
nx=ny=nz=65 収束判定値=1.0q-124倍精度-log10(残差)1510
54倍精度(n=2)
301
320
339
358377396515534553572591

0-5-10
反復回数

CGS法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12
6倍精度､8倍精度-log10(残差)15105-55160
6978
87
96105
114
123
132
141

0-10-15-20-25-30

反復回数6倍精度(n=3)8倍精度(n=4)