多项式与生成函数

1. FFT

本文暂时只讨论其在 $\text{OI}$ 领域，即加速多项式乘法的应用。

---

如果我们想要求出两个多项式之积，朴素的复杂度是 $O(n^2)$ 的。

即若 $F(x)=a_0+a_1x+…+a_Nx^n,G(x)=b_0+b_1x+…+b_mx^m$，

那么 $F(x)\cdot G(x)=\sum\limits_{i=0}^{n+m} (\sum\limits_j a_jb_{i-j}) x^i$。

---

就像空间中的点可以根据选定的基底映射到唯一的坐标上、坐标又可以映射到唯一的点上，如果我们将多项式映射到一个便于运算的“坐标”上，再将运算后的“坐标”转化回多项式，或许就可以优化这个运算过程。

• 定理：$n+1$ 个 $x$ 互不相同的点可以唯一确定一个 $n$ 次多项式。

于是我们可以将一个多项式映射到有限点集上。鉴于一些优秀的性质，我们选取的 $x$ 是复平面上的 $N(N\geq n+1)$ 个单位根 $\lbrace w_N^k\rbrace $。

• 单位根是 $z^N=1$ 的解，在复平面中均匀分布在半径为 $1$ 的圆上，根据欧拉公式有 $w_N^k=e^{i\frac{2\pi}{N}k}=\cos(\frac{2\pi}{N}k)+i\cdot\sin(\frac{2\pi}{N}k)$，其中 $i$ 是虚数单位。或者你可以理解为在单位圆的圆弧上均匀的选取 $N$ 个点。

（优秀的性质：）

首先，$w_{2N}^{2k}=w_N^k$ 是显然的。

当 $N=2^x$ 时，还能得到：

• $w_N^k=-w_N^{k+N/2}$，从而 $(w_N^k)^2=(w_N^{k+N/2})^2$。

• $(w_N^k)^2=w_N^{2k}=w_{N/2}^k$。

---

快速傅里叶变换

首先考虑怎么从多项式转化成点集。这是显然的，我们只需要求出所有 $F(w_N^k)$ 即可。

但是，朴素计算的复杂度仍然是 $O(n^2)$ 的，考虑优化。

• 利用分治的思想，如果将原函数按照次数奇偶分类：

$F(x)=a_0+a_1x+a_2x^2+a_3x^3+…$

$=(a_0+a_2x^2+…)+x(a_1+a_3x^2+…)$

$=G(x^2)+x\cdot H(x^2)$

其中 $G(x)=a_0+a_2x+…,H(x)=a_1+a_3x+…$

• 代入 $x=w_N^k$ 后：

$F(w_N^k)=G((w_N^k)^2)+w_N^k\cdot H((w_N^k)^2)$

$=G(w_{N/2}^k)+w_N^k\cdot H(w_{N/2}^k)$

• 如果代入 $x=w_N^{k+n/2}$：

$F(w_N^{k+n/2})=G(w_{N/2}^k)-w_N^k\cdot H(w_{N/2}^k)$

于是，只需要算出 $G(w_{N/2}^k)$ 和 $H(w_{N/2}^k)$，就可以知道 $F(w_N^k)$ 和 $F(w_{N}^{k+n/2})$，对 $G$ 和 $H$ 分别递归即可。

---

在此基础上，还可以进一步优化。

位逆序置换

手动模拟不断将奇偶系数拆分的过程，例如 $N=8$ 时：

• $\lbrace a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\rbrace $

• $\lbrace a_0,a_2,a_4,a_6\rbrace ,\lbrace a_1,a_3,a_5,a_7\rbrace $

• $\lbrace a_0,a_4\rbrace ,\lbrace a_2,a_6\rbrace ,\lbrace a_1,a_5\rbrace ,\lbrace a_3,a_7\rbrace $

结论是，每个位置的下标，二进制下翻转得到的数，就是系数初始位置的下标。这个变换称为位逆序置换。

具体地，当我们从前向后遍历系数下标 $i$ 时，可以同时维护一个由 $i$ 翻转的下标 $j$。每次 $i:=i+1$，而 $j$ 不断减去高位连续的 $1$，直到遇到 $0$ 后在此 $+1$。

