多倍長計算手法
平成２５年度第2四半期
目次
1. はじめに2. グリーン関数計算(続き）
3. DE(二重指数関数型積分法)の誤差評価

１．はじめに
平成２５年第1四半期では量子モンテカルロ法による物性スペクトル計算
で正しく計算できる範囲はDD形式では表現できる数値範囲の制限により,β=27,U=10までという事が明白になった。
このため今四半期はieee754-2008にデータ形式を使用して正しく計算出来るパラメータ範囲を更に広げた。結果の検証は静的かつ規則的なパスに対してG(β)=-1でチェックしている。
またループ積分においてDE（二重指数関数型積分法)の誤差の評価法が固まったので今回記載しました。

2。 グリーン関数計算(続き）

前四半期はdt  u  1fac  dt  u

条件数  efac

l2のβ  30, u  10までのケースを扱っていましたが今回はdt  u  1, dt  1(dtu)dt 

条件数  el2のケースやdt  u  1, dt  1 での場合を記述しました。

dt β, dt  u  1, dt  1でl  448における必要有効ビット数
ldd:106 td:159 qd:212 fd:265
q:113 dq:226 tq:339 qq:452 fq:565
必要有効ビット数
350
300
250
u=5u=6200

u=7u=8150

u=9u=10100
50
0102030
β4050

ieee754-2008 DQ形式結果一覧
βｌ448
448
448
448
448
448
448
448
448
448
448
448
448

u2030
4050505050505050505050

精度
9 DQ
10 TQ
10 TQ
5 TQ
6 TQ
7 TQ
8 TQ
8 TQ
9 TQ
10 16bai
9 QQ
10 QQ
10 FQ

結果
OKOK
OK
OK
OK
OK
OK
OK
4桁一致
OK
OK
OK
OK

３２倍精度限界テスト
ieee形式限界テスト
演算精度 32倍精度
beta

u
150
170
180
190
200
250
300
300
300
300
300
300
400
400
400
500
500
500
500
500
500
500
500

l=448,n=100

結果
10 OK
10 OK
10 OK
10 OK
10 OK
10 OK
5 OK
6 OK
7 OK
8 OK
9 NG
10 NG
3 OK
4 OK
5 NG
1 OK
1.5 OK
1.6 OK
1.7 OK
1.8 OK
1.9 NG(8桁一致)
2 NG
3 NG

β  500の場合はdt  1と今までと少し様子が異なる.

β=500,u=1の場合

dt  1の場合, x  dtとすると, 打ち切り誤算はO(x 2 ) e x  x  1となり, xよりちいさくないと正しい結果は得られない.
e x  x  1  xとなるxはx  1.256431209.
limit  1.256431209とおくと
dt  u  1,1 dt  limit でβ l(dtu) dtu 4

条件数  eとなる.

dt  limitの場合正しい結果はだせない.
l  448 だと
β  limit  l  562.88 ≒ 562程度まで.

betau
520
540
550
556
557
558
559
560
561
562
563
564
565
566
568
570
576
578
580
590
600

value
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.1000000000D+01
1 -0.9999999999D+00
1 -0.9999998502D+00
1 -0.1000000000D+01
1 -0.9999992758D+00
1 -0.1000000846D+01
1 -0.9999835911D+00
1 -0.9998855364D+00
1 -0.1448403101D+01
1 -0.2199689781D+01
1 0.2759895707D+09
1 -0.4289705228D+10
1 -0.1575586965D+09
1 -0.6396738371D+21
1 0.8569813910D+32

β=562.88まで正しく計算出来るとの見積もりとよく一致しています。

U=1 での実行結果一覧表(L=448,32倍精度演算)
value-9.9982000000E-01-9.9984000000E-01

550 557 559 561 563 565

β-9.9986000000E-01-9.9988000000E-01-9.9990000000E-01-9.9992000000E-01-9.9994000000E-01-9.9996000000E-01-9.9998000000E-01-1.0000000000E+00-1.0000200000E+00

value幾つかの特殊なケース
ieee 754  2008 32倍精度演算
n  40, L  112
βu
50
100
2004001000

結果
200 OK
100 OK
50 NG
25 NG
10 NG

lの値を小さくすると良くなるケース
演算精度 beta
u8倍精度
508倍精度508倍精度
50
16倍精度
100
32倍精度
250
32倍精度
300

l
101010101010

結果
448 NG
400 NG
300 OK300 OK
448 OK300 OK

別解析と対称性を仮定した計算
O(x2 )  ex  x  1  xとなるxをx0とする.
x0  1.25643120862619222963381471050114
l  448ならl x0  562.88
x  dt , u  1とすると
e x  e xe ()  x 2 x0
2βdt  とし,l  896とすると
lβ  779の場合正しく計算でき,β  780で正しく計算できない。
x2

実測結果は以下のとうり。
β  700OKβ  720OKβ  740OKβ  760OKβ  780NGβ  780β  800β  800の場合dt  776.5256125まで正しい値となっている.
NGの場合dt  778.5714286まで正しい値となっている.

O(x2 ) ex  x  1  xとなるxをx 0とする.
x 0  1.25643120862619222963381471050114
x  dt とすると  lx 0 .またe l  2 m (m : 有効ビット数 )l

これより,x0 l  m  ln(2)m  1009なら  800, l  694はこれを満足する.l  700はぎりぎり.

n
k番目の固有値をa kとするとa n  k  a k (k  1,2,...., )より
2
Snn
2

ta ke take tak 1  e    ( 1  e   1  e )e akk 1n2

( akk 1ake tak  e (  t ) ak  e tak  e (t   ) ak
(1  ek 1 ak)(1  e

分母の値はによらない。t (   m*dt ) ak
分子  e 2(  m*dt ) ak
e 2ak))l m * dtとする。
(0  m  )22
(  m*dt ) ak
e 2(   m*dt ) ak
e 2となりmのところを  mにしてもSは不変なので,よって, Sはm lで極値を取り, 最小値となる.
2

2に対して対称となる.

O(x2 ) ex  x  1  xとなるxをx 0とする.
x 0  1.25643120862619222963381471050114
x  ll    x 0 lとする.

このときの平方根計算の引数xには以下の制限が加わる。
m : 有効ビット数
1 x[m ( m )] x02ln22

 ex3x
これから, ( m ) x 0 ln 2  x
22 ln 23m ln 2 x 0x 2m  1009だとx  753.5750183
x012理論結果l  694l  694l  700l  700  816  820
  808  812x  752.5317269
x  754.3739126

OKNGOKNGx  752.0638271 OKx  753.923073 0 NGOKNG

3. DE(二重指数関数型積分法)の誤差評価
変数変換区間[ ,1   ]に対して区間幅hと分点数を定めると一次元積分での誤差gはg  exp(cN)となる.ここでgはln( N )絶対誤差なので,ここでは相対誤差と大きなオーダーの差はないと仮定する。
またDE公式を与える変換x   (t )では＝ (t 0 ),1     (t 0 )1 1となるt 0はt 0  sinh 1 ( ln( ))で与えられる。
   3.187  10 16 ,2 52    2 51とするとt 0  3.125となる。
区間幅h  0.125とすると分点数N 2は51, h  0.25とすると分点数
N1  26となる。また区間幅h  0.15625とすると分点数N 2は41,h  0.3125とすると分点数N1  21となる。
倍精度演算での評価ではこれらの区間幅, 分点数で行うのは妥当と言える。
通常解析解は不明なので, 誤差評価式のcを求めるには2つの分点数N1 , N 2での積分値y1 , y 2の差から求める必要がある。
g | y1  y 2 | exp(cN1cN 2)  exp()となるが, N1とN 2のln( N1 )ln( N 2 )大きさからg | y1  y 2 | exp(前後である.これよりc cN1)でcを求めても誤差は1%ln( N1 )ln( N1 ) 1ln( )となる。
N1gn次元の場合, 分点数( N1 , N1 ,...., N1 )のときの積分値をy1
分点数( N 2 , N 2 ,...., N 2 )のときの積分値をy 2とする.
また分点数( N1 , N1 ,...., N1 , N 2 , N1 ,..., N1 )( N 2 : i番目の値)のときの積分値をy1iとして誤差g i | y1  y1i |とすると
gngi 12i

 n g max

式g | y1  y 2 | exp(( g max  max g i )で一次元のときの評価
1i  ncN1)からcを求める事が出来る。
ln( N1 )ここでg iにばらつきがあればcの値は小さくなり, g iにばらつきがなければ, cの値は大きくなるので, cの値からある次元の分点数を増やして精度を良くするという対策が取れる。
次ぎに端点付近の影響を見る場合, DD形式の4倍精度演算ではx   (t )のコーデイングを一部修正するだけで変数変換区間[ ,1   ]で  2 1024とできるため(あくまで原理的にと言う意味で他の制限が加わる事もある),これまでと同じ方法で誤算の評価が出来る.
倍精度での変数変換区間[ 1 ,1   ](0   1   )とすると

11t1  sinh 1 ( ln( ))でt1を求める。ここで 1   (t1 ).

1t t
DEの性質より,  1  でも 1 0 は小さくて済む。またh
コーデイングの手間を考え, hの値は変更せずに  j  iのjの値をt1に対応する値に変更するだけに出来る。
例として, h  0.125の場合, N 2  51からN 2  61にすると

  3.187  10 16 から,  1  6.54  10 55となる.ソースの修正部分は
61 n) ,51n  5で約2.5倍の増加で済む。ここで変数変換区間[ ,1   ] 26  i  36  iのみとなる.演算量はn次元積分で(と[ 1 ,1   ]での積分値により,0付近の影響をみる事ができる。

5次元積分laportad,laportafの誤算評価(ケース1)
lａportaｄは2次元目と4次元目の誤差が大きい
laportafは誤差のばらつきはない。
DE誤差一覧表
N=21 & 41 h=0.3125,h=0.15625
h=0.15625,h=0.3125
係数c
誤差
laportadlaportaf
laportadlaportaf
gosa
8.715298E-06 2.189189E-07
1.69
2.22
gosa1
9.731227E-07 1.544159E-08
gosa2
7.788459E-06 6.837236E-08
gosa3
3.333854E-08 1.132839E-07
gosa4
7.613123E-06 1.506979E-08
gosa5
1.951379E-08 3.846251E-08
sum gosa 1.642756E-05 2.506301E-07

N=25 & N=49 h=0.3125,h=0.15625
誤差
laportadlaportaf
gosa
8.715298E-06 2.189189E-07
gosa1
9.731227E-07 1.544159E-08
gosa2
7.788459E-06 6.837236E-08
gosa3
3.333854E-08 1.132839E-07
gosa4
7.613123E-06 1.506979E-08
gosa5
1.951379E-08 3.846251E-08
sum gosa 1.642756E-05 2.506301E-07

laportad.laportafともに0近傍の積分への影響はない。

5次元積分laportad,laportafの誤算評価(ケース2)
lａportaｄは2次元目と4次元目の誤差が大きい
laportafは誤差のばらつきはない。
DE誤差一覧表
N=26 & 51 h=0.25,h=0.125
誤差gosa
gosa1
gosa2
gosa3
gosa4
gosa5
sum gosa

係数claportadlaportaf
laportadlaportaf
1.179577E-06 1.001327E-08
1.71
2.31
1.778147E-08 1.016066E-09
1.198222E-06 5.588154E-09
8.821643E-11 4.457539E-09
1.199953E-06 7.914422E-10
9.773612E-10 1.945273E-09
2.417022E-06 1.379847E-08

N=31 & N=61 h=0.25,h=0.125
誤差
gosa
gosa1
gosa2
gosa3
gosa4
gosa5
sum gosa

laportadlaportaf
1.179577E-06 1.001327E-08
1.778147E-08 1.016066E-09
1.198222E-06 5.588154E-09
8.821643E-11 4.457539E-09
1.199953E-06 7.914422E-10
9.773612E-10 1.945273E-09
2.417022E-06 1.379847E-08

laportad.laportafともに0近傍の積分への影響はない。