多倍長計算手法
平成２５年度第3四半期
目次
1. はじめに
2. 32倍精度演算を基本とする多倍長演算
2.1 除算2.2 平方根2.3 演算量
2.4 ヒルベルト行列
3.量子モンテカルロ法による物性スペクトル計算
3.1 32,64,96,128倍精度演算の計算結果
3.2 32,64,96,128倍精度演算性能測定結果

１．はじめに
平成２５年第2四半期では量子モンテカルロ法による物性スペクトル計算
で正しく計算できる範囲をieee754-2008データン形式の32倍精度演算まで実施しました。
このデータ形式では表現可能な数値範囲の制限から128倍精度演算まで可能であるので,これを作成して正しく計算出来るパラメータ範囲の限界を求めた。結果の検証は静的かつ規則的なパスに対してG(β)=-1でチェックしている。
128倍精度演算の作成に当たっては,32倍精度変数をつなげる方式を採用しました。
倍精度変数を4つつなげて8倍精度変数を作成した場合,消失する有効ビット数は(241-212)/241=0.12(12%)と大きいが32倍精度変数を4つつなげた場合には(4081-4036)/4081=0.011(1.1%)と消失する有効ビット数の比率が少ない事と4倍精度変数を複数個つなげた8倍精度演算
12倍精度演算,16倍精度演算がすでにあるので作業効率の良さを考慮した事によります。

2.32倍精度演算を基本とする多倍長精度演算
一般に, ある精度を持つ変数をつなげた形式で多倍長 演算を作成した場合に以下の問題が発生します。
演算精度の増加分だけ, 変数の値の範囲が制限されます。
例えば, 倍精度変数を4つつなげて,8倍精度変数と演算を定めると, 倍精度変数では10-308 ~ 10 308 の範囲で精度が保証されていますが ,8倍精度変数では,10  260 ~ 10 260 の範囲でしか精度が保証されません。倍精度変数を5つつなげて,10倍精度変数と演算を定めると, 精度が保証される範囲は10  244 ~ 10 244 となります。
これがieee754形式で多倍長演算を定めたときに量子モンテカルロ法によるスペクトル計算で発生した問題です。
32倍精度変数を複数つなげて多倍長演算を定めるとこれ以外にも問題がでてきます。
ieee754の倍精度変数は指数部のビット数と仮数部のビット数が非常に良いバランスをとっています。ieee754  2008形式ではバランスを考えれば仮数部のビット数は (52  1) * 16ビット  848ビットで 27倍精度変数となり,32倍精度変数は仮数部が160ビット多いという事になります。
また,32倍精度変数を複数つなげると, 精度が保証される範囲が１つ毎に10 303 と非常に大きな制限が加わります。
ただし,32倍精度変数を 4つつなげて128倍精度演算を定める場合と整数演算方式で128倍精度演算を作成する作業量は 2 ~ 3桁の差があります。
このため,32倍精度変数を複数つなげる影響をみるために, 除算, 平方根やヒルベルト行列を使用した連立一次方程式の解に関して調査しました。

2.1 除算
逆数がアンダーフロー が発生しない引数範囲
64倍精度 10 4627 未満96倍精度 10 4323 未満
128倍精度 10 4019 未満
で除算(逆数計算)途中で32倍精度演算でアンダーフローの値になると0としています。

整数演算方式 DD 形式の除算の精度基本となる演算精度は32倍精度

引数は 10 n とし逆数 10  n を求める。
逆数計算 精度はm倍精度で表している。
32倍精度演算はn=4931まで32倍精度分の精度を持つ。
n128倍精度 96倍精度 64倍精度
4931
0004625
64
64
64
4626
6464644627
6464644628
3232324629
32
32
32
4321
9696644322
9696644323
969664
4324
64
64
64
4325
6464
64
4017
1289664
4018
12896644019
12896644020
969664
4021
969664

2.2 平方根
平方根の計算方式には以下の2つがある。
(1)逆数平方根計算方式
1x0 a(基本演算精度)

x (1  ax n )x n 1  x n  n
2a  a  xlast21となるので , aの大きさに制限がつく 。
a2ax n  ( ax n ) x n として計算する必要がある。
2xn ≒( 2)除算方式
x0 a (基本演算精度)

