Typesetting math: 76%

快速傅里叶变换(FFT)

更好更精简更新的阅读体验(点我)

参考资料:快速傅里叶变换(FFT)详解 - 自为风月马前卒- 博客园

0. 前言

写这篇文章的时候

找不到小学数学笔记本了……

找不到初中数学笔记本了……

找不到高中数学笔记本了……

找不到脑子了……

找不到c++编译环境了……

找不到破解的几何画板了……

找不到Markdown了……

找不到……了

所以推荐在看文章之前推荐把以上东西提前准备好~

本文所有的nn都是22的整数次幂……

1. 前置知识

我再扔一次目录。觉得自己会的可以跳过……

开始进入正题……

加法

不想讲……

乘法

乘法运算的结果叫做积

很多情况下a×b=ab=aba×b=ab=ab

a×b=a++abaa×b=a++aba

同时在乘法这里讲一下累加

ri=lai=ri=lai=al+al+1++ar1+arri=lai=ri=lai=al+al+1++ar1+ar

也可以讲一下累乘

ri=lai=ri=lai=al×al+1××ar1×arri=lai=ri=lai=al×al+1××ar1×ar

乘方

乘方的结果叫做幂(power),

ab=a××abaab=a××aba

这里aa叫做底数,bb叫做指数。

根据上面的式子,可以推出

ab×ac=a××aba×a××aca=a××a(b+c)a=ab+cab×ac=a××aba×a××aca=a××a(b+c)a=ab+c

以及同理可得ab÷ac=abcab÷ac=abc

所以我们也可以知道a0=1|aRC,a0a0=1|aRC,a0

至于集合什么的……我们留到后面

为什么要强调a0a0

因为0000没有意义

我们可以知道一个00相乘是00, 就是01=001=0,然后00=01÷01=0÷000=01÷01=0÷0

小学老师告诉过我们,00除以任何非零数都是0000不能作除数,0÷00÷0没有意义

我们来设一下a÷b=ca÷b=c,即a=b×ca=b×c

如果b=0,a0b=0,a0,解关于cc方程告诉我们无解……

如果a=0,b=0a=0,b=0,解关于cc的方程告诉我们有无数个解……

所以没有意义。

所以整篇文章我们都默认底数不为00

开方

开方是乘方的逆运算。

aannnn为自然数)次方根指的是n方幂等于a的数,也就是b|bn=a,nNb|bn=a,nN。(转自百度百科)

n=2n=2时称bbaa的平方根,n=3n=3时成为立方根。所有的次方根统称方根

算术平方根

aa的算术平方根是平方为aa的非负数,即b|b2=a,b0b|b2=a,b0,也可以写成b=2a=ab=2a=a

知道这个之后我们就可以发现

ab×ab=ab+b=a2ba2b=abab×ab=ab+b=a2ba2b=ab

所以我们也可以推出mab=abmmab=abm

但是当bb不能整除mm的时候……

其实也是可以计算的

举个栗子

41.5=432=243=264=841.5=432=243=264=8

栗子吃完了,你也该懂了……

现在来说说负数次幂

其实很简单啦……

bc=mbc=m,则有ab÷ac=amab÷ac=am

m<0m<0,即b<cb<c,

am=ab÷ac=(a××aba)÷(a××aca)=1÷(a××acba)=1acbam=ab÷ac=(a××aba)÷(a××aca)=1÷(a××acba)=1acb

所以知道了吧,ab=1abab=1ab,这样转换过来就可以计算出负次数幂。

集合

我感觉我在讲数学……

转自百度百科

概念

集合是指具有某种特定性质的具体的或抽象的对象汇总而成的集体。其中,构成集合的这些对象则称为该集合的元素

表示

假设有实数x<yx<y

  1. [xy][xy] :方括号表示包括边界,即表示大于等于xx、小于等于yy的数。

  2. (xy)(xy):小括号是不包括边界,即表示大于xx、小于yy的数。

符号们

写在后面了……(因为太长了)给你们一个目录超链接去看吧……

特性

确定性

给定一个集合,任给一个元素,该元素或者属于或者不属于该集合,二者必居其一,不允许有模棱两可的情况出现。

互异性

一个集合中,任何两个元素都认为是不相同的,即每个元素只能出现一次。有时需要对同一元素出现多次的情形进行刻画,可以使用多重集,其中的元素允许出现多次。

无序性

一个集合中,每个元素的地位都是相同的,元素之间是无序的。集合上可以定义序关系,定义了序关系后,元素之间就可以按照序关系排序。但就集合本身的特性而言,元素之间没有必然的序。

分类

空集

不包含任何元素,如{x|xR,x2+1=0}{x|xR,x2+1=0},在实数范围内方程无解,所以该集合是空集,用表示。他有两个特点

  1. 空集是任意一个非空集合的真子集。
  2. 空集是任意一个集合的子集
子集

S,TS,T是两个集合,如果SS的所有元素都属于TT, 即xSxTxSxT, 则称SSTT的子集,记为STST

显然,对于任何集合SS,都有SSSS, SS,其中,符号读作包含于,表示该符号左边的集合中的元素全部是该符号右边集合的元素。

如果STST, 但TT中存在一个元素xx, xT,xSxT,xS, 即STST, 则称SSTT的一个真子集,即 STST或者STST

交集

由属于A且属于B的相同元素组成的集合,记作ABABBABA,读作“A交B”(或“B交A”),即AB={x|xA,xB}AB={x|xA,xB}.注意交集越交越少。若A包含B,则AB=BAB=AAB=BAB=A

并集

由所有属于集合A或属于集合B的元素所组成的集合,记作ABAB(或BABA),读作“A并B”(或“B并A”),即AB={x|xA,xB}AB={x|xA,xB}.注意并集越并越多,这与交集的情况正相反

相对补集

由属于A而不属于B的元素组成的集合,称为B关于A的相对补集,记作ABABA\BA\B,即AB={x|xAxB}AB={x|xAxB}

绝对补集

AA关于全集合UU的相对补集叫做AA的绝对补集,记作AA~A~A

幂集(好像不要求掌握)

设有集合A,由集合AA所有子集组成的集合,称为集合AA幂集。对于幂集有定理如下:有限集AA的幂集的基数等于22的有限集AA的基数次幂

区间

a,b(a<b)a,b(a<b)是两个相异的实数,则满足不等式a<x<ba<x<b的所有实数xx的集合称为以a,ba,b为端点的开区间,记为(a,b)={x:a<x<b}(a,b)={x:a<x<b}

满足不等式axbaxb的所有实数xx的集合成为以a,ba,b为端点的闭区间,记为[a,b]={x:axb}[a,b]={x:axb}

满足不等式a<xba<xb或不等式ax<bax<b的所有实数xx的集合称为以a,ba,b为断电的半开半闭区间,分别记为(a,b]={x:a<xb}(a,b]={x:a<xb}[a,b)={x:ax<b}[a,b)={x:ax<b}

除此之外,还有下述几类无限区间:(因为无限取不到,所以无限那半边肯定是开区间)

(a,+)={x:x>a}(,b)={x:x<b}[a,+)={x:xa}(,b]={x:xb}(,+)=R(a,+)={x:x>a}(,b)={x:x<b}[a,+)={x:xa}(,b]={x:xb}(,+)=R

相等集合

如果两个集合SSTT的元素完全相同,则称SSTT两个集合相等,记为S=TS=T 。显然有如下关系:

S=TSTTSS=TSTTS

表示方法

列举法

举法就是将集合的元素逐一列举出来的方式。例如,光学中的三原色可以用集合{红,绿,蓝}表示;由四个字母a,b,c,da,b,c,d组成的集合AA可用A=a,b,c,dA=a,b,c,d表示,如此等等。

