HDU 1402 POJ 2389 BZOJ 2179 大整数乘法 FFT

FFT果然还是看算导来得快。。

网上的讲解糊里糊涂的

不过我还是看了2小时QAQ

递归形式

function recursive-fft(a) {
n = a.length // 必为2的幂
if (n == 1)
return a
wn = Math.pow(Math.e, 2 * Math.PI * Math.i / n)
// wn还可以利用复数的三角表示法：
wn = Math.cos(2 * Math.PI / n) + Math.i * Math.sin(2 * Math.PI / n)
w = 1
a0 = [ a[0], a[2], ..., a[n - 2] ]
a1 = [ a[1], a[3], ..., a[n - 1] ]
y0 = recursive-fft(a0)
y1 = recursive-fft(a1)
for (k = 0; k < n / 2; ++k) {
y[k] = y0[k] + w * y1[k]
y[k + n / 2] = y0[k] - w * y1[k]
w *= wn
}
return y
}

如果将递归树绘出，发现各元素下标的二进制反转后的数就是其新下标，于是划分为n/2块就可以依次合并。

function bit-reverse(a) {
n = a.length
for (k = 0; k < n; ++k)
A[rev(k)] = a[k]
return A
}

function iterative-fft(a) {
a = bit-reverse(a)
n = a.length
for (m = 2; m <= n; m *= 2) {
wm = complex(cos(2 * PI / m), sin(2 * PI / m))
for (k = 0; k < n; k += m) {
w = 1;
for (j = k; j < m / 2 + k; ++j) {
t = w * a[j + m / 2];
u = a[j]
a[j] = u + t
a[j + m / 2] = u - t
w *= wm
}
}
}
return a
}

HDU 1402 BZOJ 2179是必须用FFT了。。

POJ 2389就当是强行使程序变慢吧。。

偷懒用STL了。。

#include <cstdio>
#include <cstring>
#include <complex>
#include <algorithm>
#include <cmath>
#define rep(i,j,k) for(i=j;i<k;++i)
using namespace std;
typedef complex<double> cd;
const double PI = acos(-1);
const int N = 200010;
cd x
, y
;
char a
, b
;
int s
, rev
;
void fft(cd a[], int n, int ratio) {
int m, j, k;
rep(j,0,n) if (j < rev[j]) swap(a[j], a[rev[j]]);
for (m = 2; m <= n; m *= 2) {
cd wm(cos(ratio * 2 * PI / m), sin(ratio * 2 * PI / m));
for (k = 0; k < n; k += m) {
cd w(1, 0);
rep(j, k, k + m / 2) {
cd u = a[j], t = w * a[j + m / 2];
a[j] = u + t;
a[j + m / 2] = u - t;
w *= wm;
}
}
}
}
int main() {
int la, lb, n, i, k;
while (scanf("%s%s", a, b) == 2) {
la = strlen(a), lb = strlen(b); n = 1; k = -1;
while (n < la * 2 || n < lb * 2) n *= 2, ++k;
rep(i,0,n) rev[i] = (rev[i>>1]>>1)|((i&1)<<k);
rep(i,0,la) x[i] = cd(a[la - 1 - i] - '0', 0);
rep(i,la,n) x[i] = cd(0, 0);
rep(i,0,lb) y[i] = cd(b[lb - 1 - i] - '0', 0);
rep(i,lb,n) y[i] = cd(0, 0);
fft(x, n, 1); fft(y, n, 1);
rep(i,0,n) x[i] *= y[i];
fft(x, n, -1);
rep(i,0,n) s[i] = round(x[i].real() / n);
rep(i,0,n) s[i + 1] += s[i] / 10, s[i] %= 10;
n = la + lb - 1;
while (!s
&& n) --n;
for (i = n; i >= 0; --i) putchar(s[i] + '0');
putchar('\n');
}
return 0;
}