[C++] 15576번 큰 수 곱셈 (2)

https://www.acmicpc.net/problem/15576

---

 

26/01/28

 

 

FFT의 기본을 배울 수 있는 문제이다. FFT를 배운다면 가장 먼저 이 문제를 풀어보기를 추천한다.

---

 

문제 접근 방식:

 

 

$10$진법의 정의를 생각해보자.

예를 들어, $12345242_{10}$이라고 하면, 정의에 따라 $2\cdot 10^{0}+4\cdot 10^{1}+2\cdot 10^{2}+5\cdot 10^{3}+\dots+1\cdot 10^{7}$와 같다.

그리고 이는 아래와 같은 다항식의 $x$자리에 $10$을 대입한 것과 같다.
$$2+4x+2x^{2}+5x^{3}+\dots+x^{7}$$
FFT를 활용하여 두 다항식을 빠르게 곱할 수 있다는 사실을 알고 있다면, 이제 무엇을 할지 감이 올 것이다.

먼저 위의 예시처럼 각 자리별로 다항식의 계수를 만들어준다.

이후 FFT로 두 다항식을 빠르게 곱해주자. 이때 시간 복잡도는 $\mathcal{O}(N\log N)$이 소요되며, 여기서 $10$진법을 기준으로 삼는다면 다항식의 길이 $N$이 문제에서 주어지는 $30$만 자리가 된다.

물론, $10$진법이 아니라 $100$진법, $1000$진법 등으로 끊어줄 수도 있다. 그러면 다항식의 길이 $N$또한 줄어들 것이다.

FFT로 곱해준 다항식의 결과에 이제 $10$을 대입해주자. 당연히 결과로 나온 다항식의 계수가 $10$이하라는 보장은 없다.

즉, 예를 들어 다음과 같은 결과가 나올 수 있다.
$$2+12x+21x^{2}+\dots$$
우리가 원하는 것은 $10$진법으로 표현한 모습이기 때문에 이런 경우 자리올림을 구현해주면 된다. (물론 $100$, $1000$진법 같이 다른 진법으로 구현했다면 그에 맞춰서 자리 올림을 해주면 된다.)

위의 예시라면 $12x$에서 $1$만큼의 자리올림이 일어난다.

이를 통해 결과를 만들어주면 된다.

---

아래는 내가 위의 접근 방식과 같이 작성한 C++ 코드이다. 더보기를 누르면 확인할 수 있다.

// 15576번 큰 수 곱셈 (2)
// FFT
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;
#define fastio ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
#define endl '\n'

struct FFT {
    using ld = long double;
    static constexpr ld PI = acosl(-1.0L);
private :
    // HELPER struct
    struct cd {
        ld re, im;
        cd(ld r=0, ld i=0): re(r), im(i) {}
        cd operator + (const cd& o) const { return cd(re+o.re, im+o.im); }
        cd operator - (const cd& o) const { return cd(re-o.re, im-o.im); }
        cd operator * (const cd& o) const { return cd(re*o.re - im*o.im, re*o.im + im*o.re); }
        cd& operator += (const cd& o) { re+=o.re; im+=o.im; return *this; }
        cd& operator -= (const cd& o) { re-=o.re; im-=o.im; return *this; }
        cd& operator *= (const cd& o) { *this = (*this)*o; return *this; }
    };
public : 
    void iterative_FTT(vector<cd>& a, int invert){
        int n = (int)a.size();
        // 1) bit-reversal permutation
        for (int i = 1, j = 0; i < n; i++){
            int bit = n >> 1;
            while (j & bit){
                j ^= bit;
                bit >>= 1;
            }
            j ^= bit;
            if (i < j) swap(a[i], a[j]);
        }
        // 2) butterflies by length = 2, 4, 8, ...
        for (int len = 2; len <= n; len <<= 1){
            ld ang = 2 * PI / len * (invert ? -1 : 1);
            cd w_len(cosl(ang), sinl(ang)); // primitive len-th root of unity on unit circle
            for (int i = 0; i < n; i += len){
                cd w(1, 0);
                for (int j = 0; j < len / 2; j++){
                    cd u = a[i+j];
                    cd v = w * a[i+j+len/2];
                    a[i+j] = u+v;
                    a[i+j+len/2] = u-v;
                    w *= w_len;
                }
            }
        }
        // 3) Divide by n (multiply by inverse) for inverse transform
        if (invert){
            ld inv_n = 1.0L / n;
            for (int i = 0; i < n; ++i){
                a[i].re *= inv_n;
                a[i].im *= inv_n;
            }
        }
    }
    // Input : Coefficient vector {a_0, a_1, ...}, {b_0, b_1, ...}
    // Output : Convolution of two coefficient vector
    vector<int> convolution(const vector<int>& a, const vector<int>& b) {
        if (a.empty() || b.empty()) return {};
        int S_a = (int)a.size(), S_b = (int)b.size();
        // Make vector size easy to dnc(2^n).
        int n = 1;
        while (n < S_a + S_b - 1){
            n <<= 1;
        }
        vector<cd> fa(n), fb(n);
        for (int i = 0; i < S_a; i++) fa[i] = cd((ld)a[i], 0);
        for (int i = 0; i < S_b; i++) fb[i] = cd((ld)b[i], 0);
        // FTT
        iterative_FTT(fa, 0); iterative_FTT(fb, 0);
        // Pointwise product
        for (int i = 0; i < n; i++) fa[i] *= fb[i];
        // IFTT
        iterative_FTT(fa, 1);
        vector<int> res(S_a + S_b - 1);
        for (int i = 0; i < res.size(); i++) {
            res[i] = (int)llround(fa[i].re); // rounding to nearest integer
        }
        return res;
    }
};

int main(void){
    fastio

    vector<int> a, b;
    char x; int flag = 0;
    while (cin.get(x)){
        if (x == ' ') {
            flag = 1;
            continue;
        }
        if (x == '\n') break;
        if (flag) b.push_back((int)(x-'0'));
        else a.push_back((int)(x-'0'));
    }
    FFT fft;
    auto c = fft.convolution(a, b);
    vector<int> ans;
    int carry = 0;
    for (int i = c.size()-1; i > 0; --i){
        ans.push_back((c[i]+carry) % 10);
        carry = ((c[i]+carry) / 10);
    }
    ans.push_back(c[0] + carry);
    for (int i = ans.size()-1; i >= 0; --i){
        cout << ans[i];
    }
    return 0;
}