学了好久,终于基本弄明白了
推荐两个博客:
戳我
戳我
再推荐几本书:
《ACM/ICPC算法基础训练教程》
《组合数学》(清华大学出版社)
《高中数学选修》
预备知识
复数方面
找数学老师去
\[i^{2}=-1,i为虚数的单位\]
坐标系上纵轴就是虚数轴,复数就是这上面的点
三种表示法:
\[一般:a + bi,a为实部,b为虚部\]
\[指数:e^{i\theta}*坐标系上的模长\]
\[三角:模长*(cos\theta + i sin \theta)\]
运算:
加减法:实部虚部分别相加
乘法:\[(a + bi) * (c + di) = ac + adi + bci + bdi^{2}
= ac-bd+(ad+bc)i\]
欧拉公式
\[e^{ix} = cosx + isinx(就是指数表示和三角表示)\]
\[特别的e^{i\pi} = -1\]
多项式
\[系数表示法:A(x) = \Sigma _{k=0}^{n - 1} a_kx^k\]
\[点值表示法:对于所有的x_k,求出它们对应的A(x),设为y_k\]
\[则可以用\{(x_0, y_0), (x_1, y_1), ......, (x_n-1, y_n-1)\} 表示这个多项式
并且是唯一确定的\]
单位复数根
\[n次单位复数根\omega^{n} = 1,n次单位复数根刚好有n个对应e^{\frac{2k\pi i}{n}},其中k=0到n - 1\]
三个性质:
消去引理:
\[n, d, k为正整数,则\omega^{dk}_{dn}=\omega^{k}_{n}\]
\[证明:套e^{\frac{2k\pi i}{n}} 即可\]
折半引理:
\[n为大于零的偶数,则(\omega^{k+\frac{n}{2}}_{n})^{2}=\omega^{2k+n}_{n}=\omega^{2k}_{n}\omega^{n}_{n}=(\omega^{k}_{n})^{2}\]
求和引理:
大于1的整数n,和不被n整除的非负整数k,有
\[\Sigma^{n-1}_{j=0}(\omega^{k}_{n})^{j}=0\]
证明可以用等比数列求和公式得到(很简单的,手推一遍就好)
Rader排序
其实就是二进制数位翻转
正题
DFT
对于k=0~n-1,定义:
\[y_k=A(\omega^{k}_{n}) = \Sigma^{n-1}_{j=0} a_j(\omega^{k}_{j})^j\]
\[得到的y称为a的离散傅里叶变换,记作y=DFT_n(a) (这里的y,a指的是所有的y_k, a_k,即向量y,a)\]
逆DFT
\[就是DFT的逆变换,求出向量a,记为DFT^{-1}\]
假设得到了向量y
\[对于y_k = \Sigma^{n-1}_{i=0}a_i(\omega^{k})^i\]
\[有a_k = \frac{1}{n}\Sigma^{n-1}_{i=0}y_i(\omega^{-k})^i\]
\[证明:a_k=\frac{1}{n}\Sigma^{n-1}_{i=0}y_i(\omega^{-k})^i=\frac{1}{n}\Sigma^{n-1}_{i=0}( \Sigma^{n-1}_{j=0}a_j(\omega^{k})^j)(\omega^{-k})^i=\frac{1}{n} \Sigma^{n-1}_{i=0}a_i(\Sigma^{n-1}_{j=0}(\omega^{j-k})^i)\]
\[可以用等比数列求和出上面的就是a_k(当j\ne k是括号里的为0,当j=k时为1)\]
FFT
上面已经把DFT和逆DFT搞定了,两个几乎是一样的
所以求多项式的积(卷积)可以用DFT转换成点值表示,就可以O(n),一一相乘,得到积的多项式的点值表示,最后用逆DFT得到系数表示
复杂度瓶颈在于怎样快速求解DFT(逆DFT和DFT方法一样)
FFT就是一个O(nlogn)求解DFT的方法
首先把A(x)分成奇数项和偶数项记作
\[A^{[0]}(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2} - 1}\]
\[A^{[1]}(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{\frac{n}{2} - 1}\]
\[显然A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) \]
那么
\[A(\omega^k_n)=A^{[0]}((\omega^k_n)^2) + \omega^k_n A^{[1]}((\omega^k_n)^2)=A^{[0]}(\omega^k_{\frac{n}{2}}) + \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}}) \]
\[因为\omega^{\frac{n}{2}}_{n}=\omega_{2}=e^{k\pi i}= cos k\pi + i sin k\pi = -1\]
\[所以A(\omega^{k+\frac{n}{2}}_n)=A^{[0]}(\omega^k_{\frac{n}{2}}) - \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}})\]
这称为蝴蝶操作
于是对每个y值的求解可以通过分组求出,若递归变成处理子任务,这样复杂度就成了O(nlogn)
这样不停地分组,最后就相当于Rader排序了一番,所以也可以变成非递归的
注意每次都要把多项式补成2的幂,便于FFT
递归写可能好理解一些,但不好写
还有一些东西什么的,其实记一记就好了其实自己说不清
系统的复数complex代码
# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(3e6 + 10);
const double Pi = acos(-1);
IL ll Read(){
char c = '%'; ll x = 0, z = 1;
for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
return x * z;
}
int n, m, r[_], l;
complex <double> a[_], b[_];
IL void FFT(complex <double> *P, int opt){
for(RG int i = 0; i < n; ++i) if(i < r[i]) swap(P[i], P[r[i]]); //Rader排序
for(RG int i = 1; i < n; i <<= 1){
complex <double> W(cos(Pi / i), opt * sin(Pi / i)); //旋转因子
for(RG int p = i << 1, j = 0; j < n; j += p){
complex <double> w(1, 0);
for(RG int k = 0; k < i; ++k, w *= W){
complex <double> X = P[j + k], Y = w * P[j + k + i];
P[j + k] = X + Y; P[j + k + i] = X - Y; //蝴蝶操作
}
}
}
}
int main(RG int argc, RG char *argv[]){
n = Read(); m = Read();
for(RG int i = 0; i <= n; ++i) a[i] = Read();
for(RG int i = 0; i <= m; ++i) b[i] = Read();
m += n;
for(n = 1; n <= m; n <<= 1) ++l;//补成2的幂
for(RG int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));//Rader排序预处理
FFT(a, 1); FFT(b, 1); //DFT
for(RG int i = 0; i < n; ++i) a[i] = a[i] * b[i]; //点值直接相乘
FFT(a, -1); //逆DFT
for(RG int i = 0; i <= m; ++i) printf("%d ", (int)(a[i].real() / n + 0.5));
return 0;
}
或者可以自己定义complex,用复数运算
struct Complex{
double real, image;
IL Complex(){ real = image = 0; }
IL Complex(RG double a, RG double b){ real = a; image = b; }
IL Complex operator +(RG Complex B){ return Complex(real + B.real, image + B.image); }
IL Complex operator -(RG Complex B){ return Complex(real - B.real, image - B.image); }
IL Complex operator *(RG Complex B){ return Complex(real * B.real - image * B.image, real * B.image + image * B.real); }
}