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

目次
1. はじめに
2. Hadamard積分
2.1 端点特異点の積分
2.2 一般式計算
3. 3次元反復法の反復回数
3.1 bcg法
3.2 cgs法
4. infra問題
4.1 3次元積分計算
4.2 bsgamma計算
5.量子モンテカルロ法による物性スペクトル計算
5.1 248倍精度演算結果
5.2 超多倍長計算の演算量

1.はじめにファインマンループ積分,量子モンテカルロ法による物性スペクトル計算での精度上の問題が解決を見ましたので,その結果のまとめを中心に記述しました。
また今後,行列計算における反復解法が問題になってくる事が考えられますので,その精度に関して記述しました。

これらの問題はieee754-2008形式の4倍精度及び拡張精度演算が必要になってくる場合が多く,サ-バーx5570,e5430で主に実行しています。

2.Hadamard積分
積分法は二重指数関数型積分法を使用しています。
2.1 端点特異点の積分1 dxの計算
x1111dx   dx  dx  2  dx
x (1  x )
x
1 xx
11dx
 dx  x
1 x
積分変数変換区間 [,1  ]対応する
10

１

1

1

10000

1

100tの値を[a , a ] (a  0)とする。
11dx  lim[ dx ]x (1  x )x (1  x )1 １  00 lim[   cosh( t )dt  (e  e )]aa  aaa110 x (1  x ) dx の積分結果
二重指数関数型積分  10  n , h  0 . 5 16ntmax
30 3.78125 実行結果
解析解
40 4.06250 実行結果
解析解
50 4.28125 実行結果
解析解
60 4.46875 実行結果
解析解
70 4.62500 実行結果
解析解
80 4.75000 実行結果
解析解
90 4.87500 実行結果
解析解
100 5.00000 実行結果
解析解

積分値
0.1377527175764281Q+03
0.1377527175737553Q+03
0.1825335145091959Q+03
0.1825335145056543Q+03
0.2271904668671839Q+03
0.2271904668627759Q+03
0.2740603967483030Q+03
0.2740603967429856Q+03
0.3204203634779741Q+03
0.3204203634717572Q+03
0.3630915590459988Q+030.3630915590389539Q+030.4114434511893274Q+03
0.4114434511813443Q+03
0.4662325224569613Q+030.4662325224479152Q+03

一般式計算2.2

