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

目次
1. はじめに
2. 数値積分の注意事項
3. 3次元反復法の反復回数
3.1 bcg法
3.2 cgs法
4. Hadamard積分
5.量子モンテカルロ法による物性スペクトル計算
5.1 指数部11ビットの限界について
5.2 188倍精度演算結果

1.はじめに数値積分では,積分計算そのものにデータ形式に依存するものがあり、また数学的には発散する積分値列をある範囲内でしか適用できないものがあります。この資料ではこれらをまとめた結果を記載しています。
また量子モンテカルロ法によるスペクトル計算で物理的に有効とするためにはパラメータの設定である条件がつきます。この条件下ではIeee754-2008データ形式での限界が存在します。前四半期では実行時間の関係で粒子数N=20としていましたが,本来のN=100での計算を実施し,N=20との差を求めました。

2.数値積分の注意事項
物理, 数学と計算機シュミレーションには以下の差が生じます。
端点特異点が有り、積分値が有限の場合
１

1x０dx  2

この場合, 計算機による結果の精度は, 積分区間[,1]のの値に依存します。
１1

 x dx  

この場合, 端点特異点が有り、積分値が無限ですが０1
計算機による結果は有限値で積分区間[,1]の log( )となります.(相対誤差1%以下).
1

1
1M2dx

log()の計算は2 (1  x )2
M 2  2
20 xMI

計算機による結果は, 積分区間[,1  ]のの値に依存しM2が必要となります。
2
1
1
1M2
相対誤差はI 0  dx

log()
2
 (1  x )2M 2  2( M 2  2 )  2 xM
精度の良い結果を得るには とIの値から求まります。
  0では,

11log( )2MM2
|I I|
2  0では 0
|I|
M2M2log( 2 )log( 2 )

 30
 75
M  100,   10 ,   10で相対誤差6.8  10 14log(1   M2)
2≒

  10  64 で相対誤差0.0047,   10  60 で相対誤差0.0625となります。
いくつかの計算例をこの後に記しました。

数値積分計算
１1
０ xdx  
変数変換区間[,1  2 113 ]  10 n1
相対誤差 : 積分値と log( )n値
相対誤差10 0.233407086667097208D+02 0.7842D-0220 0.465749122722682359D+02 0.7790D-02
30 0.699614108793860793D+02 0.7780D-02
40 0.927039011434535996D+02 0.7776D-02
50 0.115383602996870176D+03 0.7775D-02
60 0.139187256757241709D+03 0.7774D-02
70 0.162731950845788508D+03 0.7773D-02
80 0.184403258404195498D+03 0.7773D-02
90 0.208959620523617147D+03 0.7773D-02
100 0.233113889859605513D+03 0.7773D-02

数値積分計算
１11M2
０ xM 2  (1  x ) 2 dx  ( M 2  2 ) log( 2 )
  10  m , 変数変換区間 [  ,1   ]  10  n , h  0 . 5 6 , M  100nm=300.147365445951618924Q-01
10 0.231574255058770234Q-02
20 0.462116366109607907Q-02
30 0.694164160806739548Q-02
40 0.919820579107975093Q-02
50 0.114485346609728248Q-01
60 0.138103737135693498Q-01
70 0.147363376638424290Q-01
80 0.147363377222335175Q-01
90 0.147363377222335175Q-01
100 0.147363377222335175Q-01相対誤差
0.8429Q+00
0.6864Q+00
0.5290Q+00
0.3758Q+00
0.2231Q+00
0.6285Q-01
0.1404Q-04
0.1404Q-04
0.1404Q-04
0.1404Q-04

数値積分計算
１1
1M2dx

log()
０ xM 2  (1  x ) 2( M 2  2 )2
  10  m , 変数変換区間 [  ,1   ]  10  n , h  0 .5 7 , M  100

nm=30
0.147365445951618924Q-01
10 0.232490670862570827Q-02
20 0.463932744416145133Q-02
30 0.696889134800300170Q-02
40 0.923429795271255545Q-0250 0.114934474803299206Q-01
60 0.138645380792536845Q-01
70 0.147365445968630811Q-01
80 0.147365446342540534Q-01
90 0.147365446342540535Q-01
100 0.147365446342540535Q-01

相対誤差0.8422Q+00
0.6852Q+00
0.5271Q+00
0.3734Q+00
0.2201Q+00
0.5917Q-01
0.1154Q-09
0.2653Q-08
0.2653Q-08
0.2653Q-08

