앞서서 덧셈과 뺄셈을 구현해보았는데요, 이미 둘을 구현했다면 곱셈을 구현하는 것은 생각보다 간단합니다. 알고리즘에 대해서는 다음을 참조하시면 됩니다.
2012/10/18 - [프로그래밍/Multi precision] - BigFloat로 Pi를 구해보자-3. 곱셈 구현 + 10진수로 출력하기
제일 단순한 방법으로 곱셈을 구현하면 n길이인 정수에 대해 O(n^2) 시간에 계산을 완료할 수 있습니다. 더 빠른 알고리즘 구현은 일단 추후로 미루고ㅋ
당연하겠지만 지수부와 정수배열로 표현된 두 수의 곱은 정수배열끼리는 곱하고 지수부는 더하면 됩니다. 부호연산은 자명하니 설명이 필요없구요. 바로 코드를 보시죠.
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 |
// xl_mul : a와 b의 곱을 구하여 dest에 저장한다 void xl_mul(PXLARGE dest, PXLARGE a, PXLARGE b) { if (xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest); XL_UINT32 aLen = _XL_GetLength(a); XL_UINT32 bLen = _XL_GetLength(b); xl_init(dest, 0); //연산 결과는 a의 길이와 b의 길이의 합입니다. //처음 크기가 1이었으므로 aLen + bLen - 1를 뒤쪽에 넣으면 //전체 길이는 aLen + bLen이 되겠죠. _XL_PushBackN(dest, aLen + bLen - 1); // [0, aLen)과 [0, bLen)인 두 정수에 대해 루프를 돌립니다. for (XL_UINT32 i = 0; i < aLen; ++i) for (XL_UINT32 j = 0; j < bLen; ++j) { XL_UINT32 ad = *_XL_GetFrac(a, i); XL_UINT32 bd = *_XL_GetFrac(b, j); // 둘 중에 하나라도 0이면 결과가 0이므로 더할 필요가 없음 if (ad && bd) { // 위에서 dest의 데이터가 일렬로 배치되도록 했기에, // 두 덩어리로 나뉘는 경우를 고려할 필요가 없습니다. i_addu64(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd)); } } _XL_Shrink(dest); dest->sign = a->sign ^ b->sign; // 부호는 xor연산 dest-> exp += a-> exp + b-> exp ; // 지수부는 덧셈 연산 } |
곱셈만 구현하고 끝낼수 없죠. 곱셈한 결과를 더해주는 FMA(Fused Multiply-Add)를 구현해봅시다. 이 연산이 은근히 유용합니다.
수식으로 정리하면 위 처럼 되겠죠? 실제 구현에서는 절대값을 더해주는 경우와 빼주는 경우로 나눠서 구현을 해야합니다.
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 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 |
inline void _Abs_MulAdd(PXLARGE dest, PXLARGE a, PXLARGE b) { XL_UINT32 aLen = _XL_GetLength(a); XL_UINT32 bLen = _XL_GetLength(b); XL_INT32 tExp = a-> exp + b-> exp ; XL_INT32 tBack = tExp + aLen + bLen + 1; // carry를 대비한 여분의 +1 if (dest-> exp > tExp) _XL_PushFrontN(dest, dest-> exp - tExp); XL_INT32 dBack = dest-> exp + _XL_GetLength(dest); if (dBack < tBack) _XL_PushBackN(dest, tBack - dBack); XL_UINT32 orgDestBegin = dest->begin; if (tExp > dest-> exp ) dest->begin += tExp - dest-> exp ; if (dest->begin > dest->size) dest->begin -= dest->size; // dest 데이터가 일렬로 배치된 경우 if (dest->begin < dest->end) { for (XL_UINT32 i = 0; i < aLen; ++i) for (XL_UINT32 j = 0; j < bLen; ++j) { XL_UINT32 ad = *_XL_GetFrac(a, i); XL_UINT32 bd = *_XL_GetFrac(b, j); if (ad && bd) { i_addu64(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd)); } } } // dest 데이터가 두 덩어리로 나눠진 경우 else { for (XL_UINT32 i = 0; i < aLen; ++i) for (XL_UINT32 j = 0; j < bLen; ++j) { XL_UINT32 ad = *_XL_GetFrac(a, i); XL_UINT32 bd = *_XL_GetFrac(b, j); if (ad && bd) { if (dest->begin + i + j < dest->size) { i_addu64_2(&dest->frac[dest->begin + i + j], dest->size - dest->begin - i - j, dest->frac, dest->end, __emulu(ad, bd)); } else { i_addu64(&dest->frac[dest->begin + i + j - dest->size], dest->end + dest->size - dest->begin - i - j, __emulu(ad, bd)); } } } } dest->begin = orgDestBegin; _XL_Shrink(dest); } inline void _Abs_MulSub(PXLARGE dest, PXLARGE a, PXLARGE b) { XL_UINT32 aLen = _XL_GetLength(a); XL_UINT32 bLen = _XL_GetLength(b); XL_INT32 tExp = a-> exp + b-> exp ; XL_INT32 tBack = tExp + aLen + bLen; if (dest-> exp > tExp) _XL_PushFrontN(dest, dest-> exp - tExp); XL_INT32 dBack = dest-> exp + _XL_GetLength(dest); if (dBack < tBack) _XL_PushBackN(dest, tBack - dBack); XL_UINT32 orgDestBegin = dest->begin; if (tExp > dest-> exp ) dest->begin += tExp - dest-> exp ; if (dest->begin > dest->size) dest->begin -= dest->size; // 부호가 뒤집혔는 지 여부를 나타내는 변수 bool invert = false ; // dest 데이터가 일렬로 배치된 경우 if (dest->begin < dest->end) { for (XL_UINT32 i = 0; i < aLen; ++i) for (XL_UINT32 j = 0; j < bLen; ++j) { XL_UINT32 ad = *_XL_GetFrac(a, i); XL_UINT32 bd = *_XL_GetFrac(b, j); if (ad && bd) { // 부호가 뒤집히지 않았으면 i_subu64, 뒤집혔으면 i_addu64 호출 if ((invert ? i_addu64 : i_subu64)(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd))) { // 오버플로우가 발생한 경우 부호가 뒤집힌것 invert = true ; dest->sign = !dest->sign; if (orgDestBegin < dest->end) { i_cmplt(&dest->frac[orgDestBegin], dest->end - orgDestBegin, 0); } else { int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0); if (dest->end) i_cmplt(dest->frac, dest->end, carry); } } } } } // dest 데이터가 두 덩어리로 나눠진 경우 else { for (XL_UINT32 i = 0; i < aLen; ++i) for (XL_UINT32 j = 0; j < bLen; ++j) { XL_UINT32 ad = *_XL_GetFrac(a, i); XL_UINT32 bd = *_XL_GetFrac(b, j); if (ad && bd) { if (dest->begin + i + j < dest->size) { // 부호가 뒤집히지 않았으면 i_subu64_2, 뒤집혔으면 i_addu64_2 호출 if ((invert ? i_addu64_2 : i_subu64_2)(&dest->frac[dest->begin + i + j], dest->size - dest->begin - i - j, dest->frac, dest->end, __emulu(ad, bd))) { // 부호 뒤집기 invert = true ; dest->sign = !dest->sign; int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0); if (dest->end) i_cmplt(dest->frac, dest->end, carry); } } else { // 부호가 뒤집히지 않았으면 i_subu64, 뒤집혔으면 i_addu64 호출 if ((invert ? i_addu64 : i_subu64)(&dest->frac[dest->begin + i + j - dest->size], dest->end + dest->size - dest->begin - i - j, __emulu(ad, bd))) { // 부호 뒤집기 invert = true ; dest->sign = !dest->sign; int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0); if (dest->end) i_cmplt(dest->frac, dest->end, carry); } } } } } dest->begin = orgDestBegin; _XL_Shrink(dest); } |
진짜 FMA를 하는 함수는 위 함수를 이용해서 다음과 같이 심플하게 짤수 있겠죠.
1 2 3 4 5 6 7 8 9 10 11 |
void xl_muladd(PXLARGE dest, PXLARGE a, PXLARGE b) { if (xl_isnan(dest) || xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest); return (dest->sign == (a->sign ^ b->sign) ? _Abs_MulAdd : _Abs_MulSub)(dest, a, b); } void xl_mulsub(PXLARGE dest, PXLARGE a, PXLARGE b) { if (xl_isnan(dest) || xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest); return (dest->sign != (a->sign ^ b->sign) ? _Abs_MulAdd : _Abs_MulSub)(dest, a, b); } |
이제 곱셈이 완료되었으니, 10진수로 출력하는 함수를 만들수 있습니다. 그건 다음으로 패스!
톰-쿡 알고리즘 개요 (0) | 2014.03.30 |
---|---|
[진짜 BigFloat 구현하기] 5. 카라슈바 알고리즘으로 곱셈 성능 향상 (1) | 2014.03.11 |
[진짜 BigFloat 구현하기] 4. 시프트 연산과 나눗셈(나머지) (0) | 2014.03.02 |
[진짜 BigFloat 구현하기] 2. 덧셈과 뺄셈 구현하기 (0) | 2013.11.25 |
[진짜 BigFloat 구현하기] 1. 기본 자료구조와 수치 표현 방법 정하기 (0) | 2013.11.20 |
[진짜 BigFloat 구현하기] 0. 시작에 앞서서 (0) | 2013.11.20 |
댓글 영역