x n 1 

1a( xn )
2xna  xlast
この場合 , aの大きさには新たな制 限は加わらない。

整数演算方式 DD 形式の平方根の精度基本となる演算精度は32倍精度
n2

引数は 10 とし , 平方根 10 を求める。
n平方根計算はすべて演算精度分の精度をもつ。

2.3 演算量

1.演算量は32倍精度演算を単位として計算している。
2.平方根計算は逆数平方根算出方式(旧方式)から除算を使用する方式(新方式)に変更しています。
3.絶対値の大きい順に並べる処理では、ＩＦ文を用いない方式を取っています。
演算での注意事項
1.乗算で32倍精度変数を２つに分ける定数はr  2505  12.乗算a  bではアンダーフローを防ぐため,a, bの最初の4倍精度変数の値が0と見なされれば(a(8,1) 0, b(8,1) 0) a  b  0としている。
3.0.5をかける場合,0.5  a(8,1),0.5  a(8,2),
0.5  a(8,3),0.5  a(8,4)と最初の4倍精度のみに0.5をかける。

多倍長精度演算の演算量32倍精度演算を単位とする。

演算
加減算

乗算

除算

平方根
(新)
平方根
(旧)

精度

add,sub32 mult32
div32
64140964501289406418996722412817446643899625848128858
13864120189665196128
2950
414
64182
6396684
1681282304
460

sqrt320000002344612111000000000
11111
1

平方根計算では演算量で見ると新方式は乗算の実行回数が大きく削減されるので、実行時間では旧方式と同じか若干短くなります。

2.4 ヒルベルト行列
ヒルベルト行列の条件数
1(1  i, j  n)i  j 1(n  i  1)!(n  j  1)!
1
逆行列H ji  (1) i  j(i  j  1)[(i  1)!( j  1)!] 2 (n  i )!(n  j )!
1
1
H ijのノルムは1   ...... nが大とすると,
2nln(n)    0.57721566490153286060 (オイラー数)
ヒルベルト行列H ij 1

逆行列H ji のノルムを | H ji

1| の最大値とする。

またnは大として、スターリングの公式n!~2 e nnn

12を使用する。

逆行列の絶対値Pは変形して,(n  i )!(n  j )!i 2 j 2P
(n  i )(n  j )(i  j  1)[i! j!] 2 (n  i )!(n  j )!
i  an, j  bnとすると,P

a 2 b 2 n 2 ((1  a )n)!((1  b)n)!
((a  b)n  1)(1  a )(1  b)[(an)!(bn)!] 2 ((1  a )n)!((1  b)n)!

1 a 1 b
1  a (1 a ) n 1  b (1b ) n 1  a 2 an 1  b 2bn
1 a 1 b

()
()
() ()
2
1 bab(1  a )(1  b)((a  b)n  1)4 1  aab

Pはaとbに対して対称より, a  bとすると、
P

a2
1  a 2 (1 a ) n 1  a 4 an()()a(1  a 2 )(2an  1)4 2 1  a

1  a (1 a ) 1  a 2 a1 a2f (a )  ln[()() ]とするとf ' (a )  ln( 2 )1 aaa0.5  a  0.75ではf (a )はa  0.5  0.7071で最大となる.
i  an, j  bnが整数となる事から, a  0.7とする.
するとPmax 

0.0243 17 2.8 n 17 0.6 n
0.0243 1.530979 n( ) ( )10
1 .4 n  1 7
3
1. 4 n  1

より詳細にすると以下の様になる。
0.02433695097 17 0.6 n 17 2.8 na  0.7
P( ) ( )
1.4n  137
 101.530979068n 1.613733833log10 (1.4 n 1) (1)a  0.5 22P

1( 2n  1)4

(22  2 ( 2 2 ) n)
( 2  1) 2 2n2 2

 101.531102741n 1.596359739log10 ( 2n 1)(2)(2) /(1)  10 0.000123673n  0.017374096 log10 (1.4 n 1) log10 ( 2n 1)
10 m のm(1)n  200n  400

302.1364756
608.0304816(2)
302.1741817
608.0929301(2) /(1)

条件数比率0.03770.0624

1.091
1.155

