FFT

位逆序置换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1.0);
struct Complex
{
    double real, imag;
};
void change(Complex y[], int len)
{
    for (int i = 1, j = len / 2, k; i < len - 1; i++)
    {
        if (i < j)
            swap(y[i], y[j]);
        k = len / 2;
        while (j >= k)
        {
            j = j - k;
            k = k / 2;
        }
        if (j < k)
            j += k;
    }
}
int main()
{
    int n;
    cin >> n;
    Complex a[n];
    for (int i = 0; i < n; i++)
    {
        cin >> a[i].real;
        a[i].imag = 0;
    }
    change(a, n);
    for (int i = 0; i < n; i++)
        cout << a[i].real << " ";
    return 0;
}

DFT和IDFT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#include <bits/stdc++.h>
using namespace std;

typedef complex<double> cd;
const double PI = acos(-1.0);
void fft(vector<cd> &a, bool invert)
{
    int n = a.size();
    for (int i = 1, j = 0; i < n; ++i)
    {
        int bit = n >> 1;
        while (j >= bit)
        {
            j -= bit;
            bit >>= 1;
        }
        if (j < bit)
            j += bit;
        if (i < j)
            swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1)
    {
        double angle = 2 * PI / len * (invert ? -1 : 1);
        cd wlen(cos(angle), sin(angle));
        for (int i = 0; i < n; i += len)
        {
            cd w(1);
            for (int j = 0; j < len / 2; ++j)
            {
                cd u = a[i + j];
                cd v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert)
        for (cd &x : a)
            x /= n;
}
string format_number(double num)
{
    if (fabs(num) < 1e-8)
        num = 0.0;
    stringstream ss;
    ss << fixed << setprecision(2) << num;
    return ss.str();
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int k;
    cin >> k;
    int N = 1 << k;
    vector<int> c(N, 0);
    for (int i = 0; i < N; ++i)
        cin >> c[i];
    
    vector<cd> a(N);
    for (int i = 0; i < N; ++i)
        a[i] = cd(c[i], 0.0);
    fft(a, true);

    for (int i = 0; i < N; i++)
    {
        double real_part = a[i].real();
        double imag_part = a[i].imag();

        if (fabs(real_part) < 1e-8)
            real_part = 0.0;
        if (fabs(imag_part) < 1e-8)
            imag_part = 0.0;
        real_part = round(real_part * 100.0) / 100.0;
        imag_part = round(imag_part * 100.0) / 100.0;
        string real_str = format_number(real_part);
        string imag_str = format_number(imag_part);
        cout << real_str << " " << imag_str << "\n";
    }
    return 0;
}

乘法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <bits/stdc++.h>
#include <complex>
using namespace std;

typedef long long ll;
typedef complex<double> cd;
const double PI = acos(-1.0);

void fft(vector<cd> &a, bool invert)
{
    int n = a.size();

    for (int i = 1, j = 0; i < n; ++i)
    {
        int bit = n >> 1;
        while (j >= bit)
        {
            j -= bit;
            bit >>= 1;
        }
        if (j < bit)
            j += bit;
        if (i < j)
            swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1)
    {
        double angle = 2 * PI / len * (invert ? -1 : 1);
        cd wlen(cos(angle), sin(angle));
        for (int i = 0; i < n; i += len)
        {
            cd w(1);
            for (int j = 0; j < len / 2; ++j)
            {
                cd u = a[i + j];
                cd v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }

    if (invert)
    {
        for (cd &x : a)
            x /= n;
    }
}

string mult(string a, string b)
{
    reverse(a.begin(), a.end());
    reverse(b.begin(), b.end());

    vector<int> num1, num2;
    for (char c : a)
        num1.push_back(c - '0');
    for (char c : b)
        num2.push_back(c - '0');

    int n = 1;
    while (n < num1.size() + num2.size())
        n <<= 1;

    vector<cd> fa(n, 0), fb(n, 0);
    for (int i = 0; i < num1.size(); i++)
        fa[i] = num1[i];
    for (int i = 0; i < num2.size(); i++)
        fb[i] = num2[i];

    fft(fa, false);
    fft(fb, false);

    for (int i = 0; i < n; i++)
        fa[i] *= fb[i];

    fft(fa, true);

    vector<int> result(n, 0);
    for (int i = 0; i < n; i++)
        result[i] = round(fa[i].real());
    ll carry = 0;
    for (int i = 0; i < n; i++)
    {
        ll temp = result[i] + carry;
        result[i] = temp % 10;
        carry = temp / 10;
    }
    string ans = "";
    int i = n - 1;
    while (i > 0 && result[i] == 0)
        i--;
    for (; i >= 0; --i)
        ans += to_string(result[i]);
    return ans;
}
int main()
{
    string a, b;
    cin >> a >> b;

    if ((a.size() == 1 && a[0] == '0') || (b.size() == 1 && b[0] == '0'))
        cout << "0" << endl;
    string ans = mult(a, b);
    cout << ans << endl;
    return 0;

}

多项式相乘

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include <cmath>
#include <cstdio>
#define R register
#define I inline
using namespace std;
const int N = 4.2e6;
const double PI = acos(-1);
int n, r[N];
struct C
{ // 手写complex，比STL快一点点
    double r, i;
    I C() { r = i = 0; }
    I C(R double x, R double y)
    {
        r = x;
        i = y;
    }
    I C operator+(R C &x) { return C(r + x.r, i + x.i); }
    I C operator-(R C &x) { return C(r - x.r, i - x.i); }
    I C operator*(R C &x) { return C(r * x.r - i * x.i, r * x.i + i * x.r); }
    I void operator+=(R C &x)
    {
        r += x.r;
        i += x.i;
    }
    I void operator*=(R C &x)
    {
        R double t = r;
        r = r * x.r - i * x.i;
        i = t * x.i + i * x.r;
    }
} f[N], g[N];
I int in()
{
    R char c = getchar();
    while (c < '-')
        c = getchar();
    return c & 15;
}
I void FFT(R C *a, R int op)
{
    R C W, w, t, *a0, *a1;
    R int i, j, k;
    for (i = 0; i < n; ++i) // 根据映射关系交换元素
        if (i < r[i])
            t = a[i], a[i] = a[r[i]], a[r[i]] = t;
    for (i = 1; i < n; i <<= 1)                                               // 控制层数
        for (W = C(cos(PI / i), sin(PI / i) * op), j = 0; j < n; j += i << 1) // 控制一层中的子问题个数
            for (w = C(1, 0), a1 = i + (a0 = a + j), k = 0; k < i; ++k, ++a0, ++a1, w *= W)
                t = *a1 * w, *a1 = *a0 - t, *a0 += t; // 蝴蝶操作
}
int main()
{
    R int m, i, l = 0;
    scanf("%d%d", &n, &m);
    for (i = 0; i <= n; ++i)
        f[i].r = in();
    for (i = 0; i <= m; ++i)
        g[i].r = in();
    for (m += n, n = 1; n <= m; n <<= 1, ++l)
        ;
    for (i = 0; i < n; ++i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); // 递推求r
    FFT(f, 1);
    FFT(g, 1);
    for (i = 0; i < n; ++i)
        f[i] *= g[i];
    FFT(f, -1);
    for (i = 0; i <= m; ++i)
        printf("%.0f ", fabs(f[i].r) / n);
    puts("");
    return 0;
}