FFT

$$\rm FFT$$

好像全机房除了我以外都精通 $\rm FFT$ ，$\rm QAQ$

前言

$\rm FFT$：快速傅里叶变换，是用来做多项式乘法或者加法卷积以及其他运算的一种 $\mathcal O(n \log n)$ 的方法
$\rm FNTT及NTT$：快速傅里叶变换的优化版—>优化常数及误差
$\rm FWT$：快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题
$\rm MTT$：毛爷爷的FFT—>任意模数
$\rm FMT$：快速莫比乌斯变化-> ~~可被暴力替代~~

参考资料：
https://blog.csdn.net/enjoy_pascal/article/details/81478582
https://www.cnblogs.com/RabbitHu/p/FFT.html
https://www.luogu.com.cn/blog/Soulist/post-fft
https://www.luogu.com.cn/blog/tbr-blog/duo-xiang-shi-xue-xi-bi-ji

一、前置芝士

$\quad (Ⅰ)$ 多项式相关知识

$\quad\quad\qquad$ 多项式的定义：由若干个单项式相加组成的代数式叫做多项式。——百度百科

$\quad\quad (1)$ 多项式的表示

$\quad\quad\qquad$ 多项式有以下两种表示方法：系数表示法和点值表示法。

$\quad\qquad$系数表示法

$\qquad\qquad\quad$ 一个 $n-1$ 次 $n$ 项式可表示为 $f(x)=\sum_{i=0}^{n-1}A_ix^i$ 。
$\qquad\qquad\quad$ 于是用系数来表示 $f(x)$ 即 $f(x)=\{A_0,A_1,A_2...A_{n-1}\}$

$\qquad\quad$点值表示法

$\qquad\qquad\quad$ 把多项式看成一个函数，放到平面直角坐标系中，把 $n$ 个不同的 $x$ 代入函数中得到 $n$ 个 $y$ ，于是得到 $n$ 个不同的点。显然可以发现，把各系数堪称未知数则相当于一个 $n$ 元一次方程组，所以 $n$ 个不同的点可以唯一确定一个多项式
$\qquad\qquad\quad$ $f(x)=\{(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})\}$ ，这就是多项式的点值表示法。

$\quad\quad(2)$ 多项式乘法

$\quad\quad\qquad$对于两种不同的表示方法，多项式乘法的计算方式也有不同。
$\quad\quad\qquad$对于系数表示法来说，也就是我们平常数学上列的式子，设 $A(k),B(k),C(k)$ 分别是三个多项式 $f,g,h$ 的 $x^k$ 项的系数，可知：

\[C(k)=\sum_{i+j=k}A(i)*B(j) \]

$\quad\quad\qquad$这是很好理解的，因为 $h$ 的 $x^k$ 项的系数应当是由 $f,g$ 的所有次数和为 $k$ 的系数对相乘得来的。
$\quad\quad\qquad$而对于点值表示法，则：

\[h(x)=\{(x_0,f(x_0)*g(x_0)),(x_1,f(x_1) *g(x_1))...(x_{n-1},f(x_{n-1}) *g(x_{n-1})) \} \]

$\quad\quad\qquad$我们可以发现，如果用系数表示法来求多项式乘法的话，复杂度是 $\mathcal O(n^2)$ 的，但是点值表示法的复杂度是 $\mathcal O(n)$ 的。
$\quad\quad\qquad$这就是 $\rm FFT$ 的核心思想了：先将系数表示的两个多项式转化成点值表示法，将点值相乘后再转回系数表示法。其中，系数表示 $\to$ 点值表示的过程叫做 $\rm DFT$ ，点值表示 $\to$ 系数表示的过程叫做 $\rm IDFT$ （逆 $\rm DFT$ ）。
$\quad\quad\qquad$但是将系数表示法转化成点值表示法的复杂度也是 $\mathcal O(n^2)$ 的，具体做法是随意地选出 $n$ 个 $x$ 依次代入。这是因为没有选到合适的代入值。 $\rm FFT$ 就是用一种很巧妙的值来代入实现 $O(n \log n)$ 转化的。

$\quad (Ⅱ)$ 复数

$\quad\quad (1)$ 初识复数