列举法还包括尽管集合的元素无法一一列举,但可以将它们的变化规律表示出来的情况。如正整数集N+N+和整数集ZZ 分别可以表示为N+={1,2,3,n,}N+={1,2,3,n,}Z={0,±1,±2,,±n,}Z={0,±1,±2,,±n,}

描述法

描述法的形式为{代表元素|满足的性质}。

设集合SS是由具有某种性质PP的元素全体所构成的,则可以采用描述集合中元素公共属性的方法来表示集合:S=x|P(x)S=x|P(x)。例如,由2的平方根组成的集合B可表示为B=x|x2=2B=x|x2=2

而有理数集QQ和正实数集R+R+可以表示为Q={x|x=qp,pN+,qZ}Q={x|x=qp,pN+,qZ}R+={x:xRx>0}R+={x:xRx>0}

图像法

略……

符号法

有些集合可以用特殊符号表示,例如:

NN:非负整数集合或自然数集合0,1,2,3,0,1,2,3,

NNN+N+:正整数集合1,2,3,1,2,3,

ZZ整数集合,1,0,1,,1,0,1,

QQ有理数集合

Q+Q+:正有理数集合

QQ:负有理数集合

RR实数集合(包括有理数和无理数)

R+R+:正实数集合

RR:负实数集合

II:虚数集合

CC复数集合

:空集(不含有任何元素的集合)

(我集合好像没有讲很久吧……)

多项式

系数表示法

就是……关于xx的一个多次的式子加起来。

f(x)f(x)表示一个关于xxn1n1次多项式,那么就会有

f(x)=n1i=0ai×xif(x)=n1i=0ai×xi

我们可以发现对于固定的xx, x0,x1,,xnx0,x1,,xn都是固定的

所以表示一个多项式我们也可以去掉xx,直接表示成f(x)=a0,a1,a2an1f(x)=a0,a1,a2an1,这就是我们平时见到的多项式(?)

点值表示法

不知道什么时候我们看了一眼函数

比如说二次函数y=ax2+bx1+cx0y=ax2+bx1+cx0,对比一下xx的二次多项式f(x)=a2x2+a1x1+a0x0f(x)=a2x2+a1x1+a0x0

是不是很像?

如果对于每一个f(x)f(x)我们都把a0,a1,,an1a0,a1,,an1都看做常数的话,我们再在平面直角坐标系上画出y=f(x)y=f(x),可以发现这就是关于xx的多次函数

只要我们确定了ai(0in1)ai(0in1)nn个点,我们就确定了一个关于xxn1n1次函数

简单介绍一下DFT和IDFT

顺便说一下一些容易混淆的概念:

DFT:离散傅里叶变换,O(n2)O(n2)计算多项式乘法

FFT:快速傅里叶变换,O(nlog2n)O(nlog2n)计算多项式乘法

FNTT/NTT:快速傅里叶变换的优化版,优化常数及误差

FWT:快速沃尔什变换,利用类似FFT的东西解决一类卷积问题

MTT:毛爷爷的FFT,非常nb,任意模数

FMT:快速莫比乌斯变化,干嘛的原文没写……

究竟DFTFFT是干嘛的呢

我们通过一道裸题了解一下

题目描述

给定一个n次多项式F(x),和一个m次多项式G(x)。
请求出F(x)和G(x)的卷积。

输入

第一行2个正整数n,m。
接下来一行n+1个数字,从低到高表示F(x)的系数。
接下来一行m+1个数字,从低到高表示G(x)的系数。

输出
一行n+m+1个数字,从低到高表示F(x)∗G(x)的系数。对998244353取模。

样例输入
1 2
1 2
1 2 1
样例输出
1 4 5 2

一道很裸的模拟高精度乘法(?)(反正思路差不多,都要一个个乘)

那么时间复杂度是多少呢?

保守估计O(n2)O(n2)左右……

好了我没把提示放出来给你们看……

提示
样例解释:(2x+1)*(x^2+2x+1)=2x^3+5x^2+4x+1
故输出 1 4 5 2

两组数据,n,m<=1000000,0<=系数<=9
本题可能需要一定程度上的卡常

你还敢使用那个平方级别算法吗?

再告诉你一句话,这是个模板题

这个时候就需要我们学习一个令人窒息的操作:FFT

不过我们还是要先介绍一下傅里叶变换

前面已经说过,如果确定了ai(0in1)ai(0in1),那么我们就可以确定n1n1次函数.

所以对于两个关于xxnn次多项式,我们可以整俩nn次函数。

先把系数转为点值,然后把两个乘起来。

最后转回去。

这样每次乘法就是O(1)O(1),总共就可以O(n)O(n)暴虐多项式乘法。

是不是很棒?

算了吧……

正常的系数转点值叫做DFT(傅里叶变换)

正常的点值转系数叫做IDFT(离散傅里叶变换)

这两个东西都是O(n2)O(n2)级别的

所以使用傅里叶变换,你就把O(n2)O(n2)的东西变成了O(n2+n)O(n2+n)

那我要你有何用

那咋整。

总不能说了这么多东西白说吧……

我们可以问一下傅里叶

我们可以用不正常的转换方法

诶别用删除线。

就是用不正常的转换方法

我们知道根据点值,我们每次会有点(x,y)(x,y),只靠实数是无法优化的了。

但是我们还有虚数

还有复数。

随便用。

但在此之前还是要把我们的前置知识铺垫好。

三角函数

在讲三角函数值钱我们应该先讲一下角……

初中老师告诉过我们:

具有公共端点的两条射线组成的图形叫做角。这个公共端点叫做角的顶点,这两条射线叫做角的两条边。(静态定义)

但是初中老师好像不会教你傅里叶变换

所以我们翻到高中的定义:

一条射线绕着它的端点从一个位置旋转到另一个位置所形成的图形叫做角。所旋转射线的端点叫做角的顶点,开始位置的射线叫做角的始边,终止位置的射线叫做角的终边。(动态定义)

这里射线OAOA绕点OO旋转至射线OBOB

当然初中老师也不会跟你讲负角,不会跟你讲大于180°的角

现在我来告诉你

如果旋转方向是逆时针,那么这个角叫做正角,大小为ΘΘ,符号为++,合起来就是AOB=+ΘAOB=+Θ

(一般正号省略)

如果旋转方向是顺时针,那么这个角叫做负角,大小为ΘΘ,符号为,合起来就是AOB=ΘAOB=Θ

注意:这里的符号仅表示方向,不表示大小

如果比较大小的话,20°>5°20°>5°

舒服不?

还有一种角叫做零角,就是不旋转。

象限角

有平面直角坐标轴和一个角,角的始边与坐标轴xx的非负半轴重合,顶点和原点重合,终边在哪个象限这个角就是什么象限角。

特别地,当终边在坐标轴上时,这个角不属于任何一个象限。

别问我为什么讲角要讲圆

根据定义,圆就是在一个平面内,一个动点以一个定点为中心,以一定长度为距离旋转一周形成的封闭图形。

圆

别问我圆怎么追踪出来的

OO为定点,OAOA为定长,AA的运动轨迹形成的封闭图形就是OO

圆周角与圆心角

这圆绝对不是追踪出来的

圆心角是指在中心为O的圆中,过弧ABAB两端的半径构成的AOBAOB, 称为弧ABAB所对的圆心角。

顶点在圆上,并且两边都和圆相交的角叫做圆周角,这一定义实质上反映的是圆周角所具备的两个特征:

  1. 顶点在圆上。
  2. 两边都和圆相交。这两个条件缺一不可。

同样也有一个简单的性质,同弧或等弧所对的圆心角相等,圆心角等于同一弧所对的圆周角的两倍。百度上有非常简明易懂的证明……这里篇幅估计已经超了……

角度的弧度制

有没有在一些书上听说过ππ不是圆周率的?

这个时候就是弧度制的问题了

首先弧度制,我们就要把角放在一个圆里面。

让这个角的定点和圆心重合,对,就是让他变成圆心角。

