从N方到NlogN的转变——FFT

1.Why#

为什么我们信息学竞赛要用到FFT
因为我们要优化卷积啊
将n边为log是一个非常大的优化啊

What#

什么是FFT呢在wiki上是这么说的

快速傅里叶变换（英语：Fast Fourier Transform, FFT），是计算序列的离散傅里叶变换（DFT）或其逆变换的一种算法。傅里叶分析将信号从原始域（通常是时间或空间）转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏（大多为零）因子之积来快速计算此类变换。[1] 因此，它能够将计算DFT的复杂度从只用DFT定义计算需要的 $O(n^{2})$ ，降低到 $O(n\log n)$ ，其中 n 为数据大小。

啥你说啥我听不懂（黑人问号?）

如果简单的说就是求两个多项式的乘积

How#

其实关于FFT的具体理论很复杂，我一篇博客肯定不能接受清楚所以请有兴趣的人移步Google 或 baidu

如果简单来说FFT用的分制的方法
所以才能有log啊

将一个多项式按奇偶分开
但如果我们直接找值去计算时间复杂度还是 $O(n^2)$ 的

所以出现了一个特殊的东西他叫做n次单位根

然后利用他的性质就可以搞了

先贴两个板子

UOJ 34 多项式乘法

1
#include <cstdio>
2
#include <cstring>
3
#include <algorithm>
4
#include <cmath>
5
using namespace std;
6
const int maxn = 1 << 18 | 5;
7
namespace FFT
8
{
9
const double pi = acos(-1.0);
10
struct Complex
11
{
12
    double x, y;
13
    Complex operator+(const Complex &a)
14
    {
15
        return (Complex){x + a.x, y + a.y};
16
    }
17
    Complex operator-(const Complex &a)
18
    {
19
        return (Complex){x - a.x, y - a.y};
20
    }
21
    Complex operator*(const Complex &a)
22
    {
23
        return (Complex){x * a.x - y * a.y, x * a.y + y * a.x};
24
    }
25
    Complex operator*(const double &a)
26
    {
27
        return (Complex){x * a, y * a};
28
    }
29
    Complex Get()
30
    {
31
        return (Complex){x, -y};
32
    }
33
} c[maxn], c1[maxn], V = (Complex){0.5, 0} * (Complex){0, -0.5};
34
int rev[maxn], N, FU;
35
double Inv;
36
void FFt(Complex *a, int op)
37
{
38
    Complex w, wn, t;
39
    for (int i = 0; i < N; i++)
40
        if (i < rev[i])
41
            swap(a[i], a[rev[i]]);
42
    for (int k = 2; k <= N; k <<= 1)
43
    {
44
        wn = (Complex){cos(op * 2 * pi / k), sin(op * 2 * pi / k)};
45
        for (int j = 0; j < N; j += k)
46
        {
47
            w = (Complex){1, 0};
48
            for (int i = 0; i < (k >> 1); i++, w = w * wn)
49
            {
50
                t = a[i + j + (k >> 1)] * w;
51
                a[i + j + (k >> 1)] = a[i + j] - t;
52
                a[i + j] = a[i + j] + t;
53
            }
54
        }
55
    }
56
    if (op == -1)
57
        for (int i = 0; i < N; i++)
58
            a[i] = a[i] * Inv;
59
}
60
void FFt(int *a, int *b, int *res, int n)
61
{
62
    int j;
63
    N = 1;
64
    while (N < n)
65
        N <<= 1;
66
    FU = N - 1;
67
    Inv = 1. / N;
68
    for (int i = 0; i < N; i++)
69
        if (i & 1)
70
            rev[i] = (rev[i >> 1] >> 1) | (N >> 1);
71
        else
72
            rev[i] = (rev[i >> 1] >> 1);
73
    for (int i = 0; i < N; i++)
74
        c[i].x = a[i], c[i].y = b[i];
75
    FFt(c, 1);
76
    for (int i = 0; i < N; i++)
77
        j = (N - i) & FU, c1[i] = (c[i] + c[j].Get()) * (c[i] - c[j].Get()) * V;
78
    FFt(c1, -1);
79
    for (int i = 0; i < N; i++)
80
        res[i] = round(c1[i].x);
81
}
82
} // FFT
83

