快速傅里叶变换（FFT）学习笔记（未完待续）

最新推荐文章于 2024-05-16 10:42:21 发布

zhangjianjunab

最新推荐文章于 2024-05-16 10:42:21 发布

阅读量682

点赞数

本文链接：https://blog.csdn.net/zhangjianjunab/article/details/88592647

版权

总所周知，FFT是一个非常麻烦的算法，再加上博主语文不好，便写起来有点麻烦，但会尽力去写。 ~~要以后自己看不懂就。。。~~

注：因为最近的压力紧张，便没有继续学习FFT，这仅为目前的半成品以及一些目前已知的优化。

参考资料

讲解1

讲解2

优化

毛啸大佬的论文（不放上下载地址，因为下载论文的主要OJ太卡）

优化2

优化证明

题解

因为此篇文章比较依赖这些文章，为此先放出来。

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 (n l o g n)$ 的时间，而暴力要 $O(n^{2})$ 。

例题

普通做法

上联：打表过样例，暴力出奇迹。

下联：花中一片叶，luogu两行泪。

。。。

更高大尚的做法

于是，就有人开始尝试新的做法了！

尝试方法

说到底，FFT究竟是什么？

大概样子：

系数表达式->点值表达式->点值相乘->系数表达式

现在，就为大家一一刨析做法！

定义与一部分性质

没学过多项式的估计你还没学OI呢，~~神童？~~默认大家会一点多项式的定义。

咱们先理清我说的一些十分玄学的话。

系数表达式

设 $A (x)$ 表示一个n?1次多项式

则 $\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 + b i) + (c + d i) = (a + c) + (b + d) i$
复数相减： $(a + b i) - (c + d i) = (a - c) + (b - d) i$
复数相乘： $a+bi)+(c+di)=ac+adi+bci+dbi^{2}=(ac-bd)+(ad+bc)i$
复数的共轭： $a + b i = a - b i$

复数相乘的二维意义就是模长相乘，幅角相加。

证明来源于https://www.luogu.org/blog/command-block/fft-xue-xi-bi-ji
在这里插入图片描述

单位根

其实，FFT就是DFT加上IDFT

DFT是什么

DFT就是一个能在 $O (n l o g n)$ 内将系数表达式转成点值表达式

我们平常把系数表达式转成点值表达式不就是一个个代，然后就 $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

想象一个复数 x 。

如果 $∣ x ∣ > 1$ ，即 $x^{n}|=|x|^{n}>1$ (模长乘模长，越乘越大)

如果 $∣ x ∣ < 1$ ，即 $x^{n}|=|x|^{n}<1$ (模长乘模长，越乘越小)

所以只有模长等于1的复数才有可能成为n次单位根。

也就是说，只有单位圆上的点表示的复数才有可能成为n次单位根
我们成功缩小了探究范围，只在这一个圆上研究。

(下面提及的复数,均是在单位圆上的复数!!!)
想象一个模长为1的复数为 $X$ ，设它的幅角为(圆周 $* a$ )

(这里 $0 < = a < 1$ ,不然会产生重复)

它的n次方的幅角为 $n a ?$ 圆周

那么，如果它是n次单位根，幅角的n倍 ( $n a ?$ 圆周 ) 一定是圆周的自然数倍。

得到不定方程(听到这个名字别害怕)：

( $n a ?$ 圆周 ) = 圆周 $? k$ (k 为自然数)

等式两边同除以“圆周”

得 $n a = k$ (k 为自然数参数)

分类讨论一下：

$k = = 0$

由于 $n > 0$ (方程的次数是正整数)只有 $a = 0$ 这一种情况

$k > 0$

$k$ 可以为 $0, 1, 2, 3 . . .$

$k$ 为1时,得 $n a = 1, a = 1 / n$

$k$ 为2时,得 $n a = 2, a = 2 / n$

…

$n a = k, a = k / n$

但是， $k > 0$ 时并没有无数个单位根，他们会重叠，重叠的只算一个。

如 $k = n$ 时,得 $n a = 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次方的幅角为 $3 a ?$ 圆周

那么，如果它是3次单位根，幅角的3倍 ( $3 a ?$ 圆周 ) 一定是圆周的自然数倍。

所以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 (n l o g n)$ 的。

？？？

别急，我们一步步来，首先，我们把一个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})$ ，我们再设两个函数： $F L (x), F R (x)$ 。

有： $F(x)=FL(x^{2})+x*FR(x^{2})$ 。

