# 【多项式】FFT

## Pre-knowledge

### 多项式

#### Definition

$$f(x) = \sum_{i = 0}^{n} a_i \times x^i$$

### Operation

$$f(x)~ = ~A(x) \pm B(x)~ = ~\sum_{i = 0}^{n} (a_i \pm b_i) x^i$$

$$f(x)~=~A(x) \times B(x)~=~\sum_{i = 0}^n \sum_{j = 0}^n a_i \times b_j \times x^{i + j}$$

$$f(x) = \sum_{i = 0}^{2n} \sum_{j = 0}^{\min(i, n)} a_j \times b_{i - j} \times x^i$$

### 复数

#### Definition

$$i^2 = -1$$

$$z = a + b \times i,~~~a, b \in R$$

#### Operation

$$|z| = \sqrt{a^2 + b^2}$$

$z$ 与 $\overline z$ 满足如下关系：

$$\theta_0 + \theta_1 = \pi$$

$$|z| = |\overline z|$$

$$z_0~=~z_1 \pm z_2~=~(a_1 \pm a_2) + (b_1 \pm b_2)i$$

$$z_0~=~z_1 \times z_2~=~(a_1 + b_1i) \times (a_2 + b_2i) = a_1a_2 + a_1b_2i + a_2b_1i + b_1b_2i^2$$

$$i^2 = -1$$

$$z_0~=~(a_1a_2 - b_1b_2) + (a_1b_2 + a_2b_1) i$$

$$\theta_0 = \theta_1 + \theta_2$$

$$|z_0| = |z_1| \times |z_2|$$

$$z \times \overline z~=~(a + bi) \times (a - bi) = a^2 + b^2$$

$$z_0 = \frac{z_1}{z_2}$$

$$z_0~=~\frac{z_1 \overline {z_2}}{a_2^2 + b_2^2}$$

$$e^{i\theta} = \cos \theta + i \sin \theta$$

$$e^{i\pi} = \cos \pi + i \sin \pi$$

### 单位根

#### Definition

$$x_k~=~e^{i \frac{2 k\pi}{n}}$$

#### Proof

$$x_n^k~=~e^{i 2 k\pi}$$

$$e^{i2k\pi} = \cos 2k\pi + i \sin 2k\pi$$

$$\cos 2k\pi = 1,~~\sin 2k\pi~=~0$$

$$x_n^k = e^{i\frac{2k\pi}{n}}$$

#### Property

$n$ 个单位根在复平面上平分单位圆。他们与单位圆的 $n$ 等分线所在线段分别重合。

### 本原单位根

#### Definition

$0$ 到 $(n - 1)$ 次方的值能生成所有 $n$ 次单位根的 $n$ 次单位根称为为 $n$ 次本原单位根。

$n$ 次本原单位根可能不止一个，但是下文的“本原单位根”特指 $e^{i \frac{2\pi}{n}}$

#### Property

$$(\omega_{n}^k)^2~=~\omega_m^k$$

$$\omega_n^{m + k}~=~-\omega_{n}^k$$

#### Proof

Pre-knowledge部分结束了。

## Algorithm

### DFT

$$b_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_{n}^i$$

### FFT

$$A(x)~=~\sum_{i = 0}^{n - 1} a_i x^i~=~\sum_{i = 0}^m a_{2i} \times x^{2i} + \sum_{i = 0}^m a_{2i + 1} \times x^{2^i + 1}$$

$$A(x)~=~\sum_{i = 0}^{m - 1} a_{2i} x^{i2} + x\sum_{i = 0}^{m - 1} a_{2i + 1} x^{2i}~=~\sum_{i = 0}^{m - 1} a_{2i} (x^{2})^{i} + x\sum_{i = 0}^{m - 1} a_{2i + 1} (x^2)^{i}$$

$$A_0~=~\sum_{i = 0}^{m - 1} a_{2i} x^i$$

$$A_1~=~\sum_{i = 0}^{m - 1} a_{2i + 1}x^i$$

$$A(x)~=~A_0(x^2) + x \times A_1(x^2)$$

$$(\omega_{n}^k)^2~=~\omega_m^k$$

$$\omega_n^{m + k}~=~-\omega_{n}^k$$

$$A(\omega_n^k)~=~A_0((\omega_n^k)^2) + w_n^kA_1((\omega_n^k)^2)$$

$$A(\omega_n^k)~=~A_0(\omega_m^k) + w_n^kA_1(\omega_m^k)$$

$$A(\omega_n^{m + k})~=~A_0((\omega_n^{m + k})^2) + \omega_n^{m + k} A_1((\omega_n^{m + k})^2)$$

$$A(\omega_n^{m + k})~=~A_0((w_n^k)^2)+-\omega_n^k A_1((\omega_n^k)^2)$$

$$A(\omega_n^{m + k})~=~A_0(\omega_m^k) - w_n^kA_1(\omega_m^k)$$

### IDFT

$$b_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_n^{ik}$$

$$a_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}$$

（这里使用截图的原因是你谷的多行公式渲染有点问题，建议查看大图）。

1、当 $j = k$ 时：

$$\sum_{i = 0}^{n - 1} \omega_n^{i(j - k)}~=~n \times 1 = n$$

