FFT模板题,求A*B。

用次FFT模板需要注意的是,N应为2的幂次,不然二进制平摊反转置换会出现死循环。

取出结果值时注意精度,要加上eps才能A。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
typedef long long ll;
const double pi = acos(-1.0);
const int maxn = 50000 + 5;
const double eps = 1e-6; struct Complex {
double a, b;
Complex() {
}
Complex(double a, double b) :
a(a), b(b) {
}
Complex operator +(const Complex& t) const {
return Complex(a + t.a, b + t.b);
}
Complex operator -(const Complex& t) const {
return Complex(a - t.a, b - t.b);
}
Complex operator *(const Complex& t) const {
return Complex(a * t.a - b * t.b, a * t.b + b * t.a);
}
}; // 二进制平摊反转置换
void brc(Complex *x, int n) {
int i, j, k;
for (i = 1, j = n >> 1; i < n - 1; i++) {
if (i < j)
swap(x[i], x[j]); k = n >> 1;
while (j >= k) {
j -= k;
k >>= 1;
}
if (j < k)
j += k;
}
} // FFT,其中on==1时为DFT,on==-1时为IDFT
void FFT(Complex *x, int n, int on) {
int h, i, j, k, p;
double r;
Complex u, t;
brc(x, n);
for (h = 2; h <= n; h <<= 1) { // 控制层数
r = on * 2.0 * pi / h;
Complex wn(cos(r), sin(r));
p = h >> 1;
for (j = 0; j < n; j += h) {
Complex w(1, 0);
for (k = j; k < j + p; k++) {
u = x[k];
t = w * x[k + p];
x[k] = u + t;
x[k + p] = u - t;
w = w * wn;
}
}
}
if (on == -1) // IDFT
for (i = 0; i < n; i++)
x[i].a = x[i].a / n + eps;
} int n, ma, N;
Complex x1[maxn<<2], x2[maxn<<2];
char sa[maxn], sb[maxn];
int ans[maxn<<1]; void solve() {
int n1 = strlen(sa), n2 = strlen(sb);
int N = 1, tmpn = max(n1, n2) << 1;
// N应为2的幂次
while(N < tmpn) N <<= 1;
for(int i = 0;i < N; i++)
x1[i].a = x1[i].b = x2[i].a = x2[i].b = 0;
for(int i = 0;i < n1; i++)
x1[i].a = sa[n1-i-1] - '0';
for(int i = 0;i < n2; i++)
x2[i].a = sb[n2-i-1] - '0';
FFT(x1, N, 1); FFT(x2, N, 1);
for(int i = 0;i < N; i++)
x1[i] = x1[i]*x2[i];
FFT(x1, N, -1);
int pre = 0, top = 0;
for(int i = 0;i < n1+n2; i++) {
// 不加epsA不了~
int cur = (int)(x1[i].a + eps);
ans[++top] = (cur + pre)%10;
pre = (pre + cur)/10;
}
while(!ans[top] && top > 1) top--;
for(int i = top;i >= 1; i--)
printf("%d", ans[i]);
puts("");
} int main() {
while(scanf("%s%s", sa, &sb) != -1) {
solve();
}
return 0;
}
04-27 21:39