2. 程式人生 > >【Learning】多項式乘法與快速傅裏葉變換（FFT）

【Learning】多項式乘法與快速傅裏葉變換（FFT）

阿新 • • 發佈：2017-10-11
alt 技術 cos 相同 es2017 define ostream 意思 呵呵

簡介：

　　FFT主要運用於快速卷積，其中一個例子就是如何將兩個多項式相乘，或者高精度乘高精度的操作。

　　顯然暴搞是$O(n^2)$的復雜度，然而FFT可以將其將為$O(n lg n)$。

　　這看起來十分玄學，因為怎麽看它們的相乘操作都逃不過$O(n^2)$，FFT是如何再減少復雜度的呢？

　　講到FFT就不可避免地出現公式，但實際上它們都是比較容易理解的。

　　

全局思路

　　設兩個次數界均為$n$的多項式$\begin{aligned}A(x)&=a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}\\B(x)&=b_0x^0+b_1x^1+b_2x^2+...+b_{n-1}x^{n-1}\end{aligned}$

　　那麽我們要求$C=A*B$。

　　我們把$A$、$B$和$C$的函數圖像畫出來，考慮每一個取值點：$C(x_1)=A(x_1)*B(x_1)$，$C(x_2)=A(x_2)*B(x_2)$，$C(x_3)=A(x_3)*B(x_3)$......

　　技術分享

　　如果我們能對$A$和$B$求出它們在$x_1,x_2,...x_{2n-1}$處的值$A(x_1),A(x_2),...,A(x_{2n-1})$與$B(x_1),B(x_2),...,B(x_{2n-1})$。

　　那麽將它們一一對應地相乘，就可以得到$C(x_1),C(x_2),...,C(x_{2n-1})$。其中，$C(x)=A(x)*B(x)$。

　　接著，我們可以利用$C$的這$2n-1$個點，對$C$進行插值，求出$C$的解析式。由於兩個次數界為$n$的多項式相乘後是一個次數界為$2n-1$的多項式，因此我們需要$2n-1$個點，才能對$C$進行準確插值。

　　問題是，上面三步的時間復雜度分別為$O(n^2)$、$O(n)$和$O(n^2)$，還是沒有什麽改進。它們的名字分別是：DFT，點值乘法，逆DFT。

　　改進以後，它們分別是FFT，點值乘法，逆FFT。時間復雜度分別為$O(nlgn)$，$O(n)$，$O(nlgn)$。

　　

那就來改進吧（DFT $O(n^2)$ ---> FFT $O(nlgn)$）

　　本質老是沒有飛躍，多半是廢了

，有一個因素是取值點$x_1,x_2,...,x_{2n-1}$停留在實數範圍內，沒有太多特殊性質。

　　但如果用復數呢？

　　

　　定義$n$次單位復數根為滿足$\omega^n=1$的復數$\omega$。$n$次單位復數根恰好有$n$個：$\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$，它們的$n$次方都為1。

　　其中，$\omega_n^0=e^{2\pi i/n}$， $\omega_n^x=(\omega_n^0)^x$， $\omega_n^{n/2}=(e^{2\pi i/n})^{n/2}=e^{\pi i}=-1$。

　　（以下用$n$來代替之前提到的$2n-1$；用$A$表示一個次數界為$n$的多項式，上文提到的"$A$"和"$B$"的操作都是同理的）

　　我們把這$n$個$n$次單位復數根作為取值點，求出$A(\omega_n^0),A(\omega_n^1),...,A(\omega_n^{n-1})$。

　　這一步叫做離散傅裏葉變換(DFT)。對於$k=0,1,...,n-1$，它要求$y_k=A(\omega_n^k)=\sum\limits_{j=0}^{n-1}a_j\omega_n^{kj}$

　　然而若不利用單位復數根的性質，復雜度仍然是$O(n^2)$的。

　　

單位復數根的性質

　　這$n$個復數有神秘性質，主要用到三個：

　　1. $\omega_n^{k+n/2}=\omega_n^k*\omega_n^{n/2}=\omega_n^k*-1=-\omega_n^k$，

　　　　什麽意思呢：比如$n=8$時，$n$個單位負數根，$\omega_n^0$和$\omega_n^4$互為相反數，$\omega_n^1$和$\omega_n^5$互為相反數.....也就是$[0,n/2-1)$與$[n/2-1,n)$對應的單位復數根互為相反數。

　　2.消去引理：$\omega_{an}^{ak}=(e^{2\pi i/an})^{ak}=(e^{2\pi i/n})^k=\omega_n^k$，類似於分數約分。