2、当 $j \neq k$ 时：

$$\sum_{i = 0}^{n - 1} \omega_n^{(j - k) i}~=~\frac{1 - \omega_n^{(n - 1)(j - k)} \times \omega_n^{j - k}}{1 - \omega_{n}^{j - k}}~=~\frac{1 - \omega_n^{n(j - k)}}{1 - \omega_{n}^{j - k}}~=~\frac{1 - (\omega_n^{(j - k)})^n}{1 - \omega_{n}^{j - k}}$$

$\sum_{i = 0}^{n - 1} \omega_n^{(j - k) i}~=~\frac{1 - (\omega_n^{(j - k)})^0}{1 - \omega_{n}^{j - k}}~=~\frac{1 - 1}{1 - \omega_n^{j - k}}~=~0$

$$\sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}~=~~\sum_{j = 0}^{n - 1} a_j \times \sum_{i = 0}^{n - 1} \omega_n^{i(j - k)}~=~\sum_{j = 0}^{n - 1} a_j \times [j = k] \times n~=~a_k \times n$$

$$\text{右边}~=~\frac{1}{n} a_k \times n~=~a_k~=~\text{左边}$$

$$a_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}$$

DFT 的逆变换，称为 IDFT。我们可以通过这个式子求出这个多项式的系数表示法。

## Code

void FFT(std::complex<double> *A, int N) {
if (N == 1) {
return;
}
int M = N >> 1;
std::complex<double> A0[M], A1[M];
for (int i = 0; i < M; ++i) {
A0[i] = A[i << 1];
A1[i] = A[(i << 1) | 1];
}
FFT(A0, M); FFT(A1, M);
auto W = std::complex<double>(cos(1.0 * PI / M), sin(1.0 * PI / M)), w = std::complex<double>(1.0, 0.0);
for (int i = 0; i < M; ++i) {
A[i] = A0[i] + w * A1[i];
A[i + M] = A0[i] - w * A1[i];
w *= W;
}
}

void IFFT(std::complex<double> *A, int N) {
FFT(A, n)
std::reverse(A + 1, A + N);
}

## optimization

step 1: 0 1 2 3 4 5 6 7
step 2: 0 2 4 6,  1 3 5 7
step 3: 0 4,  2 6,  1 3,  5 7
step 4: 0,  4,  2,  6,  1,  3,  5,  7

000, 100, 010, 110, 001, 011, 101, 111

000, 001, 010, 011, 100, 110, 101, 111

0, 1, 2, 3, 4, 5, 6, 7

000, 100

000, 001

010, 110

000, 100, 010, 110

000, 001, 010, 011

void MakeRev(const int N) {
int d = N >> 1, p = 0;
tax[p++] = 0;
tax[p++] = d;
for (int w = 2; w <= N; w <<= 1) {
d >>= 1;
for (int p0 = 0; p0 < w; ++p0) {
tax[p++] = tax[p0] | d;
}
}
}

for (int i = 1; i < N; ++i) if (tax[i] > i) {
std::swap(A[i], A[tax[i]]);
}

## Final Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <complex>
#include <iostream>
#include <algorithm>

const int maxn = 6000006;
const double PI = acos(-1);

int n, m, k, sm;
int tax[maxn];
std::complex<double> F[maxn], G[maxn];

void MakeRev(const int N);
void FFT(std::complex<double> *A, int N);

int main() {
freopen("1.in", "r", stdin);
qr(n); qr(m);
for (int i = 0, x; i <= n; ++i) {
x = 0; qr(x); F[i] = x;
}
for (int i = 0, x; i <= m; ++i) {
x = 0; qr(x); G[i] = x;
}
sm = n + m; k = 1;
while (k <= sm) k <<= 1;
MakeRev(k);
FFT(F, k);
FFT(G, k);
for (int i = 0; i < k; ++i) {
F[i] *= G[i];
}
FFT(F, k);
std::reverse(F + 1, F + k);
for (int i = 0; i < sm; ++i) {
qw(int((F[i].real()) / k + 0.5), ' ', true);
}
qw(int(F[sm].real() / k + 0.5), '\n', true);
return 0;
}

void FFT(std::complex<double> *A, int N) {
for (int i = 1; i < N; ++i) if (tax[i] > i) {
std::swap(A[i], A[tax[i]]);
}
for (int len = 2, M = 1; len <= N; M = len, len <<= 1) {
std::complex<double> W(cos(PI / M), sin(PI / M)), w(1.0, 0.0);
for (auto L = 0, R = len - 1; R <= N; L += len, R += len) {
auto w0 = w;
for (auto p = L, lim = L + M; p < lim; ++p) {
auto x = A[p] + w0 * A[p + M], y = A[p] - w0 * A[p + M];
A[p] = x; A[p + M] = y;
w0 *= W;
}
}
}
}

void MakeRev(const int N) {
int d = N >> 1, p = 0;
tax[p++] = 0;
tax[p++] = d;
for (int w = 2; w <= N; w <<= 1) {
d >>= 1;
for (int p0 = 0; p0 < w; ++p0) {
tax[p++] = tax[p0] | d;
}
}
}