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

目次
1. はじめに
2. 反復解法
2.1 対称行列
2.2 非対称行列
2.2.1 ｂｃｇ法
2.2.2 ｃｇｓ法
3. 乱数ルーチン
3.1 手法概略
3.2 並列化コーデイング例
3.2.1 乗算型合同式法
3.2.2 混合（線形)型合同式法
4. ループ積分(数値積分)における特異点処理

１．はじめに数学と計算機で実行した結果は必ずしも一致しません。これは,１．計算機で使用される変数は連続でなく離散的であること。
２．浮動小数点演算では有限のビット数しかもっていないので演算毎に誤差が発生。
３．また表現可能な数値範囲が有限であること。
４．極限操作が有限回の操作で行うこと。などが挙げられます。
このため、必要とされる精度を持った結果を得るには、それに必要なビット数が問題毎にあります。
一般に使用されている演算精度は倍精度ですが、これではビット数が不足すると多倍長計算が必要となります。
もう一つの問題に極限操作があり,加速法の使用により,収束範囲が数学より広くなる場合があり,結果の妥当性のチェックが要求される場合が発生します。これは必要なビット数との相互作用がある場合があるため、ここでも多倍長計算が必要になってきます。

2.反復解法
2.1 対称行列
3次元ポアソン方程式を差分法で解く際に現れる係数行列をA,Ax  bでx T  (1,1,....,1,1)になる様にbを設定.x T の初期値をx0  (0,0,...,0,0)としてサイズ129  129  129
T

収束判定値は残差10 12とする。
３次元反復法SR16000/M1

129*129*129 収束判定値 1.0e-12

解法
bcg
cgs
scg
bicgs
cgsmilgpbicgcgsilucgsiluq

倍精度

実行時間(秒)
反復回数 single32core
64smp
618
24.2743.571
2.8641416.649
2.721
2.232
618
16.2551.975
1.799
16.339
2.393
2.382
反復回数
402395420
26.12518.98519.739
反復回数
27729230824.897
3.094
2.707
反復回数
434381414
精度single32core64smp
倍精度
944.418 644.799 635.057
反復回数
10000
10000
10000
4倍精度
307.614 159.012 158.055
反復回数
275
275
275

(注) 10000 は収束しないのでこの回数で打ち切ったもの

この様に比較的収束し易い問題でも４倍精度が必要になる事があります。

2.2 非対称行列

反復回数の調査は,SR16000/M1とx5570を使用して行っています。

2.2.1 ｂｃｇ 法
演算精度と反復回数は以下の様になっています。これより,４倍精度演算が良い事がわかります。
bcg 法

反復回数
3500
30002911

2500

2000
反復回数

1500
1000
500

1075

346

173

163

060 68 76 81 86 94 102113177241

仮数部のﾋﾞｯﾄ数

2.2.2 cgs 法
演算精度と反復回数は以下の様になっています。
これより,６倍精度演算が良い事がわかります。
cgs 法

反復回数
3500
3000

2911

2500
2000
反復回数

1500
1000
587

500

364
151

128

0
97

99

102 106 145 177 212

仮数部のﾋﾞｯﾄ数

3.乱数ルーチン
3.1 手法概略
疑似乱数に関する概略説明(1)混合型合同式法(mixedcongruencemethod )x n 1  x n   (mod M )（ア）とMは互いに素（イ）PがMの素因数ならば  1 (mod M )（ウ） 4がMの素因数ならば  1 (mod 4)ならば周期M .
これは, M  2 s ならば, が奇数で  1 (mod 4)
すなわち,   2 a  1, a  2はすべて周期2 s .
系列相関係数x n , x n 1の間の相関係数 ≒(1 6
(1  )  K
 MMK )M
M(2)乗算型合同式法
x n 1  x n

(mod M )M  2 35 ,   2 7  1
M  2 47 ,   2 9  1
M  2 s ; x0が奇数,   5(mod 8),または  3 (mod 8)ならば最大周期2 s  2.

  5 2 k 1 ,   7 4 k 1 ,   1313

