快速傅里葉變換FFT& 數論變換NTT

時間  2019-11-06
標籤 快速 傅里葉 變換 fft 數論 ntt 简体版
原文   原文鏈接

相關知識

時間域上的函數f(t)通過傅里葉變換（Fourier Transform）變成頻率域上的F(w)，也就是用一些不一樣頻率正弦曲線的加 權疊加獲得時間域上的信號。html

\[ F(\omega)=\mathcal{F}[f(t)]=\int\limits_{-\infty}^\infty f(t)e^{-iwt}dt \]c++

傅里葉逆變換是將頻率域上的F(w)變成時間域上的函數f(t)，通常稱\(f(t)\)爲原函數，稱\(F(w)\)爲象函數。原函數和象函數構成一個傅里葉變換對。web

\[ f(t)=\mathcal{F^{-1}}[F(w)]=\frac 1 {2\pi}\int\limits_{-\infty}^\infty F(w)e^{iwt}dw \]算法

離散傅里葉變換（DFT）是傅里葉變換在時域和頻域都呈現離散的形式。數組

\[ x_n=\sum_{k=0}^{N-1}X_ke^{\frac{2\pi i}{N} kn},n=0,...,N-1 \]函數

固然也就有離散傅里葉逆變換（IDFT）：ui

\[ x_n=\frac 1 N\sum_{k=0}^{N-1}X_ke^{-\frac{2\pi i}{N} kn},n=0,...,N-1 \]spa

快速傅里葉變換FFT

快速傅里葉變換（FFT）利用分治的思想，將 DFT 分解成更小的 DFT進行計算，能夠在\(O(n\log n)\)時間複雜度內完成 DFT 的算法。可用來加速卷積和多項式乘法。code

多項式乘法

對於\(C(x)=A(x)B(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{i}a_jb_{i-j}x^i\)，（A的最高項次數加上B的最高項次數爲n-1）直接計算是\(O(n^2)\)的。咱們將它們變成點值表達式orm

\[ A(x)：{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}\\\\ B(x)：{(x_0,y_0'),(x_1,y_1'),...,(x_{n-1},y_{n-1}')}\\\\ C(x)：{(x_0,y_0y_0'),(x_1,y_1y_1'),...,(x_{n-1},y_{n-1}y_{n-1}')} \]

那麼就能夠\(O(n)\)計算出C(x)。

將其轉爲點值表達式須要帶入n個不一樣的x，利用霍納法則，帶入一個x須要\(O(n)\)時間：

\[ A(x_0)=a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0(a_{n-1}))...)) \]

從點值表達式轉換爲係數表達式，用插值法，將n個不一樣的值帶入A，惟一肯定了一個係數表達式（證實可見《算法導論》）。

利用 FFT 能夠在\(O(n\log n)\) 內完成係數到點值，以及從點值到係數表達式。

n次單位複數根

n次單位複數根\(w\)是知足\(w^n=1\)的複數。 主n次單位根爲\(w_n=e^{\frac {2\pi } ni}\)，其它全部n次單位複數根都是\(w_n\)的冪次。

咱們知道複數的指數形式的定義是：

\[ e^{ui}=\cos(u)+i\sin(u) \]

因此n個複數根均勻分佈在以複平面原點爲圓心的單位圓上。

單位根的性質

  1. 消去引理 \(w^{dk}_{dn}=w^k_n\)，帶入定義便可證實。推論:\(w^{n/2}_n=w_2=-1\)

  2. 折半引理\((w^{k+n/2}_n)^2=(w_n^k)^2\)，展開便可證實。推論:\(w^{k+n/2}_n=-w^{k}_n\)

爲了方便分治，咱們將n補到2的次冪（多項式後面加0）。將n個n次單位根帶入A(x)的過程就是對係數向量進行 DFT 。

\[ A(w_n^k)=\sum_{j=0}^{n-1}a_jw^{kj}_n,k=0,1,..,n-1 \]

考慮按指數奇偶分類:

\[ A_{even}(x)=a_0+a_2x+a_4x^2+...\\ A_{odd}(x)=a_1+a_3x+a_5x^2+...\\ \]

