高速根号計算 (fast sqrt algorithm)

　概要: C言語のsqrt(float)より精度は若干劣るものの，2倍以上速いsqrtのalgorithm．
　ググって見つけた物が，非常に面白かったのでまとめておく．
　精度より速度が求められる場面で活躍する（シェーダ-とか）．
　精度はあまりでないのでシミュレーションなどには使えない．

　参考文献
　[1] David Eberly, Fast Inverse Square Root (Revisited), http://www.geometrictools.com/Documentation/ FastInverseSqrt.pdf, 1/2002-.

実装

//---Algorithm float(IEEE754)用------
inline float t_sqrtF( const float& x )
{
　　float xHalf = 0.5f * x;
　　int   tmp   = 0x5F3759DF - ( *(int*)&x >> 1 ); //initial guess
　　float xRes  = *(float*)&tmp;

　　xRes *= ( 1.5f - ( xHalf * xRes * xRes ) );
　　//xRes *= ( 1.5f - ( xHalf * xRes * xRes ) );//コメントアウトを外すと精度が上がる
　　return xRes * x;
}

//---Algorithm double(IEEE754)用------
inline double t_sqrtD( const double &x) 
{
　　double         xHalf = 0.5 * x;
　　long long int  tmp   = 0x5FE6EB50C7B537AAl - ( *(long long int*)&x >> 1);//initial guess
　　double         xRes  = * (double*)&tmp;

　　xRes *= ( 1.5 - ( xHalf * xRes * xRes ) );
　　//xRes *= ( 1.5 - ( xHalf * xRes * xRes ) );//コメントアウトを外すと精度が上がる
　　return xRes * x;
}

説明

　　　この二つのアルゴリズムはそれぞれ，float(IEEE 754), double(IEEE754)用のsqrt(x)を求める関数である．
　　　[1]によれば，2002/1/9にcom.graphics.algorithmに投稿された物らしい.
　　　(doubleバージョンは私がweb上の情報を元に適当に書き換えた物で，Visual Studio 2010で動作を確認した．)
　　　内容自体は非常に単純で， $\frac{1}{\sqrt{x}}$ をニュートン法で求め， $\frac{x}{\sqrt{x}} = \sqrt{x}$ を計算している．

　　　ニュートン法で1/√a を求めるには， $f(x) = x^{-2} - a$ と置き， $f(x) = 0$ となる x を逐次的に求める．
　　　ニュートン法のイテレーションは　
　　　　　 $\displaystyle x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$
　　　
　　　なので，次のようにxnを繰り返し更新すればよい．
　
　　　　　 $\displaystyle x_{k+1} = x_k - \frac{x_k^{-2}-a}{-2x_k^{-3}} = x_k (1.5-0.5ax_k^2)$
　
　　　ここで重要なのは，このイテレーションには，足し算と掛け算しか含まれておらず，速度低下の原因になる割り算がない点である．
　　　また，上のアルゴリズムをみると，このニュートン法のイテレーションを一度しかまわしていない．
　　　しかし，例えばsqrt(2.0)を計算すると，以下の通り，小数点第2位まで正解している．

   sqrt(2.0 ) = 1.4142135623730950(解析解)
　 sqrt(2.0f) = 1.4142135381698608(VCのsqrt())
　 sqrt(2.0 ) = 1.4142135623730951(VCのsqrt())
t_sqrtF(2.0f) = 1.4138600826263428
t_sqrtD(2.0 ) = 1.4138593015909278

　　　ちなみに，2回イテレーションを回すと，以下の通り，少数第6位まで同じ値になり，

   sqrt(2.0 ) = 1.4142135623730950(解析解)
　 sqrt(2.0f) = 1.4142135381698608(VCのsqrt())
　 sqrt(2.0 ) = 1.4142135623730951(VCのsqrt())
t_sqrtF(2.0f) = 1.4142132997512817
t_sqrtD(2.0 ) = 1.4142134292706141

　　　3回イテレーションを回すと，以下の通り，floatでは誤差のためか精度向上が見られないが，doubleでは大きな精度向上がみられた．

   sqrt(2.0 ) = 1.4142135623730950(解析解)
　 sqrt(2.0f) = 1.4142135381698608(VCのsqrt())
　 sqrt(2.0 ) = 1.4142135623730951(VCのsqrt())
t_sqrtF(2.0f) = 1.4142136573791504
t_sqrtD(2.0 ) = 1.4142135623730765