弧度制是用弧长和半径之比来度量对应的圆心角角度的方式。用符号radrad表示,读作弧度,一般忽略。

所以以后看见A=1A=1不是一度的意思,是一弧度的意思。因为它没带单位

那么一弧度是多少度呢?

都说了是用弧长和半径的比来度量了……

那么当弧长等于半径的时候,这个弧所对的圆心角就是1弧度

百度是这样说的:

等于半径长的圆弧所对的圆心角叫做1弧度的角。

由于圆弧长短与圆半径之比,不因为圆的大小而改变,所以弧度数也是一个与圆的半径无关的量。角度以弧度给出时,通常不写弧度单位。

然后我们再推算一下

设在OO中,AOB=1rad=x°AOB=1rad=x°,半径为rr,根据弧长等于半径可以得到AB=r=2πr×x360AB=r=2πr×x360,可以得到等式r=πrx180r=πrx180

化简之后得到x=180πx=180π

也及时1(rad)=180°π1(rad)=180°π

那么180°=π(rad)180°=π(rad)就是移个项的事而已……

为什么我们要用弧度制呢?

为了方便计算。

弧度制下的角的大小,是一个实数,可以和其他数字比较大小。而角度制下只是一个角度。

还有一个很重要的原因:作为程序猿,我相信你不知道(我才知道)(看好了):

当然,这里计算器没有错,C++也没有错

只是他们的单位不一样而已

计算器计算的是sin30°sin30°,C++ 计算的是sin30(rad)sin30(rad)

如果我们修改一下,在C++里面输入sin(3.1415926)sin(3.1415926),或者更精准一点,输入sin(π)sin(π)

这里拓展一个小知识,acos()是cos()的反函数,即acos(cos(x))=xacos(cos(x))=x

后面我们会告诉你cos(π(rad))=1cos(π(rad))=1,这里你只需要知道acos(1)=πacos(1)=π就可以了

舒服了吗?

C++里面是弧度制计算的。

所以我们要学习弧度制。

函数

普通的函数概念就不在这里讲了……

我们只来补充一下函数三要素:定义域值域对应法则

看一下对于函数y=f(x)y=f(x)

定义域

指的就是自变量xx的取值范围

如果f(x)=2xf(x)=2x,在实数范围内,这个函数的定义域(xx的取值范围)就是区间[0,+)[0,+)

如果f(x)=xf(x)=x,在实数范围内,这个函数的定义域就是区间(,+)(,+)

值域

函数经典定义中,随自变量改变的因变量改变的取值范围叫做这个函数的值域,在函数现代定义中是指定义域中所有元素在某个对应法则下对应的所有的象所组成的集合。如:f(x)=x,那么f(x)的取值范围就是函数f(x)的值域。

百度百科的解释是不是太长了?

没关系我们来举个栗子。

如果f(x)=x2f(x)=x2,很明显这个函数的定义域是所有实数,那么根据平方非负,我们可以知道这个函数的值域就是区间[0,+)[0,+)

如果f(x)=xf(x)=x,这个函数的定义域是(,+)(,+),值域是(,+)(,+)

对应法则

对于函数y=f(x)y=f(x)来说,ff指的就是一个对应法则。

若函数的定义域为AA,值域为BB,

y=f(x)y=f(x)表明,xA,!yB使得y=f(x)xA,!yB使y=f(x)

即对于定义域中所有xx的值,在对应法则ff的作用下,可以得到值域中唯一的yy的值

反函数

设有函数y=f(x)y=f(x),定义域是DD, 值域是f(D)f(D)

如果有一个对应法则gg,使得yf(D),!xD使得x=g(y)yf(D),!xD使x=g(y)

那么按照法则gg得到了一个定义域是f(D)f(D),值域是DD的函数,我们把gg成为ff的反函数

记为g=f1g=f1, x=f1(y),yf(D)x=f1(y),yf(D)

由此易得,gg的反函数是ff,也就是fff1f1互为反函数。

同样我们也很容易得到f(g(x))=xf(g(x))=x

而且y=f(x)y=f(x)y=f1(x)y=f1(x)关于直线y=xy=x对称

证明:设点(a,b)(a,b)y=f(x)y=f(x)上的一点,则b=f(a)b=f(a)

然后我们用反函数代进去,y=f1(x)y=f1(x), f1(b)=af1(b)=a,所以点(b,a)(b,a)在函数y=f1(x)y=f1(x)

易得(感性理解一下就好)(a,b)(a,b)(b,a)(b,a)关于直线y=xy=x对称

a,ba,b的任意性可知,y=f(x)y=f(x)y=f1(x)y=f1(x)关于直线y=xy=x对称


初中的三角函数

在讲三角函数之前,我应该先贴一幅图

初中学的三角函数应该仅限于锐角三角形

设三角形ABCABCA,B,CA,B,C所对的边分别是a,b,ca,b,c,那么sinC=ca,cosC=ba,tanC=sinCcosC=cbsinC=ca,cosC=ba,tanC=sinCcosC=cb

除了商数关系tanC=sinCcosCtanC=sinCcosC之外,还有平方关系sin2C+cos2C=1sin2C+cos2C=1

特别的,A=π(rad)2A=π(rad)2时,sinA=aa=1sinA=aa=1,虽然根据常规的方法我们不知道哪条是AA的邻边,但是根据平方关系我们可以解得cos2A=0cos2A=0, 即cosA=0cosA=0,这个时候根据商数关系就没有tanAtanA

然后我们开始讲任意角的三角函数

任意角的三角函数

在平面直角坐标系中作单位圆(就是半径为1的圆啦),AOCAOC的始边OCOCxx轴非负半轴重合,顶点OO与直面直角坐标系顶点OO重合,终边OAOA交单位圆于点AA

易得OAOA = 1

我们可以脑补一下当AOCAOC 为锐角的时候,这个时候sinAOC=ABAO=AB1=AB=yAsinAOC=ABAO=AB1=AB=yA, 同理可得cosAOC=xacosAOC=xa

其实你已经得出了任意角三角函数的一个巧妙的计算式了。

对于任意角AOCAOC, 令OAOA交单位圆于点AA,都有sinAOC=yA,cosAOC=xAsinAOC=yA,cosAOC=xA

可以求得tanAOC=yAxA,AOCπ2,cotAOC=xAyA,AOCπtanAOC=yAxA,AOCπ2,cotAOC=xAyA,AOCπ

诱导公式

这里直接给出公式,推理方面嘛……自己画个单位圆,再画两个角,用终边和单位圆的交点的坐标推理一下就知道了

角度 αα α+2kπ(kZ)α+2kπ(kZ) αα π+απ+α παπα 2πα2πα π2±απ2±α
sinsin sinαsinα sinαsinα sinαsinα sinαsinα sinαsinα sinαsinα cosαcosα
coscos cosαcosα cosαcosα cosαcosα cosαcosα cosαcosα cosαcosα sinαsinα
tantan tanαtanα tanαtanα tanαtanα tanαtanα tanαtanα tanαtanα cotαcotα
cotcot cotαcotα cotαcotα cotαcotα cotαcotα cotαcotα cotαcotα tanαtanα

三角函数的反函数

sin1=asinsin1=asin, cos1=acoscos1=acos

知道这两个应该够了……

补:三角形的面积公式

小学老师是不是告诉过你们

SABC=AB×CD÷2SABC=AB×CD÷2

相信不用我说你也知道作CDABCDABDD

但是老师们基本不会告诉你们SABC=AB×AC×sinBAC÷2SABC=AB×AC×sinBAC÷2

我们可以知道SABC=AC×BD2SABC=AC×BD2

而且我们也知道BD=sinA×ABBD=sinA×AB,

代回去:SABC=AC×sinA×AB2=sinA×AB×AC2SABC=AC×sinA×AB2=sinA×AB×AC2

所以三角形面积也可以等于任意两条边的乘积再乘上夹角的正弦的一半