$F L (x), F R (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 \geq 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 (n l o g n)$ 的到了点值表达式

IDFT

IDFT就是一个可以将点值表达式转成系数表达式的东西，时间复杂度 $O (n l o g n)$

如何转？拉格朗日插值法~~其实我不会拉格朗日插值法~~？高斯消元？

貌似都可以，但是都很慢。~~争做世界上最快的男人~~

于是，又有大佬想出了方法，因为证明过程太过牛逼~~作者太懒~~，所以放上一名大佬的证明

在这里插入图片描述

因此，DFT与IDFT大家都会了，就到了高潮重点了。

蝴蝶迭代优化

因为FFT的分治过程是用DFS，而DFS会严重拖慢速度，于是有大佬想到了蝴蝶迭代，使得 FFT拜托了DFS的“控制”

因为~~作者懒~~作者语文不好，为了更好的食用，放上大佬的博客摘录。

在这里插入图片描述

这个也挺好证的：

0 1 2 3|4 5 6 7 不难发现，左边的一层最高位为0，右边的一层最高位为1
0 2 4 6|1 3 5 7 不难发现，左边的一层最低位为0，右边的一层最低位为1

然后不断分下去皆是如此。

单位根求法

不是，单位根怎么求呀，在此之前，请大家记一下三角函数。

尤其是正弦函数(cos)与余弦函数(sin)，在C++的cmath库内有函数。

在这里插入图片描述
B点便是我们要求的了，那么他的坐标就是 $(A C, B C)$ 怎么求？我们知道AB长为1，同时

于是，大部分我们都准备好了，可以开始实现了。 $∠ A C B = 90 °, ∠ B A C$ 也可以求出来，那么利用三角函数易得知 $C B = s i n (∠ B A C) * A B, A C = c o s (∠ B A C) * A B$ ，同时 $A B = 1$ ，所以 $C B = s i n (∠ B A C), A C = c o s (∠ B A C)$ 。

但是C++里面的三角函数采用弧度制下一个整圆不再是 $360$ °,而是 $2 π$ 了。

所以 $∠ B A C =$ $2π\over n$ ，而不是 $∠ B A C =$ $360°\over n$ 了

当然，万事皆有烦人之事， $w_{n}^{1}$ 在 $(0, 1)$ ，自己画图就知道三角函数这时候就很惆怅了，有两种方法：1.霸王硬上弓，一般不影响结果。2.特判，更稳，但要记。

实现、细节与小优化

就到了我们的实现细节了。

细节

π

$π$ 等于多少？

当然，你可以去背，再这打上：3.1415926535897932

也可以用acos(-1)，精度很好，不过最后一位错了。

反正我是背的，当然，听说用acos(-1)更好。

补齐？

我们之前也讲了，DFT的分治基于的是 $n$ （ $n$ 次多项式）能写成 $2^{i}-1$ 的形式，也就是多项式中有 $2^{i}$ 项，但是题目中并不能保证，我们可以手动补齐，计算离 $2 n + 1$ （因为结果的最大次为 $2 n$ ）最近且比他大的2的幂次。

如： $n = 5$ 因补齐到15。

合并

通过看大佬的博客，相信大家也会有所了解。

IDFT的代码和DFT的代码可以合并，但是需要注意的是IDFT中的单位根都是 $w_{n}^{-1}$ 什么的，于是我们的 $y$ 坐标要标成负数，当然，我们只需要多传一个 $t y p e$ 值就行了。（注意，IDFT的结果要除以多项式的项数 $(l i m i t)$ ）

预处理蝴蝶迭代？

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，怕被卡就用NTT吧（来自作者第二次更新时的提醒，当然作者自己也决定不复习这个优化，等到以后什么时候有用的时候再复习吧）。

这题

普通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其实就是我们的 $l i m i t$ 了。

不过其实我十分怀疑这个优化适用于系数为复数的情况，因为把复数带进去貌似也是可以证明的，以后试试吧。

而且这个优化应该满足可加性质。

放上代码：

#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;
}

论如何将两次DFT变成一次

