# Convolution

$\boldsymbol{h(x) = \int _{- \infty} ^{\infty}g(\tau) \cdot f(x - \tau)} \rm{d\tau}$

$A(x) \cdot B(x) =\sum\limits_{i = 0}^{n} \sum\limits_{j=0}^{i}{a_jb_{i-j}}$

# Dot Method

$f(x)\Longleftrightarrow{(x_0,y_0),(x_1,y_1),(x_2,y_2)....(x_n,y_n)}\\ \forall k,y_k = f(x_k)$

## Advanced Trick Point $\color{red}{1}$ Multiplication

$A(x) = {(x_0,y_0),(x_1,y_1)....(x_{2n},y_{2n})}\\ B(x) = {(x_0,y_0'),(x_1,y_1')....(x_{2n},y_{2n})}$

$A(x)B(x) = {(x_0,y_0y_0'),(x_1,y_1y_1')\cdots(x_{2n},y_{2n}y_{2n}')}$

$\longrightarrow$ 点乘法(时间复杂度 $\Theta(n)$ )

$$\longrightarrow$$ 将得出来的 $C(x)$ 的点值表达式再转换成系数表达式(复杂度 $\Theta(n^2)$ )。

# Base of Optimization

## Advanced Trick Point $\color{red}{2}$ Unit Complex Root

$n$ 次单位复根是满足 $\omega^n = 1$ 的复数 $\omega$ ，其中我们可以由复数的运算法则(辐角相加，模相乘)很简单地得出 $n$ 次单位根有 $n$ 个 这个结论——亦或者是用代数基本定理证，都可以。

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

$e^{2 \pi i} = 1 = \omega^n \Longleftrightarrow \omega = e^{\frac{2\pi i}{n}}$

### $$\frak{Elimination~Lemma}$$

$\omega_{dn}^{dk} = \omega_n^k$

$\mathcal{Proof.}$

$\omega_{dn}^{dk} = \omega^{\frac{2\pi dk}{dn}} = e^{\frac{2\pi ik}{n}} = \omega_n^k$

### $\frak{Binary~Lemma}$

$\mathcal{Proof.}$

$(\omega _n^k)^2 = \omega^{2k}_n=\omega_{\frac{n}{2}}^k$

$(\omega_n^{k + \frac{n}{2}})^2 = \omega_n^{2k + n} \Longrightarrow \omega_n^{2k} \cdot \omega_n^n \Longrightarrow \omega_n^{2k} = (\omega_n^k)^2$

### $\frak{Sum~Lemma}$

$\sum\limits_{j =0}^{n-1}{(\omega_n^k)^j} = 0$

$\sum\limits_{j = 0}^{n}{x^j} = \frac{x^{j +1} -1}{x -1}$

$\sum\limits_{j =0}^{n-1}{(\omega_n^k)^j} = \frac{(\omega_n^k)^n -1}{\omega_n^k -1} \Longrightarrow \frac{(\omega_n^n)^k -1}{\omega_n^k -1} = \frac{(1)^k -1}{\omega_n^k -1}$

# DFT $$\to$$ FFT

## DFT

$A(x) = \sum\limits_{i =0}^{n - 1}{a_ix^i}$

$\omega_n^0,\omega_n^1,\omega_n^2 \cdots \omega_n^{n-1}$

$y_k = A(\omega_n^k)=\sum\limits_{j =0}^{n -1}{a_j \cdot \omega_n^{kj}}$

$\vec{y} = {y_0, y_1, y_2 \cdots y_{n-1}}$

## FFT 优化 DFT

$A(x) = a_0+a_1x+ a_2x^2 \cdots +a_{n-1}x^{n-1}$

$A^{[0]}(x) = a_0 + a_2x+a_4x^2 \cdots +a_{n-2}x^{\frac{n}{2} - 1} \\ A^{[1]}(x) = a_1 + a_3x+a_5x^2 \cdots +a_{n-1}x^{\frac{n}{2} - 1}$

$A(x) = A^{[0]}(x^2)+xA^{[1]}(x^2)$

$\omega_n^0,\omega_n^1,\omega_n^2 \cdots ,\omega_n^{n-1}$

$A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A^{[1]}(\omega_n^{2k+n}) \Longrightarrow A^{[0]}(\omega_n^{2k}\cdot \omega_n^{n})-\omega_n^{k}A^{[1]}(\omega_n^{2k}\cdot \omega_n^{n})$

$A(\omega_n^{k+\frac{n}{2}})=A^{[0]}(\omega_n^{2k})-\omega_n^{k}A^{[1]}(\omega_n^{2k})$

int Lim = 1, N, M ;
function FFT(int lenth, complex *A, int flag){
IF (Lim == 1) return ;
complex A0[lenth >> 1], A1[lenth >> 1] ;//分成两部分
for(int j : 0 to lenth by_grow 2) A0[j >> 1] = A[j], A1[j >> 1] = A[j + 1] ;
FFT(lenth >> 1, A0, flag) ;
FFT(lenth >> 1, A1, flag) ;
complex Wn = unit(,) , w = (1, 0) ;
//Wn是单位根，w用来枚举幂，即我们令主次单位根不变，由于其余单位根都是其整次幂，所以可以用一个w来记录到第几次幂
/*此处求单位根的时候会用到我们的参数flag……嗯没错就用这一次，并且flag的值域为(-1, 1)……是的，只会有两个值*/
for(int j : 0 to (lenth >> 1) by_grow 1 with w = w * Wn){
A[i] = A0[i] + A1[i] * w ;//应用公式，下同
A[i + (lenth >> 1)] = A0[i] - A1[i] * w ; //顺便求出另一半，由折半引理可显然。
}
}
function Main{
input(N), input(M) ;
for(i : 0 to N by_grow 1) => input(A) ;
for(i : 0 to M by_grow 1) => input(B) ;
while(Lim < N + M) Lim <<= 1 ;//Lim为结果多项式的长度（暂时），化为2的幂的便于分治（二分）
FFT(Lim, A, 1) ;//两遍FFT表示从系数化为点值
FFT(Lim, B, 1) ;
for(i : 0 to Lim by_grow 2) => A[i] *= B[i] ;//点乘法，此处需要重定义乘号，因为每一项现在表示的是一个点，有x和y两个属性qwq
}



## FFT 优化 IDFT

$\left|\begin{array}{ccccc}1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_{n} & \omega_{n}^{2} & \cdots & w_{n}^{n-1} \\ 1 & \omega_{n}^{2} & \omega_{n}^{4} & \cdots & \omega_{n}^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_{n}^{n-1} & \omega_{n}^{2(n-1)} & \cdots & \omega_{n}^{(n-1)(n-1)}\end{array}\right| \quad\left|\begin{array}{c}a_{0} \\ a_{1} \\ a_{2} \\ \vdots \\ a_{n-1}\end{array}\right|=\left|\begin{array}{c}y_{0} \\ y_{1} \\ y_{2} \\ \vdots \\ y_{n-1}\end{array}\right|$

$V'V = \sum\limits^ {n-1}_{k=0}{(\dfrac{\omega_n^{-ki}}{n})} \cdot {\omega_n^{kj}} = \frac{1}{n} \sum\limits^ {n-1}_{k=0}{\omega_n^{k(j-i)}}$

$\mathrm{IDFT}_n(y) = \frac{1}{n}\sum\limits_{k = 0}^{n-1}{y_k\omega_n^{-kj}},j\in [0,n-1]$

void FFT(int Lim,complex *A,int flag){
if(Lim == 1) return ;
complex A0[Lim >> 1], A1[Lim >> 1] ;
for(int i = 0; i <= Lim ; i += 2)
A0[i >> 1] = A[i], A1[i >> 1] = A[i+1] ;
FFT(Lim >> 1, A0, flag) ;
FFT(Lim >> 1, A1, flag) ;
complex unit = (complex){cos(2.0 * Pi / Lim) , flag * sin(2.0 * Pi / Lim)}, w = complex(1, 0) ;//欧拉公式
for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
A[i] = A0[i] + w * A1[i] ;
A[i + (Lim>>1)] = A0[i] - w*A1[i];
}
}

int main(){
......................
FFT(A, 1), FFT(B, 1) ;
for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
FFT(A, -1) ;
......................
}


# Iterative Optimization

## Advanced Trick Point $\color{red}{3}$ The Butterfly Operation

$emmm$ 先上一个不是特别卡常数的优化。我们观察之前的代码中，有这么一步：

   for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
a[i] = A0[i] + w * A1[i] ;
a[i + (Lim>>1)] = A0[i] - w*A1[i];
}


   for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