那麼

\[ A(x)=A_{even}(x^2)+xA_{odd}(x^2)\\ A(w_n^k)=A_{even}(w_{n}^{2k})+w_n^kA_{odd}(w_n^{2k})\\ \]

當\(k< n/2\)時，由消去引理

\[ A_{even}(w_n^{2k})+w_n^kA_{odd}(w_n^{2k})=A_{even}(w_{n/2}^k)+w_n^kA_{odd}(w_{n/2}^k) \]

也就是至關於對兩個長度爲n/2的多項式\(A_{even}(x)\),\(A_{odd}(x)\)進行 DFT。

指數不小於\(n/2\)的部分與\(k+n/2\)一一對應，由折半引理

\[ A_{even}((w_n^{k+n/2})^2)+w_n^{k+n/2}A_{odd}((w_n^{k+n/2})^2)=A_{even}(w_{n/2}^k)-w_n^kA_{odd}(w_{n/2}^k) \]

因此這個過程是能夠遞歸的，且\(T(n)=2T(n/2)+O(n)\)。

遞歸的FFT（c++代碼）

typedef complex<double> CD;
const double pi = acos(-1);
CD tmp[N],epsilon[N];
void init_epsilon(int n){
    for(int i = 0; i < n; ++i){
        epsilon[i] = CD(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n)); 
        arti_epsilon[i] = conj(epsilon[i]);
    }
}
void recursive_fft(int n, CD* A,int offset, int step, CD* w){
    if(n==1)return;
    int m=n>>1;
    recursive_fft(m,A,offset,step<<1,w);
    recursive_fft(m,A,offset+step,step<<1,w);
    for(int k=0;k<m;++k){
        int pos=2*step*k;
        tmp[k]  =A[pos+offset]+w[k*step]*A[pos+offset+step];
        tmp[k+m]=A[pos+offset]-w[k*step]*A[pos+offset+step];
    }
    for(int i=0;i<n;++i)
        A[i*step+offset]=tmp[i];
}

迭代實現

可是遞歸須要比較大的空間，如何實現迭代的寫法？

觀察遞歸的過程，第一步：

0（000）2（010）4（100）6（110），1（001）3（011）5（101）7（111）

第二步：

0 （000）4（100），2（010） 6（110），1（001）5（101）3（011）7（111）

將下標的二進制翻轉過來就是：

000，001，010，011，100，101，110，111

對應了0，1，2，3，...

用reverse(i)表示i的二進制翻轉後的數。就至關於二進制從高位到低位的+1，是從左往右遇到第一個0，就改成1，左邊的1改成0。

int reverse(int x){
    for(int l=1<<bit_length;(x^=l)<l;l>>=l);
    return x;
}

咱們獲得遞歸最後一步的數組：

0，4，2，6，1，5，3，7

就能夠從下到上迭代了。