首先，我们要求出 $A(x)=a_0+a_1x+a_2x^2+a_3x^3+...,a_nx^nB(x)=b_0+b_1x+b_2x^2+b_3x^3+...b_nx^n$
处理这两个的DFT，我们想到了我们可以把 $B (x)$ 放在虚数部，那么就是 $P (x) = A (x) + B (x) i$ ，我们设函数 $A (x)$ 代入 $w_{n}^{k}$ 后的式子为 $A^{'} (k)$ （其他也是一样）。

那么很明显 $P^{'} (x) = A^{'} (x) + B^{'} (x) i$ ，因为我们拆开之后： $P'(x)=\sum\limits_{j=0}^{n}(a_{j}+b_{j}i)w_{n}^{xj}=\sum\limits_{j=0}^{n}a_{j}w_{n}^{xj}+b_{j}w_{n}^{xj}i=A'(x)+B'(x)$

然后我们设 $Q (x) = A (x) - B (x)$ ，那么 $Q^{'} (x) = A^{'} (x) - B^{'} (x)$ 。（注意：以上证明系数可以为复数）

那么我们只要知道 $P^{'}, Q^{'}$ ，那么就可以用公式： $A'(x)=\frac{Q'(x)+P'(x)}{2},B'(x)=\frac{P'(x)-Q'(x)}{2i}$ ，算出来了。

但是关键是如何用 $P^{'}$ 算 $Q^{'}$ 呢？

看证明！！！来自毛老爷的证明（温馨提示： $c o s (- x) = c o s (x)$ ）
在这里插入图片描述

注意，这里的 $2 L - k$ 是在 $\%2L$ 的意义下， $0$ 还是有点特殊的。

但是如果我们有了 $A^{'}, B^{'}$ ，怎么一次转回去呢？自然就是 $P^{'} (x) = A^{'} (x) + B^{'} (x) i$ ，然后跑一下IDFT，然后实数就是 $A^{'}$ 的IDFT，虚数就是 $B^{'}$ 的IDFT。

等会！ $A^{'}, B^{'}$ 是复数域的，那么是不是代表复数也可以用此优化？很可惜，不能用，即使能用，也不是这个证明。

因为我们可以假设 $A, B$ 就是复数域的。

那么 $A(x)=\sum\limits_{j=0}^{n}(a_{j}+c_ji)x^j,B(x)=\sum\limits_{j=0}^{n}(b_{j}+d_ji)x^j$ 。

那么 $P(x)=\sum\limits_{j=0}^{n}((a_{j}-d_{j})+(b_{j}+c_{j})i)x^{j}$ ，也就是归到原本应有的部分。

而 $Q(x)=\sum\limits_{j=0}^{n}((a_{j}+d_{j})+(-b_{j}+c_{j})i)x^{j}$

那么化出来就已经不是一样的了。

我们化完之后应该是 $Q'(k)=conj(\sum\limits_{j=0}^{n}((a_{j}+d_{j})+(b_{j}-c_{j})i)w_{n}^{kj})=conj(\sum\limits_{j=0}^{n}(a_{j}-c_{j}i)w_{n}^{kj}+(b_{j}-d_{j}i)w_{n}^{kj})$ ，结果我们求到的里面也是共轭。

但是为什么可以 $I D F T$ 求复数域的。

因为我们是利用了 $I D F T$ 的本质，把一个 $D F T$ 过的转回去，当然这里也就体现了 $D F T$ 支持复数系数。

而且 $I D F T$ 求完后不是用公式而是直接看实数虚数域的。

这个优化在MTT（或者叫拆系数FFT）中优化比较明显，而且代码量也比较可观。

利用这个优化可以实现MTT中的 $7$ 次变 $4$ 次优化。

论如何将一次DFT变成0.5次

这个其实还是比较简单的，如果你会了上面那个的话。

对于 $f(w_{n}^{k})=f_0(w_{n/2}^{k})+w_{n}^{k}f_1(w_{n/2}^{k}),f(w_{n}^{k+n/2})=f_0(w_{n/2}^{k})-w_{n}^{k}f_1(w_{n/2}^{k})$ ，那么我们是不是可以算出 $f_0,f_1$ ，然后推出 $f$ 呢？

而且长度减半。

那么对于 $I D F T$ 呢？

我们逆推回去，根据上面两条，我们在知道了 $f$ 之后，可以逆推出 $f_0,f_1$ ，然后求 $f_0,f_1$ 的IDFT，用上述方法合并就行了，不过对于 $w_{n}^k$ ，我们需要求出 $w_n^{n-k}$ 然后乘起来。

我的代码比较丑，不要嫌弃，因为只是随手打的，没有美不美观之说。

其实这个也可以实现3次变1.5次优化，这道题目就是模版题。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define  N  210000
#define  NN  310000
using  namespace  std;
const  double  pai=acos(-1.0);
int  limit,l,r[N];
struct  node
{
    double  x,y;
    node(double  xx=0,double  yy=0){x=xx;y=yy;}
}a[N],b[N],c[N],aa[NN],bb[NN],wn[NN];double  msin[30],mcos[30];
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  node  operator!(node  x){return  node(x.x,-x.y);}
inline  node  swap(node  &x,node  &y){node  z=x;x=y;y=z;}
void  fft_prepare(int  len)
{
    limit=1;l=0;
    while(limit<=len)limit<<=1,l++;
    limit>>=1;l--;
    for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int  i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
}
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  L=0,R=mid<<1;L<limit;L+=R)
        {
            for(int  k=0;k<mid;k++)
            {
                node  x=A[L+k],y=wn[k]*A[L+k+mid];
                A[L+k]=x+y;A[L+k+mid]=x-y;
            }
        }
    }
}
void  chuli(node  x[],node  y[])
{
    for(int  i=0;i<limit;i++)//求出f
    {
        int  j=(limit-1)&(limit-i);
        node  xx=(x[i]+!x[j])*node(0.5,0),yy=(x[i]-!x[j])*node(0,-0.5)*wn[i];
        y[i]=xx+yy;y[i+limit]=xx-yy;
    }
}
int  n,m,nn;
inline  int  mymax(int  x,int  y){return  x>y?x:y;}
int  main()
{
    wn[0]=node(1,0);
    scanf("%d%d",&n,&m);nn=n+m;
    fft_prepare(nn);
    for(int  i=0;i<=n;i++)
    {
        if(i&1)scanf("%lf",&a[i/2].y);
        else  scanf("%lf",&a[i/2].x);
    }
    for(int  i=0;i<=m;i++)
    {
        if(i&1)scanf("%lf",&b[i/2].y);
        else  scanf("%lf",&b[i/2].x);
    }
    fft(a,1);fft(b,1);
    wn[1]=node(cos(pai/(1<<l)),sin(pai/(1<<l)));
    int  ed=limit<<1;
    for(int  k=2;k<ed;k++)wn[k]=(wn[k>>1]*wn[(k>>1)+(k&1)]);//处理单位根
    chuli(a,aa);chuli(b,bb);
    for(int  i=0;i<ed;i++)bb[i]=bb[i]*aa[i];
    for(int  i=0;i<limit;i++)//利用f推回f0,f1
    {
        node  x=bb[i],y=bb[i+limit];
        bb[i]=(x+y)*node(0.5,0);bb[i+limit]=(x-y)*node(0.5,0)*wn[((limit<<1)-i)&((limit<<1)-1)];
    }
    for(int  i=0;i<limit;i++)c[i]=bb[i]+bb[i+limit]*node(0,1);
    fft(c,-1);
    for(int  i=0;i<=n+m;i++)
    {
        if(i&1)printf("%d ",(int)(c[i/2].y/limit+0.5));
        else  printf("%d ",(int)(c[i/2].x/limit+0.5));
    }
    return  0;
}

FFT的更多小优化

想不到吧，还有优化。

我老是忘记区分DFT和IDFT

有些人在 $F F T$ 的时候老是忘记打 $* t y p e$ 。

对于 $n - 1$ 次多项式而言

但是又想想看， $w_{n}^{-k}=w_{n}^{n-k}$ ，所以我们是不是只要对原数列做些操作，是不是就可以了？

我们发现 $w_{n}^{-kj}$ 是对于第 $j$ 项的，而如果不做任何操作 $w_n^{kj}=w_{n}^{k(j-n)}$ 会乘到 $j$ 项，而正确应该乘到 $j - n$ 项才对，所以我们需要对于 $0, n$ 翻转一下，不对哦， $0, n$ 是一个位置哦，所以应该是 $1, n - 1$ 翻转，其实就是第 $n / 2$ 位不动，然后整体绕 $n / 2$ 翻转。

而这个，你又不用特地打个IDFT。~~不是写的更多了吗？~~

精度和时间两开花，我的单位根太SB了吧

相信许多人都有这样的情况，就是如果每次DFT三角函数爆算单位根，时间开花，如果打倍增乘法优化，精度开花，那么有没有两全其美的方法呢？

有一种牺牲一点点时间的方法就是换成 $l o n g$ $d o u b l e$ ，那么大部分卡精题目都可以过。

但是其实我们可以预处理单位根，把所有单位根处理出来，而且对于每个长度的单位根所需要的个数都是 $m i d$ 的，所以仔细想想就可以发现对于长度 $l i m i t$ ，长度 $l i m i t$ 的DFT所需要的单位根个数为 $l i m i t$ 个。

当然，有的神人甚至做到了 $l i m i t / 2$ 个，为什么他的这么少，因为我们的单位根次数是翻倍的，即一次，二次，四次，八次单位根等等，那么只要处理出 $l$ 次单位根，对于第 $i$ 个 $l - 1$ 次单位根，就是 $i * 2$ 的 $l$ 次单位根。

当然预处理单位根一般和 $l o n g$ $d o u b l e$ 时间差不多，但是如果 $D F T$ 次数上来了，那么预处理就是质的飞跃了。

MTT（拆系数FFT）

例题

有时候对于模数 $1 e 9$ 级别，而且项数在 $1 e 5$ 级别，那么乘完后的级别就到了 $1 e 23$ 级别了，连 $l o n g$ $l o n g$ 都死了的那种！！！~~不是int128就可以了吗~~不给用！！！！

那么我们设一个常数 $k$ ，然后对于 $A (x)$ 系数拆成 $B (x), C (x)$ ，满足 $A(x)=k*B(x)+C(x)(c_{j}<k)$ ，其实也就是对于 $A (x)$ 的系数拆成 $a_{j}/k,a_{j}\%k$ 两个部分，那么对于 $A (x) D (x)$ ，就可以拆成 $B(x)E(x)*k^2+(B(x)F(x)+C(x)E(x))*k+C(x)F(x)$ ，那么我们就可以对于 $B, C, E, F$ 求 $D F T$ ，然后对于那四个乘积求 $I D F T$ ，八次DFT。

但是其实中间乘以 $k$ 的可以合并。

那么就是七次，考虑 $3$ 次变 $1.5$ 次优化，很明显我们可以相对的 $2$ 次DFT求出所有的DFT，并且两两相乘，实现 $4$ 次DFT，当然如果这个优化支持可加性质的话，可以实现 $3.5$ 次。

而如果使用 $2$ 变 $1$ 优化的话，由于 $7$ 会剩一个，所以就是 $4$ 次，也就是 $8$ 个两两合并的。（其实对其中 $6$ 个用这个优化， $1$ 个用下个优化也行）。

如果用 $1$ 次变 $0.5$ 次，就是 $3.5$ 次了。

这里给出的是最纯朴的拆系数FFT，也就是 $7$ 次。

另外对于 $k$ 的选择，选择 $\sqrt{mod}$ 是最好的，比较平均。

//这道题目用三角函数预处理单位根时间会快不少，大概会变成现在的二分之一
#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  310000
using  namespace  std;
typedef  long  long  LL;
LL  mod=1e9+7,fuck=32768/*文明*/;
long  double  pai=3.1415926535897932;
int  limit=1,l,r[N],pre;
struct  node
{
	long  double  x,y;
	node(long  double  xx=0,long  double  yy=0){x=xx;y=yy;}
}a1[N],a2[N],b1[N],b2[N],lin[N],wn[N];long  double  msin[30],mcos[30];
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  L=0,R=mid<<1;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  n,m;
inline  void  chuli(node  A[],LL  k)
{
//	k%=mod;
	for(int  i=0;i<limit;i++)A[i].x=(((LL(A[i].x/limit+0.5))%mod)*k)%mod;
}
int  main()
{
//	freopen("1.in","r",stdin);
//	freopen("1.out","w",stdout);
	wn[0]=node(1,0);
	scanf("%d%d%lld",&n,&m,&mod);
	for(int  i=0;i<=n;i++){int  x;scanf("%d",&x);x%=mod;a1[i].x=x/fuck;a2[i].x=x%fuck;}
	for(int  i=0;i<=m;i++){int  x;scanf("%d",&x);x%=mod;b1[i].x=x/fuck;b2[i].x=x%fuck;}
	while(limit<=n+m)limit<<=1,l++;
	for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int  i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
	fft(a1,1);fft(a2,1);fft(b1,1);fft(b2,1);
	for(int  i=0;i<limit;i++)lin[i]=a1[i]*b1[i],a1[i]=a1[i]*b2[i]+a2[i]*b1[i],a2[i]=a2[i]*b2[i];
	fft(lin,-1);fft(a1,-1);fft(a2,-1);
	chuli(lin,fuck*fuck);
	chuli(a1,fuck);
	chuli(a2,1);
	for(int  i=0;i<n+m;i++)printf("%lld ",((LL)(lin[i].x)+(LL)(a1[i].x)+(LL)a2[i].x)%mod);
	printf("%lld\n",((LL)(lin[n+m].x)+(LL)(a1[n+m].x)+(LL)a2[n+m].x)%mod);
	return  0;
}

生成函数，FFT的得力干员

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了。

实战！

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;
}