数値積分計算
１11M2dx log( 2 )2 (1  x ) 2( M 2  2 )
０ xM
  10  m , 変数変換区間 [  ,1   ]  10  n , h  0 .5 8 , M  100

nm=300.147365445951618924Q-01 相対誤差
10 0.232949769142978581Q-02 0.8419Q+00
20 0.464842710263577845Q-02 0.6846Q+00
30 0.698254290616501867Q-02 0.5262Q+00
40 0.925237939730574551Q-02 0.3721Q+00
50 0.115159479054159220Q-01 0.2185Q+00
60 0.138916715932203850Q-01 0.5733Q-01
70 0.147365445665596625Q-01 0.1941Q-08
80 0.147365445951618909Q-01 0.9727Q-16
90 0.147365445951618910Q-01 0.9649Q-16
100 0.147365445951618910Q-01 0.9649Q-16

数値積分計算
１11M2
０ xM 2  (1  x ) 2 dx  ( M 2  2 ) log( 2 )
  10  m , 変数変換区間 [  ,1   ]  10  n , h  0 .5 9 , M  100nm=300.147365445951618924Q-01 相対誤差
10 0.233179540879295303Q-02 0.8418Q+00
20 0.465298137363800405Q-02 0.6843Q+00
30 0.698937535733912137Q-02 0.5257Q+00
40 0.926142896060555719Q-02 0.3715Q+0050 0.115272091218842446Q-01 0.2178Q+00
60 0.139052511237161374Q-01 0.5641Q-01
70 0.147365445704516375Q-01 0.1677Q-08
80 0.147365445951618924Q-01 0.6627Q-18
90 0.147365445951618924Q-01 0.1458Q-28
100 0.147365445951618924Q-01 0.1470Q-30

数値積分計算
１11M2
０ xM 2  (1  x )2 dx  ( M 2  2 ) log( 2 )
  10  m , 変数変換区間 [  ,1   ]  10  n , h  0 .5 10 , M  100nm=300.147365445951618924Q-01 相対誤差
10 0.233294482396574749Q-02 0.8417Q+00
20 0.465525961958223562Q-02 0.6841Q+00
30 0.699279325095223009Q-02 0.5255Q+00
40 0.926595595250995594Q-02 0.3712Q+00
50 0.115328424811043277Q-01 0.2174Q+00
60 0.139120440781553936Q-01 0.5595Q-01
70 0.147365445722666186Q-01 0.1554Q-08
80 0.147365445951618924Q-01 0.6081Q-18
90 0.147365445951618924Q-01 0.1329Q-28
100 0.147365445951618924Q-01 0.1960Q-31以上の例からh  0.5 6 ,   10 30 の場合  10 100 程度必要になります。これは一次元の例ですが多次元や被積分150の場合,dd形式(指数部11ビット)
では,
10程度
nDしか取れない場合が多 々発生します。

関数が

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  129cgs(conjugategradientsquare)法R  100, nx  ny  nz  113
収束判定値は共役残差10 12
初期値u ( x , y, z)  0

反復回数10000まででは、bcg法は倍精度演算,ｃｇｓ法は4倍精度演算まで収束しません。

3.1 bcg法

bcg 法反復回数一覧
演算精度 SR16000 x5570
E5-2670 Phi5110P
3倍精度
482
4倍精度
316
308
308
314
5倍精度
317
6倍精度
317
7倍精度
317
8倍精度
317

性能面からは4倍精度演算が最も適していると言えます。

3.2 cgs法

cgs 法反復回数一覧表
E5-2670 Phi5110P
演算精度 SR16000 x5570
368368
373
5倍精度
240
240
237
281
6倍精度
217
217
217
7倍精度
217
217
217
217
8倍精度

(注 )
1.赤字はその計算機で最 も性能がよかったもの 。
2.6倍精度の有効ビット数 はSR16000は159ビット, 他は177ビット
3.E5  2670, Phi 5110Pでは内積部分に無限精度演算を適 用しています。

性能面からは6倍精度演算が最も適していると言えます。

4.Hadamard積分
ループ積分では内部特異点がある場合と端点特異点がある場合があるが, 両方含む場合の計算は 以下の資料をもとに問題を作成し精度の検証を行った。
数理解析研究所講究録第791巻
1992年 206  219
Hadamard有限部分積分に対するDE公式
緒方秀教, 杉原正顯, 森正武 : 東京大学工学部物理工学科
P.215
11