三角函数的和角公式

正弦

当然sin(α+β)=sinαcosβ+sinβcosαsin(α+β)=sinαcosβ+sinβcosα也是一个很重要的共识

(偷图我们最擅长了)

我们不难发现SABC=SABD=SADCSABC=SABD=SADC

再根据上面推的三角形面积公式:

AB×AC×sin(α+β)2=AB×AD×sinα2+AD×AC×sinβ2AB×AC×sin(α+β)=AD×(AB×sinα+AC×sinβ)AD=cosα×AB=cosβ×ACAB=ADcosα,AC=ACcosβAD2cosαcosβ×sin(α+β)=AD×(AD×sinαcosα+AD×sinβcosβ)sin(α+β)cosαcosβ=sinαcosα+sinβcosβ=sinαcosβ+sinβcosαcosαcosβsin(α+β)=sinαcosβ+sinβcosαAB×AC×sin(α+β)2=AB×AD×sinα2+AD×AC×sinβ2AB×AC×sin(α+β)=AD×(AB×sinα+AC×sinβ)AD=cosα×AB=cosβ×ACAB=ADcosα,AC=ACcosβAD2cosαcosβ×sin(α+β)=AD×(AD×sinαcosα+AD×sinβcosβ)sin(α+β)cosαcosβ=sinαcosα+sinβcosβ=sinαcosβ+sinβcosαcosαcosβsin(α+β)=sinαcosβ+sinβcosα

余弦

(这图有点小)

作单位圆OOxx正半轴于AA,作POA=αPOA=α, ROA=βROA=β连接PRPR ,将PORPOR绕点OO旋转ββAOQAOQ

A(1,0),A(1,0),P(cosα,sinα)P(cosα,sinα), R(cos(β)=cosβ,sin(beta)=sinβ)R(cos(β)=cosβ,sin(beta)=sinβ), AOQ=α+βAOQ=α+β, Q(cos(α+β),sin(α+β))Q(cos(α+β),sin(α+β))

因为旋转所以|QA|=|PR|,QA2=PR2|QA|=|PR|,QA2=PR2

[cos(α+β)1]2+sin2(α+β)=[cosαcos(β)]2+[sinαsin(β)]2cos2(α+β)+12cos(α+β)+sin2(α+β)=cos2α+cos2β2cosαcosβ+sin2α+sin2β+2sinαsinβ[cos2(α+β)+sin2(α+β)]+12cos(α+β)=(sinα+cosα)+(sinβ+cosβ)2cosαcosβ+2sinαsinβ1+12cos(α+β)=1+12(cosαcosβsinαsinβ)cos(α+β)=cosαcosβsinαsinβ[cos(α+β)1]2+sin2(α+β)=[cosαcos(β)]2+[sinαsin(β)]2cos2(α+β)+12cos(α+β)+sin2(α+β)=cos2α+cos2β2cosαcosβ+sin2α+sin2β+2sinαsinβ[cos2(α+β)+sin2(α+β)]+12cos(α+β)=(sinα+cosα)+(sinβ+cosβ)2cosαcosβ+2sinαsinβ1+12cos(α+β)=1+12(cosαcosβsinαsinβ)cos(α+β)=cosαcosβsinαsinβ

这些和角公式都要记住……后面推出来有用的……

复数

向量

我们经常学的应该叫做标量。

有大小没有方向

比如说质量,体积。

但是我们也接触过一些有方向有大小的量

他们叫做矢量。

比如说力,速度。

注:数学里面计算的那个应该叫做速率?

向量也叫矢量,在二维及以上维度既有大小又有方向的量为矢量(矢量定义)

矢量可以用一个点指向另一个点的箭头表示。

一个平面向量也可以用一个数对(x,y)(x,y)表示,表示这个向量从(0,0)(0,0)指向(x,y)(x,y)

当然向量不知平面向量,今天我们的FFT只用二维向量就可以了。所以我们只关注二维向量。

两个向量相等,当且仅当他们的方向,大小都相等。所以起点不是(0,0)(0,0)的矢量也可以用数对表示,把他平移到起点为(0,0)(0,0)就可以表示了

向量和常数(标量)相乘,大小相乘,方向不变。

所以a//bλ,使得b=λaa//bλ,使b=λa

向量的模

|a||a|就是取aa的长度……说完了

向量加法

向量加法满足平行四边形法则三角形法则

三角形法则

(图太难画我就不画了)

假设有三个点:饭堂,学校,考场。

aa从学校指向饭堂,bb从饭堂指向考场, cc从学校指向考场

a+ba+b意味着你先从学校出发,去吃了个饭,然后去慷慨赴死考试

不看过程的话,其实你就是从学校出发去了考场。中途经过了什么地方没有半毛钱关系……

所以a+b=ca+b=c

(等于你没吃饭就上去了)

平行四边形法则

我们可以发现三角形法则的适用范围

当矢量aabb首尾相接的时候才能用

但是阴险狡诈的出题人怎么可能会让你有首尾相接的矢量

现在我们重新定义三个矢量

aa表示向北偏东37°方向走37米,bb表示向正东走4米

那么a+ba+b,就执行这两个指令就可以了……

a+b=ca+b=c, 那么cc指向的点应该就是向北偏东37°走37米,再向东走4米的那个点。

因为先走aa再走bb和先走bb再走aa没有区别……

所以a+b=b+aa+b=b+a

好了回到正题。

向量平移是不会改变大小和方向的

所以我们可以把aa 平移一下。

使得bbaa首尾相接

这个时候如果我们把aaaa的终点连接起来……其实就是一个平行四边形

cc就是平行四边形的一条对角线。

搬运某大佬的图:

矢量乘法

数乘

也叫标量乘法,是矢量和标量相乘

就是多走几个矢量而已……方向不变,长度相乘,前面讲过了。

这里注意,b=λab=λa, λ=0λ=0bb是零向量,零向量没有方向。 λ=1λ=1b=ab=aλ=1λ=1bbaa互为反向量

数量积

也叫点积,是向量和向量的乘积。这里注意符号要用点乘,不能用叉乘

ab=|a||b|cosθab=|a||b|cosθ

θθ为向量aabb的夹角

当然这个也可以是aabb方向上的投影长度。这个是后话……

注意点乘出来的是标量,比如说物理中的做功就是用力的向量乘位移的向量,W=FsW=Fs

向量积

不说了,怕你们混淆。

(其实我也混淆了)

有兴趣的小伙伴自行百度~

复数

是时候回到我们的复数了

虚数

初中老师是不是告诉你们没有2121这东西?

但是高中老师会告诉你们有,而且这东西还很重要

定义虚数单位iii2=1i2=1

那么ii是多少?i=21i=21

初中有没有学过化简二次根式?

2a2b=a2b2a2b=a2b

那么我们也可以推出27=21×7=27i2=27i27=21×7=27i2=27i

这样就变成了一个实数和ii相乘

当然,ii, 2727这种东西,在数轴上面是表示不出来的

所以相对应的,我们把它叫做虚数

2121就是虚数单位

复数

(有没有发现标题里面复数就有好多个……)

类似于有理数可以和无理数相加,实数也可以和虚数相加

18+24=18+2i18+24=18+2i

类似于这样的数就叫做复数。

而且Python/C++里面好像有复数这一概念……

#include <complex>

是可以编译通过的……

x = 18 + 4j
print (type(x))

你会发现它输出了

<class 'complex'>

当然Python里面的虚数单位叫做j……(我也不知道为啥蒙了好久)

共轭复数

相信我,轭念(e4)

共轭复数,指两个实部相等,虚部互为相反数的复数互为共轭复数。

3+2i3+2i32i32i就是共轭复数

特别的,3(+0i)3(+0i)3(0i)3(0i)(就是和自己)互为共轭复数

复数zz的共轭复数记作¯z¯¯¯z

