제곱근(sqrt) 함수를 구현해보자.

우리가 프로그래밍을 하면서 자주 사용하는 여러 수학 함수들이 있다. 할일 없는 사람들은 이 함수들이 어떻게 작동하는지 정말 궁금해할지도 모른다. (물론 나도 그런 할 일 없는사람이다.) 한 번 내가 알고 있는 것들을 글로 써보며 생각을 정리해보고자 한다.

수학함수를 구현할때 가장 중요한 점은 정밀도와 속도이다.
늘상 있는 일이지만 프로그래머는 정확성과 속도 사이를 잘 조율해야한다. 필요이상의 정밀도는 시간 낭비일뿐이다.

먼저 수학함수 중에서 가장 인기가 높은 제곱근 함수 (sqrt함수)
거리 구할때를 비롯하여 안 쓰이는 곳이 없는 인기 좋은 함수이다.
한번 float sqrt(float x)를 구현해보자.

부동소수점 데이터는 꼴로 표현된다. (음수일 경우와 0일 경우는 무시했다.)
여기서 두 가지 경우로 나누어보자.

첫번째 n=2k 인 경우

가 된다.

두번째 n=2k+1 인 경우

가 된다.
는 상수이므로, 컴파일 타임에 결정되므로 크게 상관하지 않아도 좋다.
여기서 중요한것은 실질적으로 우리가 계산해야할 값은 라는 것이다.
k는 n을 2로 나눈값이므로 시프트 연산만으로 구현할 수 있기 때문이다.

를 어떻게 구할까?
이를 예견한것인지는 몰라도 뛰어난 수학자들이 몇 세기전에 좋은 해결방법을 만들어 놓았다.
그 이름은 바로 테일러 급수.


단 ( |x| < 1 이어야한다.)
자세한 것은 위키백과 테일러 급수 참조

라고 하면 가 된다.
따라서 위 테일러 전개식의 수렴조건도 만족한다. 적당한 부분까지 테일러 전개식에 대입하여 를 구해내면된다.
게다가 부동소수점 데이터의 특성상 곧바로 x값을 구할수 있다.

32비트 부동소수점 데이터를 예로 들자면
0 00000000 00000000000000000000000
부호  지수부  가수부
그런데 가수부는 항상 가장 첫째자리의 1을 뺀 나머지 부분을 저장한다.
예를 들어 1.25 * 2^8 같은 경우
부호는 0 (양수이므로)
지수부는 10001000 (바이어스로 128을 더해주어서 136)
가수부는 01000000000000000000000 (1.25 이므로 1.01₍₂₎이나, 앞의 1을 빼버린다.)

그리고 위 식을

이렇게 묶어내면 연산횟수를 줄일 수 있다.

이 정도면 꽤 최적화된 제곱근 계산하는 함수를 구현할수 있을것이다.

float sqrt(float x)
{

	DWORD *px=(DWORD*)(&x);
	//음수일 경우 NAN을 돌려준다.
	if((*px)&0x80000000)return NAN;
	//0일 경우 0을 돌려준다.
	if((*px)==0)return 0.f;
	//지수부를 구한다.
	long exp=(((*px)&0x7F800000)>>23)-128;
	//가수부를 구한다.
	__int64 signif=(*px)&0x7FFFFF,temp;
	//1-(5/8)x
	temp=(1<<24-(5*(signif>>3)));
	//1-(1/2)x(1-(5/8)x)
	temp=1<<24-(((signif>>1)*temp)>>23);
	//1-(1/4)x(1-(1/2)x(1-(5/8)x))
	temp=1<<24-(((signif>>2)*temp)>>23);
	//(1/2)x(1-(1/4)x(1-(1/2)x(1-(5/8)x)))
	temp=1<<24-(((signif>>1)*temp)>>23);
	float ret;
	*(DWORD*)&ret=(exp/2+128)<<23 | (signif&0x7FFFFF);
	//지수부가 홀수인경우
	if(exp&1)return 1.4142136f*ret;
	else return ret;
}