[진짜 BigFloat 구현하기] 4. 시프트 연산과 나눗셈(나머지)

앞서서 곱셈을 구현했으니, 이제 나눗셈을 구현할 차례겠죠?

근데 문제는 나눗셈을 구현하는게 생각보다 까다롭다는거지요. 

일반적으로 우리가 학교에서 배우는 나눗셈 방법은 장제법(long division)이라고 하는데,

이게 얼마나 어려운지는 초등학생들이 더 잘 압니다.

http://bab2min.tistory.com/54

자세한 알고리즘에 대한 설명은 위 페이지로 대신합니다.

우리가 나눗셈을 할때 나누는 수랑 나눠지는 수 자리를 맞춰서 계산하죠. 자리 맞추기 위해 숫자들을 왼쪽으로 옮기고 오른쪽으로 옮기고 하는데, 2진법에서 이것과 근본적으로 같은 연산이 비트 시프트 연산입니다.

다른말로 하면 나머지를 구하는 나눗셈 연산을 위해서는 먼저 시프트 연산을 구현할 필요가 있다는 거지요.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

void xl_setnan(PXLARGE pxl); // pxl을 NaN으로 설정하는 함수 int xl_isnan(PCXLARGE pxl); // pxl이 NaN인지 확인하는 함수   int xl_cmp(PCXLARGE a, PCXLARGE b); // a > b 이면 1, a == b이면 0, a < b이면 -1을 반환하는 함수 int xl_cmpzero(PCXLARGE a); // a > 0 이면 1, a == 0이면 0, a < 0이면 -1을 반환하는 함수   void xl_shl(PXLARGE dest, XL_UINT32 bit); // dest를 왼쪽으로 bit비트만큼 시프트하는 함수 void xl_shr(PXLARGE dest, XL_UINT32 bit); // dest를 오른쪽으로 bit비트만큼 시프트하는 함수   void xl_mod(PXLARGE dest, PXLARGE divisor, PXLARGE quot); // 오늘의 목표! // dest 나누기 divisor를 수행해서 // 몫은 quot에 두고, 나머지를 dest에 넣어주는 함수입니다.

오늘 구현하게 될 함수들 목록입니다. 많아보이지만 별거없어요.

NaN이라는 개념이 등장했어요~! 난!이라는 게 아니라 Not A Number의 약자로 연산이 불가능하거나 연산결과가 실수 범위에서 나오지 않는 경우를 표현하기 위해 주로 부동소수점에서 사용되는 값입니다. 덧셈, 뺄셈, 곱셈 연산에서는 문제가 없었지만, 나눗셈이 등장하는 순간 0으로 나누기라는 함정이 생기게 되죠. 0으로 나누는게 불가능하다는 건 다들 알고계실테니 설명하지 않겟구요, 0으로 나눈 결과를 어떻게 표현하기는 해야하니 NaN으로 세팅해주는거죠.

일단 NaN관련 함수를 구현합니다.

1 2 3 4 5 6 7 8 9 10 11