容易得到zׯz=(x+yi)(xyi)=x2y2i2=x2+y2zׯ¯¯z=(x+yi)(xyi)=x2y2i2=x2+y2

容易得到共轭复数关于实轴对称

复平面

不难发现,对于一个复数z=a+biz=a+bi

只要确定了a,ba,b就确定以唯一的复数zz

这就和平面直角坐标系上确定(x,y)(x,y)就可以确定一个点A(x,y)A(x,y)非常像

所以我们整一个复平面:

在复平面内,xx轴上的点都表示实数,yy轴除了原点之外都是纯虚数。

所以一般把xx成为实轴, yy称为虚轴

Z(a,b)Z(a,b)就可以表示复数z=a+biz=a+bi

然后我们发现,Z(a,b)Z(a,b)也可以表示向量OZOZ

所以复数z=a+biz=a+bi也可以用向量OZOZ来表示

模长

复数的模……和向量的模差不多

OZOZ的长度就是|OZ|=|z||OZ|=|z|

不难发现|a+bi|=2a2+b2|a+bi|=2a2+b2

辐角

z0z0的时候,向量OZOZxx的夹角就叫做辐角。

当然也可以理解为以OZOZ为终边的角的度数

这里注意一下区分正角和负角

有了辐角和模长,我们就可以表示一个向量了。

|OA|=r,AOB=θ|OA|=r,AOB=θ,根据三角函数,AB=r×sinθAB=r×sinθ, OB=r×cosθOB=r×cosθ

根据三角形法则,OA=OB+BA=r×cosθ+r×sinθ×iOA=OB+BA=r×cosθ+r×sinθ×i

化简一下,若aa的模长为|a||a|, 辐角为θθ,则a=r(cosθ+i×sinθ)a=r(cosθ+i×sinθ)

同样滴,aa对应的复数AA 也可以用模长和辐角表示

复数乘法

矢量乘法我们没讲过……qwqqwq

但是复数是可以相乘的

因为他们是数。

我们可以用乘法分配律

设复数A,BA,B, A=a+biA=a+bi, b=c+dib=c+di

A×B=(a+bi)×(c+di)=ac+adi+bci+bdi2=ac+adi+bcibd=(acbd)+(bc+ad)iA×B=(a+bi)×(c+di)=ac+adi+bci+bdi2=ac+adi+bcibd=(acbd)+(bc+ad)i

这是用普通方法表示

那么用模长和辐角表示呢?

设复数A,BA,B, A=rA(cosθA+sinθAi),B=rB(cosθB+sinθBi)A=rA(cosθA+sinθAi),B=rB(cosθB+sinθBi)

(前方高能)

A×B=rArB×(cosθA+sinθAi)(cosθB+sinθBi)=rArB×(cosθAcosθB+i×sin θAcosθB+i×sinθBcosθA+sinθAsinθB×i2)=rArB×[cosθAcosθBsinθAsinθB+(sinθBcosθA+sinθAcosθB)i]=rarB×(cos(θA+θB)+sin(θA+θB)i)A×B=rArB×(cosθA+sinθAi)(cosθB+sinθBi)=rArB×(cosθAcosθB+i×sin θAcosθB+i×sinθBcosθA+sinθAsinθB×i2)=rArB×[cosθAcosθBsinθAsinθB+(sinθBcosθA+sinθAcosθB)i]=rarB×(cos(θA+θB)+sin(θA+θB)i)

(根据三角函数和角公式,上面应该有,没有就去百度去)
所以两个复数相乘的结果就是:模长相乘,辐角相加

单位根

在复平面内,其中一个顶点为1(+0i)1(+0i)nn边型的所有顶点称为nn次单位根。

特别的,11次单位根为11, 22次单位根为±1±1

看一下百度百科给的定义

数学上,n次单位根是n次幂为1的复数。它们位于复平面的单位圆上,构成正n边形的顶点,其中一个顶点是1。

如图,ABCABC是正三角形就不证了……

这里OO是单位圆,A(1,0)A(1,0),那么OA,OB,OCOA,OB,OC表示的复数都是三次单位根。

类似的,你也可以求出其他单位根……

我们来研究一下nn次单位根有什么性质

用数学方法,我们可以求出A=1A=1, B=0.5+232iB=0.5+232i, C=0.5232iC=0.5232i

这样我们试验了很多性质,发现好像都不太对……

就差An,Bn,CnAn,Bn,Cn没有试了……

但是写成和的形式……您慢慢拆吧……

我们把它写成模长和辐角表示的形式

我们来回忆一下一个向量怎么写成这种形式

化简一下,若aa的模长为|a||a|, 辐角为θθ,则a=r(cosθ+i×sinθ)a=r(cosθ+i×sinθ)(不错我就是复制的)

AA是单位圆和实轴的交点,BB是辐角为正且最小的单位根,CC是辐角为负且最小的单位根。

这里注意CC的辐角最小指大小,和符号没有关系

A=1×(cosAOA+isinAOA)A=1×(cosAOA+isinAOA)

B=1×(cosAOB+isinAOB)B=1×(cosAOB+isinAOB)

C=1×(cosAOC+isinAOC)C=1×(cosAOC+isinAOC)注意这里AOCAOC应该是一个负角

根据模长相乘,辐角相加,可得

An=1n×[cos(nAOA)+isin(nAOA)]An=1n×[cos(nAOA)+isin(nAOA)]

Bn=1n×[cos(nAOB)+isin(nAOB)]Bn=1n×[cos(nAOB)+isin(nAOB)]

Cn=1n×[cos(nAOC)+isin(nAOC)]Cn=1n×[cos(nAOC)+isin(nAOC)]

我们再看一下单位根的定义,是正nn边形的定点。

那么在正多边形中,中心角该都是相等的。

PS:正多边形的外接圆(或内切圆)的圆心叫做正多边形的中心,外接圆的半径叫做正多边形的半径,内切圆的半径叫做正多边形的边心距.正多边形各边所对的外接圆的圆心角都相等.正多边形每一边所对的外接圆的圆心角叫做正多边形的中心角.正n边形的每个中心角都等于360°/n

也就是AOB=AOC=2π(rad)nAOB=AOC=2π(rad)n

那么nAOB=2π=360=0,nAOC=2π=0nAOB=2π=360=0,nAOC=2π=0

然后sin0=0,cos0=1sin0=0,cos0=1

所以An=Bn=Cn=1An=Bn=Cn=1

那么推广一下,如果是任意一个单位根呢?

我们先设辐角最小且为正的向量对应的单位根为ω1nω1n

因为中心角相等,所以第二个单位根的辐角就是第一个单位根的两倍,模长相乘之后还是1,

因此,我们就可以用ω1n×ω1n=ω2nω1n×ω1n=ω2n表示第二个nn次单位根。

同样滴,我们可以用ωin(0in)ωin(0in)表示第iinn次单位根。

注意ω0n=ωnn=1ω0n=ωnn=1

又因为中心角相等,所以ωinωin的辐角就应该是i×2πni×2πn

模长是1

所以win=1i×(cos2πin+isin2πin)win=1i×(cos2πin+isin2πin)

然后(ωin)n=1in×[cos(2πin×n)+isin(2πin×n)]=cos2πi+sin2πi=cos2π+sin2π=cos0+sin0=1(ωin)n=1in×[cos(2πin×n)+isin(2πin×n)]=cos2πi+sin2πi=cos2π+sin2π=cos0+sin0=1

所以单位根的定义就是:{z|zn=1},nN{z|zn=1},nN
(就是所有正整数次幂能达到1的复数)

