# 快速傅里叶变换（FFT）

## 虚数

\begin{aligned} (a+bi) + (c+di) &= (a+c) + (b+d)i \\ (a+bi)(c+di) &= (ac - bd) + (ad + bc)i \\ \cfrac {(a+bi)}{(c+di)} &= \cfrac {ac + bd}{c^2 + d^2} + \cfrac {bc - ad}{c^2 + d^2} \end{aligned}

### 重要性质

$e^{ix} = \cos x + i \sin x$

\begin{aligned} e^x &= 1 + x + \cfrac {x^2}{2!} + \cfrac {x^3}{3!} + \dots \\ \cos x &= 1 - \cfrac {x^2}{2!} + \cfrac {x^4}{4!} - \cfrac {x^6}{6!} + \dots \\ \sin x &= x - \cfrac {x^3}{3!} + \cfrac {x^5}{5!} - \cfrac {x^7}{7!} + \dots \\ \end{aligned}

$e^{ix} = 1 + ix - \cfrac {x^2}{2!} - \cfrac {ix^3}{3!} + \dots$

### 代码实现

struct Complex {
double real, imag;
Complex() : real(0), imag(0) {}
Complex(double re, double im) : real(re), imag(im) {};
inline Complex operator + (const Complex & b) {
return Complex(real + b.real, imag + b.imag);
}
inline Complex operator - (const Complex & b) {
return Complex(real - b.real, imag - b.imag);
}
inline Complex operator * (const Complex & b) {
return Complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
}
};


## 单位根

### 三个引理

• 消去定理：$$\omega_{dn}^{dk} = \omega_n^k$$

• 折半引理：$$(\omega_n^{k+\frac n2})^2 = (w_n^k)^2 = \omega_{\frac n2}^k$$

\begin{aligned} \omega_n^{k+\frac n2} &= \omega_n^k\omega_n^{\frac n2} = -\omega_n^k \\ (\omega_n^k)^2 &= \omega_n^{2k} = \omega_{\frac n2}^k \end{aligned}

• 求和引理：$$\sum_{i=0}^{n-1} (\omega_n^k)^i = 0$$

$\sum_{i=0}^{n-1} (\omega_n^k)^i = \cfrac {(w_n^k)^n - 1}{w_n^k - 1} = \cfrac {(w_n^n)^k - 1}{w_n^k - 1} = \cfrac {1^k - 1}{w_n^k - 1} = 0$

## 多项式

(OI中)一般形式$$F(x) = a_0 + a_1x + a_2 x^2 + \cdots + a_n x^n$$

\begin{aligned} A(x) &= \sum_{i = 0}^n a_i x^i \\ B(x) &= \sum_{i = 0}^n b_i x^i \\ \end{aligned}

• 加法：

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

• 乘法

\begin{aligned} c_i &= \sum_{j = 0}^{i} a_j b_{i - j} \\ A(x)B(x) &= \sum_{i = 0}^{2n} c_i x^i \\ \end{aligned}

• 系数表示

它将一个多项式表示成由其系数构成的向量的形式

例如 $$A = [a_0, a_1, a_2, \dots, a_n]^T$$

加法即为 $$A_1+ A_2$$，直接相加即可。时间复杂度 $$O(n)$$

乘法则做向量卷积，为 $$A_1 \otimes A_2$$。一般来说，时间复杂度为 $$O(n^2)$$

如果给定 $$x$$ 求值，则可以使用霍纳法则或者秦九昭算法。时间复杂度为 $$O(n)$$

• 点值表示

用至少 $$n$$ 个多项式上的点来表示

