2.7 나눗셈(divide)
큰 수에 대한 나눗셈을 생각해 보자.
곱셈의 경우는 앞 장에서 설명한 바와 같이, m자리와 n자리의 두 수를 곱하면 (m+n)자리의 결과값이 나온다.
나눗셈에서는 (m+n)자리의 수를 n자리 정수로 나누면, m자리 이하의 몫이 나온다.
나눗셈에 대한 알고리즘도 덧셈,뺄셈,곱셈과 마찬가지로 손으로 공책에 써가며 계산하는 손계산 방법을 생각해보자.
먼저, 익숙한 십진수에 대한 예를 가지고 나눗셈에 대한 기본 원리를 살펴 보겠다.
| 7 | 1 | 2 | 1 | 0 | 0 | 0 |
| 2 | ||||||
| 7 | 1 | 2 | 1 | 0 | 0 | 0 |
| 1 | 4 | 2 | ||||
| 6 | 8 | |||||
| 2 | 9 | |||||
| 7 | 1 | 2 | 1 | 0 | 0 | 0 |
| 1 | 4 | 2 | ||||
| 6 | 8 | 0 | ||||
| 6 | 3 | 9 | ||||
| 4 | 1 | |||||
| 2 | 9 | 5 | ||||
| 7 | 1 | 2 | 1 | 0 | 0 | 0 |
| 1 | 4 | 2 | ||||
| 6 | 8 | 0 | ||||
| 6 | 3 | 9 | ||||
| 4 | 1 | 0 | ||||
| 3 | 5 | 5 | ||||
| 6 | 5 |
위는 손계산으로 십진수 21000을 71로 나누는 과정이다. 나눈 결과 몫은 295가 나오고 나머지는 65이다.
피제수를 u, 제수를 v, 몫을 q, 나머지를 r이라 했을 때, 위 예에서 u=21000, v=71, q=295, r=65 이다.
일반화하면 u = (um+n-1 um+n-2 ... u1 u0)b이고 v=(vn-1 vn-2 ... v1 v0)b일 때, 몫 q=(qm-1 qm-2 ... q1 q0)이고 r은 (0<= r < q)이고, u = vq + r로 표현된다. 알고리즘은 u와 v가 주어졌을 때, 이와 같은 조건을 만족시키는 q와 r을 구하는 것이된다.
그렇다면 위의 손계산을 어떻게 알고리즘화 할 것인가? 덧셈, 뺄셈, 곱셈과 달리 약간의 어림짐작이 필요하다. 즉, 위 예에서 몫의 첫번째 수가 '2'가 됨을 첫번째 두 자리 21 < 71이므로 210까지 확장해서 피제수를 고려하고, 이 210/71을 만족하는 것이 '2'라는 것을 머릿속으로 알아내는 과정이 손계산에서는 일어난다. 이것을 일반화해서 알고리즘화해야 한다.
위의 나눗셈과정을 자세히 보면 5자리의 수 21000를 2자리수 71로 나누는데, 각 과정은 먼저 210을 71로 나눠서 2라는 몫을 알아내고, 그 다음 210에서 (71 x 2)한 142을 빼서 나온 68을, 10을 곱한 후 그 다음 자릿수를 더한 값인 680로 놓고, 이 680에 대해서 다시 71로 나누는 과정을 반복함을 알 수 있다.
즉, 각 단계는 (n+1)자리 피제수 u를 n자리의 제수 v로 나누는 것으로 생각할 수 있고, 이 때 0 <= u/v < b가 된다. (십진수에서 b = 10)
그럼 (un un-1 ... u1 u0)를 (vn-1 vn-2 ... v1 v0)로 나누는 qm-1의 값을 어떻게 구할까?
여기서 사용할 수 있는 좋은 정리가 있다.
q' = ⌊(unb+un-1)/vn-1⌋일 때, 만약 vn-1 >= ⌊b/2⌋이면 q' - 2 <= q <= q' 이다.
위 예제에서 q' = (2x10+1)/7 = 3이고 vn-1 = 7 >= 5 따라서, q는 1,2,3중 하나가 된다.
먼저 3으로 가정하면 210 - (71 x 3) = 210 - 213 = -3 < 0 으로 부적합.
1을 감소시킨 2로 가정하면 210 - (71 x 2) = 210 - 142 = 68 > 0 으로 적합
이러한 과정으로 손계산을 해보면 다음과 같겠다.
vn-1 >= ⌊b/2⌋ 이기만 하다면 많아도 3번 이내에 q값들을 찾아낼 수 있게 되었다.
그러나, 문제는 vn-1 >= ⌊b/2⌋이어야 한다. 이것은 강제로 vn-1을 ⌊b/2⌋보다 크게 만들면 해결된다. 즉, vn-1에 어떤 값 p를 곱해서 ⌊b/2⌋보다 크게만들고, 같은 값 p를 u에다 곱해주는 것이다. 그렇게 해도 q는 같은 값을 가지게 된다. 왜냐하면 q = u/v = up/vp 이기 때문
이제 위와같은 원리를 이용해서 나눗셈에 대한 알고리즘을 생각해 보면 다음과 같다.
알고리즘 (음이 아닌 정수들의 나눗셈)
음의 아닌 정수 u = (um+n-1 um+n-2 ... u1 u0)b이고 v=(vn-1 vn-2 ... v1 v0)b가 주어졌으며 vn-1 ≠ 0 이고, n > 1, b는 기수일 때, 이 알고리즘은 몫 ⌊u/v⌋ = q = (qm-1 qm-2 ... q1 q0)b와 나머지 u mod v = r = (rn-1 rn-2 ... r1 r0)b를 구한다.
S1. 정규화
vn-1 >= ⌊b/2⌋이 되도록하는 p를 구하고, u와 v에 곱한다.
S2. 루프(j)
j ← m
S3. q'을 계산
S4. 곱하고 뺀다.
S5 j에 대한 루프
if(j >= 0) goto S2
END
S1. 정규화
vn-1 >= ⌊b/2⌋이 되도록하는 p를 구하고, u와 v에 곱한다.
p를 구하는 쉬운 방법은 vn-1이 ⌊b/2⌋될때까지 2를 곱하는 것이다. 즉, vn-1 x 2d >= ⌊b/2⌋되는 d를 구하면 되는 것이고, ⌊b/2⌋보다 커질 때까지 vn-1을 왼쪽으로 shift시키는 것으로 쉽게 구현될 수 있다.
S2. 루프(j)
j의 초기값은 m이고, 루프에서는 (uj+n ... uj+1 uj)b를 (vn-1 ... v1 v0)b로 나누어서 한 자리 몫 qj를 얻는다.
S3. q'을 계산
q' ← ⌊(uj+nb + uj+n-1) / vn-1⌋로 놓고, 나머지 r'을 (uj+nb + uj+n-1) mod vn-1로 둔다.
이제 q'vn-2 > br' + uj+n-2를 판정한다. 만일 참이면 q'을 1감소시킨다.
S4. 곱하고 뺀다
(uj+n ... uj+1 uj)b를 (uj+n ... uj+1 uj)b - q'(0vn-1 ... v1 v0)b로 대체한다.
S5. j에 대한 루프
j >= 0에 대해서 루프 진행
n=1일 때의 나눗셈
제수 v가 n=1인 한 자리 수인 경우는 나눗셈이 간단해진다. 이 때는 q'을 유추할 필요 없이 바로 q의 계산이 가능하다.
십진수로 예를 들면 (2345 / 3)같은 경우로, 처음 ⌊23/3⌋에 해당하는 7을 구하고, 3x7=21을 23에서 뺀값인 2와 v의 그 다음 자릿수인 4를 결합해서 24로 만들고, 다시 ⌊24/3⌋=8을 구하는 식으로 진행이 된다.
| 3 | 2 | 3 | 4 | 5 | |
| 7 | |||||
| 3 | 2 | 3 | 4 | 5 | |
| 2 | 1 | ||||
| 2 | |||||
| 7 | 8 | ||||
| 3 | 2 | 3 | 4 | 5 | |
| 2 | 1 | ||||
| 2 | 4 | ||||
| 2 | 4 | ||||
| 7 | 8 | 1 | |||
| 3 | 2 | 3 | 4 | 5 | |
| 2 | 1 | ||||
| 2 | 4 | ||||
| 2 | 4 | ||||
| 5 | |||||
| 3 | |||||
| 2 |
코딩
나눗셈을 구현하는 함수를 3개로 구분하도록 하자.
- divide_int: 나누는 제수의 길이가 1 즉, INT형일때 나눗셈을 계산하는 함수
- divide_positive: 양의 정수에 대한 계산을 담당. 여기서 제수가 1인 경우는 divide_int이용 계산
- divide: 양수, 음수 모두 사용할 수 있는 함수. 절대값으로 수를 바꾼후 divide_positive로 계산 후, 부호 붙여주는 구조
divide_int
/************************************************************************** divide BIGINT number by INT. A/b = Qb + R Input: @A = dividend. should be positive BIGINT value @b = divisor. should be positive INT value Output: @Q = quotient @R = Remainder Returns: OK : no error E_DivideByZero: in case of divisor(B) is zero Comment: ***************************************************************************/ int divide_int(BIGINT A, INT b, BIGINT Q, BIGINT R){ register INT *rPtr; INT *qPtr; INT r[(MAX_INT_LEN << 1) + 2]; INT tmp; DINT rhat; unsigned int d = 0; int shift; copy_big(r,A); //divide by zero -> Error if(b == 0) return E_DivideByZero; shift = sizeof(INT) * 8; tmp = 0; rPtr = r + SIZE(r); qPtr = Q + SIZE(r); for (; rPtr >= (r+1); rPtr--, qPtr--) { *qPtr = (INT)((rhat = ((((DINT)tmp) << shift) + (DINT)*rPtr)) / b); tmp = (INT)(rhat - (DINT)b * (DINT)*qPtr); } Q[0] = SIZE(r); REMOVE_ZERO(Q); R[0]=1; R[1]=tmp; return OK; }
divide_positive
제수와 피제수가 0 혹은 1 등 특수한 경우의 수를 고려하여 처리한다.
A/B에 대해 몫이 Q, 나머지가 R인 경우
Case1. zero로 나누는 경우는 에러 리턴. if(B==0) Error
Case2. 0을 나누는 경우는(0/B), 몫과 나머지가 zero
Case3. A와 B가 같은 경우, Q=1 R=0
Case4. A < B의 경우 Q=0 R=A
Case 5. B가 한 자리 수인 경우, divide_int함수 이용
Case 6. 그 외의 경우 divide_positive알고리즘 이용
/************************************************************************** divide BIGINT number by BIGINT. A/B = qB + r Input: @A = dividend. should be positive BIGINT value @B = divisor. should be positive BIGINT value Output: @Q = quotient @R = Remainder Returns: OK : no error E_DivideByZero: in case of divisor(B) is zero Comment: ***************************************************************************/ int divide_positive(BIGINT A, BIGINT B, BIGINT Q, BIGINT R){ register INT *rPtr, *bPtr; INT *qPtr, *endB, *startR, *endR; BIGINT b; INT r[(MAX_INT_LEN << 1) + 2]; DINT HALF_BASE = (BASE >> 1), shiftVal; INT b_n1,b_n2; INT r_i0, r_i1, r_i2; INT qhat; DINT rhat,right,left,t,carry; unsigned int d = 0; int i; int shift; copy_big(r,A); copy_big(b,B); //Case1. divide by zero -> Error if(*b == 0) return E_DivideByZero; //Case2. dividend(=A) is zero -> q=0 r=0 // divide함수에서 처리 //Case3. A == B -> q=1 r=0 i = compare_big(r,b); if(i == 0){ Q[0] = 1; Q[1] = 1; R[0] = 0; return OK; }else if(i == -1){ //Case4. A < B -> q=0 r=A Q[0] = 0; copy_big(R, A); return OK; } //Case5. B's length = 1 if(SIZE(b) == 1){ return divide_int(A,B[1],Q,R); } //Case6. B's length > 1 //S1. 정규화 // make b_n1 sufficently large number such as > HALF_BASE shift = sizeof(INT) * 8; endB = b + SIZE(b); b_n1 = *endB; d=0; while(b_n1 < HALF_BASE){ d++; b_n1 <<= 1; } shiftVal = (int)((sizeof(INT) * 8) - d); if(d>0){ b_n1 += (INT)((DINT)(*(endB -1)) >> shiftVal); //b[n-1] b[n-2]를 d만큼 shift한 값 if(SIZE(b) > 2){ b_n2 = (INT)((*(endB -1) << d) + ((DINT)(*(endB -2)) >> shiftVal)) ; }else{ b_n2 = (INT)(*(endB -1) << d); } }else{ b_n2 = (INT)(*(endB - 1)); } endR = r + SIZE(r)+1; //r[2+ 2*MAX_INT_LEN] startR = (r + SIZE(r)) - SIZE(b) + 1; // [0] [1=startR] ...[endR]=0 [endR+1]... [1+2MAX] *endR = 0; qPtr = Q + SIZE(r) - SIZE(b) + 1; //S2. 루프 while(startR >= (r+1)){ //S3. q'계산 r_i0 = (INT)((*endR << d) + ((DINT)(*(endR-1)) >> shiftVal)); r_i1 = (INT)((*(endR-1) << d) + ((DINT)(*(endR-2)) >> shiftVal)); if((endR-3) > r) r_i2 = (INT)((*(endR-2) << d) + ((DINT)(*(endR-3)) >> shiftVal)); else r_i2 = (INT)(*(endR-2) << d); // if(r_i0 != b_n1){ qhat = (INT)((rhat = ((DINT)r_i0 << shift) + (DINT)r_i1) / b_n1); right = ((rhat = (rhat - (DINT)b_n1 * qhat)) << shift) + r_i2; if((left = (DINT)b_n2 * qhat) > right){ //예상했던 q가 너무 크면 qhat--; if ((rhat + b_n1) < BASE) { if ((left - b_n2) > (right + ((DINT)b_n1 << shift))) { qhat--; } } } }else{ qhat = MAX_INT; //ffff or ffffffff right = ((DINT)(rhat = (DINT)b_n1 + (DINT)r_i1) << shift) + r_i2; if (rhat < BASE) { if ((left = (DINT)b_n2 * qhat) > right) { qhat--; if ((rhat + b_n1) < BASE) { if ((left - b_n2) > (right + ((DINT)b_n1 << shift))) { qhat--; } } } } } //S4. 곱하고 뺀다. t = BASE; carry = 0; bPtr = b+1; rPtr = startR; for (; bPtr <= endB; bPtr++, rPtr++) { if (t >= BASE){ *rPtr = (INT)(t = ((DINT)(*rPtr) + BASE - (DINT)(INT)(carry = (DINT)(*bPtr) * qhat + (DINT)(INT)(carry >> shift)))); }else{ *rPtr = (INT)(t = ((DINT)(*rPtr) + MAX_INT - (DINT)(INT)(carry = (DINT)(*bPtr) * qhat + (DINT)(INT)(carry >> shift)))); } } if (t >= BASE){ *rPtr = (INT)(t = ((DINT)(*rPtr) + BASE - (DINT)(INT)(carry >> shift))); }else{ *rPtr = (INT)(t = ((DINT)(*rPtr) + MAX_INT -(DINT)(INT)(carry >> shift))); } // *qPtr = qhat; if (t < BASE) { carry = 0; bPtr = b+1; rPtr = startR; for (; bPtr <= endB; bPtr++, rPtr++){ *rPtr = (INT)(carry = ((DINT)(*rPtr) + (DINT)(*bPtr) + (DINT)(INT)(carry >> shift))); } *rPtr += (INT)(carry >> shift); (*qPtr)--; } endR--; startR--; qPtr--; } Q[0] = SIZE(r) - SIZE(b) + 1; REMOVE_ZERO(Q); r[0] = SIZE(b); copy_big(R,r); return OK; }
divide
양수, 음수 등 모든 경우의 수에 대해 나눗셈을 수행할 수 있는 함수이다.
A/B에 대해서 몫이 Q, 나머지가 R이라 할 때, 각 경우의 수를 보면,
1. 0으로 나누는 경우(B=0), 에러
2. 0을 나누는 경우, Q=0 R=0
3. A와 B의 부호가 같은 경우 Q > 0, R >= 0
4. A와 B의 부호가 다른 경우 Q < 0, R >= 0
/************************************************************************** divide BIGINT number by BIGINT. A/B = qB + r Input: @A = dividend. should be positive BIGINT value @B = divisor. should be positive BIGINT value Output: @Q = quotient @R = Remainder Returns: OK : no error E_DivideByZero: in case of divisor(B) is zero Comment: ***************************************************************************/ int divide(BIGINT A, BIGINT B, BIGINT Q, BIGINT R){ int resultSign, ret; BIGINT absA, absB; //0. if divide by zero, return E_DivideByZero if(*B == 0){ *Q = 0; *R = 0; return E_DivideByZero; } //1. divide zero then Q=0 R=0 if(*A == 0){ *Q = 0; *R = 0; return OK; } //2. 절대값(=양수)로 변환 후 계산 resultSign = isPositive(A) * isPositive(B); copy_big(absA,A); abs(absA); copy_big(absB,B); abs(absB); ret = divide_positive(absA,absB,Q,R); //3. 부호 결정 if(ret == OK){ if(resultSign == -1) //Q만 음수로 바꿔줌. R은 이미 양수 return negative(Q); else return OK; }else{ return ret; } }
테스트
다음과 같은 테스트 케이스를 생각할 수 있다.
1. divide by zero -> error
2. 0 / b = 0
3. a / 1 = a
4. a / -1 = -a
5. -a / 1 = -a
6. -a / -1 = a
7. a / b : q=0, r=a (if a < b)
8. a / b (if b's length = 1, r=0)
9. a / b (if b's length = 1, r!=0)
10. a / b (if b's length > 1)
위 케이스에 따라 구현된 테스트 함수는 다음과 같다.
int test_divide(){ BIGINT a,b,q,r; BIGINT expected_q, expected_r; int result=-1; int ret=-1; int i; //1. divide by zero -> error hexstr2bigint("112233445566",a); b[0]=0; ret = divide(a,b,q,r); if(ret != E_DivideByZero) return -1; //2. 0 / b = 0 a[0] = 0; hexstr2bigint("112233445566",b); ret = divide(a,b,q,r); if(ret != OK) return -1; expected_q[0]=0; result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0]=0; result = compare_big(r,expected_r); if(result != 0) return -1; //3. a / 1 = a hexstr2bigint("112233445566",a); hexstr2bigint("01",b); ret = divide(a,b,q,r); if(ret != OK) return -1; hexstr2bigint("112233445566",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0]=0; result = compare_big(r,expected_r); if(result != 0) return -1; //4. a / -1 = -a hexstr2bigint("112233445566",a); hexstr2bigint("-1",b); ret = divide(a,b,q,r); if(ret != OK) return -1; hexstr2bigint("-112233445566",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0]=0; result = compare_big(r,expected_r); if(result != 0) return -1; //5. -a / 1 = -a hexstr2bigint("-112233445566",a); hexstr2bigint("01",b); ret = divide(a,b,q,r); if(ret != OK) return -1; hexstr2bigint("-112233445566",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0]=0; result = compare_big(r,expected_r); if(result != 0) return -1; //6. -a / -1 = a hexstr2bigint("-112233445566",a); hexstr2bigint("-1",b); ret = divide(a,b,q,r); if(ret != OK) return -1; hexstr2bigint("112233445566",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0]=0; result = compare_big(r,expected_r); if(result != 0) return -1; //7. a / b : q=0, r=a (if a < b) hexstr2bigint("0100",a); hexstr2bigint("112233445566",b); ret = divide(a,b,q,r); ret = divide(a,b,q,r); if(ret != OK) return -1; expected_q[0] = 0; result = compare_big(q,expected_q); if(result != 0) return -1; hexstr2bigint("0100",expected_r); result = compare_big(r,expected_r); if(result != 0) return -1; //8. a / b (if b's length = 1, r=0) hexstr2bigint("01258BF25308",a); hexstr2bigint("1122",b); ret = divide(a,b,q,r); hexstr2bigint("11223344",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; expected_r[0] = 0; result = compare_big(r,expected_r); if(result != 0) return -1; //9. a / b (if b's length = 1, r!=0) hexstr2bigint("112233445566",a); hexstr2bigint("0100",b); ret = divide(a,b,q,r); hexstr2bigint("1122334455",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; hexstr2bigint("66",expected_r); result = compare_big(r,expected_r); if(result != 0) return -1; //10. a / b (if b's length > 1) hexstr2bigint("11223344556677",a); hexstr2bigint("1122334455",b); ret = divide(a,b,q,r); hexstr2bigint("10000",expected_q); result = compare_big(q,expected_q); if(result != 0) return -1; hexstr2bigint("6677",expected_r); result = compare_big(r,expected_r); if(result != 0) return -1; return OK; }
프로젝트 파일
007.Divide.zip
'암호화프로그래밍 > (Old)C-BIGINT' 카테고리의 다른 글
| 015. 최대공약수(gcd) (0) | 2013.08.09 |
|---|---|
| 014. 나머지(mod) (0) | 2013.08.08 |
| 012. 제곱(square) (0) | 2013.07.31 |
| 011. 곱셈(multiply) (0) | 2013.07.30 |
| 010. 감소(decrese) (0) | 2013.07.30 |