Third

【题意】
n（n<=40000）个物品，可以用1/2/3个不同的物品组成不同的价值，求每种价值有多少种方案(顺序不同算一种)
【输入】
第一行是整数N，表示有N个物品
接下来n行升序输入N个数字Ai(Ai≤40000)，表示每个物品的价值。
【输出】
若干行，按升序对于所有可能的总价值输出一行x y，x为价值，y为方案数。
【样例输入】
4 

4 

5 

6 

8 

【样例输出】
4 1
5 1
6 1
8 1
9 1
10 1
11 1
12 1
13 1
14 1
15 1
17 1
18 1
19 1

这道题目就是生成函数了，不过还有容斥原理（顺序不同算一种，且每个物品只能选一次）。

我们先考虑一下： $(a, b, c) (a \neq = b \neq = c)$ 有几种情况？容易知道有六种情况。
$(a, b, c) (a = b \neq = c)$ 有几种情况？容易知道有三种情况。
$(a, b, c) (a = b = c)$ 有几种情况？你不废话，就一种。

那么如何容斥？

我的方法仔细讲讲：

对于选两种的方案，其实就是选了一种为 $0$ 的，而选了一种的后面可以直接加上去。

那么我们就只要知道序列加上一个 $0$ 后选三个的方案数。

然后我们就对生成函数乘个三次方。

但是我们又发现了一个事情，就是对于 $(a, b, c)$ 方案，可以有 $6$ 种组合方法，所以要除以 $6$ ，当然还是存在一些特判的情况，所以我们先不除 $6$ 。

就是当选了两种相同的时候，我们需要把这个列成两个生成函数相乘，一个生成函数就是对于值为 $y$ 的数字， $x^y$ 的系数 $+ 1$ ，对于另外一个，就是 $x^{2y}$ 的系数 $+ 1$ 。

然后一乘就知道了，不过这个要乘以 $3$ 再减过去。

然后就是对于 $3$ 个同样的情况，那么第一次FFT这种情况每种情况只统计了一次，而在第二次 $F F T$ 也出现了一次，但是乘了 $3$ ，所以就是 $- 2$ 次，我们最后 $+ 2$ 就行了。

然后把值除以 $6$ ，并且别忘了 $+ 1$ ，就是一个的情况。

不过需要注意的是这里的相同不是指权值相同，是指位置。

代码：

#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;//最后除以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;
}

Fourth

【题意】
请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ，并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。
【输入】
第一行一个整数N,接下来N行，第i+2..i+N-1行，每行两个数，依次表示a[i],b[i] (0 < = i < N)。
【输出】
输出N行，每行一个整数，第i行输出C[i-1]。
【样例输入】
5
3 1
2 4
1 1
2 4
1 4
【样例输出】
24
12
10
6
1

