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

目次
1. はじめに2. 整数演算方式による128倍精度演算の作成
2.1 概略
2.2 加減算
2.3 乗算2.4 作成作業に関して
3. 128倍精度演算の精度
3.1 平方根計算
3.2 ヒルベルト行列
4. 128倍精度演算の性能
5.量子モンテカルロ法による物性スペクトル計算
5.1 大きなβの場合に関して
5.2 パラメータ領域拡大結果

1.はじめに量子モンテカルロ法による物性スペクトル計算で正しく計算出来るパラメータ範囲の限界を求めるために,32倍精度変数をつなげる方式で,64,96,128倍精度演算を作成して32倍精度演算でのパラメータ範囲の限界値  250
U  10, N  100, L  448の条件でx5570で実行した所,64倍精度演算で 2日弱,96倍精度演算で1週間,128倍精度演算で 2週間かかっています。 最終的には 64倍精度演算で正しく計算出来るパラメータ範囲は  600, U  10, N  100
L  448まで拡大しています。
ieee754  2008形式では表現できる数値範囲から理論的には128倍精度演算まで実行可能ですが,1ケースのテストで 2週間もかかるため,128倍精度演算を整数演算方式で行う様にしました。
ここで問題となるのはソースステップ数の多さです。
原因としては8倍精度演算等の低精度演算では性能などの問題、32倍精度より低精度演算では繰り返しパターンがほとんどないためdo文を使用せず, 演算をすべて展開するという方式を取っていた事によります。
128倍精度演算となると繰り返しパターンも多くなり, do文を使用せず演算を展開してもその性能上の効果もなく, あまり展開が多いと性能上も逆効果となりますので, do文を使用してソースステップ数の削減を図りました。

整数演算方式では正整数演算a  b, a  bが負の値ならない様にするため, 整数型4倍バイト配列は各要素 30ビット使用する,整数型8倍バイト配列は各要素 60ビット使用する様にしています。
また8倍精度変数の仮数部のビット数は 240ビットと低精度演算ではもっとも扱いやすくな っています。このため,8  15  k (kは正整数）
精度演算が同じパターンがでる事になります。また 4n倍精度のほうが無駄な処理が少なくなるため8  60  k倍精度演算が適していることがわかります。k  2とすれば8  60  k  128となり,128倍精度演算が適している事がわかります。
最終的には,1つのサブルーチンが350ステップ内で収まり,加減算, 乗算のルーチンでは共通する部分がありため、デバッグ作業等も含めても,1人で扱える作業量となりました。
ヒルベルト行列の次元数800(条件数101200 )で10進8桁まで一致しました。
これから,128倍精度演算から188倍精度演算,248倍精度演算にするのには定数を変更するだけで良く、実際デバッグ ( N  100のヒルベルト行列の誤差を求める)を含め,188倍精度は 2時間、これで慣れた事により248倍精度演算は1時間で作業が終了しました。
この事から8  60  k (k  0,1,2,3,4,5,6,7,8)倍精度演算を作成し,加減算, 乗算, 除算, 平方根計算, 判定( gt , eq, lt )の合計のステップ数
はコメント , 定義文等を含め, k  0で1000ステップ, k  0で1500ステップとなっています。今回は128倍精度演算を中心に必要におおじて他の精度演算にかんしても記述しました。

2.整数演算方式による128倍精度演算の作成
2.1概略
除算,平方根計算は以下の様に反復法を用いる事で加減算,乗算の作成に帰着されます。
平方根計算
f ( x )  x 2  a, f ' ( x )  2 xx a 1ax n 1  x n  n
 ( xn )
2 xn
2xn
21aa( xn )
2xnxn

最終段階x n 1 

初期値x 0 a ( 4倍精度).除算a1 a  で計算します。
bb
逆数計算は
c

1
1
 b, f ' ( x )   2xx
1bxnx n 1  x n  x n  x n (1  bx n )  x n ( 2  bx n )1f ( x) xn2で行います。
1( 4倍精度)b(注) x n ,1  bx n が小さい場合, x n (1  bx n )の計算でアンダー
初期値ｘ 0 

フロー発 生の可能性があるのでx n 1  x n ( 2  bx n )として います。

2.2 加減算
加減算は引数を整数型8バイトの60ビットのみ使用する配列に展開します。準備す る配列要素の数は 4080  60  1  69となります。引数を正の整数の演算になる様に引数を組み替えます。
（この処理は加算及び減算に共通です.。.)この処理は低精度演算のものを定数を変更するだけですので前処理サブルーチンに分けます。（ソースステップ数約150)
隠れビットと丸めよう のビットのために、先 頭配列ix (69)の値は加算 2  ix (69)  3

