FFT 模板

$\DeclareMathOperator{\rev}{rev}$

FFT（Fast Fourier Transform），确切地说应该称之为 FDFT（Fast Discrete Fourier Transform），因为 FFT 是为了解决 DFT 问题而设计的一种快速算法。在深入讨论之前，有必要特别指出这一点。

DFT 问题：

给定复数域上的 $n-1$ 次多项式 $A(x)$ 的系数表示（coefficient representation）$(a_0, a_1,\dots, a_{n-1})$，求 $A(x)$ 的某个点-值表示（point-value representation）：

\[((x_0, y_0), (x_1, y_1), (x_2, y_2), \dots, (x_{n-1}, y_{n-1}))\]

 

 

FFT 的非递归（迭代）实现 Version I 手写 Complex 类（《算法导论》）

 

#include <bits/stdc++.h>
#define rep(i, l, r) for(int i=l; i<r; i++)
using namespace std;
const double PI(acos(-1));

struct Complex{
    double r, i;
    Complex(double r, double i):r(r), i(i){}
    Complex(int n):r(cos(2*PI/n)), i(sin(2*PI/n)){}    //!!error-prone
    Complex():r(0), i(0){}    //default constructor
    Complex &operator*=(const Complex &a){
        double R=r*a.r-i*a.i, I=r*a.i+a.r*i;
        r=R, i=I;
        return *this;
    }
    Complex operator+(const Complex a){
        return Complex(r+a.r, i+a.i);
    }
    Complex operator-(const Complex a){
        return Complex(r-a.r, i-a.i);
    }
    Complex operator*(const Complex a){
        return Complex(r*a.r-i*a.i, r*a.i+a.r*i);
    }
    void out(){
        cout<<r<<' '<<i<<endl;
    }
};

const int N(1<<17);
int ans[N];
Complex a[N], b[N];
char s[N], t[N];

void bit_reverse_swap(Complex *a, int n){
    for(int i=1, j=n>>1, k; i<n-1; i++){
        if(i < j) swap(a[i],a[j]);
        //tricky
        for(k=n>>1; j>=k; j-=k, k>>=1);    //inspect the highest "1"
        j+=k;
    }
}

void FFT(Complex* a, int n, int t){
    bit_reverse_swap(a, n);
    for(int i=2; i<=n; i<<=1){
        Complex wi(i*t);
        for(int j=0; j<n; j+=i){
            Complex w(1, 0);
            for(int k=j, h=i>>1; k<j+h; k++){
                Complex t=w*a[k+h], u=a[k];
                a[k]=u+t;
                a[k+h]=u-t;
                w*=wi;
            }
        }
    }
    if(t==-1) rep(i, 0, n) a[i].r/=n;    //!!error-prone
}

int trans(int x){
    int i=0;
    for(; x>1<<i; i++);
    return 1<<i;
}

int main(){
    for(; ~scanf("%s%s", s, t); ){
        int n=strlen(s), m=strlen(t), l=trans(n+m-1);
        rep(i, 0, n) a[i]=Complex(s[n-1-i]-'0', 0);
        rep(i, n, l) a[i]=Complex(0, 0);
        rep(i, 0, m) b[i]=Complex(t[m-1-i]-'0', 0);
        rep(i, m, l) b[i]=Complex(0, 0);

        FFT(a, l, 1), FFT(b, l, 1);
        rep(i, 0, l) a[i]*=b[i];
        FFT(a, l, -1);
        rep(i, 0, l) ans[i]=(int)(a[i].r+0.5); ans[l]=0;    //error-prone
        rep(i, 0, l) ans[i+1]+=ans[i]/10, ans[i]%=10;
        int c=l;
        for(;c && !ans[c]; --c);
        for(; ~c; putchar(ans[c--]+'0'));    //error-prone
        puts("");
    }
    return 0;
}

 

 

 

Comment:

1. 此代码是为 HDU1402 写的。代码中，凡注释 error-prone 处，都应特别小心。我犯的最傻逼的错误是第9行，应当是2*PI，我写成PI了。

2. FFT的数值稳定性（精度）问题，还有待考虑。(UPD)多次做多项式乘法时，精度损失较快，这时将 double 换成 long double 可缓解精度损失。

---