$\quad\quad\qquad$初中我们学过实数，那时候觉得好像所有数都是实数，完全无法想象虚数是什么样的。
$\quad\quad\qquad$定义 $i^2=-1$ ，这在实数上是不可能成立的，但是它是虚数的定义之类的东西。其中， $i$ 称为虚数单位。这样，我们可以用 $b*i$ 的形式来表示一个虚数。那么现在负数的平方根我们也能表示了， $\sqrt{-7}=\sqrt{7}i$ 。
$\quad\quad\qquad$复数就是既有实数部分又有虚数部分的数，通常表示为 $a+b*i$ 。其中 $a$ 称为实部， $b$ 称为虚部。这样的一个形式像不像一个函数？~~（像个鬼）（你说像就像好吧）~~ 所以我们可以像平面直角坐标系一样表示复数。

$\quad\quad\qquad$~~（很显然这张图是 kuai 的）~~
$\quad\quad\qquad$这样的横轴表示实部，纵轴表示虚部的坐标系我们把它叫做复平面。这样一个复数就可以表示为复平面上的一个点。比如说图上那个点表示的就是 $2+3i$ 。而这个复数的模就是连接点到原点的线段的长度即 $\sqrt{13}$ 。
$\quad\quad\qquad$（感觉讲得不是很好，如果有需要可以自己去查找一下其他资料。）

$\quad\quad (2)$ 复数的运算

$\quad\quad\qquad$既然是数，那么复数肯定也可以四则运算。~~（似乎挺对的）~~
$\quad\quad\qquad$比如说两个复数 $z_1=a+b_i,z_2=c+d_i$ ，那么 $z_1+z_2=(a+c)+(b+d)i$ （实部虚部分别相加）， $z_1*z_2=(ac-bd)+(ad+bc)i$ （各项分别相乘）。
$\quad\quad\qquad$表示在复平面上，复数相加也满足平行四边形法则（这个自己动手画一画吧）~~（ ta 某曾经曰过：“我不是懒，只是怕东西太多页面卡”）~~
$\quad\quad\qquad$另外还有一个很重要的性质：复数相乘等效于复平面上表示复数的两条线段模长相乘，幅角相加。（幅角就是从横轴正半轴转到这条线段所经过的角度。）这个性质在下面会要用到。

$\quad\quad (3)$ 简介单位根

$\quad\quad\qquad$这就是之前所说的很妙的值。

$\quad\quad\qquad$数学上， $\rm n$ 次单位根是 $\rm n$ 次幂为 $\rm 1$ 的复数。它们位于复平面的单位圆上，构成正 $n$ 边形的顶点，其中一个顶点是 $1$ 。——百度百科

$\quad\quad\qquad$我们在复平面上画一个单位圆，圆上点表示的数都可以经过若干次方得到 $\rm 1$ （这是因为复数相乘模长相乘辐角相加）。同时根据这一个性质，圆上所有的点都可以表示成一个复数的次方。这个复数就是单位根， $\rm n$ 次单位根表示为 $\omega_{n}^{1}$ 。
$\quad\quad\qquad$由于我们要取 $n$ 个数来求点值，所以我们可以将单位圆分成均等的 $n$ 份，然后取这 $n$ 个数。

