상세 컨텐츠

본문 제목

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

프로그래밍/테크닉

by ∫2tdt=t²+c 2008. 11. 26. 00:45

본문

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

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

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

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

첫번째 n=2k 인 경우
가 된다.

두번째 n=2k+1 인 경우
가 된다.
http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=%5Csqrt%7B2%7D는 상수이므로, 컴파일 타임에 결정되므로 크게 상관하지 않아도 좋다.
여기서 중요한것은 실질적으로 우리가 계산해야할 값은 http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=%5Csqrt%7B%5Calpha%7D라는 것이다.
k는 n을 2로 나눈값이므로 시프트 연산만으로 구현할 수 있기 때문이다.

http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=%5Csqrt%7B%5Calpha%7D를 어떻게 구할까?
이를 예견한것인지는 몰라도 뛰어난 수학자들이 몇 세기전에 좋은 해결방법을 만들어 놓았다.
그 이름은 바로 테일러 급수.

http://www.sitmo.com/gg/latex/latex2png.2.php?z=150&eq=%5Csqrt%7B1%20%2B%20x%7D%20%3D%201%20%2B%20%5Ctextstyle%20%5Cfrac%7B1%7D%7B2%7Dx%20-%20%5Cfrac%7B1%7D%7B8%7Dx%5E2%20%2B%20%5Cfrac%7B1%7D%7B16%7D%20x%5E3%20-%20%5Cfrac%7B5%7D%7B128%7D%20x%5E4%20%2B%20%5Cdots%5C!
단 ( |x| < 1 이어야한다.)
자세한 것은 위키백과 테일러 급수 참조

http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=%5Calpha%3D1%2Bx라고 하면 http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=0%5Cle%20x%3C1가 된다.
따라서 위 테일러 전개식의 수렴조건도 만족한다. 적당한 부분까지 테일러 전개식에 대입하여 http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=%5Csqrt%7B%5Calpha%7D를 구해내면된다.
게다가 부동소수점 데이터의 특성상 곧바로 x값을 구할수 있다.

그리고 위 식을
http://www.sitmo.com/gg/latex/latex2png.2.php?z=100&eq=1%2B%5Cfrac%7B1%7D%7B2%7Dx(1-%5Cfrac%7B1%7D%7B4%7Dx(1-%5Cfrac%7B1%7D%7B2%7Dx(1-%5Cfrac%7B5%7D%7B8%7Dx(%5Cldots))))
이렇게 묶어내면 연산횟수를 줄일 수 있다.

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


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; }


관련글 더보기

댓글 영역