Comment:

1. bit_reverse_swap()函数是对算法导论上的bit_reverse_copy()的改进，将下标互为bit-reverse 的两元素互换位置，就免去了 copy 所需的空间。

2.bit_revrese_copy()不太好懂，需要一点解释：

void bit_reverse_swap(Complex *a, int n){
    for(int i=1, j=n>>1, k; i<n-1; i++){
        if(i < j) swap(a[i],a[j]);
        //tricky
        for(k=n>>1; j>=k; j-=k, k>>=1);    //inspect the highest "1"
        j+=k;
    }
}

将$i$的bit-reverse记作$\rev(i)$。

(i). 由于 $\rev(0)=1, \rev(n-1)=n-1$（$n$ 是 $2$ 的幂），所以第 2 行的主循环可令 $i$ 从 $1$ 循环到 $n-2$。同时 $j$ 从 $\rev(1)= \frac{n}{2}$，“循环”到 $\rev(n-2)$ 。

(ii). 第 3 行的判断 if(i < j) 避免了重复交换

(iii).第 5 行的循环的作用就是将 $j$ 从 $\rev(i)$ 变成 $\rev(i+1)$：

首先应当注意到，$i$ 的最低位恰是 $rev(i)$ 的最高位。若 $\rev(i)$ 的最高位是 $0$ 那么 $\rev(i+1)$ 就是 $\rev(i)+ \frac{n}{2}$，否则，$i$ 加上 $1$ 后，最低位将变成 $0$，并且向高一位进 $1$ 。相应的，$\rev(i+1)$ 的最高位应置 $0$（即代码中的 j-=k），并且向低一位"进“ $1$（对应代码中的 k>>=1）。这样从高位往低位检查，遇到 $1$（对应代码中的条件 j>=k）就进位，遇到 $0$ 就退出循环。
3. 我写代码时把第 58 行的 == 写成了 =，结果DEBUG 一个多小时。。。

Version II: 用 C++ 标准库中的 complex<double> 类，代码短一些，但也会慢一些：

#include <bits/stdc++.h>
#define rep(i, l, r) for(int i=l; i<r; i++)
using namespace std;
const double PI(acos(-1));
typedef complex<double> C;

const int N(1<<17);
int ans[N];
C a[N], b[N];
char s[N], t[N];

void bit_reverse_swap(C *a, int n){
    for(int i=1, j=n>>1, k; i<n-1; i++){
        if(i < j) swap(a[i],a[j]);
        //tricky
        for(k=n>>1; j>=k; j-=k, k>>=1);    //inspect the highest "1"
        j+=k;
    }
}

void FFT(C* a, int n, int t){
    bit_reverse_swap(a, n);
    for(int i=2; i<=n; i<<=1){
        C wi(cos(t*2*PI/i), sin(t*2*PI/i));
        for(int j=0; j<n; j+=i){
            C w(1);
            for(int k=j, h=i>>1; k<j+h; k++){
                C t=w*a[k+h], u=a[k];
                a[k]=u+t;
                a[k+h]=u-t;
                w*=wi;
            }
        }
    }
    if(t==-1) rep(i, 0, n) a[i]/=n;    //!!error-prone: typo ==/=
}

int trans(int x){
    int i=0;
    for(; x>1<<i; i++);
    return 1<<i;
}

int main(){
    for(; ~scanf("%s%s", s, t); ){
        int n=strlen(s), m=strlen(t), l=trans(n+m-1);
        rep(i, 0, n) a[i]=C(s[n-1-i]-'0');
        rep(i, n, l) a[i]=C(0);
        rep(i, 0, m) b[i]=C(t[m-1-i]-'0');
        rep(i, m, l) b[i]=C(0);

        FFT(a, l, 1), FFT(b, l, 1);
        rep(i, 0, l) a[i]*=b[i];
        FFT(a, l, -1);
        rep(i, 0, l) ans[i]=(int)(a[i].real()+0.5); ans[l]=0;    //error-prone
        rep(i, 0, l) ans[i+1]+=ans[i]/10, ans[i]%=10;
        int p=l;
        for(;p && !ans[p]; --p);
        for(; ~p; putchar(ans[p--]+'0'));    //error-prone
        puts("");
    }
    return 0;
}