【算法】FFT从敲门到跑路

请注意,本文编写于 246 天前,最后修改于 114 天前,其中某些信息可能已经过时。

在此特感谢function2 深入浅出的讲解,本文出自其讲义

多项式的数域

对于数域$F$,若$a_0 , a_1 , a_2 , \cdots , a_n \in F$

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

为数域$F$上的一个多项式

其数域$F$上的所有多项式记作$F[x]$

同样的

$Q[x]$为有理数域的所有多项式

$R[x]$为实数域的所有多项式

$C[x]$为复数域的所有多项式

多项式的度数

多项式最高非零项的次数成为多项式的度数,用$deg$表示

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

其中$a_i\neq 0$

记作$degf(x)=n$

首一多项式(略)

加法与数乘运算(略)

多项式的线性空间

线性空间,即向量空间

满足线性性,封闭性

线性性: 量与量之间按比例呈线性的关系

封闭性: 其两个集合内元素进行二元运算结果必然在集合内

显然多项式加法,数乘运算满足线性性

记$F_n[x]$表示数域$F$上度数不超过$n$的全体多项式的集合,显然它们构成一个线性空间

多项式与向量满足一一对应的双射关系

$$f(x)=\sum_{i = 0} ^ {n}a_ix^i \Leftrightarrow \vec{a}=(a_0 , a_1 , \cdots , a_n)$$

多项式的系数表达

用一个$n+1$维的向量来表示$F_n[x]$中的全体多项式

多项式乘法

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

$$g(x)=\sum_{i=0}^{n}b_ix_i$$

$$h(x) = f(x)*g(x)=\sum_{i=0}^{2*n}\sum_{j=0}^{i}a_jb_{i-j}x^i$$

向量的卷积

多项式乘法中将多项式用向量表示其系数关系就是向量的卷积
$$f(x)=\vec{a}$$

$$g(x)=\vec{b}$$

$$h(x)=\vec{c}=\vec{a}\otimes\vec{b}$$

$$c_i=\sum_{j=0}^{i}a_j*b_{i-j}$$

多项式的幂

$$f^{k}(x)=\prod_{i=1}^{k}f(x)$$

多项式的复合

$$(f \circ g)(x)=f(g(x))=\sum_{i=0}^{n}a_ig^{i}(x)$$

多项式的分治乘法

设$n$为$2$的自然数次幂,即$\exists k\in N,n=2^k$

设$f,g\in F_{2n-1}[x]$(一共有2n项)

将其分为两半

$$f(x)=f_0(x)+x^nf_1(x)$$

$$g(x)=g_0(x)+x^ng_1(x)$$

显然的,$f_0,f_1,g_0,g_1\in F_{n-1}[x]$

$$(f \times g)(x) = (f_0 \times g_0)(x) +x_n(f_0 \times g_1 +f_1 \times g_0)(x) + x_{2n}(f_1 \times g_1)(x)$$

观察中间项,$(f_0 \times g_1 +f_1 \times g_0)(x)=((f_0+f_1)\times(g_0 + g_1) - f_0 * g_0 - f_1 * g_1)(x)$

那么我们就只需要处理$(f_0 \times g_0)(x) , (f_1 \times g_1)(x) , ((f_0 +f_1) \times (g_0 + g_1))(x)$,规模均为原问题的一半

所以$T(n)=3T(\frac{n}{2})+O(n)=O(n^{log_23})\approx O(n^{1.585})$

多项式的点值表达

$$degf(x) = n$$
$$x_0 , x_1 , x_2 , \cdots , x_n \in F$$
$$y_0 = f(x_0) , y_1 = f(x_1) , y_2 = f(x_2) , \cdots , y_n = f(x_n)$$
则$(x_0 , y_0) , (x_1 , y_1) , (x_2 , y_2) , \cdots ,(x_n , y_n)$为多项式$f(x)$的点值表达

点值表达与系数表达满足双射关系

有范德蒙德矩阵$V$

