FFT & NTT 学习 模板

参考资料：

算导第 30 章

http://www.gatevin.moe/acm/fft%E7%AE%97%E6%B3%95%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/

http://blog.csdn.net/acdreamers/article/details/39026505

先来两道裸的

hdu 1402 (DFT)

Code:

#include <bits/stdc++.h>#define PI acos(-1.0)#define maxn 140010#define maxl 50010using namespace std;inline int max(int a, int b) { return a > b ? a : b; }struct Complex {    double real, image;    Complex(double _real = 0, double _image = 0) : real(_real), image(_image) {}};Complex operator + (const Complex& c1, const Complex& c2) { return Complex(c1.real + c2.real, c1.image + c2.image); }Complex operator - (const Complex& c1, const Complex& c2) { return Complex(c1.real - c2.real, c1.image - c2.image); }Complex operator * (const Complex& c1, const Complex& c2) { return Complex(c1.real * c2.real - c1.image * c2.image, c1.real * c2.image + c2.real * c1.image); }ostream& operator << (ostream& out, Complex c) { out << c.real << "+ i * " << c.image << endl;}Complex a[maxn], b[maxn], A[maxn];char s1[maxl], s2[maxl];int ans[maxn];int rev(int x, int len) {    int ret = 0, mask = 1;    for (int i = 0; i < len; ++i) {        ret <<= 1;        if (mask & x) ret |= 1;        mask <<= 1;    }    return ret;}void dft(Complex* a, int len, int D) {    int h = (int)((double)log(len) / log(2) + 0.5);    for (int i = 0; i < len; ++i) {        A[rev(i, h)] = a[i];    }//    for (int i = 0; i < len; ++i) cout << A[i]; cout << "\n";    for (int s = 1; s <= h; ++s) {        int m = 1 << s;        Complex wm = Complex(cos(D * 2 * PI / m), sin(D * 2 * PI / m));        for (int k = 0; k < len; k += m) {            Complex w = Complex(1, 0);            for (int j = 0; j < (m >> 1); ++j) {                Complex temp = w * A[k + j + (m >> 1)];                Complex ori = A[k + j];                A[k + j] = ori + temp;                A[k + j + (m >> 1)] = ori - temp;                w = w * wm;            }        }    }    if (D == -1) {        for (int i = 0; i < len; ++i) A[i].real /= len, A[i].image /= len;    }    for (int i = 0; i < len; ++i) a[i] = A[i];}void work() {    int len1 = strlen(s1), len2 = strlen(s2), len = len1 + len2;    int n = 1;    while (n < len) n <<= 1;    len = n;    memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b));    for (int i = 0; i < len1; ++i) a[i] = Complex(s1[len1 - 1 - i] - '0', 0);    for (int i = 0; i < len2; ++i) b[i] = Complex(s2[len2 - 1 - i] - '0', 0);    dft(a, len, 1); dft(b, len, 1);    for (int i = 0; i < len; ++i) a[i] = a[i] * b[i];    dft(a, len, -1);    memset(ans, 0, sizeof(ans));    for (int i = 0; i < len; ++i) ans[i] = (int)(a[i].real + 0.5);    for (int i = 0; i < len; ++i) {        ans[i + 1] += ans[i] / 10;        ans[i] %= 10;    }    if (ans[len]) ++len;    int i = len - 1;    while (i >= 0 && !ans[i]) --i;    if (i == -1) { printf("0\n"); return; }    for (; i >= 0; --i) printf("%d", ans[i]);    printf("\n");}int main() {    freopen("in.txt", "r", stdin);    while (scanf("%s%s", s1, s2) != EOF) work();    return 0;}

51nod 1028 (NTT)

Code:

#include <bits/stdc++.h>#define maxn 300010typedef long long LL;using namespace std;LL a[maxn], b[maxn], A[maxn], wn[22];char s1[maxn], s2[maxn];const int N = 1 << 18;const int P = (479 << 21) + 1;const int G = 3;const int NUM = 20;LL poww(LL a, LL b) {    LL ret = 1;    while (b) {        if (b & 1) ret = ret * a % P;        a = a * a % P;        b >>= 1;    }    return ret;}LL pre() {    for (int i = 0; i < NUM; ++i) {        int t = 1 << i;        wn[i] = poww(G, (P - 1) / t);    }}int geta(char* s, LL* a) {    int len = strlen(s);    for (int i = 0; i < len; ++i) a[i] = s[len - 1 - i] - '0';    return len;}int rev(int x, int len) {    int ret = 0, mask = 1;    for (int i = 1; (1 << i) <= len; ++i) {        ret <<= 1;        if (mask & x) ret |= 1;        mask <<= 1;    }    return ret;}void ntt(LL* a, int len, int N) {    int id = 0;    memset(A, 0, sizeof(A));    for (int i = 0; i < len; ++i) A[rev(i, len)] = a[i];    for (int i = 1; (1 << i) <= len; ++i) {        int m = 1 << i;        ++id;        for (int j = 0; j < len; j += m) {            LL w = 1;            for (int k = 0; k < (m >> 1); ++k) {                LL temp = w * A[j + k + (m >> 1)] % P;                LL ori = A[j + k] % P;                A[j + k] = (ori + temp) % P;                A[j + k + (m >> 1)] = (ori - temp + P) % P;                w = w * wn[id] % P;            }        }    }//    for (int i = 0; i < len; ++i) printf("%lld", A[i]); printf("\n");    if (N == -1) {        for (int i = 1, j = len - 1; i < j; ++i, --j) swap(A[i], A[j]);        LL inv = poww(len, P - 2);        for (int i = 0; i < len; ++i) A[i] = A[i] * inv % P;    }    for (int i = 0; i < len; ++i) a[i] = A[i];}void work() {    memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b));    int len1 = geta(s1, a), len2 = geta(s2, b);    int len = len1 + len2, len0 = 1;    while (len0 < len) len0 <<= 1;    ntt(a, len0, 1); ntt(b, len0, 1);    for (int i = 0; i < len0; ++i) a[i] = a[i] * b[i] % P;    ntt(a, len0, -1);    for (int i = 0; i < len; ++i) {        a[i + 1] += a[i] / 10;        a[i] %= 10;    }    if (a[len]) ++len;    int p = len - 1;    while (p >= 0 && a[p] == 0) --p;    if (p == -1) { printf("0\n"); return; }    for (int i = p; i >= 0; --i) printf("%lld", a[i]); printf("\n");}int main() {    pre();    while (scanf("%s%s", s1, s2) != EOF) work();    return 0;}