傅里叶变换(FFT)学习笔记

——傅里叶变换在信息学竞赛中主要作用是利用卷积思想,化乘为加,快速计算多项式乘法。 我太蒟了,看了$F(x)$篇博文,还是没看懂。 关于多项式,有了$O(nlogn)$乘法,就有了全世(jia)界(tong)! $update(2019.6.27):$更新了一些证明和式子。 $update(2019.10.28):$发现模板题时限缩小,开始修锅(去冗余)。 当年写这篇文章的时候比较菜(现在仍然很菜),所以讲的比较含混不清,现在更新和精简了一些内容,希望能对大家有所帮助~ 源码删去了17KB,更新了9KB内容。 文章中的代码已经进一步优化,可以通过模板题。 由于文章前后内容有所变化,已将讨论区清空。 ------------ **警告** : 本文对数学水平有一定要求,如果发现看不懂,可以先存着。 另外,如果你掌握**和式**的话,学习速度将会大幅度提高。 此外,由于前面的内容太掉逼格,**建议水平高的选手直接跳到“单位根”**。 -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # -1.前置知识 1.线段树(或者基本分治)：[题目](https://www.luogu.org/problemnew/show/P3372) 2.基本数论：[题目](https://www.luogu.org/problemnew/show/P3811) 3.码力：[题目](https://www.luogu.org/problemnew/show/P1074) 4.基本数学：[题目](https://www.luogu.org/problemnew/show/P2261) 然后在你刷完这些题之后,发现题和下面的内容并没什么关系。 附上一句,这是$\color{purple}\textsf{省选}$内容哦！(然而并没有什么影响) $\color{Black}\colorbox{lightgreen}{总结一下}$ 这里不会满屏都是三角函数!!,没有$e^i$和什么欧拉定理!! 不需要你会高数,只要用过**平面直角坐标系**,就可以看。 ~~这篇文章源码**40KB**,写的时候挑战洛谷Blog(浏览器)性能极限~~ -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 0.概(che)论(dan) **FFT**全称(Fast Fourier Transformation) 中文名：**快速傅里叶离散变换** 作用 : 以$O(nlogn)$的复杂度计算多项式乘法(你以为呢?`Of course more than that!`)。 大家在学多项式乘法的时候,是不是觉得拆括号特别麻烦。 比如说计算$(x^2+3x-1)*(2x^2+x-5)$ - 原式$=x^2(2x^2+x-5)+3x(2x^2+x-5)-(2x^2+x-5)$ $\quad\ \ \ =(2x^4+x^3-5x^2)+(6x^3+3x^2-15x)-(2x^2+x-5)$ $\quad\ \ \ =2x^4+x^3-5x^2+6x^3+3x^2-15x-2x^2-x+5$ $\quad\ \ \ =2x^4+(x^3+6x^3)+(3x^2-2x^2-5x^2)+(-x-15x)+5$ $\quad\ \ \ =2x^4+7x^3-4x^2-16x+5$ 很复杂... 现在我们就要用计算机**把两个多项式乘起来**(所谓多项式问题)。 **P.S:**下面提及的多项式均**只有一元**,你可以认为只有$x$这一元。 为了表述方便,下面定义一些**记号**(请仔细阅读)： - $F(x)$表示一个多项式,简单的叫做“多项式$F$”。 $F(x)$是简便写法,比如我们设$F(x)=x^2+x+1$,以后我们提起$x^2+x+1$就不用那么啰嗦,直接说$F(x)$就好啦。 你也可以把它理解为函数,上面的$F(x)=x^2+x+1$ 那么$F(5)=5^2+5+1=31$,非常和谐。 - 系数(参数) 这里有一个**n次多项式**$F(x)$ 差不多长成这样：$a*x^n+b*x^{n-1}+c*x^{n-2}+...$ 其中$a,b,c...$是**参数**,因为他们在系数的位置上,所以也叫系数. 用不同的字母来表示系数很烦,我们就用数组。 通常把$F(x)$的$n$次项系数称作$F[n]$。 也即$F(x)=\sum^{n}_{i=0}F[i]x^i$ - 乘法的本质(卷积) 现在如果让你把两个$n$次多项式$F(x)$和$P(x)$乘起来,你会怎么写程序？ 简单啦！ 设$C=A*B$(多项式). $C[k]=\sum\limits_{i=0}^kA[i]B[k-i]$ **理解**:想要乘出$x^k$,需要$A$的$x^i$项和$B$的$x^{k-i}$项。 也就是说$A[i]B[k-i]$乘起来之后,他们后面的未知数就变成了$x^k$,所以要往$C[k]$贡献。 也可以写做等价的$C[k]=\sum\limits_{i+j==k}A[i]B[j]$(不懂没关系,看代码你就懂了) 形为$C[k]=\sum\limits_{i⊕j==k}A[i]B[j]$的式子称为卷积,其中$⊕$为某种运算。 那么你可能观察到了,**多项式乘法就是加法卷积**。 (后面我们还会学习到其他的卷积) 多项式乘法拥有交换律结合律分配率什么的,就不多说了…… 亮出[模板题](https://www.luogu.org/problemnew/show/P3803) 暴力卷积Code(用于上述模板题,44'): ```cpp #include<algorithm> #include<cstdio> #define Maxn 1000500 using namespace std; inline int read() { register int X=0; register char ch=0; while(ch<48||ch>57)ch=getchar(); while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar(); return X; } int n,m; long long f[Maxn],g[Maxn],s[Maxn]; void mul(long long *s,long long *f,long long *g) { for (int k=0;k<n+m-1;k++) for (int i=0;i<=k;i++) s[k]+=f[i]*g[k-i]; } int main() { scanf("%d%d",&n,&m);n++;m++; for (int i=0;i<n;i++)f[i]=read(); for (int i=0;i<m;i++)g[i]=read(); mul(s,f,g); for (int i=0;i<n+m-1;i++)printf("%lld ",s[i]); return 0; } ``` $\color{Black}\colorbox{lightgreen}{总结一下}$ 要用计算机**把两个多项式乘起来**就是所谓多项式问题。 数学家为了偷懒发明了**麻烦的记号**。 $\color{Black}\colorbox{yellow}{习题}$ - 验算上面的式子,你同意作者的化简过程吗？ -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 1.DFT & IDFT的思想 首先说一下多项式的**点值表达**。 如果把多项式看作函数,画出**图像**： 比如一个一次函数(多项式)$F(x)=2x+1$ 如果一开始不告诉你解析式,只跟你说这个**一次**函数**经过点(0,1)和(1,3)**,你能搞出解析式吗？  **待定系数法**(两点确定一条直线)即可。 我们这里不讲待定系数法,我们只讲待定系数法的**推论** 在平面直角坐标系中,$n+1$ 个点就能唯一确定一个 $n$ 次多项式。 如果对这方面的知识感兴趣,可以查看 [从拉插到快速插值求值](https://www.luogu.com.cn/blog/command-block/zong-la-cha-dao-kuai-su-cha-zhi-qiu-zhi) 所以我们可以**用n+1个点值(有序数对)来描述一个多项式**。 比如两个$n$次多项式$F(x)$和$P(x)$ 令数列$X=\{x_0,x_1,x_2,...,x_n\}$ 把数列$X=\{x_0,x_1,...,x_n\}$代入多项式$F(x)$,得到的$n+1$个点分别为$(x_0,y_0)(x_1,y_1)...(x_n,y_n)$ 同样的,多项式$P(x)$得到的的$n+1$个点分别为$(x_0 , y_0^*)(x_1 , y_1^*)...(x_n , y_n^*)$ 比如说$F(x)=x^2/2+x-1;$(绿色)$\ \ P(x)=-x^2/4+2x$(红色) 取$X=\{-1,0,1\}$(这里次数$n=2$) 函数图像与得到的点如下。  三个**红点**的函数值分别为$Fy=\{-2.25;\ 0;\ 1.75\}$ 三个**绿点**的函数值分别为$Py=\{-1.5;\ -1;\ 0.5\}$(均为从左到右) 让你求两个多项式之积$W(x)=F(x)*P(x)$的**点值表达**,你会怎么干？ 点值表达有一个好处,给你两个多项式点值表达,让你求两个多项式之积的点值表达,直接将对应项乘起来就好了。 (点值的乘法对应着整个多项式的乘法,也就是**浓缩了整个多项式**) 比如上文$F(x)$和$P(x)$这两个多项式的第一个有序数对$(-1,-2.25)$和$(-1,-1.5)$ 就相当于告诉我们$F(-1)=-2.25,G(-1)=-1.5$ 那么$F(-1)*G(-1)$呢？ 废话,肯定等于$-2.25*-1.5$ 所以$W(-1)=F(-1)*G(-1)=-2.25*-1.5=3.375$ $W(0)=F(0)*G(0)=0*-1=0$ $W(1)=F(1)*G(1)=1.75*0.5=0.875$ 所以W在$X$处的函数值为$Wy=\{3.375;\ 0;\ 0.875\}$ **验算**一下？ $W(x)=F(x)*P(x)=(x^2/2+x-1)(-x^2/4+2x)$ $=-x^4/8+x^3-x^3/4+2x^2+x^2/4-2x$ $=-x^4/8+3x^3/4+9x^2/4-2x$ 代入$X=\{-1,0,1\}$ 得到$W(-1)=-1/8-3/4+9/4+2=3.375$ $W(0)=0$ $W(1)=-1/8+3/4+9/4-2=0.875$ 总结一下就是$F(x)G(x)=(F*G)(x)$ 好像那里不对,W的**次数**明摆着是**2n**。 在平面直角坐标系中,给你(n+1)个点就能确定一个n次多项式。 所以**需要2n+1个点**才行,**还差n个**,肿么办？ 反正点都是你自己找的,再多来几个不就好了…… 现在$F(x)$和$P(x)$两个多项式各有$2n+1$个点 所以W的点值表达为$(0, Fy_0*Py_0)(1,Fy_1*Py_1)...(4,Fy_{2n}*Py_{2n})$ 只需要**做2n+1次乘法**即可。 看起来很好用,其实只是几个点而已,和**系数表达**差着十万八千里。 于是我们就想,能不能把系数表达**转换**为点值表达呢？ 算法流程： 把系数表达转换为点值表达->点值表达相乘->把点值表达转换为系数表达。 **这就是 _FFT_ 的算法流程。** “把**系数**表达转换为**点值**表达”的算法叫做**DFT** “把**点值**表达转换为**系数**表达”的算法叫做**IDFT**(DFT的逆运算) P.S: - 从一个多项式的系数表达确定其点值表达的过程称为**求值**(毕竟求点值表达的过程就是取了 n 个 x 然后扔进了多项式求了 n 个值出来); - 而求值运算的逆运算(也就是从一个多项式的点值表达确定其系数表达)被称为**插值**. $\color{Black}\colorbox{lightgreen}{总结一下}$ 在平面直角坐标系中,给你(n+1)个点就能确定一个n次多项式(函数)。 点值表示相乘(点乘)远快于暴力卷积。 -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 2.单位根(复数)及其性质 (如果你不知道$i$的故事,请百度[“虚数$i$”](http://baijiahao.baidu.com/s?id=1581985849864403512&wfr=spider&for=pc)) ### P.S：学会了复数计算之后推荐给自己出点计算题来做 DFT有一个妙处,就是**代入什么由你自己决定**,只要点值个数足够。 我们这些蒟蒻第一感觉都会选择人畜无害的正整数,或者有理数什么的。 但是这些东西虽然在普通人看来计算简单,但并没有什么能够加以利用的优秀性质。 事实证明,找一些毒瘤东西代入进去是个好想法。 上古之时,有一位**dalao**横空出世,他就是**傅里叶**！ 好像是个搞**复数**的专家,有一天他突发奇想： 把**单位根**$\omega_{n+1}^0$~$\omega_{n+1}^n$代入了多项式想求点值表达(~~什么鬼~~) 然后$\omega_{n+1}^0$~$\omega_{n+1}^n$有一些性质(~~蛤玩意儿~~) 然后分治,就$O(nlogn)$了(好像只有最后一个分句看得懂) ------------ 好吧先来介绍复数。 **复数**,即形为$a+bi$的数,这里$\large i^2=-1$ 比如说数轴原来是一条直线：  某个实数一定可以用数轴上一个点来代表。 多么和谐自然。 然而**复数**来了。 数轴原来是一条横着的直线,现在又多了一条竖着的,变成了一个坐标系。  原来的那条直线(横)称为实轴 新的那条直线(竖)称为虚轴 每个复数一定可以用这个坐标系上一个点来代表。 比如蓝点是$-1+2i$,红点是$1+i$ ### 其实形为$a+bi$的复数也可以理解为一个点(a,b) 精妙之处在于,**复数之间能做运算**。  不懂的话就多琢磨下~~我写完这些之后就懂了~~ 我们来介绍下**复数运算**： - 复数相加：实部和虚部分别相加。 如：$(3-2i)+(5+i)=(3+5)+(-2i+i)=8-i$ $(a+bi)+(c+di)=(a+c)+(b+d)i$ - 复数相减：取个相反数再相加。 - 复数相乘：像**一次多项式**一样相乘。 如：$(3-2i)*(5+i)=3*(5+i)-2i*(5+i)$ $\qquad\qquad\qquad\qquad\quad=15+3i-10i-2i^2=17-7i$ **P.S:$i^2=-1$** $(a+bi)*(c+di)=a(c+di)*bi(c+di)=ac+adi+bci+bdi^2$ $\qquad\qquad\qquad\qquad=ac-bd+(ad+bc)i$ - 复数的共轭：实部相同,虚部相反(根据实轴对称)。 如$3-2i$的共轭是$3+2i$ **P.S**:一个复数与自己的共轭相乘一定会得到实数：$(a-bi)(a+bi)=a^2+b^2$ - 复数相除：分子分母同乘分母的共轭 将分母实数化,在分子分母同时乘上分母的共轭,即可把分母变成实数。 然后直觉化简。 $\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{a c + bd+ (bc-ad)i }{c^2+d^2}=\dfrac{a c + bd }{c^2+d^2}+\dfrac{bc-ad}{c^2+d^2}i$ ------------ 而且这些运算在几何(也就是上面的坐标系)里面也有神奇的性质。 ### 复数相乘  这里画的是复数$(3+2i)$与$(1+4i)$相乘的结果$(-5+14i)$ 看出来什么没有？ 好吧,看见那些从$A,B,C$点连向原点的细线了么。 以你做几何题的直觉？？ 没错！如果设原点为点$O$,**数5表示的那个点为P**(忘记画点P)。 有性质 : $∠POC=∠BOA$ , $OB*OC=OA$ 我们把一个复数表示的点到原点的距离叫做这个复数的**模长**,也就是这里的$OB;OC;OA$,复数$x=(a+bi)$的模长称作$|x|$。 也即 : **复数相乘时,模长相乘,幅角相加！** 原点到该点的射线与实轴正方向射线组成的角 **（逆时针旋转）** 的角度乘坐这个复数的**幅角**,复数$a+bi$的幅角称作$arg(a+bi)$。 (不清楚的看图)  接下来是**证明**(最好跟着把图画出来)： 为保证一般性,下面涉及到的都是字母。 - 模长相乘 设两个复数为$a+bi$和$c+di$ 那么它们表示的点就是$(a,b)$和$(c,d)$,分别称这两个点为$B$和$C$ 这两个复数的积为$ac-bd+(ad+bc)i$；表示的点是$(ac-bd,ad+bc)$,称这个点为$A$ - 使用勾股定理可得$BO=\sqrt{a^2+b^2};\ CO=\sqrt{c^2+d^2}$ $AO=\sqrt{(ac-bd)^2+(ad+bc)^2};=\sqrt{a^2c^2 + a^2d^2 + b^2c^2 + b^2d^2}$ $\overline{BO}*\overline{CO}=\sqrt{a^2+b^2}*\sqrt{c^2+d^2}=\sqrt{(c^2+d^2)*(a^2+b^2)}=\sqrt{a^2c^2 + a^2d^2 + b^2c^2 + b^2d^2}=AO$ 证毕。 - 幅角相加 第二个规律需要相似+爆算的,不想了解请跳过。  那么,如果能证明途中两个三角形相似,就能证明幅角相加定理了。 我们已经证明了$OA=OC*OB$(模长相乘),相似比例就是$OB$。 $OB:O1$也是$OB$,我们只要再证明$AB:1C=OB$即可,也即$AB^2=OB^21C^2$。 $AB^2=(ac-bd-c)^2+(ab+bc-d)^2=a^2c^2+a^2d^2-2ac^2-2ad^2+b^2c^2+b^2d^2+c^2+d^2$ $OB^21C^2=(c^2+d^2)((a-1)^2+b^2)=a^2c^2+a^2d^2-2ac^2-2ad^2+b^2c^2+b^2d^2+c^2+d^2=AB^2$ 证毕。这玩意特别重要。 ------------ 接下来我们介绍**单位根**。 看起来很高大上,这是傅里叶快速变换的**重要**内容。 - n次单位根(**n为正整数**)是**n次幂为1**的复数。 换句话说,就是方程$\large x^n=1$的复数解。 到这里你也许会有一点不懂,我们来一起推导一下。 学过三角函数的话应该对**单位圆**很熟悉。 不熟悉的话也没事,单位圆就是： ### 圆心在原点且半径为1的圆(如图)  我们在复平面上放一个单位圆,在这个单位圆上的点表示的复数都有一个性质： **模长为1**(圆上每一点到圆心的距离都相等) 可还记得？n次单位根是**n次幂为1**的复数,1的**模长就是1**。 考虑一个复数$x$。 如果$|x|>1$,即$|x^n|=|x|^n>1$(模长乘模长,越乘越大) 如果$|x|<1$,即$|x^n|=|x|^n<1$(模长乘模长,越乘越小) 所以只有**模长等于1**的复数才有可能成为n次单位根。 也就是说,只有单位圆上的点表示的复数才有可能成为n次单位根(必要性) 我们成功缩小了探究范围,**只在这一个圆上研究**。 (下面提及的复数,均是在单位圆上的复数!!!) 接下来,我们要**找**单位根。 这些单位根的模长都是$1$,只有幅角存在差别,**下面我们就不考虑模长**。 容易找到 : **幅角**为$0,\dfrac{1}{n}\text{圆周},\dfrac{2}{n}\text{圆周},...,\dfrac{n-1}{n}\text{圆周}\ $的复数,都是单位根,共n个。 也就是说,假如一个复数,其幅角为$\dfrac{a}{n}(0\leq a<n,a∈Z)$,那么这就是一个单位根。 容易知道这玩意的$n$次方幅角是$a$倍圆周,那么就和1重合。 根据代数基本定理 : **n次方程在复数域内有且只有n个根**。 所以这些不重不漏,就是所有的单位根了。 还不懂也没关系,我们来实践一下。 设$n=3;$ 有一个模长为1的复数为3次单位根,它的幅角为(**圆周$*a$**) 它的3次方的幅角为$3a*$圆周 那么,如果它是3次单位根,幅角的3倍$(3a*$圆周$)$一定是圆周的自然数倍。 **所以3a是自然数。** a=0时,它的幅角为0(其实这个复数就是1) a=1/3时,它的幅角为(圆周/3) a=2/3时,它的幅角为(2圆周/3)  有**幅角**为$0,\dfrac{1}{n}\text{圆周},\dfrac{2}{n}\text{圆周},...,\dfrac{n-1}{n}\text{圆周}\ $的单位根,共n个。 每一个都比上一个"多转了"$\dfrac{1}{n}\text{圆周}$。 所以n次单位根**n等分**单位圆。(重要)。 怎么称呼这些单位根呢 “$\omega$”是单位根的符号,偶尔也写作$w$(~~颜文字请走开~~,LaTeX是`$\omega$`) 我们**从1开始**(1一定是单位根),沿着单位圆**逆时针**把单位根标上号。 $\omega_{n}^0=1$ $\omega_{n}^1$为第二个n次单位根 ... $\omega_{n}^{n-1}$为第n个n次单位根  比如上图,最右边的点就是$\omega_{3}^0$,左上角的点就是$\omega_{3}^1$,左下角的点就是$\omega_{3}^2$. **P.S** 虽然我们只承认$\omega_{n}^0,\omega_{n}^1,\omega_{n}^2$~$\omega_{n}^{n-1}$, 但是也有$\omega_{n}^{k}$的$k>=n$或$k<0$的情况([带正负的角](https://baike.baidu.com/item/%E8%A7%92/3555592?fr=aladdin#1))。 有$\large\omega_{n}^k=\omega_{n}^{k\%n}$。 随意称呼,切了书写方便,一般不取模。 ------------ 关于单位根的基本定义相信你已经Get到了。单位根的世界,就是一个单位圆。 下面还有一些性质(类比**切蛋糕**记忆) - -1.对任意的n,$\ \ \large{\omega_{n}^{0}=1}$ 显而易见,不懂的把上面的图看一看。 - 0.$\ \large{\omega_{n}^k=(\omega_{n}^1)^k}$ $\omega_{n}^k=\omega_{n}^1*\omega_{n}^1*...*\omega_{n}^1\text{(共k个)}$ 每乘上一个$\omega_{n}^1$代表“**把幅角逆时针转动**$\dfrac{1}{n}\text{圆周}$” 转了k次,即为$\dfrac{k}{n}\text{圆周}$。 - 1.两个n次单位根$\ \large{\omega_{n}^j*\omega_{n}^k=\omega_{n}^{j+k}}$ 由(0)可得$\ \large{\omega_{n}^j*\omega_{n}^k=(\omega_{n}^1)^j*(\omega_{n}^1)^k=(\omega_{n}^1)^{j+k}}=\omega_{n}^{j+k}$ 感性理解 : 我们把一个圆**等分**成n份(想想你切蛋糕的时候) 由(1)可得: - $(ω^k_n)^{-1}=(ω^k_n)^{-1}=ω_n^{-k}=ω_n^{n-k}$ 根据$\omega_{n}^j*\omega_{n}^k=\omega_{n}^{j+k}$ $ω_n^{n-k}*ω^k_n=\omega_{n}^{n-k+k}=\omega_{n}^{n}=\omega_{n}^{0}=1$ - $(ω^k_n)^j=ω_n^{kj}$ 根据$\omega_{n}^j*\omega_{n}^k=\omega_{n}^{j+k}$ $(ω^k_n)^j=\begin{matrix}\text{共j个}\\\overbrace{(ω^k_n)*(ω^k_n)*...*(ω^k_n)}\end{matrix}==ω_n^{kj}$ - 任何一个n次单位根都可以写成$ω_n^1$的整幂 显而易见,$ω_n^k=(ω_n^1)^k$ - 2.$\large{ω^{2k}_{2n}=ω^k_n}$ $ω^k_n$相当于我们把一个圆(蛋糕)**等分**成**n**份然后取**k**份 类似的,$ω^{2k}_{2n}$相当于我们把一个圆(蛋糕)**等分**成**2n**份然后取**2k**份 显然是等价的,只不过多切了几刀而已。 同理,有$ω^{pk}_{pn}=ω^k_n$ - 3.如果n是偶数,$\ \large{ω^{(k+n/2)}_{n}=-ω^k_n}$ 根据$\omega_{n}^{(j+k)} = \omega_{n}^j*\omega_{n}^k$ $ω^{(k+n/2)}_{n}=ω^k_n*ω^{n/2}_{n}$ 根据(1),把“乘$\omega_{n}^{n/2}$”这个操作理解为把幅角逆时针转动($\large\frac{(\frac{n}{2})}{n}=\frac{1}{2}$圆周)。 也就是旋转半个圆周,也就是**关于原点对称**,也就是**取相反数**(建议**画个图**) 好啦,关于单位根想必你已经掌握啦(~~终于！！~~)。 _**坐稳了,前方高能！！！**_ $\color{Black}\colorbox{lightgreen}{总结一下}$ 复数把数轴扩展到了**复平面**上,复数可以对应复平面上一个点。 复数也有**四则运算**。复数相乘时,**模长**相乘,**幅角**相加。 n次单位根(n为正整数)是**n次幂为1**的复数。 把单位根画到单位圆上之后,就能整出一些**性质**。 $\color{Black}\colorbox{yellow}{习题}$ - 画出五次单位根并计算它们的具体值 -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 3.DFT的加速版本 我们来讲讲FFT的**分治方式** 傅里叶一巴掌把这个多项式拍成碎片,然后**按奇偶**拼成两块 考虑**n-1次**(也就是说有n项)多项式F(x) 的每一项按照下标(次数)的奇偶分成两部分： $F(x)=(F[0]+F[2]x^2+...+F[n-2]x^{n-2})+(F[1]x+F[3]x^3+...+F[n-1]x^{n-1})$ 这里**保证n是2的整幂**,不会出现分得不均匀的情况。 又设**两个有n/2项**的多项式$FL(x)$和$FR(x)$, 他们的系数依赖于$F(x)$(具体看式子) $FL(x)=F[0]+F[2]x+...+F[n-2]x^{n/2-1}$ $FR(x)=F[1]+F[3]x+...+F[n-1]x^{n/2-1}$ (又把多项式忘了的,赶紧回去看) 则可以得出$F(x)=FL(x^2)+xFR(x^2)$(**非常重要**) 我们要把$ω^k_n$代入$F(x)$： - $k<n/2$,代入$ω^k_n$ $F(ω^k_n)=FL((ω^k_n)^2)+ω^k_nFR((ω^k_n)^2)$ 因为$(ω^k_n)^2=ω^k_{n/2}$（不理解这个的请见第二节） $\large{F(ω^k_n)=Fl(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})}$ - $k<n/2$,代入$ω^{k+n/2}_n$ $F(ω^{k+n/2}_n)=FL((ω^{k+n/2}_n)^2)+ω^{k+n/2}_nFR((ω^{k+n/2}_n)^2)$ (下一步用到了$(ω^k_n)^j=ω_n^{kj}$) $\qquad\qquad\ \ \ =FL(ω^{2k+n}_n)+ω^{k+n/2}_nFR(ω^{2k+n}_n)$ (下一步用到了$\omega_{n}^k=\omega_{n}^{k\%n}$) $\qquad\qquad\ \ \ =FL(ω^{2k}_n)+ω^{k+n/2}_nFR(ω^{2k}_n)$ (下一步用到了$\omega_{2n}^{2k}=\omega_{n}^{k}$) $\qquad\qquad\ \ \ =FL(ω^{k}_{n/2})+ω^{k+n/2}_nFR(ω^{k}_{n/2})$ (下一步用到了$\omega_{n}^{k+n/2}=-\omega_{n}^k$) $\qquad\qquad\ \ \ =FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})$ 最后得到$\large{F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})}$ **比对一下**上一个式子$F(ω^k_n)=Fl(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})$ 等式右边不就是只有一个**正负**号的区别吗？ 这两条式子有什么用呢？ 到这里,如果我们知道两个多项式$FL(x)$和$FR(x)$分别在$ω^0_{n/2},ω^1_{n/2},ω^2_{n/2},...,ω^{n/2-1}_{n/2}$的点值表示, 套用$F(ω^k_n)=Fl(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})$可以$O(n)$求出$F(x)$在$ω^0_{n},ω^1_{n},ω^2_{n},...,ω^{n/2-1}_{n}$处的点值表示。 套用$F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})$可以$O(n)$求出$F(x)$在$ω^{n/2}_{n},ω^{n/2+1}_{n},ω^{n/2+2}_{n},...,ω^{n-1}_{n}$处的点值表示。 所以如果我们知道两个多项式$FL(x)$和$FR(x)$分别在$ω^0_{n/2},ω^1_{n/2},ω^2_{n/2},...,ω^{n/2-1}_{n/2}$的点值表示,就可以$O(n)$求出$F(x)$在$ω^0_{n},ω^1_{n},ω^2_{n},...,ω^{n-1}_{n}$处的点值表示。 我们就成功的获得了$F(x)$的点值标示。 (懵逼不怕,具体见代码) 好像有哪里不对？ **上文:** 如果我们知道两个多项式$FL(x)$和$FR(x)$分别在$ω^0_{n/2},ω^1_{n/2},ω^2_{n/2},...,ω^{n/2-1}_{n/2}$的点值表示...就可以$O(n)$求出$F(x)$在$ω^0_{n},ω^1_{n},ω^2_{n},...,ω^{n-1}_{n}$处的点值表示。 问题在于我们不知道啊？ 对$FL(x)$和$FR(x)$暴力代入？ NO~NO~NO~ 上面的工作,其实是一个**分治**的过程。 $FL(x)$和$FR(x)$这两个多项式,是规模为原多项式一般的**子问题**,他们的**性质完全相同**！ 这个可以一直**分治**下去,直到多项式只剩下一个项为止(此时什么操作也不用做)。  ### 实践出真知,我们现在就用代码来表现丰富多彩的数学吧！ $\color{Black}\colorbox{lightgreen}{总结一下}$ FFT**分治**的过程,就是根据： $\large{F(ω^k_n)=Fl(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})}$ $\large{F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})}$ 这两个式子,实现子问题的合并。 记忆方法:$\large F(x)=FL(x^2)+xFR(x^2)$,将单位根分上下半圆代入讨论。 $\color{Black}\colorbox{yellow}{习题}$ - 默写上述两个式子的证明过程 -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 4.DFT的代码实现 还有一个小细节。 上文有一句话：“保证**n是2的整幂**,不会出现分得不均匀的情况。” 实际应用中,n不一定是2的正整数次幂。 我们可以补项,在最高次强行添加一些**系数为0**的项(类似于高精度补0)。不影响我们的计算结果,却占了位置。(具体见代码) 讲完了这些我们可以开始写$DFT$了 ## 1.复数结构体 根据上文复数的四则运算重载。 `CP`是 $complex$ 的简称。 ```cpp #include <iostream> using namespace std; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} CP operator / (CP const &B) const { double t=B.x*B.x+B.y*B.y; return CP((x*B.x+y*B.y)/t, (y*B.x-x*B.y)/t); } }a,b; int main() { cin>>a.x>>a.y>>b.x>>b.y; CP c; c=a+b; cout<<'('<<c.x<<','<<c.y<<")\n"; c=a-b; cout<<'('<<c.x<<','<<c.y<<")\n"; c=a*b; cout<<'('<<c.x<<','<<c.y<<")\n"; c=a/b; cout<<'('<<c.x<<','<<c.y<<")\n"; return 0; } ``` (别用STL,会被卡常) ------------ ### 2.预处理单位根 之前我们一直在说“把**单位根**代入...” 现在我们来学习如何求单位根。 前面我们说过$\ \large{\omega_{n}^k=(\omega_{n}^1)^k}$ 也就是说只要求出$ω^1_n$然后把它乘n次,就能得到$\{ω^0_n,ω^1_n,ω^2_n,ω^3_n,...,ω^{n-1}_n\}$ 我们怎么求$ω^1_n$呢？  点A就是$ω^1_n$,我们需要求出它的横坐标$OB$和纵坐标$AB$ 在$RT$三角形$OAB$中,我们已知$∠AOB=\dfrac{1}{n}$圆,$∠OBA=90$°,$OA=1$。 其他什么也不知道。我们需要求$OB$和$AB$。 [百度百科](https://baike.baidu.com/item/%E4%B8%89%E8%A7%92%E5%87%BD%E6%95%B0/1652457?fr=aladdin#2)  | 英文缩写 | 表达式 | 语言描述 | 函数名称 | | :----: | :----: | :----: | :----: | | sin | a/c | ∠A的对边比斜边 | 正弦函数 | | cos | b/c | ∠A的邻边比斜边 | 余弦函数 | | tan | a/b | ∠A的对边比邻边 | 正切函数 | 观察上图,需要求$OB$和$AB$,它们俩分别是$∠O$的邻边和对边. $sin(∠O)=∠O$的对边比斜边; $cos(∠O)=∠O$的邻边比斜边; 所以$OB=cos(∠O)*AO=cos(∠O);\ AB=sin(∠O)*AO=sin(∠O)$ （**模长**$AO=1$） 所以$ω^1_n$的坐标为$(cos(∠O),sin(∠O))$ C++的三角函数采用**弧度制**。一个整圆不再是$360$°,**而是$2π$** 所以$∠O=\dfrac{1}{n}$圆$=\dfrac{2π}{n}$,$ω^1_n$的坐标为$\left(cos(\dfrac{2π}{n}),sin(\dfrac{2π}{n})\right)$。 只要求出$ω^1_n$然后把它乘n次,就能得到$\{ω^0_n,ω^1_n,ω^2_n,ω^3_n,...,ω^{n-1}_n\}$ 代码： ```cpp #include<cstdio> #include<cmath> #define Maxn 1000500 //用这句话能得到得到精确的π,可以当做结论来记 const double Pi=acos(-1); using namespace std; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }w[Maxn]; //w长得是不是很像ω? int n; int main() { scanf("%d",&n); CP sav(cos(2*Pi/n),sin(2*Pi/n)),buf(1,0); for (int i=0;i<n;i++){ w[i]=buf; buf=buf*sav; } for (int i=0;i<n;i++) printf("w[%d][n]=(%.4lf,%.4lf)\n",i,w[i].x,w[i].y); //由于精度问题会出现-0.0000的情况,将就看吧 return 0; } ``` ------------ ### 3.递归实现DFT(简单版) 具体见代码,多说无益。 ```cpp #include <cstdio> #include <cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }f[Maxn<<1],sav[Maxn<<1]; void dft(CP *f,int len) { if (len==1)return ;//边界 //指针的使用比较巧妙 CP *fl=f,*fr=f+len/2; for (int k=0;k<len;k++)sav[k]=f[k]; for (int k=0;k<len/2;k++)//分奇偶打乱 {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} dft(fl,len/2); dft(fr,len/2);//处理子问题 //由于每次使用的单位根次数不同(len次单位根),所以要重新求。 CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); for (int k=0;k<len/2;k++){ //这里buf = (len次单位根的第k个) sav[k]=fl[k]+buf*fr[k];//(1) sav[k+len/2]=fl[k]-buf*fr[k];//(2) //这两条语句具体见Tips的式子 buf=buf*tG;//得到下一个单位根。 }for (int k=0;k<len;k++)f[k]=sav[k]; } int main() { scanf("%d",&n); for (int i=0;i<n;i++)scanf("%lf",&f[i].x); //一开始都是实数,虚部为0 for(m=1;m<n;m<<=1); //把长度补到2的幂,不必担心高次项的系数,因为默认为0 dft(f,m); for(int i=0;i<m;++i) printf("(%.4f,%.4f)\n",f[i].x,f[i].y); return 0; } ``` **Tips** : (1)$F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})$ (2)$F(ω^k_n)=FL(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})$ 现在问题来了,DFT输出的都是些杂乱的**点值表达**,所以解决不了问题。 上文说过,**IDFT**是DFT的逆(~~CP~~),她可以把点值还原成多项式,最终完成乘法。 现在到了讲IDFT的时候了!!! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 5.IDFT理论与FFT初步实现 **IDFT和DFT**就是两句话的区别,这个结论实在太好记了。 **结论** : 把**DFT**中的$ω^1_n$换成$ω^{-1}_n$,做完之后除以$n$即可。(记忆) 当然,证明还是要有的,而且很有价值,从这里可以延伸出 [单位根反演](https://www.luogu.com.cn/blog/command-block/dan-wei-gen-fan-yan-xiao-ji)。 我们知道$DFT$就是求点值,我们不需要再去证明点值的相关性质,我们只要弄出一个$DFT$的逆运算就好了。 还是多项式$F(x)$,设我们变换之后,得到的点值数列为$G$。 即$G=DFT(F)$。 脑补代入的过程,可得$\large G[k]=\sum\limits_{i=0}^{n-1}(ω^k_n)^iF[i]$ 那么**结论**就是$\large n*F[k]=\sum\limits_{i=0}^{n-1}(ω^{-k}_n)^iG[i]$ 通过这个式子我们可以把点值还原。 证明: 右边$=\sum\limits_{i=0}^{n-1}(ω^{-k}_n)^iG[i]$ 代入可得: $=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}ω^{ij}_nω^{-ik}_nF[j]$ $=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}ω^{i(j-k)}_nF[j]$ 分类讨论: - ①$j==k$ 那么贡献就是$\sum\limits_{i=0}^{n-1}ω^{0}_nF[k]=nF[k]$ - ②$j≠k$,设$p=j-k$ 那么贡献就是$\sum\limits_{i=0}^{n-1}ω^{ip}_nF[k+p]$ $=ω^{p}_n(\sum\limits_{i=0}^{n-1}ω^{i}_n)F[k+p]$ 等比数列求和可以得到: $(\sum\limits_{i=0}^{n-1}ω^{i}_n)=\dfrac{ω_n^{np}-1}{ω_n^{p}-1}=\dfrac{ω_n^{0}-1}{ω_n^{p}-1}=0$ 所以这部分贡献为0,我们就成功地证明了上述结论。 > 如果你想知道的多一些,我可以告诉你这个东西本质是单位根反演。 > 或者可以理解成范德蒙德矩阵求逆。 $nF[k]=\sum\limits_{i=0}^{n-1}(ω^{-k}_n)^iG[i]$ 就是相当于把$G$数列**当做系数**,再**代入**一遍求值。 不同的是,这次代入的是$\{ω^{0}_n,ω^{-1}_n,ω^{-2}_n,...,ω^{-n+1}_n\}$ 这里就要求找到一个东西能求出$\{ω^{0}_n,ω^{-1}_n,ω^{-2}_n,...,ω^{-n+1}_n\}$ $ω^{-1}_n$就可以。 我悄悄告诉你$ω^{-1}_n$长这个样子  得$(ω^{-1}_n)^j=ω_n^{-j}$ 也就是说只要求出$ω^{-1}_n$然后**把它乘n次**,就能得到$\{ω^{0}_n,ω^{-1}_n,ω^{-2}_n,...,ω^{-n+1}_n\}$ 它的各方面和$ω^{1}_n$都很像,其实就是$\left(cos(\dfrac{2π}{n}),-sin(\dfrac{2π}{n})\right)$ [题目Link](https://www.luogu.org/problemnew/show/P3803) 至此我们已经写出了第一个版本的FFT(蛤？还有几个？) 代码很好写,和DFT不同的地方很少： ```cpp #include <cstdio> #include <cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }f[Maxn<<1],p[Maxn<<1],sav[Maxn<<1]; void dft(CP *f,int len) { if (len==1)return ;//边界 //指针的使用比较巧妙 CP *fl=f,*fr=f+len/2; for (int k=0;k<len;k++)sav[k]=f[k]; for (int k=0;k<len/2;k++)//分奇偶打乱 {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} dft(fl,len/2); dft(fr,len/2);//处理子问题 //由于每次使用的单位根次数不同(len次单位根),所以要重新求。 CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); for (int k=0;k<len/2;k++){ //这里buf = (len次单位根的第k个) sav[k]=fl[k]+buf*fr[k];//(1) sav[k+len/2]=fl[k]-buf*fr[k];//(2) //这两条语句具体见Tips的式子 buf=buf*tG;//得到下一个单位根。 }for (int k=0;k<len;k++)f[k]=sav[k]; } void idft(CP *f,int len) { if (len==1)return ;//边界 //指针的使用比较巧妙 CP *fl=f,*fr=f+len/2; for (int k=0;k<len;k++)sav[k]=f[k]; for (int k=0;k<len/2;k++)//分奇偶打乱 {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} idft(fl,len/2); idft(fr,len/2);//处理子问题 CP tG(cos(2*Pi/len), - sin(2*Pi/len)),buf(1,0); //注意这 ↑ 个负号! for (int k=0;k<len/2;k++){ //这里buf = (len次反单位根的第k个) sav[k]=fl[k]+buf*fr[k]; sav[k+len/2]=fl[k]-buf*fr[k]; buf=buf*tG;//得到下一个反单位根。 }for (int k=0;k<len;k++)f[k]=sav[k]; } int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)scanf("%lf",&f[i].x); for (int i=0;i<=m;i++)scanf("%lf",&p[i].x); for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂 dft(f,n);dft(p,n);//DFT for(int i=0;i<n;++i)f[i]=f[i]*p[i];//点值直接乘 idft(f,n);//IDFT for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49)); //注意结果除以n return 0; } ``` [提交记录](https://www.luogu.org/record/25868084) 可以看到因为常数过大而TLE了,看来我们需要一个常数更小的写法。 先别着急卡常,前文我们说过,IDFT和DFT只有一个符号的区别。 那么我们何不减少一下代码量呢: ```cpp #include <cstdio> #include <cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }f[Maxn<<1],p[Maxn<<1],sav[Maxn<<1]; //flag=1 -> DFT flag=0 -> IDFT void fft(CP *f,int len,bool flag) { if (len==1)return ; CP *fl=f,*fr=f+len/2; for (int k=0;k<len;k++)sav[k]=f[k]; for (int k=0;k<len/2;k++) {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} fft(fl,len/2,flag); fft(fr,len/2,flag); CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); if (!flag)tG.y*=-1; for (int k=0;k<len/2;k++){ sav[k]=fl[k]+buf*fr[k]; sav[k+len/2]=fl[k]-buf*fr[k]; buf=buf*tG; }for (int k=0;k<len;k++)f[k]=sav[k]; } int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)scanf("%lf",&f[i].x); for (int i=0;i<=m;i++)scanf("%lf",&p[i].x); for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂 fft(f,n,1);fft(p,n,1);//DFT for(int i=0;i<n;++i)f[i]=f[i]*p[i]; fft(f,n,0); for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49)); return 0; } ``` -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- ------------ # 6.FFT的精细实现 卡常的方法分为以下几类： 0.硬件优化：浮点加速器,CPU等等。 1.编译优化：编译开关,inline,register 之类的。 2.封装优化：把代码的封装拆开,封装得过死会使效率降低(stl)。 3.编写优化：细节优化,手写栈之类。 4.时空互易: 空间换时间以达到时空均衡(常数层面上)。 硬件优化没什么指望,编译优化能力有限又太玄学。 咱们这份代码如果硬要把封装拆开就会显得十分丑陋。 只剩下**编写优化**,**时空互易**两条路子了。 首先观察这个程序片段: ```cpp for (int k=0;k<len/2;k++){ sav[k]=fl[k]+buf*fr[k]; sav[k+len/2]=fl[k]-buf*fr[k]; buf=buf*tG; } ``` 复数的乘法代价是很高的,而`buf*fr[k]`我们计算了两次! 换成下面的操作,常数就小多了: ```cpp for (int k=0;k<len/2;k++){ CP tt=buf*fr[k]; sav[k]=fl[k]+tt; sav[k+len/2]=fl[k]-tt; buf=buf*tG; } ``` 代码中我们使用了大量数组拷贝,这会影响性能。我们能否通过某些手段规避数组拷贝呢? 这个分治究竟在干什么? 分治过程 : 先按照某种顺序把东西分开排列,然后再合并上去。 - “**分开**”的优化 具体的来说： ``` 原来的递归版(数组下标,先偶后奇,从0开始)： 0 1 2 3 4 5 6 7 第1层 0 2 4 6|1 3 5 7 第2层 0 4|2 6|1 5|3 7 第3层 0|4|2|6|1|5|3|7 第4层 ``` 我们要想办法求出第4层的数组状况,然后才能往上合并。 （这个东西叫做蝴蝶变换来着） 具体怎么求呢?多多观察可得,最后的序列是原序列的**二进制反转**。 比如$6=(110)_2$反过来就是$(011)_2=3$ “6”确实排在原来“3”的位置。 如何得到二进制翻转后的数列呢？ $O(nlogn)$倒是可以,但是这就违背了我们卡常的本愿啦…… 可以用递推$O(n)$搞定。 ``` for(int i=0;i<limit;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?limit>>1:0); //tr[i]是i的二进制翻转 ``` 这玩意是怎么得到的呢？ 我们可以把一个二进制数的反转拆成两部分来看： ``` ~~~~~ x 其他部分 最后一位 翻转的话就相当于 x 连接上 其他部分的反转 其他部分的反转=r[i>>1]>>1 然后判一下最后一位,如果是1则加上Limit>>1 ``` 可以类比数位DP理解。 - “**合并**”的优化 观察合并的代码片段。 ```cpp CP *fl=f,*fr=f+len/2; for (int k=0;k<len/2;k++){ CP tt=buf*fr[k]; sav[k]=fl[k]+tt; sav[k+len/2]=fl[k]-tt; buf=buf*tG; }for (int k=0;k<len;k++)f[k]=sav[k]; ``` 注意到指针`fl[k] <--> f[k] . . . fr[k] <--> f[k+len/2]` 那么我们完全可以用如下代码来代替: ```cpp fft(f,len/2,flag); fft(f+len/2,len/2,flag); for (int k=0;k<len/2;k++){ CP tt=buf*f[k+len/2]; f[k+len/2]=f[k]-tt; f[k]=f[k]+tt;//注意赋值顺序! buf=buf*tG; } ``` 这就规避了所有的数组拷贝。 可以写出如下代码: ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} }f[Maxn<<1],p[Maxn<<1]; //flag=1 -> DFT flag=0 -> IDFT void fft(CP *f,int len,bool flag) { if (len==1)return ; fft(f,len/2,flag); fft(f+len/2,len/2,flag); CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); if (!flag)tG.y*=-1; for (int k=0;k<len/2;k++){ CP tt=buf*f[k+len/2]; f[k+len/2]=f[k]-tt; f[k]=f[k]+tt;//注意赋值顺序! buf=buf*tG; } } int tr[Maxn<<1]; int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)scanf("%lf",&f[i].x); for (int i=0;i<=m;i++)scanf("%lf",&p[i].x); for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂 for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0); for (int i=0;i<n;i++) if (i<tr[i])swap(f[i],f[tr[i]]); for (int i=0;i<n;i++) if (i<tr[i])swap(p[i],p[tr[i]]); fft(f,n,1);fft(p,n,1);//DFT for(int i=0;i<n;++i)f[i]=f[i]*p[i]; for (int i=0;i<n;i++) if (i<tr[i])swap(f[i],f[tr[i]]); fft(f,n,0); for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49)); return 0; } ``` [提交记录](https://www.luogu.org/record/25870016) 可以看到还是TLE一个点QAQ 三角函数求值很慢,我们却使用了$O(nlogn)$次三角函数求值,可以考虑预处理来减小常数。 不过一种更有效的方式是:迭代实现(**嘴巴讲不清楚的,看代码秒懂**)。 还有一点,`/2`可以用位运算优化,我们把`/2`统一换成`>>1`。 迭代版代码: ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} }f[Maxn<<1],p[Maxn<<1]; int tr[Maxn<<1]; void fft(CP *f,bool flag) { for (int i=0;i<n;i++) if (i<tr[i])swap(f[i],f[tr[i]]); //枚举区间长度 for(int p=2;p<=n;p<<=1){ int len=p>>1;//待合并的长度 CP tG(cos(2*Pi/p),sin(2*Pi/p)); if(!flag)tG.y*=-1;//p次单位根 for(int k=0;k<n;k+=p){//枚举起始点 CP buf(1,0);//遍历一个子问题并合并 for(int l=k;l<k+len;l++){ CP tt=buf*f[len+l]; f[len+l]=f[l]-tt; f[l]=f[l]+tt; buf=buf*tG; } } } } int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)scanf("%lf",&f[i].x); for (int i=0;i<=m;i++)scanf("%lf",&p[i].x); for(m+=n,n=1;n<=m;n<<=1); for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0); fft(f,1);fft(p,1);//DFT for(int i=0;i<n;++i)f[i]=f[i]*p[i]; fft(f,0); for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49)); return 0; } ``` 怎么觉得比递归版好写？ 这个版本是最经典的,建议背诵全文。 三句话记忆: ```cpp //F(x)=FL(x^2)+x*FR(x^2) //F(W^k)=FL(w^k)+W^k*FR(w^k) //F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k) ``` [评测记录](https://www.luogu.org/record/25876484) 可以看到已经成功地通过了这道题。 ## 附 : "三次变两次"优化 根据$(a+bi)*(c+di)=a(c+di)*bi(c+di)=ac-bd+adi+bci$ 假设我们需要求$F(x)*G(x)$ 设**复多项式**$P(x)=F(x)+G(x)i$ 也就是实部为$F(x)$,虚部为$G(x)$. 则$P(x)^2=(F(x)+G(x)i)^2=F(x)^2-G(x)^2+2F(x)G(x)i$ 发现$P(x)^2$的虚部为 $2F(x)G(x)i$ 也就是说求出$P(x)^2$之后,把它的**虚部**除以2即可. **Code:** 为了卡常数还写了个快读…… ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define Maxn 1350000 using namespace std; const double Pi=acos(-1); inline int read() { register char ch=0; while(ch<48||ch>57)ch=getchar(); return ch-'0'; } int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} }f[Maxn<<1];//只用了一个复数数组 int tr[Maxn<<1]; void fft(CP *f,bool flag) { for (int i=0;i<n;i++) if (i<tr[i])swap(f[i],f[tr[i]]); for(int p=2;p<=n;p<<=1){ int len=p>>1; CP tG(cos(2*Pi/p),sin(2*Pi/p)); if(!flag)tG.y*=-1; for(int k=0;k<n;k+=p){ CP buf(1,0); for(int l=k;l<k+len;l++){ CP tt=buf*f[len+l]; f[len+l]=f[l]-tt; f[l]=f[l]+tt; buf=buf*tG; } } } } int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)f[i].x=read(); for (int i=0;i<=m;i++)f[i].y=read(); for(m+=n,n=1;n<=m;n<<=1); for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0); fft(f,1); for(int i=0;i<n;++i)f[i]=f[i]*f[i]; fft(f,0); for(int i=0;i<=m;++i) printf("%d ",(int)(f[i].y/n/2+0.49)); return 0; } ``` 相信经过上面的学习,这份代码大家都能够理解,我就不再赘述了. [评测记录](https://www.luogu.org/record/25877049) $\quad$( $2020.5.12$ 更新 : [评测记录](https://www.luogu.com.cn/record/33544825)) 快了一点吧！虽然本蒟蒻的常数还是很大QwQ - **Warning** : 如果需要卷积的两个多项式值域相差太大,就会卡精度,慎用。 比如说$F(x)$的系数均在$[10^{-5},10^{-6}]$之间,$G(x)$的系数均在$[10^5,10^6]$之间。 直接做$FFT$,涉及的精度跨度上限是${10^{12}}$。 假如使用三次变两次优化,由于平方项的存在,涉及的精度跨度上限是${10^{24}}$,严重掉精度。 下方例题第二题直接铁头硬上会WA。 至于修正方法,可以把两个多项式**数乘**到相同的值域范围,这样子就不会产生平方项的大误差了。 注意最后要除回去。 ## 附 : 点值的性质 我们的点值代表了整个多项式,加法和乘法都是对应的。 这可以用于某些(初等)卡常,比如我们如果要计算$A*B+B*C+A*C$ 如果完全封装,每次乘法就要计算$3$次FFT,一共就是$9$次FFT。 化为$A*(B+C)+B*C$就是两次乘法,需要$6$次FFT。 我们也不必对$A*(B+C)$,$B*C$分别IDFT,可以加在一起统一做,这样一共就是$5$次FFT。 又注意到$DFT(B+C)=DFT(B)+DFT(C)$,又可以省去一次。 采用"三次变两次"优化,构造$H1=(A+(B+C)i);H2=(B+Ci)$ 平方之后也不需要分别IDFT,这样做是$3$次FFT,只折合原来的一次暴力。 $\color{Black}\colorbox{lightgreen}{总结一下}$ 卡常有什么好总结的？去~~吊打集训队~~问[**wys**](https://www.luogu.org/problemnew/show/P4604)吧。 $\color{Black}\colorbox{yellow}{习题}$ - [A*B Problem升级版](https://www.luogu.org/problem/P1919) - [力](https://www.luogu.org/problemnew/show/P3338) ------------ # 5.后话 那么讲到这里,想必FFT的大部分内容你已经掌握了,到了说再见的时候了…… **再见啦~(挥手)** 欲知后事如何,请看下集:[NTT与多项式全家桶](https://www.luogu.org/blog/command-block/ntt-yu-duo-xiang-shi-quan-jia-tong)! (看完你会了解更多的多项式工业)