离散傅里叶变换是一种快速进行多项式乘法的方法，朴素的多项式乘法时间复杂度一般为 $O(n^2)$ ,通过离散傅里叶变换，时间复杂度可以优化到 $O(n\log n)$

多项式的表示方法

系数表示法

我一直以来学习的多项式表示方法都是系数表示法，简单来说就是 $f_n(x)=a_0x^0+a_1x^1+a_2x^2……a_nx^n\\$ 其中 $f$ 是 $n$ 阶多项式

$g_m(x)=b_0x^0+b_1x^1+b_2x^2……b_mx^m\\$ 其中 $g$ 是 $m$ 阶多项式

那么对于这两个多项式的乘积，就是 $h_{n+m}(x)=\sum_{d=0}^{n+m}x^d\sum_{i=0}^{n}a_ib_{d-i}$
$*$ 高中知识，用 $f$ 中每一项乘 $g$ 中每一项，然后合并同类项，显然复杂度是 $O(n^2)$ 的

点值表示法

众所周知 $n+1$ 个点可以确定一个 $n$ 阶的多项式，所以一个 $n$ 阶多项式可以表示成 $f_n(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),……(x_3,f(x_3))\}\\$

那么 $f$ 和 $g$ 相乘就可以写成 $fg(x)=\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),(x_2,f(x_2)g(x_2)),……(x_n,f(x_n)g(x_n))\}\\$

只需要 $O(n)$ 的复杂度就可以完成

系数表示法转点值表示法

朴素的方法就是一个个算，每算一次就要 $O(n)$ 的复杂度，算 $n$ 次，需要 $O(n^2)$ 的复杂度，得不偿失
而离散傅里叶变换就是一种能够在 $O(n\log n)$ 的时间复杂度内将系数表示法转换为点值表示法的算法
考虑我们想要的其实就是能够快速得算出每一个点值的大小，如果 $x=\{1,-1\}$ 就可以方便得算出答案，但是显然实数范围内只有这两个特殊得数，对于阶数足够高的题，起不了多大的作用，所以我们就去复数范围内找答案

单位复根

在负数里，一个负数可以表示为复平面上的一个向量，两个复数相乘可以表示为向量的模相乘幅角相加
也就是说，如果有一个复数 $\omega_n$ 所表示的向量恰好落在以(0,0)为圆心的单位圆上，幅角为 $\frac{2\pi}{n}$ ,那么 $(\omega_n)^n$ 就是 $1$ 了。
同时 $((\omega_n)^2)^n$ 也是 $1$ ,以此类推， $\omega_n$ 的整数次幂的 $n$ 次幂
其中 $\omega_n$ 就被称为单位复根
图片说明
如图，红色的向量就是 $n=8$ 的单位复根，绿色的是 $n=4$ 的单位复根，蓝色的是 $n=2$ 的
单位复根可以表示为一个复数 $\omega_n=a+bi=\cos{\frac{2\pi}{n}}+\sin{\frac{2\pi}{n}}i$ 同时根据欧拉公式，又 $=e^{\frac{2\pi i}{n}}$

单位复根的性质

对于任意的正整数 $n$ 和整数 $k$
$\omega_n^n=1\\$

$\omega_n^{2k}=\omega_{\frac{n}{2}}^{k}\\$

$\omega_n^{k+\frac{n}{2}}=-\omega_n^k\\$

快速傅里叶变换

考虑一个 $8$ 项的多项式 $f(x)=a_0x^0+a_1x^1+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7\\$
按下标的奇偶划分 $f(x)=a_0x^0+a_2x^2+a_4x^4+a_6x^6+a_1x^1+a_3x^3+a_5x^5+a_7x^7\\$
对奇数项提取公因子 $x$
$f(x)=(a_0x^0+a_2x^2+a_4x^4+a_6x^6)+x(a_1x^0+a_3x^2+a_5x^4+a_7x^6)\\$
整理一下

