多倍長計算手法
平成２４年度第１四半期

目次
1. DD形式の利点
2 . 精度保証

1.

DD形式の利点
浮動小数点演算では、近接する２つの数a, bに対しa  bの演算で桁落ちが発生します。
例えば, x  1   , (1    0)で1  xを計算する場合があげられます。
倍精度演算と単精度変数を２つつなげた擬倍精度演算を比べた場合、有効ビット数は倍精度変数のほうが多いですが,1  xを計算した場合の桁落ちも、倍精度変数のほうが大きくなります。これは倍精度変数では1  2 53と１の間の実数xを表す事ができませんが、擬倍精度変数ｘはx  a  b(a, b単精度)a  1, b  0で表す事が可能になる事でわかります。
x  1  10  ii= 11i= 12i= 13i= 14i= 15i= 16i= 17の場合の1  xの値は以下の様になっています。

倍精度演算
0.1000000082740371D-10
0.9999778782798785D-12
0.1000310945187266D-12
0.9992007221626409D-14
0.9992007221626409D-15
0.1110223024625157D-15
0.0000000000000000D+00

擬倍精度演算i= 11 0.9999999960041972D-11
i= 12 0.9999999960041972D-12
i= 13 0.9999999824516700D-13
i= 14 0.9999999824516700D-14
i= 15 0.1000000003627494D-14
i= 16 0.1000000016862384D-15
i= 17 0.9999999837751590D-17同様の事は、多倍長演算でのDD形式とieee形式の演算を比較する場合も同じ様になります。
すなわち、計算において、DD形式で表現できる数値範囲で済む場合は、精度上DD形式が有利になります。
例としては,1x 1(1  x)  1 dx  B ( ,  ) [ ,   0]

ieee形式が有利

1DD形式が有利

011  x  x   dx ≒ 4 ln(  )22(  10 60 )0があります。
これをより一般的な多次元数値積分に適用した例を以下に示します。
(0,1)の浮動小数点数xには,   x  1  となる1    0が存在します。変数変換区間[0,1]の二重指数関数型積分es
法を使用する場合, x  1  はs  sinh(t ), x  sから
2e  escosh(t )  4s 2   21≒ e 。重み係数はw  s≒
≒ 2 s ≒  ln( )。
s 22s(e  e )e
11tの値はt  sinh 1 ( ln( ))となります。具体的な数値としては以下
12sの様になります。

  10 30 の場合, t  3.784,   10 300 の場合, t  6.086となります。
2tとなり、演算量的にh
実際の計算ではh  10 3より大きな値をとる必要があります。
刻み幅をhとすると分点数NはN 

n次元有界領域においてD  0, またはD  0の場合,
X  ( x1 , x2 ,..., xn )  xi  1   (iは1からnのなかの数)でw1w2 ...wn
D ( X )の絶対値E, 積分値の絶対値Fとすると、hn
E / Fの値より使用する演算精度で得られる結果の誤差を推定する事ができます。

INFRABOX
1 1 x 1 x  yI   0001dzdydx
D2D   sxy  tz (1  x  y  z )  ( x  y )2 (1  x  y )(1  x  y  z ) mes  0, t  0, | s |, | t | me2

 z (1  x  y ) m f2

 22

変数変換をして整理すると
x (1  x )dzdydx
D20 0 01 1 1I  D   sx 2 y (1  y )  x2  (1  x ) 2 [tz 2  ( t  m fP   sx 2 y (1  y )  x2  (1  x ) 2 me
2 me ) z  me ]222Q   sx y (1  y )  x  (1  x ) [ me2222 t( t  mf2

 me2tx (1  x )x (1  x )x (1  x )とすると,。以上から、
P2D2Q22)2 ]x (1  x )x (1  x )0 0 P 2 dydx, 0 0 Q 2 dydxでの分点での値と解析近似解
1 11 1を比較すれば誤差が推定できます。二重指数関数型積分法
ではある分点を境に急速に分点での絶対値が減少しますので誤差推定は少ない点で計算できるという利点もあります。
いま、刻み幅h, x  1   , y h 2 2 ln(値は,1)(2sP2) ln(s2)

2s

として二重指数関数型でのP  22( t  m f )1sln( 2 ) ln()2222(  s )( t  m f )me m f2

これと, I 2との比(相対誤差, Iを２で割っているのは, 二重指数関数型積分の値の対称性より)はs  500 2 , t  150 2 , m f  150, me  0.0005,   10 30h  0.56  596 / 1024

  10 34 (ieee 4倍精度)で6.6207  10 -8  10 32 ( DD 4倍精度)で6.2313  10 - 4A( z )  m f とし,2x (1  x )
dydx
Q20 0

1 1I  Q   sx 2 y (1  y )  x2  m f2で計算しても,   10 32 ( DD 4倍精度)で1.3799  10 - 4となり10進6桁の精度を得るにはDD形式の４倍精度では困難な事がわかります。

実測結果は以下の様になりました。

DD形式４倍精度演算結果
h=0.5q0**6*593.0q0/1024.0q0
do i=1,nx=-595+it=x*h
y=pi*0.5q0*sinh(t)x30(i)=exp(y)/(exp(y)+exp(-y))
gw30(i)=cosh(t)*pi/(exp(y)+exp(-y))**2
end do