84
int main(int argc, char const *argv[])
85
{
86
    static int a[maxn],b[maxn];
87
    int n1, n2, m;
88
    scanf("%d%d", &n1, &n2);
89
    m = n1 + n2 + 1;
90
    for (int i = 0; i <= n1; i++)
91
        scanf("%d", &a[i]);
92
    for (int i = 0; i <= n2; i++)
93
        scanf("%d", &b[i]);
94
    FFT::FFt(a, b, a, m);
95
    for (int i = 0; i < m; i++)
96
        printf("%d ", a[i]);
97
    return 0;
98
}

COGS 1473. 很强的乘法问题

1
#include <cstring>
2
#include <cstdio>
3
#include <algorithm>
4
#include <cmath>
5
const double pi = acos(-1.);
6
struct Complex
7
{
8
    double x, y;
9
    Complex() { ; }
10
    Complex(double a, double b) : x(a), y(b) {}
11
    Complex operator+(const Complex &a) { return Complex(x + a.x, y + a.y); }
12
    Complex operator-(const Complex &a) { return Complex(x - a.x, y - a.y); }
13
    Complex operator*(const Complex &a) { return Complex(x * a.x - y * a.y, x * a.y + y * a.x); }
14
    Complex operator*(const double a) { return Complex(x * a, y * a); }
15
    Complex Get() { return Complex(x, -y); }
16
} A[150005 << 3], B[150005 << 3];
17
int rev[150005 << 3], N, FU;
18
double INV;
19
void FFt(Complex *a, int op)
20
{
21
    Complex w, wn, t;
22
    for (int i = 0; i < N; i++)
23
        if (i < rev[i])
24
            std::swap(a[i], a[rev[i]]);
25
    for (int k = 2; k <= N; k <<= 1)
26
    {
27
        wn = Complex(cos(op * 2 * pi / k), sin(op * 2 * pi / k));
28
        for (int j = 0; j < N; j += k)
29
        {
30
            w = Complex(1, 0);
31
            for (int i = 0; i < (k >> 1); i++, w = w * wn)
32
            {
33
                t = a[i + j + (k >> 1)] * w;
34
                a[i + j + (k >> 1)] = a[i + j] - t;
35
                a[i + j] = a[i + j] + t;
36
            }
37
        }
38
    }
39
    if (op == -1)
40
        for (int i = 0; i < N; i++)
41
            a[i] = a[i] * INV;
42
}
43
void FFt(const int *a, const int *b, int *res, int n)
44
{
45
    N = 1;
46
    while (N < n)
47
        N <<= 1;
48
    INV = 1. / N;
49
    for (int i = 0; i < N; i++)
50
        if (i & 1)
51
            rev[i] = (rev[i >> 1] >> 1) | (N >> 1);
52
        else
53
            rev[i] = (rev[i >> 1] >> 1);
54
    for (int i = 0; i < N; i++)
55
        A[i].x = a[i], B[i].x = b[i];
56
    FFt(A, 1), FFt(B, 1);
57
    for (int i = 0; i < N; i++)
58
        A[i] = A[i] * B[i];
59
    FFt(A, -1);
60
    for (int i = 0; i < N; i++)
61
        res[i] = round(A[i].x);
62
}
63
char s[150005 << 2];
64
struct BigNum
65
{
66
    int len;
67
    int a[1000005];
68
    void read()
69
    {
70
        scanf("%s", s);
71
        len = strlen(s);
72
        for (int i = len - 1, j = 0; i >= 0; i--, j++)
73
            a[j] = s[i] - '0';
74
    }
75
    BigNum operator*(const BigNum &b)
76
    {
77
        BigNum ans;
78
        ans.len = len + b.len - 1;
79
        FFt(a, b.a, ans.a, ans.len + 1);
80
        for (int i = 0; i <= ans.len + 2; i++)
81
        {
82
            if (ans.a[i] > 9)
83
            {
84
                ans.a[i + 1] += ans.a[i] / 10;
85
                ans.a[i] %= 10;
86
            }
87
        }
88
        while (ans.a[ans.len])
89
            ans.len++;
90
        return ans;
91
    }
92
    void print()
93
    {
94
        for (int i = len - 1; i >= 0; i--)
95
            printf("%d", a[i]);
96
        printf("\n");
97
    }
98
} a, b;
99

100
int main()
101
{
102
    freopen("bettermul.in","r",stdin);
103
    freopen("bettermul.out","w",stdout);
104
    a.read();
105
    b.read();
106
    (a * b).print();
107
    //while (1)
108
        ;
109
}