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

今回はパラメータ―の設定と精度に関してまとめて記述しました.
ループ積分と呼ばれる数値積分計算では,質量0の光子や質量が非常に小さい事はわかっているが,その値は不明なニュートリノに対して赤外発散を防ぐため微小量を与えて計算しています.
この設定する微少量の値により,結果の精度及び反復に要する時間が大きく作用したり,誤った値を得る事があります.
ここでは典型的な3つのケースで説明します.
(1)Infra box(2)W-vertex(3)bsgamma planar

(1)Infra box
1 1 x 1 x  y

1I     2 dzdydx
0 0 0 DD  sxy  tz (1  x  y  z )  ( x  y)2  (1  x  y  z )(1  x  y)me z (1  x  y)m f2

解析近似解
s  0, t  0, s,t  me  2の場合
2I 1
 s ( t  m f )2ln(s2( t  m f ) 22) ln

2me m f

最大相対誤差ln(s2 ( t  m f 2 )( t  m f ) 22) m ln(e2me m f2)22

要求する精度の結果を得るのに必要な仮数部のビット数は以下の様にして推定できます。
s  0では inf radの最小値は,box( x, y, z )  (1,0,0), (0,1,0)で2となります.
x  1   , y  0, zは0      1とすると,d  2   2 me となります.いま要求される相対誤差を10  n2とします.すると   me２22

10  n  (1 )よりが求まりますが,22

1dx  02 .
2d1 0 d1

これを 0とします. 1

1D.Eではx  1  の場合かかる重みは,  ((ln ) 2   2 ) 2となります.



1

また重みをw( x)としますと,  0  に対し  w( x)dx 1

n

1 2
1
 (ln   )   0 .
2211011   0  x  1ではd   (1 )ですから. ((ln ) 2   2 ) 2   02となるが表現できれば, D.E.で要求される精度の結果が得られる事に2なります.
1これから   ((ln( ))   )me
２2222



110  n  (1 )を解けば良い事になります.
2式を整理すると ((ln( ))   ) 
222



2

２22me 10 n.

  10 30 , n  6, me  0.0005とすると, x  ln( )2 x  ln( x 2   2 )  137.4619584  0  1.93  10 32.有効ビット106ビットとすると2 106  1.23  10 32.
このことから桁落ち等を工夫すれば, DD形式４倍精度でぎりぎり要求される精度の結果が得られる事になります.

s  500 2 , t  150 2 , m f  150, me  0.0005,   10 30 の計算を行うとします.計算はs   s  iとして複素数計算として実行します.
結果の精度は解析近似式を複素数計算して得られたものと比較します.式から虚数部分の値はに依存しない事がわかります.
また一般に演算に必要な精度はに依存します.
演算に必要な仮数部のビット数は log 21 要求精度の10進桁数.

例えば要求精度を10進6桁とすると, log 2 10 30  6  106となリます.
実数部の計算には4倍精度演算が必要で, 虚数部の計算はに依存しないので倍精度演算でも可能です.またD   D  iで計算すると虚数部の値は非常に収束性の良い列がえられますが, の初期値の取り方により3つの異なる値が得られます.

反復の最初の値をとします.

  me 2

ex : 25 0.142647572794761143D - 09 0.137D - 15m e    2

ex : 15 0.371536892867592250D - 08 0.266D - 25

2  ex : 15 0.743073784900021472D - 08 0.841D - 192

正しい値は２  の場合ですが, 結局4倍精度演算が必要という事になってしまいます.
s   s  i (iとしたのは符号を合わせるため)としてを10 5 程度にしますと, 倍精度演算で十分な精度が得られます.
このように, 計算の特徴に合わせて計算を行う必要があります.

(2)W-vertex

D   xys  ( x  y )( x  y  1)m 2  ( x  y ) M w  (1  x  y )ma2M w : W  boson

80.3477 2 , m : electron

20.0005

ma : neurino
この様にM w  m, maの場合, s  0では
1  4M w / s  1 2
11dydx(ln()) とm, maの値に依存しない事
0 0 D2
s1  4M w / s  11 1 x2になります.
計算手法の検討のため, s  0の場合を考えます.

以降２ページにs<0,s=0の場合の実行結果を示し,その後,s=0の解析式の導入方法を示します.
s  0の解析式の検討で数学公式に出てくる, 式からは,ln(1ma2)の項がでて, 発散するように見えるが,うまく切り抜けられる事を示しています.ma  0として計算すると演算量も削減され,誤差も簡単に見積もれる事がわかります.