s=
t=
me=
mf=
ramda=-250000.00000000000000000000000000000-22500.000000000000000000000000000000
0.49999999999999999999999999999999978E-003
150.00000000000000000000000000000000
0.99999999999999999999999999999999605E-030

result= 0.35615337207937480457621671164174208E-006

これを改善するために、
x30(i )  exp( y ) /(exp( y )  exp( y ))をx30(i )  1.0q0 /(1.0q0  exp(2.0q0 * y ))とプログラムを変更します。

ieee形式４倍精度
s= -250000.000000000000000000000000000
t= -22500.0000000000000000000000000000
fme= 5.000000000000000000000000000000000E-0004
fmf= 150.000000000000000000000000000000
ramda= 1.000000000000000000000000000000000E-0030
size=1024
result= 3.561732905959858107980654331577447E-0007

DD形式４倍精度
size=1024
s= -250000.000000000
t= -22500.0000000000
mf**2= 22500.0000000000
me**2= 2.500000000000000E-007
ramda**2= 1.000000000000000E-060
result= 0.356173662449695233D-06

変数変換区間を[1,1]で行う必要がある場合はx30(i )  (exp( y )  exp( y )) /(exp( y )  exp( y ))を
y  0なら
x30(i )  (1.0q 0  exp(2.0q 0 * y )) /(1.0q 0  exp(2.0q 0 * y ))y  0ならx30(i )  (1.0q 0  exp(2.0q 0 * y )) /(1.0q 0  exp(2.0q 0 * y ))と修正するだけで, 精度改善が実施できます。

N  4096まで使用すると,   10 100 でも
解析近似解

 0.111864900351 10 5に対し変数変換区間[0,1]  0.111864900547  10 5
変数変換区間[1,1]  0.111864900549  10 5と非常に良い結果となり, N  8192まで使用
すると,   10 150 でも(2の計算があるのでDD形式では限界ともいえます。）
解析近似解  0.166327423367  10 5に対し,変数変換区間[0,1]で0.166327137485  10 -5とieee形式では32倍精度でないと良い精度の計算ができないものまで計算できます。

ただし、コーデイングでは,アンダーフロー,オーバーフローを防ぐ必要があります。

gw30(i1)*gw30(i2)*gw30(i3)/d**2を(gw30(i1)/d)*(gw30(i2)/d)*gw30(i3)に修正が必要です。1に近いxを取る様にすると、gw30(i1)*gw30(i2)*gw30(i3）、d**2 でアンダーフローが発生する場合があるためです。

2.

精度保証

誤差を含まない変数配列a, bに対する, 総和、内積nnk 1k 1s  s   a (i ), s  s   a (i ) * b(i )の計算でsの精度をa, bと同じ精度を保証する方式を考えます。一般に浮動小数点演算をマイクロ命令レベルで見ると、FPGAなどでわかるように, 整数演算、シフト演算、論理演算、マスク演算で行われています。これをFORTRANやCで行う(整数演算方式）事を考えます。
通常の演算では、浮動小数点演算方式と整数演算方式を比べると, 乗算処理はあまり差がなく、加減算の処理に幾分差がでます。精度を保証する場合は、乗算処理は、整数演算方式が有利となり, 加減算処理
では, 浮動小数点演算方式で行う事が非常に困難になります。
一般に0以外の浮動小数点数は2 n (1.F) (0  F  1)で表せます
ので,0以外の浮動小数点数の絶対値最小の数2 mに対し2 2 mを単位とし, 絶対値最大の数の二乗分の整数テーブル(内積演算途中でのアンダーフロー, オーバーフローを防ぐため)を作成しそれを使って加減算処理をします。
テーブルは正の数値用IPと負の数値用IMの２つを持ち, a (i ), a (i ) * b(i )が正の場合はIPに, 負の場合は絶対値をIMに加算します。
0の場合は何もする必要がありません。
すべてのデータの処理が終了した後, IP  IMならIP  IM ,IP  IMなら, IM  IPを計算し, 所要の精度の浮動小数点数に変換します。
( IP  IMの場合は先頭ビットを１にします。）
これを, T 2 Kで6倍精度,8倍精度の通常の内積演算と精度保証付きの内積演算を行った結果は以下の様になっています。
ここで, 通常の内積演算は浮動小数点演算方式で, 精度保証付きの内積演算は整数演算方式で行っています。

T2k 内積演算性能比較

内積演算 N=4,000,000*MPI数 演算量

精度

６倍精度

８倍精度

ＭＰＩ数

8*MPI数 MFLOP

DD形式（通常演算） IEEE形式（精度保証）
実行時間 性能
実行時間 性能(秒）(MFLOPs) (秒）(MFLOPs)
640.909 563.26
0.807 634.45
1280.902 1135.25
0.812 1261.08
256
0.904 2265.49
0.848 2415.09
512
0.904 4338.75
0.811 5050.55
64
1.85 276.76
1.312 390.24
128
1.85 553.51
1.344761.9
256
1.862 1099.89
1.331 1538.69
512
1.867 2193.89
1.331 3077.39