1. No category

整数演算方式による超多倍長演算
日時：平成26年3月7日(金) 1１:00 ～１１：３０
会場：工学院大学 新宿校舎 28階 第4会議室
高エネルギー加速器研究機構
濱口信行
1.はじめに
2.整数演算方式の多倍長演算の作成方法
2.1 加減算
2.2 乗算
2.3 具体例
3.ヒルベルト行列Hを用いた連立一次方程式の求解の例
4.量子モンテカルロ法による物性スペクトル計算への適用例
5.おわりに
1.はじめに
整数演算方式とは浮動小数点演算を符号部,指数部,仮数部に
分け整数演算、論理演算およびシフト演算を使用して行う方式
で、ＦＯＲＴＲＡＮ言語で作成しています。
倍精度変数を複数個つなげて,既存の倍精度演算を使って4倍
精度,8倍精度浮動小数点演算が行われているのに、わざわざ
整数演算方式を取る理由とは？
現在の倍精度のデータ形式は1985年に制定されたもので、
当時の最高性能の計算機が1GFLOPsに届くかどうか,メモリは
仮想空間が16MBから2GBに移行中で,TOP500の前身の
リンパックレポートではチューニングを競うサイズは
1000元(所要メモリ8MB)の時代です。
その時代から約３０年近く経っていますのでプログラムの規模も、
内容も大きくかつ複雑になり、当時は充分な仕様でもやはり問題
がでてきます。
問題は２つあり,仮数部のビット数不足,指数部のビット数不足
（表現できる数値範囲の狭さ）がありますが取られている対策
は仮数部のビット数不足のみです。4倍精度データに関し
2008年に新しく規定が改定されているにも関わらず
(指数部のビット数が11ビットから15ビットに拡張)多くの
計算機でハードウエア命令がサポートされていません。
この様な理由から整数演算方式を使用するようにしました。
作成自体が難しいと考えられるかもしれませんが、私自身
特にFORTRAN言語の教育は学生時代の教養課程で森口繁一
先生のJIS FORTRAN7000を受けた程度でメーカー時代は
特にありませんので皆さんの基礎知識があれば充分こなせる
ものと思います。
私自身が既存の方式で経験した問題は以下の様なものです。
1.ループ積分で仮想光子(本来質量0のものに計算上微小な質量
を与える)の 質量の値を小さくすると,精度の良い結果が
得られない。
2.全球スペクトルモデルで高解像度になると、ルジャンドル陪関数の
テーブル 作成や極付近の計算精度に支障をきたした。
3.量子モンテカルロ法による物性スペクトル計算を解析の際
（プログラムは高エネルギー加速器研究機構 岩野氏より提供
されました）
日常温度やそれ以下の温度で計算する必要があるのに、357Kと
いう高温で表現できる数値範囲の制限によりオーバーフローが
発生した。
更に細かいところでは,ヒルベルト行列の条件数を求めるのにサイズ
Nが200近くなると計算出来なかった等々があります。
2.整数演算方式の多倍長演算の作成方法
加減算,乗算,除算,平方根計算を考えた場合,
除算は逆数と乗算に分け,逆数計算の初期値を
平方根計算は初期値を既存の4倍精度演算を
使用して反復法で求める事ができます。
この事から整数演算方式の多倍長演算作成は
加減算,乗算の作成に帰着されます。
加減算,乗算は符号部と指数部と仮数部に
分けそれぞれ整数演算、論理演算,シフト演算を使用します。
注意する点は、以下の１点のみです。
整数演算では,正または0の演算にして
結果も正または0になる様にする必要があります。
多倍長整数演算(大きな整数の演算 )は, 既に市販の書物
にも書かれています。 今回は整数型 4バイト変数 (i 4)の30
ビット, 整数型8バイト変数 (i8)の60ビットを使用するように
しています。これから 仮数部のビット数が 60の倍数となる
精度演算作成が適している事になります。
作成に必要な共通処理としては
1.引数を30ビット,60ビット配列に展開する 処理。
2.30ビット,60ビット配列に展開され た結果から浮動小数点
数を作成する。
3.丸め処理
があり,
加減算ではc  a  bをic   (ia  ib)に分ける
処理と桁合わせ処 理
更に減算では ia  ibとなる様にするという条件と、
桁落ち数を求める処理が加わります。
2.1加減算
浮動小数点数は符号部1ビット,指数部15ビット,
仮数部mビットとします。
まず多倍長整数演算を行うための場合分けと,指数部の差を
求めます。
多倍長整数演算では丸めを考慮し,加算m+1ビット,減算m+2ビット
演算を行います。
２つの浮動小数点数a,bの加減算は,正整数(a>b) のa+b,a-bに
帰着されます。また,a=0,b=0なら,符号変換等で演算できます。
a=bの場合は,指数部の更新及び,0とする事ができます。
加減算の場合分け
a,bの条件
a>b>0
b>a.>0
a>0>b,a>|b|
a>0>b,a<|b|
b>0>a,b>|a|
b>0>a,b<|a|
0>a>b
0>b>a
ia
ia=a
ia=b
ia=a
ia=-b
ia=b
ia=-a
ia=-b
ia=-a
ib
ib=b
ib=a
ib=-b
ib=a
ib=-a
ib=b
ib=-a
ib=-b
c=a+b
ia+ib
ia+ib
ia-ib
ー(ia-ib)
ia-ib
ー(ia-ib)
ー(ia+ib)
ー(ia+ib)
c=a-b
ia-ib
ー(ia-ib)
ia+ib
ia+ib
ー(ia+ib)
ー(ia+ib)
ia-ib
ー(ia-ib)
ここで引数を60ビットのみ使用する整数型8バイトの配列に入れます。
多倍長整数演算は8バイト整数演算を行います。
c  a ' b' でc  2 m  2 なら桁上がり1。
c  2 m  2 なら桁上がりなし 。
c  a ' b' でc  2 m  2 なら桁下がりなし 。
2 m 1  c  2 m  2 なら桁下がり1.
2 m  c  2 m １ なら桁下がり 2
c  2 m .なら何桁桁下がりがあるか計算する
必要があ ります。(左から最初の1が現れるビット位置を求めます)
2.2.乗算
ここで引数を30ビット使用する整数型4バイトの配列に入れ
ます。
a , bは 2 m  仮数部として多倍長整数乗算を行います。
4バイト整数を8バイト整数に代入して8バイト整数乗算
を行います。（整数乗算の結果が負にならないため)
c  a  b　でc  2 2 m 1 なら桁上がり1.c  2 2 m 1 なら桁上がりなし。
指数部の計算は符号部を除いた15ビットの整数加算で行います。
符号部はa , bの符号部のexclusive　orを取ります。
2.3 具体例
まず使用したメモを紹介します。規則性がみてとれると思います。
超多倍長精度演算
作成メモ
精度
仮数部ビット数
R16宣言
I8宣言
I4宣言
i4*i4宣言
i4*i4先頭
i4*i4最下位
用途
精度の保持
呼ぶ側の宣言
加減算用
乗算用
整数乗算
桁上がりの有無
丸め位置
8
240
a(2)
ix(5)
ix(9)
iz(18)
iz(17)
iz(9)
68
2160
a(17)
ix(37)
ix(73)
iz(146)
iz(145)
iz(73)
精度　8  60  k　(k  0,1,2,3,4,5,6,7,8)
仮数部のビット数　精度  32  16
R16宣言数  精度 / 4　kとk  1の差15
I8宣言数  仮数部のビット数 / 60  1　kとk  1の差32
I4宣言数  仮数部のビット数 / 30  1　kとk  1の差64
I4  I4宣言数  I4宣言数  2　kとk  1の差128
I4  I4先頭  I4  I4宣言数  1　kとk  1の差128
I4  I4最下位  I4宣言数　kとk  1の差64
128
4080
a(32)
ix(69)
ix(137)
iz(274)
iz(273)
iz(137)
188
6000
a(47)
ix(101)
ix(201)
iz(402)
iz(401)
iz(201)
248
7920
a(62)
ix(133)
ix(265)
iz(530)
iz(529)
iz(265)
308
9840
a(77)
ix(165)
ix(329)
iz(658)
iz(657)
iz(329)
368
428
488
11760
13680
15600
a(92)
a(107)
a(122)
ix(197) ix(229) ix(261)
ix(393) ix(457) ix(521)
iz(786) iz(914) iz(1042)
iz(785) iz(913) iz(1041)
iz(393) iz(457) iz(521)
次ぎに128倍精度演算作成時のメモを紹介します。
加算の配列代入
integer*4 a(128) =>integer*8 ix(69)
ix(69)
ix(68)
ix(67)
ix(66)
ix(65)
ix(64)
ix(63)
ix(62)
ix(61)
ix(60)
ix(59)
2+a(128)(15:15)
a(128)(0:14) a(127)(17:31)
a(126)(0:18) a(125)(21:31)
a(124)(0:22) a(123)(25:31)
a(122)(0:26) a(121)(29:31)
a(120)(1:30)
a(119)(0:2) a(118)(5:31)
a(117)(0:6) a(116)(9:31)
a(115)(0:10) a(114)(13:31)
a(113)(0:14) a(112)(17:31)
a(111)(0:18) a(110)(21:31)
a(127)(0:16) a(126)(19:31)
a(125)(0:20) a(124)(23:31)
a(123)(0:24) a(122)(27:31)
a(121)(0:28) a(120)(31:31)
a(120)(0:0) a(119)(3:31)
a(118)(0:4) a(117)(7:31)
a(116)(0:8) a(115)(11:31)
a(114)(0:12) a(113)(15:31)
a(112)(0:16) a(111)(19:31)
a(110)(0:20) a(109)(23:31)
ix(66)からix(59)までのパターンを7回繰り返す。
ix(58) a(109)(0,22)
ix(50) a(94)(0,22)
ix(42) a(79)(0,22)
ix(34) a(64)(0,22)
ix(26) a(49)(0,22)
ix(18) a(34)(0,22)
ix(10) a(19)(0,22)
あまり部分
ix(2) a(4)(0:22) a(3)(25:31) a(3)(0:24) a(2)(27:31)
ix(1) a(2)(0:26) a(1)(29:31) a(1)(0:28)
8倍精度演算では繰り返し部分を削除し前後をつなげれば良い。
ix(5)
ix(4)
ix(3)
ix(2)
ix(1)
2+a(8)(15:15)
a(8)(0:14) a(7)(17:31)
a(6)(0:18) a(5)(21:31)
a(4)(0:22) a(3)(25:31)
a(2)(0:26) a(1)(29:31)
a(7)(0:16) a(6)(19:31)
a(5)(0:20) a(4)(23:31)
a(3)(0:24) a(2)(27:31)
a(1)(0:28)
減算の配列代入
integer*4 a(128) =>integer*8 ix(69)
ix(69)
ix(68)
ix(67)
ix(66)
ix(65)
ix(64)
ix(63)
ix(62)
ix(61)
ix(60)
ix(59)
4+a(128)(14:15)
a(128)(0:13) a(127)(16:31)
a(126)(0:17) a(125)(20:31)
a(124)(0:21) a(123)(24:31)
a(122)(0:25) a(121)(28:31)
a(120)(0:29)
a(119)(0:1) a(118)(4:31)
a(117)(0:5) a(116)(8:31)
a(115)(0:9) a(114)(12:31)
a(113)(0:13) a(112)(16:31)
a(111)(0:17) a(110)(20:31)
a(127)(0:15) a(126)(18:31)
a(125)(0:19) a(124)(22:31)
a(123)(0:23) a(122)(26:31)
a(121)(0:27) a(120)(30:31)
a(119)(2:31)
a(118)(0:3) a(117)(6:31)
a(116)(0:7) a(115)(10:31)
a(114)(0:11) a(113)(14:31)
a(112)(0:15) a(111)(18:31)
a(110)(0:19) a(109)(22:31)
ix(66)からix(59)までのパターンを7回繰り返す。
ix(58) a(109)(0,21)
ix(50) a(94)(0,21)
ix(42) a(79)(0,21)
ix(34) a(64)(0,21)
ix(26) a(49)(0,21)
ix(18) a(34)(0,21)
ix(10) a(19)(0,21)
あまり部分
ix(2) a(4)(0:21) a(3)(24:31) a(3)(0:23) a(2)(26:31)
ix(1) a(2)(0:25) a(1)(28:31) a(1)(0:27)
8倍精度演算では繰り返し部分を削除し前後をつなげれば良い。
ix(5)
ix(4)
ix(3)
ix(2)
ix(1)
4+a(8)(14:15)
a(8)(0:13) a(7)(16:31)
a(6)(0:17) a(5)(20:31)
a(4)(0:21) a(3)(24:31)
a(2)(0:25) a(1)(28:31)
a(7)(0:15) a(6)(18:31)
a(5)(0:19) a(4)(22:31)
a(3)(0:23) a(2)(26:31)
a(1)(0:27)
乗算の配列代入
integer*4 a(128) =>integer*4 ix(137)
ix(137)
ix(136)
ix(135)
ix(134)
ix(133)
ix(132)
ix(131)
ix(130)
ix(129)
ix(128)
ix(127)
ix(126)
ix(125)
ix(124)
ix(123)
ix(122)
ix(121)
1
a(128)(0:15)
a(127)(0:17)
a(126)(0:19)
a(125)(0:21)
a(124)(0:23)
a(123)(0:25)
a(122)(0:27)
a(121)(0:29)
a(120)(2:30)
a(120)(0:1)
a(119)(0:3)
a(118)(0:5)
a(117)(0:7)
a(116)(0:9)
a(115)(0:11)
a(114)(0:13)
a(127)(18:31)
a(126)(20:31)
a(125)(22:31)
a(124)(24:31)
a(123)(26:31)
a(122)(28:31)
a(121)(30:31)
a(119)(4:31)
a(118)(6:31)
a(117)(8:31)
a(116)(10:31)
a(115)(12:31)
a(114)(14:31)
a(113)(16:31)
ix(136)からix(121)までのパターンを7回繰り返す。
ix(120)
ix(104)
ix(88)
ix(72)
ix(56)
ix(40)
ix(24)
a(113)(0:15)
a(98)(0:15)
a(83)(0:15)
a(68)(0:15)
a(53)(0:15)
a(38)(0:15)
a(23)(0:15)
あまり部分
ix(8)
ix(7)
ix(6)
ix(5)
ix(4)
ix(3)
ix(2)
ix(1)
a(8)(0:15)
a(7)(0:17)
a(6)(0:19)
a(5)(0:21)
a(4)(0:23)
a(3)(0:25)
a(2)(0:27)
a(1)(0:29)
a(7)(18:31)
a(6)(20:31)
a(5)(22:31)
a(4)(24:31)
a(3)(26:31)
a(2)(28:31)
a(1)(30:31)
8倍精度演算では繰り返し部分を削除し,ix(9)=1とあまり
の部分となります。
桁上がりと丸め位置
加減算：integer*8 iz(69)
加算:iz(69)が4以上 桁上がりあり。
仮数部iz(1)(2)~iz(69)(1) 丸め位置 iz(1)(1)
加算:iz(69)が3 以下 桁上がりなし。
仮数部iz(1)(1)~iz(69)(0) 丸め位置 iz(1)(0)
減算:iz(69)が4以上
仮数部iz(1)(2)~iz(69)(1)
減算:iz(69)が2または3
仮数部iz(1)(1)~iz(69)(0)
減算:iz(69)が1
仮数部iz(1)(0)~iz(68)(59)
桁下がりなし。
丸め位置 iz(1)(1)
桁下がりあり。1ビット
丸め位置 iz(1)(0)
桁下がりあり。2ビット
丸めなし。
別処理が必要
減算:iz(69)が0 桁下がりあり。
仮数部iz(1)(1)~iz(68)(59) 可変
丸めなし
乗算
integer*4 ix(137),iy(137),iz(274)
ix(i)*iy(j)
iz(i+j) 上位30ビット
iz(I+j-1) 下位30ビット
iz(274)は常に0となる。
乗算:iz(273)が2以上 桁上がりあり。
仮数部iz(137)(1)~iz(273)(0)
丸め位置 iz(137)(0)
乗算:iz(273)が2未満 桁上がりなし。
仮数部iz(137)(0)~iz(272)(29)
丸め位置 iz(136)(29)
桁合わせ
i=l/60
J=l-60*i
i=0 nooperation
i≧1
do m=1,69-i
iy(m)=iy(m+i)
end do
do m=69-i+1,69
iy(m)=0
end do
ｌ：指数部の差から得られたシフト数
iand の引数はinteger*8で,ishftの引数はinteger*4で
J=0 nooperation
j ≧1
do m=1,68-i
iy(m)=ishft(iy(m),-j)+ishft(iand(iy(m+1),2**j-1),60-j)
end do
iy(69-i)=ishft(iy(69-i),-j)
減算でiz(69)=0の場合
iz (l ) 1  l  68　で最初に1となるビット位置m
桁下がり数　2  (68  l )  60  (60  m)
l  68　nooperatio n
l  67
do　n  l ,1,1
iz (n  68  l )  iz ( n)
end　do
do　n  1,68  l
iz (n)  0
end　do
iz(68)=iand(iz(68),2**m-1)
do n=68,2,-1
iz(n)=ishft(iand(iz(n),2**m-1),60-m)
+iand((ishft(iz(n-1),-m),2**(60-m)-1)
end do
iz(1)=ishft(iand(iz(1),2**m-1),60-m)
整数型8バイトの詰め込み
integer*8
iz(68)=>iz(64)
以下のマスクを準備しておく。
integer*8 ip4,ip8,ip12,ip16,ip20,ip24,ip32,
c
ip36,ip40,ip44,ip48,ip52,ip56,ip60
data
data
:
:
data
data
ip4/z000000000000000f/
ip8/z00000000000000ff/
ip56/z00ffffffffffffff/
ip60/z0fffffffffffffff/
乗算の場合は結果の仮数部をinteger*4 izz(136)に入
れた後、2要素組みでinteger*8 iz(68)に入れておきます。
iz(1)=ior(iz(1),ishft(iand(iz(2),2**4-1),60))
iz(2)=ior(ishft(iz(2),-4),ishft(iand(iz(3),2**8-1),56))
iz(3)=ior(ishft(iz(3),-8),ishft(iand(iz(4),2**12-1),52))
iz(4)=ior(ishft(iz(4),-12),ishft(iand(iz(5),2**16-1),48))
iz(5)=ior(ishft(iz(5),-16),ishft(iand(iz(6),2**20-1),44))
iz(6)=ior(ishft(iz(6),-20),ishft(iand(iz(7),2**24-1),40))
iz(7)=ior(ishft(iz(7),-24),ishft(iand(iz(8),2**28-1),36))
iz(8)=ior(ishft(iz(8),-28),ishft(iand(iz(9),2**32-1),32))
iz(9)=ior(ishft(iz(9),-32),ishft(iand(iz(10),2**36-1),28))
iz(10)=ior(ishft(iz(10),-36),ishft(iand(iz(11),2**40-1),24))
iz(11)=ior(ishft(iz(11),-40),ishft(iand(iz(12),2**44-1),20))
iz(12)=ior(ishft(iz(12),-44),ishft(iand(iz(13),2**48-1),16))
iz(13)=ior(ishft(iz(13),-48),ishft(iand(iz(14),2**52-1),12))
iz(14)=ior(ishft(iz(14),-52),ishft(iand(iz(15),2**56-1),8))
iz(15)=ior(ishft(iz(15),-56),ishft(iand(iz(16),2**60-1),4))
iz(1)～iz(15)のパターンを3回繰り返す。
iz(61)=ior(iz(65),ishft(iand(iz(66),2**4-1),60))
iz(62)=ior(ishft(iz(66),-4),ishft(iand(iz(67),2**8-1),56))
iz(63)=ior(ishft(iz(67),-8),ishft(iand(iz(68),2**12-1),52))
iz(64)=iand(ishft(iz(68),-12),2**48-1)
準備したマスクを
使用する
3.ヒルベルト行列Hを用いた連立一次方程式の求解の例
1
ヒルベルト行列H ij 
(1  i, j  n )
i  j 1
を用いてHx  b　x  (1,1,......,1)
T
となる様にbを定めこの方程式を解き
得られたHy  bから max
(
abs
(
y
(
i
)

1
))
を
1i  n
T
求めます。
nは大と仮定して条件数を以下の近似式
として結果をチェックしています。　ヒルベルト行列の条件数
1
ヒルベルト行列H ij 
(1  i, j  n )
i  j 1
(n  i  1)!(n  j  1)!
1
i j
逆行列H ji  (1)
(i  j  1)[(i  1)!( j  1)!] 2 (n  i)!(n  j)!
1
1
H ijのノルムは1   ......  nが大とすると,
2
n
ln(n )  　  0.57721566490153286060　(オイラー数 )
1
逆行列H ji のノルムを | H ji
1
| の最大値とします。
またnは大として、スターリングの公式n!~
2 e  n n
逆行列の絶対値Pは変形して,
(n  i)!(n  j)!i 2 j2
P
(n  i)(n  j)(i  j  1)[i! j!] 2 (n  i)!(n  j)!
n
1
2
を使用します。
i  an, j  bnとすると,
a 2 b 2 n 2 ((1  a )n )!((1  b)n )!
P
((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 2 bn
1 a 1 b

(
)
(
)
(
) (
)
2
1 b
a
b
(1  a )(1  b)((a  b)n  1)4 1  a
ab
Pはaとbに対して対称より, a  bとすると、
a2
1  a 2 (1a ) n 1  a 4 an
P
(
)
(
)
2
2
a
(1  a )(2an  1)4 1  a
1  a (1a ) 1  a 2 a
1 a2
f (a )  ln[(
)
(
) ]とするとf ' (a )  ln( 2 )
1 a
a
a
0.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.4n  1 7
3
1.4n  1
ヒルベルト行列の条件数
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
1450
7361
1500
7615
128倍精度
有効ビット数 4081
188 倍精度
有効ビット数
6001
248 倍精度
有効ビット数　7921
実測結果は以下の様になっています。
ヒルベルト行列の精度
N
最大誤差
400
500
600
700
800
1000
1100
1400
1.014E-621
1.681E-467
1.321E-314
1.349E-161
7.603E-009
1.013E-280
4.211E-127
1.013E-245
精度
128倍精度
128倍精度
128倍精度
128倍精度
128倍精度
188倍精度
188倍精度
248倍精度
最大誤差のビット数  有効ビット数  条件数のビット数
となってい ます。
248倍精度演算ではN  1400が使用できるメモリ容量の
上限となっ ています。
超多倍長演算での精度
N=100でのヒルベルト行列計算
精度
308
368
428
488
488
条件数0.127E+152
ビット数501
有効ビット数
理論値
2.80E-2812
2.52E-3390
2.66E-3968
2.80E-4546 (注)
2.80E-4546
9841
11761
13681
15601
15601
誤差
実測
4.73E-2814
3.38E-3392
1.61E-3971
5.36E-4198
1.00E-4547
(注)逆数計算で4倍精度の初期値から7回反復した場合の値
308,368,428倍精度は7回反復,488倍精度は8回反復が
必要
ヒルベルト行列のN=400の実行時間一覧表
方式
浮動型
浮動型
浮動型
整数型
整数型
整数型
整数型
精度
64倍精度
96倍精度
128倍精度
68倍精度
128倍精度
188倍精度
248倍精度
x5570
540.603
1551.696
3467.115
179.059
652.896
1235.567
2125.342
単位: 秒
e5430
743.938
2165.918
4280.500
277.818
918.142
1945.560
3337.811
32倍精度変数をつなげる方式に比べ,
整数演算方式の性能が非常に良くなっています。
超多倍長精度での性能
N=100でのヒルベルト行列計算
精度
実行時間　(秒)
x5570
e5430
308
50.298
81.123
368
70.468 115.572
428
95.133 156.090
488
123.216 205.238
488
123.253 205.404
反復回数7回
反復回数の増加の影響はほとんどありません。
4.量子モンテカルロ法による物性スペクトル計算への適用例
量子モンテカルロ法による物性スペクトル計算を数値計算的に言えば
サイズ NのL個の行列 B1 , B2 ,......,B Lから L  1個の行列を定め
G m  ( Bm Bm 1 .......B1 B0 )[I  B L B L 1 ......B1 ] 1 (0  m  L, B0  I )
L  1個のグリーン関数 Tr (G m ) (Tr : 行列のトレース )の値を求める事
となります。
行列積B m B m 1 .......B1をそのまま計算すると桁落ちが大きいので、
行列積計算をQDR分解を用いて計算します。ＱＤＲ分解とは,
行列Vをグラムシュミットの直交化でV=QR’
(Q:直交行列,R’：上三角行列)
にして,R’の対角要素からなる行列ＤとR’=DRとなるR
（上三角行列で対角要素は１）からV=QDRと分解する事
を言います。
計算の方法は以下の様になります。
( S. R. White et al, PRB 40, 506 (1989) が提案）
B m B m 1 ......B1を
B1  Q1D1R1　'
'
''
B 2 Q1  Q 2 D 2 R 2 B 2 B1  Q 2 D 2 R 2 D1R1  Q 2 D 2 D1 R 2 R1  Q 2 D 2 R 2
''
:
:
''
B m Q m 1  Q m D m R m
''
B m B m 1 ......B1  QDR
で計算します。
また[I  B L B L 1 ........B1 ] 1 は
B L B L 1 ........B1  Q1 D1 R 1とするQDR分解で
D1  Q1 R 1
T
1
 Q 2 D 2 R 2 を求めI  B L B L 1 ........B1  Q1Q 2 D 2 R 2 R 1
とQDR分解が出来ます。このとき
[I  B L B L 1 ........B1 ] 1  R 1 D 1Q 1  R 1 D 1Q T と演算量も削減されます。
Dは対角行列なので通常はベクトル(1次元配列)に収められます。
今回扱ったケースは静的かつ規則的なパス
Si
m
 ( 1) m に対するものでパラメ ータt 0 , N, L,
, Uのうち, t 0  1, U  10, N  100, L  448と固定。
(ここで, t 0 は電子のホッピング , Uは電子間クーロン
斥力に関連するパラメ ータでNは粒子数を表わして
います。)
G m : グリーン関数
U
L 1

), E ( m ) | Tr ( G m ) | として
L
2
L
1
E (0)  E ( L)  1, E (
) 
で結果の検証を行いました。
2
P
P  e Q (Q 
 の値は10000 /  が絶対温度(K )での温度を表わしています。
以前は倍精度変数を複数個つなげる方式では
  28でオーバーフロー, 整数演算方式で
32倍精度演算で  250,32倍精度変数を 2つつなげた
64倍精度演算で  600まで正しく計算できていました。
今まで扱った量子モンテカルロ法による物性スペクトル計算
では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  1
(
))とすると,
L
2
1
結果の絶対値の最小値は ,  が小さいので
P
Dの最大数はP 2 となり, グラムシュミット法で
P  e Q (Q 
平方数の和の平方根を求めるため,
プログラムに現れる最 大数はP 4 となります。
今回は扱う が大きいので条件数は / Lの大きさに依存
する度合いが大きくなります。
これは, B L ......B1  QDRとし ( I  QDR ) 1 を求める過程で
D  Q T R 1のQDR分解の際オーバーフロ ーが発生する事に
現れます。
正確には, b 
a1  a 2
2
2
 ......  a n の計算の際
2
a i  a i の計算でオーバーフロ ーが発生します。
これを防ぐには, a max  絶対値最大値 (a i , i  1,2,...., n )
として, b | a max |  c1  c 2
2
2
 .....  c n (c i 
2
ai
)
a max
とします。これで更に 大きな で計算出来る様になりましたが,
 / Lの値も大きくなりますので、正しい答えので る
 の最大数における, D' ' の最大数Ｓ,
Se
 UL 
としてその影響を調査しました。
調査ではN  20, U  10, L  448と固定しています。
これは, 実行時間短縮とN  20で正しい結果が得られる場合は
N  100でも正しい結果が得られている事によります。
調査結果
条件　N  20, U  10, L  448
精度
68倍精度
128倍精度
188倍精度
248倍精度
308倍精度
β の最大数
732
1620
2590
3609
3776
最大値
2.518D+1287
7.778D+2441
6.198D+3593
6.580D+4745
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.00000000000000000000000000000
8.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.1879571429D+04
0.1888000000D+04
0.1896428571D+04
0.3776000000D+04
-0.1000000000D+01
-0.7381849773-890
-0.1520746834-893
-0.7381849773-890
-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=
2
QDR = 8.395584469058393694708873673477994E+4931
I+QDR= 3.111402625840683560090283433225110E+4770
0.0000000000D+00
0.1880069196D+04
0.1888500000D+04
0.1896930804D+04
0.3777000000D+04
-0.4999999997D+00
-0.4586860821-890
-0.7614617438-894
-0.2814546226-890
-0.5000000000D+00
パラメータ領域拡大結果
測定結果一覧
x5570
測定条件　N  100, U  10, L  448
と絶対温度Kとの関係
絶対温度(K )  10000 / 
β
精度
730
1600
2500
10000/3
3770
68倍精度
128倍精度
188倍精度
248倍精度
308倍精度
絶対値最小値
実測値
0.824D-384
0.170D-582
0.773D-728
0.292D-840
0.156D-893
理論値
0.269D-383
0.214D-582
0.462D-728
0.978D-841
0.397D-894
実行時間
(秒)
57642
185097
403155
673631
975735
この条件では,32倍精度  250(絶対温度40K )
から,68倍精度で13.7 K ,128倍精度で 6.26K
188倍精度で 4K ,248倍精度で3K ,308倍精度
で 2.65Kまでパラメータ領域が 拡大されました。
これ以上の演算精度を用いても改善はされま せんの
でこの値が限界と言えます。
5.おわりに
整数演算方式を用いる事により,量子モンテカルロ法による
物性スペクトル計算において既存の方式では357Kまでしか
計算できなかったが、2.65Kまで計算可能とパラメータ領域を
拡大できました。また、ソースステップ数も8倍精度演算で1000
ステップ,その他の精度演算で1500ステップと作成作業
及びデバッグ作業まで一人で扱える範囲となりました。
今後さらにパラメータ領域を広げたり複雑なケースを扱うには
指数部の拡張が必要であり、現在、修正量ができるだけ少なく、
効率的に作業が行える方式を検討中です。
なお、これらの情報の更新は３か月毎にHPCI戦略5の
高性能計算の扉 http://www.jicfus.jp/field5/jp/promotion
で行っています。
最後にプログラムを提供してくださいました高エネルギー加速器
研究機構物構研 岩野薫氏、計算環境を整えていただいた
高エネルギー加速器研究機構石川正准教授、湯浅富久子准教授
に感謝いたします。