还能说什么，SB题。

设 $d_{i}=a_{n-i-1}$ ，然后一看，卷积！！！

看出过程：

$c_{i}=\sum\limits_{j=i}^{n-1}d_{n-j-1}b_{j-i}$

把项数一加 $(n - 1) - i$ ，所以最后逆着输出即可。

              


?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
 
#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  310000
using  namespace  std;
const  double  pai=acos(-1.0);//数据掉精逼我也掉精。。。
int  limit=1,l,r[N],n,a1[N],a2[N];
struct  node
{
    double  x,y;
    node(double  xx=0,double  yy=0){x=xx;y=yy;}
}a[N],b[N],c[N],wn[N];double  msin[30],mcos[30];
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[1];
        for(int  L=0,R=mid<<1;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);
    scanf("%d",&n);
    for(int  i=1;i<=n;i++)scanf("%d%d",&a1[i],&a2[i]);
    while(limit<n+n)limit<<=1,l++;
    for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int  i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
    for(int  i=1;i<=n;i++)a[i].x=a1[n-i+1],b[i-1].x=a2[i];
    fft(a,1);
    fft(b,1);
    for(int  i=0;i<limit;i++)c[i]=a[i]*b[i];
    fft(c,-1);
    for(int  i=n;i>=1;i--)printf("%d\n",int(c[i].x/limit+0.5));
    return  0;
}

Fifth

【题意】
有T组数据
每组数据给出n条线段，问任意取三条，可以组成三角形的概率 

【输入】
开头一行输入T（T<=10）
下来T组数据，每组数据第一行输入一个n（3<=n<=10^5）
第二行输入n个数，表示n条线段
线段长度（1<=n<=10^5） 

【输出】
每组数据输出一个数p
表示可以组成三角形的概率
保留七位小数 

【样例输入】
2 
4 
1 3 3 4 
4 
2 3 3 4 

【样例输出】
0.5000000
1.0000000

初三再学FFT做这道题目感觉了如指掌。

太水了。

很明显对于数对 $a, b, c (a > = b > = c)$ ，如果 $b + c < = a$ 就是不成立的，所以我们只要统计这个数对的数量就行了。

那么我们先用生成函数处理出 $b + c$ 的方案数。（别忘了一个数字选两个的情况要减）

然后做一遍前缀和，对于每个数字都去看有多少种方案小于等于它，别忘了除 $2$ ，而我们也可以利用 $n * (n - 1) * (n - 2) / 6$ 算出总方案数。

然后就可以算期望了。

#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  1100000
#define  pai  3.1415926535897932
using  namespace  std;
typedef  long  long  LL;
int  limit,l,r[N];
struct  node
{
    double  x,y;
    node(double  xx=0,double  yy=0){x=xx;y=yy;}
}a[N],b[N],wn[N];double  msin[30],mcos[30];
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  L=0,R=mid<<1;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  n,pp[N];
LL  ans[N];
inline  int  mymax(int  x,int  y){return  x>y?x:y;}
int  main()
{
    wn[0]=node(1,0);
    for(int  i=0;i<30;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
    int  T;scanf("%d",&T);
    while(T--)
    {
        int  maxnum=0;
        memset(a,0,sizeof(a));
        scanf("%d",&n);
        for(int  i=1;i<=n;i++)
        {
            scanf("%d",&pp[i]);maxnum=mymax(maxnum,pp[i]);
            a[pp[i]].x+=1.0;
        }
        limit=1;l=0;
        while(limit<=maxnum+maxnum)limit<<=1,l++;
        for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        fft(a,1);
        for(int  i=0;i<limit;i++)b[i]=a[i]*a[i];
        fft(b,-1);
        for(int  i=0;i<limit;i++)b[i].x=b[i].x/limit;
        for(int  i=1;i<=n;i++)b[pp[i]+pp[i]].x-=1.0;//即这个数字选了两次的情况
        maxnum<<=1;
        double  gailv=(LL)n*(n-1)*(n-2)/6;
        for(int  i=1;i<=maxnum;i++)ans[i]=ans[i-1]+(LL)(b[i].x+0.5);
        double  cnt=0;
        for(int  i=1;i<=n;i++)cnt+=ans[pp[i]];//统计小于等于他的情况。
        printf("%.7lf\n",(gailv-cnt/2)/gailv);
    }
    return  0;
}

Sixth

【题意】
给出一串数a[i]，定义数组b经过累加变换a得到。
如a｛1,2,3,4｝
则第一个b是｛1,3,6,10｝
此后的b则为累加变换b得到
如b｛1,3,6,10｝
则第二个b｛1,4,10,20｝
我们想知道第k次变换后的b数组
【输入】
第1行，2个数N和K，中间用空格分隔，N表示数组的长度，K表示处理的次数（2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9)
【输出】
共N行，每行一个数，对应经过K次处理后的结果。输出结果 Mod 10^9 + 7。
【样例输入】
4 1
1
2
3
4
【样例输出】
1
3
6
10

