HDU 1402 A * B Problem Plus (FFT模板题)

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;
}