题目

http://acm.hdu.edu.cn/showproblem.php?pid=1402

知识点

FFT 快速傅立叶变换

算法导论(第2版)第30章的《多项式与快速傅里叶变换》 

多项式有系数表示法和点值表示法。 两个次数界为n的多项式A(x)和B(x)相乘,输入输出均采用系数表示法。(假定n为2的幂)

1)使次数界增加一倍:A(x)和B(x)扩充为次数界为2n的多项式,并构造起系数表示

2)求值:两次应用2n阶FFT,计算出A(x)和B(x)的长度为2n的点值表示

3)点乘:计算多项式C(x)=A(x)*B(x)的点值表示

4)插值:对2n个点值对应用一次FFT计算出其逆DFT,就可以构造出多项式C(x)的系数表示

第1步和第3步需要时间为O(n),第2步和第4步运用FFT需要时间为O(nlgn)。 第2步求取n个点的值所需时间为O(n^2),如果我们巧妙地选取x(k)则其运行时间变为O(nlgn),FFT主要利用了单位复根(w^n=1的w就是n次单位复根,他们是e^(2PI*k/n),k=0,1,……n-1)的特殊性质。

代码

/*
 * Author: Haipz
 * School: HDU
 * File Name: 1402.cpp
 */
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cctype>
#include <cfloat>
#include <cstdlib>
#include <cstring>
#include <climits>
#include <iostream>
#include <vector>
#include <string>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
using namespace std;
const int N = 2e5, L = 5e4;
const double pi = acos(-1.0);
char a[L + 5], b[L + 5];
struct Complex {
    double r, i;
    Complex(double real = 0.0, double image = 0.0) {
        r = real, i = image;
    }
    Complex operator + (const Complex o) {
        return Complex(r + o.r, i + o.i);
    }
    Complex operator - (const Complex o) {
        return Complex(r - o.r, i - o.i);
    }
    Complex operator * (const Complex o) {
        return Complex(r*o.r - i*o.i, r*o.i + i*o.r);
    }
} x1[N + 5], x2[N + 5];
int sum[N + 5];
void brc(Complex *y, int l) {
    for (int i = 1, j = l/2; i < l - 1; ++i) {
        if (i < j) swap(y[i], y[j]);
        int k = l/2;
        while (j >= k) j -= k, k /= 2;
        if (j < k) j += k;
    }
    return;
}
void fft(Complex *y, int l, double on) {
    brc(y, l);
    for (int i = 2; i <= l; i *= 2) {
        Complex wn(cos(on*2*pi/i), sin(on*2*pi/i));
        for (int j = 0; j < l; j += i) {
            Complex w(1.0, 0.0);
            for (int k = j; k < j + i/2; ++k) {
                Complex u = y[k], t = w*y[k + i/2];
                y[k] = u + t, y[k + i/2] = u - t, w = w*wn;
            }
        }
    }
    if (on == -1)
        for (int i = 0; i < l; ++i) y[i].r /= l;
    return;
}
int main() {
    while (scanf("%s%s", a, b) != EOF) {
        int l1 = strlen(a), l2 = strlen(b), l = 1;
        while (l < l1*2 || l < l2*2) l *= 2;
        for (int i = 0; i < l; ++i) x1[i].r = x1[i].i = x2[i].r = x2[i].i= 0.0;
        for (int i = 0; i < l1; ++i) x1[i].r = a[l1 - i - 1] - '0', x1[i].i = 0.0;
        for (int i = 0; i < l2; ++i) x2[i].r = b[l2 - i - 1] - '0', x2[i].i = 0.0;
        fft(x1, l, 1), fft(x2, l, 1);
        for (int i = 0; i < l; ++i) x1[i] = x1[i]*x2[i];
        fft(x1, l, -1);
        for (int i = 0; i < l; ++i) sum[i] = x1[i].r + 0.5;
        for (int i = 0; i < l; ++i) sum[i + 1] += sum[i]/10, sum[i] %= 10;
        l = l1 + l2 - 1;
        while (sum[l] <= 0 && l > 0) --l;
        for (int i = l; i >= 0; --i) printf("%d", sum[i]);
        printf("\n");
    }
    return 0;
}
转载保留版权:http://haipz.com/blog/i/1486 - 海胖博客