性质

  • ωkn=cos2πkn+isin2πknωkn=cos2πkn+isin2πkn,这个上面讲过不说了。
  • ω2k2n=ωknω2k2n=ωkn
    感性理解:ω2k2nω2k2n把圆分成2n2n份,取了2k2k份,就是一整个圆的2k2n=kn2k2n=kn ωkn把圆分成n份,取了k份,就是一整个圆的kn
    理性分析:ω2k2n=cos2π×2k2n+isin2π×2k2n=cos2πkn+icos2πkn=ωkn
  • ωk+n2n=ωkn
    感性理解:多转了半圈,关于原点对称,实部和虚部都要取反。
    理性分析:ωn2n=cos2πn2n+isin2πn2n=cosπnn+isinπnn=cosπ+isinπ=1+0=1;
    ωkn×ωn2n=ωk+n2n,这个是同底数幂相乘……前面说过的
    代入ωn2n=1 得到 ωkn=ωk+n2n
  • ω0n=ωnn=1

好了前方高能,进入正题。

2. 快速傅里叶变换(FFT)

*注:看不懂的可以坚持看下去,或许看到后面你就懂了呢?

讲完了前置知识是不是很开心

那么我们讲单位根是干嘛的呢

如果我们设x=ω1n,那么x2就是第二个n次单位根,x3就是第三个n次单位根……

然后我们就可以用单位根表示一个多项式。

两个多项式相乘……就是单位根相乘……

(好像挺快的)

(看下去再说吧)

为了方便讨论,下文默认n2的整数次幂

前面我们提到过,一个n次多项式可以被n+1个点唯一确定

这些点ai(0in)就是i次项系数

我们设xn次多项式f(x)={ai|0in1}

那么f(x)的值就应该是f(x)=n1i=0aixi

按照次数的奇偶性分类一下f(x)=n22i=0a2ix2i+n22i=0a2i+1x2i+1

注意这里的i不是虚数单位,不要混淆了。

如果我们设f1(x)=n22i=0a2ixi=a0+a2x+a4x2++an2xn22

如果我们设f2(x)=n22i=0a2i+1xi=a1+a3x+a5x2++an1xn22

我们能不能用f1,f2表示出f呢?

肯定是可以的
f1包含了所有系数下标为偶数的单项式,f2包含了所有系数下标为奇数的单项式。

我们可以把x的指数扩大到和系数的下标一样

观察f1,每个系数的下标都是指数的两倍,考虑代入x2

会有f1(x2)=a0x0+a2x2+a4x4++an2xn2

同样的,f2(x2)=a1x0+a3x2++an1xn2

还差了1咋办?

给他乘上呗~

x×f2(x2)=a1x1+a3x3++an1xn1

根据前置知识里面的加法,可以得到f1(x2)+xf2(x2)=f(x)

根据前文,我们把单位根代入回来。

当然那么多单位根,为了确保普遍性,我们代入ωknk的取值我就不多说了……

f(ωkn)=f1(ω2kn)+ωknf2(ω2kn)

没了?怎么可能

还记不记得我们有一条公式ωk+n2n=ωkn

如果k<n2,也就是k+n2<n,那么我们也可以有

f(ωk+n2n)=f1(ω2×(k+n2)n)+ωk+n2nf2(ω2×(k+n2)n)=f1(ω2k+nn)+(ωkn)f2(ω2k+nn)=f1(ω2kn×ωnn)ωknf2(ω2kn×ωnn)=f1(ω2kn)ωknf2(w2kn)

对比一下式子(1)(2) (编号在最右边)

不难发现, 对于确定的 ωkn, f(ωkn)f(ωk+n2n)(k<n2)只有一个常数项 ωknf2(ω2kn) 的符号不同……

所以我们在使用一些花里胡哨的算法枚举f(ωkn), 也就是计算f1(ω2kn)ωknf2(w2kn)的时候,也可以用O(1)的复杂度求出f(ωk+n2n)

这里k的取值范围是0k<n2,那么n2k+n2<n

尝试把小于号变成小于等于。

0k<n2,n2k+n2n1

这就意味着我们只需要求f(ωin)(0i<n2)就可以求出f(ωjn)(0jn1)

这样我们就通过二分把问题缩小了一半

然后我们发现,我们要求的f(ωin)(0i<n2)也可以用二分

也就是说我们可以递归下去……

递归边界就是当i=0的时候,多项式只剩下一个常数项,直接返回就好了

所以FFT相当于二分过程……时间复杂度是O(nlogn)

3. 快速傅里叶逆变换(IFFT)

但是我们FFT求出来的东西

能直接输出吗?

我们平时是很少用复数来表示一个多项式的……

我们借助一个偏门的方法,快速求出了多项式。

结果发现结果你看不懂……

那有啥用?

所以我们要通过快速傅里叶变换的逆变换把多项式转换回来。

前方高能

听不懂的没关系,把结果记下来也可以愉快地写IFFT,反正我也没听懂

假设我们现在通过FFT得到了一些点值{yi|0in1}, 我们要把它变换回多项式系数{ai|0in1}

设有向量ci(0in1)满足

ck=n1i=0yi(ωkn)i

好像又到了推公式的时间……

ck=n1i=0yi(ωkn)i=n1i=0[n1j=0aj(ωin)j(ωkn)i]=n1i=0[n1j=0aj(ωjn)i(ωkn)i]=n1i=0[n1j=0aj(ωijn)(ωikn)]=n1i=0[n1j=0aj(ωijikn)]=n1i=0[n1j=0aj(ωjkn)i]=n1i=0n1j=0aj(ωjkn)i=n1j=0aj[n1i=0(ωjkn)i]

我们设有数列A=x0,x1,,xn1,不难发现则是一个等比数列,公比为x

那么这个等比数列的和S(x)就是

S(x)=n1i=0xixS(x)=ni=1xixS(x)S(x)=(x1)S(x)=xnx0S(x)=xn1x1

如果,这个数列的公比是ωkn

那么

S(ωkn)=(ωkn)n1ωkn1=ωknnωkn1=(ωnn)k1ωkn1=1k1ωkn1=0ωkn1

这里ωkn1是分母不0

所以当ωkn1的时候需要特别讨论

ωkn1,就说明这个数列的公比是1,而且k=0k=n

那么很清楚了吧?这个时候S(ωkn)=S(1)=nx0=n

很明显,S(ωkn)=n1i=0ωin=0

考虑回前面的式子(3)

ck=n1j=0aj[n1i=0(ωjkn)i]

里面有一个n1i=0(ωjkn)i

根据上面讨论的结果,如果jk=0 或者jk=n的时候,(4)的结果是n,否则就是0

当然, 这里jn1,k0,jk<n1<n,所以jk=0,j=k

代入式子(3)

