• FFT


    \[\text{FFT}(快速傅里叶变换) \]

    多项式

    定义:多项式是指由变量、系数以及它们之间的加、减、乘、幂运算(非负整数次方)得到的表达式。

    • :多项式中的每个单项式叫做多项式的项。
    • 多项式的次数:这些单项式中的最高项次数,就是这个多项式的次数。
      一个 \(n\) 次多项式可以表示为 \(\sum\limits_{i=0}^na_ix^i\),其中 \(a_i\) 为系数。

    复数

    定义:我们把形如 \(z=a+bi\)\(a,b\) 均为实数)的数称为 复数,其中 \(a\) 称为实部,\(b\) 称为虚部,\(i\) 称为虚数单位,(\(i^2=-1\)
    \(z_1=a+bi\),\(z_2=c+di\)

    \[z_1\pm z_2=\left(a\pm c\right)+\left(b\pm d\right)i \]

    \[z_1\times z_2=\left(ac-bd\right)+\left(ad+bc\right)i \]

    \[z_1 \div z_2=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i \]

    单位根

    定义:记方程 \(\omega^n=1\left(n\in N^{+}\right)\)\(n\) 个复数解为 \(n\) 次单位根。
    可以解出 \(\omega=\mathrm{e}^{\frac{2k\pi i}{n}}\)
    定义 主次单位根\(\omega_n=\mathrm{e}^{\frac{2\pi i}{n}}\)
    对于 所有单位根,记作:\(\omega_n^k=\mathrm{e}^{\frac{2k\pi i}{n}}\)

    多项式乘法卷积

    对于 \(n\) 次多项式 \(A,B\)\(A=\sum\limits_{i=0}^na_ix^i\),\(B=\sum\limits_{i=0}^nb_ix^i\)
    \(A\)\(B\) 的乘法卷积为 \(C\)\(C=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{n}a_ib_jx^{i+j}\)

    多项式的点值表示

    众所周知,一个 \(n\) 次函数可以由 \(n+1\) 个点对 \(\left(x_i,y_i\right)\) 唯一确定(其中 \(x_i\) 互不相同)
    因此可以用 \(n+1\) 个点对表示出一个 \(n\) 次的多项式。
    \(A=\left(\left(x_1,y_1\right),\left(x_2,y_2\right),\dots,\left(x_n,y_n\right)\right)\)\(B=\left(\left(x_1,y'_1\right),\left(x_2,y_2'\right),\dots,\left(x_n,y_n'\right)\right)\)
    \(A\times B=\left(\left(x_1,y_1y_1'\right),\left(x_2,y_2y_2'\right),\dots,\left(x_n,y_ny_n'\right)\right)\)
    所以如果多项式的点值已知,那么多项式的乘法是 \(O\left(n\right)\) 的。
    若把点对的数量扩大至 \(2n+1\),则可以唯一确定 \(C\)

    DFT & IDFT

    DFT(离散傅里叶变换) 是系数数列 \(a_0,a_1,\cdots,a_n\) 转为点值数列的过程,IDFT(逆离散傅里叶变换)DFT 的逆运算。
    注意到单位根存在以下性质

    \[\omega_{2n}^{k}+\omega_{2n}^{n+k}=0 \]

    \[\left(\omega_{2n}^k\right)^2=\omega_{n}^k \]

    考虑用单位根的性质优化 DFT。
    对于多项式 \(A\left(x\right)=\sum_{i=0}^{2n-1}a_ix^i\),其中 \(2n=2^k\left(k\in N^{+}\right)\),我们将 \(\omega_{2n}^0,\omega_{2n}^1,\cdots,\omega_{2n}^{{2n}-1}\) 代入公式计算点值。
    现在重定义两个多项式

    \[A_0\left(x\right)=a_0+a_2x+a_4x^2+\cdots+a_{2n}x^n \]

    \[A_1\left(x\right)=a_1+a_3x+a_5x^2+\cdots+a_{2n-1}x^n \]

    显然

    \[A\left(x\right)=A_0\left(x^2\right)+xA_1\left(x^2\right) \]

    将单位复根代入上式

    \[\begin{aligned}A\left(\omega_{2n}^k\right) & =A_0\left(\left(\omega_{2n}^k\right)^2\right)+\omega_{2n}^kA_1\left(\left(\omega_{2n}^k\right)^2\right)\\&=A_0\left(\omega_n^k\right)+\omega_{2n}^kA_1\left(\omega_n^k\right)\\A\left(\omega_{2n}^{n+k}\right)&=A_0\left(\left(\omega_{2n}^{n+k}\right)^2\right)+\omega_{2n}^{n+k}A_1\left(\left(\omega_{2n}^{n+k}\right)^2\right)\\&=A_0\left(\omega_n^k\right)-\omega_{2n}^kA_1\left(\omega_n^k\right)\end{aligned} \]

    发现对于 \(k\in\left[0,1,\cdots,n-1\right]\) \(A\left(\omega_{2n}^k\right)\)\(A\left(\omega_{2n}^{n+k}\right)\) 是可以在一起计算的。
    于是有以下伪代码

    function FFT(complex A[], int Length)
        if Length == 1 return
        complex A0[Length / 2], A1[Length / 2]
        for int i = 0 to Length / 2 - 1 with i += 1
        	A0[i] = A[i * 2]
            A1[i] = A[i * 2 + 1]
    	FFT(A0, Length / 2)
    	FFT(A1, Length / 2)
    	complex wn = (cos(2 * Pi / Length), sin(2 * Pi / Length))
    	complex w = (1, 0)
    	for int i = 0 to Length / 2 - 1 with i += 1, w *= wn
    		A[i] = A0[i] + A1[i] * w
    		A[i + Length / 2] = A0[i] - A1[i] * w
    

    IDFT 是 DFT 的逆运算,令 \(DFT\left(\left\{a_i\right\}\right)=\left\{b_i\right\}\),已知

    \[b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki} \]

    存在结论

    \[a_k=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki} \]

    证明:将前式代入后式

    \[\begin{aligned}a_k&=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\\&=\frac 1 n\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_n^{-ki}\omega_{n}^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i\left(j-k\right)}\end{aligned} \]

    考虑 \(\sum_{i=0}^{n-1}\omega_n^{i\left(j-k\right)}\)
    \(j=k\)

    \[\sum_{i=0}^{n-1}\omega_n^{i\left(j-k\right)}=\sum1=n \]

    \(j\neq k\)

    \[\begin{aligned}\sum_{i=0}^{n-1}\omega_n^{i\left(j-k\right)}&=\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i\\&=\frac{1-\left(\omega_n^{j-k}\right)^n}{1-\omega_n^{j-k}}\\&=\frac{1-\left(\omega_n^n\right)^{j-k}}{1-\omega_n^{j-k}}=\frac{1-1}{1-\omega_n^{j-k}}=0\end{aligned} \]

    所以,

    \[\begin{aligned}a_k&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i\left(j-k\right)}\\&=\frac 1 n\cdot na_k=a_k\end{aligned} \]

    发现 DFT 和 IDFT 的公式形式一样,所以可以用一个函数实现。

    蝴蝶操作

    考虑 FFT 的分治过程,以 \(n=16\) 为例

    \[\left\{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15\right\} \]

    \[\left\{0,2,4,6,8,10,12,14\right\},\left\{1,3,5,7,9,11,13,15\right\} \]

    \[\left\{0,4,8,12\right\},\left\{2,6,10,14\right\},\left\{1,5,9,13\right\},\left\{3,7,11,15\right\} \]

    \[\left\{0,8\right\},\left\{4,12\right\},\left\{2,10\right\},\left\{6,14\right\},\left\{1,9\right\},\left\{5,13\right\},\left\{3,11\right\},\left\{7,15\right\} \]

    其二进制下表示,

    \[\left\{0000,1000\right\},\left\{0100,1100\right\},\left\{0010,1010\right\},\left\{0110,1110\right\}, \]

    \[\left\{0001,1001\right\},\left\{0101,1101\right\},\left\{0011,1011\right\},\left\{0111,1111\right\} \]

    观察发现,若干次蝴蝶操作(奇偶分治)后,其数位二进制翻转后是连续的。
    对于二进制翻转,可用递推计算,rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0)
    考虑用蝴蝶定理使 FFT 的过程避免递归,可以先将 \(\left\{a_i\right\}\)\(rev\) 序重新排列,在分治树上从下往上,

    function FFT(complex A[], int Length, int on)
    	for int i = 0 to Length - 1 with i += 1
    		if i < Rev[i]
    			swap(A[i], A[Rev[i]])
    	for int k = 1 to n - 1 with k *= 2
    		complex wn = (cos(Pi / k), sin(on * Pi / k))
    		for int i = 0 to n with i += k * 2
    			complex w = (1, 0)
    			for int j = i to i + k - 1 with j += 1
    				complex x = a[j]
    				complex y = a[j + k]
    				a[j] = x + y
    				a[j + k] = x - y
    	if on == -1
    		for int i = 0 to n - 1 with i += 1
    			a[i] /= n
    

    然后枚举当前合并的长度 \(2k\),枚举合并区间开始位置 \(i\),枚举区间中的元素 \(a_j\)

    struct node{
    	double x,y;
    }A[MAXN],B[MAXN];
    node operator+ (node p,node q){return (node){p.x+q.x,p.y+q.y};}
    node operator- (node p,node q){return (node){p.x-q.x,p.y-q.y};}
    node operator* (node p,node q){return (node){p.x*q.x-p.y*q.y,p.x*q.y+p.y*q.x};}
    void FFT(node *X,int flag){
        for (int i=0;i<len;i++) if (i<rev[i]) swap(X[i],X[rev[i]]);
        for (int l=1;l<len;l<<=1){
            node T=(node){cos(Pi/l),flag*sin(Pi/l)};
            for (int s=0;s<len;s+=(l<<1)){
                node t=(node){1,0};
                for (int k=0;k<l;k++,t=t*T){
                    node Nx=X[s+k],Ny=t*X[s+k+l];
                    X[s+k]=Nx+Ny,X[s+k+l]=Nx-Ny;
                }
            }
        }
        for (register int i=0;i<len;i++) A[i].x/=len;
        return;
    }
    int main(){
        n=qr(),m=qr();
        for (int i=0;i<=n;i++) A[i].x=qr();
        for (int i=0;i<=m;i++) B[i].x=qr();
        while (len<=n+m) len<<=1,++L;
        for (int i=0;i<len;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
        FFT(A,1);FFT(B,1);
        for (int i=0;i<=len;i++) A[i]=A[i]*B[i];
        FFT(A,-1);
        for (int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x+0.5));
        return 0;
    }
    
    
  • 相关阅读:
    培训课程大纲
    十个心理细节
    海马记忆训练
    手把手教你_怎么找android应用的包名和启动activity
    LoaderManager使用具体解释(四)---实例:AppListLoader
    strtok函数
    猫猫学iOS 之微博项目实战(2)微博主框架-自己定义导航控制器NavigationController
    OpenCV实践之路——Python的安装和使用
    状态模式
    一个有意思的Ruby Webdriver超时问题的解决过程
  • 原文地址:https://www.cnblogs.com/bestlxm/p/12156371.html
Copyright © 2020-2023  润新知