빠른 1/sqrt(n) 계산방법

http://en.wikipedia.org/wiki/Fast_inverse_square_root

사실 애니메이션 본TM계산시 엄청나게 많이 쓰는데 당연한 얘기지만 각종 노말라이즈에 다 들어가는 수식이기 때문이가.

그런데 생각보다 속도가 느리다... 이를 최적화하는 기법. 

최초 구현된 것은 ID의 Quake 3 Arena에서 존카멕 형이 만든 것. 아래의 소스 코드이다.

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;
 
	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
 
	return y;
}

아래는 언리얼4엔진에서 사용하고 있는 방법. SSE 를 사용하고 있다.

static FORCEINLINE float InvSqrt( float F )

{

static const __m128 fThree = _mm_set_ss( 3.0f );

static const __m128 fOneHalf = _mm_set_ss( 0.5f );

__m128 Y0, X0, Temp;

float temp;

Y0 = _mm_set_ss( F );

X0 = _mm_rsqrt_ss( Y0 ); // 1/sqrt estimate (12 bits)

// Newton-Raphson iteration (X1 = 0.5*X0*(3-(Y*X0)*X0))

Temp = _mm_mul_ss( _mm_mul_ss(Y0, X0), X0 ); // (Y*X0)*X0

Temp = _mm_sub_ss( fThree, Temp ); // (3-(Y*X0)*X0)

Temp = _mm_mul_ss( X0, Temp ); // X0*(3-(Y*X0)*X0)

Temp = _mm_mul_ss( fOneHalf, Temp ); // 0.5*X0*(3-(Y*X0)*X0)

_mm_store_ss( &temp, Temp );

return temp;

}