一开始我们一看，不就是构造一个 $1, 1, 1, 1, 1 . . .$ 的多项式，然后乘 $k$ 次，套快速幂不就行了。

只不过每次快速幂需要只保留前面的一些数字，所以每次都需要IDFT，所以时间复杂度 $n l o g n l o g k$ ，结果模数爆炸，又改成MTT，直接T成了屎。

上网看了 $R o s e$ 巨佬的博客：

在这里插入图片描述
下面列的是 $a_0$ 的系数方阵（ $a_1$ 的就是整体右移一位）。

发现其实就是斜着的杨辉三角形，对于杨辉三角形的第 $i$ 行 $j$ 项就是 $C_(i-1)^{j-1}$ 。

不难发现如果对于最后一行的第一个数字是杨辉三角形中的第 $i$ 行第一个，那么第二个就是第 $i + 1$ 行第二个，第三个就是 $i + 2$ 行第三个。

那么不难发现设 $c_i$ 为最后一行第 $i + 1$ 个的系数，那么 $c_0=1,c_i=\frac{c_{i-1}(i+k-1)}{i}$ ，而对于除 $i$ 求逆元就行了。

然后我们就发现： $a_{i}=a_{j}*c_{i-j}$ ，卷积即可。

当然使用MTT啦。

时间复杂度： $O (n l o g n)$

#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  210000
using  namespace  std;
typedef  long  long  LL;
LL  mod=1e9+7,c[N],fuck=32768/*文明*/;
const  long  double  pai=3.1415926535897932;
int  limit=1,l,r[N],pre;
struct  node
{
    long  double  x,y;
    node(long  double  xx=0,long  double  yy=0){x=xx;y=yy;}
}a1[N],a2[N],b1[N],b2[N],lin[N],wn[N];long  double  msin[30],mcos[30];
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  L=0,R=mid<<1;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;
            }
        }
    }
}
LL  ksm(LL  x,LL  k)
{
    LL  ans=1;
    while(k>1)
    {
        if(k&1)ans=(ans*x)%mod;
        x=(x*x)%mod;
        k>>=1;
    }
    ans=(ans*x)%mod;
    return  ans;
}
int  n,m;
inline  void  chuli(node  A[],LL  k)
{
//  k%=mod;
    for(int  i=0;i<limit;i++)A[i].x=(((LL(A[i].x/limit+0.5))%mod)*k)%mod;
}
int  main()
{
//  freopen("1.in","r",stdin);
//  freopen("1.out","w",stdout);
    wn[0]=node(1,0);
    scanf("%d%d",&n,&m);
    for(int  i=1;i<=n;i++){int  x;scanf("%d",&x);x%=mod;a1[i-1].x=x/fuck;a2[i-1].x=x%fuck;}
    while(limit<=n+n-2)limit<<=1,l++;
     
    for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int  i=0;i<l;i++)msin[i]=sin(pai/(1<<i)),mcos[i]=cos(pai/(1<<i));
    if(m==0)for(int  i=0;i<n;i++)printf("%d\n",int(a1[i].x*fuck+a2[i].x)%mod);
    else
    {
        c[0]=1;for(int  i=1;i<n;i++){c[i]=(((c[i-1]*(LL)(i+m-1)%mod)%mod)*ksm(i,mod-2))%mod;}
        for(int  i=0;i<n;i++)b1[i].x=c[i]/fuck,b2[i].x=c[i]%fuck;
        fft(a1,1);fft(a2,1);fft(b1,1);fft(b2,1);
        for(int  i=0;i<limit;i++)lin[i]=a1[i]*b1[i],a1[i]=a1[i]*b2[i]+a2[i]*b1[i],a2[i]=a2[i]*b2[i];
        fft(lin,-1);fft(a1,-1);fft(a2,-1);
        chuli(lin,fuck*fuck);
        chuli(a1,fuck);
        chuli(a2,1);
        for(int  i=0;i<n;i++)printf("%lld\n",((LL)(lin[i].x)+(LL)(a1[i].x)+(LL)a2[i].x)%mod);
    }
    return  0;
}

n阶等差也可以这么搞。

Last

【题意】
给定两个字符串a（长度<=10^5）和b（长度<=a），求b在a中出现了几次以及分别在哪些位置出现。
b中会存在“？”字符，这个字符可以匹配所有字母
【输入】
两行两个字符串，分别代表a和b
【输出】
第一行一个正整数m，表示b在a中出现了几次
接下来m行正整数，分别代表b每次在a中出现的开始位置。按照从小到大的顺序输出，a下标从0开始。
【样例输入】
abc
a
【样例输出】
1
0

聪明的孩子一看就发现这是个构造。~~反正我是没看出来~~

然后我们构造 $s, t$ 数组。

对于 $s [i] = a [i], t [m - i] = b [i]$ （当然 $m - i$ 就代表当 $i = m$ 时为 $0$ ）

