[[PukiWiki/1.4/Manual/Plugin]]

*TEST [#s0532dca]
#math( y = \sqrt{x} )
#math( z = \frac{1}{y} )
#math( y = ax + b )

*逆数平方根 [#l0a5570b]
以下、基本的に
[[逆数と平方根を求める高次収束アルゴリズム:http://www.finetune.jp/~lyuka/technote/fract/sqrt.html]]
をベースに数式を整形したものです。
また、SSEでの実装については
[[advanced optimization SSE:http://homepage1.nifty.com/herumi/adv/adv50.html]]
も参考になるでしょう。

除算と平方根は一般的に加減乗算よりも計算コストが一桁高く、
数値計算においてもこの部分を高速化することが重要になることはよくあります。
よく知られている方法として、逆数（&math(1/x);）ないし逆数平方根（&math(1/\sqrt{x});）
の適当な初期値推定から、
加減乗算のみを用いて精度を改善するというものがあります。

最近はハードウェアで高速にこの初期値推定を行う命令をそなえたプロセッサも増えてきました
（SSEやHPC-ACE）。ただし、初期値の精度（有効ビット数）はハードウェア依存、
最終的に必要とされる精度はアプリケーション依存です。

ここでは、「Newton-Raphson法を倍精度に収束するまで闇雲に繰り返すよりもう少し効率のいい方法はないか」という話題を取り扱っていきます。

**逆数 [#m5a0b70b]
#math( y_n \sim 1 / x, )
#math( y_n \simeq 1 / x, )
として、二次収束：
#math( y_{n+1}  = 2 y_n - x y_n^2. )
 add, mul, mulsub
と覚えるとよいでしょう（２倍を加算で実現すれば定数が不要）。

より高次のものは、
#math(h_n \equiv 1 - x y_n)
としてて、
として、
三次収束：
#math( y_ {n+1} = (1 + h_n(1  + h_n)) y_n, )
四次収束は因数分解できて、
四次収束は因数分解できて：
#math( y_{n+1} = (1 + h_n)(1 + h_n^2) y_n, )
八次収束：
#math( y_{n+1} = (1 + h_n)(1 + h_n^2)(1 + h_n^4) y_n, )
FMA dual-issueの計算機を考えると、
 h    = 1.0 - x*y,  hp1  = 2.0 - x*y;
 h2   = h*h,        tmp1 = hp1*y;
 h2p1 = 1.0 + h2,   h4p1 = 1.0 + h2*h2;
 tmp2 = h2p1*h4p1;
 y    = tmp1*tmp2;
のように実装でき、スループットベースでは5 cycle程度の消費となる。

#math( y_n = \frac1x (1 - h_n), )
と言う事を考えれば上記の公式はすんなりと納得出来ると思う。

**逆数平方根 [#ac26a679]
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS