카라추바 곱셈 알고리즘

알고리즘 문제해결전략 183p에 소개된 카라추바 곰셈 알고리즘입니다.

 

2N자리수인 a와 b를 다음과 같이 정의한다고 하자.

 

$$a = a_{1} \times 10^{N}+a_{0}$$

$$ b = b_{1} \times 10^{N}+b_{0}$$

 

a x b 는 다음과 같이 표현할 수 있다.

 

$$a \times b \\ = ( a_{1}\times 10^{N} +a_{0}) \times (b_{1} \times 10^{N} + b_{0}) \\ = a_{1} \times b_{1} \times 10^{2N} + (a_{1} \times b_{0} + a_{0} \times b_{1})\times10^{N}+a_{0}\times b_{0}$$

 

그렇다면 계산을 위해 다음과 같이 곱셈을 기준으로 4가지 조각으로 나눌 수 있다.

$$a_{1} \times b_{1} $$

$$ a_{1} \times b_{0} $$

$$ a_{0} \times b_{1} $$

$$ a_{0} \times b_{0} $$

 

분할정복으로 4가지로 나눠 재귀적으로 계산한다면 시간복잡도는 다음과 같다.

$$O^{\log_{2}4} = O^{2}$$

 

카라추바 곱셈 알고리즘은 이 4조각을 3조각으로 바꾸어 계산하는 것이다.

(곱셈을 4개로 나누어져 있는 것을 3개로 나누어 계산하도록 한다.)

 

$$ z_{0} = a_{0} \times b_{0} $$

$$ z_{2} = a_{1} \times b_{1} $$

$$ z_{1} = (a_{0} + a_{1}) \times (b_{0} + b_{1}) - z_{0} - z_{1} $$

$$ ( z_{1} = a_{0} \times b_{1} + a_{1} \times b_{0} ) $$

$$ \therefore a \times b = z_{0} \times 10^{2N} + z_{1} \times 10^{N} + z_{2} $$

 

위와 같이 계산하면 곱셈을 3번만 사용하므로 최종 시간복잡도는 다음과 같습니다.

$$O^{\log_{2}3}$$

 

 

#include <bits/stdc++.h>
using namespace std;
#define swap(a,b) (a)^=(b)^=(a)^=(b)
#define endl '\n'
typedef long long lld;
 
void addto(vector<int>& a, vector<int>& b, int k)
{
    a.resize(max(a.size()+1, b.size()+k+1), 0);
    int carry = 0;
    for(int i=0;i<b.size();i++) a[i+k] += b[i];
    for(int i=0;i<a.size();i++)
    {
        a[i] += carry;
        carry = a[i]/10;
        a[i] %= 10;   
    }
    while(a.size()>=1 && a.back()==0) a.pop_back();
}
 
void subfrom(vector<int>& a, vector<int>& b)
{
    int bor = 0;
    for(int i=0;i<a.size();i++)
    {
        if(i<b.size()) a[i] -= b[i]+bor;
        else a[i] -= bor;
        if(a[i]<0)
        {
            a[i] += 10;
            bor = 1;
        }
        else bor = 0;
    }
}
 
vector<int> multiple(vector<int>& a, vector<int>& b)
{
    // a.size() >= b.size()
    vector<int> res(a.size()+b.size()+2, 0);
    for(int i=0;i<b.size();i++)
    {
        for(int j=0;j<a.size();j++)
        {
            res[i+j] += a[j]*b[i];
        }
    }
    int carry = 0;
    for(int i=0;i<res.size();i++)
    {
        res[i] += carry;
        carry = res[i]/10;
        res[i] %= 10;
    }
    while(res.back()==0) res.pop_back();
    return res;
}
 
 
vector<int> karatsuba(vector<int>& a, vector<int>& b)
{
    int alen = a.size(), blen = b.size();
    if(alen<blen) return karatsuba(b, a);
    if(alen==0 || blen==0) return vector<int>();
    else if(alen<=230)
        return multiple(a, b);
 
    int half = alen/2;
    vector<int> a0(a.begin(), a.begin()+half);
    vector<int> a1(a.begin()+half, a.end());
    vector<int> b0(b.begin(), b.begin()+min((int)b.size(), half));
    vector<int> b1(b.begin()+min((int)b.size(), half), b.end());
 
    vector<int> z2 = karatsuba(a1, b1);
    vector<int> z0 = karatsuba(a0, b0);
    addto(a0, a1, 0);
    addto(b0, b1, 0);
    vector<int> z1 = karatsuba(a0, b0);
    subfrom(z1, z0);
    subfrom(z1, z2);
 
    vector<int> ret;
    addto(ret, z0, 0);
    addto(ret, z1, half);
    addto(ret, z2, half+half);
    return ret;
 
}
 
int main()
{
    ios_base::sync_with_stdio(NULL);
    cin.tie(NULL);
    cout.tie(NULL);
    time_t start = clock();
    string ai, bi;
    cin>>ai>>bi;
    vector<int> a(ai.length(), 0), b(bi.length(), 0);
    for(int i=0;i<ai.length();i++)
        a[i] = ai[ai.length()-i-1]-'0';
    for(int i=0;i<bi.length();i++)
        b[i] = bi[bi.length()-i-1]-'0';
    vector<int> ans = karatsuba(a, b);
    for(int i=ans.size()-1;i>=0;i--)
        cout<<ans[i];
    cout<<endl;
    cout<<((double)clock()-start)/1000;
    return 0;
}