int temp = w * A1[i] ;
a[i] = A0[i] + t ;
a[i + (Lim>>1)] = A0[i] - t ;
}


$qwq$ 这分明是骗小孩子的啦……如果单单这一步就可以卡常数的话，那这个世界该多么美好 $\mathcal{QAQ}$ 。

step 1 成对地取出儿子节点，用蝴蝶操作计算出其 DFT。
step 2 用这一步的 DFT 替换之前的。
step 3 直到我们迭代到根节点为止，否则返回 step 1。

for(j = 1; j < Lim; j <<= 1){//枚举区间长度，从小区间到大区间依次合并。
node T(cos(Pi / j), flag * sin(Pi / j)) ;
for(k = 0; k < Lim; k += (j << 1) ){//两段区间两段区间的枚举，用于合并
node t(1, 0) ;
for(l = 0 ; l < j; l ++, t = t * T){//枚举k所枚举的两个区间内的值，并进行蝴蝶操作。
node Nx = J[k + l], Ny = t * J[k + j + l] ;
J[k + l] = Nx + Ny ; J[k + j + l] = Nx - Ny ;//一次蝴蝶操作
}
}
}


## Trick Point $\color{red}{4}$ The Butterfly Law

for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;


#include <cmath>
#include <cstdio>
#include <iostream>
#define il inline

