« 逆数の高精度化 | メイン | テクスチャマッピングの高速化 その3 »

2008年02月23日

x86 SIMD Technique:: 平方根の逆数の高精度化

    

平方根の逆数を求める命令として、_mm_rsqrt_ss / _mm_rsqrt_ps がある。
でも、これは 11bit 精度しかない ( 11bit 精度だと、10 進数で3桁程度 ) 。
と言うことで、11bit 精度で得られた平方根の逆数の近似値をニュートン-ラフソン法 ( Newton-Raphson ) を用いて高精度化する。
SSE のソースコードは以下。

// 22bit 精度で平方根の逆数を求める。
// 1/√x = 3/2 * x0 - 1/2 * a * x0^3 を用いる
inline __m128 m128_rsqrt_22bit_ps( const __m128& a ) {
  static const __declspec(align(16)) float valm05[4] = {-0.5f,-0.5f,-0.5f,-0.5f};
  static const __declspec(align(16)) float val15[4] = {1.5f,1.5f,1.5f,1.5f};
  static const __m128 Q_m0p5 = *(__m128*)valm05;
  static const __m128 Q_1p5 = *(__m128*)val15;
  __m128 xm0 = a;
  __m128 xm1 = _mm_rsqrt_ps( xm0 );
  __m128 xm2 = _mm_mul_ps( xm0, xm1 );
  xm2 = _mm_mul_ps( xm2, xm1 );
  xm2 = _mm_mul_ps( xm2, xm1 );
  xm2 = _mm_mul_ps( xm2, Q_m0p5 );
  xm0 = _mm_mul_ps( xm1, Q_1p5 );
  return _mm_add_ps( xm0, xm2 );
}
inline __m128 m128_rsqrt_22bit_ss( const __m128& a ) {
  // 必要ないけど、4つセット
  static const __declspec(align(16)) float valm05[4] = {-0.5f,-0.5f,-0.5f,-0.5f};
  static const __declspec(align(16)) float val15[4] = {1.5f,1.5f,1.5f,1.5f};
  static const __m128 Q_m0p5 = *(__m128*)valm05;
  static const __m128 Q_1p5 = *(__m128*)val15;
  __m128 xm0 = a;
  __m128 xm1 = _mm_rsqrt_ss( xm0 );
  __m128 xm2 = _mm_mul_ss( xm0, xm1 );
  xm2 = _mm_mul_ss( xm2, xm1 );
  xm2 = _mm_mul_ss( xm2, xm1 );
  xm2 = _mm_mul_ss( xm2, Q_m0p5 );
  xm0 = _mm_mul_ss( xm1, Q_1p5 );
  return _mm_add_ss( xm0, xm2 );
}
inline float rsqrt_sse( float a ) {
  float  ret;
  __m128 xm0 = _mm_set_ss(a);
  _mm_store_ss( &ret, m128_rsqrt_22bit_ss( xm0 ) );
  return ret;
}

これと逆数の高精度化を用いれば、平方根が得られるが、平方根は素直に _mm_sqrt_ss / _mm_sqrt_ps を使った方が速い様子。
また、標準関数の実装にもよると思うが、以下のようにして SSE の平方根を使う方が速いこともあるようだ ( VC2005 では以下のほうが速かった ) 。
精度も 23bit なので、困ることはないと思う。

inline float sqrt_sse( float a ) {
  float  ret;
  __m128 xm0 = _mm_set_ss(a);
  _mm_store_ss( &ret, _mm_sqrt_ss( xm0 ) );
  return ret;
}

平方根の逆数は距離で割る時など、初めから逆数になっているのでかけるだけで済む。

追記:
1/√a * a = √a なので、上の平方根の逆数に逆数云々ではなく、単純に平方根の逆数に元の値を掛ければ平方根が求まる。
こうすれば上の高精度化した平方根の逆数に元の数を掛けた方が _mm_sqrt_ss / _mm_sqrt_ps よりも倍近く速い。
2 ~ 1千万までの平方根を求める時間を計測すると、VC 2005 の sqrt、_mm_sqrt_ss 使用版、上の高精度化版*元の数 では、だいたい 3 : 2 : 1 ぐらいの速度比になった。
でも、状況によってはコンパイラの最適化との兼ね合いか、 _mm_sqrt_ss / _mm_sqrt_ps を使う方が速いこともあるよう。詳しくは調べてない。

さらに追記:
( 2 ~ 1千万 ) / 100.0 で平方根を求める時間を計測すると、45 : 33 : 35 ぐらいになった。
小数の時は _mm_sqrt_ss / _mm_sqrt_ps の方が速いようだ。
_mm_rsqrt_ss / _mm_rsqrt_ps はテーブルを引くらしいのでその辺りの兼ね合いか?



投稿者 Takenori : 2008年02月23日 20:59




comments powered by Disqus
Total : Today : Yesterday : なかのひと