一般形式如 $$\{(x_0, A(x_0)), (x_1, A(x_1), \dots, (x_n, A(x_n))\}$$

进行运算是，一般要保证两个多项式在同一位置取值相同，即 $$x_i$$ 相同

加法运算直接将两点坐标相加即可，时间复杂度为 $$O(n)$$

乘法运算只需要将两点坐标相乘即可。时间复杂度为 $$O(n)$$，太好了！

如果我们需要 $$A(x)$$ ，这个过程叫做插值，可以通过拉格朗日插值公式进行计算，复杂度为 $$O(n^2)$$，这里不展开讲述。

## 离散傅里叶变换(DFT)

DFT(Discrete Fourier Transform) 是快速傅里叶变换（FFT）的基础，也是快速数论变换（NTT）的基础

$h(x) = \sum_{i = 0}^{n-1}c_i x^i$

$[\ h(\omega^0), h(\omega^1), h(\omega^2), \dots, h(\omega^{n-1})\ ]^T$

$\begin{bmatrix} h(\omega^0) \\ h(\omega^1) \\ h(\omega^2) \\ \vdots \\ h(\omega^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ 1 & w_n^2 & w_n^{2\times2} & w_n^{2\times3} & \cdots & w_n^{2\times(n-1)}\\ \vdots & \vdots & \vdots & \vdots & \ddots & \cdots \\ 1 & w_n^{(n-1)\times2} & w_n^{(n-1)\times3} & w_n^{(n-1)\times2} & \cdots & w_n^{(n-1)\times(n-1)} \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{n-1} \end{bmatrix}$

### 重要性质（由点值到多项式）

\begin{aligned} g(\omega^{-k}) &= h(\omega^0) \omega^{-0k} + h(\omega_1) \omega^{-k} + h(w_2) \omega^{-2k} + \cdots + h(\omega^{n-1}) \omega^{-(n-1)k} \\ &= \sum_{i=0}^{n-1} h(\omega^i)\omega^{-ik} \\ &= \sum_{i=0}^{n-1} \omega^{-ik} \sum_{j = 0}^{n-1} c_j \omega^{ij} \\ &= \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} \omega^{(j-k)i} c_j \\ &= \sum_{i=0}^{n-1} c_j \sum_{j=0}^{n-1} \omega^{(j-k)i} \end{aligned}

• $$j = k$$，那么此时 $$\sum_{j=0}^{n-1} \omega^{(j-k)y} = \sum_{j=0}^{n-1} \omega^0 = n$$。对于 $$c_k$$ 做出的贡献为 $$n$$

• $$j \ne k$$，将 $$k$$ 看作常数，那么此时 $$\sum_{j=0}^{n-1} \omega^{(j-k)i} = \sum_{j=0}^{n-1} \omega^{j}$$。依据上文中求和引理，其值为 $$0$$，也就是对 $$c_j$$ 做出的贡献为 $$0$$

• 将两个多项式改写为向量的形式，并分别对其做一次离散傅里叶变换 (DFT)

• 将变换过后的两个序列相乘（点对点相乘）得出一个新的序列

• 我们再对此序列做一次逆傅里叶变换，也就是将序列变为 $$<g(\omega^{-0}), g(\omega^{-1}),g(\omega^{-2}), \dots, g(\omega^{-(n-1)})>$$

也就是 $$<n \times c_0, n \times c_1, \dots, n \times c_{n-1}>$$

最后对于每一项除以 $$n$$ 即可。

## 快速傅里叶变换（FFT）

FFT -> Fast-Fast-TLE

\begin{aligned} f_{even}(x) &= c_0 + c_2 x + c_4 x^2 + \dots + c_{n-2} x^{\frac n2 - 1} \\ f_{odd}(x) &= c_1 + c_3 x + c_5 x^2 + \dots + c_{n-1} x^{\frac n2 - 1} \\ \end{aligned}

$f(x) = f_{even}(x^2) + x f_{odd}(x^2)$

\begin{aligned} f(\omega_n^k) &= f_{even}(\omega_{n}^{2k}) + \omega_n^kf_{odd}(\omega_n^{2k}) \\ f(\omega_n^{k+\frac n2}) &= f_{even}(\omega_n^{2k+n}) + \omega_n^{k+\frac n2}f_{odd}(\omega_n^{2k+n}) \end{aligned}

$\omega_n^{\frac n2} = -1$

\begin{aligned} f(\omega_n^k) &= f_{even}(\omega_{\frac n2}^k) + \omega_n^kf_{odd}(\omega_{\frac n2}^k) \\ f(\omega_n^{k+\frac n2}) &= f_{even}(\omega_{\frac n2}^k) - \omega_n^kf_{odd}(\omega_{\frac n2}^k) \\ \end{aligned}

#include <iostream>
#include <algorithm>
#include <vector>

struct Complex {
double real, imag;
Complex() : real(0), imag(0) {}
Complex(double re, double im) : real(re), imag(im) {};
inline Complex operator + (const Complex & b) { return Complex(real + b.real, imag + b.imag); }
inline Complex operator - (const Complex & b) { return Complex(real - b.real, imag - b.imag); }
inline Complex operator * (const Complex & b) { return Complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real); }
};
typedef std::vector<Complex> Vector;

const double PI = acos(-1);

void FFT(Vector &v, int n, int inv) {
if (n == 1) return; // 递归边界，只有一个元素，不做变换

// 奇偶变化为两个向量
int mid = n >> 1;
Vector even(mid), odd(mid);
for (int i(0); i < n; i += 2) {
even[i >> 1] = v[i], odd[i >> 1] = v[i + 1];
}
// 递归操作
FFT(even, mid, inv), FFT(odd, mid, inv);

// 进行合并操作
// 定义基本 omega
Complex omega(cos(PI * 2 / n), inv * sin(PI * 2 / n));
// 当前旋转因子
Complex w(1, 0);
for (int i(0); i < mid; ++i, w = w * omega) {
v[i] = even[i] + w * odd[i];
v[i + mid] = even[i] - w * odd[i];
}
}

int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(0), std::cout.tie(0);

int n, m;
std::cin >> n >> m;

// 获取最终的长度，必须是 2 的次幂，且比两个向量卷起来要长
int O(1);
while (O <= m + n) O <<= 1;
// std::cout << PI << " " << O << std::endl;

Vector A(O), B(O);
for (int i(0); i <= n; ++i) std::cin >> A[i].real;
for (int i(0); i <= m; ++i) std::cin >> B[i].real;

FFT(A, O, 1);
FFT(B, O, 1);

// 我们单点相乘，然后进行逆变换，求出每一项的系数
for (int i(0); i < O; ++i) A[i] = A[i] * B[i];
FFT(A, O, -1);

