雪舞喵和FFT=。=

[UPD1]2017.01.12

于是今天雪舞喵听了一上午的fft,终于把迭代版预处理搞会了=。=然后就来更新一下

那个迭代版预处理其实就是观察一下预处理前后序列的性质=。=比如说原来一个8长度的序列, 下标分别是

0 1 2 3 4 5 6 7

递归预处理之后呢 重新排列原序列的下标 变成了这样

0 4 2 6 1 5 3 7

考虑一下把它们写成二进制形式

000 001 010 011 100 101 110 111 ==>
000 100 010 110 001 101 011 111

不就是把二进制形式下的高位变成低位,低位变成高位么=。=

换句话说 在二进制下倒过来读

这个很简单啊 只要先预处理出0~n-1的数字二进制反转后是多少,然后再对原序列应用这种变换就好了喵。

而且0~n-1的数字二进制反转可以线性预处理。首先0反转之后是0,1反转之后是n>>1

接下来 对于一个二进制数x(比如,0100101100),设它的反转为rev(x),我们考虑rev(x>>1)(是0110100100),只要把它>>1在把x的最低位填到它的最高位就可以啦

最后附上新(O(n))的狗带函数代码

 

顺便一提,原来那个递归狗带在codevs上交超大整数乘法是合计543ms,最慢175ms。改成这个之后变成合计337ms,最慢106ms了。


喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵(分割线)喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵喵


(^=。=^)于是经过昨天一天的艰苦奋斗,今天总算是把fft这个噩梦解决了=。=

以下是雪舞喵的个人理解,如果有不清楚的地方,原因一定是你是神犇,我是苣蒻,我们思维不是一个级别的。

所以呢,为什么要学FFT呢,因为我们要用它计算DFT,为什么又要计算DFT呢,因为我们要算两个序列的卷积。然而暴力算卷积是 \(O(n^2)\) 的……喵?卷积是什么?

嗯……大概就是这样,给定两个长度为n的序列A和B, 则它们的卷积为一个长度为n的序列C,并满足以下关系:

\[C_{k}=\sum_{i=0}^{k}A_{i}B_{k-i}\]

那么假如现在有两个n项的多项式 \(A(x)=\sum_{i=0}^{n}a_ix^i, B(x)=\sum_{i=0}^{n}b_ix^i\)

我们得到它们的系数序列 \(A=[a_1, a_2, ..., a_n], B=[b_1, b_2, ..., b_n]\) ,然后手动推一下(或者脑补一下)两个多项式相乘得到的系数序列C,可以发现C序列正好是A,B序列的卷积(只要在A,B序列后各补上n个0)。

因此,多项式乘法就是多项式系数序列的卷积。

然而暴力算卷积是 \(O(n^2)\) 的……为了优化,我们需要知道一个定理,叫做卷积定理。

这个定理的内容是这样的,设F(A)为A序列进行DFT后的结果,C序列是A和B序列的卷积,那么通过定理我们可以推出F(C)和F(A),F(B)的关系 \(F(C)_k = F(A)_kF(B)_k\) ,这个关系就简单多了啊,这样,我们如果知道了F(A)和F(B),就能在O(n)的时间内就能算出来F(C)。

至于证明雪舞喵一向非常头疼,如果需要请自行百度=。=

于是下一步来了,怎么计算A和B的DFT呢?

根据定义, \(F(A)_{k} = \sum_{i=0}^{n-1}\omega _{n}^{ik}x_i\)

这个 \(\omega _n\) 是一个叫做n次单位根的东西,他等于 \(sin(2\pi /n)+cos(2\pi /n)i\) ,定义就是它的n次幂等于1。

这是个复数啊,因此FFT中的所有计算(+,-,*)都是复数计算。