S<0の場合の結果
s= -250000.00000000000000000000000000000
me2= 0.26112099999999999999999999999999780E-006
mw2= 6455.7528952900000000000000000000004
size=1024
kai= 0.549449295429470167Q-04
ma2:10^{-n}= 12.0 result= 0.549449295437675428Q-04
ma2:10^{-n}= 14.0 result= 0.549449295437683377Q-04
ma2:10^{-n}= 16.0 result= 0.549449295437683468Q-04
ma2:10^{-n}= 18.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 20.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 22.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 24.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 26.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 28.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 30.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 32.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 34.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 36.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 38.0 result= 0.549449295437683469Q-04
ma2:10^{-n}= 40.0 result= 0.549449295437683469Q-04

S=0の場合の結果
me2= 0.26112099999999999999999999999999780E-006
mw2= 6455.7528952900000000000000000000004
size=1024
kai= 0.154900600478343235Q-03
ma2:10^{-n}= 12.0 result= 0.154900600478342386Q-03
ma2:10^{-n}= 14.0 result= 0.154900600478343226Q-03
ma2:10^{-n}= 16.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 18.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 20.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 22.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 24.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 26.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 28.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 30.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 32.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 34.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 36.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 38.0 result= 0.154900600478343235Q-03
ma2:10^{-n}= 40.0 result= 0.154900600478343235Q-03

s0D  ( x  y )( x  y  1)m 2  ( x  y ) M W  (1  x  y )ma2

1 1 x

2