Hadamard積分は
111dx  log( )]x0n 1lim[ n 1lim[   01111dx ( ) n 1 ]n  0n 1 0 xで定義されます .。
今二重指数関数型積分 で分点が[,1  ]にありそれに対応する ｔの値を[ a , a ](a  0)とします。

1 1 1 111I1   dx , I 2   2 dx , I n   n dx ( n  3)とします。
x
 x
 x
1
 1  e  sinh( a )1 n  1では 


1 a

1
1dx  2  dx    cosh( t )dt より
x (1  x )
x
aaa

1
I1 
 cosh( t )dt    cosh( t )dt   sinh(a )2 a0aI 2    cosh( t )e   sinh( t ) dt  e  sinh( a )  e   sinh( a )a1I1  log( )   sinh(a )  log(1  e  sinh( a ) )  e   sinh( a )1I 2   1  e   sinh( a )n 21
 ( i 1) sinh( a )I n  e   ( i 1) sinh( a ) )n  2 C i (ei 0 i  1n 1111( ) n 1 (1   n 1C i (e i sinh( a ) ))n 1 n 1i 1

よって I n n 21111
  ( i 1) sinh( a )( ) n 1  )n  2 C i (en 1 n  1 i 0 i  1  2 53 なら a  3.15
Hadamarde   sinh( a )  10 15.89 でI1  01( n  2)n 1となり, 算法での値とも一致す る。
HadamardIn  3.3次元反復法の反復回数
問題は 3次元ポアソン方程式
uu  R fx
領域[0,1]  [0,1]  [0,1]解析解 u ( x , y, z)  e xyz  sin( x )  sin( y)  sin( z)
非対称問題。
解法
bcg(biconjugategradient )法R  100, nx  ny  nz  129
cgs(conjugategradientsquare)法R  100, nx  ny  nz  65
収束判定値は共役残差 10 12
初期値 u ( x , y, z)  0

3.1 bcg法
bcg法反復回数
精度
56ビット
60ビット
64ビット
68ビット
72ビット
76ビット
3倍精度
84ビット
88ビット
92ビット
96ビット
100ビット
104ビット108ビット
4倍精度

反復回数 正規化
実測/正規化
50001 none
none
4315 node
none
1579
931
1.696
938
826
1.136
784
738
1.062
593
664
0.893
482
600
0.803
439
545
0.806
396
497
0.797
368
455
0.809
336
418
0.804
325
386
0.842
319
357
0.894
313
331
0.946
308
308
1

(注)50001回は50000回反復で収束しなかった事を示しています。
正規化
P有効ビット数=精度のビット数+1

正規化反復回数308  (113)2P

性能的には3倍精度が最も効率が良いと言えます。

3.2 cgs法
e5430 cgs 反復回数
倍精度は10000回で収束せず。

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

反復回数
ieee
ieee smp
dd392
448
226
155 none
130
131
128
128 none
128
128

670
144
128

実行時間(秒)
精度
4倍精度
5倍精度
6倍精度
7倍精度
8倍精度

ieee
159.2898
221.2924
144.4976
171.5329
213.6745

ieee smp
105.5012
75.494268.869469.804461.7321dd
95.3938none
39.1221none
70.0613

q4sum,q5sum 使用
精度
4倍精度
5倍精度

反復回数
実行時間
378
160.5706
155
157.2171(1)ieee形式,dd形式ともに6倍精度演算が最も効率が良い事を示しています。
(2)q4sum,q5sumの様に総和演算を無限精度演算を使用すると,5倍精度の収束が良くなり,ieeeでのsmp実行と同じ反復回数となっています。

dd形式cgs反復回数
ビット数

反復回数
954931962155973506982811
99125010011031011080
10210271038151047021056301066704 infra 問題
4.1 3次元積分計算
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4I  







000001d
D2d  dx 5 dx 4 dx 3 dx 2 dx 1
D  C( x 3 M2 x 6 2 )

 ( x 1  x 2  x 3 )( x 4  x 5 ) 2
 ( x 1  x 2 ) 2 (1  x 1  x 2 ) 2 x 3 ( x 1  x 2 )( x 4  x 5 )C  x 3 (1  x 1  x 2  x 3 )  ( x 1  x 2 )(1  x 1  x 2 )x 6  1  x1  x 2  x 3  x 4  x 5
で  ｎ  10  n での積分値 I n の一次差分 ,二次差分が10進10桁一致する事を検証。

3次元化
1 1 1

Jdxdydz
2D0 0 0I  J  (1  x ) 2 yz
D  [(1  x  xy(1  y)] [ x (1  y) M 2  (1  x )(1  z )2 ] [ xy 2 (1  xy)  (1  x ) 2 z 2
 2 x (1  x ) y(1  y) z ]M=0の場合をcase1,M=91.19/0.1057の場合をcase2,M=125.5/0.1057の場合をcase3としています。

ncase1 2次差分一覧表
N=728
N=1456
N=2912
16 5.3018776181 5.3018981105 5.3018981105
17 5.3020175196 5.3018981105 5.3018981105
18 5.3020700159 5.3018981101 5.3018981105
19 5.3016335118 5.3018981100 5.3018981105
20 5.3023935051 5.3018981077 5.3018981105
21 5.3029475759 5.3018981216 5.3018981105
22 5.2971188923 5.3018980753 5.3018981105
23 5.3041142583 5.3018982336 5.3018981105
24 5.3093273173 5.3018981250 5.3018981105
25 5.2943816325 5.3018976075 5.3018981105
26 5.2794768170 5.3018979849 5.3018981105
27 5.3558839326 5.3019004772 5.3018981105
28 5.3091134086 5.3018977130 5.3018981105

Case1 の場合は1
積分値は log( 2 )の 2次式になるので
二次差分を求めていま す。

 n  10  nnn

case2 一次差分一覧表
N=728
N=1456
N=2912
16 0.000038733241268 0.000038733248507 0.000038733248507
17 0.000038733234064 0.000038733248507 0.000038733248507
18 0.000038733266331 0.000038733248507 0.000038733248507
19 0.000038733326864 0.000038733248507 0.000038733248507
20 0.000038733174086 0.000038733248506 0.000038733248507
21 0.000038733353236 0.000038733248506 0.000038733248507
22 0.000038733544343 0.000038733248510 0.000038733248507
23 0.000038731850666 0.000038733248499 0.000038733248507
24 0.000038732625311 0.000038733248540 0.000038733248507
25 0.000038734471892 0.000038733248546 0.000038733248507
26 0.000038731896322 0.000038733248424 0.000038733248507
27 0.000038724804488 0.000038733248425 0.000038733248507
28 0.000038738845469 0.000038733249112 0.000038733248507
29 0.000038740941884 0.000038733249063 0.000038733248507

case3 一次差分一覧表
N=728
N=1456
N=2912
16 0.000021493099355 0.000021493102821 0.000021493102821
17 0.000021493093286 0.000021493102821 0.000021493102821
18 0.000021493115638 0.000021493102821 0.000021493102821
19 0.000021493142483 0.000021493102821 0.000021493102821
20 0.000021493082925 0.000021493102821 0.000021493102821
21 0.000021493117428 0.000021493102821 0.000021493102821
22 0.000021493302447 0.000021493102823 0.000021493102821
23 0.000021492401183 0.000021493102817 0.000021493102821
24 0.000021492629664 0.000021493102839 0.000021493102821
25 0.000021493704654 0.000021493102844 0.000021493102821
26 0.000021492781172 0.000021493102769 0.000021493102821
27 0.000021488800060 0.000021493102792 0.000021493102821
28 0.000021495565878 0.000021493103126 0.000021493102821
29 0.000021496356431 0.000021493103198 0.000021493102821

1)の一次式2となるので一次差分を 求め
積分値は log(ています。  n  10  n

4.2 bsgamma計算
bsgammaplanarD  c( x 1 ( p 1  m 1 )  x 2 ( p 2  m 2 )  x 3 m 3  x 6 m 6 )222222

 c1 (x 5 p 2  x 4 p1  x 4 x 5 (p1  p 2 ))222222

 c 2 (x 2 p 2  x 1 p1  x 1 x 2 (p1  p 2 ))222222 2x 2 x 3 x 5 p 2  2x 1 x 3 x 4 p1  x 3 (x 2 x 4  x 1 x 5 )(p1  p 2 )2222x 6  1  x1  x 2  x 3  x 4  x 5c1  x 1  x 2  x 3c 2  1  x1  x 2c  x 3 (1  x 1  x 2  x 3 )  (x 1  x 2 )(1  x 1  x 2 )  x 3 (c 2  x 3 )  (1  c 2 )c 2
p1  m b , p 2  m s , m1  m 2  1.5, m 3  80.3477, m b  4.7, m s  0.0942222m6  
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4I 0 00001dx dx dx dx dx
2 5 4 3 2 1D

これを4次元化して計算しています。

nbsgamma 一次差分結果
N=728
N=1024
N=1456
16 0.000881676036303 0.000881676024318 0.000881676024308
17 0.000881675957819 0.000881676024256 0.000881676024308
18 0.000881676127056 0.000881676024379 0.00088167602430819 0.000881676192593 0.000881676024096 0.000881676024308
20 0.000881675929067 0.000881676025850 0.000881676024307
21 0.000881675933258 0.000881676026869 0.000881676024308
22 0.000881675979610 0.000881676024874 0.000881676024305

1
積分値は log( 2 )の一次式となりますの で
一次差分を求めていま す。
 n  10 n5.量子モンテカルロ法による物性スペクトル計算

5.1 248倍精度演算結果
x5570 L=20β ,U=10精度
68倍精度
128倍精度188倍精度248倍精度

実行時間(秒)

有効ビット数 β
2161
4081
6001
7921

n=20100180
250300n=1002790
25460916956
1528008
50123
4375879
1002229036094

絶対値最小値一覧精度
68倍精度
128倍精度188倍精度248倍精度n=20
0.330D-307
0.514D-553
0.476D-768
0.126D-921

n=100
0.330D-307
0.512D-553
0.468D-768
0.121D-921

N=20とN=100での実測値では大きな差は見られません。

5.2 超多倍長計算の演算量
演算量計算
ここで演算量は実行時間に比例すると仮定しています。
演算量I(n)I(100)  I( 20) m
行列積m=3精度
68倍精度
128倍精度188倍精度248倍精度

m
2.805
2.797
2.777
2.797

68倍精度演算のLで正規化した演算量
I(P) P有効ビット数

I( P )  (乗算P) m I ( 2161 )
2161

m=2精度
128倍精度188倍精度248倍精度

m
1.894
1.888
1.902