/*

**对于不知道什么是复数的,雪舞喵也科普一下吧

**首先呢,在实数域内,对一个负数开平方根是无意义的,这个就比较麻烦,总是需要乱七八糟讨论。

**数学家们和雪舞喵一样非常懒,觉的这讨论太麻烦了,就让他有意义吧,然后就发明了复数。

**一个复数由两个实数组成,称为实部和虚部,假设它的实部为a,虚部为b,那么它形式上表示为a+bi

**这个i是什么东西呢,可以理解为一个常量,定义就是 \(i^2 = -1\)

**因此他就类似一个二项式,那么两个复数(a1, b1)+(a2, b2)(用一种比较方便的表示方法)就等于(a1+a2, b1+b2),减法同理

**乘法(a1, b1)*(a2, b2) = (a1*a2-b1*b2, a1*b2+a2*b1)

*/

然而,这公式喵了个咪的不是 \(O(n^2)\) 的么!

不要急,接下来才是真正的FFT。

刚才说过FFT是个算法,他是用来计算DFT的,并且可以在O(nlogn)的时间内算出一个序列的DFT

不过他有两个小缺点。

小缺点待会再说,先研究一下他是怎么回事=。=

对于一个序列A,我们把他拆偶数下标元素组成的偶子序列E和奇数下标元素组成的奇子序列O,那么如果我们算出F(E)和F(O),就可以在线性时间内将两个DFT合并成为F(A)。这个证明起来比较容易

假定n为偶数,则 \(F(A)_{k} = \sum_{i = 0}^{n-1} \omega _{n}^{ik}x_i = \sum_{i = 2t}^{n-2}\omega _{n}^{ik}x_i+\sum_{i = 2t+1}^{n-1}\omega _n^{ik}x_i\)

\[=\sum_{t = 0}^{n/2}\omega _{n}^{2tk}x_{2t}+\sum_{t=0}^{n/2}\omega _n^{2tk+k}x_{2t+1}=\sum_{t=0}^{n-1}\omega _{n/2}^{tk}x_{2t}+\omega _n^k\sum_{t=0}^{n-1}\omega_{n/2}^{tk}x_{2t+1}\]

\[=F(E)_k+\omega _n^kF(O)_k\]

证毕。

以上证明用到了单位根的一些奇怪的性质,比如说 \(\omega _n^{2t} = \omega _{n/2}^t\) 之类的,它还有一条性质 \(\omega _n^k = -\omega _n^{n/2+k}\)

这性质有啥用呢,用来计算 \(F(A)_{k+n/2}\) 。求出的两个子序列长度只有原序列的一半,因此如果一次计算序列的一个元素的话需要扫描两个子序列两遍,这样就不能在一个大数组中搞了,本着节省空间的原则,我们同时计算出 \(F(A)_k\)\(F(A)_{k+n/2}\) 就可以覆盖 \(F(E)_k\)\(F(O)_k\) ,就可以不断的在一个数组里来回搞了。

由于要求n为偶数,而且我们需要递归解决问题,因此我们必须先把序列补到2的幂,这是第一个小缺点。

怎么补呢,在低位方向不断补0就可以了。(低位在后的话就往后插,低位在前就往前插)

这样我们就可以使用分治法了,概括一下FFT的流程:

1.划分问题:把给定的长度为2的幂的序列划分为偶子序列和奇子序列两部分。

2.递归求解:递归的计算偶子序列和奇子序列的DFT。

3.合并问题:把两个子序列线性合并。

显然FFT的复杂度是O(nlogn)

边界就是一个元素的序列DFT后还是它本身。

然而在实现上还有一个细节需要注意,就是我们需要对这个序列进行预处理。

为什么呢,因为我们如果正常的计算这个序列,比如说计算一个长度为8的序列,在算第0位时可以同时算出第4位,然而 \(F(E)_0\) 保存在第0位, \(F(O)_0\) 却保存在第1位,不能正好覆盖原来的计算结果。为了避免这个问题,我们操作一下这个序列,只要把偶子序列放在序列的左半边,奇子序列放在序列的右半边,在递归的预处理左半边序列和右半边序列。这样就可以避免这个问题啦。

同时这样处理之后我们还可以把FFT的主过程写成迭代版的。