$$ \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{pmatrix} \begin{pmatrix} a_0\\ a_1\\ \vdots\\ a_n \end{pmatrix} = \begin{pmatrix} y_0\\ y_1\\ \vdots\\ y_n \end{pmatrix} $$

范德蒙德矩阵的行列式为$\prod_{j < i} (x_i - x_j)$

由于$x_i \neq x_j$,所以行列式值非零,有唯一解

在$x$确定时,$a,y$有一一对应的双射关系

多项式的点值乘法(略)

多项式的插值

给出多项式的点值表达,求多项式的系数表达

  1. 高斯消元$O(n^3)$
  2. 拉格朗日插值法$O(n^2)$

高斯消元(略)

拉格朗日插值法

$$L(x)=\sum_{i=0}^{n}\frac{\prod_{j \neq i}x - x_j}{\prod_{j \neq i}x_i - x_j}y_i$$
不难发现$L(x_i) = y_i$

所以$L,f$有相同的系数表达

自然数幂级数和

自然数幂级数和记作$S$
$$S(n,k)=\sum_{i = 0} ^ {n}i^k$$
$$S(n,0) = n + 1$$

在这里我们将$0^0$看为$1$

$$S(n,1) = \frac{n(n+1)}{2}$$
$$S(n,2) = \frac{n(n+1)(2n+1)}{6}$$

合理推测

$S(n,k)$为$n$的$k+1$次多项式

为了计算$S(n,k)$

选择$x_i=i(i = 0,\cdots , k),y_i = S(x_i , k)$

插值求解即可

证明

$$ (n+1)^{k+1} - n^{k+1} = \sum_{i=0}^{k} \begin{pmatrix} k+1\\ i \end{pmatrix} n^i $$

$$ n^{k+1} - (n - 1) ^ {k+1} = \sum_{i=0}^{k} \begin{pmatrix} k+1\\ i \end{pmatrix} (n - 1)^i $$

$$\vdots$$

$$ 1^{k+1}-0^{k+1}=\sum_{i=0}^{k} \begin{pmatrix} k+1\\ i \end{pmatrix} 0^i $$

累加得

$$ (n+1)^{k+1}=\sum_{i=0}^{k} \begin{pmatrix} k+1\\ i \end{pmatrix} S(n,i) $$

显然可得$S(n,k)$的度数

欧拉公式

$$e^{i\theta}=cos(\theta)+isin(\theta)$$
几何意义上是单位圆上角度为$\theta$的点

Circle.png
Circle.png

单位复数根

考虑构造一个$x$,使得$x^n=1$

由欧拉公式可得$e^{2k\pi i}=1$

则$x^n = e^{2k\pi i} , x = e^{\frac{2k\pi i}{n}}$

当$k = 1$ 时,令$\omega_n=e^{\frac{2\pi i}{n}}$,称为主$n$次单位根

$w_n^k = e^{\frac{2k\pi i}{n}} = cos(\frac{2k\pi}{n}) + isin(\frac{2k\pi}{n})$

由三角函数的周期性$T=2\pi$,不难得到$w_n^k$的周期性

$w_n^k = w_n^{k\mod n}$

$w_n^n = w_n^0 = 1$

当$k=0,1,2,\cdots,n-1$时,得到$w_n^0,w_n^1,w_n^2,\cdots,w_n^{n-1}$,称为$n$次单位复数根

单位复数根的性质

消去引理:对于常数$d>0$,有$w_{dn}^{dk} = w_n^k$
折半引理:若$n$为偶数,则$n$次单位复数根的平方的集合就是$\frac{n}{2}$次单位复数根的集合。特别的,每个$\frac{n}{2}$次单位复数根出现两次
求和引理:

$$ \sum_{i=0}^{n-1}(w_n^k)^i= \begin{cases} n &n\mid k\\ 0 &n\nmid k \end{cases} $$

离散傅里叶变换(DFT)

多项式的系数表达式$f=(a_0,a_1,a_2,\cdots,a_{n-1})$,求出在$x_0 = w_n^0 , x_1 = w_n^1 , x_2 = w_n^2 , \cdots , x_{n-1} = w_n^{n-1}$初的点值