３．２ 並列化コーデイング例
3.2.1 乗算型合同式法
x n 1  x n (mod M )x0  iseed (乱数ルーチンの種)
一般にM  2 s ; x0が奇数,   5(mod 8),または  3 (mod 8)ならば最大周期2 s 2となる.事と混合型（線形)合同式法にみられる系列相関係数x n , x n 1が大きくなる事のない  1313 , M  2 53 (2 59 の方が良いと考えられますが, 拡張倍精度が必要となるためにM  2 53にしています.)オリジナルソースをinline 展開したソース
fave=(double) (ib)*r ;for (i=1;i<68719476736;i++) {ic=ia*ib ;ib=ic&mask ;g05faf=(double) ( ib) *r ;fave=fave+g05faf ;}fave=fave*t ;での結果 fave= 0.50000054390657433512
に対し 修正した結果
1cpu
fave= 0.50000054390657433512
8smp fave= 0.50000054390660353398となり,1ｃｐｕでは結果は最後まで一致しています.

プログラムでは, x n 1  x n (mod M )を変形しx n  n x 0 (mod M )として先にn (mod M )を計算してテ  ブルにします.
long int i, j, k ;ilist[0]  1 ;ib  1 ;for(i  1; i  n; i   ) {ic  ia * ib ;ib  ic & mask ;ilist[i]  ib ;}ic  ia * ib ;ia  ic & mask ;jlist[0]  1 ;ib  1 ;for(i  1; i  n; i   ) {ic  ia * ib ;ib  ic & mask ;jlist[i]  ib ;}ic  ia * ib ;ia  ic & mask ;klist[0]  1 ;ib  1 ;for(i  1; i  n; i   ) {ic  ia * ib ;ib  ic & mask ;klist[i]  ib ;}実行部分のソースは以下のようになります。
#pragma omp parallel private(i,j,k,ia,ib,ic,g05faf)
reduction(+:fave) ¥ num_threads(8)#pragma omp forfor (l=0;l<68719476736;l++) {k=l>>24 ;j=(l-(long int)(k<<24))>>12 ;i=l-(long int)(k<<24)-(long int)(j<<12) ;ib=klist[k] ;ia=jlist[j] ;ic=(ia*ib)&mask ;ia=ilist[i];ib=(ia*ic)&mask ;ia=(ib*iseed)&mask ;g05faf=(double) ( ia) *r ;fave=fave+g05faf ;}fave=fave*t ;テーブル作成には、iseed(種)は全く影響しません。並列実行時は参照のみとなります。
N=2^36での実行時間(秒)は,オリジナル(ｉｎｌｉｎｅ化ソース)
修正ソース
修正ソース(smp)

e5430
235
415
72(8smp)コンパイルオプション icc -Os -lmt2k
248
589
44(16smp)

同じ問題をＭＰＩ（フラットＭＰＩ）で実行した結果は以下の様になっています.

乗算合同法
T2K
BG/Q ３２ノード
mpi数 実行時間(秒） mpi数
実行時間（秒)
16 40.302161
1 2576.929578
32 20.273622
512
5.090237
64 10.1492
1024
3.278346
128 5.073721
2048
2.882218
256 2.550183 7.297119
並列化効果がはっきりでています。

拡張倍精度乱数列rn  arn 1(mod 2 59 )a  1313 , iseed (種)  511

N  2 36 での平均値
倍精度

0.50000207337461433088

拡張倍精度
４倍精度
実行時間(秒)0.500002073371655848301
0.500002073371655848463279703253193(e54301cpu )

倍精度

234.624821

拡張倍精度

217.760483

４倍精度

3240.429380

条件
倍精度 icc -Os -lm
拡張倍精度 icc -Os –fp-model extended -lm
4倍精度 ifort -Os

精度、性能面では拡張倍精度演算が倍精度演算より良い事を示している。

3.2.2 混合(線形)型合同式法
rn  arn 1  c (mod M )nrn  a r0  (  a m 1 )c (mod M )nm 0n

乗算合同法のテーブルに加えて,  a m 1のテーブルを作成する.
m 0l  i  j  4096  k  4096 2(0  i  4095,0  j  4095,0  k  4095, l  0)
乗算合同法のテーブルをilist , jlist , klistとして, 追加するテーブル
を, ilist1, jlist1, klist1とする.
iilist1(i )   a m(0  i  4095)m 0jjlist1( j )   (a 4096 ) m(0  j  4095)