ヒルベルト行列Hを使用した連立一次方程式Hx=b
では,bの要素に誤差があると条件数が大きいので結果に大きな影響がでます。
その影響をみるためx  (1,1,....,1) T となる様に, bを同じ精度の演算で計算して, Hx  bからxをもとめ, 得たx T と(1,1,....,1) T の値を比較しました。(bには演算精度におおじた誤差が入る.)

これまでのヒルベルト行列実測結果
連立一次方程式は部分軸選択付きのLU分解法を使用しています。
このためbに全く誤差のない値（たとえばbの要素すべて1など)を入れると,得られた解xを使用してHxとｂを比較するとほとんど誤差は出ません.
ヒルベルト行列最大誤差一覧表
次元数倍精度
４倍精度
８倍精度
2051.616
1.44E-07 2.13E-47
40
182.10762.034 8.72E-16
60 1971.569
62.804251.13280642.074
300.734
252.653100
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

詳細なテスト結果
ヒルベルト行列の精度結果のビット数が113以上は誤差なしとして0としている。
有効ビット数 32倍精度 1009,64倍精度 2018
演算精度 32倍精度64倍精度
結果
条件数
条件数
結果
条件数
条件数
次元数
（ビット数) （ビット数) (定義式) 次元数（ビット数) （ビット数) (定義式)
170
0853
857
370
01870
1874
172
0864
867
3720
1880
1884
1740874
878
3740
1890
1894
176
0884
888
3760
1900
1905
1780
894
898
378
106
1910
1915
180105904908
380
97
1920
1925
182
101914918
382
861931
1935
184
8592492838476
1941
1945
186
76
935
939
386
76
1951
1955
188
74
945
949
388
61
1961
1966
190
54
955
959
390
47
1971
1976
192
46
965
969
392
39
1981
1986
194
37
975
979
394
25
1992
1996
196
27985
989
396
16
2002
2006
198
15
996
1000
398
5
2012
2016
200
6
1006
1010
400
-6
2022
2027
202-8
1016
1020
204
-5
1026
1030
206
-71036
1040
208
-6
1046
1050
210
-7
1057
1061
(注)マイナスのビット数は最初の１桁目より結果が異なる事を示す。

ヒルベルト行列Aを使用して,Ax  bでx  (1,1,...,1) T となる様に
nbi   a ij として, bを定めますが条件数はj 1

このbを作成する際の誤差に大きく影響します。
例えばAを32倍精度で求め, その他の演算を全て64,96,128倍精度で行うと条件数の影響は全くでなくなり、解 x  (1,1,...,1) T が得られます。

簡単な例

