[진짜 BigFloat 구현하기] 5. 카라슈바 알고리즘으로 곱셈 성능 향상

곱셈은 앞서(http://bab2min.tistory.com/405)에서 이미 구현했습니다. 다 좋은데 시간복잡도가 O(n^2)이라 숫자가 커지면 커질수록 너무 느려진다는 단점이 있죠. 하지만 이런 문제를 겪을줄 알았는지 이미 일찍이 수학자들이 여러 가지 빠른 곱셈 알고리즘을 발명해 놓았습니다. 여기에서는 곱셈 알고리즘 중 두번째로 구현이 간단한(그 댓가로 뒤에서 두번째로 빠른) 카라슈바 알고리즘(Karatsuba Algorithm)을 구현해보도록 할것입니다.

카라슈바 알고리즘의 개요는 다음과 같습니다. 2자리수 간의 곱을 계산하려면 일반적인 곱셈법으로는 한 자릿수 곱셈이 총 4번이 필요한데, 카라슈바는 한 자릿수 곱셈을 3번만 하고, 덧셈/뺄셈으로 나머지 부분의 결과를 구해내는 아이디어를 떠올립니다.

(편의상 여러 자리 수를 [a, b, c]와 같은 형태로 표시하도록 하겠습니다.)

[a, b] * [c, d] 를 계산하면 [ac, ad bc, bd]의 결과가 나와야 하는거지요.

이를 위해서는 ac, ad, bc, bd를 구해야하므로 총 4번의 곱셈이 필요합니다.

하지만 ad bc = (a b)(c d) - ac - bd 라는것에 주목!

ac, bd는 원래 구해야하는 값이니깐 기억해두었다가, (a b)(c d)를 계산한뒤 빼주면

ad bc를 구할수 있겠죠. 이렇게 되면 곱셈이 3번으로 줄어들게 됩니다. 그 대신 약간의 덧셈 뺄셈이 들어가게 되겠죠.

이 방법을 사용하면 두 2N자리의 곱셈은 3번의 N자리 곱셈과 약간의 덧뺄셈으로 바뀔수 있겠죠. 마찬가지로 N자리끼리의 곱셈은 N/2자리끼리의 곱셈 3번으로 처리할수 있을 것이구요. 이를 재귀적으로 활용하면 원래는 N^2번의 곱셈이 필요했던 연산을 N^(log3)번의 곱셈으로 계산할수 있게 됩니다!

어떤 알고리즘인지 이해가 되셨다면, 바로 구현으로 들어가봅시다.

/* 앞서 구현했었던, 곱셈 알고리즘이에요. 혼란을 피하기 위해 이름을 _XL_MulN으로 바꾸었습니다. 크기가 작은 수의 경우는 그냥 이 알고리즘으로 계산하는게 가장 빠릅니다. */ void _XL_MulN(PXLARGE dest, PCXLARGE a, PCXLARGE b) { XL_UINT32 aLen = _XL_GetLength(a); XL_UINT32 bLen = _XL_GetLength(b); xl_init(dest, 0); _XL_PushBackN(dest, aLen bLen - 1); 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)); } } _XL_Shrink(dest); dest->exp = a->exp b->exp; }

예전에 구현했던 곱셈 연산함수입니다.

/* 해당 XLarge의 하위 N단위를 버리고 상위부분만 가져오는 함수입니다. */ inline void _XL_FracUpper(PXLARGE dest, PCXLARGE src, XL_UINT32 n) { dest->frac = src->frac; dest->size = src->size; dest->exp = 0; dest->sign = 0; dest->begin = (src->begin n) % src->size; dest->end = src->end; } /* 해당 XLarge의 하위 N단위만 가져오는 함수입니다. */ inline void _XL_FracLower(PXLARGE dest, PCXLARGE src, XL_UINT32 n) { dest->frac = src->frac; dest->size = src->size; dest->exp = 0; dest->sign = 0; dest->begin = src->begin; dest->end = (src->begin n) % src->size; }

이 두 함수를 구현한 것은 긴 정수배열을 반으로 잘라서 곱셈에 사용해야하기 때문이죠. 자주 사용하게 될테니, 미리 구현해두었습니다.

