日本物理学会2013年秋季大会（於 高知大学朝倉campus）

講演21aSB-6 (2013年9月21日)

高スピン間の回転行列の数値評価における著しい桁落の回避方法

田嶋直樹（福井大工）
1.回転演算子の角運動量固有状態基底での表現行列（D関数）をWignerの公式で数値的に求めると、角運動量jが大きいとき著しい桁落ちが起きる。
「j=1/2ないし1との角運動量合成に関連する漸化式」で求める場合も同程度の桁落ちが起きる。
2. D関数をフーリエ（有限）級数で表せば非常に大きなjでも桁落ちはほとんど起きないことを示す。
3.本講演では、フーリエ 展開係数は数式処理で求めてデータファイルに書き出しておき、数値計算プログラムはそれを読み込んで使う。数値計算プ
ログラム内で係数を計算できる桁落ちの小さい公式を考案すべく研究中。

回転行列 (D 関数) とは

Euler 角 ϕ, θ, ψ で指定された回転を行う演算子
R̂(ϕ, θ, ψ) = e−iϕJz e−iθJy e−iψJzˆˆˆの「角運動量の固有状態で構成された基底」での行列要素である。
jDmk= ⟨jm|R̂(ϕ, θ, ψ)|jk⟩ = e−imϕe−ikψ djmk (θ),

djmk (θ) = ⟨jm|e−iθJy |jk⟩ˆ

基底の位相の規約により Jˆy の行列要素は純虚数なので、dimk (θ) は実数値をとる。
j• Dmkの表式は Wigner による定義であり、Bohr-Mottelson はその複素共役で定義する。

しかし、djmk (θ) に関しては両者の定義で違いは生じない。

• R̂ はユニタリ演算子であるから回転行列は規格化された状態間の重なりである。
∴ |dimk (θ)| ≤ 1

Wigner の公式 [1]：djmk (θ) =∑n(−1)ntn(j, m, k; θ),n : integer in [max(0, k − m, ), min(j − m, j + k)],√()()(j + m)!(j − m)!(j + k)!(j − k)!
θ 2j+k−m−2nθ m−k+2ntn(j, m, k; θ) =· cos 2· − sin 2,(j − m − n)!(j + k − n)!(n + m − k)!n!

この公式は j が大きい場合に著しい桁落ちが起きる。
例えば、j が整数、θ = π2 , m = k = 0 の場合、
(2π)1j!

 ,tn j, 0, 0;= j
