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

저번 연재에서는 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까지 정적분하고 루트를 씌우면 되겠지요. 평균오차를 최소화시키기 위한 간단한 일차식을 찾아보면 다음을 구할수 있습니다.

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

생각보다 어려워서 더보기로 빼놨습니다.

0.5~1범위에서 1/x과 제일 근사하는 일차식을 찾아보도록 합시다.

이 두 식의 오차가 최소가 되도록하는 a,b를 구하는게 목적이 되겠지요. 오차는 다음과 같이 계산할 수 있습니다.

이 값이 최소가 되도록하여야 하므로, 전개하여 계산해봅시다.

(복잡고 귀찮은 중간 계산과정은 생략ㅋ)

이 식을 최소가 되도록하는 a,b를 찾아야하므로, 미분을 해봅시다.

E는 a,b에 관한 이차함수이므로, 극소값은 미분값이 0이 되는 지점에서 형성될 것입니다. 그러므로 다음과 같은 연립방정식을 세울수 있으므로

이 식을 풀면

을 얻을 수 있습니다. 유리수로 근사시키면 위 결과를 얻을수 있습니다.

....어렵죠

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

코드를 볼까요?

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를 뺌)이지요. 이 함수 구현은 의외로 간단합니다.

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진수로 출력하기)와 같은 방법으로 구현할 수 있습니다.

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

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

#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 방법의 속도를 훨씬 끌어올릴만한 최적화가 아직 남아있습니다. 다음시간에 계속...)