　　　計算時間は，2回イテレーションを回しても倍程度速かった(Visual C++, Intel Core i7).

初期値 - Initial Guess

　　　このアルゴリズムで非常に重要かつ面白いのは，初期解を求める以下の行である．

　　int   tmp   = 0x5F3759DF - ( *(int*)&x >> 1 ); //initial guess

　　　この一行に，実はかなり深い意味が隠されている. 以下に説明を書くが，より詳しくは[1]を参照．

　準備

　　　IEEE 754 float型は，32bitの浮動小数点数で，1bitの符号ビット，8bitの指数ビット，23bitの仮数ビットからなる．
　　　例えば，
　　　　　　c = 0101 1111 0011 0000 0000 0000 0000 0000
　　　の表す値を考える．

　　　符号 : 0 = 正
　　　指数 : 101 1111 0 = 128 + 32 + 16 + 8 + 4 + 2 = 190
　　　仮数 : 011 0000 0000 0000 0000 0000 = 0/2 + 1/2^2+ 1/2^3+ ... + 1/2^23= 0.375

　　　指数部には，-127の下駄をはかせ，仮数部には1.0を足す必要がある事を考慮すると，cは次のようになる．
　　　　　　c = (1.0 + 0.375) * 2^(190-127)

　初期の計算法(指数部分)

　　　1/√aをニュートン法で求める場合，欲しい初期値は1/√aそのものである．
　　　この近似値を，√や/を使わないで求めたい．

　　　aはfloatで与えられるので，その指数部の値をE，仮数部をTとすると......
　　　　　 $\displaystyle a = 2 ^{E-127}\time(1.0 + T)$
　　　すると......
　　　　　 $\displaystyle \frac{1}{\sqrt{a}} = a ^{-1/2} = 2^{-E/2 + 63} \time \frac{1}{\sqrt{1.0+T}}$
　　　となる.
　　　求めたい数 -- 1/√aの仮数部E'は，- E/2 + 63 + 127 - 127 = (190 - E/2) - 127であるので，

　　　　　　E' = 190 - E/2

　　　となる．つまり，指数部分に関しては、aを1bit右にシフトさせ(a>>1)，これを190から引けばよい．つまり，

　　　　　　1/√a ≒ 0101 111 0*** **** **** **** **** **** - (a>>1) = 0x5f****** - (a>>1)

　　　*は仮数部分．指数部分に関しては，ビット演算のみで正確な近似解が得られる．

　初期の計算法(仮数部分)

　　　仮数部分23bitは，テーラー展開などで近似するのでなく，力技で求める．
　　　つまり、前もって全通り（2^24通り）試しておいて、一番誤差の少ないものを初期値とする．

　　　入力されうるすべての数 a を考える
　　　　　　符号部分と指数部分の上位7bitを固定し，下位24ビット全ての組み合わせを考慮した2^24通りのa=0x3f ** ** ** を考える．
　　　　　　ただし, 　1/2 <= a < 2.0 となる

　　　初期解の引かれる数 c のすべての通りを考える
　　　　　　符号部分と指数部分の上位7bitを固定し、2^24通りの c = 0x 5f ** ** ** を用意する．

　　　今，cを適当に一つ選ぶと，全ての入力aに対する最大誤差が計算できる．

　　　　　　 $\displaystyle \max_{a} |sqrt(a) - (c-(a>>1))|$

　　　全てのcにおいて上の最大誤差を計算し，それを最小化するcを選び出せばよい.

　　　　　　 $\displaystyle \min_{c} \max_{a} |sqrt(a) - (c-(a>>1))|$
　　　この計算により，c_minmax = 0x537642fが得られるが[1]，この値はオリジナル c_orig=0x5f3759dfと異なる．
　　　[1]によると，c_minmaxを用いた方が初期値のエラーは小さいが，c_origを用いた方がニュートン法を一ステップ回した後のエラーが小さくなる．
　　　c_origは，ニュートン法1ステップ込みの最大誤差の最小化を行って求めたものだと考えられる．

　　　解析: 0x5f3759df - (a>>1)という初期値は，指数部分の桁落ちがある場合/無い場合，仮数部分の引き算が正の場合/負の場合で，4通りの異なる直線を描く(ありえない1通りがあるので実際は3本)．
　　　つまり，3本のpiece-wise-linearなグラフになり，これが非常によく1/√aにフィットする([1]参照)．