減算 4  ix (69)  7となります。

引数の大きい方を配列ix, 小さい方を配列iyにいれます。
そして指数部の差分だけ配列iyをシフトして, 桁合わせをします。
ここまでは加算及び減算に共通です。ビット 位置が1つことなるのでサブルーチンとしての抜出はしていません。
多倍長整数加算または乗算を行い結果を配列izにいれます。
この段階では整数型8バイト配列izは各要素60ビットしか使用いません.ここで指数部と丸め処理をします。
加算の場合iz (69)  4 桁上がりあり, iz (69)  3

桁上がりなし減算の場合iz (69)  4 桁下がりなし,2  iz (69)  3 桁下がり1ビット
このケースの場合, ともに丸め処理後, 指数部の値と仮数部を配列izの68要素にいれます。処理 がことなるのは, 減算でiz (69)  1の場合は桁下がり2ビットで丸め処理は有りません。iz (69)  0の場合はiz (68), iz (67).....と調べて最初に1があるビット位置を求める処理がある事です。これはiz (l )  0, iz (l )  2 m  k (1  k  2 m  1)を実数にするとmが指数部の値に一致する事を使用します。これで指数部の値と仮数部を配列izの68要素にいれます。（ iz (69)  0の場合は前につめ、下位の配列要素には 0をいれます。）。この後は加減算, 乗算共通で60ビット 68要素配列を64ビット 64要素配列に詰め,128倍精度浮動小数点数を作成します。とくにここもサブルーチン化して抜き出さなくてむ、コメントなども含めても約350ステップ数内に収まっています。

2.3 乗算
2つの引数から積の符号 と指数部の値を計算して引数を4080ビットの正整数にしま す。これを整数型 4バイトの30ビットのみ使用する配列に展開します。準備する配列要素の数は4080  30  1  137となります。二つの引 数を配列ix, iyにいれ多倍長正整数乗算を行い, 結果を整数型４倍バイトの30ビット使用する配列izにいれます.(izの要素数は 274)。
ix (137), iy (137)には隠れビット1をいれますので , 桁上がりがあるかどうかはiz ( 273)の値で求まり, これで丸めビットの位置がわかります。
iz ( 273)  2だと桁上がりがあり, iz ( 272)  1だと桁上がりが有りません。
ここで整数型8バイト配列izz（配列要素数69)を用意して, iz   izzにいれ, 加減算と同じ処理をして,128倍精度浮動小数点数を作成します。これらを含めて300ステップ以下で作成できています。

2.4 作成作業に関して整数演算方式の作成の作業量は連続して使用できるビット数によります。例えば8192ビット使用できるとすれば、ソースステップ数は128倍精度演算は20-30程度で行う事が可能と言えます。
また丸め処理は行わない(例えば切り捨て)とすればやはりソースステップ数が大きく削減されます。FPGA等で作成する場合もほとんど同じ様になります。

3.128倍精度演算の精度

量子モンテカルロ法による物性スペクトル計算で演算精度に影響があるのは、逆行列を求める計算とグラムシュミット法の計算です。
そこでヒルベルト行列を使用した連立一次方程式の求解と平方根計算の精度を求めました。

また128倍精度演算では,量子モンテカルロ法による物性スペクトル計算では表現可能な数値範囲の制限にかかる可能性があり,問題発生時にその原因が有効ビット数不足なのか、表現可能な数値範囲の制限によるものかの切り分けのため、ヒルベルト行列では188倍精度演算,248倍精度演算,308倍精度演算,368倍精度演算,428倍精度演算,488倍精度演算でも実行しています。

3.1 平方根計算
テスト内容は以下のものです。

a  10 (n  4931,4932)
n1d
ab ac dsqrt div limit test
a= 1.000000000000000000000000000000000E+4931
b= 3.162277660168379331998893544432718E-2466
c= 9.999999999999999999999999999999999E-4932
d= 3.162277660168379331998893544432718E+2465

a  10 4931 では、除算, 平方根は正しく計算できています。

sqrt div limit test
a= 9.999999999999999999999999999999999E+4931
b= 0.000000000000000000000000000000000E+0000
c= 0.000000000000000000000000000000000E+0000
d= 9.999999999999999999999999999999999E+2465

a  10 4932 の場合, 除算が正しく実行されていません。
これは指数部の値が 0の場合には 0としている事によります。
加減算, 乗算などのアンダーフ ロー処理を容易にするために行っています。このため最小値の10の指数は
 16382  log10 2  4931.473389......となります3.2 ヒルベルト行列

ヒルベルト行列のサイズNと条件数には近似的には以下の様になっています。
ヒルベルト行列の条件数n条件数(ビット数表示)
400
2022
450
2276
500
2530
550
2785
600
3039
650
3293
700
3547
750
3801
800
4056
850
4310
900
4564
950
4818
1000
5072
1050
5327
1100
5581
1150
5835
1200
6089
1250
6344
1300
6598
1350
6852
1400
7106
14507361
1500
7615

128倍精度
有効ビット数 4081

188 倍精度
有効ビット数 6001

248 倍精度
有効ビット数7921

実測結果は以下の様になっています。
ヒルベルト行列の精度

N

最大誤差
400
500
600
700
800
1000
1100
1400

精度

1.014E-621 128倍精度
1.681E-467 128倍精度
1.321E-314 128倍精度
1.349E-161 128倍精度
7.603E-009 128倍精度
1.013E-280 188倍精度
4.211E-127 188倍精度
1.013E-245 248倍精度

最大誤差のビット数  有効ビット数  条件数のビット数となっています。
248倍精度演算ではN  1400が使用できるメモリ容量の上限となっています。

より高精度演算での精度

N=100でのヒルベルト行列計算
条件数0.127E+152
精度

ビット数501

有効ビット数 誤差
実測
理論値
308
9841 4.73E-2814 2.80E-2812
36811761 3.38E-3392 2.52E-3390
428
13681 1.61E-3971 2.66E-3968
488
15601 5.36E-4198 2.80E-4546 (注)488
15601 1.00E-4547 2.80E-4546(注)逆数計算で4倍精度の初期値から7回反復した場合の値
308,368,428倍精度は7回反復,488倍精度は8回反復が必要

488倍精度以上の精度を作成してもこれ以上精度向上はないと言えます。

4.128倍精度演算の性能比較のため行列サイズN=400のヒルベルト行列による連立一次方程式の求解で実施しました。
ヒルベルト行列のN=400の実行時間一覧表

方式
浮動型
浮動型
浮動型
整数型
整数型
整数型
整数型

精度
64倍精度
96倍精度
128倍精度
68倍精度
128倍精度188倍精度248倍精度

単位: 秒x5570
e5430
540.603 743.938
1551.696 2165.918
3467.115 4280.500
179.059 277.818
652.896 918.142
1235.567 1945.560
2125.342 3337.811

32倍精度変数をつなげる方式に比べ,整数演算方式の性能が非常に良くなっています。

他の精度での性能
N=100でのヒルベルト行列計算

精度

実行時間 (秒)
x5570
e5430
308
50.29881.123368
70.468115.572428
95.133156.090488
123.216205.238
反復回数7回
488
123.253205.404

N=400 での8倍精度でのヒルベルト行列計算
Ax=b でＡの要素は4倍精度演算で作成

ソース
従来
新規

実行時間 (秒)
x5570
e5430
4.686
5.926
4.885
6.354

1.反復回数の増加の影響はほとんどない.
2.8倍精度演算ではソースステップ数の大幅削減をしても性能低下は5%程度に収まっています。

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

5.1 大きなβの場合に関して量子モンテカルロ法による物性スペクトル計算
ではB L ......B1  QDRとし(I  QDR ) 1 を求めるのに,D  Q T R 1  Q' D' R ' とQDR分解し,Q(D  Q T R 1 )R  QQ' D' R ' R  Q' ' D' ' R ' '  I  QDRから逆行列を求めます。
U L  1P  e (Q ())とすると,
L2Q

1
U  Lなら結果の絶対値の最小値は ,P
Dの最大数はP 2 となり, グラムシュミット法で平方数の和の平方根を求めるため,プログラムに現れる最 大数はP 4 となります。

U  L, では条件数は / Lの大きさ依存する度合いが大きくなる。
これは, B L ......B1  QDRとし(I  QDR ) 1 を求める際に発生します。
D  Q T R 1  Q' D' R ' とQDR分解し,Q(D  Q T R 1 )R  QQ' D' R ' R  Q' ' D' ' R ' '  I  QDRから逆行列を求めます.
このD  Q T R 1のQDR分解の際オーバーフロ ーが発生します.
正確には, b  a 1  a 2  ......  a n の計算の際
222a i  a i の計算で発生します。
オーバーフローを防ぐ には,a max  絶対値最大値(a i , i  1,2,...., n )として, b | a max |  c1  c 2  .....  c n222(c i ai)
a maxとします。ここで正しい答えのでる最大数における, D' ' の最大数Ｓ,S  e UL  としてその影響を調査しました。
調査ではn  20, U  10, L  448と固定しています。

