command_block
2018-08-19 08:59:28
——傅里叶变换在信息学竞赛中主要作用是利用卷积思想,化乘为加,快速计算多项式乘法。
我太蒟了,看了
关于多项式,有了
这里不会满屏都是三角函数!!,没有
不需要你会高数,只要用过平面直角坐标系,就可以看。
这篇文章源码40KB,写的时候挑战洛谷Blog(浏览器)性能极限
FFT全称(Fast Fourier Transformation)
中文名：快速傅里叶离散变换
作用 : 以Of course more than that!)。
大家在学多项式乘法的时候,是不是觉得拆括号特别麻烦。
比如说计算
原式
很复杂...
现在我们就要用计算机把两个多项式乘起来(所谓多项式问题)。
P.S:下面提及的多项式均只有一元,你可以认为只有
为了表述方便,下面定义一些记号(请仔细阅读)：
那么
系数(参数)
这里有一个n次多项式
差不多长成这样：
其中
用不同的字母来表示系数很烦,我们就用数组。
通常把
也即
乘法的本质(卷积)
现在如果让你把两个
简单啦！
设
理解:想要乘出
也就是说
也可以写做等价的
形为
那么你可能观察到了,多项式乘法就是加法卷积。
(后面我们还会学习到其他的卷积)
多项式乘法拥有交换律结合律分配率什么的,就不多说了……
亮出模板题
暴力卷积Code(用于上述模板题,44'):
#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;
}要用计算机把两个多项式乘起来就是所谓多项式问题。
数学家为了偷懒发明了麻烦的记号。
首先说一下多项式的点值表达。
如果把多项式看作函数,画出图像：
比如一个一次函数(多项式)
如果一开始不告诉你解析式,只跟你说这个一次函数经过点(0,1)和(1,3),你能搞出解析式吗？
待定系数法(两点确定一条直线)即可。
我们这里不讲待定系数法,我们只讲待定系数法的推论
在平面直角坐标系中,
如果对这方面的知识感兴趣,可以查看 从拉插到快速插值求值
所以我们可以用n+1个点值(有序数对)来描述一个多项式。
比如两个
令数列
把数列
同样的,多项式
比如说
取
函数图像与得到的点如下。
三个红点的函数值分别为
三个绿点的函数值分别为
让你求两个多项式之积
点值表达有一个好处,给你两个多项式点值表达,让你求两个多项式之积的点值表达,直接将对应项乘起来就好了。
(点值的乘法对应着整个多项式的乘法,也就是浓缩了整个多项式)
比如上文
就相当于告诉我们
那么
废话,肯定等于
所以
所以W在
验算一下？
代入
得到
总结一下就是
好像那里不对,W的次数明摆着是2n。
在平面直角坐标系中,给你(n+1)个点就能确定一个n次多项式。
所以需要2n+1个点才行,还差n个,肿么办？
反正点都是你自己找的,再多来几个不就好了……
现在
所以W的点值表达为
只需要做2n+1次乘法即可。
看起来很好用,其实只是几个点而已,和系数表达差着十万八千里。
于是我们就想,能不能把系数表达转换为点值表达呢？
算法流程：
把系数表达转换为点值表达->点值表达相乘->把点值表达转换为系数表达。
这就是 FFT 的算法流程。
“把系数表达转换为点值表达”的算法叫做DFT
“把点值表达转换为系数表达”的算法叫做IDFT(DFT的逆运算)
P.S:
从一个多项式的系数表达确定其点值表达的过程称为求值(毕竟求点值表达的过程就是取了 n 个 x 然后扔进了多项式求了 n 个值出来);
而求值运算的逆运算(也就是从一个多项式的点值表达确定其系数表达)被称为插值.
在平面直角坐标系中,给你(n+1)个点就能确定一个n次多项式(函数)。
点值表示相乘(点乘)远快于暴力卷积。
(如果你不知道
DFT有一个妙处,就是代入什么由你自己决定,只要点值个数足够。
我们这些蒟蒻第一感觉都会选择人畜无害的正整数,或者有理数什么的。
但是这些东西虽然在普通人看来计算简单,但并没有什么能够加以利用的优秀性质。
事实证明,找一些毒瘤东西代入进去是个好想法。
上古之时,有一位dalao横空出世,他就是傅里叶！
好像是个搞复数的专家,有一天他突发奇想：
把单位根什么鬼)
然后蛤玩意儿)
然后分治,就
好吧先来介绍复数。
复数,即形为
比如说数轴原来是一条直线：
某个实数一定可以用数轴上一个点来代表。
多么和谐自然。
然而复数来了。
数轴原来是一条横着的直线,现在又多了一条竖着的,变成了一个坐标系。
原来的那条直线(横)称为实轴
新的那条直线(竖)称为虚轴
每个复数一定可以用这个坐标系上一个点来代表。
比如蓝点是
精妙之处在于,复数之间能做运算。
不懂的话就多琢磨下我写完这些之后就懂了
我们来介绍下复数运算：
复数相加：实部和虚部分别相加。
如：
复数相减：取个相反数再相加。
复数相乘：像一次多项式一样相乘。
如：
P.S:
复数的共轭：实部相同,虚部相反(根据实轴对称)。
如
P.S:一个复数与自己的共轭相乘一定会得到实数：
复数相除：分子分母同乘分母的共轭
将分母实数化,在分子分母同时乘上分母的共轭,即可把分母变成实数。
然后直觉化简。
而且这些运算在几何(也就是上面的坐标系)里面也有神奇的性质。
这里画的是复数
看出来什么没有？
好吧,看见那些从
以你做几何题的直觉？？
没错！如果设原点为点
有性质 :
我们把一个复数表示的点到原点的距离叫做这个复数的模长,也就是这里的
也即 : 复数相乘时,模长相乘,幅角相加！
原点到该点的射线与实轴正方向射线组成的角 （逆时针旋转） 的角度乘坐这个复数的幅角,复数
(不清楚的看图)
接下来是证明(最好跟着把图画出来)：
为保证一般性,下面涉及到的都是字母。
模长相乘
设两个复数为
那么它们表示的点就是
这两个复数的积为
证毕。
幅角相加
第二个规律需要相似+爆算的,不想了解请跳过。
那么,如果能证明途中两个三角形相似,就能证明幅角相加定理了。
我们已经证明了
证毕。这玩意特别重要。
接下来我们介绍单位根。
看起来很高大上,这是傅里叶快速变换的重要内容。
换句话说,就是方程
到这里你也许会有一点不懂,我们来一起推导一下。
学过三角函数的话应该对单位圆很熟悉。
不熟悉的话也没事,单位圆就是：
我们在复平面上放一个单位圆,在这个单位圆上的点表示的复数都有一个性质：
模长为1(圆上每一点到圆心的距离都相等)
可还记得？n次单位根是n次幂为1的复数,1的模长就是1。
考虑一个复数
如果
如果
所以只有模长等于1的复数才有可能成为n次单位根。
也就是说,只有单位圆上的点表示的复数才有可能成为n次单位根(必要性)
我们成功缩小了探究范围,只在这一个圆上研究。
(下面提及的复数,均是在单位圆上的复数!!!)
接下来,我们要找单位根。
这些单位根的模长都是
容易找到 : 幅角为
也就是说,假如一个复数,其幅角为
容易知道这玩意的
根据代数基本定理 : n次方程在复数域内有且只有n个根。
所以这些不重不漏,就是所有的单位根了。
还不懂也没关系,我们来实践一下。
设
有一个模长为1的复数为3次单位根,它的幅角为(*圆周$a$**)
它的3次方的幅角为
那么,如果它是3次单位根,幅角的3倍
所以3a是自然数。
a=0时,它的幅角为0(其实这个复数就是1)
a=1/3时,它的幅角为(圆周/3)
a=2/3时,它的幅角为(2圆周/3)
有幅角为
每一个都比上一个"多转了"
所以n次单位根n等分单位圆。(重要)。
怎么称呼这些单位根呢
“颜文字请走开,LaTeX是$\omega$)
我们从1开始(1一定是单位根),沿着单位圆逆时针把单位根标上号。
显而易见,不懂的把上面的图看一看。
0.
每乘上一个
转了k次,即为
1.两个n次单位根
由(0)可得
感性理解 : 我们把一个圆等分成n份(想想你切蛋糕的时候)
由(1)可得:
根据
根据
显而易见,
2.
3.如果n是偶数,
根据
根据(1),把“乘
也就是旋转半个圆周,也就是关于原点对称,也就是取相反数(建议画个图)
好啦,关于单位根想必你已经掌握啦(终于！！)。
坐稳了,前方高能！！！
复数把数轴扩展到了复平面上,复数可以对应复平面上一个点。
复数也有四则运算。复数相乘时,模长相乘,幅角相加。
n次单位根(n为正整数)是n次幂为1的复数。
把单位根画到单位圆上之后,就能整出一些性质。
我们来讲讲FFT的分治方式
傅里叶一巴掌把这个多项式拍成碎片,然后按奇偶拼成两块
考虑n-1次(也就是说有n项)多项式F(x)
的每一项按照下标(次数)的奇偶分成两部分：
这里保证n是2的整幂,不会出现分得不均匀的情况。
又设两个有n/2项的多项式
他们的系数依赖于
(又把多项式忘了的,赶紧回去看)
则可以得出
我们要把
因为
(下一步用到了
(下一步用到了
(下一步用到了
(下一步用到了
最后得到
比对一下上一个式子
等式右边不就是只有一个正负号的区别吗？
这两条式子有什么用呢？
到这里,如果我们知道两个多项式
套用
套用
所以如果我们知道两个多项式
我们就成功的获得了
(懵逼不怕,具体见代码)
好像有哪里不对？
上文: 如果我们知道两个多项式
问题在于我们不知道啊？
对
NO~NO~NO~
上面的工作,其实是一个分治的过程。
FFT分治的过程,就是根据：
这两个式子,实现子问题的合并。
记忆方法:
还有一个小细节。
上文有一句话：“保证n是2的整幂,不会出现分得不均匀的情况。”
实际应用中,n不一定是2的正整数次幂。
我们可以补项,在最高次强行添加一些系数为0的项(类似于高精度补0)。不影响我们的计算结果,却占了位置。(具体见代码)
讲完了这些我们可以开始写
根据上文复数的四则运算重载。
CP是
#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,会被卡常)
之前我们一直在说“把单位根代入...”
现在我们来学习如何求单位根。
前面我们说过
也就是说只要求出
我们怎么求
点A就是
在
其他什么也不知道。我们需要求
百度百科
| 英文缩写 | 表达式 | 语言描述 | 函数名称 |
|---|---|---|---|
| sin | a/c | ∠A的对边比斜边 | 正弦函数 |
| cos | b/c | ∠A的邻边比斜边 | 余弦函数 |
| tan | a/b | ∠A的对边比邻边 | 正切函数 |
观察上图,需要求
（模长
所以
C++的三角函数采用弧度制。一个整圆不再是
所以
只要求出
代码：
#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;
}具体见代码,多说无益。
#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)
(2)
现在问题来了,DFT输出的都是些杂乱的点值表达,所以解决不了问题。
上文说过,IDFT是DFT的逆(CP),她可以把点值还原成多项式,最终完成乘法。
现在到了讲IDFT的时候了!!!
IDFT和DFT就是两句话的区别,这个结论实在太好记了。
结论 : 把DFT中的
当然,证明还是要有的,而且很有价值,从这里可以延伸出 单位根反演。
我们知道
还是多项式
即
脑补代入的过程,可得
那么结论就是
通过这个式子我们可以把点值还原。
证明:
右边
代入可得:
分类讨论:
那么贡献就是
那么贡献就是
等比数列求和可以得到:
所以这部分贡献为0,我们就成功地证明了上述结论。
如果你想知道的多一些,我可以告诉你这个东西本质是单位根反演。
或者可以理解成范德蒙德矩阵求逆。
就是相当于把
不同的是,这次代入的是
这里就要求找到一个东西能求出
也就是说只要求出
它的各方面和
题目Link
至此我们已经写出了第一个版本的FFT(蛤？还有几个？)
代码很好写,和DFT不同的地方很少：
#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;
}提交记录
可以看到因为常数过大而TLE了,看来我们需要一个常数更小的写法。
先别着急卡常,前文我们说过,IDFT和DFT只有一个符号的区别。
那么我们何不减少一下代码量呢:
#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;
}卡常的方法分为以下几类：
0.硬件优化：浮点加速器,CPU等等。
1.编译优化：编译开关,inline,register 之类的。
2.封装优化：把代码的封装拆开,封装得过死会使效率降低(stl)。
3.编写优化：细节优化,手写栈之类。
4.时空互易: 空间换时间以达到时空均衡(常数层面上)。
硬件优化没什么指望,编译优化能力有限又太玄学。
咱们这份代码如果硬要把封装拆开就会显得十分丑陋。
只剩下编写优化,时空互易两条路子了。
首先观察这个程序片段:
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]我们计算了两次!
换成下面的操作,常数就小多了:
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”确实排在原来“3”的位置。
如何得到二进制翻转后的数列呢？
假设我们需要求
设复多项式
则
发现
也就是说求出
Code:
为了卡常数还写了个快读……
#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;
}相信经过上面的学习,这份代码大家都能够理解,我就不再赘述了.
评测记录
快了一点吧！虽然本蒟蒻的常数还是很大QwQ
比如说
直接做
假如使用三次变两次优化,由于平方项的存在,涉及的精度跨度上限是
下方例题第二题直接铁头硬上会WA。
至于修正方法,可以把两个多项式数乘到相同的值域范围,这样子就不会产生平方项的大误差了。
注意最后要除回去。
我们的点值代表了整个多项式,加法和乘法都是对应的。
这可以用于某些(初等)卡常,比如我们如果要计算
如果完全封装,每次乘法就要计算
化为
我们也不必对
又注意到
采用"三次变两次"优化,构造
平方之后也不需要分别IDFT,这样做是
卡常有什么好总结的？去吊打集训队问wys吧。
A*B Problem升级版
力
那么讲到这里,想必FFT的大部分内容你已经掌握了,到了说再见的时候了……
再见啦~(挥手)
欲知后事如何,请看下集:NTT与多项式全家桶! (看完你会了解更多的多项式工业)