void bit_reverse(CD* A,int n){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(A[i],A[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
}
void fft(CD* A, int n, CD* w){
    bit_reverse(A,n);
    for(int i=2;i<=n;i<<=1)//自下到上，i爲每一層的步長，或者說子問題長度
      for(int j=0,m=i>>1;j<n;j+=i)//j爲偏移量，或者說這一層的每個子問題的起點
        for(int k=0;k<m;++k){//k=0..i/2，計算出子問題的第k個和第k+i/2個的值
          CD b=w[n/i*k]*A[j+m+k];
          A[j+m+k]=A[j+k]-b;
          A[j+k]+=b;
        }
}

離散傅里葉逆變換（IDFT）

將點值轉回係數表達式，至關於解方程組\(\vec{y}=V_n\vec{a}\)，其中\(\vec {a}\)是係數向量。\(V_n\)是以下範德蒙德矩陣：

\[ \begin{bmatrix} 1&1&1&1&\cdots &1\\ 1&w;_n&w;_n^2&w;_n^3&\cdots &w;_n^{n-1}\\ 1&w;_n^2&w;_n^4&w;_n^6&\cdots &w;_n^{2(n-1)}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ 1&w;_n^{n-1}&w;_n^{2(n-1)}&w;_n^{3(n-1)}&\cdots &w;_n^{(n-1)(n-1)}\\ \end{bmatrix} \]

那麼\(V_n^{-1}\vec{y}=\vec{a}\)。

可令\([V_n^{-1}]_{kj}=w^{-kj}_n/n\)，只需證實\(V_n^{-1}V_n=I_n\)：

\[ [V_n^{-1}V_n]_{jj'}=\sum_{k=0}^{n-1}(w_n^{-kj}/n)(w_n^{kj'})=\sum_{k=0}^{n-1}w_n^{k(j'-j)}/n = \left\{ \begin{aligned} 0,j'\neq j \\ 1,j'=j \end{aligned} \right. \]

因此只要用\(w_n^{-1}\)替換\(w_n\)，並將結果每一個元素除以n，就能夠計算出\(DDF_n^{-1}\)。

FFT解決高精度乘法，c++代碼

51 Nod 1028 大數乘法 V2

注意FFT由於用了cos和sin，以及是浮點數計算，因此會有精度偏差，故加了0.5。

#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l;i<r;++i)
#define per(i,l,r) for(int i=r-1;i>=l;--i)
#define SZ(x) ((int)(x).size())

typedef double dd;
typedef complex<dd> CD;
const dd PI=acos(-1.0);
const int L=18,N=1<<L;

CD eps[N],inv_eps[N],f[N],g[N];
void init_eps(int p){
    rep(i,0,p)eps[i]=CD(cos(PI*i*2/p),sin(PI*i*2/p)),inv_eps[i]=conj(eps[i]);
}
void fft(CD p[], int n, CD w[]){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(p[i],p[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=n;i<<=1)
        for(int j=0,m=i>>1;j<n;j+=i)
            rep(k,0,m){
                CD b=w[n/i*k]*p[j+m+k];
                p[j+m+k]=p[j+k]-b;
                p[j+k]+=b;
            }
}
int ans[N];
int main(){
    string a,b;
    cin>>a>>b;
    int n=max(SZ(a),SZ(b)),p=1;
    while(p<n)p<<=1;p<<=1;
    rep(i,0,p)f[i]=g[i]=0;
    n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
    n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';

    init_eps(p);
    fft(f,p,eps);fft(g,p,eps);
    rep(i,0,p)f[i]*=g[i];
    fft(f,p,inv_eps);

    int t=0;
    rep(i,0,p){
        ans[i]=t+(f[i].real()+0.5)/p;
        if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
        else t=0;
    }
    bool flag=0;
    per(i,0,p)if(ans[i]||flag){
        printf("%d",ans[i]);flag=1;
    }
    if(flag==0)puts("0");
    return 0;
}

數論變換（NTT）

爲了解決FFT產生的精度偏差，就有了一種在模意義下的方法，運算過程只有整數。它就是NTT（Number-Theoretic Transform）。

原根

設m是正整數，a是整數，則使得同餘式\(a^r \equiv 1 \mod m\) 成立的最小正整數 \(r\) 叫作 \(a\) 對於模 \(m\) 的指數。

若a模m的指數等於φ(m)（歐拉函數），則稱a爲模m的一個原根。

原根的性質

  1. 對於質數p來講原根 \(g\) 知足 \(g^0,g^1,g^2,…,g^{p-1}\) 構成模 \(p\) 的簡化剩餘系。

咱們令\(p=c\cdot 2^k+1\) ，那麼對於2的冪n有 \(n|(p-1)\)。設

\[ g_n=g^{\frac {p-1}n} \]

考慮FFT中，須要用到單位根的如下性質：

  1. \(w_n^k(0\le k < n)\)互不相同，保證點值表示的合法；

  2. \(\omega_{n} ^ {2k} = \omega_{n/2} ^k\)，用於分治；

  3. \(\omega_n ^ { k + \frac{n}{2} } = -\omega_n ^ k\)，用於分治；

  4. 僅當\(k \neq 0\)時，\(\sum_{j=0}^{n-1}(w_n^k)^j=0\)，用於逆變換。

原根是否也有以上性質呢？

  1. 令\(p=c\cdot 2^k+1\)，取\(n=2^m\)，則(g_n^k(0\le k<n)=g^{\frac {(p-1)k}="" n}="g^0,g^{c\cdot" 2^{k-m}},g^{2c\cdot="" 2^{k-m}}…,g^{(p-1){c\cdot="" 2^{k-m}}})<="" span="">，互不相等。</n)=g^{\frac>


    這裏有個小知識點，若 \(a_1,a_2,..a_p\) 是模p的徹底剩餘系，那麼 \(aa_1+b,aa_2+b,..aa_p+b\)也是模p的徹底剩餘系，證實：若 \(aa_j+ba\equiv a_i+b \mod p\)則 \(a_i \equiv a_j \mod p\)，矛盾。

  2. \(g_n^{2k}=g^{\frac {2k(p-1)} n}=g^{\frac{k(p-1)}{n/2}}=g_{n/2}^k\)

  3. \(g_n^n=g^{p-1}\equiv 1\mod p\Rightarrow g^{\frac {p-1}2}\equiv\sqrt 1 \mod p\)，由於\(g^k\)互不相同，因此\(g^{\frac {p-1}2}\equiv -1\mod p\)，因而\(g^{k+\frac n 2}_n=g_n^k\cdot g_n^{\frac n 2}=g_n^k\cdot g^{\frac {p-1} 2}=-g_n^k\)。

  4. 當\(k \neq 0\)時，\(\sum_{j=0}^{n-1}(g_n^k)^j=\frac{1-(g_n^k)^n}{1-g_n^k}\)，由3得，\(g_n^n=1\)，因此\(1-(g_n^k)^n=1-(g_n^n)^k=0\)。不然等於n。

所以在FFT中用原根代替單位根就是NTT了。若是要求模數任意，只要再用 中國剩餘定理（CRT）合併便可。

NTT的c++代碼

仍然是上面高精度乘法那題

#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l;i<r;++i)
#define per(i,l,r) for(int i=r-1;i>=l;--i)
#define SZ(x) ((int)(x).size())

typedef long long LL;
const int L=18,N=1<<L;

const LL C = 479;
const LL P = (C << 21) + 1;
const LL G = 3;

LL qpow(LL a, LL b, LL m){
    LL ans = 1;
    for(a%=m;b;b>>=1,a=a*a%m)if(b&amp;1)ans=ans*a%m;
    return ans;
}

LL eps[N],inv_eps[N],f[N],g[N];
void init_eps(int n){
    LL t=(P-1)/n, invG=qpow(G,P-2,P);
    rep(i,0,n) eps[i]=qpow(G,t*i,P),inv_eps[i]=qpow(invG,t*i,P);
}
void fft(LL p[], int n, LL w[]){
    for(int i=0,j=0;i<n;++i){
        if(i>j)swap(p[i],p[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=n;i<<=1)
        for(int j=0,m=i>>1;j<n;j+=i)
            rep(k,0,m){
                LL b=w[n/i*k]*p[j+m+k]%P;
                p[j+m+k]=(p[j+k]-b+P)%P;
                p[j+k]=(p[j+k]+b)%P;
            }
}
LL ans[N];
int main(){
    string a,b;
    cin>>a>>b;
    int n=max(SZ(a),SZ(b)),p=1;
    while(p<n)p<<=1;p<<=1;
    rep(i,0,p)f[i]=g[i]=0;
    n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
    n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';

    init_eps(p);
    fft(f,p,eps);fft(g,p,eps);
    rep(i,0,p)f[i]=f[i]*g[i]%P;
    fft(f,p,inv_eps);

    int t=0;
    LL invp=qpow(p,P-2,P);
    rep(i,0,p){
        ans[i]=(t+f[i]*invp%P)%P;
        if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
        else t=0;
    }
    bool flag=0;
    per(i,0,p)if(ans[i]||flag){
        printf("%lld",ans[i]);flag=1;
    }
    if(flag==0)puts("0");
    return 0;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程