22 (j − n)!n!となり、j が偶数なら n = 2j で最大値
2vuu 21  j! π)u ∼ ttj/2 j, 0, 0;= j
2j22 (j/2)!
πj(をとる。
• この巨大な値をとりうる項 tn の差し引きとして絶対値が 1 以下である djmk (θ) を求めるようとしても、j ∼ 54 (114) で倍精度（四倍精度）浮動小数点実数の精度は完全に失われてしまう。
• Wigner の公式の基となったような (角運動量合成 j + 21 = j ⊗ 21 , j + 1 = j ⊗ 1 等にかかわる）漸化式を直接用いても桁落ち問題に大幅な改善はないことが文献 [2] からわかる。
[1] e.g., see M.E. Rose, Elementary Theory of Angular Momentum, 1957.
[2] “Fast and accurate determination of the Wigner rotation matrices in the fast multipole method”,
H. Dachsel, J. Chem. Phys., 124, 144115 (2006)

( )

π
d60
0,0 2

=

012345678910111213141516171819202122232425262728
2930313233

0.000000000000000-0.000000000000003
0.000000000002717-0.000000001015688
0.000000206248120-0.000025871764215
0.002173946854217-0.129372020957098
5.678218857320146-189.55436778016883*
4930.309105962191***-101865.89061905353****1698472.2456690802*****-23155503.278234087******
260971973.171526*********-2454296423.2486625********
19413868191.71305**********-130052764080.12619**********742183829580.7201***********-3626626801607.7295***********15240899133756.482************-55295779170091.55*************
173770413466341.4**************-474337385719086.94*************1127374793488593.8**************-2337724371777948.***************4236260880810631.***************-6717582411820424.***************9330927610296482.***************-11361319706234956.***************12131364708546436.***************-11361319706234956.***************9330927610296482.***************-6717582411820424.***************

344236260880810631.***************
35-2337724371777948.***************
361127374793488593.8**************
37-474337385719086.94*************
38173770413466341.4**************
39-55295779170091.55*************
4015240899133756.482************
41-3626626801607.7295***********
42742183829580.7201***********
43-130052764080.12619**********
4419413868191.71305**********
45-2454296423.2486625********
46260971973.171526*********
47-23155503.278234087******
481698472.2456690802*****
49-101865.89061905353****504930.309105962191***
51-189.55436778016883*
52
5.678218857320146
53-0.129372020957098
540.002173946854217
55-0.000025871764215
56
0.000000206248120
57-0.000000001015688
58
0.000000000002717
59-0.000000000000003
60
0.000000000000000
total = d^{60}_{0,0}(pi/2)= 0.102578173008570

フーリエ（有限）級数による djkm の表式
Wigner の公式に現れる cosn 2θ sinm 2θ の形をした各項は、x = 2θ として、三角関数の恒等式

sin2 x = − 12 cos 2x + 12sin3 x = − 34 sin 3x + 34 sin x······
cos2 x = 12 cos 2x + 12cos3 x = 14 cos 3x + 34 cos x······
sin µx cos νx = 12 sin(µ + ν)x + 12 sin(µ − ν)xkθを利用すれば {sin kθ,cos
22 |0 ≤ k ≤ n + m} の線形結合の形に表せることがわかる。
詳細に考察すると下記の一般式を得る。

 even jjmk cosdmk (θ) = cn nθ
form−k: odd  sin n∑0evenwhere n = j, j − 1, · · · 
for2j:
 1 
 odd 
2cjmk=n∫ 4π


1cosevenjd(θ)
nθdθform−k:
 sin 
 odd 2νπ 0 mk 






2 · · · for n = 0 ∩ m − k:even jwhere ν = andWigner
formula→dmk
 1 · · · other wise

具体例

d7
2
1 ,− 1227 θ sin2
(θ) = 144 144(= −d72
1 ,− 322( )( )35 sin(θ) = 48

√15



7t2)

−

cos2 2θ

− 5 sin( )cosθ
2(+125t2( )( )

sin5 2θ)(+ 15 sin64

( )sin6 2θ

48

cos4 2θ3t2( )−)t2sin
8−( )

cos6 2θ

sin 2θ 
36( )− 9 sin 2θ
( )sin4 2θ

cos3 2θ( )( )3+( )

sin2 2θ 



24( )
( )
( )
( )
√√
√√
7t
5t
3t
7 15 cos 2 − 3 15 cos 2 − 15 cos 2 − 3 15 cos 2θ
= −64
( )

9 cosθ
212( )

cos5 2θ( )sin7 2θ

cos

3( )

t2( )

sin5 2θ( )cos5 2θ( )

sin3 2θ −2404872√√√
37 2 sin (4 t) − 7 2 sin (3 t) − 2 2 sin (2 t) − 3 2 sin t + sin(0 · θ)=32d41,−2(θ) = 45 · 2 2 −d42,0(θ) = 288
= −

7

√
√



( )

2 θ
 cos
2


10 ( )

sin6 2θ

96+( )−

cos4 2θ

10 cos (4 θ) +0 · cos(3θ) − 4√( )sin4 2θ

36( )+cos6 2θ( )

sin2 2θ 9610 cos (2 θ) +0 · cos(2θ) − 364√

10 cos(0θ)

本研究では現在までに• フーリエ級数展開の係数 cjmkを数式処理で求めた。
n• 数式処理ソフト Maxima を使用した。
これまでの約 20 年間は、もっぱら京大基研の Mathematica にお世話になりました。

• 人間の側の手間を惜しまなければ、数式は不要で、任意長の分数計算だけで計算できる。
• 求まった係数は倍精度実数に精度を落としてデータファイルに保存した。
• jmax 以下の全て（整数＋半整数）の j の djmk を表すための係数の総個数は約

j4
2.

8byte の倍精度実数なら jmax = 50 で 25MB, jmax = 100 で 400MB, jmax = 200 で 6.4GB• その係数データを読み込んで djmk (θ) の値を倍精度実数で数値的に求めるプログラムを作成した。
• 以下では、このプログラム および「Wigner 公式通りに倍精度実数で数値的に総和を計算して求めるプログラム」の誤差を「数式処理で求めた正確な値」と比較して求め、グラフにして示す。
各 j 値に |m|, |k| ≤ j をみたす (2j + 1)2 個の (m, k) 対があるが、対称性

djmk (θ) = (−1)m−k djkm(θ) = (−1)m−k dj−m,−k (θ) = dj−k,−m(θ)
を考慮して 約 14 の領域 {(m, k)|0 ≤ m ≤ j, |k| ≤ m} だけを示す。

cos nθ, sin nθ の計算方法
三角関数の加法定理により計算すると計算量が少なく精度も高い:

n が整数のときnn− 1 2cos nθ, − sin nθ   cos θ, − sin θ  = sin nθ, cos nθsin θ, cos θn が半整数 (=半奇数) のとき

cos nθ, − sin nθ   cos θ, − sin θ = sin θ, cos θsin nθ, cos nθ

cos 2θ ,sin 2θ ,− sin 2θ ,cos 2θ

cos θ = cos2 2θ − sin2 2θ sin θ = 2 sin θ cos θ22

一方、nθ の値を計算してからそれを関数 sin, cos に渡す計算方法では、sin,cos を n 回ずつ呼ぶため計算時間が長くなることに加えて、|nθ| > π4 ならば、 π2 の整数倍を差し引いてそれを [− π4 , π4 ] に還元する際に桁落ちが生じる。
加法定理の適用を n 回繰り返して生じる数値誤差は、この桁落ちより小さい。

まとめ1.回転演算子の角運動量固有状態基底での表現行列（D関数）をWignerの公式で数値的に求めると、角運動量jが大きいとき著しい桁落ちが起きる。
j=1/2ないし1との角運動量合成に関連する漸化式で求める場合も桁落ちに大幅な改善は見られない。
2. D関数をフーリエ（有限）級数で表せば非常に大きなjでも桁落ちはほとんど起きないことを示した。
3.本講演では、フーリエ 展開係数は数式処理で求めてデータファイルに書き出しておき、数値計算プログラムはそれを読み込んで使った。現在、数値計算プログラム内で係数を計算できる桁落ちの小さい公式を考案すべく研究を継続中である。
4. 関数値が非常に小さい場合に限り、Wigner公式のほうが誤差が小くなることがある。あらゆる用途を想定して、最終的には引数によってはWigner公式に切り換えるようプログラムを改良すべきである。