/* 카라슈바 알고리즘으로 곱해볼까요 */ void _XL_MulK(PXLARGE dest, PCXLARGE a, PCXLARGE b) { XL_UINT32 aLen = _XL_GetLength(a); XL_UINT32 bLen = _XL_GetLength(b); // a나 b의 길이가 threshold보다 작으면 일반 곱셈으로 넘겨버립니다. if(MIN(aLen, bLen) < K_THRESHOLD) { return _XL_MulN(dest, a, b); } /* 총 3가지 경우로 나눠볼수가 있습니다. 1. a가 b보다 2배 이상 길 때 2. b가 a보다 2배 이상 길 때 3. 나머지 경우 3의 경우는 적당한 크기로 a, b를 잘라서 카라슈바 알고리즘을 적용하면 됩니다. 하지만 1, 2번의 경우는 긴 쪽을 반으로 잘라도 짧은쪽과 길이가 비슷해서 카라슈바 알고리즘을 적용해봤자 이득이 없으니, 적당히 넘겨줍니다. */ if((aLen 1) / 2 4 >= bLen) // 1번 경우 { XL_UINT32 half = (aLen 1) / 2; XLarge a0, a1; _XL_FracLower(&a0, a, half); _XL_FracUpper(&a1, a, half); // a를 절반으로 잘라서 상위 하위를 나눕니다. PXLARGE z1 = xl_create(aLen 1); // 넉넉한 크기의 임시변수 생성 _XL_MulK(dest, &a0, b); // dest = a0 * b _XL_MulK(z1, &a1, b); // z1 = a1 * b z1->exp = half; // z1은 한 단위만큼 높은 자릿수여야하니 // half만큼 지수부를 더해줍니다. xl_add(dest, z1); // dest = z1 ; 즉 dest = (a0 a1 * (2^half)) * b xl_delete(z1); // 임시변수 삭제 dest->exp = a->exp; // a의 원래 지수부만큼을 더 더해줘 곱셈을 마무리합니다. } else if(aLen <= (bLen 1) / 2 4) // 2번 경우. 1번과는 a랑 b랑 바뀐거 빼고는 다 동일. { XL_UINT32 half = (bLen 1) / 2; XLarge b0, b1; _XL_FracLower(&b0, b, half); _XL_FracUpper(&b1, b, half); PXLARGE z1 = xl_create(bLen 1); _XL_MulK(dest, a, &b0); _XL_MulK(z1, a, &b1); z1->exp = half; xl_add(dest, z1); xl_delete(z1); dest->exp = b->exp; } else // 3번 경우. { XL_UINT32 half = (MAX(aLen, bLen) 1) / 2; // 둘 중 긴 쪽에 맞춰서 절반 크기를 구합니다. XLarge a0, a1, b0, b1; _XL_FracLower(&a0, a, half); _XL_FracLower(&b0, b, half); _XL_FracUpper(&a1, a, half); _XL_FracUpper(&b1, b, half); // a, b를 상위부분, 하위부분으로 쪼개구요. PXLARGE z0 = xl_create(2*half); PXLARGE z1 = xl_create(2*half); PXLARGE z2 = xl_create(2*half); // 넉넉한 크기로 임시변수 3개를 만듭니다. xl_move(z0, &a0); xl_add(z0, &a1); // z0 = a0 a1 xl_move(z2, &b0); xl_add(z2, &b1); // z2 = b0 b1 _XL_MulK(z1, z0, z2); // z1 = (a0 a1) * (b0 b1) _XL_MulK(z0, &a0, &b0); // z0 = a0 * b0 _XL_MulK(z2, &a1, &b1); // z2 = a1 * b1 xl_sub(z1, z0); xl_sub(z1, z2); // z1 = (a0 a1) * (b0 b1) - z0 - z2 z2->exp = 2 * half; z1->exp = half; xl_move(dest, z0); xl_add(dest, z1); xl_add(dest, z2); xl_delete(z2); xl_delete(z1); xl_delete(z0); dest->exp = a->exp b->exp; // a, b의 원래 지수만큼을 더 더해줘서 곱셈을 마무리합니다. } _XL_Shrink(dest); }

실제 카라슈바 알고리즘 구현!

void xl_mul(PXLARGE dest, PCXLARGE a, PCXLARGE b) { if(xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest); _XL_MulK(dest, a, b); dest->sign = a->sign ^ b->sign; }

그리고 마무리로 xl_mul 함수를 이렇게 고쳤습니다. 이젠 xl_mul 함수를 호출하면 크기에 따라 적절하게 일반 알고리즘 혹은 카라슈바 알고리즘이 호출되겠죠?

집컴퓨터에서 테스트해보니 K_THRESHOLD는 대략 40~48 정도로 잡는게 적당할듯합니다. 그 정도 크기 이상부터 카라슈바 알고리즘이 더 빨라지게 됩니다. (구현이 썩 좋지는 않은가 봅니다. CPU에 따라 다르지만, 보통 2^320 정도부터는 카라슈바가 더 빠르다고 하는데...)