目录
总所周知，FFT是一个非常麻烦的算法，再加上博主语文不好，便写起来有点麻烦，但会尽力去写。要以后自己看不懂就。。。
注：因为最近的压力紧张，便没有继续学习FFT，这仅为目前的半成品以及一些目前已知的优化。
参考资料
毛啸大佬的论文（不放上下载地址，因为下载论文的主要OJ太卡）
因为此篇文章比较依赖这些文章，为此先放出来。
FFT
吹水
FFT到底是干什么的？法法塔清小怪，不不不，人家的全名叫快速傅里叶变换。听说与傅里叶没有半毛钱关系。
FFT能解决两个多项式相乘。
例如：
\(a_{2}x^{2}+a_{1}x^{1}+a_{0}\)与\(b_{2}x^{2}+b_{1}x^{1}+b_{0}\)相乘，便可以算出\(x^{4},x^{3},x^{2},x^{1},x^{0}\)的系数，在\(O(nlogn)\)的时间，而暴力要\(O(n^{2})\)。
例题
普通做法
上联：打表过样例，暴力出奇迹。
下联：花中一片叶，luogu两行泪。
。。。
更高大尚的做法
于是，就有人开始尝试新的做法了！
说到底，FFT究竟是什么？
大概样子：
系数表达式->点值表达式->点值相乘->系数表达式
现在，就为大家一一刨析做法！
定义与一部分性质
没学过多项式的估计你还没学OI呢，神童？默认大家会一点多项式的定义。
咱们先理清我说的一些十分玄学的话。
系数表达式
设\(A(x)\)表示一个n−1次多项式
则\(A(x)=p = \sum_{i=0}^{n} ai∗xi\)
例如：\(A(3)=2+3∗x+x^{2}\)
这就是个系数表达式。
---以上摘自(http://www.cnblogs.com/zwfymqz/p/8244902.html#_label1_1)
点值表达式
这个就比较玄学，其实每一个n次的系数表达式都可以用n+1个点表是，将n+1个不同的值代入系数表达式，得出了n+1个点值\((x,y)\)，而这n+1个点值就是这个系数表达式的点值表达式。
性质1
为什么n+1个点可以表示n次表达式，年轻人，我跟你共呀。
我们列一个n次方程：\(a_{n}x^{n}+a_{n-1}x^{n-1}+...+a_{1}x^{1}+a_{0}=y\)，将n+1个x,y带进去，列出一堆a的连七八糟的方程，用高斯消元便可以解出这个方程了。
点值相乘？？？
首先，系数相乘要\(n^{2}\)，但是我们想想看点值相乘有没有可能乘出系数相乘后的表达式的点值？如果可以的话，岂不是\(O(n)\)乘完？
明显，传说有个大佬想到了这一点，于是轻微的证了一下：
\(C(x)=A(x)*B(x)\)
所以
\(C(0.5)=A(0.5)*B(0.5)\)
这与点值有什么关系？
把0.5带到A(x)与B(x)不就是点值了，相乘不就是A(0.5)*B(0.5)了想乘了？
证毕。
但是，由于乘积后的多项式是2n次多项式，因此原本的两个多项式也需要准备2n+1个点来得出最后的多项式。
只可惜，点值相乘确实牛逼，但是转成点值并且如何转回来更难
卷积
两个多项式的乘积就是他们的卷积。
复数
数轴上的点是实数，但是二维坐标系中的点如何用一个表达式表示呢？一个巨人一巴掌拍成了数轴就可以了。。。
我们看一张图：
看，这个点的x坐标可以用实数表示，但y坐标。。。
还是(x,y)友好
于是就有大佬发明了个单位i，就是传说中的复数单位，i有个性质，\(i^{2}=-1\)。
于是这个点的表示就是4+5i
于是大佬好人做到底，送佛送到西
大佬们又发明了复数运算。
仅介绍与FFT有关的运算。
复数相加：\((a+bi)+(c+di)=(a+c)+(b+d)i\)
复数相减：\((a+bi)-(c+di)=(a-c)+(b-d)i\)
复数相乘：\((a+bi)+(c+di)=ac+adi+bci+dbi^{2}=(ac-bd)+(ad+bc)i\)
复数的共轭：\(a+bi=a-bi\)
---以下来源于https://www.luogu.org/blog/command-block/fft-xue-xi-bi-ji
单位根
其实，FFT就是DFT加上IDFT
DFT是什么
DFT就是一个能在\(O(nlogn)\)内将系数表达式转成点值表达式
我们平常把系数表达式转成点值表达式不就是一个个代，然后就\(O(n^{2})\)了
但是，其实有一个神奇的东西，估计大家很熟悉，叫单位根。
放入大佬（https://www.luogu.org/blog/command-block/fft-xue-xi-bi-ji）的讲解：
n次单位根(n为正整数)是n次幂为1的复数！
到这里你也许会有一点不懂，我们来一起推导一下。
换句话说，就是方程\(x^{n}=1\)的复数解。
学过三角函数(解三角形)的童鞋应该对单位圆很熟悉
不熟悉的话也没事，单位圆就是：
圆心在原点且半径为1的圆(如图)
我们在复平面上放一个单位圆，在这个单位圆上的点表示的复数都有一个性质：
模长为1(圆上每一点到圆心的距离都相等)
可还记得？n次单位根是n次幂为1的复数
1的模长就是1
想象一个复数 xx 。
如果 \(|x|>1\) ，即\(|x^{n}|=|x|^{n}>1\) (模长乘模长，越乘越大)
如果 \(|x|<1\)， 即\(|x^{n}|=|x|^{n}<1\) (模长乘模长，越乘越小)
所以只有模长等于1的复数才有可能成为n次单位根。
所以只有模长等于1的复数才有可能成为n次单位根。
也就是说，只有单位圆上的点表示的复数才有可能成为n次单位根
我们成功缩小了探究范围，只在这一个圆上研究。
(下面提及的复数,均是在单位圆上的复数!!!)
想象一个模长为1的复数为 \(X\) ，设它的幅角为(圆周 \(*a\) )
(这里 \(0<=a<1\) ,不然会产生重复)
它的n次方的幅角为 \(na∗\) 圆周
那么，如果它是n次单位根，幅角的n倍 ( \(na∗\) 圆周 ) 一定是圆周的自然数倍。
得到不定方程(听到这个名字别害怕)：
(\(na∗\) 圆周 ) = 圆周 \(∗k\) (k 为自然数)
等式两边同除以“圆周”
得 \(na=k\) (k 为自然数参数)
分类讨论一下：
\(k==0\)
由于 \(n>0\) (方程的次数是正整数)只有 \(a=0\) 这一种情况
\(k>0\)
\(k\) 可以为 \(0,1,2,3...\)
\(k\)为1时,得 \(na=1,a=1/n\)
\(k\)为2时,得 \(na=2,a=2/n\)
...
\(na=k, a=k/n\)
但是， \(k>0\) 时并没有无数个单位根，他们会重叠，重叠的只算一个。
如 \(k=n\) 时,得 \(na=n,a=n/n=1\)
这就相当于绕了单位圆一圈又回到了起点，与k为0时没什么不同。
前面说 \(0<=a<1\) ,到这里就该打住。
所以本质不同的单位根是 \(0<k<n\) 的 \(n-1\) 个单位根,共\(n\)个。
综上所述， \(k==0\) 时必有一个单位根， \(k>0\) 时有\(n-1\)个单位根，所以\(n\)次单位根一共有\(n\)个。
另一种理解方法:
根据代数基本定理,n次方程在复数域内有且只有n个根。
上面我们验证了,幅角为 \(0,\dfrac{1}{n}\text{圆周},\dfrac{2}{n}\text{圆周},...,1\text{圆周}\)的单位根,共n个。
所以这些就是所有的单位根了。
还不懂也没关系，我们来实践一下。
设 n=3;n=3;
有一个模长为1的复数为3次单位根，它的幅角为(圆周 \(∗a\) )
它的3次方的幅角为 \(3a∗\) 圆周
那么，如果它是3次单位根，幅角的3倍 (\(3a∗\) 圆周 ) 一定是圆周的自然数倍。
所以3a是自然数。
a=0时，它的幅角为0(其实这个复数就是1)
a=1/3时，它的幅角为(圆周/3)
a=2/3时，它的幅角为(2圆周/3)
可以看出,\(n\)次单位根\(n\)等分单位圆。
证明:
有幅角为 \(0,\dfrac{1}{n}\text{圆周},\dfrac{2}{n}\text{圆周},...,1\text{圆周}\)的单位根,共n个。
每一个都比上一个"多转了" \(\dfrac{1}{n}\text{圆周}\)。
所以n次单位根n等分单位圆。
(非常重要!)。
怎么称呼这些单位根呢
“\(\omega\)” 是单位根的符号
我们从1开始(1一定是单位根)，沿着单位圆逆时针把单位根标上号。
\(w_{n}^{0}=1\)
\(w_{n}^{1}\)为第二个n次单位根
...
\(w_{n}^{n-1}\)为第二个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\) 的情况(学过带正负的角吗)。
所以 \(\omega_{n}^k=\omega_{n}^{k\%n}\)。
随意称呼,切了书写方便,一般不取模。
关于单位根的基本定义相信你已经Get到了。
单位根的世界，就是一个单位圆。
摘抄完毕。（我不是博主，我只是博客的搬运工）
当然单位根还有一系列要证的，再次，我形容一个抽象的概念帮助大家更好的食用。
我们可以把n次单位根平均切成n块，例如四次单位根：
而且我们可以证明它，不是说了吗，\(w_{n}^{i}\)的n次方是1吗，也就是图中被圈的位置，我们把被圈上的单位根叫\(w_{n}^{0}\)（注：\(w_{n}^{0}=1\)），然后逆时针标上\(w_{n}^{1},w_{n}^{2},w_{n}^{3}...\)
根据图中，我们可以知道，\(w_{n}^{1}\)的幅角为\(1*\)\(360\over n\)度，\(w_{n}^{2}\)的幅角为\(2*\)\(360\over n\)度...\(w_{n}^{i}\)的幅角为\(i*\)\(360\over n\)度，那么\(w_{n}^{i}\)的n次方的幅角就是\(w_{n}^{i}\)的幅角乘以n，容易知道为360的倍数，所以就在\(0,1\)的位置，所以可以知道，单位根连接原点的这几条线均分单位圆，用大佬的话来讲，就是将“蛋糕”分成了n份。
性质1：\(w^{i}_{n}*w^{j}_{n}=w_{n}^{i+j}\)
我们想呢，取i份蛋糕再取j份蛋糕不就是i+j份蛋糕？（复数乘法幅角相加）
性质2：\((w^{i}_{n})^{k}=w_{n}^{ik}\)
我们想呀，\(w_{n}^{i}\)就是取出“蛋糕”中的i份，而k次方就是k次取出其中的i份，不就是ki份了吗？
性质3：\(w_{kn}^{ki}=w_{n}^{i}\)
显而易见，切kn刀拿ki份等价于n刀i份。
性质4：\(w_{n}^{k+n/2}=-w_{n}^{k}\)(\(n\)被2整除)
这个比较抽象（貌似其他也很抽象），就是你现在有一个单位根，而加上了所谓的\(n\over2\)块蛋糕后，就相当于取了目前单位根的负数（旋转180度），给上两张图让大家理解：
自己画几个图便可以理解。
于是，在千辛万苦之下，我们懂了单位根。
DFT
首先，我们不说一些BB的话了。
首先，DFT的本质就是分治同时利用数学的奇妙性质实现\(O(nlogn)\)的。
？？？
别急，我们一步步来，首先，我们把一个n次多项式（n可以写成\(2^i-1\)的形式，也就是说里面有\(2^{i}\)个项）按奇偶分开。
如：\(a_{0}+a_{1}x^{1}+a_{2}x^{2}+a_{3}x^{3}\)可以分成\(a_{0}+a_{2}x^{2}\)和\(a_{1}x^{1}+a_{3}x^{3}\)
我们设\(F(x)=a_{0}+a_{1}x^{1}+a_{2}x^{2}+a_{3}x^{3}=(a_{0}+a_{2}x^{2})+(a_{1}x^{1}+a_{3}x^{3})\)，我们再设两个函数：\(FL(x),FR(x)\)。
有：\(F(x)=FL(x^{2})+x*FR(x^{2})\)。
\(FL(x),FR(x)\)是依赖于\(F(x)\)的。
\(FL(x)=a_{0}+a_{2}x+.......\)
\(FR(x)=a_{1}+a_{3}x+.......\)
那么，我们就可以分治啦。
不是不是呀，我们只是分治了项数呀，但是要找的点数却不变，所以分治不顶用呀！
急什么。
不如重点。
性质1：\(F(w_{n}^{i})=FL(w_{n/2}^{i})+w_{n}^{i}*FR(w_{n/2}^{i})\) \((i<n/2)\)
我们怎么证？
\(F(w_{n}^{i})=FL((w_{n}^{i})^{2})+w_{n}^{i}*FR((w_{n}^{i})^{2})\) \((i<n/2)\)
因为\((w_{n}^{i})^{k}=w_{n}^{ki}\)
\(F(w_{n}^{i})=FL(w_{n}^{2i})+w_{n}^{i}*FR(w_{n}^{2i})\) \((i<n/2)\)
因为\(2|n\)并且\(w_{2n}^{2i}=w_{n}^{i}\)
\(F(w_{n}^{i})=FL(w_{n/2}^{i})+w_{n}^{i}*FR(w_{n/2}^{i})\) \((i≥n/2)\)
这样我们就利用\(n/2\)次单位根推出了目前的\(w_{n}^{0},w_{n}^{1},w_{n}^{2}......,w_{n}^{n/2-1}\)
但是另一半呢？
性质2：\(F(w_{n}^{i+n/2})=FL(w_{n/2}^{i})-w_{n}^{i}*FR(w_{n/2}^{i})\) \((i<n/2)\)
\(F(w_{n}^{i+n/2})=FL((w_{n}^{i+n/2})^{2})+w_{n}^{i+n/2}*FR((w_{n}^{i+n/2})^{2})\) \((i<n/2)\)
因为\((w_{n}^{i})^{k}=w_{n}^{ki}\)
\(F(w_{n}^{i+n/2})=FL(w_{n}^{2(i+n/2)}))+w_{n}^{i+n/2}*FR(w_{n}^{2(i+n/2)})\) \((i<n/2)\)
因为\(2|n\)并且\(w_{2n}^{2i}=w_{n}^{i}\)
\(F(w_{n}^{i+n/2})=FL(w_{n/2}^{i+n/2}))+w_{n}^{i+n/2}*FR(w_{n/2}^{i+n/2})\) \((i<n/2)\)
因为\(w_{n}^{i}=w_{n}^{i\% n}\)，\((i<n/2,i+n/2>n/2)\)
\(F(w_{n}^{i+n/2})=FL(w_{n/2}^{i}))+w_{n}^{i+n/2}*FR(w_{n/2}^{i})\) \((i<n/2)\)
因为\(w_{n}^{k+n/2}=-w_{n}^{k}\)
\(F(w_{n}^{i+n/2})=FL(w_{n/2}^{i}))-w_{n}^{i}*FR(w_{n/2}^{i})\) \((i<n/2)\)
我们就完美的解决了另外一半的点，而且我们发现这些只有正负号的差距，于是这两条式子就可以帮助我们一步步分治了。调试也难了。。。
于是，我们就再\(O(nlogn)\)的到了点值表达式
IDFT
IDFT就是一个可以将点值表达式转成系数表达式的东西，时间复杂度\(O(nlogn)\)
如何转？拉格朗日插值法其实我不会拉格朗日插值法？高斯消元？
貌似都可以，但是都很慢。争做世界上最快的男人
于是，又有大佬想出了方法，因为证明过程太过牛逼作者太懒，所以放上一名大佬的证明
因此，DFT与IDFT大家都会了，就到了高潮重点了。
蝴蝶迭代优化
因为FFT的分治过程是用DFS，而DFS会严重拖慢速度，于是有大佬想到了蝴蝶迭代，使得 FFT拜托了DFS的“控制”
因为作者懒作者语文不好，为了更好的食用，放上大佬的博客摘录。
单位根求法
不是，单位根怎么求呀，在此之前，请大家记一下三角函数。
尤其是正弦函数(cos)与余弦函数(sin)，在C++的cmath库内有函数。
B点便是我们要求的了，那么他的坐标就是\((AC,BC)\)怎么求？我们知道AB长为1，同时
于是，大部分我们都准备好了，可以开始实现了。\(∠ACB=90°,∠BAC\)也可以求出来，那么利用三角函数易得知\(CB=sin(∠BAC)*AB,AC=cos(∠BAC)*AB\)，同时\(AB=1\)，所以\(CB=sin(∠BAC),AC=cos(∠BAC)\)。
但是C++里面的三角函数采用弧度制下一个整圆不再是 \(360\) °,而是 \(2π\) 了。
所以\(∠BAC=\)\(2π\over n\)，而不是\(∠BAC=\)\(360°\over n\)了
当然，万事皆有烦人之事，\(w_{n}^{1}\)在\((0,1)\)，自己画图就知道三角函数这时候就很惆怅了，有两种方法：1.霸王硬上弓，一般不影响结果。2.特判，更稳，但要记。
实现、细节与小优化
就到了我们的实现细节了。
细节
π
\(π\)等于多少？
当然，你可以去背，再这打上：3.1415926535897932
也可以用acos(-1)，精度很好，不过最后一位错了。
反正我是背的，当然，听说用acos(-1)更好。
补齐？
我们之前也讲了，DFT的分治基于的是\(n\)（\(n\)次多项式）能写成\(2^{i}-1\)的形式，也就是多项式中有\(2^{i}\)项，但是题目中并不能保证，我们可以手动补齐，计算离\(2n+1\)（因为结果的最大次为\(2n\)）最近且比他大的2的幂次。
如：\(n=5\)因补齐到15。
合并
通过看大佬的博客，相信大家也会有所了解。
IDFT的代码和DFT的代码可以合并，但是需要注意的是IDFT中的单位根都是\(w_{n}^{-1}\)什么的，于是我们的\(y\)坐标要标成负数，当然，我们只需要多传一个\(type\)值就行了。（注意，IDFT的结果要除以多项式的项数\((limit)\)）
预处理蝴蝶迭代？
for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((n&1)<<(l-1));
//l代表的是离2n+1最近且比他大的2的幂次是2的几次方。复数都去哪了？
公式都推出来了，怎么可能不对？
机房大佬：都推出来了，怕什么，你这是不自信的体现，反正原本都是实数，推回去应该也是实数呀。
小优化
优化不能少呀。
作为毒瘤选手为了爆破毒瘤出题人的癖好，我们一定要学会小优化。
我们可以预处理出cos和sin的值，毕竟我预处理了以后我的速度就猛如虎了。
实现
强势代码来一波
我没开O2
// luogu-judger-enable-o2
#include<cstdio>
#include<cstring>
#include<cmath>
#define pai 3.1415926535897932
#define N 2500000
using namespace std;
inline int getx()
{
int f=1,x=0;char c=getchar();
while(c>'9' || c<'0')c=='-'?f=-1:0,c=getchar();
while(c>='0' && c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*f;
}//我怎么可能会用快读呢？
int n,m,limit/*这就是传说中离2n+1最...的2的幂次*/,l/*之前说过了*/,r[N]/*预处理蝴蝶迭代*/;
double mcos[30],msin[30];//预处理cos和sin
struct node
{
double x,y;
node(double xx=0,double yy=0){x=xx;y=yy;}//一个特别神奇的操作，在node里面打上=代表如果没有传参数就等于这个数
}a[N],b[N];node wn[1200000];
inline node operator+(node x,node y){return node(x.x+y.x,x.y+y.y);}
inline node operator-(node x,node y){return node(x.x-y.x,x.y-y.y);}
inline node operator*(node x,node y){return node(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
inline void swap(node &x,node &y){node z=x;x=y;y=z;}
//复数的一些操作
void fft(node A[],int type)
{
for(int i=0;i<limit;i++)
{
if(i<r[i])swap(A[i],A[r[i]]);
}//预处理基层
for(int id=0;id<l;id++)
{
int mid=(1<<id);
wn[1]=node(mcos[id],msin[id]*type);
for(int k=2;k<mid;k++)wn[k]=wn[(k>>1)]*wn[(k>>1)+(k&1)]/*用这种方法可以尽量减少每个单位根乘的次数，减少误差*/;
for(int R=mid<<1,L=0;L<limit;L+=R)
{
for(int k=0;k<mid;k++)
{
node x=A[L+k],y=wn[k]*A[L+mid+k];
A[L+k]=x+y;A[L+mid+k]=x-y;//之前推的两公式
}
}
}
}
int main()
{
wn[0]=node(1,0);
n=getx();m=getx();
for(int i=0;i<=n;i++)a[i].x=getx();
for(int i=0;i<=m;i++)b[i].x=getx();
l=0;limit=1;
while(limit<=n+m)limit<<=1,l++;//传说中......作者已经死机
for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<l;i++)mcos[i]=cos(pai/(1<<i)),msin[i]=sin(pai/(1<<i));
fft(a,1);
fft(b,1);
for(int i=0;i<=limit;i++)a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<n+m;i++)printf("%d ",(int)(a[i].x/limit+0.5));//int没有自带四舍五入
printf("%d\n",(int)((a[n+m].x/limit+0.5)));
return 0;
}于是，到此，普通FFT结束了。
超~毒瘤优化。
在大佬的博客里面有三次变两次的优化。
然而，今天着重讲这个大佬的博客。
其实就是毛潇大佬的3次变1.5次优化。注：此方法掉精度会比原FFT多！
普通FFT。
3->1.5优化：
所以：优化莫乱用，父母两行泪。
开始讲优化吧：
在我们观察发现，FFT的复数部分是由没有到有到没有，于是就有大佬想了，如果开始就有了会怎样？
于是大佬把奇数项放实数，偶数项放复数。
如：\(F(x)=(a_{n}^{0}+a_{n}^{1}i)+(a_{n}^{2}+a_{n}^{3}i)x+(a_{n}^{4}+a_{n}^{5}i)x^{2}+...\)
以下摘自这个大佬
这里主要讲解了复数拆开的可优化性。
于是我们开始自推自导（跟在字母后的i为复数单位）：
我们有函数\(A(x),B(x),C(x)\)，同时n=a/2,m=b/2。（a,b为A，B的最高次，且a，b为奇数，当然，不是可以补0，也可以完成证明）
\(A(x)=\sum\limits_{i=0}^{n}(a_{2i}+a_{2i+1}i)x^{i},B(x)=\sum\limits_{i=0}^{m}(b_{2i}+b_{2i+1}i)x^{i}\)
而\(C(x)\)为目标函数。
注：如\(a_{-1}\)之类的越界数字等于0
\(C(x)=\sum\limits_{i=0}^{n}(\sum\limits_{j=0}^{m}((a_{2i}b_{2j}+a_{2i+1}b_{2j-1})+(a_{2i+1}b_{2j}+a_{2i}b_{2j+1})i)x^{i+j})+a_{2i-1}b^{2m+1}x^{i+m}\)（将目标函数奇偶分后的东西）
而：\(A(x)*B(x)=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}((a_{2i}b_{2j}-a_{2i+1}b_{2j+1})+(a_{2i+1}b_{2j}+a_{2i}b_{2j+1})i)x^{i+j}\)
我们可以相减看一下有什么方法可以把\(A(x)*B(x)\)可以变成\(C(x)\)。注：N=n+m+1
我们现在一定很好奇为什么N要多+1，考虑多项式其实要带入项数个单位根的，于是+1就是项数
\[C(x)-A(x)*B(x)=\sum\limits_{i=0}^{n}(\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j-1}+a_{2i+1}b_{2j+1})x^{i})+a_{2i-1}b_{2m+1}x^{i+m}\]
然后我们尝试将\(w_{N}^{k}(0≤k≤N-1)\)带入
\[=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j-1})w_{N}^{k(i+j)}+\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}+a_{2i-1}b_{2m+1}w_{N}^{k(i+m)}\]
\[=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j+1)}+\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}\]
\[=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j+1)}+(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}\]
\[=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}*w_{N}^{k}+(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}\]
\[=(1+w_{N}^{k})\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}(a_{2i+1}b_{2j+1})w_{N}^{k(i+j)}\]
\[=(1+w_{N}^{k})\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}a_{2i+1}*w_{N}^{ki}*b_{2j+1}*w_{N}^{kj}\]
在上文我们知道\(X_{i}-\overline{X_{N-i}}=2i\sum\limits^{N-1}_{j=0}b_{j}w_{N}^{ij}\)（\(b_{j}\)代表的是复数部，同时N其实代表的是X多项式的项数，而不是最高次，且这里的N不是上面的N），所以我们也可以得知\(A_{i}-\overline{A_{N-i}}=2i\sum\limits^{n}_{j=0}a_{2j}x_{n+1}^{ij}\)（注：\(A_{i}\)表示将\(w_{N}^{i}\)带入\(A(x)\)得到的点值）
同时，从FFT的实现细节我们有知道\(a=b=\)乘完后的最高次（补零办到的），这与N-1是同一个道理，所以：
\(A_{i}-\overline{A_{N-i}}=2i\sum\limits^{n}_{j=0}a_{2j}x_{N}^{ij}\)
所以
\[=(1+w_{N}^{k})\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}a_{2i+1}*w_{N}^{ki}*b_{2j+1}*w_{N}^{kj}\]
\[=-\frac{(1+w_{N}^{k})(A_{i}-\overline{A_{N-i}})(B_{i}-\overline{B_{N-i}})}{4}\]
于是我们便可以得出来了，但是这种方法掉精度比较严重。
同时，N其实就是我们的\(limit\)了。
放上代码：
#include<cstdio>
#include<cstring>
#include<cmath>
#define pai 3.1415926535897932
#define N 1100000
using namespace std;
inline int get()
{
int x=0,f=1;char c=getchar();
while(c>'9' || c<'0')c=='-'?f=-1:0,c=getchar();
while(c<='9' && c>='0')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*f;
}
inline int log2(int x)
{
char result=0;
x&0xffff0000?result+=16,x>>=16:0;
x&0x0000ff00?result+=8,x>>=8:0;
x&0x000000f0?result+=4,x>>=4:0;
x&0x0000000c?result+=2,x>>=2:0;
x&0x00000002?result+=1,x>>=1:0;
return result;
}
//小工具区
struct dian
{
double x,y;
dian(double xx=0,double yy=0){x=xx;y=yy;}
}a[N],b[N],c[N],wn[600000];double msin[30],mcos[30];
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator*(dian x,dian y){return dian(x.x*y.x-x.y*y.y,x.y*y.x+x.x*y.y);}
inline dian operator!(dian x){return dian(x.x,-x.y);}
inline void dwap(dian &x,dian &y){dian tt=x;x=y;y=tt;}
int r[N],l,limit,n,m;
//定义
void fft(dian A[],int type)
{
for(int i=0;i<limit;i++)
{
if(i<r[i])dwap(A[i],A[r[i]]);
}
//调整好蝴蝶迭代。
for(int id=0;id<l;id++)
{
int mid=(1<<id);
wn[1]=dian(mcos[id],msin[id]*type);for(int k=2;k<mid;k++)wn[k]=wn[(k>>1)]*wn[(k>>1)+(k&1)];
for(int L=0,R=(mid<<1);L<limit;L+=R)
{
for(int k=0;k<mid;k++)//开始算
{
dian x=A[L+k],y=A[k+mid+L]*wn[k];
A[L+k]=x+y;A[k+mid+L]=x-y;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m);
l=log2((n+m)>>1)+1;limit=(1<<l);//处理好边界
for(int i=0;i<=n;i++)
{
double x=get();
(i&1)?a[i>>1].y=x:a[i>>1].x=x;
}
for(int i=0;i<=m;i++)
{
double x=get();
(i&1)?b[i>>1].y=x:b[i>>1].x=x;
}
//输入
for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
//蝴蝶
wn[0]=dian(1,0);
for(int i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
//单位根
fft(a,1);
fft(b,1);
c[0]=a[0]*b[0]-((dian(1,0)+wn[0])*(a[0]-!a[0])*(b[0]-!b[0]))*dian(0.25,0);//wn[limit]=wn[0]
for(int i=1;i<(limit>>1);i++)
{
int j=limit-i;
c[i]=a[i]*b[i]-((a[i]-!a[j])*(b[i]-!b[j])*(dian(1,0)+wn[i]))*dian(0.25,0);//小数不能用>>2;
}
for(int i=(limit>>1);i<limit;i++)
{
int j=limit-i;
c[i]=a[i]*b[i]-((a[i]-!a[j])*(b[i]-!b[j])*(dian(1,0)-wn[i^(limit>>1)]))*dian(0.25,0);//小数不能用>>2;
}//为什么分成两个循环？因为在FFT里面我们只处理了一半的单位根。
fft(c,-1);
//处理完毕
printf("%d",int(c[0].x/limit+0.5));//防止double的精度误差
for(int i=1;i<=n+m;i++)
{
if(i&1)printf(" %d",int(c[i>>1].y/limit+0.5));
else printf(" %d",int(c[i>>1].x/limit+0.5));
}
printf("\n");
return 0;
}实战！
First
许多人就是为了这个学习的FFT吧，其实非常简单，只要乘完以后在进位就行了。
#include<cstdio>
#include<cstring>
#include<cmath>
#define N 80000
#define pai 3.1415926535897932
using namespace std;
inline int log2(int x)
{
char result=0;
x&0xffff0000?result+=16,x>>=16:0;
x&0x0000ff00?result+=8,x>>=8:0;
x&0x000000f0?result+=4,x>>=4:0;
x&0x0000000c?result+=2,x>>=2:0;
x&0x00000002?result+=1,x>>=1:0;
return result;
}
//小工具
struct dian
{
double x,y;
dian(double xx=0,double yy=0){x=xx;y=yy;}
}a[N],b[N],c[N],wn[40000];double msin[30],mcos[30];
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator*(dian x,dian y){return dian(x.x*y.x-x.y*y.y,x.y*y.x+x.x*y.y);}
inline dian operator!(dian x){return dian(x.x,-x.y);}
inline void swap(dian &x,dian &y){dian z=x;x=y;y=z;}
int n,m,l,r[N],limit,res[130000];
char st[N];
//定义区
void fft(dian *A,int type)
{
for(int i=0;i<limit;i++)
{
if(i<r[i])swap(A[i],A[r[i]]);
}
for(int id=0;id<l;id++)
{
int mid=(1<<id);
wn[1]=dian(mcos[id],msin[id]*type);for(int k=2;k<mid;k++)wn[k]=wn[(k>>1)]*wn[(k>>1)+(k&1)];
for(int L=0,R=(mid<<1);L<limit;L+=R)
{
for(int k=0;k<mid;k++)
{
dian x=A[L+k],y=A[L+mid+k]*wn[k];
A[L+k]=x+y;A[L+mid+k]=x-y;
}
}
}
}
int main()
{
// freopen("hehe.cpp","r",stdin);
// freopen("heho.txt","w",stdout);
scanf("%d",&n);m=--n;
scanf("%s",st+1);
for(int j=1;j<=n+1;j++)
{
int i=n-j+1;
(i&1)?a[i>>1].y=st[j]-'0':a[i>>1].x=st[j]-'0';
}
scanf("%s",st+1);
for(int j=1;j<=m+1;j++)
{
int i=m-j+1;
(i&1)?b[i>>1].y=st[j]-'0':b[i>>1].x=st[j]-'0';
}
//数日
l=log2(n+m>>1)+1;limit=(1<<l);//有可能n+m=0的情况下l会多，但是不怕，多就多，耗不了几ms
//FFT的2的l次方的限制
wn[0]=dian(1,0);for(int i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
//单位根
for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
//蝴蝶。
fft(a,1);
fft(b,1);
c[0]=a[0]*b[0]-(a[0]-!a[0])*(b[0]-!b[0])*(dian(1,0)+wn[0])*dian(0.25,0);
for(int i=1;i<(limit>>1);i++)
{
int j=limit-i;
c[i]=a[i]*b[i]-(a[i]-!a[j])*(b[i]-!b[j])*(dian(1,0)+wn[i])*dian(0.25,0);
}
for(int i=(limit>>1);i<limit;i++)
{
int j=limit-i;
c[i]=a[i]*b[i]-(a[i]-!a[j])*(b[i]-!b[j])*(dian(1,0)-wn[i^(limit>>1)])*dian(0.25,0);
}
fft(c,-1);
//FFT
for(int i=0;i<=n+m;i++)
{
i&1?res[i]+=int(c[i>>1].y/limit+0.5):res[i]+=int(c[i>>1].x/limit+0.5);
res[i+1]+=res[i]/10;
res[i]%=10;
}
while(res[n+m+1])
{
n++;
res[n+m+1]+=res[n+m]/10;
res[n+m]%=10;
}
while(!res[n+m])n--;
//进位
for(int i=n+m;i>=0;i--)printf("%d",res[i]);
printf("\n");
//输出
return 0;
}Second
原式等于：
\(E_{j}=\frac{F_{j}}{q_{j}}=\sum_{i<j} \frac{q_{i}}{(j-i)^{2}}-\sum_{i>j} \frac{q_{i}}{(i-j)^{2}}\)
于是我们设\(A_{i}=\frac{1}{i^{2}},B_{i}=qi\)
我们很容易看出\(E_{j}=\frac{F_{j}}{q_{j}}=\sum_{i<j} A_{i}B_{j-i}-\sum_{i>j} A_{i}B_{i-j}\)
前面的我们发现是卷积，但是后面的呢？我们也可以卷呀，我们设\(C_{i}=A_{n-i+1}\)
原式为\(E_{j}=\frac{F_{j}}{q_{j}}=\sum_{i<j} A_{i}B_{j-i}-\sum_{i>j} C_{n-i+1}B_{i-j}\)
我们可以看出，后面的那个式子也可以优化了，因为\((n-i+1)+(i-j)=n-j+1\)。
所以，我们可以处理出\(A\)与\(B\)的卷积为\(D\)、\(A\)与\(B\)的卷积为\(H\)（第0次可以手动补齐），然后\(E_{j}=D_{j}-H_{n-j+1}\)
重要提示：只能用普通的FFT
代码：
#include<cstdio>
#include<cstring>
#include<cmath>
#define N 280000
#define pai 3.1415926535897932
using namespace std;
inline int log2(int x)
{
char result=0;
x&0xffff0000?result+=16,x>>=16:0;
x&0x0000ff00?result+=8,x>>=8:0;
x&0x000000f0?result+=4,x>>=4:0;
x&0x0000000c?result+=2,x>>=2:0;
x&0x00000002?result+=1,x>>=1:0;
return result;
}
//Log运算
struct dian
{
double x,y;
dian(double xx=0,double yy=0){x=xx;y=yy;}
}x1[N]/*正数组*/,x2[N]/*倒数组*/,x3[N]/*处理1/i^2*/,wn[140000];double msin[30],mcos[30];
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator*(dian x,dian y){return dian(x.x*y.x-x.y*y.y,x.y*y.x+x.x*y.y);}
inline dian operator!(dian x){return dian(x.x,-x.y);}
inline void dwap(dian &x,dian &y){dian z=x;x=y;y=z;}
int limit,l,r[N],n;
//数组的定义
void fft(dian *A,int type)//FFT主要的函数
{
for(int i=0;i<limit;i++)
{
if(i<r[i])dwap(A[i],A[r[i]]);
}
for(int id=0;id<l;id++)
{
int mid=(1<<id);
wn[1]=dian(cos(pai/(1<<id)),sin(pai/(1<<id))*type);for(int k=2;k<mid;k++)wn[k]=wn[k-1]*wn[1];
//单位根处理完毕
for(int R=(mid<<1),L=0;L<limit;L+=R)
{
for(int k=0;k<mid;k++)
{
dian x=A[L+k],y=A[L+mid+k]*wn[k];
A[L+k]=x+y;A[L+mid+k]=x-y;
}
}
}
}
int main()
{
// freopen("testdata.in","r",stdin);
// freopen("my.out","w",stdout);
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
scanf("%lf",&x1[i].x);
x2[n-i+1].x=x1[i].x;x3[i].x=1.0/(i*1.0)/(i*1.0);
}
//处理所有的数组
l=log2(n+n/*n+n>>1*/)+1;limit=(1<<l);
//处理了2^i
wn[0]=dian(1,0);
//单位根处理
for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
//蝴蝶。
fft(x1,1);fft(x2,1);fft(x3,1);
for(int i=0;i<limit;i++)x1[i]=x1[i]*x3[i],x2[i]=x2[i]*x3[i];
fft(x1,-1);fft(x2,-1);
for(int i=0;i<=limit;i++)x1[i].x/=limit,x2[i].x/=limit;
for(int i=1;i<=n;i++)printf("%.4lf\n",x1[i].x-x2[n-i+1].x);
return 0;
}温馨插入：生成函数
FFT有一大趴的题目都是生成函数，生成函数是什么？
据网上的大佬说：生成函数分好几类。
但是作者只会一点点比较简单的。重点
我们引出一道十分简单的题目，有5个数字，分别为\(1,2,3,4,5\)（一下的都是可以重复选的）
那么，如何求选两个数字后所能得到的所有数字呢？
暴力背包来一波！
于是我们可以求出两个数字的。
但是可以不选呢？
我们可以多加一个为0的数字，也许我们看起来不重要，毕竟背包也可以无视他，但是后面会起到重要作用。
于是，有一个大佬问了，那如果我让你求每个数字的方案数呢？
还是背包呀？
漏了一句，我们把数列扩大到万级别的数列，每个数字不一定按顺序。
不是，这怎么做，这人想想都只能\(n^{2}\)吧。数学原本就很不符合常人逻辑
但是，我们设一个多项式：\(x^{0}+x^{1}+x^{2}+...\)
而这个序列的系数是什么？\(x_{i}\)的系数表示我们原本的系列中有多少为\(i\)的数，没有为则\(0\)。
这个序列有什么用？
拿\(1,2\)举例（可以不选）：\((x^{0}+x^{1}+x^{2})(x^{0}+x^{1}+x^{2})\)
自乘一下，得出：
\(x^{0}+2x^{1}+3x^{2}+2x^{3}+x^{4}\)
我们又发现，\(x^{i}\)的系数表示的就是选两次后可以组成\(i\)的方案数（这都被你发现了。），那么，以此类推：乘三次的话就是选三次的。
为什么？
我们分析一下，乘完以后是\(x^{i}\)的不就是\(x^{a}*x^{b}=x^{a+b}=x^{i}\)吗？所以这个过程就模拟了加的时候的过程，而且系数也可以对应的乘起来加起来。
至此，简单的生成函数OK了。
Third
这道题目就是生成函数了，不过还有容斥原理（顺序不同算一种，且每个物品只能选一次）。
我们先考虑一下：\((a,b,c)(a≠b≠c)\)有几种情况？容易知道有六种情况。
\((a,b,c)(a=b≠c)\)有几种情况？容易知道有三种情况。
\((a,b,c)(a=b=c)\)有几种情况？你不废话，就一种。
那么如何容斥？
我的方法：
添加一个0（不选），自乘两次，容易知道包括了选0、1、2、3次的结果，设这个多项式为\(B\)。
但是十分麻烦的是，我们要用容斥排除掉不是的。
我们先考虑减掉两次重复选同一个的情况，怎么办呢？
我们可以强行把每个数字的平方做成一个新的系列，也就是\((a_{1})^{2}+(a_{2})^{2}+...\)，这个数列表示每个数字重复选两次的情况，这次不同了，不是自乘的生成函数，而是把新的数列和原数列用生成函数乘一下，就可以得到每一个重复选同一个数两次的情况，把这个多项式乘为\(C\)，然后\(B=B-3*C\)。
现在再看重复选三个数的情况：
不对！
好像C里面也包括了重复选三个数的情况！怎么办？这样B里面重复选三个数字的部分就成了-2，就不是0了，于是，我们提前在这些位置加上2。
大功告成！
而且选两次的只不过包括了一个0。也满足一个数没有重复选的情况。（注：一个数重复选只是单纯指一个位置上的数，不包括其他和他值相同的数。）
不对！（这作者怎么老是一惊一乍。）
选1次的为什么没有呢？
因为选一次包括了两个0呀！
所以我们只要加上就行了。
真！大功告成。
代码：
#include<cstdio>
#include<cstring>
#include<cmath>
#define N 150000
#define pai 3.1415926535897932
using namespace std;
inline int log2(int x)
{
char result=0;
x&0xffff0000?(x>>=16),result+=16:0;
x&0x0000ff00?(x>>=8),result+=8:0;
x&0x000000f0?(x>>=4),result+=4:0;
x&0x0000000c?(x>>=2),result+=2:0;
x&0x00000002?(x>>=1),result+=1:0;
return result;
}
inline int getx()
{
int x=0,f=1;char c=getchar();
while(c>'9' || c<'0')c=='-'?f=-1:0,c=getchar();
while(c<='9' && c>='0')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*f;
}
inline int mymax(int x,int y){return x>y?x:y;}
int r[N],limit,l,n,m,res[N],a[N];
struct dian
{
double x,y;
dian(double xx=0,double yy=0){x=xx;y=yy;}
}aa[N],bb[N],cc[N],dd[N],wn[80000];double msin[30],mcos[30];
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator*(dian x,dian y){return dian(x.x*y.x-x.y*y.y,x.y*y.x+x.x*y.y);}
inline dian operator!(dian x){return dian(x.x,-x.y);}
inline void dwap(dian &x,dian &y){dian tt=x;x=y;y=tt;}
inline int red(dian *A,int id){return int(((id&1)?(A[id>>1].y):(A[id>>1].x))/limit+0.5);}
void fft(dian *A,int type)
{
for(int i=0;i<limit;i++)
{
if(i<r[i])dwap(A[i],A[r[i]]);
}
for(int id=0;id<l;id++)
{
int mid=(1<<id);
wn[1]=dian(mcos[id],msin[id]*type);for(int k=2;k<mid;k++)wn[k]=wn[k>>1]*wn[(k>>1)+(k&1)];
for(int R=(mid<<1),L=0;L<limit;L+=R)
{
for(int k=0;k<mid;k++)
{
dian x=A[L+k],y=A[L+mid+k]*wn[k];
A[L+k]=x+y;A[L+mid+k]=x-y;
}
}
}
}
inline void work(dian *A,dian *B,dian *C)
{
C[0]=A[0]*B[0]-(dian(1,0)+wn[0])*(A[0]-!A[0])*(B[0]-!B[0])*dian(0.25,0);
for(int i=1;i<(limit>>1);i++)
{
int j=limit-i;
C[i]=A[i]*B[i]-(dian(1,0)+wn[i])*(A[i]-!A[j])*(B[i]-!B[j])*dian(0.25,0);
}
for(int i=(limit>>1);i<limit;i++)
{
int j=limit-i;
C[i]=A[i]*B[i]-(dian(1,0)-wn[i^(limit>>1)])*(A[i]-!A[j])*(B[i]-!B[j])*dian(0.25,0);
}
}
int main()
{
n=getx();
aa[0].x++;
for(int i=1;i<=n;i++)
{
a[i]=getx();res[a[i]*3]+=2;/*容斥*/
m=mymax(a[i],m);
(a[i]&1)?aa[a[i]>>1].y++:aa[a[i]>>1].x++;
}
m=(m<<1)+m;
l=log2(m>>1)+1;limit=(1<<l);for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
wn[0]=dian{1,0};msin[1]=1;mcos[1]=0;for(int i=2;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
fft(aa,1);
work(aa,aa,bb);
work(aa,bb,cc);
fft(cc,-1);
memset(bb,0,sizeof(bb));bb[0].x++;//重复利用b数组
for(int i=1;i<=n;i++)bb[a[i]].x++;
fft(bb,1);
work(aa,bb,dd);
fft(dd,-1);
for(int i=0;i<=m;i++)res[i]+=red(cc,i)-red(dd,i)*3,res[i]/=6;
for(int i=1;i<=n;i++)res[a[i]]++;
for(int i=0;i<=m;i++)
{
if(res[i])printf("%d %d\n",i,res[i]);
}
return 0;
}未完待续。。。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
