1402大数相乘 FFT算法学习！
网上找的，慢慢研究下
View Code
1 #include <cmath>
2 #include <complex>
3 #include <cstring>
4 using namespace std;
5
6 const double PI = acos(-1);
7 typedef complex<double> cp;
8 typedef long long int64;
9
10 const int N = 1 << 16;
11 int64 a[N], b[N], c[N << 1];
12
13 void bit_reverse_copy(cp a[], int n, cp b[])
14 {
15 int i, j, k, u, m;
16 for (u = 1, m = 0; u < n; u <<= 1, ++m);
17 for (i = 0; i < n; ++i)
18 {
19 j = i; k = 0;
20 for (u = 0; u < m; ++u, j >>= 1)
21 k = (k << 1) | (j & 1);
22 b[k] = a[i];
23 }
24 }
25
26 void FFT(cp _x[], int n, bool flag)
27 {
28 static cp x[N << 1];
29 bit_reverse_copy(_x, n, x);
30 int i, j, k, kk, p, m;
31 for (i = 1, m = 0; i < n; i <<= 1, ++m);
32 double alpha = 2 * PI;
33 if (flag) alpha = -alpha;
34 for (i = 0, k = 2; i < m; ++i, k <<= 1)
35 {
36 cp wn = cp(cos(alpha / k), sin(alpha / k));
37 for (j = 0; j < n; j += k)
38 {
39 cp w = 1, u, t;
40 kk = k >> 1;
41 for (p = 0; p < kk; ++p)
42 {
43 t = w * x[j + p + kk];
44 u = x[j + p];
45 x[j + p] = u + t;
46 x[j + p + kk] = u - t;
47 w *= wn;
48 }
49 }
50 }
51 memcpy(_x, x, sizeof(cp) * n);
52 }
53
54 void polynomial_multiply(int64 a[], int na, int64 b[], int nb, int64 c[], int &nc)
55 {
56 int i, n;
57 i = max(na, nb) << 1;
58 for (n = 1; n < i; n <<= 1);
59 static cp x[N << 1], y[N << 1];
60 for (i = 0; i < na; ++i) x[i] = a[i];
61 for (; i < n; ++i) x[i] = 0;
62 FFT(x, n, 0);
63 for (i = 0; i < nb; ++i) y[i] = b[i];
64 for (; i < n; ++i) y[i] = 0;
65 FFT(y, n, 0);
66 for (i = 0; i < n; ++i) x[i] *= y[i];
67 FFT(x, n, 1);
68 for (i = 0; i < n; ++i)
69 {
70 c[i] =(int64)(x[i].real() / n + 0.5);
71 }
72 for (nc = na + nb - 1; nc > 1 && !c[nc - 1]; --nc);
73 }
74
75 const int LEN = 5, MOD = 100000;
76 void convert(char *s, int64 a[], int &n)
77 {
78 int len = strlen(s), i, j, k;
79 for (n = 0, i = len - LEN; i >= 0; i -= LEN)
80 {
81 for (j = k = 0; j < LEN; ++j)
82 k = k * 10 + (s[i + j] - '0');
83 a[n++] = k;
84 }
85 i += LEN;
86 if (i)
87 {
88 for (j = k = 0; j < i; ++j)
89 k = k * 10 + (s[j] - '0');
90 a[n++] = k;
91 }
92 }
93
94 void print(int64 a[], int n)
95 {
96 printf("%I64d", a[--n]);
97 while (n) printf("%05I64d", a[--n]);
98 puts("");
99 }
100
101 char buf[N + 10];
102
103 int main()
104 {
105 int na, nb, nc;
106
107 while (scanf("%s", buf) != EOF)
108 {
109 bool sign = false;
110 if (buf[0] == '-')
111 {
112 sign = !sign;
113 convert(buf + 1, a, na);
114 }
115 else convert(buf, a, na);
116 scanf("%s", buf);
117 if (buf[0] == '-')
118 {
119 sign = !sign;
120 convert(buf + 1, b, nb);
121 }
122 else convert(buf, b, nb);
123 polynomial_multiply(a, na, b, nb, c, nc);
124 int64 t1, t2;
125 t1 = 0;
126 for (int i = 0; i < nc; ++i)
127 {
128 t2 = t1 + c[i];
129 c[i] = t2 % MOD;
130 t1 = t2 / MOD;
131 }
132 for (; t1; t1 /= MOD) c[nc++] = t1 % MOD;
133 if (sign) putchar('-');
134 print(c, nc);
135 }
136 return 0;
137 }