不过如果 $b [i] = ?$ ，那么 $t [m - i] = 0$ 。

那么我们只要求出 $f[i+j]+=(s[i]-t[j])^2*t[j]$ 。

如果 $f [k] = 0$ ，那么 $[k - m + 1, k]$ 与 $b$ 相匹配。

为什么？

因为如果 $s [i] = t [j]$ ，那么 $s[i]-t[j])^2*t[j]=0$ ，如果 $t [j] = 0$ ，那么也是 $0$ ，而且因为非负，只要有一个 $> 0$ ，式子就绝对不是 $0$ 。

但是要怎么算呢？我们只要算出 $s[i]^2t[j]$ 和 $s[i]s[j]^2$ 就行了。

而 $t[j]^3$ 完全可以独立处理。

#include<cstdio>
#include<cstring>
#include<cmath>
#define  N  310000
#define  pai  3.1415926535897932
using  namespace  std;
int  limit=1,l,r[N];
struct  node
{
    double  x,y;
    node(double  xx=0,double  yy=0){x=xx;y=yy;}
}a[N],b[N],c[N],d[N],wn[N];double  msin[30],mcos[30];
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  zwap(node  &x,node  &y){node  z=x;x=y;y=z;}
void  fft(node  A[],int  type)
{
    for(int  i=1;i<limit;i++)
    {
        if(i<r[i])zwap(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  L=0,R=mid<<1;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;
            }
        }
    }
}
char  s1[N],s2[N];
int  n,m;
int  ans=0,sta[N],top;
int  main()
{
    wn[0]=node(1,0);
    scanf("%s%s",s1+1,s2+1);n=strlen(s1+1);m=strlen(s2+1);
    for(int  i=1;i<=n;i++)a[i].x=s1[i]-'a'+1,d[i]=a[i].x*a[i].x;
    for(int  i=1;i<=m;i++)
    {
        int  x=(s2[i]=='?'?0:(s2[i]-'a'+1));
        b[m-i].x=x;c[m-i].x=x*x,ans+=x*x*x;//独立处理
    }
    while(limit<=n+m-1)limit<<=1,l++;
    for(int  i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    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);fft(c,1);fft(d,1);
    for(int  i=0;i<limit;i++)a[i]=a[i]*c[i],b[i]=b[i]*d[i];
    fft(a,-1);fft(b,-1);
    for(int  i=m;i<=n;i++)
    {
        int  x=int(b[i].x/limit+0.5)-int(a[i].x/limit+0.5)*2+ans;
        if(!x)sta[++top]=i-m;
    }
    printf("%d\n",top);
    for(int  i=1;i<=top;i++)printf("%d\n",sta[i]);
    return  0;
}

高维FFT

自己YY出来了QAQ（这种方法特别通用（不是我发明的，早就有人说过了），但是常数大！~~FFT哪有常数小的~~）。

首先，我们考虑二维（高维也一直是这样）。

也就是快速计算这个式子（ $0 \leq i \leq n, 0 \leq j \leq m$ ）：

$f(i,j)=\sum\limits_{q=0}^{i}\sum\limits_{p=0}^{j}f(q,p)f(i-q,j-p)$

也就是差不多是这么一个图像：

在这里插入图片描述

其实我们仔细思考是可以用一维FFT来做的，主要就是化编号。

原本绿色矩阵的 $(q, p)$ 在一维中的编号就是 $q * 2 m + j$ ，简单来说就是每一行扩大一倍，然后扩大的范围的通通为 $0$ ，然后一行拼到前面的一行。

然后就按FFT固定套路，多扩展一倍然后一维FFT走起，仔细一想我们的那个式子刚好乘起来的编号就是原本二维FFT乘完的编号转换过来的。

关键是这个做法还很容易扩展到高维，二维的时间复杂度为： $O(nm\log{nm})$

坑

3次变1.5次优化中试试系数是否能为复数（FFT支持复数系数，并且在2次DFT变1次中的IDFT优化已经体现），因为里面没涉及到什么系数是复数就不行的操作，而且还要实验这个优化对于 $P (x) = A (x) B (x) + C (x) D (x)$ ，处理了 $A (x) B (x), C (x) D (x)$ ，是否可以直接加起来然后一遍IDFT。（当然如果有巨佬否定了或者有实验肯定了，请在评论区告诉我，谢谢）。

zhangjianjunab

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
1
评论
快速傅里叶变换（FFT）学习笔记（未完待续）

总所周知，FFT是一个非常麻烦的算法，再加上博主语文不好，便写起来有点麻烦，但会尽力去写。要以后自己看不懂就。。。注：因为最近的压力紧张，便没有继续学习FFT，这仅为目前的半成品以及一些目前已知的优化。目录参考资料FFT吹水例题普通做法更高大尚的做法定义与一部分性质系数表达式点值表达式性质1点值相乘？？？卷积复数DFT参考资料讲解1讲解2优化毛啸大佬的论文（不放上下载地址，因为下载论...
复制链接

扫一扫