using namespace std ;
int N, M, K ;
const int MAXN = 3000100 ;
const double Pi = acos(-1.0) ;
int i, j, k, l, Lim = 1, L, R[MAXN] ;
struct node{
double x, y ;
node (double xx = 0, double yy = 0){
x = xx, y = yy ;
}
}A[MAXN], B[MAXN] ;
node operator * (node J, node Q){
return node(J.x * Q.x - J.y * Q.y , J.x * Q.y + J.y * Q.x);
}
node operator + (node J, node Q){
return node(J.x + Q.x , J.y + Q.y);
}
node operator - (node J, node Q){
return node(J.x - Q.x , J.y - Q.y );
}

il int qr(){
int k = 0, f = 1 ;
char c = getchar() ;
while(!isdigit(c)){if(c == '-') f = -1 ;c = getchar() ;}
while(isdigit(c)) k = (k << 1) + (k << 3) + c - 48 ,c = getchar() ;
return k * f ;
}
void FFT(node *J, int flag){
for(i = 0; i < Lim; i ++)
if(i < R[i]) swap(J[i], J[R[i]]) ;//前面的if保证只换一次
for(j = 1; j < Lim; j <<= 1){
node T(cos(Pi / j), flag * sin(Pi / j)) ;
for(k = 0; k < Lim; k += (j << 1) ){
node t(1, 0) ;
for(l = 0 ; l < j; l ++, t = t * T){
node Nx = J[k + l], Ny = t * J[k + j + l] ;
J[k + l] = Nx + Ny ;
J[k + j + l] = Nx - Ny ;
}
}
}
}
int main(){
N = qr(), M = qr() ;
for(i = 0; i <= N; i ++) A[i].x = qr() ;
for(i = 0; i <= M; i ++) B[i].x = qr() ;
while(Lim <= N + M) Lim <<= 1, L ++ ;
for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;
FFT(A, 1), FFT(B, 1) ;
for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
FFT(A, -1) ;
for(i = 0; i <= N + M; i ++)
printf("%d ", (int)(A[i].x / Lim + 0.5)) ;//我们推过的公式里面有一个1/n这一项，最后输出的时候添上即可
return 0 ;
}


# Afterword

1、好多自己当初不理解的地方在代码里就只有半行qaq

2、三个引理中，只有消去引理跟算法的实现没有关系——消去引理主要是用来证明其他引理的

## $\mathfrak{writter:pks}$

posted @ 2018-07-01 20:17  皎月半洒花  阅读(15541)  评论(14编辑  收藏  举报