1F( x )
f .p. dx , F( z )  (1  z ) 4 (1  z ) 42
1 ( x   )5( 46)3
厳密値   (1   ) 4 (1   ) 42を変形し,( 47)端点( x  0orx  1)と内部特異点を持つ例題を作成してDE  eps算法を検証しました。
厳密値の計算は以下の とうり。
F1 ( z )  z

1n(1  z )1n
11F2 ( z )  z n (1  z ) n

1 z n1 t() ,c z1z
t()n ,c 
1 z
1 
1とすると,
F1 ( x )11dx
dt
22  n(x)
t

c00

1I1  

F2 ( x )11dx
dt
22  n(x)
(1)
t
c00

1I2  1dt  0からI1  I 2  0
c0 tn  2の場合 n111 1n  3の場合  ndt  2  3dt  2ct c
3 3 30 t
c3 0c


1
11
 1n  4の場合  ndt dt  
1  44 34ct c0 t
2 02cc厳密解の計算例
11(1  x ) 4 x 4Ｉ2  dx( x  ) 2011xt44t 34t () ,x  4, dx 1 xt 1( t 4  1) 2(1   ) t 4  x よりt4 1
4t 4
14t 4
I2  dt
dt
422 
42((1)
t)
(1)
(tc)
00)1 (c 

4t 444cdt   4dt   4
dt
42
 c)
 c)
 c) 20 (t0 (t0 (tE   (t41t
31dt 
dtから244

4c t  c c)
4c ( t  c )t
1
1| tt  
dt   4
dt
0 44t ctct
c
00E 

14(1  x ) x( x  ) 20

1I1  14dxでも,1 x 414t 3) ,x  4, dx   4xt 1( t  1) 21t (  ( t 4  1)  1よりt4 1

4t 4
14t 4I1  dt
dt
44 ) 2
2  c) 20 (1   t0 (tx (c 1 )となり,


I1 , I 2 はともにEを計算すれば良い事になります。


1
11dt [ 24c2 c 0 t 0 tE  

1414t  c y, dt  c dyからE  cdt  0 t

1


2 c012

c1
14c (y2dt ]  
2 c t0dy   1)

4c

これから
5I1  

3

3

5



1

 4

 (1   ) 4 , I 2    4 (1   ) 4 となります。
443
4.

12

cdt

より一般のnに対して厳密解を求める場合は以下の計算を使用できます。
1
dt


cot() ( n  2)
0 t n  1nn
1n 



1 
1c1
ndt
dy

ccot() ( n  2)n
0 t n  c
c 0 y 1nny

t, dy 111nndt,tc
c(y 1)
1cncn
テストしたのは, 以下の8ケースで良い結果が得られています。
例題
11

1F( x )f .p. dx , F(z)  (1  z) 4 (1  z) 42
1 ( x   )
5

3


厳密値   (1   ) 4 (1   ) 4
2
ex : 25 - 0.1570797130571282Q  01 0.954Q - 05
ramda  0.99999999999999999999999999999999605E - 030
kai  - 1.5707963267948966192313216916389764

例題2111F( x )dx, F(z)  (1  z) 4 (1  z) 4f .p. 
2
1 ( x   )
5

3

53


厳密値   (1  ) 4 (1  ) 4
2
1変数変換t  ( x  1)で
2
1 1
1F( x )
dx, F(z)  (1  z) 4 z 4
2  f .p.20 (2x  1  )
厳密値   (1  ) 4 (1  ) 4
2

ex : 35 - 0.1570796327089221Q  01 0.603Q - 09
ramda  0.99999999999999999999999999999999605E - 030
kai  - 1.5707963267948966192313216916389764

例題3
1

1

1F( x )
dx, F(z)  (1  z) 4 z 4f .p.20 ( x  )
5

3


厳密値    4 (1  ) 4
4

ex : 35 - 0.1511499660473267Q  02 0.100Q - 04
ramda  0.99999999999999999999999999999999692E - 001
kai  - 15.114994701951815421617318813684932
例題41

11


F( x )
dx, F(z)  (1  z) 4 z 4f .p.20 ( x  )
3

5


厳密値    4 (1  ) 44
ex : 35 - 0.1511499660473267Q  02 0.100Q - 04
ramda  0.89999999999999999999999999999999877
kai  - 15.114994701951815421617318813685031

例題5111F( x )33
f .p.
dx,
F(z)(1z)z
20 ( x  )
厳密値  43 (1  )

2
33 3
ex : 35 - 0.1397353864686661Q  02 0.107Q - 04
ramda  0.99999999999999999999999999999999692E - 001
kai  - 13.973535984464536977830852668014869例題61

1

1F( x )
33
f .p.
dx,F(z)(1z)z
20 ( x  )
厳密値  23

 (1  )

4
33 3
ex : 35 - 0.1397353864686661Q  02 0.107Q - 04
ramda  0.89999999999999999999999999999999877
kai  - 13.973535984464536977830852668014277

例題71

1

1F( x )
f .p.
dx, F(z)  (1  z) 2 z 2
20 ( x  )
厳密値  0ex : 35 0.1404624934040962Q - 08 0.723Q - 08
ramda  0.99999999999999999999999999999999692E - 001
kai  0.00000000000000000000000000000000000E  000
例題8111F( x )
f .p.
dx, F(z)  (1  z) 2 z 2
20 ( x  )

厳密値  0ex : 35 0.1404624933843097Q - 08 0.723Q - 08
ramda  0.89999999999999999999999999999999877
kai  0.00000000000000000000000000000000000E  000

5.量子モンテカルロ法による物性スペクトル計算

量子モンテカルロ法では, 安定した計算を行うにはLはL  20  の関係が満たされている事が必要になります。
U  10ではプログラムを実行 中に変数がとる値の最大値はe 200となり, 指数部が11ビットだと  50.14,15ビットだと  802.978となります。
また粒子数Nは100程度が必要となります。

5.1 指数部が11ビットの限界について
量子モンテカルロ法  可変, L  20  , U  10
倍精度の最大値 1.7976931348623157E  0308
l=980
beta = 49.0000000000000000000000000000000

ispin =1

QDR = 1.785864139172194426883476026500380E + 0303
I + QDR = 1.785864139172194426883476026500380E + 0303

ispin =2
QDR = 1.016926695260681160826110539840909E + 0304
I + QDR = 1.016926695260681160826110539840909E + 0304
l=1000beta = 50.0000000000000000000000000000000
spin =1
QDR = 2.862712242794756891557352959861312E + 0309
I + QDR = 2.862712242794756891557352959861312E + 0309

ispin =2
QDR = 1.651604970572476540194914368437567E + 0310
I + QDR = 1.651604970572476540194914368437567E + 0310

  49.6

OK,   49.7

NGl =992

beta =49.6000000000000000000000000000000

ispin =1

QDR = 9.436165589130513054255944449201070E + 0306
I + QDR = 9.436165589130513054255944449201070E + 0306
ispin =2

QDR = 5.415677874136262454106581048120274E + 0307
I + QDR = 5.415677874136262454106581048120274E + 0307
l =994
beta = 49.7000000000000000000000000000000
ispin =QDR

1

= 3.938127827940264707461183253762617E + 0307

I + QDR = 3.938127827940264707461183253762617E + 0307
ispin =2
QDR = 2.263160356072587757351389081936904E + 0308
I + QDR = 2.263160356072587757351389081936904E + 0308
l =beta =
ispin =

996
49.8000000000000000000000000000000
1QDR = 1.643557573940390154190429173363875E + 0308
I + QDR = 1.643557573940390154190429173363875E + 0308
ispin =2
QDR = 9.457544390246214860747340892604545E + 0308
I + QDR = 9.457544390246214860747340892604545E + 0308
l =beta =

998
49.9000000000000000000000000000000

ispin =1
QDR = 6.859317195871621637905579118618770E + 0308
I + QDR = 6.859317195871621637905579118618770E + 0308
ispin =2
QDR = 3.952228149418581143027425919256247E + 0309
I + QDR = 3.952228149418581143027425919256247E + 0309

5.2 188倍精度演算結果
188倍精度演算で粒子数N=100で実行出来ましたのでこれまでの結果と合わせてまとめて記載しました。
量子モンテカルロ法
精度調査一覧
演算精度 β68倍精度
128倍精度188倍精度L
100
180250

最小値
理論値2000 0.567D-307
3600 0.120D-552
5000 0.130D-767

実測値(N=20) 実測値(N=100)
0.330D-307 0.330D-307
0.514D-553 0.512D-553
0.476D-768 0.468D-768

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