$\quad\quad\qquad$就比如说这张图中，数字标号 $k$ 分别表示这个点是 $\omega_{n}^{k}$ 。同时，很显然每个 $\omega$ 都可以通过三角函数求出来的，$\omega_{n}^{k}=cos^k(\frac{2\pi}{n})+sin^k(\frac{2\pi}{n})*i$。
$\quad\quad\qquad$在信息学中，我们通常用一个结构体来表示复数，分别存储实部和虚部。

    struct node {
        double x, y;
    }a[maxn], b[maxn];
    il node operator + (node a, node b) {return (node){a.x + b.x, a.y + b.y};}
    il node operator - (node a, node b) {return (node){a.x - b.x, a.y - b.y};}
    il node operator * (node a, node b) {return (node){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}

$\quad\quad (4)$ 单位根的性质

$\quad\qquad ①$ $\omega_{2n}^{2k}=\omega_{n}^{k}$

$\quad\quad\qquad$这个性质在图上挺好理解的。划成 $2n$ 份的单位圆上的第 $2k$ 个点和划成 $n$ 份的单位圆上的第 $k$ 个点显然就是同一个点了。

$\quad\qquad ②$ $\omega_{n}^{k+n/2}=-\omega_{n}^{k}$

$\quad\quad\qquad$这个性质也可以在图上得到很好的体现。

$\quad\qquad ③$ $\sum_{j=0}^{n-1} (\omega_{n}^{k})^j=n*[k==0]$

$\quad\quad\qquad$若 $k=0$ ，则显然答案为 $n$。
$\quad\quad\qquad$若 $k \ne 0$ ，设

\[S=\sum_{j=0}^{n-1} (\omega_{n}^{k})^j \]

$\quad\quad\qquad$根据等必须列求和的方法，则

\[\omega_{n}^{k}*S=\sum_{j=1}^{n} (\omega_{n}^{k})^j \]

$\quad\quad\qquad$两式相减，则

\[S=\dfrac{(\omega_{n}^{k})^n-(\omega_{n}^{k})^0}{\omega_{n}^{k}-1} \]

$\quad\quad\qquad$显然可见，分母为 $0$ 分子不为 $0$ 。

$\quad\qquad$这三个性质在后面都很有用。

二、正文

$\quad$以下均保证 $n$ 为 $2$ 的次幂。

$\quad\quad \rm DFT$

$\quad\qquad$对于

\[F(x)=A_0*x^0+A_1*x^1+...+A_{n-1}*x^{n-1} \]

$\quad\qquad$我们按照下标奇偶分成

\[F1(x)=A_0*x^0+A_2 *x^1+...+A_{n-2} *x^{\frac{n}{2}-1} \]

\[F2(x)=A_1*x^0+A_3 *x^1+...+A_{n-1} *x^{\frac{n}{2}-1} \]

$\quad\qquad$那么：

\[F(x)=F1(x^2)+x *F2(x^2) \]

$\quad\qquad$把单位根代入：

\[F(\omega_{n}^{x})=F1(\omega_{n}^{2x})+\omega_{n}^{x}*F2(\omega_{n}^{2x}) \]

$\quad\qquad$根据性质 ①

\[F(\omega_{n}^{x})=F1(\omega_{n/2}^{x})+\omega_{n}^{x}*F2(\omega_{n/2}^{x}) \]

$\quad\qquad$当 $x < n/2$ 时，我们可以直接用 $F1,F2$ 的 $\frac{n}{2}$ 个值就能求出 $F$ 的 $n$ 个值。
$\quad\qquad$当 $x \ge n/2$ 时，则设 $x=x'+n/2$

\[F(\omega_{n}^{x})=F1((\omega_{n}^{x'+n/2})^2)+\omega_{n}^{x'+n/2}*F2((\omega_{n/2}^{x'+n/2})^2) \]

\[F(\omega_{n}^{x})=F1((-\omega_{n}^{x'})^2)-\omega_{n}^{x'}*F2((-\omega_{n}^{x'})^2) \]

\[F(\omega_{n}^{x})=F1(\omega_{n/2}^{x'})-\omega_{n}^{x'}*F2(\omega_{n/2}^{x'}) \]

$\quad\qquad$其实和 $k < n/2$ 的情况差不多。
$\quad\qquad$于是我们只需要知道 $F1,F2$ 的 $\frac{n}{2}$ 个值就可以求出 $F$ 的 $n$ 个值。
$\quad\qquad$复杂度 $T(x)=x+2*T(n/2)=x \log n$

$\quad\quad \rm IDFT$

$\quad\qquad$现在我们已知了 $n$ 个点的点值 $a_i$，求系数 $A_i$ 。

\[a_0=A_0*(\omega_{n}^{0})^0+A_1*(\omega_{n}^{0})^1+ \cdots+(\omega_{n}^{0})^{n-1} \]

\[a_1=A_0*(\omega_{n}^{1})^0+A_1*(\omega_{n}^{1})^1+ \cdots+(\omega_{n}^{1})^{n-1} \]

\[\vdots \]

\[a_{n-1}=A_{n-1}*(\omega_{n}^{n-1})^0+A_1*(\omega_{n}^{n-1})^1+ \cdots+(\omega_{n}^{n-1})^{n-1} \]

$\quad\qquad$其实相当于一个矩阵形式。

\[\left[ \begin{matrix} A_0 \\ A_1 \\ \vdots \\ A_{n-1}\\ \end{matrix} \right]* \left[ \begin{matrix} (\omega_{n}^{0})^0&(\omega_{n}^{0})^1 & \cdots&(\omega_{n}^{0})^{n-1}\\ (\omega_{n}^{1})^0&(\omega_{n}^{1})^1 & \cdots&(\omega_{n}^{1})^{n-1}\\ \vdots& \vdots & \vdots & \vdots\\ (\omega_{n}^{n-1})^0&(\omega_{n}^{n-1})^1 & \cdots&(\omega_{n}^{n-1})^{n-1}\\ \end{matrix} \right]= \left[ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\ \end{matrix} \right] \]

$\quad\qquad$设三个矩阵分别为 $A,V,a$ ，那么我们要求 $A$ 相当于求 $a*V^{-1}$ 。
$\quad\qquad$令 $V^{-1}(i,j)=(\omega_{n}^{-i})^j$ ，可知 $V*V^{-1}(i,j)=\sum_{k=0}^{n-1}(\omega_{n}^{j-i})^k$。根据性质三，当且仅当 $i=j$ 的时候 $V*V^{-1}(i,j)$ 不为 $0$ ，且其一定为 $n$ 。所以 $V$ 的逆 $V_{-1}$ 实际是 $\frac{1}{n}(\omega_{n}^{-i})^j$
$\quad\qquad$然后我们发现 $A=a*V^{-1}$ 其实可以看成 $a'=A'*V'$ ，即将 $A$ 视为新的点值矩阵， $V^{-1}$ 视为新的单位根矩阵， $a$ 视为新的系数矩阵，就可以再套一次 $\rm DFT$ 的板子。
$\quad\qquad \rm Warning:$ 最后的答案一定要除以 $n$ 即 $lim$;

Code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define Mid ((l+r)/2)
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define drep(i,a,b) for(int i=(a);i>=(b);--i)
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout);
const int maxn=1e5+5,mod=1e9+7,inf=0x3f3f3f3f;
int n,m,Q,K,T;
int read(){
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
const double pi=acos(-1.0);
struct node{double x,y;}a[maxn],b[maxn];
node operator + (node a,node b){return (node){a.x+b.x,a.y+b.y};}//定义复数运算
node operator - (node a,node b){return (node){a.x-b.x,a.y-b.y};}
node operator * (node a,node b){return (node){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
int lim;
void FFT(node *x,int f,int len){
    if(len==1)return;
    node f1[len>>1],f2[len>>1];
    for(int i=0;i<len;i+=2)f1[i/2]=x[i],f2[i/2]=x[i+1];
    FFT(f1,f,(len>>1)),FFT(f2,f,(len>>1));
    node w=(node){1,0},bas=(node){cos(2.0*pi/len),f*sin(2.0*pi/len)};
    for(int i=0;i<(len>>1);++i,w=w*bas){
        x[i]=f1[i]+w*f2[i],x[i+(len>>1)]=f1[i]-w*f2[i];
    }return;
}
signed main(){
    //file(a);
    n=read();m=read();
    rep(i,0,n)a[i].x=read();
    rep(i,0,m)b[i].x=read();
    while((1<<lim)<=n+m)++lim;//保证长度是 2 的幂次，否则二分不完整无法 FFT 。
    //如果 n+m 不是 2 的幂次，可以在后面补 0 
    FFT(a,1,(1<<lim));FFT(b,1,(1<<lim));//处理出点值
    rep(i,0,(1<<lim))a[i]=a[i]*b[i];//点值相乘
    FFT(a,-1,(1<<lim));//IDFT，转回系数表示
    rep(i,0,n+m)printf("%d ",int(a[i].x/(1<<lim)+0.5));//向上取整
    puts("");

    return 0;
}
</details>

可以发现，多项式乘法的系数表示法和加法卷积是一样的，所以 $\rm FFT$ 也可以用来做加法卷积。

posted @ 2020-07-07 19:49 maskey 阅读(142) 评论(0) 收藏举报

刷新页面返回顶部

maskey

FFT

$$\rm FFT$$

前言

一、前置芝士

\(\quad (Ⅰ)\) 多项式相关知识

\(\quad\quad (1)\) 多项式的表示

\(\quad\qquad\)系数表示法

\(\qquad\quad\)点值表示法

\(\quad\quad(2)\) 多项式乘法

\(\quad (Ⅱ)\) 复数

\(\quad\quad (1)\) 初识复数

\(\quad\quad (2)\) 复数的运算

\(\quad\quad (3)\) 简介单位根

\(\quad\quad (4)\) 单位根的性质

二、正文

\(\quad\quad \rm DFT\)

\(\quad\quad \rm IDFT\)

公告