m 0kklist1(k )   (a 40964096 ) m(0  k  4095)

m 04095A   a  ilist1(4095), B 
m40964096 1
m

m 0l

am 1 a  A  jlist1(4096)

m 0 B  klist1(k  1)  klist (k )  [ A  jlist1( j  1)  jlist ( j )  ilist1(i  1)]

m 1(mod M )k  0のときは, klist1(k  1)  0, j  0のときはjlist1( j  1)  0,i  0のときは, ilist1(i  1)  0とする.

オリジナルソースをinline化したものは以下の様になります。
for (i=1;i<68719476736;i++) {ie=ia*ib+ic ;
ib=ie&mask ;g05faf=(double) ( ib) *r ;fave=fave+g05faf ;}fave=fave*t ;N=2^36での実行時間(秒)は,オリジナル(ｉｎｌｉｎｅ化ソース)
修正ソース
修正ソース(smp)

e5430
262
708
120(8smp)

t2k2561145
82(16smp)となっています.
コンパイルオプションはicc -Os –lm です.
具体的なテーブル作成部分と実行部分のソースを以降に記述しました。

追加するテーブル部分のコーデイング
ilist1[0]=1 ;for(i=1;i<n;i++) {ie=ia*ilist1[i-1]+1 ;ib=ie&mask ;ilist1[i]=ib ;}
ih=jlist[1] ;jlist1[0]=1 ;for(i=1;i<n;i++) {ie=ih*jlist1[i-1]+1 ;ib=ie&mask ;jlist1[i]=ib ;}
ih=klist[1] ;klist1[0]=1 ;for(i=1;i<n;i++) {ie=ih*klist1[i-1]+1 ;ib=ie&mask ;klist1[i]=ib ;}結果にかんしては、オリジナルの
fave= 0.50000176689742203973
に対し,シングルジョブでは
fave= 0.50000176689742203973
8ｓｍｐでは以下の様になっています.
fave= 0.50000176689738828895

実行部分のソースは以下の様になります。

#pragma omp parallel private(i,j,k,ib,ih,ie,ik1,ik2,ij1,ij2,ii,g05faf)
reduction(+:fave) num_threads(8)#pragma omp forfor (l=1;l<68719476736;l++) {k=l>>24 ;j=(l-(long int)(k<<24))>>12 ;i=l-(long int)(k<<24)-(long int)(j<<12) ;ib=klist[k] ;ih=jlist[j] ;ie=(ih*ib)&mask ;ih=ilist[i];ib=(ih*ie)&mask ;次ぎに続きます。

if(k==0) {ik1=0 ;
ik2=1 ;}
else {
ik1=klist1[k-1] ;ik2=klist[k] ;}if(j==0) {ij1=0 ;
ij2=1 ;}
else {
ij1=jlist1[j-1] ;ij2=jlist[j] ;}if(i==0) {ii=0 ;}
else {
ii=ilist1[i-1] ;}ie=b*ik1+ik2*(a*ij1+ij2*ii) ;ih=(ib*iseed+ie*ic)&mask ;g05faf=(double) ( ih) *r ;fave=fave+g05faf ;}fave=fave*t ;ＭＰＩでの値は以下の様になっています。

混合法
T2K
BG/Q 32ノード
mpi数 実行時間（秒）mpi数
実行時間（秒)
16 75.631308
1 6335.624106
32 38.844884
512 12.06219
64 21.414003
1024 6.877529
128 10.84745
2048 5.518025
256 7.297119
乗算合同法と同じく並列化効果がはっきりでています。

4.ループ積分(数値積分)における特異点処理
ループ積分では複数の数値積分を行い, その和を最終値とする処理を行っています。
nS   Ai
i 1

Ai : 個々の数値積分の値, S : 最終値

ここで, Sは有限値となりますが, 個々のAiでは発散する場合があります。
一般に発散する被積分関数Dに対しては、微小量iを加え,    0としたときの値で発散する積分の有限部分を取り出す操作を行います.ただし計算機では数学的な   0の操作は出来ないので,

 n   0となる点列を取り, D  i nの値の極限値を有限回の操作で求める事になります。

数学的に求まる値と計算機で求まる値
0  a  1,   0110 x  a dx :計算機で求まる値  主値積分の値