調査結果条件 N  20, U  10, L  448
精度
β の最大数 最大値
α68倍精度732 2.518D+1287
128倍精度
1620 7.778D+2441
188倍精度2590 6.198D+3593
248倍精度
3609 6.580D+4745
308倍精度3776 6.591D+4930

1.64
2.09
2.43
2.72
2.76

368倍精度では308倍精度と同じく、β=3776では正しく計算出来ていて,β=3777では結果が不正となっている事から308倍精度の結果が限界を示していると言えます。

368倍精度演算で正しく計算されている例
n=20
l=448
beta, dt= 3776.000000000000000000000000000008.42857142857142857142857142857143
u, t0= 10.0000000000000000000000000000000
1.00000000000000000000000000000000
fac= 9.18072515031978755045201686537202
For denominator: mstep=
1
For green func: mstep2=
1
ispin=1
QDR = 6.993199813124998895732748089737576E+4922
I+QDR= 6.993199813124998895732748089737576E+4922

ispin=2
QDR = 6.591017770912797832925133116465994E+4930
I+QDR= 6.591017770912797832925133116465994E+4930

0.0000000000D+00 -0.1000000000D+01
0.1879571429D+04 -0.7381849773-890
0.1888000000D+04 -0.1520746834-893
0.1896428571D+04 -0.7381849773-890
0.3776000000D+04 -0.1000000000D+01

368倍精度演算でも正しく計算されてない例
n=
20
l=
448
beta, dt= 3777.00000000000000000000000000000
8.43080357142857142857142857142857
u, t0= 10.0000000000000000000000000000000
1.00000000000000000000000000000000
fac= 9.18194073789880841037203642516443
For denominator: mstep=
1
For green func: mstep2=1

ispin=1
QDR = 8.886250347365217688103022558980503E+4923
I+QDR= 8.886250347365217688103022558980503E+4923

ispin=2QDR = 8.395584469058393694708873673477994E+4931
I+QDR= 3.111402625840683560090283433225110E+4770
0.0000000000D+00 -0.4999999997D+00
0.1880069196D+04 -0.4586860821-890
0.1888500000D+04 -0.7614617438-894
0.1896930804D+04 -0.2814546226-890
0.3777000000D+04 -0.5000000000D+00

N=20,U=10,L=448と固定すると正しく計算出来るβの最大値は3776でα=2.76であった。
ここでN=20,U=10でＬを小さくした場合正しく計算できた場合は以下の表の様になっています。
n=20,U=10 での実行結果一覧
βL3776
4000
4250
4500

最大値
448 6.571D+4930
320 2.857D+4931
208 6.161D+4919
132 1.077D+4927√β ＵＬα1.730D+1786
6.014D+1553
1.780D+1291
2.934D+1058

2.76
3.18
3.814.66

この表から が3776以上になると,  が 250増加するとLは0.64倍となる。  3776と4000では が 224の増加
なのでLは 250 / 224  0.64  0.714倍となっています。
１ , U 1 , L1 ,  1が既知として1   2 の２ を選ぶと, L 2が推定され実行するU 2 を決めれば (通常  U 1 )とすると
 ２  1

１ U 1 L1
から 2 が求まり実行可能かどうかがわかる。
2 U 2L2

  4750, L  84,   5.686

最大値  10 4932.6

  5000, L  54,   6.912

最大値  10 4932.5と計算が困難になる事がわかります。

5.2 パラメータ領域拡大結果
測定結果一覧

x5570

測定条件 N  100, U  10, L  448
と絶対温度Kとの関係
絶対温度(K)  10000 / β

精度

絶対値最小値
実測値
理論値
730 68倍精度 0.824D-384 0.269D-383
1600 128倍精度 0.170D-582 0.214D-582
2500 188倍精度 0.773D-728 0.462D-728
10000/3 248倍精度 0.292D-840 0.978D-841
3770 308倍精度 0.156D-893 0.397D-894

実行時間(秒)
57642
185097403155673631975735

この条件では,32倍精度  250(絶対温度40K )から,68倍精度で13.7 K,128倍精度で 6.26K
188倍精度で 4K,248倍精度で3K,308倍精度
で 2.65Kまでパラメータ領域が拡大されました。
これ以上の演算精度を用いても改善はされませんのでこの値が限界と言えます。