解題區►解題區 - Zerojudge

a577. 高精度乘法

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution

Problem

這題十分直接，就是要你計算兩個很大的數字的乘積。

Sample Input

1
2
3
4
5
0
5
11
12

Sample Output

1
2
0
132

Solution

快速傅立葉 FFT 對我來說也是一知半解，但是我能知道他計算旋積可以在 O(n log n) 完成 n 個答案結果。

旋積計算就是對於維度為 n 的兩個向量，相互作內積，其中一個向量不斷地作 rotate shift right 的操作，因此會有 n 個內積結果。

如果要套用在這一題，相當於操作多項式乘法的計算，舉個例子來說

123 x 456

相當於下面兩個向量作旋積

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
(3, 2, 1, 0, 0, 0, 0, 0)
(0, 0, 0, 0, 0, 4, 5, 6)
dot = 0
(3, 2, 1, 0, 0, 0, 0, 0)
(0, 0, 0, 0, 4, 5, 6, 0)
dot = 0
(3, 2, 1, 0, 0, 0, 0, 0)
(0, 0, 0, 4, 5, 6, 0, 0)
dot = 0
(3, 2, 1, 0, 0, 0, 0, 0)
(0, 0, 4, 5, 6, 0, 0, 0)
dot = 4
(3, 2, 1, 0, 0, 0, 0, 0)
(0, 4, 5, 6, 0, 0, 0, 0)
dot = 13

如此類推，不過 FFT 只能使用 double 運算，因此在儲存精度上不是很好拿捏，誤差大概在 1e-1 ~ 1e-2 之間。

感謝 liouzhou_101 的誤差指導

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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <string.h>
#include <algorithm>
#include <queue>
#include <stack>
#include <math.h>
#include <iostream>
#include <map>
#include <complex>
#include <cmath>
using namespace std;
#define for if (0); else for
using namespace std;
const int MaxFastBits = 16;
int **gFFTBitTable = 0;
int NumberOfBitsNeeded(int PowerOfTwo) {
	for (int i = 0;; ++i) {
		if (PowerOfTwo & (1 << i)) {
			return i;
		}
	}
}
int ReverseBits(int index, int NumBits) {
	int ret = 0;
	for (int i = 0; i < NumBits; ++i, index >>= 1) {
		ret = (ret << 1) | (index & 1);
	}
	return ret;
}
void InitFFT() {
    gFFTBitTable = new int *[MaxFastBits];
    for (int i = 1, length = 2; i <= MaxFastBits; ++i, length <<= 1) {
        gFFTBitTable[i - 1] = new int[length];
        for (int j = 0; j < length; ++j) {
            gFFTBitTable[i - 1][j] = ReverseBits(j, i);
		}
    }
}
inline int FastReverseBits(int i, int NumBits) {
    return NumBits <= MaxFastBits ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits);
}
void FFT(bool InverseTransform, vector<complex<double> >& In, vector<complex<double> >& Out) {
    if (!gFFTBitTable) { InitFFT(); }
    // simultaneous data copy and bit-reversal ordering into outputs
	int NumSamples = In.size();
    int NumBits = NumberOfBitsNeeded(NumSamples);
    for (int i = 0; i < NumSamples; ++i) {
		Out[FastReverseBits(i, NumBits)] = In[i];
    }
    // the FFT process
    double angle_numerator = acos(-1.) * (InverseTransform ? -2 : 2);
    for (int BlockEnd = 1, BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
        double delta_angle = angle_numerator / BlockSize;
        double sin1 = sin(-delta_angle);
        double cos1 = cos(-delta_angle);
        double sin2 = sin(-delta_angle * 2);
        double cos2 = cos(-delta_angle * 2);
        for (int i = 0; i < NumSamples; i += BlockSize) {
			complex<double> a1(cos1, sin1), a2(cos2, sin2);
            for (int j = i, n = 0; n < BlockEnd; ++j, ++n) {
				complex<double> a0(2 * cos1 * a1.real() - a2.real(), 2 * cos1 * a1.imag() - a2.imag());
				a2 = a1;
				a1 = a0;
				complex<double> a = a0 * Out[j + BlockEnd];
				Out[j + BlockEnd] = Out[j] - a;
				Out[j] += a;
            }
        }
        BlockEnd = BlockSize;
    }
    // normalize if inverse transform
    if (InverseTransform) {
        for (int i = 0; i < NumSamples; ++i) {
			Out[i] /= NumSamples;
        }
    }
}
vector<double> convolution(vector<double> a, vector<double> b) {
	int n = a.size();
	vector<complex<double> > s(n), d1(n), d2(n), y(n);
    for (int i = 0; i < n; ++i) {
        s[i] = complex<double>(a[i], 0);
	}
    FFT(false, s, d1);
    s[0] = complex<double>(b[0], 0);
    for (int i = 1; i < n; ++i) {
		s[i] = complex<double>(b[n - i], 0);
	}
    FFT(false, s, d2);
    for (int i = 0; i < n; ++i) {
		y[i] = d1[i] * d2[i];
    }
    FFT(true, y, s);
	vector<double> ret(n);
	for (int i = 0; i < n; ++i) {
		ret[i] = s[i].real();
	}
	return ret;
}
struct Polynomial {
    vector<double> v;
    Polynomial operator*(const Polynomial &other) const {
        int n = (int) max(v.size(), other.v.size()) * 2, m;
        for (m = 2; m < n; m <<= 1);
        vector<double> na, nb;
        na.resize(m, 0), nb.resize(m, 0);
        for (int i = 0; i < v.size(); i++)
            na[i] = v[i];
        for (int i = 0, j = m - 1; i < other.v.size(); i++, j--)
            nb[j] = other.v[i];
        Polynomial ret;
        ret.v = convolution(na, nb);
        for (int i = 1; i < ret.v.size(); i++)
            ret.v[i - 1] = ret.v[i];
        ret.v[ret.v.size() - 1] = 0;
        return ret;
    }
};
char sa[1<<18], sb[1<<18];
long long ret[1<<19];
int main() {
    while (scanf("%s %s", sa, sb) == 2) {
        Polynomial a, b, c;
        for (int i = (int)strlen(sa) - 1; i >= 0; i--)
            a.v.push_back(sa[i] - '0');
        
        for (int i = (int)strlen(sb) - 1; i >= 0; i--)
            b.v.push_back(sb[i] - '0');
        c = a * b;
        
        memset(ret, 0, sizeof(ret));
        int f = 0;
        double eps = 1.5e-1;
        for (int i = 0; i < c.v.size(); i++)
            ret[i] = (long long) (c.v[i] + eps);
        int n = (int)c.v.size();
        for (int i = 0; i < n; i++) {
            if (ret[i] >= 10) {
                ret[i + 1] += ret[i]/10;
                ret[i] %= 10;
            }
        }
        for (int i = n; i >= 0; i--) {
            if (ret[i])
                f = 1;
            if (f)
                printf("%lld", ret[i]);
        }
        if (!f)
            printf("0");
        puts("");
    }
    return 0;
}
/*
 0
 5
 
 11
 12
 */