for(int i = 1, j = len >> 1, k; i < len && i < j; ++ i ){
    swap(a[i], a[j]);
    k = len >> 1;
    while(j >= k){
        j -= k;
        k >>= 1;
    }
    j += k;
}

---

蝶形运算优化

位逆序置换后，$G(w_{N/2}^k)$ 存在位置 $k$ 上，$H(w_{N/2}^k)$ 存在位置 $k+n/2$ 上。而 $F(w_N^k)$ 需要存到位置 $k$ 上，$F(w_N^{k+n/2})$ 需要存到 $k+n/2$ 上，于是不需要额外开空间，直接覆写即可。

---

快速傅里叶逆变换

已知 $F(w_N^k)=\sum_i a_i\cdot w_N^{k\cdot i}$，下证 $a_i=\frac{1}{n}\sum_k F(w_N^k)\cdot w_N^{-k\cdot i}$。

• 左式代入右式得：

$a_i=\sum_k w_N^{-k\cdot i} \cdot (\sum_j a_j\cdot w_N^{k\cdot j})$

$=\sum_k\sum_j a_j\cdot w_N^{k\cdot(j-i)}$

$=\sum_j a_j \sum_k (w_N^{j-i})^k$

设 $S(w_N^a)=\sum_i (w_N^a)^i$。

• 当 $a\equiv0\mod N$ 时，$S(w_N^a)=n$。

• 当 $a\not\equiv0\mod N$ 时，错位相减解得：

$S(w_N^a)=\frac{(w_N^a)^n-(w_N^a)^0}{w_N^a-1}=0$。

所以 $a_i=\frac{1}{n}\sum a_i$。

于是对单位根取倒数，对函数值跑一遍 $\text{FFT}$，再除以 $N$ 即得到系数。

---

Code

洛谷 P3803 FFT 模板

可以在此基础上增加预处理原根等优化。

#include<iostream>
#include<algorithm>
#include<string.h>
#include<complex>
#include<vector>
#include<queue>
#include<set>
#include<map>

using namespace std;

typedef long long ll;
typedef pair<int, int> pii;

const int N = 3e6 + 9;
int n, m;
complex<double> F[N], G[N];

inline int read(){
	int x = 0, f = 1;
	char ch = getchar();
	while(ch < '0' || ch > '9'){
		if(ch == '-') f = -1;
		ch = getchar();
	}
	while(ch >= '0' && ch <= '9'){
		x = x * 10 + (ch - '0');
		ch = getchar();
	}
	return x * f;
}

//---

void DFT(complex<double>* A, int lim, int op){
    for(int i = 1, j = lim >> 1, k; i < lim; ++ i ){
        if(i < j) swap(A[i], A[j]);

        k = lim >> 1;
        while(j >= k && k > 0){
            j -= k;
            k >>= 1;
        }
        j += k;
    }

    for(int len = 2; len <= lim; len <<= 1 ){
        double theta = op * acos(-1.0) * 2 / len;
        for(int i = 0; i < lim; i += len ){
            for(int j = 0; j < (len >> 1); ++ j ){
                complex<double> w(cos(theta * j), sin(theta * j));
                complex<double> u = A[i + j];
                complex<double> v = A[i + j + (len >> 1)] * w;
                A[i + j] = u + v;
                A[i + j + (len >> 1)] = u - v;
            }
        }
    }
}

int main(){
    /* freopen("a.in", "r", stdin);
    freopen("a.out", "w", stdout); */

    scanf("%d %d", &n, &m);
    for(int i = 0; i <= n; ++ i )
        F[i].real(read());
    for(int i = 0; i <= m; ++ i )
        G[i].real(read());
    
    int x = 1;
    while(x <= n + m) x <<= 1;

    DFT(F, x, 1), DFT(G, x, 1);
    for(int i = 0; i < x; ++ i ) F[i] *= G[i];
    DFT(F, x, -1);

    for(int i = 0; i <= n + m; ++ i )
        printf("%lld ", (ll)round(F[i].real() / x));
    return 0;
}