FFT

前置知识

两种多项式表示法：

系数表示法：$A(x)=\sum_{i=0}^{n-1}a_ix^i$ ，其中 $\{a_{n-1}\}$ 就是系数。

点值表示法：代入 $n$ 个不同的 $x$ 得到对应的 $A(x)$ 记为 $(x,y)$ ，这 $n$ 个点 $(x_0,y_0),(x_1,y_1)\cdots(x_{n-1},y_{n-1})$ 唯一确定了 $A(x)$ 。

由系数表示法变为点值表示法称为求值，点值表示法变为系数表示法称为插值。

点值表示法好处：快速求 $A(x)\times B(x)$ 的点值表达式——只需要将相同 $x$ 的 $y$ 乘起来。

单位复数根：

用 $n$ 个模为 $1$ 的复数 $n$ 等分复平面，令幅角最小的是 $\omega_n$ ( $n$ 次单位复数根)，那么其他 $n-1$ 个就是 $\omega_n^2,\omega_n^3\cdots\omega_n^n$ （两个复数相乘，模长相乘，幅角相加），画图yy一下可以发现 $\omega_n^0=\omega_n^n=1$ 。

因为是 $n$ 等分，所以 $\omega_n^k=(cos{2k\pi\over n},sin{2k\pi\over n})$ ，由此可见 $\omega_{2n}^{2k}=\omega_n^k,\omega_n^{k+{n\over 2}}=-\omega_n^k$ 。

快速傅里叶变换

因为点值表示法可以快速进行多项式乘法，所以我们想将系数表示法转换成点值表示法，再进行多项式相乘，最后将点值表示法转换成系数表示法。两个过程分别为离散傅里叶变换和傅里叶逆变换。

ps：接下来都默认 $n$ 是 $2^k$ （如果不够那么高位补 $0$ 即可）。

离散傅里叶变换：

如果将 $A(x)$ 的系数奇偶拆分，我们可以得到：

$$A_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{{n\over 2}-1}\\ A_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{{n\over 2}-1}\\ A(x)=A_0(x^2)+xA_1(x^2)$$
因为要代入 $n$ 个不同的 $x$ ，我们会联想到 $n$ 次单位复数根，又因为如果 $n$ 是偶数，则 $(\omega_n^{k})^2=\omega_{n\over 2}^{k}$ ，也就是说 $(\omega_{n}^{k})^2$ 的取值只有 $n\over 2$ 种。于是我们惊奇的发现只要代入 $n$ 个 $n$ 次单位复数根， $A_0(x^2)$ 和 $A_1(x^2)$ 就变成了与 $A(x)$ 同样的问题（这就是为什么 $n$ 要是 $2^k$ ）！

所以可以用 $O(nlog_2n)$ 的复杂度来愉快的离散傅里叶变换。

傅里叶逆变换：

只需要用 $\omega_n^{-k}$ 再做一次离散傅里叶变换，证明可以参考Menci大佬的博客QAQ，我太弱了就只记了结论：

$$a_k={1\over n}\sum_{i=0}^{n-1}a_i(\omega_n^{-k})^{i}$$

实现

离散傅里叶变换的递归写法很显然，但是常数巨大啊。实际上有迭代写法，我们观察每个 递归到最后的位置：

0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7
二进制：
 0   1   2   3   4   5   6   7
000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111
 0   4   2   6   1   5   3   7

惊奇的发现竟然是二进制反转QAQ，这样我们就可以从底向上迭代处理了。

模板

UOJ34，烂大街的FFT模板题。

#include<cstdio>
#include<cctype>
#include<cmath>
#include<algorithm>
#define fr first
#define sc second
#define mp make_pair
using namespace std;
typedef long double DB;typedef pair<DB,DB> C;
const int maxn=262144;const DB pi=acos(-1);

int n,m,R[maxn+5];C a[maxn+5],b[maxn+5];

#define Eoln(x) ((x)==10||(x)==13||(x)==EOF)
inline char readc(){
    static char buf[100000],*l=buf,*r=buf;
    if (l==r) r=(l=buf)+fread(buf,1,100000,stdin);
    if (l==r) return EOF;return *l++;
}
inline int readi(int &x){
    int tot=0,f=1;char ch=readc(),lst='+';
    while (!isdigit(ch)) {if (ch==EOF) return EOF;lst=ch;ch=readc();}
    if (lst=='-') f=-f;
    while (isdigit(ch)) tot=(tot<<3)+(tot<<1)+ch-48,ch=readc();
    return x=tot*f,Eoln(ch);
}
inline int Rev(int x,int len){
    static int buf[31];for (int i=0;i<len;i++) buf[i]=x&1,x>>=1;
    for (int i=0;i<len;i++) x=x<<1|buf[i];return x;
}
C operator + (const C &a,const C &b) {return mp(a.fr+b.fr,a.sc+b.sc);}
C operator - (const C &a,const C &b) {return mp(a.fr-b.fr,a.sc-b.sc);}
C operator * (const C &a,const C &b) {return mp(a.fr*b.fr-a.sc*b.sc,a.fr*b.sc+a.sc*b.fr);}
inline void FFT(C *a,int n,int f){
    for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1){
        C w=mp(1,0),wn=mp(cos(pi/k),sin(f*pi/k)),x,y;
        for (int i=0;i<n;i+=(k<<1),w=mp(1,0))
            for (int j=0;j<k;j++,w=w*wn)
                x=a[i+j],y=w*a[i+j+k],a[i+j]=x+y,a[i+j+k]=x-y;
    }
}
int main(){
    freopen("FFT.in","r",stdin);
    freopen("FFT.out","w",stdout);
    readi(n);readi(m);
    for (int i=0,x;i<=n;i++) readi(x),a[i]=mp(x,0);
    for (int i=0,x;i<=m;i++) readi(x),b[i]=mp(x,0);
    m+=n;int len=0;for (n=1;n<=m;n<<=1) len++;for (int i=0;i<n;i++) R[i]=Rev(i,len);
    FFT(a,n,1);FFT(b,n,1);for (int i=0;i<n;i++) a[i]=a[i]*b[i];
    FFT(a,n,-1);for (int i=0;i<=m;i++) printf("%d ",(int)(a[i].fr/n+0.1));
    return 0;
}