// 最后进行输出
// 记得两个东西卷起来之后是 n + m 次的
for (int i(0); i <= n + m; ++i) {
// 这里是向上取整？
std::cout << (long long)(A[i].real / O + 0.5) << ' ';
} std::cout << std::flush;
return 0;
}


### 蝶形优化

for (int i(0); i < mid; ++i, w = w * omega) {
v[i] = even[i] + w * odd[i];
v[i + mid] = even[i] - w * odd[i];
}


$$rev[(abcd)_2] = (dcba)_2 = (rev[(0abc)_2]>>1)\ or\ ((d\&1)<<3)$$

dp[x] = (dp[x>>1] >> 1) | ((x&1) << (log2(n) - 1))

#include <complex>
#include <iostream>
#include <algorithm>
#include <vector>

struct Complex {
double real, imag;
Complex() : real(0), imag(0) {}
Complex(double re, double im) : real(re), imag(im) {};
inline Complex operator + (const Complex & b) { return Complex(real + b.real, imag + b.imag); }
inline Complex operator - (const Complex & b) { return Complex(real - b.real, imag - b.imag); }
inline Complex operator * (const Complex & b) { return Complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real); }
};

typedef std::vector<Complex> Vector;

const double PI = acos(-1);
int O(1), logO(0);

void FFT(std::vector<int> &rev, Vector &v, int inv) {
for (int i(0); i < O; ++i) {
if (i < rev[i]) std::swap(v[i], v[rev[i]]);
}

// 第 log(k) 次合并，一共logO次
// 合并之后区间的长度为 k
for (int k(1); k < O; k <<= 1) {
Complex omega(cos(PI / k), inv * sin(PI / k));
for (int i(0); i < O; i += (k<<1)) { // 处理行
Complex w(1, 0);
for (int j = 0; j < k; ++j, w = w * omega) {
Complex s = v[i + j], t = v[i + j + k] * w;
v[i + j] = s + t, v[i + j + k] = s - t;
}
}
}

if (inv == -1) for (int i(0); i < O; ++i) v[i].real /= O;
}

int main() {
int n, m;
std::cin >> n >> m;

while (O <= n + m) O <<= 1, ++logO;

Vector A(O), B(O);
for (int i(0); i <= n; ++i) std::cin >> A[i].real;
for (int i(0); i <= m; ++i) std::cin >> B[i].real;

std::vector<int> rev(O);
for (int i(0); i < O; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (logO - 1));

FFT(rev, A, 1), FFT(rev, B, 1);

for (int i(0); i < O; ++i) A[i] = A[i] * B[i];
FFT(rev, A, -1);

for (int i(0); i <= n + m; ++i) {
std::cout << (long long)(A[i].real + 0.1) << ' '; // 向上取整
} std::cout << std::flush;
return 0;
}


## 快速数论变换（NTT）

\begin{aligned} \forall x \in [1, \varphi(p)), g^x \not\equiv 1 \pmod p \\ g^{\varphi(p)} \equiv 1 \pmod p \end{aligned}

$$v = \varphi(p)$$，考虑：

$w_n = g^{\frac vn}$

• 消去$$w_{dn}^{dk} = w_n^k$$

证明：考虑展开即可：$$w_{dn}^{dk} = g^{\frac v{dn} dk} = g^{\frac nk}$$

• 折半$$(w_n^{k+\frac n2})^2 = w_n^{2k + n} = w_n^{2k} = w_{\frac n2}^k$$

• 求和$$\sum_{i=0}^{n-1} (\omega_n^k)^i = 0$$

证明：也考虑等比数列：

\begin{aligned} \sum_{i=0}^{n-1} (\omega_n^k)^i &= \cfrac {(w_n^k)^n - 1}{w_n^k - 1} \\ &= \cfrac {(w_n^n)^k - 1}{w_n^k - 1} \\ &= \cfrac {1^k - 1}{w_n^k - 1} = 0 \end{aligned}

这过程怎么一模一样……

const int MOD = 998244353, g = 3, ig = 332748118;
int O(1), logO(0);

typedef vector<long long> vec;
typedef vector<int> ivec;

void NTT(ivec &rev, vec &v, int inv) {
for (int i(0); i < O; ++i) {
if (rev[i] < i) std::swap(v[i], v[rev[i]]);
}

// 第 log(k) 次合并
for (int k(1); k < O; k <<= 1) {
long long omega = qpow((long long)(~inv ? g : ig), (MOD - 1) / (k << 1), MOD);
for (int i(0); i < O; i += (k << 1)) {
long long w(1);
for (int j(0); j < k; ++j, w = w * omega % MOD) {
long long s = v[i + j], t = w * v[i + j + k];
v[i + j] = (s + t) % MOD, v[i + j + k] = ((s - t) % MOD + MOD) % MOD;
}
}
}

if (inv == -1) {
for (int iv(qpow((long long)O, MOD - 2, MOD)), i(0); i < O; ++i) {
v[i] = v[i] * iv % MOD;
}
}
}


posted @ 2023-02-10 15:28  jeefy  阅读(160)  评论(1编辑  收藏  举报