BigFloat로 Pi를 구해보자-6. 나눗셈 구현(Newton-Raphson 방법)

Posted by 적분 ∫2tdt=t²+c
2012.10.28 22:57 프로그래밍/Multi precision



저번 연재에서는 Goldschmidt의 방법을 이용해서 나눗셈을 구현했었는데요, 이번에는 Newton-Raphson의 방법을 이용해서 구현해보려고 합니다. 그리고 그 둘 중에서 어떤 놈이 빠른지 확인해보는거지요.


Newton-Raphson 방법에서는 초기값이 중요합니다. 초기값을 잘못 잡으면 아예 수렴 방향이 바뀌어서 엉뚱한 곳으로 수렴할 수가 있거든요. (기억이 나지 않는다면 BigFloat로 Pi를 구해보자-4.나눗셈은 어려워 (나눗셈 구현 알고리즘) 참고)



초기값을 적당히 잡아서 i를 키워나가면 Xi는 1/D에 수렴하겠지요. 그러면 초기값은 어떻게 잡을까요? D의 범위를 0.5 < D < 1로 잡았으니, 1/D의 범위는 1 < 1/D < 2 가 될겁니다.


그냥 단순하게 1차 선형으로 근사를 시켜보면, D = 0.5일때 1/D = 2이고, D = 1일때 1/D = 1이므로

라고 할 수 있겠네요. 평균오차는 ((1/D) - (3-2D))^2를 0.5에서 1까지 정적분하고 루트를 씌우면 되겠지요. 평균오차를 최소화시키기 위한 간단한 일차식을 찾아보면 다음을 구할수 있습니다.

(위와 같은 근사가 어떻게 구해질까요?)


자, 이제 다 됐습니다. x_0를 1/D의 근사값으로 넣고 계산을 시작하면 되겠네요.




코드를 볼까요?



먼저 코드를 보면 아직 구현하지 않은 muladd와 mulsub가 등장하는데요, 이는 Fused Multiply-Add 연산을 흉내내는 함수입니다. 즉 muladd(a, b, c)는 a+=b*c (a에 b*c를 합함), mulsub(a, b, c)는 a-=b*c(a에 b*c를 뺌)이지요. 이 함수 구현은 의외로 간단합니다.



RR_mul함수 구현과 거의 동일합니다. 즉 추가적인 연산비용없이 곱셈+덧셈을 한번에 수행할수 있지요.

다만 mulsub구현은 조금 복잡한데, 곱셈 결과를 빼나가다 보면 dest의 부호가 바뀌는 경우가 있기때문에 그렇지요. 여기에 필요한 i_subu64 함수는 i_addu64함수(3. 곱셈 구현 + 10진수로 출력하기)와 같은 방법으로 구현할 수 있습니다.




자 이렇게 두 가지 방법으로 나눗셈 함수를 구현해봤습니다. 그러면 어떤게 나은지 성능 테스트를 해봐야지요??





결과는 Newton-Raphson방법의 승리네요!

더욱 놀라운건 두 방법의 오차가... 0 ! 제로입니다. 아주 완벽하네요ㅋ

(하지만 Newton_Raphson 방법의 속도를 훨씬 끌어올릴만한 최적화가 아직 남아있습니다. 다음시간에 계속...)

Tags
이 댓글을 비밀 댓글로