111 1A 1 A1 1 11
 (1  ）(2   ) 24
条件数 ≒ (1  ） 
解 x  (1,1) T となる様にbを定める場合の値でb0  (2,2  2 ) T , b1  (2   ,2  2 )となる場合があります。
Ax  b0 の解はx  (2,0) T , Ax  b1の解はx  (1,1) T3.量子モンテカルロ法による物性スペクトル計算結果のチェックとして, L  1個の計算結果の絶対値をE(i)(i  0,1,....., L)としたとき、L
E(0) 1, E(L) 1, E( )で絶対値が最小値をとる2という条件とした。
n  20とすると,
1βU L  1P
 e (P )となる。
L
L
2E( )
2β  L,βU  Lの場合,βU βL条件数は, e (P () )となり注意を要する。
Lβ
またβ  Lとなると誤差がO(( )2 )から正しく計算Lする条件はさらに厳しくなる。正しく計算されるとP

1βU L  1P
 e (P )となる。
L
L
2E( )
2

3.1 32,64,96，128倍精度演算の結果の精度
N=20,L=448 結果一覧
32倍精度演算βu

実測値
理論値
10 3.12E-231 4.79E-231u

実測値
理論値10 0.202E-325 0.190E-325
10 0.193E-356 0.155E-356
5 2.45E-273 3.03E-273
5 4.13E-292 4.63E-292
3 3.01E-240 4.41E-240
2 5.69E-207 9.90E-207

u

実測値
理論値
8 0.228E-344 0.194E-344
7 0.228E-344 0.194E-344
4 3.41E-277 4.13E-277
3 3.63E-253 4.98E-253

u

実測値
理論値
9 0.312E-365 0.239E-365
7 0.227E-344 0.194E-344
8 0.407E-368 0.308E-368
4 3.41E-277 4.13E-277
5 0.961E-309 0.985E-309
3 3.63E-253 4.98E-253
4 4.13E-292 4.63E-29225064倍精度演算
β500
600
700
800
900
1000
96倍精度演算
β700
800
900
1000
128倍精度演算
β700
800
800
900
900
1000
1000

実測値と理論値は良く一致している。！

N  30, L  448 での計算式
11L 1Pm
 e ( P  ( xU  (ln( x)) ) )
絶対値最小値m2
  200x  1mL
200

N=30,L=448 の実行結果
64倍精度演算
βu
500
600
700
800
900
1000

実測値
理論値
10 0.630E-328 0.818E-328
10 0.438E-360 0.107E-360
5 2.39E-278 1.96E-278
5 1.21E-298 2.14E-298
3 1.42E-248 5.57E-248
2 2.40E-217 7.29E-217

実測値と理論値は良く一致している。！

32倍精度と64倍精度演算の結果比較N  100, L  448で正しく計算出来る

 ,U32倍精度演算   250, U  10

  300, U  8
64倍精度演算   500, U  10
  600, U  10
64倍精度と96倍精度,128倍精度演算結果で差がでるものは意外に少ない。

32倍精度変数を複数つなげる事による引数範囲の縮小化による。
この後差の出たケースの結果をしめしました。

N  12,   1000,U  10, L  160
0.9750000000D+030.9812500000D+030.9875000000D+030.9937500000D+030.1000000000D+040.1185494857+131
0.9203004545+139
0.1006367236+149
0.1112616274+157
0.1036389481+166

64bai elapse= 118.255023000000

sec

0.9750000000D+03 -0.3077877772D-14
0.9812500000D+03 -0.8349313422D-11
0.9875000000D+03 -0.2264972890D-07
0.9937500000D+03 -0.6191951556D-04
0.1000000000D+04 -0.1000000000D+01
96bai elapse= 334.062215000000

sec

0.9750000000D+03 -0.3077877772D-14
0.9812500000D+03 -0.8349313422D-11
0.9875000000D+03 -0.2264972890D-07
0.9937500000D+03 -0.6191951556D-04
0.1000000000D+04 -0.1000000000D+01
128bai elapse= 652.552798000000

secN  12,   1200, U  10, L  100

0.1152000000D+04 -0.1556030131D-19
0.1164000000D+04 -0.8901756069D-15
0.1176000000D+04 -0.5092527422D-10
0.1188000000D+04 -0.2913410882D-05
0.1200000000D+04 -0.1000000000D+01
96bai elapse= 194.011505000000

sec

0.1152000000D+04 -0.1556030131D-19
0.1164000000D+04 -0.8901756069D-15
0.1176000000D+04 -0.5092527422D-10
0.1188000000D+04 -0.2913410882D-05
0.1200000000D+04 -0.1000000000D+01
128bai elapse= 442.914668000000

secN  8,   1500,U  10, L  96

0.1437500000D+040.1453125000D+040.1468750000D+040.1484375000D+040.1500000000D+040.2067081954+228
0.8661847662+246
0.3114430095+265
0.2193652260+284
0.5614536671+302

96bai elapse= 65.1950880000000

sec

0.1437500000D+04 -0.4821874620D-22
0.1453125000D+04 -0.1293888751D-16
0.1468750000D+04 -0.3471985966D-11
0.1484375000D+04 -0.9316632940D-06
0.1500000000D+04 -0.1005946110D+01
128bai elapse= 121.441539000000

sec3.2 32,64,96,128倍精度演算性能測定結果

測定条件
N=20,L=448
X5570 -O2 –fp-model strict で実行

32倍精度
64倍精度
96倍精度
128倍精度

600秒
6000秒
17052秒
34407秒

110
28.4
57.3

(注)96倍精度はもとになった6倍精度の減算に無駄な演算があったので修正しています。
6倍精度は,x5570,e5430,SR16000/M1,BG/Q,T2Kで性能測定をして、性能が向上したのを確認しています。