void xl_setnan(PXLARGE pxl) {     pxl->sign = (XL_UINT32)-1; // 단순히 sign값을 -1인 경우를 NaN으로 하기로 약속하죠. }   int xl_isnan(PCXLARGE pxl) {     return pxl->sign == (XL_UINT32)-1; // 구현 참 쉽죠? }

먼저 시프트연산보다 앞서서 비교함수를 구현해보죠. 나눗셈 구현에는 비교연산도 필수적이니까요.

하나는 두 XLarge를 비교하는 함수구요, 다른 하나는 한 XLarge랑 0을 비교하는 함수입니다.

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

int xl_cmp(PCXLARGE a, PCXLARGE b) {     if(xl_isnan(a) || xl_isnan(b)) return -1; // 둘 중에 하나라도 NaN이면 무조건 -1을 반환합니다.     if(a->sign == 0)     {         if(b->sign == 0) return _Abs_Cmp(a, b); // 둘 다 양수라면 절대값을 비교합니다.         else return 1; // b가 음수라면 당연히 a가 크겠죠     }     else     {         if(b->sign == 0) return -1; // a가 음수인데 b가 양수면, 당연히 b가 크죠         else return -_Abs_Cmp(a, b); // 둘다 음수인 경우는 절대값을 비교하고 결과를 뒤집어 줍니다. // 음수면 절대값이 작을수록 큰거니까요     } }   int xl_cmpzero(PCXLARGE a) {     if(xl_isnan(a)) return -1; // NaN이면 -1 반환 끝.     if(*_XL_GetFracBack(a) == 0 && _XL_GetLength(a) == 1) return 0; // 값이 0일경우 전체 길이가 1이고, 가수부가 다 0이 되도록 약속해놨잖아요!     if(a->sign == 0) return 1; // 부호가 0이면 양수     return -1; // 아니면 음수 }

뭔가 없는거 같네요. 아 절대값 비교를 하는 _Abs_Cmp를 구현안했죠. 구현해보입시다.

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

inline int i_cmps2(const XL_UINT32* dest1, XL_UINT32 destSep, const XL_UINT32* dest2, const XL_UINT32* src1, XL_UINT32 srcSep, const XL_UINT32* src2, XL_UINT32 len) // 이 함수는 길이가 같은 두 정수 배열을 비교해줍니다. // 자료구조가 환형큐라서 중간에 끊길수 있기때문에 // dest배열은 dest1과 dest1의 길이인 destSep, dest2 // src배열은 src1과 src1의 길이인 srcSep, src2를 인수로 받습니다. // dest2와 src2의 길이를 받지 않는 이유는 전체 길이 len을 받으면 // destSep과 srcSep을 빼서 dest2, src2의 길이를 알아낼수 있기 때문이죠 { /* 당황하지 말고 침착하게 생각해보면 총 4가지의 경우가 나온다는걸 알수있어요. 비교할 대상이 1. 모두 dest1과 src1에 있다 (즉 len이 짧아서 destSep, srcSep을 넘지 못하는 경우) : d1s1 2. dest는 dest1에는 다 있지만 src는 src1과 src2에 나눠져있다. (len이 destSep은 넘지 않지만, srcSep은 넘는 경우) : d1s2 3. dest는 dest1과 dest2에 나눠져있지만 src는 src1에 다 있다. (len이 destSep은 넘지만, srcSep은 넘지 않는 경우) : d2s1 4. dest1, 2와 src1, 2에 모두 나눠져 있다. (len이 destSep도 넘고, srcSep도 넘는 경우)   일단 4번의 상태라고 가정을 하고 비교를 해서 jump해나가도록 합니다. */       XL_UINT32 i = len - 1;       if(i < destSep && i < srcSep) goto d1s1; // 1번 경우죠. d1s1로 goto!     if(i < destSep && i >= srcSep) goto d1s2; // 2번 경우죠. d1s2로 goto!     if(i >= destSep && i < srcSep) goto d2s1; // 3번 경우죠. d2s1로 goto! // 나머지는 4번 경우인거에요.       for(; i != destSep-1 && i != srcSep-1; --i)     {         XL_UINT32 dV = dest2[i - destSep];         XL_UINT32 sV = src2[i - srcSep];           if(dV < sV) return -1;         if(dV > sV) return 1;     } // 제일 윗 자리에서부터 비교해서 내려옵니다. // dest2, src2 비교가 끝날때까지도 return이 안된건 걔네들이 모두 같았다는거죠. // i값이 줄어들었으니 다시한번 점프를 해봅시다.     if(i < destSep && i < srcSep) goto d1s1; // 1번 경우죠. d1s1로 goto!     if(i < destSep) goto d1s2; // 2번 경우     if(i < srcSep) goto d2s1; // 3번 경우   d2s1: //3번 경우의 처리과정입니다.     for(; i != (XL_UINT32)(-1) && i != destSep-1; --i)     {         XL_UINT32 dV = dest2[i - destSep];         XL_UINT32 sV = src1[i];           if(dV < sV) return -1;         if(dV > sV) return 1;     } // dest2와 src1을 비교해서 내려옵니다. // 이 값들도 다 같으면, 다음에는 1번경우로 건너뜁니다.     goto d1s1;   d1s2: //2번 경우의 처리과정입니다.     for(; i != (XL_UINT32)(-1) && i != srcSep-1; --i)     {         XL_UINT32 dV = dest1[i];         XL_UINT32 sV = src2[i - srcSep];           if(dV < sV) return -1;         if(dV > sV) return 1;     } // dest1와 src2을 비교해서 내려옵니다. // 이 값들도 다 같으면, 다음에는 1번경우로 건너뜁니다.     goto d1s1;   d1s1:     for(; i != (XL_UINT32)(-1); --i)     {         XL_UINT32 dV = dest1[i];         XL_UINT32 sV = src1[i];           if(dV < sV) return -1;         if(dV > sV) return 1;     } // dest1와 src1을 비교해서 내려옵니다. // 이 값들도 다 같으면 두 배열이 아예 똑같은거죠.     return 0; }

(조금 복잡해보이죠? 하지만 5분만 생각해보면 당연하다는 걸 이해하실수 있을거에요.

goto를 적절히 사용하면 로직이 단순해지는 경우도 있습니다.)

이제 이 함수를 이용해서 _Abs_Cmp를 구현합니다.

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

inline int _Abs_Cmp(PCXLARGE a, PCXLARGE b) // 절대값을 비교하는 함수이므로 부호따위는 쿨하게 무시합니다. { // 지수부에 전체 자릿수 길이를 더하면 가장 높은 자릿수를 구할수 있죠. // 가장 높은 자릿수를 먼저 비교해봅시다.     XL_UINT32 aLen = _XL_GetLength(a);     XL_UINT32 bLen = _XL_GetLength(b);     if(a->exp   aLen < b->exp   bLen) return -1; // b가 자릿수가 더 높으니 b가 큰거     if(a->exp   aLen > b->exp   bLen) return 1; // a가 자릿수가 더 높으니 a가 큰거       // 끈질긴 녀석. 여기까지 왔다는건 최고 자릿수가 같다는거죠.     XL_INT32 mExp = MAX(a->exp, b->exp); // 일단 a와 b의 가수부 배열이 끊겨 있지 않다고 가정하며 시작합니다. // 둘 중에 지수부가 큰 녀석에 맞춰놓고 비교합시다. // 지수부가 큰데 최고 자릿수가 같다는건, 전체 길이가 짧다는거죠. // 즉 전체 길이를 구하려면 둘의 길이 중 짧은 걸 골라야한다는겁니다.       XL_UINT32* a1 = &a->frac[a->begin   mExp - a->exp];     XL_UINT32* a2 = nullptr;     XL_UINT32 aSep = MIN(aLen, bLen); // 짧은게 전체길이     XL_UINT32* b1 = &b->frac[b->begin   mExp - b->exp];     XL_UINT32* b2 = nullptr;     XL_UINT32 bSep = MIN(aLen, bLen); // 짧은게 전체길이   // a가 끊겨 있는 경우 처리     if(a1 > a->frac   a->size)     {         a1 -= a->size;     }     else if(a->begin > a->end)     {         a2 = a->frac;         aSep = a->size - (a->begin   mExp - a->exp);     }   // b가 끊겨 있는 경우 처리     if(b1 > b->frac   b->size)     {         b1 -= b->size;     }     else if(b->begin > b->end)     {         b2 = b->frac;         bSep = b->size - (b->begin   mExp - b->exp);     }   // 자, 여기까지해서 a1, aSep, a2, b1, bSep, b2를 모두 채웠습니다. // 이제 앞서 구현한 비교하는 함수에 넣고 비교를 해보죠.     int r = i_cmps2(a1, aSep, a2, b1, bSep, b2, MIN(aLen, bLen));     if(r) return r; // 이마저도 결과가 같게 나왔다면 최종적으로 a와 b의 지수부를 비교해봅니다. // 앞서 말했듯이 이미 여기까지 왔다는건 최고 자릿수가 같다는 거고 // 지수부가 높은 쪽이 전체 길이가 짧은거고, 그렇기에 전체 값이 작은겁니다. // 예로 123 * 10^2 과 1234 * 10^1 을 생각하면 됩니다.     if(a->exp > b->exp) return -1;     if(a->exp < b->exp) return 1; // 지수부까지도 같으면 완전 같은걸로 인정.     return 0; }

환형큐로 부동소수점 포맷을 구현하다보니 비교함수쪽이 매우 복잡하죠? 다시 정리하자면 비교 프로세스는 다음과 같습니다.

1. 최고 자릿수를 비교합니다. 최고 자릿수가 큰 쪽이 더 큰 수

2. 최고 자릿수가 같으면, 전체 길이가 짧은 쪽에 맞춰서 정수 배열을 비교합니다. 정수 배열이 큰 쪽이 더 큰수

3. 정수 배열도 같으면, 지수부를 비교합니다. 지수부가 작은 쪽이 더 큰 수 

4. 지수부도 같으면 두 수는 같은거.

꺄울! 수고많으셨습니다. 이제야... 시프트함수를 구현합니다. 먼저 도우미 함수 보고 가실게요.

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

inline XL_UINT32 i_shrs(XL_UINT32* dest, XL_UINT32 len, XL_UINT32 shift, XL_UINT32 fill) // len길이의 dest 배열을 오른쪽으로 shift비트만큼 시프트시킵니다. 최상위의 빈 자리는 fill로 채웁니다. // 최하위에서 오른쪽으로 밀려난 값이 결과값으로 반환됩니다 {     if(!shift) return 0;       XL_UINT32 rem = dest[0] << (32 - shift);     dest[0] = dest[0] >> shift;     for(XL_UINT32 i = 1; i < len;   i)     {         dest[i - 1] |= dest[i] << (32 - shift);         dest[i] = dest[i] >> shift;     }     dest[len - 1] |= fill;       return rem; }   inline XL_UINT32 i_shls(XL_UINT32* dest, XL_UINT32 len, XL_UINT32 shift, XL_UINT32 fill) // len길이의 dest 배열을 왼쪽으로 shift비트만큼 시프트시킵니다. 최하위의 빈 자리는 fill로 채웁니다. // 최상위에서 왼쪽으로 밀려난 값이 결과값으로 반환됩니다 {       if(!shift) return 0;       XL_UINT32 rem = dest[len - 1] >> (32 - shift);     dest[len - 1] = dest[len - 1] << shift;     for(XL_UINT32 i = len - 2; i != (XL_UINT32)(-1); --i)     {         dest[i   1] |= dest[i] >> (32-shift);         dest[i] = dest[i] << shift;     }     dest[0] |= fill;       return rem; }

직관적인 코드라 설명은 별도로 필요가 없을거같네요! 아주아주 중요한 주의사항은 shift값이 32보다 작아야한다는 겁니다. 어차피 정수가 32비트 단위로 저장되기 때문에 32비트 이상 시프트하는 건 정수 단위로 움직이는것 32미만의 비트 시프트로 분할해서 생각할 수 있기 때문이죠. 그래서 아예 32비트 이상 시프트하는것은 고려하지도 않았습니다.

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

void xl_shl(PXLARGE dest, XL_UINT32 bit) {     if(xl_isnan(dest)) return;     dest->exp  = bit >> 5; // 32의 배수로 시프트하는 경우는 지수부를 1씩 더해주기만 하면됩니다. // 32비트 미만 시프트하는 경우만 실제로 시프트 연산을 해주면 되겠죠     if(bit & 31)     {         _XL_PushBack(dest); //시프트 연산하면 혹시나 밀려나는 값이 생길수 있으니, 보관할수 있게 공간을 마련해줍니다.         if(dest->begin < dest->end)         {             i_shls(&dest->frac[dest->begin], dest->end - dest->begin, bit & 31, 0); // 배열 한 덩어리인 경우는 그냥 시프트하면 끝.         }         else         {             XL_UINT32 carry = i_shls(&dest->frac[dest->begin], dest->size - dest->begin, bit & 31, 0); // 두 덩어리로 쪼개져있는 경우는 높은쪽을 먼저 시프트하고, 밀려난 비트들을 carry에 담아두었다가             if(dest->end) i_shls(dest->frac, dest->end, bit & 31, carry); // 낮은쪽 시프트 연산시에  보충해줍니다.         }     }     _XL_Shrink(dest); }   void xl_shr(PXLARGE dest, XL_UINT32 bit) // 위와 모두 동일하고 동작 방향만 반대입니다. {     if(xl_isnan(dest)) return;     dest->exp -= bit >> 5;     if(bit & 31)     {         _XL_PushFront(dest);         if(dest->begin < dest->end)         {             i_shrs(&dest->frac[dest->begin], dest->end - dest->begin, bit & 31, 0);         }         else         {             XL_UINT32 carry = 0;             if(dest->end) carry = i_shrs(dest->frac, dest->end, bit & 31, 0);             i_shrs(&dest->frac[dest->begin], dest->size - dest->begin, bit & 31, carry);         }     }     _XL_Shrink(dest); }

자, 시프트 연산 구현을 마쳤으니, 이제 본격적으로 나머지를 구하는 나눗셈 연산을 구현해봅시다!

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

void xl_mod(PXLARGE dest, PXLARGE divisor, PXLARGE quot) {     if(xl_isnan(dest) || xl_isnan(divisor)) return xl_setnan(quot);     XL_UINT32 bsDivisor = i_bsu(*_XL_GetFracBack(divisor)); // 최고 자릿수에 들어가있는 정수의 최상위 비트를 구합니다. // 최고 자릿수의 모든 비트가 0이라는건 값이 0이라는 말이죠     if(bsDivisor == (XL_UINT32)-1)     {         xl_setnan(quot); // 과감하게 NaN을 설정하고 끝.         return;     }       XL_UINT32 destLen = _XL_GetLength(dest);     XL_UINT32 divLen = _XL_GetLength(divisor);     XL_UINT32 expBias = divisor->exp - dest->exp; // 나누기를 무한정 반복할수는 없으니 // 나누는 수의 지수부에 맞춰서 나누는 걸로 약속해요 // 그렇게 되면 몫이 자동적으로 정수만 나오게 됩니다.     if(dest->exp > divisor->exp) // 나눠지는 수가 나누는 수보다 지수부가 큰 경우 // 아래쪽으로 확장해서 미래를 대비합니다.     {         expBias = 0;         _XL_PushFrontN(dest, dest->exp - divisor->exp);         destLen = _XL_GetLength(dest);     }   // 나누는 수가 더 크면 몫은 0 // 나머지는 dest그대로가 됩니다.     if(dest->exp   destLen < divisor->exp   divLen)     {         xl_init(quot, 0);         return;     }     XL_INT32 bit = ((dest->exp   destLen) - (divisor->exp   divLen)) << 5;     bit  = i_bsu(*_XL_GetFracBack(dest)) - bsDivisor; // 나눠지는 수와 나누는 수의 최상위비트 자릿수의 차이를 구합니다.     if(bit < 0) // 나눠지는 수의 최상위 비트가 더 작다는건 나누는수가 더 크다는거고 // 역시나 몫이 0이되겠죠     {         xl_init(quot, 0);         return;     }   // 이제 나눗셈을 시작할거에요 // 나누는수를 왼쪽으로 시프트해서 나눠지는 수랑 최상위 비트를 같게만듭니다.     xl_shl(divisor, bit); // 몫은 0으로 초기화하고, 넉넉하게 여유공간을 잡아둡니다.     xl_init(quot, 0);     _XL_PushBackN(quot, bit >> 5);   // 최상위비트부터 차례로 1씩 내려오면서     for(XL_UINT32 i = bit;;--i)     { // 나누는수와 나눠지는 수를 비교합니다. // 나눠지는수가 더 크다는것은 // 그 자리에서 한번 뺄 수 있다는걸 의미하죠         if(_Abs_Cmp(dest, divisor) >= 0)         {             _Abs_Sub(dest, divisor); // 그럼 빼줘야지.             *_XL_GetFrac(quot, i >> 5) |= 1 << (i & 31); // 한 번 뺐다는 것은 그 자릿수의 몫이 1이 된다는거죠.         }         if(i == 0) break; // 탈출조건         xl_shr(divisor, 1); // 한 자리 내려가기 위해, 나누는수를 오른쪽으로 1비트 시프트합니다     } // 몫의 부호는 나누는수와 나눠지는 수의 부호의 곱으로 하고 // 나머지의 부호는 나눠지는 수를 그대로 따라가도록 합니다.     quot->sign = dest->sign ^ divisor->sign;     _XL_Shrink(quot); }

산수시간에 배웠던 나눗셈 방법과 본질적으로 동일합니다. 차이점이 있다면 이진법을 쓴다는거죠. 몫의 각 자리가 0 아니면 1이 되기때문에 나눠지는수와 나누는수를 비교해서 나눠지는수가 크다면 그 자리의 몫은 1, 아니면 0으로 간단하게 결정할 수 있습니다.

이방법은 사실 코드만 봐도 저렇게 긴데, 당연히 빠른 알고리즘은 아닙니다. 뺄셈이 O(n)시간에 수행된다고 했었는데, 그런 뺄셈을 최악의 경우 32*n번 해야하니, 전체 시간복잡도는 O(n^2)이 됩니다. 물론 곱셈도 O(n^2)이었는데, 나눗셈이 곱셈보다도 훨씬 느립니다.

속도를 더 빠르게 할 방법이 없을까요?

없다면 이 글을 안 썼겠죠. 분할정복(Divide & Conquer)의 방법을 사용해서 시간복잡도를 비약적으로 낮추는 방법이 있습니다. 여기까지 가기 전에 먼저 곱셈 알고리즘부터 개선하고 가는게 좋겠지요?

다음 글에서는 O(n^2)의 벽을 깨버린 곱셈 알고리즘을 소개하는 시간을 갖도록하겠습니다.