1x
1
dydx
0 m 2 x 2  (M 2  m 2  m 2 ) x  m 2 dx
0 0 DaaWa  m 2 , b  M w  m 2  ma , c  ma
222(2ax  b  b 2  4ac ) 21b1x
2 ax 2  bx  c dx  2a ln(| ax  bx  c |)  2a b 2  4ac ln | ax 2  bx  c |2acbb2ac 1 2,b 2  4ac  b bbb 2  4ac b(1  2ac )2b1(b  c) 2c1abc1x))  (  2 ) ln(ln(dx I  2c)cb
a(a2ca2bc
bx
ax
0(b  c) 2c
1 abc))  2 ln( ln(( a  b  c )cbcab1 ab)c  0, I  ln(ba2MW11 ab)) ≒ 2 ln(解析近似解I ≒ ln(2bamMW  m22meMW  m22m a1(b  c) 2c)|ln(|≒|)
ln() ||誤差 | 2 ln(242222( a  b  c )cbma
MWM W ma
( M W  m 2  ma ) 2222m 2  ma なら
21 1 x211
MW1dxx1
0 0 D dydx  0 m 2 x 2  (M 2  m 2 ) x dx  0 m 2 x  (M 2  m 2 )  m 2 ln( M 2  m 2 )WWW(3) bsgammaplanarD  c( x1 ( p1  m1 )  x2 ( p 2  m2 )  x3 m3  x6 m6 )222222

 c1 ( x5 p 2  x4 p1  x4 x5 ( p1  p 2 ))222222

 c2 ( x2 p 2  x1 p1  x1 x2 ( p1  p 2 ))222222 2 x2 x3 x5 p 2  2 x1 x3 x4 p1  x3 ( x2 x4  x1 x5 )( p1  p 2 )2222x6  1  x1  x2  x3  x4  x5
c1  x1  x2  x3
c2  1  x1  x2
c  x3 (1  x1  x2  x3 )  ( x1  x2 )(1  x1  x2 )  x3 (c2  x3 )  (1  c2 )c2
p1  mb , p 2  ms , m1  m2  1.5, m3  80.3477, mb  4.7, ms  0.0942222m6  
1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x4

I  

0 00001
dx dx dx dx dx
2 5 4 3 2 1
D1-xの形の計算があるため桁落ちに注意が必要になります.

式を整理すると下記の様になります.
D  x6 ( x3 ( x1 p1  x 2 p 2 )  c2 )22

 ( x 4 p1  x5 p 2 )[ x3 ( x1  x 2 )  ( x1  x 2  x3 )( x 4  x5 )]22

 c( x1 m1  x 2 m2  x3 m3 )22

2D0  D  cx6 2  x6 ( x3 ( x1 p1  x 2 p 2 )22

 ( x 4 p1  x5 p 2 )[ x3 ( x1  x 2 )  ( x1  x 2  x3 )( x 4  x5 )]22

 c( x1 m1  x 2 m2  x3 m3 )とおくと,
222D0   x3 ( x 4  x5 )( x1 p1  x 2 p 2 )  ( x 4 p1  x5 p 2 )2222( x3 ( x1  x 2 )  ( x1  x 2  x3 )( x 4  x5 )) c( x1 m1  x 2 m2 )22

 x3 (cm3  (1  x1  x 2  x3 )( x1 p1  x 2 p 2 ))222cm3  ( x1  x 2 )(1  x1  x 2 )m3  ( x1  x 2 )(1  x1  x 2  x3 )m3
222cm3  (1  x1  x 2  x3 )( x1 p1  x 2 p 2 )  (1  x1  x 2  x3 )222(( x1  x 2 )m3  x1 p1  x 2 p 2 )222

 (1  x1  x 2  x3 )( x1 (m3  p1 )  x 2 (m3  p 2 ))  0
2222

以上からD  D0  0.
よってD  0ならD0  0.D0  0ならc  0より, D  0式の形から, 定義式どうり計算するより, 式を変形した形で計算した方が, 精度, 演算量からみて, 良い事がわかります.

性能モニターでの演算量

bsgamma planar

演算量

originn=128n=256n=512n=1024

14844GFLOP475008GFLOP
15200TFLOP
486PFLOP

tune
n=128n=256n=512n=10246489GFLOP
207648GFLOP
6645TFLOP
213PFLOP

15200/6645=2.29
演算量は式の変形のみで２倍以上の差があります.

計算する前にbsgamma

planar

でM Ｗがない場合で, 発散は起きない事.誤差は２になる事を確認しました.
bsgamma

でM Ｗがない場合の式はplanar

以下の様に変わります.
D   x(1  x  y ) p1  y (1  x  y ) p 2  xm1  ym2  (1  x  y )22222(m1  m2  m 2 )221 1 x

1dydxとなります.
D0 0I  22p1p2  0とすると, I (SP()
SP())
2222mm( p1  p 2 )1

この式からI  0.02593114637754185125591974634710566
これを4倍精度DQでD-  D  iの算法で計算した結果は下記のとうりで, 発散は起きていません.
eps = 0.13822538D - 05
0.2593108797935867D - 01 0.240D - 25000ex 0.2593114637772012D - 01 0.841D - 07 0.25931D - 01
ex 0.2593114637754155D - 01 0.361D - 12 0.25931D - 01
ex 0.2593114637754185D - 01 0.642D - 15 0.25931D - 01
ex 0.2593114637754185D - 01 0.221D - 19 0.25931D - 01
ex 0.2593114637754185D - 01 0.145D - 23 0.25931D - 01
ex 0.2593114637754185D - 01 0.128D - 24 0.25931D - 01
ex 0.2593114637754185D - 01 0.128D - 25 0.25931D - 01
ex : 15 0.2593114637754185D - 01 0.501D - 27に値をいれると, 誤差が2になる事が以下の結果でわかりました.

  10 5

ex : 15 0.2593114641144115D - 01 0.116D - 11  10 7

ex : 15 0.2593114637754597D - 01 0.493D - 16

  10 -8

ex : 15 0.2593114637754185D - 01 0.177D - 16

  10 -10

ex : 15 0.2593114637754185D - 01 0.255D - 20

  10 -15

ex : 15 0.2593114637754185D - 01 0.495D - 27

  10 -30

ex : 15 0.2593114637754185D - 01 0.501D - 27

(注)  10 -15 ,10 30 で誤差に変化がないのは, プログラムで相対誤差の許容値を10  24 に設定している事によります.

このあと,BG/Qで必要な分点数を定め,計算はSR16000/M1で実行したものを採用しています.
これは,必要な分点数を求めるには演算量が多くなる
事と,この計算は精度に鋭敏な事によります.
分点数N=1024 BG/Q 512ノード 28537秒
分点数N=512 SR16000/M1 8ノード 8486秒
演算量は分点数の５乗に比例.

積分計算は、二重指数関数型積分で行い, 変数変換区間は[10 150 ,1  10 150 ]で実施しています.また, 一般にこの種の積分は1  xの形の式がでますので, 桁落ちを考慮してDD形式の4倍精度を使用する必要があります.

また区間幅が結果の精度に影響しますので,分点数を変化させ,必要な分点数をBG/Qでもとめました。(λ=0)
結果
result=result=result=result=result=result=result=result=result=.331405208498465261D-01.488169643516346941D-01.460458451277416842D-01.497962422755366272D-01.544217811073604471D-01.565118842794349136D-01.558519811521311674D-01.566275466678280369D-01.565795177564476502D-01

分点数
128
192
256
320
384
512
640
768
1024

SR16000/M1 8ノードでいくつかのλの値とλ=0の場合の値は以下の様になっています.

value
0.0600000000000000
0.0500000000000000
0.0400000000000000
0.0300000000000000

value

0.0200000000000000
0.0100000000000000

n   log10 (2 )0.0000000000000000
40 60 80 100 120 140 160nResult=0.0566344454366303263