快速傅里叶变换FFT& 数论变换NTT

介绍快速傅里叶变换算法及数论变换算法原理及 C++实现。

相关知识

时间域上的函数f(t)经过傅里叶变换(Fourier Transform)变成频率域上的F(w),也就是用一些不同频率正弦曲线的加 权叠加得到时间域上的信号。

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

傅里叶逆变换是将频率域上的F(w)变成时间域上的函数f(t),一般称\(f(t)\)为原函数,称\(F(w)\)为象函数。原函数和象函数构成一个傅里叶变换对。

\[ 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):

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

快速傅里叶变换FFT

快速傅里叶变换(FFT)利用分治的思想,将 DFT 分解成更小的 DFT进行计算,可以在\(O(n\log n)\)时间复杂度内完成 DFT 的算法。可用来加速卷积和多项式乘法。

多项式乘法

对于\(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)\)的。我们将它们变成点值表达式

\[ 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<="" 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;
}

(来自之前的博客,内容最后更新于 2017-9-16)

Previous Post Next Post

Blog Comments powered by Disqus.