고정 헤더 영역
상세 컨텐츠
본문
저번 연재에서는 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의 근사값으로 넣고 계산을 시작하면 되겠네요.
코드를 볼까요?
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 |
int RR_inv(PRR dest, PRR a){ if(RR_iszero(a)) return -1; PRR t = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); PRR na = RR_create(dest->size + 1); RRU32 w = i_bscs(a->frac, a->size);//a를 normalize하여 na에 저장. if(w < 32) { RR_shrcopy(na, a, 32 - w); } else if(w > 32) { RR_shlcopy(na, a, w - 32); } else { RR_move(na, a); } RR_initu(dest, 2, 0); dest->frac[1] = 3537031890; // dest = 2.823529 RR c1882; static RRU32 c1882d[2] = {1, 3789677025}; c1882.size = 2; c1882.alloced = 2; c1882.sign = 0; c1882.frac = c1882d; RR_mulsub(dest, na, &c1882); // dest = 2.823529 - na * 1.882353// 1/na의 초기값을 계산합니다. RRU32 dsize = dest->size; while(1) { RR_initu(t, 1, 0); RR_mulsub(t, na, dest); // t = 1 - na*dest if(i_cmpzero(t->frac, dsize))break; RR_mul(u, t, dest); // u = dest*(1 - na*dest) i_roundup(u->frac, u->alloced); RR_add(dest, u); // dest += u } if(w < 32) { RR_shr(dest, 32 - w); } else if(w > 32) { RR_shl(dest, w - 32); }// a를 시프트한만큼 결과값도 시프트해주어야// 제대로된 결과가 나옴. RR_delete(na); RR_delete(u); RR_delete(t); return 0;} |
먼저 코드를 보면 아직 구현하지 않은 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를 뺌)이지요. 이 함수 구현은 의외로 간단합니다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
inline void i_RR_muladd(PRR dest, PRR a, PRR b){ for(RRU32 j = b->size - 1; j != (RRU32)(-1); --j) { for(RRU32 i = a->size - 1; i != (RRU32)(-1); --i) { if(a->frac[i] && b->frac[j]) { i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b->frac[j]), i + j); } } }}//3번째 연재했던 RR_mul함수와 동일ㅋinline void i_RR_mulsub(PRR dest, PRR a, PRR b){ bool invert = false; for(RRU32 j = b->size - 1; j != (RRU32)(-1); --j) { for(RRU32 i = a->size - 1; i != (RRU32)(-1); --i) { if(a->frac[i] && b->frac[j]) { if(invert) {//부호가 반전된 이후로는 빼면 안되고 더해야함. i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b->frac[j]), i + j); } else { if(i_subu64(dest->frac, dest->alloced, __emulu(a->frac[i], b->frac[j]), i + j))//만약 언더플로우가 일어난 경우는 부호가 뒤집혀야 되므로 { dest->sign = ! dest->sign; i_cmplt(dest->frac, dest->size); invert = true; } } } } }}void RR_muladd(PRR dest, PRR a, PRR b){//dest와 a*b의 부호가 같으면 덧셈 if(dest->sign == (a->sign ^ b->sign)) { i_RR_muladd(dest, a, b); }//다르면 뺄셈 수행 else { i_RR_mulsub(dest, a, b); }}void RR_mulsub(PRR dest, PRR a, PRR b){//dest와 a*b의 부호가 같으면 뺄셈 if(dest->sign == (a->sign ^ b->sign)) { i_RR_mulsub(dest, a, b); }//다르면 덧셈 수행 else { i_RR_muladd(dest, a, b); }} |
RR_mul함수 구현과 거의 동일합니다. 즉 추가적인 연산비용없이 곱셈+덧셈을 한번에 수행할수 있지요.
다만 mulsub구현은 조금 복잡한데, 곱셈 결과를 빼나가다 보면 dest의 부호가 바뀌는 경우가 있기때문에 그렇지요. 여기에 필요한 i_subu64 함수는 i_addu64함수(3. 곱셈 구현 + 10진수로 출력하기)와 같은 방법으로 구현할 수 있습니다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
inline int i_subu32(RRU32* dest, RRU32 len, RRU32 src, RRU32 pos){ int carry = 0; RRU32 i = pos; if(dest[i] < src) carry = 1; dest[i--] -= src; for(; i != (RRU32)(-1); --i) { if(carry) { carry = dest[i] == 0; dest[i]--; } } return carry;}inline int i_subu64(RRU32* dest, RRU32 len, RRU64 src, RRU32 pos){ if(len < pos) return 0; if(len == pos) { return i_subu32(dest, len, (RRU32)(src >> 32), pos - 1); } int carry = 0; RRU32 i = pos; carry = dest[i] < (RRU32)(src & 0xFFFFFFFF); dest[i--] -= (RRU32)(src & 0xFFFFFFFF); if(i == (RRU32)(-1)) return carry; if(carry) { dest[i]--; carry = dest[i] == (RRU32)(-1) || dest[i] < (RRU32)(src >> 32); } else { carry = dest[i] < (RRU32)(src >> 32); } dest[i--] -= (RRU32)(src >> 32); for(; i != (RRU32)(-1); --i) { if(carry) { carry = dest[i] == 0; dest[i]--; } } return carry;} |
자 이렇게 두 가지 방법으로 나눗셈 함수를 구현해봤습니다. 그러면 어떤게 나은지 성능 테스트를 해봐야지요??
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
#include "stdafx.h"#include "RealReal.h"#include <Windows.h>#include <MMSystem.h>#pragma comment(lib, "winmm.lib")int _tmain(int argc, _TCHAR* argv[]){ DWORD t; PRR a,b,c,d; a = RR_create(128); b = RR_create(1024); c = RR_create(1024); while(1) { RR_init(a, 1234567890); t = timeGetTime(); for(int i = 0; i < 32; ++i) { RR_inv(b, a); } printf("Newton-Raphson's : %dms\n", timeGetTime() - t); RR_print(b, stdout, 64); printf("\n"); RR_init(a, 1234567890); t = timeGetTime(); for(int i = 0; i < 32; ++i) { RR_inv2(c, a); } printf("Goldschmidt's : %dms\n", timeGetTime() - t); RR_print(c, stdout, 64); printf("\nerror: %d %d\n\n", b->frac[b->size-2] - c->frac[c->size-2], b->frac[b->size-1] - c->frac[c->size-1]); } printf("\n");//메모리 해제는 귀찮아서 안했음ㅋ while(1); return 0;} |
결과는 Newton-Raphson방법의 승리네요!
더욱 놀라운건 두 방법의 오차가... 0 ! 제로입니다. 아주 완벽하네요ㅋ
(하지만 Newton_Raphson 방법의 속도를 훨씬 끌어올릴만한 최적화가 아직 남아있습니다. 다음시간에 계속...)
'프로그래밍 > Multi precision' 카테고리의 다른 글
| BigFloat로 Pi를 구해보자-8. 어셈블리어로 최적화하기 (0) | 2012.11.05 |
|---|---|
| BigFloat로 Pi를 구해보자-부. 나눗셈 구현 코드 (0) | 2012.11.04 |
| BigFloat로 Pi를 구해보자-7. 나눗셈 최적화와 큰O 표기법(Big-O Notation) (0) | 2012.11.01 |
| BigFloat로 Pi를 구해보자-5. 나눗셈 구현 (Goldschmidt division) (0) | 2012.10.24 |
| BigFloat로 Pi를 구해보자-4. 나눗셈은 어려워 (나눗셈 구현 알고리즘) (0) | 2012.10.24 |
| BigFloat로 Pi를 구해보자-3. 곱셈 구현 + 10진수로 출력하기 (0) | 2012.10.18 |
댓글 영역