算法导论30（多项式与快速傅里叶变换）

$$\begin{array}{*{20}{l}} {A(x) = {a_0} + {a_1}x + \cdots + {a_{n - 1}}{x^{n - 1}} = \sum\limits_{j = 0}^{n - 1} {{a_j}{x^j}} }\\ {({a_0} + {a_1}x + \cdots + {a_{n - 1}}{x^{n - 1}})({b_0} + {b_1}x + \cdots + {b_{n - 1}}{x^{n - 1}})}\\ { = {a_0}{b_0} + ({a_0}{b_1} + {a_1}{b_0})x + \cdots + {a_{n - 1}}{b_{n - 1}}{x^{2n - 2}}}\\ { = \sum\limits_{j = 0}^{2n - 2} {(\sum\limits_{k = 0}^j {{a_k}{b_{j - k}}} ){x^j}} ({a_j} = {b_j} = 0,j \ge n)}\\ \end{array}$$

30.1 多项式的表达

系数表达：$({a_0},{a_1}, \cdots ,{a_{n - 1}})$

点值表达：$\{ ({x_0},{y_0}),({x_1},{y_1}), \cdots ,({x_{n - 1}},{y_{n - 1}})\} ({y_k} = A({x_k}),0 \le k \le n - 1)$

求值（系数表达转换为点值表达，${\rm O}({n^2})$）：
霍纳准则：$A({x_0}) = {a_0} + {x_0}({a_1} + \cdots + {x_0}({a_{n - 2}} + {x_0}{a_{n - 1}}) \cdots )$

插值（点值表达转换为系数表达，${\rm O}({n^2})$）：
拉格朗日公式：$A(x) = \sum\limits_{k = 0}^{n - 1} {{y_k}} \frac{{\prod\limits_{j \ne k} {(x - {x_j})} }}{{\prod\limits_{j \ne k} {({x_k} - {x_j})} }}$

30.2 DFT与FFT

$$\begin{array}{l} {w^n} = 1\\ w = w_n^0,w_n^1, \cdots ,w_n^{n - 1}({w_n} = {e^{2\pi i/n}})\\ {y_k} = A(w_n^k) = \sum\limits_{j = 0}^{n - 1} {{a_j}w_n^{kj}} (0 \le k \le n - 1) \end{array}$$

$y = ({y_0},{y_1}, \cdots ,{y_{n - 1}})$是$a = ({a_0},{a_1}, \cdots ,{a_{n - 1}})$的离散傅里叶变换（DFT）。

$$\begin{array}{l} n = {2^k}\\ {A^{[0]}}(x) = {a_0} + {a_2}x + \cdots + {a_{n - 2}}{x^{n/2 - 1}}\\ {A^{[1]}}(x) = {a_1} + {a_3}x + \cdots + {a_{n - 1}}{x^{n/2 - 1}}\\ A(x) = {A^{[0]}}({x^2}) + x{A^{[1]}}({x^2}) \end{array}$$

$$\begin{array}{l} 0 \le k \le n/2 - 1\\ {y_k} = y_k^{[0]} + w_n^ky_k^{[1]}\\ = {A^{[0]}}(w_{n/2}^k) + w_n^k{A^{[1]}}(w_{n/2}^k)\\ = {A^{[0]}}(w_n^{2k}) + w_n^k{A^{[1]}}(w_n^{2k})\\ = A(w_n^k)\\ {y_{k + n/2}} = y_k^{[0]} - w_n^ky_k^{[1]}\\ = {A^{[0]}}(w_{n/2}^k) - w_n^k{A^{[1]}}(w_{n/2}^k)\\ = {A^{[0]}}(w_n^{2k + n}) + w_n^{k + n/2}{A^{[1]}}(w_n^{2k + n})\\ = A(w_n^{k + n/2}) \end{array}$$

const double PI=3.141592654;

class Complex
{
public:
    Complex(double x=0,double y=0):a(x),b(y){}
    friend Complex operator+(Complex x,Complex y)
    {
        return Complex(x.a+y.a,x.b+y.b);
    }
    friend Complex operator-(Complex x,Complex y)
    {
        return Complex(x.a-y.a,x.b-y.b);
    }
    friend Complex operator*(Complex x,Complex y)
    {
        return Complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
    }
    double a,b;
};

//递归
vector<Complex> recursiveFFT(vector<Complex> a)
{
    int n=a.size();
    if(n==1)return a;
    Complex wn(cos(2*PI/n),sin(2*PI/n)),w(1,0);
    vector<Complex> a0,a1,y0,y1,y;
    for(int i=0;i<n;i+=2)a0.push_back(a[i]);
    for(int i=1;i<n;i+=2)a1.push_back(a[i]);
    y0=recursiveFFT(a0);
    y1=recursiveFFT(a1);
    for(int k=0;k<=n/2-1;++k)
    {
        y[k]=y0[k]+w*y1[k];
        y[k+n/2]=y0[k]-w*y1[k];
        w=w*wn;
    }
    return y;
}

$T(n) = 2T(n/2) + \Theta (n) = \Theta (n\log n)$

30.3 高效FFT实现

//蝴蝶操作
for(int k=0;k<=n/2-1;++k)
{
    Complex t=w*y1[k];
    y[k]=y0[k]+t;
    y[k+n/2]=y0[k]-t;
    w=w*wn;
}

//位逆序置换
int rev(int k,int n)
{
    int k0=0,m=log((double)n)/log((double)2);
    for(int j=0;j<m;++j)
    {
        if(k&(1<<j))k0+=(1<<(m-1-j));
    }
    return k0;
}

void bitReverseCopy(vector<Complex> a,vector<Complex> &A)
{
    int n=a.size();
    for(int k=0;k<=n-1;++k)A[rev(k,n)]=a[k];
}

//迭代
vector<Complex> iterativeFFT(vector<Complex> a)
{
    vector<Complex> A;
    bitReverseCopy(a,A);
    int n=a.size();
    for(int s=1;(1<<s)<=n;++s)
    {
        int m=(1<<s);
        Complex wm(cos(2*PI/m),sin(2*PI/m));
        for(int k=0;k<=n-1;k+=m)
        {
            Complex w(1,0);
            for(int j=0;j<=m/2-1;++j)
            {
                Complex u=A[k+j],t=w*A[k+j+m/2];
                A[k+j]=u+t;
                A[k+j+m/2]=u-t;
                w=w*wm;
            }
        }
    }
    return A;
}