　　3.折半引理：$$\begin{aligned}(\omega_n^k)^2&=\omega_n^{2k}=\omega_{n/2}^k\\(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}*\omega_n^n&=\omega_n^{2k}=\omega_{n/2}^k\end{aligned}$$

　　　　什麽意思呢？如果把$n$個單位根分成兩組$\omega_n^0...\omega_n^{n/2-1}$ 和 $\omega_n^{n/2}...\omega_n^{n-1}$，兩兩對應位置的單位根的平方是相同的。

　　　　如$n==8$時：

　　　　$(\omega_8^0)^2=(\omega_8^4)^2=\omega_{4}^0\\(\omega_8^1)^2=(\omega_8^5)^2=\omega_{4}^1\\(\omega_8^2)^2=(\omega_8^6)^2=\omega_{4}^2\\(\omega_8^3)^2=(\omega_8^7)^2=\omega_{4}^3$

　　　　也就是$n$個$n$次單位復數根的平方的集合，等於$n/2$個$n/2$次單位復數根的集合。

多項式的拆分

　　我們回來看一下$A$可以如何拆分：記$A$的系數為$a_0,a_1,...,a_{n-1}$。

　　如果我們設$\begin{aligned}A_0(x)&=a_0x^0+a_2x^1+a_4x^2+...+a_{n-2}x^{n/2}\\A_1(x)&=a_1x^0+a_3x^1+a_5x^2+...+a_{n-1}x^{n/2}\end{aligned}$，也就是將$A$的系數奇偶分組，成為兩個次數界為$n/2$的多項式。

　　那麽有$$A(x)=A_0(x^2)+x*A_1(x^2)$$

　　我們求的是$A(\omega_n^0),A(\omega_n^1),...,A(\omega_n^{n-1})$，那麽轉換一下就變成求

$$\begin{aligned}
A(\omega_n^0)&=A_0((\omega_n^0)^2)+\omega_n^0*A_1((\omega_n^0)^2)\\
A(\omega_n^1)&=A_0((\omega_n^1)^2)+\omega_n^1*A_1((\omega_n^1)^2)\\
&...\\
A(\omega_n^{n-1})&=A_0((\omega_n^{n-1})^2)+\omega_n^{n-1}*A_1((\omega_n^{n-1})^2)\\
\end{aligned}$$

　　求解$A_0$和$A_1$在$n$個單位復數根，我們用遞歸實現。

　　我們發現代入$A_0$和$A_1$的參數是一個單位復數根的平方，這意味著代入$A_0$和$A_1$的單位復數根並沒有$n$個。根據折半引理，代入$A_0$和$A_1$的總共只有$n/2$個不同的數：$\omega_{n/2}^0,\omega_{n/2}^1,...,\omega_{n/2}^{n/2-1}$，因為$(\omega_n^k)^2=(\omega_n^{k+n/2})^2$。

　　我們像上面把單位復數根分為$[0,n/2)$和$[n/2,n)$兩組，觀察$A(\omega_n^k)$和$A(\omega_n^{k+n/2})$，也就是相對的兩個單位復數根的代入：

　　　　

\begin{aligned}
A(\omega_n^k)&=A_0(\omega_{n/2}^k)+\omega_n^k*A_1(\omega_{n/2}^k)\\
A(\omega_n^{k+n/2})&=A_0(\omega_{n/2}^k)+\omega_n^{k+n/2}*A_1(\omega_{n/2}^k)\\
&=A_0(\omega_{n/2}^k)-\omega_n^k*A_1(\omega_{n/2}^k)
\end{aligned}

　　它們長得好像！

　　這下可好，我們只需要遞歸求解$A_0(\omega_n^0...\omega_n^{n/2-1})$和$A_1(\omega_n^0...\omega_n^{n/2-1})$，就可以求出$A(\omega_n^0...\omega_n^{n-1})$了。

　　時間復雜度下降的原因就在於，用$n/2$次的遞歸得到的數據，可以求出右半邊的數值。

點值乘法 （呵呵 $O(n)$ 不優化了吧這個）

　　對於$A$和$B$都進行DFT後，我們對$n$個點值直接相乘，得到$C$的$n$個點值。

IFFT （IDFT $O(n^2)$ ---> IFFT $O(nlgn)$）

　　如果我們知道$C$的$n$個點值，如何知道$C$的解析式呢？

　　我們看一下DFT的矩陣形式：$y=V_na$，分別與下式對應：