虽然市面上流行版本的FFT代码这个预处理也是迭代版的,用了一些奇怪的位运算,把空间复杂度彻底压到O(n)了。

然而我并不会奇怪的位运算,就写了一个递归版预处理。(这样的话空间复杂度其实也是O(n),因为n+n/2+n/4+n/8+......=2n,不过new的时间常数稍大)(见下面代码)

接下来研究一下实现多项式乘法。

首先得到系数序列。如果给定的两个序列不等长,需要在较短的序列高位方向补0补到一样长。

接下来向后补到序列长度的二倍。

再向后补到2的幂。(这个后事实上是看你怎么保存这个序列,如果高位在前就向低位方向,如果低位在前就向高位方向,这两部补位只要保证算卷积的时候如果按顺序计算,必须先计算的是保存数的部分,在计算保存0的部分就可以啦)

然后计算两个序列的DFT。

根据卷积定理,将DFT对应项相乘。

计算结果序列的IDFT(即DFT的逆运算,定义是 \(F^{-1}(A)_k = \frac{1}{n}\sum_{i=0}^{n-1}\omega _{n}^{-ik}x_i\) ,实现上只需要把 \(2\pi\) 改成 \(-2\pi\) 就可以啦。

所以上代码了(我以大数乘法为例)

从代码中可以看出来,由于使用了double类型,精度误差是不小的,这是第二个小缺点。

然而这第二个小缺点是可以避免的,那就是使用NTT。

NTT是使用整数(而不是会掉精度的复数)在模一个素数的意义下进行的(因此为了保证结果正确,这个素数要尽量大)。

其实实现上和FFT几乎没区别,就是把复数运算换成带模数的整数运算,单位根换成一个整数而已。

这样的话我们需要选择一个数 \(a_n\) 来代替原来使用的 \(\omega _n\) ,选择什么样的数呢,只要满足单位根的3条性质就可以啦。

刚才已经说过两条了,其实第3条也说过,就是它的定义,需要保证 \(a_n^n=1, a_n^p \ne 1 (p \ne kn, k \in \mathbb{Z})\)

然后经过一系列复杂的推理证明 \(a_n=g^{\frac{M-1}{n}}\) 是可行的(其中M是那个素数,g是M的原根。注意NTT中所有运算全是在模M的意义下进行的,以下部分就省略了)。

/*

**关于原根:原根的定义是满足方程 \(g^{x}\equiv 1(mod\,M)\) 的最小整数解为 \(x=\varphi(M)\) 的所有g。对于一个素数就是x=n-1

**由于一个数n他的最小原根是loglogloglogloglogn级别的(6个log),所以可以直接从2枚举

**判断他是不是原根怎么搞呢,首先我们对M-1进行质因数分解,得到它的所有不同的质因数p1, p2, p3……pk

**只要判断出 \(g^{\frac{M-1}{p_i}}\not\equiv 1(mod\,M)\) 对于 \(i \in [1, k]\) 均成立即可。证明较复杂请自行百度。

**事实上我们可以直接暴力判断1到n-1

**因为这个素数肯定是写代码的时候就知道是多少(要么是给定的,要么是自己选择的)

*/

然而想让这个 \(a_n\) 可行我们对M还有一些要求。

首先M是个素数。

其次M-1必须能除开n,即 \(M=k2^p+1(k, p \in \mathbb{Z})\)

然后M还要尽量大(但是M的平方还不能爆long long)

因此如果你看一道题它的模数特别的奇怪(比如998244353之类的)(他等于119*2^23+1)而不是常规的1e9+7,1e8+7啥的,就大概需要FFT(然而知道这个并没有什么卵用,还是不会)

好了这样终于不用担心掉精度啦,只要把刚才的代码稍加修改即可~

其他的只需要把INTT中的1/n换成n的逆元,g换成g的逆元即可

好了,到这就差不多了吧,然而依然没有什么卵用,原因是FFT的题依然不会喵喵喵~

Yukimai(^=。=^)

仅有1条评论 发表评论

  1. shorn /

    前排资瓷

发表评论