$\\ \left.\begin{cases} f(x)=g(x^2)+xh(x^2)\\ g(x)=a_0x^0+a_2x^1+a_4x^2+a_6x^3\\ h(x)=a_1x^0+a_3x^1+a_5x^2+a_7x^3\\ \end{cases}\right.$
将 $x=\omega_n^k$ 带入 $\\ f(\omega_n^k)=g((\omega_n^{k})^2)+\omega_n^kh((\omega_n^{k})^2)\\ =g(\omega_n^{2k})+\omega_n^kh(\omega_n^{2k})\\ =g(\omega_{\frac{n}{2}}^{k})+\omega_n^kh(\omega_{\frac{n}{2}}^{k})\\$
可以发现 $g,h,f$ 非常相似，仅仅只是系数的不同，所以可以递归得分治来解决这个问题

我们假设 $f[i]$ 表示 $f_{[1,n]}(\omega_n^i)$ 的值
在我们求 $f[i]$ 时，先递归的求解 $f_{[1,\frac{n}{2}]}(\omega_\frac{n}{2}^i)$ 和 $f_{[\frac{n}{2}+1,n]}(\omega_\frac{n}{2}^i)$ 的值，则 $\\ f_{[1,n]}(\omega_n^i)=f_{[1,\frac{n}{2}]}(\omega_\frac{n}{2}^i)+\omega_n^i f_{[\frac{n}{2}+1,n]}(\omega_\frac{n}{2}^i)$
同时根据 $-\omega_n^i=\omega_n^{i+\frac{n}{2}}$
$f[i+\frac{n}{2}]=f_{[1,n]}(-\omega_n^i)\\ =f_{[1,\frac{n}{2}]}\left((-\omega_n^i)^2\right)-\omega_n^i f_{[\frac{n}{2}+1,n]}\left((-\omega_n^i)^2\right)\\ =f_{[1,\frac{n}{2}]}(\omega_\frac{n}{2}^i)-\omega_n^i f_{[\frac{n}{2}+1,n]}(\omega_\frac{n}{2}^i)$
所以我们只需要遍历一遍 $[1,\frac n2]$ 就可以得到所有的答案了

递归版离散傅里叶变换模板

void DFT(Comp f[], int n, int rev)//f[i]表示x=(w_n)^i的值
{
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  for (int i = 0; i < n; ++i) {  //偶数放左边，奇数放右边
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  DFT(g, n / 2, rev);
  DFT(h, n / 2, rev);  //递归DFT
  Comp cur(1, 0), wn(cos(2 * M_PI / n), sin(2 * M_PI * rev / n));
  for (int k = 0; k < n / 2; ++k) {
    tmp[k] = g[k] + cur * h[k];
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur = cur * wn;//只需要每次都乘w_n就可以得到每一个单位复根
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

蝴蝶变换

列表发现，序列中每一项在递归终点的位置为它下标二进制表示反转之后的值
例如在 $n=8$ 的离散傅里叶变换中， $1$ 的二进制表示是 $001_2$ ，在递归终点时，它的位置为 $4$ ,即 $100_2$
如下图所示
图片说明
我们可以动态规划来求解每个数转换后下标的位置

dp实现

记 $pos[i]$ 为原下标为 $i$ 的系数最终的位置
显然 $pos[0]=0$
对于任意的 $i$ ，他的二进制表示为 $[b_n……b_2b_1b_0]_2\\$

则 $\left\lfloor\frac i2\right\rfloor$ 的二进制表示为 $[0b_n……b_2b_1]_2\\$

翻转后，也就是 $pos[i]$ 的二进制表示为 $[b_0b_1b_2……b_n]_2\\$

则 $pos\left[\left\lfloor\frac i2\right\rfloor\right]$ 的二进制表示为 $[b_1b_2……b_n0]_2\\$

可以发现 $pos[i]=\left(pos\left[\left\lfloor\frac i2\right\rfloor\right]>>1\right)|(i\&1)<<\left(\left\lceil\log(i)\right\rceil\right)$

for (int i = 0; i < len; ++i)
{
    pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (l - 1));
}

于是fft的代码变为

void change(comp y[],int len)
{
    for (int i = 0; i < len; ++i)
    {
        if(i<pos[i])
            swap(y[i], y[pos[i]]);
    }
}

void DFT(comp y[],int len,int rev)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        comp wn(cos(2 * PI / h), sin(rev * 2 * PI / h));
        for (int j = 0; j < len; j += h)
        {
            comp w(1,0);
            for (int k = j; k < j + h / 2; k++)
            {
                comp u = y[k];
                comp t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w * wn;
            }
        }
    }
    if (rev == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
}

数学-离散傅里叶变换