$$\begin{bmatrix}
y_0\\y_1\\y_2\\y_3\\.\\.\\y_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
1&1&1&1&...&1\\
1&\omega_n&\omega_n^2&\omega_n^3&...&\omega_n^{n-1}\\
1&\omega_n^2&\omega_n^4&\omega_n^6&...&\omega_n^{2(n-1)}\\
1&\omega_n^3&\omega_n^6&\omega_n^9&...&\omega_n^{3(n-1)}\\
...&...&...&...&...&...\\
1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&...&\omega_n^{(n-1)(n-1)}
\end{bmatrix}
*
\begin{bmatrix}
a_0\\a_1\\a_2\\a_3\\.\\.\\a_{n-1}
\end{bmatrix}$$

　　我們所求的是$a$，而$a=yV_n^{-1}$，求出$V_n$的逆矩陣即萬事大吉了。

　　定理：對於$V_n$，$(k,j)$處的元素為$\omega_n^{kj}$。

　　　　 而對於$V_n^{-1}$，$(k,j)$處的元素為$\omega_n^{-kj}/n$。

　　　　 如果想簡單證明的話，將$V_n^{-1}$寫出來算一算就好。（可以參見算導）

　　那麽$a_j=\frac{1}{n}\sum\limits_{k=0}^{n-1}y_k\omega_n^{-kj}$。

　　看回上面DFT的算式表達，我們發現它們長得幾乎一樣：IFFT的表達，僅僅是多了一個$\frac{1}{n}$，以及單位復數根的指數取負數。

　　這就非常棒了：IFFT的程序其實和FFT一樣，只不過單位復數根替換一下，算完以後，每一個數值都除去$n$即可，具體參見代碼解釋。

END

　　FFT的應用，主要是將問題轉化成如DFT式子的形式，用FFT來進行加速或計算的操作。

　　附上遞歸版代碼和非遞歸版代碼：

技術分享 
#include <cstdio>
#include <vector>
#include <cmath>
#define max(a,b) ((a)>(b)?(a):(b))
using namespace std;
const int N=50010;
const double Pi=3.14159265358979323846;
struct Comp{//手寫了一個復數類 
    double a,b; 
    Comp(){a=b=0.0;}
    Comp(double x,double y){a=x;b=y;}
    friend Comp operator + (Comp x,Comp y){
        return Comp(x.a+y.a,x.b+y.b);
    }
    friend Comp operator - (Comp x,Comp y){
        return Comp(x.a-y.a,x.b-y.b);
    }
    friend Comp operator * (Comp x,Comp y){
        return Comp(x.a*y.a-x.b*y.b,x.b*y.a+x.a*y.b);
    }
};  
typedef vector<Comp> vc; 
int A,B,type,len;
vc a,b,c;
vc fft(vc u,int flag){//flag標識是否為逆FFT 
    int n=u.size();
    if(n==1) return u;//規模為1時，只有一個常數項的多項式的FFT就為這個常數，可以直接返回了 
    Comp w_n=Comp(cos(2*Pi/n),sin(2*Pi/n)),w=Comp(1,0);//算出單位復數根的底w_n；w是用來叠代的，減少計算次數 
    if(flag) w_n.b*=-1.0;//逆FFT與FFT的不同 
    vc a0,a1,v;
    a0.clear(); a1.clear(); v.clear();
    for(int i=0;i<n;i++){//系數按奇偶分組 
        if(i&1) a1.push_back(u[i]);
        else    a0.push_back(u[i]);
        v.push_back(Comp(0,0));
    }
    //遞歸求解A0和A1 
    a0=fft(a0,flag);
    a1=fft(a1,flag);
    //用一半的數據，綜合算出全部的結果，w在此處不斷乘上w_n，就保證它是w_n的k次方 
    for(int k=0;k<=n/2-1;k++){
        v[k]=a0[k]+w*a1[k];
        v[k+n/2]=a0[k]-w*a1[k];
        w=w*w_n;
    }
    return v;
}
int main(){
//原題：求兩個多項式相乘後的系數（系數都為整數）
    scanf("%d%d%d",&A,&B,&type);
    A++; B++;
    for(int i=0,x;i<A;i++) scanf("%d",&x),a.push_back(Comp(x,0));
    for(int i=0,x;i<B;i++) scanf("%d",&x),b.push_back(Comp(x,0));
    len=1;//算出高位補齊len（上文提到的至少需要2n-1個點），並把兩個多項式的次數都擴展到len 
    //代碼裏的len指的是上文提到的n 
    while(len<(max(A,B)*2)) len<<=1;
    for(int i=A;i<len;i++) a.push_back(Comp(0,0));
    for(int i=B;i<len;i++) b.push_back(Comp(0,0));
    //求兩個多項式在n個單位復數根的值O(nlgn)
    a=fft(a,0);
    b=fft(b,0);
    //點值乘法 O(n)
    for(int i=0;i<len;i++) c.push_back(a[i]*b[i]);
    //對點值乘法的結果進行逆FFT O(nlgn)
    c=fft(c,1);
    for(int i=0;i<A+B-1;i++) printf("%d ",(int)(c[i].a/len+0.5));//除去len，四舍五入（這題是整數） 
    return 0;
}
遞歸版FFT 技術分享 
#include <cstdio>
#include <iostream>
#include <cmath>
#define max(a,b) ((a)>(b)?(a):(b))
using namespace std;
const int N=50010;
const double Pi=3.14159265358979323846;
struct Comp{
    double a,b;    
    Comp(){a=b=0.0;}
    Comp(double x,double y){a=x;b=y;}
    friend Comp operator + (Comp x,Comp y){return Comp(x.a+y.a,x.b+y.b);}
    friend Comp operator - (Comp x,Comp y){return Comp(x.a-y.a,x.b-y.b);}
    friend Comp operator * (Comp x,Comp y){return Comp(x.a*y.a-x.b*y.b,x.b*y.a+x.a*y.b);}
}a[N*4],b[N*4];    
int A,B,type,n;
inline int rev(int x){
    int ret=0;
    for(int i=1;i<n;i<<=1,x>>=1)
        ret=(ret<<1|(x&1));
    return ret;
}
void fft(Comp *a,int f){
    int lg=log2(n),len;    
    Comp w,w_n,u,v;
    for(int i=0,t;i<n;i++)
        if(i<(t=rev(i))) swap(a[i],a[t]);
    for(int i=1;i<=lg;i++){
        len=1<<i;
        w_n=Comp(cos(2*Pi/len),sin(2*Pi/len)*f);
        for(int j=0;j<n;j+=len){
            w=Comp(1,0);
            for(int k=0;k<=len/2-1;k++){
                u=a[j+k]; v=w*a[j+k+len/2];
                a[j+k]=u+v; a[j+k+len/2]=u-v;
                w=w*w_n;
            }
        }
    }
}
int main(){
    scanf("%d%d%d",&A,&B,&type);
    for(int i=0,x;i<A;i++) scanf("%lf",&a[i].a);
    for(int i=0,x;i<B;i++) scanf("%lf",&b[i].a);
    for(n=1;n<A+B;n<<=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<A+B-1;i++) 
        printf("%d\n",(int)(a[i].a/n+0.5));
    return 0;
}
非遞歸版FFT（常數小）