ck=n1j=0(aj×{0jknj=k)

因为我们知道0kn1,0jn1,所以在求和枚举的过程中,!j=k

因为只有当j=k的时候这个单项式才是nak,其他时候都是0,

所以我们就得出了:

ck=nakak=ckn

根据这个,就可以用IFFT把点转回系数

然后就会有人(我)问

这个ck怎么求呢?

我们观察一下它的定义式

ck=n1i=0yi(ωkn)i=(n1i=0yi)×[n1i=0(ωin)k]=(n1i=0yi)×[n1i=0(ωnin)k]

可以发现,对于每一个ck,前面的乘数都是一个固定的常数

只要求后半部分就可以了


总结一下

一张图总结(引用一下远航之曲大佬的图)

诶诶诶别走啊

不看完你知道怎么写代码吗?

不知道怎么写代码学这么多干嘛?

代码实现

递归实现

按照上面说的方法,二分过程求f(ωin)

我可以直接贴代码了!

(好了忘给题目了……)传送门

这里提醒一下,虽然C++有complex库,但还是推荐手写吧……花不了多久的……

还有就是acos(1.0)=π(rad)

数组的话……需要开大一点

具体多大呢?看完代码我们再来推一下

#include <cmath>
#include <string>
#include <cstdio>
#include <cstring>
using std::memset;
const int mod = 998244353;
// #define Cdebug

const double Pi = acos(-1.0);
class complex
{
private:

public:
    complex(const double _real = 0, const double _image = 0) : real(_real), image(_image) {};
    virtual ~complex() {}
    inline double get_real() const { return real; }
    inline double get_image() const { return image; }
    inline void set_real(const double& _real) { real = _real; }
    inline void set_image(const double& _image) { image = _image; }
    complex operator = (const complex&);
    friend const complex operator + (const complex&, const complex&);
    friend const complex operator - (const complex&, const complex&);
    friend const complex operator * (const complex&, const complex&);
    #ifdef Cdebug
    void print();
    #endif // Cdebug
protected:
    double real, image;
};

template <typename T>
const T read()
{
    int x = 0;
    char ch = getchar();
    while(ch < '0' or ch > '9')
    {
        ch = getchar();
    }
    while(ch >= '0' and ch <= '9')
    {
        x = (x << 1) + (x << 3) + (ch ^ 48);
        ch = getchar();
    }
    return x;
}

template <typename T>
void read(T& x)
{
    x = 0;
    char ch = getchar();
    while(ch < '0' or ch > '9')
    {
        ch = getchar();
    }
    while(ch >= '0' and ch <= '9')
    {
        x = (x << 1) + (x << 3) + (ch ^ 48);
        ch = getchar();
    }
}

template <typename T>
void write(const T x)
{
    if(x < 0)
    {
        putchar('-');
        if(x > 9) write(x / 10);
        putchar(x % 10 + 48);
    }else
    {
        if(x > 9) write(x / 10);
        putchar(x % 10 + 48);
    }
}

void FFT(int, complex[], const double = 1.0);

complex a[2097153], b[2097153];
int n, m, limit;
int main()
{
    read(n);
    read(m);
    for(register int i = 0; i <= n; i++) a[i].set_real(read<int>());
    for(register int i = 0; i <= m; i++) b[i].set_real(read<int>());
    limit = 1;
    while(limit <= n + m) limit <<= 1;
    FFT(limit, a);
    FFT(limit, b);
    for(register int i = 0; i <= limit; i++)
        a[i] = a[i] * b[i];
    FFT(limit, a, -1);
    for(int i = 0; i <= n + m; i++)
    {
        #ifdef Cdebug
        a[i].print();
        putchar(10);
        #endif // Cdebug
        printf("%d ", (int)(a[i].get_real() / limit + 0.5) % mod);
    }
    return 0;
}

complex complex::operator = (const complex& a)
{
    this -> real = a.real;
    this -> image = a.image;
    return *this;
}

const complex operator + (const complex& a, const complex& b)
{
    return complex(a.real + b.real, a.image + b.image);
}

const complex operator - (const complex& a, const complex& b)
{
    return complex(a.real - b.real, a.image - b.image);
}

const complex operator * (const complex& a, const complex& b)
{
    return complex(a.real * b.real - a.image * b.image,
                   a.image * b.real + a.real * b.image);
}

void FFT(int limit, complex a[], const double I)
{
    if(limit == 1) return ;
    #ifdef Cdebug
    for(register int i = 0; i <= limit; i++)
    {
        printf("a%d:", i); a[i].print();
        putchar(10);
    }
    #endif // Cdebug
    complex a1[(limit >> 1) + 1], a2[(limit >> 1) + 1];
    for(register int i = 0; i <= limit; i += 2)
    {
        a1[i >> 1] = a[i];
        a2[i >> 1] = a[i + 1];
    }
    FFT(limit >> 1, a1, I);
    FFT(limit >> 1, a2, I);
    register const complex omega(cos(2.0 * Pi / limit), I * sin(2.0 * Pi / limit));
    register complex root(1, 0);
    register const int limits = limit >> 1;
    for(register int i = 0; i < limits; i++, root = root * omega)
    {
        a[i] = a1[i] + root * a2[i];
        a[i + limits] = a1[i] - root * a2[i];
    }
    return ;
}

#ifdef Cdebug
void complex::print()
{
    printf("real:%lf image:%lf %lf", real, image, real / limit);
}
#endif // Cdebug

我们看见洛谷的数据范围1n,m,106

我们要发扬卡空间的优良传统,开数组最多只开大一个

是不是以为开到106+1就够了?

这是开到了106+20的:

事实上,开到3×106是可以满分的

但是我们来卡一下空间

看一下主函数里面的这句:

while(limit <= n + m) limit <<= 1;

结合一下数据范围:n+m2×106

那么limit最大取到多少呢?

log22000000=20.0+

所以limit最大会取到221=2097152

观察代码,我们会发现limit就是数组下标上限

所以我们的数组开到2097153就可以精准卡空间。

迭代优化

以为这就完了?

都不知道谁想出来的迭代优化。

我们知道我们二分的过程中是按照数组下标的奇偶性来分类的

那么我们来画一张二分的图,观察一下下标的二进制

有没有发现,后序列每一个数的二进制就是原序列每一个数的二进制反转

A=A0,A1,,An,n=2k1,kN为原序列,B=B0,B1,,Bn,n=2k为后序列易得A0=0,B0=0

那么我们知道A的通项公式:Ai=i因为原序列就是按顺序的下标。

所以A1=1=(¯00k101)2,根据规律,B1=(¯100k10)2

目前还发现不了什么规律,我们继续看下去

A0=(00k0)2B0=(00k0)2A1=(00k101)2B1=(100k10)2A2=(00k2010)2B2=(0100k20)2A3=(00k2011)2B3=(1100k20)2A4=(00k30100)2B4=(00100k30)2A5=(00k30101)2B5=(10100k30)2

发现规律没有?

没有

我们要求的是数列B,但是BA之间好像没有关系

那么如果我们单独研究序列B呢?

B0=(00k0)2B0=(00k0)2B1=(100k10)2B1=(100k10)2B2=(0100k20)2B3=(1100k20)2B2=(0100k20)2B4=(00100k30)2B5=(10100k30)2

观察一下BiBi2之间的关系

(我都把他们放在一行了)

首先我们要知道在C++中有右移运算符,a >> k相当于a整除2k,而且也很形象,二进制下把所有数字向右移k位,最右边k位舍弃

所以整除2的整数次幂就直接二进制下右移就可以了

好了有了知识铺垫,我们发现,B2iB2i+1,iN,只有最高位不同

B2i的最高位是0, B2i+1的最高位是1

很好,距离规律很近了

再观察一下BiB2i的区别

观察到B2i=Bi>>1

所以我们就得到了B的通项公式

Bi={0i=0Bi2/2+2k1×(imod2)i>0

为什么大于0的时候后面要乘以i2的余数呢?这是因为如果i是奇数,那么他的翻转的最高位就应该是1,否则就是0
所以我们就要加上这一句
写成代码就是
B[i]=(B[i>>1]>>1)|((i&1)<<(limit-1));
得到了B,我们再根据

使用B不断向上合并,就可以得到结果了

那么我们有了原序列,怎么得到后序列呢

注意B是不能修改的,还要留着下次用,我们只能把A改成后序列

我们可以每次比较,如果Ai=i<Bi,那么就交换AiABi,这样就可以得到后序列了

详见代码:

#include <bits/stdc++.h>
using namespace std;

complex <double> a[2097153], b[2097153];
int n, m, limit, l, r[2097153];
double c;
const double Pi = acos(-1.0);

void FFT(const int, complex<double>[], const int = 1);

int main()
{
    scanf("%d%d", &n, &m);
    for(register int i = 0; i <= n; ++i)
    {
        scanf("%lf", &c);
        a[i].real(c);
    }
    for(register int i = 0; i <= m; ++i)
    {
        scanf("%lf", &c);
        b[i].real(c);
    }
    limit = 1;
    while(limit <= n + m) limit <<= 1, ++l;
    for(register int i = 1; i <= limit; ++i)
    {
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    }
    FFT(limit, a);
    FFT(limit, b);
    for(register int i = 0; i <= limit; ++i) a[i] *= b[i];
    FFT(limit, a, -1);
    for(register int i = 0; i <= n + m; ++i)
    {
        printf("%d ", (int)(a[i].real() / limit + 0.5));
    }
    return 0;
}

void FFT(const int limit, complex<double> a[], const int flag)
{
    for(register int i = 0; i < limit; ++i)
    {
        if(i < r[i]) swap(a[i], a[r[i]]);
    }
    for(register int mid = 1; mid < limit; mid <<= 1)
    {
        // omega(cos(Pi * 2 / (mid * 2)), flag * sin(Pi * 2 / (mid * 2)))
        complex<double> root(cos(Pi / mid), flag * sin(Pi / mid));
        for(register int R = mid <<1, j = 0; j < limit; j += R)
        {
            complex<double> omega(1, 0);
            for(register int k = 0; k < mid; ++k, omega *= root)
            {
                complex<double> x = a[j + k], y = omega * a[j + mid + k];
                a[j + k] = x + y;
                a[j + mid + k] = x - y;
            }
        }
    }
}

补充:符号们

当然也要给你们一个目录可以快速回到之前读到的地方

参考资料:还是百度百科

基本数学符号
符号种类 符号 名字 解说 例子 读作
逻辑符号 实质蕴含 A ⇒ B 意味着如果A为真,B也为真。如果A为假,则对B没有任何影响 x = 2 ⇒ x² = 4 为真,但 x² = 4 ⇒ x = 2 一般为假(因为 x 可以是 −2)。 蕴含;如果..那么
可能意味着同 ⇒ 一样的意思(这个符号也可以指示函数的域和陪域;参见数学符号表)。
可能意味着同 ⇒ 一样的意思(这个符号也可以指示超集)。
实质等价 A ⇔ B 意味着 A 为真如果 B 为真,和 A 为假如果 B 为假。 x + 5 = y +2 ⇔ x + 3 = y 当且仅当,iff
¬ 逻辑否定 ¬A 为真,当且仅当 A 为假。 ¬(¬A) ⇔ A
/ 命题逻辑 穿过其他算符的斜线 等同于在它前面放置的"¬"。 x ≠ y ⇔ ¬(x = y)
逻辑合取 如果 A 与 B 二者都为真,则陈述 A ∧ B 为真;否则为假。 n < 4 ∧ n >2 ⇔ n = 3(当 n 是自 然数的时候)。
逻辑折取 如果 A 或 B有一个为真陈述 或二者均为真陈述,则 A ∨ B 为真;如果二者都为假,则 陈述为假。 n < 3 ∨ n > 3 ⇔ n ≠ 3
xor 陈述 A ⊕ B 为真,当且仅当A, B真假性不同。A ⊻ B 意思相同。 (¬A) ⊕ A 总是真,A ⊕ A 总是假。 异或
全称量词 ∀x,P(x)为真,意味着对于所有存在的x,都有P(x)为真 ∀ n, n² ∈ N(n ∈ R) 对于所有;对于任何;对于每个;任意的
存在量词 ∃ x, P(x)为真,意味着存在至少一个x使P(x)为真。 ∃ n ∈ N(n 是偶数)。 存在着
∃! 唯一量词 ∃! x,P(x),意味着精确的有一个x使P(x)为真。 ∃! n ∈ N(n + 5 = 2n) 精确的存在一个
:= 定义 x := y 或 x ≡ y 意味着 x 被定义为 y 的另一个名字(但要注意 ≡ 也可以意味着其他东西,比如恒等)。 cosh x := (1/2)(exp x + exp (−x)) (别问,我也看不懂) 被定义为
:⇔ P :⇔ Q 意味着 P 被定义为逻辑运算上等价于 Q。 A XOR B :⇔ (A ∨ B) ∧ ¬(A ∧ B)
() 优先组合 优先进行括号内的运算 (8/4)/2 = 2/2 = 1, 而 8/(4/2) = 8/2 = 4。 (好像没有)
推论 x ├ y 意味着 y 推导自 x。 A → B ├ ¬B → ¬A 推论;推导
数量符号 π 圆周率 圆的周长和直径的比值 [paɪ]
e 自然常数 [iː]
φ 黄金分割数 把一条线段分割为两部分,使较大部分与全长的比值等于较小部分与较大的比值,则这个比值即为黄金分割。 [faɪ]
i 虚数单位 规定 i²=-1 [ai]
(好像找不到) 毕达哥拉斯常数 根号2 自己算去 Runtime Error
运算符号 四则运算符号 +、-、×或·、÷或/ 懒……
集合运算符号 交集 设A,B是两个集合,由所有属于集合A且属于集合B的元素所组成的集合,叫做集合A与集合B的交集(intersection)。即:A∩B= {x|x∈A∧x∈B}。
并集 若A和B是集合,则A和B并集是有所有A的元素和所有B的元素,而没有其他元素的集合。A和B的并集通常写作 "A∪B",读作“A并B”,用符号语言表示,即:A∪B={x|x∈A,或x∈B}
对数 log 对数(logarithms) 求幂的逆运算.如果a的x次方等于N(a>0,且a≠1),那么数x叫做以a为底N的对数(logarithm),记作x=loga N。其中,a叫做对数的底数,N叫做真数。 [lɒɡ]
lg 以10为底的对数 lg a = log10 a,计算机里面的log其实就是lg lg……吧
ln 以无理数e为底的对数 ln a = loge a ln……吧
lb 以2为底的对数,多用于计算机 lb n = log2 n lb……吧
lim 极限对数 某一个函数中的某一个变量,此变量在变大(或者变小)的永远变化的过程中,逐渐向某一个确定的数值A不断地逼近而“永远不能够重合到A”(“永远不能够等于A,但是取等于A‘已经足够取得高精度计算结果)的过程。 limits……吧
: 懒……
|| 绝对值符号
能用到的基本都在这了……

补充:快读快写

正常的输入输出是不是

cin >> a;
cout << a;

如果是要用文件读入的话

freopen("", "", );
fopen("", "");
fstream fin("");

我们在这里使用

fstream fin("");

对输入输出进行优化


有没有听教练说过iostream比cstdio慢?

为什么呢?

首先我们要知道他们的工作原理。

cstdio的话……我找不到代码了

他们是不断调用底层程序,而且需要有格式化字符串(就是前面你写的一大堆)来识别变量类型。相对来说会快。

但是iostream的话,他们使用的是重定义运算符。而且你有没有发现iostream有很多花里胡哨的操作?

他们所依靠的operator<<函数进行的是格式化输入,这意味着,每次你调用的时候它们都必须做大量工作。它们必须建立和销毁sentry对象(为每个operator<<调用进行建立和清除活动的特殊的iostream对象),

它们必须检查可能影响它们行为的流标志(比如skipws),它们必须进行全面的读取错误检查,而且如果它们遇到问题,它们必须检查流的异常掩码来决定是否该抛出一个异常。

如果进行格式化输入,这些都是重要的活动,但如果你需要的只是从输入流中抓取下一个字符,这样做就过度了。

(不错,复制自CSDN

那咋整?

上面是不是说,他们必须检查很多东西,而且要不断调用IO,所以慢。

那我们就优化一下,使用缓存输入输出。到最后一次性输出,看你还能怎么慢。

posted @   IdanSuce  阅读(602)  评论(0编辑  收藏  举报
编辑推荐:
· 开发者必知的日志记录最佳实践
· SQL Server 2025 AI相关能力初探
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
阅读排行:
· 阿里最新开源QwQ-32B,效果媲美deepseek-r1满血版,部署成本又又又降低了!
· 开源Multi-agent AI智能体框架aevatar.ai,欢迎大家贡献代码
· Manus重磅发布:全球首款通用AI代理技术深度解析与实战指南
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· AI技术革命,工作效率10个最佳AI工具
点击右上角即可分享
微信分享提示