$$y_k = f(x_k) = f(w_n^k)$$

输出向量$y=(y_0 , y_1 , y_2 , \cdots , y_{n-1})$

代入范德蒙德矩阵得:

$$ \begin{pmatrix} 1 & w_n^{0} & {w_n^{0}}^2 & \cdots & {w_n^{0}}^{n-1}\\ 1 & w_n^{1} & {w_n^{1}}^2 & \cdots & {w_n^{1}}^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{(n-1)} & {w_n^{(n-1)}}^2 & \cdots & {w_n^{(n-1)}}^{n-1}\\ \end{pmatrix} \begin{pmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{pmatrix} $$

此过程称为长度为$n$的离散傅里叶变换,记作$y = DFT_n(f)$

注意: 多项式的系数表达式的度数为$n-1$,而单位根为$n$次单位复数根

快速傅里叶变换(FFT)

$f = (a_0 , a_1 , a_ 2 , \cdots , a_{n-1})$

将奇数项,偶数项分别提出,得到构造出两个度为$\frac{n}{2}$的系数表示(假设$n$为偶数)

$f_0 = (a_0 , a_2 , a_4 , \cdots , a_{n-2})$

$f_1 = (a_1 , a_3 , a_5 , \cdots , a_{n-1})$

那么$f(x) = f_0(x^2) + x*f(1)(x^2)$

设$k < \frac{n}{2}$

$$f(w_n^k) = f_0(w_2^{2k}) + w_n^k f_1(w_n^{2k}) = f_0(w_{\frac{n}{2}}^k) + w_n^k f_1(w_{\frac{n}{2}}^k)$$

$$f(w_n^{k+\frac{n}{2}}) = f_0(w_n^{2k+n}) + w_n^{k+\frac{n}{2}} f_1(w_n^{2k+n}) = f_0(w_{\frac{n}{2}}^k) - w_n^k f_1(w_{\frac{n}{2}}^k$$

令$y_k = f(w_n^k) , y_k^{[0]} = f_0(w_{\frac{n}{2}}^k) , y_k^{[1]} = f_1(w_{\frac{n}{2}}^k)$

那么

$$y_k = y_k^{[0]} + w_n^ky_k^{[1]}$$

$$y_{k+\frac{n}{2}} = y_k^{[0]} - w_n^ky_k^{[1]}$$

这一系列操作称为蝴蝶操作,$w_n^k$称为旋转因子

每次递归求得$y_k^{[0]} , y_k^{[1]}$,问题规模变为$\frac{n}{2})$,合并得到$y_k$

$T(n)=2T(\frac{n}{2}) + O(n) = O(n\log n)$

逆离散傅里叶变换(IDFT)

多项式的点值表达式$(w_n^0 , y_0) , (w_n^1 , y_1) , (w_n^2 , y_2) , \cdots , (w_n^{n-1} , y_{n-1})$,转化为系数表达式

尝试求得范德蒙德矩阵的逆矩阵$V^{-1}$

即求$V_n V_n^{-1} = I_n$

即$v_i^{T}v_j^{-1} = \begin{cases} 1 &i=j\\0 &i\neq j\end{cases}$,其中$v_i^{T}$为$V$的行向量,$T$为转置,$v_j^{-1}$为$V^{-1}$的列向量

$v_i^{T} = (w_n^{0i} , w_n^{1i} , w_n^{2i} , \cdots , w_n^{(n-1)i})$

联系求和定理

$$ \sum_{i=0}^{n-1}(w_n^k)^i= \begin{cases} n &k = 0\\ 0 &k \neq 0 \end{cases} $$

那么构造$v_j^{-1} = \frac{1}{n}(w_n^{-0j} , w_n^{-1j} , w_n^{-2j} , \cdots , w_n^{-(n-1)j})^{T}$

验证:

$$v_i^{T}v_j^{-1} = \frac{1}{n}\sum_{k = 0}^{n - 1}(w_n^{i-j})^k = \begin{cases} n &i=j\\ 0 &i\neq j\end{cases}$$

将$Va = y$同乘$V^{-1}$得$a = yV^{-1}$,即:

$$ \frac{1}{n} \begin{pmatrix} 1 & w_n^{-0} & {w_n^{-0}}^2 & \cdots & {w_n^{-0}}^{n-1}\\ 1 & w_n^{-1} & {w_n^{-1}}^2 & \cdots & {w_n^{-1}}^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{-(n-1)} & {w_n^{-(n-1)}}^2 & \cdots & {w_n^{-(n-1)}}^{n-1}\\ \end{pmatrix} \begin{pmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{pmatrix} = \begin{pmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{pmatrix} $$

可见,$a$是系数表达$y = (y_0 , y_1 , y_2 , \cdots , y_{n-1})$的多项式在$w_n^0 , w_n^1 , w_n^2 , \cdots , w_n^{n-1}$处取的的点值除以$n$

在这些位置求出点值的过程被称为逆离散傅里叶变换,记作$DFT_n^{-1}(y)$,只需要对前面的DFT稍加修改即可实现

卷积定理

设$a,b$为$n$维向量,$a,b$的卷积$c = a\otimes b$满足:

$$c_i = \sum_{j = 0} ^{i} a_j * b_{i-j}$$

通过多项式乘法与其系数表达卷积的关系可以得到:

$$c = DFT_{2n}^{-1}(DFT_{2n}(a) * DFT_{2n}(b))$$

FFT的迭代实现

FFT.png
FFT.png

对于$n = 2^k$,考虑第$i$次分治(第$i$层)

二进制第$i$位为$0$,放入左侧,最终位置的第$k-i$位为$0$

二进制第$i$位为$1$,放入右侧,最终位置的第$k-i$位为$1$

不难发现,最终位置的二进制为原来位置的逆序

FFT题目:

BZOJ 3160 万径人踪灭

BZOJ 3527

BZOJ 4259 残缺的字符串

HDU 4609 3-idiots

codechef Prime Distance On Tree

板子

#include <bits/stdc++.h>
using namespace std;
#define N 3000010
const double pi=acos(-1);
struct Complex {
    double x;double y;
    Complex (){};
    Complex (double _x,double _y){x=_x;y=_y;}
    Complex operator + (const Complex o) {return Complex(x+o.x,y+o.y);}
    Complex operator - (const Complex o) {return Complex(x-o.x,y-o.y);}
    Complex operator * (const Complex o) {return Complex(x*o.x-y*o.y,x*o.y+y*o.x);}
}a[N],b[N];
int n,m,l;
int r[N],c[N];
char s[60010];
int read() {
    int ans=0,flag=1;
    char ch=getchar();
    while((ch>'9' || ch<'0') && ch!='-') ch=getchar();
    if(ch=='-') flag=-1,ch=getchar();
    while(ch>='0' && ch<='9') ans=ans*10+ch-'0',ch=getchar();
    return ans*flag;
}
void fft(Complex *now,int f) {
    for(int i=0;i<n;i++)
        if(i<r[i]) swap(now[i],now[r[i]]);
    for(int i=1;i<n;i<<=1) {
        Complex wn(cos(pi/i),f*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)) {
            Complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn) {
                Complex x=now[j+k],y=w*now[j+k+i];
                now[j+k]=x+y;now[j+k+i]=x-y;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<n;i++)
            now[i].x/=n;
}
int main() {
    n=read();m=read();
    for(int i=0;i<=n;i++) a[i].x=read();
    for(int i=0;i<=m;i++) b[i].x=read();
    m+=n;
    for(n=1;n<=m;n<<=1) l++;
    for(int i=0;i<n;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    fft(a,1);fft(b,1);
    for(int i=0;i<=n;i++) a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=m;i++)
        printf("%d ",(int)(a[i].x+0.5));
    return 0;
}
Comments

添加新评论