【Learning】多項式乘法與快速傅裏葉變換（FFT）

相關推薦

【Learning】多項式乘法與快速傅裏葉變換（FFT）

alt 技術 cos 相同 es2017 define ostream 意思 呵呵 簡介： 　　FFT主要運用於快速卷積，其中一個例子就是如何將兩個多項式相乘，或者高精度乘高精度的操作。 　　顯然暴搞是$O(n^2)$的復雜度，然而FFT可以將其將為$O(n lg

【模板】快速傅裏葉變換（FFT）（BZOJ2179）

turn const 一個 algorithm har string using tchar out Description 給出兩個n位10進制整數x和y，你需要計算x*y。 Input 第一行一個正整數n。 第二行描述一個位數為n的正整數x。 第三行描述一個位數為

OI中的快速傅裏葉變換（FFT）

ali 系列 我認 交流 cpp 實現 AS 需要 部分 快速傅裏葉變換（FFT） 前言 關於這篇文章 ? ? 非常高興能有機會來探討快速傅裏葉變換，也就是大家熟知的 \(FFT\) 在 \(OI\) 中的運用。以前了解過一次 \(FFT\) ，現在過了幾個月，數學和 \

【知識總結】快速傅立葉變換（FFT）

這可能是我第五次學FFT了……菜哭qwq 先給出一些個人認為非常優秀的參考資料： 一小時學會快速傅立葉變換（Fast Fourier Transform） - 知乎 小學生都能看懂的FFT！！！ - 胡小兔 - 部落格園 快速傅立葉變換（FFT）用於計算兩個\(n\)次多項式相乘，能把複雜度從樸素的\

【演算法】快速傅立葉變換（FFT）的遞迴實現

FFT是數字訊號處理中的重要演算法，在matlab中可以直接呼叫fft()函式，本文是C++版的FFT演算法，用遞迴方式進行實現。 //================================== //Signal_Process_FFT_v1.cpp //====

【模板】快速傅立葉變換（FFT）（BZOJ2179）

Description 給出兩個n位10進位制整數x和y，你需要計算x*y。 Input 第一行一個正整數n。 第二行描述一個位數為n的正整數x。 第三行描述一個位數為n的正整數y。 Output 輸出一行，即x*y的結果。 Code #include<cstdio> #

離散傅立葉變換（DFT）和快速傅立葉變換（FFT）原理與實現

目錄 1、影象變換 2、離散傅立葉變換（Discrete Fourier Transform） 3、DFT性質 4、DFT與數字影象處理 5、FFT-快速傅立葉變換 6、DFT與FFT的演算法實現 1. 影象變換 — —數學領域中有很多種變換，如傅立葉變換、拉普拉斯變

快速傅立葉變換（FFT）及其應用

有一個 swap max read mes turn scan 原本 color 在信息學競賽中FFT只有一個用處那就是加速多項式的乘法 多項式乘法原本的時間復雜度是O（n^2）的，然後經過FFT之後可以優化為O（nlogn） FFT就是將系數表示法轉化成點值表示法相乘，再

快速傅立葉變換（FFT）（學習筆記）

學習了一波 F F T FFT

11.頻域裡的卷積——介紹，傅立葉變換和卷積，快速傅立葉變換（FFT）_1

目錄 介紹 FFT 介紹 我們將繼續討論頻率分析以及如何用頻率分量的概念來研究影象。如果你還記得上次我們講過的基於頻率的影象分解的概念。我們通過給你們看這張照片來回憶它（如圖）。這是著名的Dali圖片，當你在那裡允許高頻影象時，你會看到一個女人在欣賞地中海之類的東

我所理解的快速傅立葉變換（FFT）

轉載自：https://blog.csdn.net/shenziheng1/article/details/52891807 1.歷史放在最前頭 首先FFT是離散傅立葉變換(DFT)的快速演算法，那麼說到FFT，我們自然要先講清楚傅立葉變換。先來看看傅

快速傅立葉變換（FFT）學習

首先，在寫這篇部落格之前，我還沒有完全學會FFT。 先把會的部分打好，加深一下記憶（也可以說是做筆記吧）。 初三了，還不會FFT，要退役嘍…… 多項式乘法 點開這篇部落格之前，你就應該知道，FFT是用來求多項式乘法的。 什麼是多項式，什麼是多項式乘法？ 不講。初一內容。 如

卷積和快速傅立葉變換（FFT）的實現

卷積運算 卷積可以說是影象處理中最基本的操作。線性濾波通過不同的卷積核，可以產生很多不同的效果。假如有一個要處理的二維影象，通過二維的濾波矩陣（卷積核），對於影象的每一個畫素點，計算它的領域畫素和濾波器矩陣的對應元素的乘積，然後累加，作為該畫素位置的值。關於影

快速傅立葉變換（FFT）（筆記）

今天看了一下午FFT，總算看懂了。怕以後忘記，留個圖片。 首先以下圖片來自 https://www.cnblogs.com/zwfymqz/p/8244902.html?mType=Group https://www.cnblogs.com/wangyh1008/p/9325715.ht

影象處理 之 一維快速傅立葉變換（FFT）

# -*- coding: utf-8 -*- """ Created on Sun Jul 8 21:05:51 2018 @author: Diko """ import numpy def

快速傅立葉變換（FFT）和數論變換（NTT）模板

#include<bits/stdc++.h> #define ll long long #define dob complex<double> const int N=120010; const double pi=acos(-1.0); usin

快速傅立葉變換（FFT）C語言函式

/********************************************************************* 快速福利葉變換C函式 函式簡介：此函式是通用的快速傅立葉變換C語言函式，移植性強，以下部分不依 賴硬體。此函式採用聯合體

傅裏葉變換（FFT）的多相濾波結構實現

class hold es2017 要求 sca 接收機 .html mar try 作者：桂。 時間：2017-09-25 14:53:01 鏈接：http://www.cnblogs.com/xingshansi/p/7591868.html 前言

【bzoj2179】FFT快速傅裏葉變換（優化高精度乘法）

efi swa 快速傅裏葉變換 得到 urn lse can sin %s #include<bits/stdc++.h> using namespace std; #define pi acos(-1) typedef complex<double&g

洛谷P1919 A*B problem 快速傅裏葉變換模板 [FFT]

sca 快速 https -o ora stream ros for AR 　　題目傳送門 A*B problem 題目描述 給出兩個n位10進制整數x和y，你需要計算x*y。 輸入輸出格式 輸入格式： 第一行一個正整數n。 第二行描述一個位數為n的正整數