110 ( x  a) 2 dx :計算機で求まる値  Hadamard積分の値

説明は次ページ以降に記述します。

11 x  a dx
0a lim[   00111 adx  dx]  lim (ln(1  a )  ln a )  ln()  0xaxaaa 1xa1(1  a ) 2   21 alim[dx]  lim[ ln()]  ln()22
  0  ( x  a ) 2   2
  0 2aa 0111

 ( x  a) dx
20a 

111112lim[ dx  dx]  lim[ ( ) ]22
  0  0a 1 a
0 ( x  a)a  ( x  a )有限部分  (11)a 1 a( x  a) 2   2
 (1  a )a
11dx]  lim[ 2]  ( )
22 2222  00a1
a((xa))(1a)a
01lim[ S n  (112)
a 1 an( 0   1  .....   n   0)

エイトケン加速法、及びその一般形の  算法では(11)に収束する.
a 1 aS n  a  b n(|  | 1)

エイトケン加速法

S n  S n  S n 1  bn 1 (  1)
2 S n  S n  S n 1  bn  2 (  1) 2( S n ) 2b 2 2 n  2 (  1) 2nTn  S n  a  b a2n22 Snb (  1)

特異点が端点の場合の例にちなむ)

H ( HadamardbH  f ( x)dx

区間[a, b]での発散積分の有限部分.

a

積分のCauchyの主値は発散積分の有限部分の一例b  算法は端点特異点の場合H  f ( x)dxを求めている事になる.

af ( x) : [a   , b]で有界, 積分可能で, aの近傍でLaurent展開が可能。
mc0crf ( x)  h( x )x  a r 1 ( x  a ) r

[ r  1](h( x)は[a, b]で積分可能)bb

H  f ( x)dx  lim[  f ( x)dx  c0 loga  0a b1

mc r 1 r 1( ) ]1r 1 r

mcr
1 r 1()r 1  r  1 b  a  h( x)dx  c0 log(b  a )  a1 x
1 11dx

()
3
1 2
20 x

1H

h( x)  0, c0  0, c1  1, c 2  1, 1  2,  2  3

1 x
1 1 11dx  322 2 nn
n x
1an  ( n 1)2nの  算法の結果は

ex - 0.364048765700024415D + 04 0.104D + 10 - 0.36405D + 04
ex 0.500000000000000000D + 00 0.637D + 04 0.50000D + 00
ex 0.500000000000000000D + 00 0.147D - 25 0.50000D + 00
ex 0.500000000000000000D + 00 0.235D - 26 0.50000D + 00
ex 0.500000000000000000D + 00 0.861D - 28 0.50000D + 00
ex 0.500000000000000000D + 00 0.319D - 29 0.50000D + 00
ex 0.500000000000000000D + 00 0.281D - 29 0.50000D + 00
ie = 150.500000000000000000D + 00

0.476667661235840874D - 31

ε-算法の結果はHadamard積分の値に一致しています。

より一般的な場合の例
inf ra box s  0,   0, t  150 2 ,m f  150, me  0.0005の計算結果
1 x1dx
32x01Hとした場合の解析解  44.4444444507677456465204143093617

の初期値は1.2 -90 ,  n   n 1 / 1.2 で15回反復計算.
DEeps = 0.58230485Q - 08
result = - 0.602472064922940173Q + 03
ex 0.443856621401209096Q + 02 0.755Q + 03 0.44386Q + 02
ex 0.444445514441291720Q + 02 0.122Q + 00 0.44445Q + 02
ex 0.444444440470920105Q + 02 0.461Q - 03 0.44444Q + 02
ex 0.444444444236695617Q + 02 0.348Q - 05 0.44444Q + 02
ex 0.444444446858620597Q + 02 0.352Q - 06 0.44444Q + 02
ex 0.444444444677995797Q + 02 0.451Q - 06 0.44444Q + 02
ex 0.444444444079332656Q + 02 0.901Q - 06 0.44444Q + 02
ex : 15 0.4444444468586206Q + 02 0.225Q - 06

ループ積分においては被積分関数Dに対し,微小量iεを加え,ε-算法により積分値を求める
方法は非常に有用な手法となっています。