最近ぼちぼち色々な計算機で FP16 が使えるようになってきたので、
それで多体問題、少なくとも重力多体問題ができないかという雑談を
何人かでしていた時にでてきた方法について。
仮数長は 10ビットもあるので、GRAPE-1, 1A, 3, 5, 8 あたりに比べて長く、
相互作用計算にはいってからの精度には特に問題はありませんが、問題は
-
座標表現、積算には精度不足
-
途中で指数がオーバーフローないしアンダーフローする
ということです。座標表現は、実効的に固定小数点にして、その上位と下位を
別々のワードに格納し、別々に引き算してから加算することで精度を確保する
ことはできなくはありません。但し、2語使っても20ビットなので、3語くらい
使う必要がありそうです。
とはいえ、自分は原点近くにいて、遠くのほうは相対精度が落ちてよいとすれば、遠くは「桁を合わせない」という方法でなんかできる気がします。
積算は複数語のアキュムレータを使う方法が古くから知られています。
問題は途中でのオーバーフローです。距離のダイナミックレンジ自体は30ビッ
ト程度はほしいので、途中で必要な距離の3乗の逆数を計算すると90ビットに
なり、 FP16 で持つ5ビットの指数では表現できなくなります。
ここで、ツリー法を使い、質量の空間分布が2ないし3次元でほぼ一様な場合だ
けをまず考えると、相互作用を計算する超粒子の質量は距離の 2 ないし3乗に
比例します。まず、色々上手くいきそうな 2次元の場合を考えると、
距離と質量の平方根が比例するので、距離を質量の平方根でわったものは
ほぼ一定になります。従って、普通なら(dx 等は3次元ベクトルで内積があるとして)
dx = xj-xi
r2= dx*dx
r=sqrt(r2)
mr3inv=mj/(r*r*r)
fi += mr3inv*dx
と書くところを
q=1/sqrt(mj)
dxhq = q*(xjh-xih)
dxlq = q*(xjh-xil)
dxq= dxhq+dxlq
r2q2= dxq*dxq
rq=sqrt(r2q2)
mr3invqinv=1/(rq*rq*rq)
fi += mr3invqinv*dxq
としてやると、 dxhq や dxlq はダイナミックレンジが小さく、基本的に1の
オーダーになっているので、指数のオーバーフローやアンダーフローが
ほぼ起きなくなることが期待できます。最初の、質量の逆数の平方根をとると
ころは、ツリー構築を作った時に1度だけすればいいので、計算量の増加は無視できます。
3次元の場合には、一様分布の場合遠くの超粒子からの力が近くの粒子からの
力より大きくなりますが、これは距離に比例する程度なのでやはり上のような
修正を行うことでダイナミックレンジを大幅に小さくできます。
3次元で TreePM とか使うなら距離にはカットオフがあるのでこれで問題ない
でしょう。開放境界の3次元問題の時には、結局、トータルの力が決まるあた
りの特徴的な距離スケールがあるので、そこでオーバーフローが起きないよう
にスケールを決めることができるはずです。